An ALMA survey of the SCUBA-2 Cosmology Legacy Survey UKIDSS/UDS field: High-resolution dust continuum morphologies and the link between sub-millimetre galaxies and spheroid formation

We present an analysis of the morphology and profiles of the dust continuum emission in 153 bright sub-millimetre galaxies (SMGs) detected with ALMA at S/N ratios of $>8$ in high-resolution $0.18''$ ($\sim1$kpc) 870$\mu$m maps. We measure sizes, shapes and light profiles for the rest-frame far-infrared emission from these luminous star-forming systems and derive a median effective radius ($R_e$) of $0.10''\pm0.04''$ for our sample with a median flux of $S_{870}=5.6\pm0.2$mJy. We find that the apparent axial ratio ($b/a$) distribution of the SMGs peaks at $b/a\sim0.63\pm0.24$ and is best described by triaxial morphologies, while their emission profiles are best fit by a Sersic model with $n\simeq1.0\pm0.1$, similar to exponential discs. This combination of triaxiality and $n\sim1$ Sersic index are characteristic of bars and we suggest that the bulk of the 870$\mu$m dust continuum emission in the central $\sim2$kpc of these galaxies arises from bar-like structures. By stacking our 870$\mu$m maps we recover faint extended dust continuum emission on $\sim4$kpc scales which contributes $13\pm1$% of the total 870$\mu$m emission. The scale of this extended emission is similar to that seen for the molecular gas and rest-frame optical light in these systems, suggesting that it represents an extended dust and gas disc at radii larger than the more active bar component. Including this component in our estimated size of the sources we derive a typical effective radius of $\simeq0.15''\pm0.05''$ or $1.2\pm0.4$kpc. Our results suggest that kpc-scale bars are ubiquitous features of high star-formation rate systems at $z\gg1$, while these systems also contain fainter and more extended gas and stellar envelopes. We suggest that these features, seen some $10-12$Gyrs ago, represent the formation phase of the earliest galactic-scale components: stellar bulges.

The high star-formation rates are similar to those measured for local ultra-luminous infrared galaxies (ULIRGs, e.g. Sanders & Mirabel 1996;Tacconi et al. 2008;Engel et al. 2010;Riechers et al. 2011;Bothwell et al. 2013). The intense star-formation activity in local ULIRGs is believed to be triggered and fuelled by mergers, resulting in irregular morphologies at UV/optical wavelengths, with single and double nuclei and tidal tails (e.g. Clements & Baker 1996;Farrah et al. 2001;Surace et al. 2001;Veilleux 2002;Psychogyios et al. 2016). Theoretical models provide support for this suggestion: hydrodynamical simulations of mergers can result in remnants with a central starburst event building a bulge. After 1 Gyr the merger remnant comprises a central bulge with in situ star formation and an extended disc/envelope dominated by stars formed before the merger (Hopkins et al. 2013). It has been similarly suggested that major and minor mergers may also be the trigger for the activity in the high-redshift SMG population (e.g. McAlpine et al. 2019).
The spatial extents of local (U)LIRGs has been shown to vary strongly depending upon the observed wavelength: with the highest star-formation rate (U)LIRGs displaying the most extended emission in the optical, while at the same time showing the most compact emission in the mid-infrared, which is thought to trace the on-going star formation (Chen et al. 2010;Psychogyios et al. 2016). Optical depths effects are a likely explanation for these varying trends, and this suggest that the physical size measured in the optical is highly dependent on the geometry of the dust distribution (Calzetti et al. 2007;Psychogyios et al. 2016). Comparisons of the rest-frame optical and far-infrared sizes of highredshift SMGs have suggested similar behaviour, with much more extended optical sizes, compared to those derived from interferometric observations in the sub-millimetre, which trace the bulk of the star-formation activity visible in the rest-frame far-infrared waveband (Simpson et al. 2015b;Ikarashi et al. 2015;Hodge et al. 2016;Lutz et al. 2016).
However, while there are similarities, there are also apparent differences between the observed properties of SMGs and those of comparably strongly star-forming ULIRGs in the local Universe. One notable difference being the large spatial extent of the star-formation activity in the high redshift sources, which was hinted at in early interferometric studies (Chapman et al. 2004;Sakamoto et al. 2008;Ivison et al. 2012). This has now been clearly demonstrated by ALMA: while the the typical extent of the starburst seen in local ULIRGs is of the order of a few 100's pc to ∼ 1 kpc, the rest-frame far-infrared emission in high-redshift SMGs arises from a region with an effective FWHM of ∼ 2-3 kpc (e.g. Simpson et al. 2015b;Hodge et al. 2016). There are also hints that the dust continuum morphologies of some high-redshift SMGs show features which are not found in the typically more complex local counterparts. Thus recent high-resolution (0. 03-0. 3) studies with ALMA have found that the dust continuum emission in SMGs has a disc-like brightness profile (Simpson et al. 2015a;Hodge et al. 2016;Ikarashi et al. 2017;Gullberg et al. 2018). While a study of six SMGs at z 2.5 by Hodge et al. (2019) at 0. 07 resolution (∼ 0.5 kpc) has shown spatially resolved 870 µm dust continuum morphologies with "clump-like" structures bracketing elongated nuclear emission, reminiscent of bars and rings (Kormendy 2013). The sizes of the "bars" and "rings" are in the ratio of 1.9 ± 0.3, consistent with that expected for Lindblad resonances. If these are indeed bars and rings then analytic theory and numerical simulations (e.g. Binney & Tremaine 1987;Lynden-Bell 1996;Athanassoula 2003) have shown that the ring is formed by gas outside the point of co-rotation being driven outward, by angular momentum transfer, collect into a ring near the outer Lindblad resonance. At radii inside the point of co-rotation, however, the gas falls inwards to the centre creating the bar. A bar is a means to drive gas from the outer part of the galaxy towards the centre, as the incoming gas is robbed of its energy, due to shocks. Gas can be funnelled inwards by the bar over an extended period, so maintaining the star formation in the central region. Simulations have suggested, however, that a bar can also cause quenching of star formation in the central region by sweeping up the gas within the co-rotational radius within a few rotations, which is then consumed in a vigorous burst of star formation (Gavazzi et al. 2015).
Hence, while SMGs and ULIRGs share some similar physical characteristics, the difference in the extent of their star formation, and current view of the dust continuum morphologies leaves open the possibility that the star formation activity in the two populations are not driven by the same processes. Indeed, alternative theories have been proposed for how cold gas fuels the star formation in highredshift starburst galaxies and how the star formation is triggered, through accretion from the cosmic web (Bournaud & Elmegreen 2009;Dekel et al. 2009). In this scenario these galaxies rapidly accrete gas from the cosmic web and disc instabilities cause clumps to migrate to the nucleus where they form a bulge. These theories predict star-formation rates similar to the gas accretion rate of ∼ 100 M yr −1 and a resulting morphology of a gas disc twice the size of the nuclear star forming bulge (Dekel et al. 2009).
To improve our understanding of the dust continuum structures of strongly star-forming galaxies at high redshift and so throw light on their possible formation and triggering mechanisms, we have analysed the morphologies of a much larger sample of SMGs to those studied to date. In this paper we present the result of this high-resolution (0. 18) continuum study of the 870 µm morphology of 153 SMGs from the AS2UDS ALMA survey of sub-millimetre sources in the SCUBA-2 Cosmology Legacy Survey UDS field (Stach et al. 2019). This large sample of uniformly selected SMGs, with integrated continuum signal-to-noise ratios of ≥ 8, provides a statistically robust constraint on the sizes and, for the first time, the shapes of this high redshift population. By selecting only the highest resolution observations and applying a conservative signal to noise cut, we seek to go beyond measuring crude sizes for the SMGs and instead derive constraints on their profiles and axial ratios for large statistical samples. These can then be used to investigate the physical nature of the dust continuum emission in these systems. Our observations resolve the 870 µm dust emission in these sources and so provide reliable measures of the shape and profile parameters, such as the effective radius (R eff ), axis ratio distribution of the population and typical Sérsic indices.
An outline of the structure of this paper is as follows: in §2 we present the observations used for in our analysis. We analyse these and describe our basic results in §3. We then interpret and discuss these in §4, before giving our conclusions in §5. We adopt a standard concordance, flat ΛCDM cosmology of H 0 = 71 km s −1 Mpc −1 , Ω Λ = 0.73, and Ω M = 0.27 (Spergel et al. 2007).

SAMPLE SELECTION
Our sample is drawn from an ALMA follow-up study, called AS2UDS (Stach et al. 2019), of the sub-millimetre sources discovered in the SCUBA-2 Cosmology Legacy Survey map of the Ultra Deep Survey field (S2CLS Geach et al. 2017). Details of the AS2UDS observations, data reduction and catalog are given in Stach et al. (2019), although we briefly summarise these here. Using ALMA in Cycles 1, 3, 4 and 5 we targeted a complete sample of 716 single-dish SCUBA-2 850 µm sources with observed flux densities of S 850 > 3.6 mJy (corresponding to > 4σ detection significance in the SCUBA-2 map). For all the observations, the central frequency of the receivers was tuned to 344 GHz and the FWHM of the ALMA primary beam was 17. 3 (encompassing the FWHM of the SCUBA-2 beam of 14. 7).
To reduce the data, we used the Common Astronomy Software Application (casa, McMullin et al. 2007) version 4.5.3 using the standard ALMA calibration scripts. The data were imaged using the clean, algorithm in casa with natural weighting (ROBUST = 2). We cleaned the images to the 1.5σ level. Due to configuration differences during these cycles, the FWHM of the naturally weighted synthesised beam varies from 0. 18 to 0. 35 (with a small number of repeat observations obtained at 0. 7 in Cycle 5 to test if flux was being resolved out of the higher resolution maps, Stach et al. 2019). Hence to construct the catalogue, all of the maps were tapered to 0. 5 FWHM. The noise in these tapered maps varies between 0.09-0.34 mJy beam −1 (see Stach et al. 2019 for more information about the AS2UDS data reduction).
The final AS2UDS catalogue contains 706 SMGs that are detected at > 4.3σ (2% false positive rate) with a median flux density of S 870 ∼ 3.7 mJy (Stach et al. 2019). We note that in the tapered maps, on average we recover the majority of the single-dish flux for sources with S 870 3.5 mJy (Stach et al. 2019).
For this morphological study of the dust emission in SMGs, we concentrate on the Cycle 3 observations where a subset of 507 SMGs from the AS2UDS survey were detected in maps at a native resolution of 0. 18 FWHM. For relatively high resolution observations similar to these, Simpson et al. (2015b) showed that for a signal-to-noise ratio of S/N > 8, the uncertainties on the resulting size measurements of sources are 35% (Simpson et al. 2015b). We therefore select all 153 SMGs which were observed in Cycle 3 (0. 18 FWHM) and are detected with S/N > 8 in the 0. 5 tapered maps, and these form the sample for the remainder of our analysis. This selection should ensure we can measure robust sizes and shapes for these sources and that we are sensitive to a broad range in source sizes. This sample of 153 SMGs has a median flux density of S 870 = 5.6 ± 0.2 mJy, roughly ∼ 50% brighter than the full sample, with a range in flux density of 2.9-11.9 mJy, spanning the bulk of the SMG population which has been studied with ALMA. This means that though the exposure time of 40 seconds per source is short, we reach similar S/N levels as studies of fainter SMGs but with longer exposure times (e.g. Tadaki et al. 2017).
Our ALMA survey was carried out in the ∼ 1 degree 2 UDS field, part of which was observed with Hubble Space Telescope; (HST) by the Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey (CANDELS). (Grogin et al. 2011) This provides F606W, F814W, F125W, and F160W-band observations for 47 of the AS2UDS SMGs that lie in the central region of the UDS. We compare the optical/near-infrared and 870 µm dust continuum morphologies of some of these galaxies in Figure 1. This shows HST I JH-band colour thumbnails of eight ALMA SMGs, overlaid with the dust continuum emission from ALMA. At the median redshift of our sample (z ∼ 3), the observed I JH bands samples the rest-frame mid-UV to B-band, and as Fig. 1 shows, the rest-frame UV/optical morphologies display a range of structures on arcsecond-scales from disc-like to apparently multi-component mergers, point sources and SMGs which are undetectable in even the reddest HST filters. In contrast, on average, the 870 µm continuum appears much more compact than the rest-frame UV/optical emission, although generally the emission in the two wavebands is centred in the same position.

Multi-wavelength data sets and physical
properties of the sample Before we assess the dust continuum sizes of our SMG sample, for context we review the physical properties of our high-resolution sub-sample from AS2UDS and place them in context of the parent population of 706 ALMA SMGs in this field. In particular, (Dudzevičiūtė et al. 2019) use the extensive multiwavelength imaging of the UDS to estimate the photometric redshifts and physical properties of the complete sample, including inferring their stellar masses, starformation rates and dust masses. To achieve this, (Dudzevičiūtė et al. 2019) exploit 22-band photometry (or limits) for each SMG 1 , building on the UDS DR11 K-band selected catalogue of Almaini et al. (in preparation), and fit the spectral energy distribution, including deriving the photometric redshift estimates and uncertainties, using the high-redshift version of magphys (da Cunha et al. 2015;Battisti et al. 2019) From the analysis of the multi-wavelength SEDs, (Dudzevičiūtė et al. 2019) determine that the median redshift of the full sample of 706 ALMA SMGs is z = 2.61±0.08, with a quartile range of 1.8-3.4. The median star-formation rate determined for the full parent sample is 235 M yr −1 , the median dust mass is 6.7 × 10 8 M and the stellar mass is (1.3 ± 0.1) × 10 11 M (Dudzevičiūtė et al. 2019). The corresponding values for the sub-sample of 153 SMGs in our high-resolution sample are z = 2.9±0.1, with a quartile range 2.5-3.5, a median star-formation rate of 380 M yr −1 , a dust mass of 1.1 × 10 9 M and stellar mass of (1.3 ± 0.1) × 10 11 M . As expected, our high-S/N SMG sample, which are roughly 50% brighter in S 870 than the full sample, also exhibit correspondingly higher dust masses and star-formation rates and due to the correlation between observed S 870 and redshift reported by Stach et al. (2019) means that lie at somewhat higher redshifts than the full sample.We will return to a discussion of the trends of dust continuum structure with flux and star-formation rate in §4.

ANALYSIS AND RESULTS
In this section, we first assess the spatial extent of the dust continuum emission in our high-resolution observations of SMGs through measurement in both the the uv-amplitude plane, and fitting models to the image plane maps. We then derive the Sérsic profiles and axis ratios for the continuum emission.

Sizes measurements from the visibility plane
The spatial extent of the 870 µm continuum emission of our SMGs can be derived by measuring the amplitude as a function of uv-distance. We apply this approach to each SMG by first aligning the phase centre of our visibilities with the source positions from Stach et al. (2019) using the casa task fixvis and then radially averaging the amplitudes in 75 kλ bins (the choice of 75 kλ bins is arbitrary, although this binning minimises the scatter). In Figure 2 we show the real part of the amplitude as a function of uv-distance for the uv-range out to 1200 kλ for the same eight galaxies shown in Figure 1. The error bars on the amplitudes are given by the error on the mean in each bin. In this figure, we also include the total flux measurements from the 0. 5 resolution uv-tapered maps. Figure 2 shows that in all cases, the amplitude declines as a function of increasing uv-distance. This is a clear indication that the emission from the source is resolved in these observations. Fitting Gaussian light profiles we derive a median FWHM size of 0. 25 ± 0. 03.
A subset of our sample were observed at 1.1 mm with ∼ 0. 7 resolution using ALMA by Ikarashi et al. (2017) who studied a sample of millimetre sources selected from the 1.1-mm AzTEC map of the UDS field. Unsurprisingly the sources in this bright 1.1-mm sample overlap with brighter 870 µm sources in the S2CLS map and as a result 65 of the 69 sources in Ikarashi et al. (2017) are also included in AS2UDS, of which 30 are in our high-resolution 0. 18 subsample. We compare the ratio of the estimated sizes from uv-fits to the 870 µm and the lower-resolution (but S/N > 10) 1.1 mm observations for these 30 sources and derive a median ratio of FWHM from Gaussian fits of 0.95 ± 0.05. This provides strong independent confirmation of the reliability of our derived sizes using a completely independent observations, reduction and analysis method.
One noteworthy feature of Figure 2 is that it is clear that the amplitude does not converge to zero at large uvdistances in many cases. Indeed, in 119 SMGs (out of 153) the amplitude is non-zero (at > 3σ) at 1200 kλ. Naively this would suggest that a large fraction (∼ 80%) of sources  Fig. 1. The amplitudes are extracted by radially averaging the visibilities in 75 kλ bins over the full frequency range and the the total flux densities recovered in the maps uv-tapered to 0. 5 resolution are plotted as a square. We overlay half-Gaussian fits to the continuum emission as a dashed lines, and a Sérsic fit with n = 1 by the solid curve. The 870 µm dust continuum of the SMGs are all resolved in our observations. We note that the Gaussian fits frequently show an apparent compact or unresolved component, indicated by a non-zero flux at large uv-distances. However, this is an artefact of the fact that the profiles are poorly described by a Gaussian, while an n = 1 Sérsic better reproduces both the compact and extended emission. The 870 µm dust continuum sizes of the sources are listed in Table 2. Figure 3. Visibility amplitudes as a function of the uv-distance for three model profiles with n = 0.5, 1 and 2 and radius of 0. 15, observed with the same array configuration set-up as our observations, compared to the median composite profile of our 153 SMGs. This illustrate that for a source with an intrinsic Gaussian profile (n = 0.5) the amplitude quickly converges to zero, while for higher n the amplitude at the largest uv distances remains non-zero out to at least 1,000 kλ, mimicking the signal of a point-source component. This implies that for marginally-resolved galaxies with profiles steeper than n > 0.5 the inner part of the galaxy will appear unresolved, incorrectly suggesting the presence of a pointsource component.
contain an unresolved component (compared to a Gaussian model) comprising on average 13 ± 1% (or typically S 870 = 0.63 ± 0.05 mJy) of the emission.
Large uv-distances correspond to small physical scales, and so this "compact" emission seen in Figure 2 must arise on 0. 18 scales. One option is that the our observed uv-amplitude profile comprises a luminous, extended (Gaussian-like) starburst with a ∼ 13% contribution from a central point-source (Tadaki et al. 2017). However, it is also possible that the "compact" emission instead arises from a light profile which is more centrally concentrated than a Gaussian (which is generally a poor description of the light or mass profiles of resolved galaxies). To investigate how different light profiles should appear in the uv-amplitude plane, in Figure 3 we show the simulated uv-amplitude-distance as a function of Sérsic index (n) with n = 0.5, 1 and 2 and the median composite profile of the 153 SMGs in our sample. These profiles were created using the casa simulation tool with the same configuration as our observations (and hence the same synthesised beam FWHM). For brightness profiles with Sérsic index n > 0.5 the uv-amplitude profile does not converge to zero by a uv-distance of ∼ 1200 kλ -a consequence of the steep central light profile which gives rise to apparently compact emission. Hence, fitting a Gaussian model to a marginally resolved source with an intrinsic Sérsic n = 1 profile you would conclude that a second compact emission component was implied by the non-zero amplitude at large uv-distance. To illustrate the difference in Fig. 2 we overlay the best-fit models with n = 0.5 and n = 1 to the eight galaxies shown. The χ 2 distributions of the indi- Effective radii (R eff ) values measured in the visibility plane (with n = 1) compared to the circularised R eff measured using Sérsic fits with fixed n = 1 in the image plane. The R eff measurements from the two methods show a one-to-one correlation, with a scatter partly due to the difference in profile. Profiles close to n ∼ 1 show better agreement than profiles with n ∼ 2-5. B: Comparison of the circularised R eff (R circ ) measurements from the n = 1 and free Sérsic fits. For n 1.6 the radii for the n = 1 and free Sérsic fits agree within 1σ, while at n 1.6 the free Sérsic fit returns up to ∼ 2.5 times higher R circ than the n = 1 fit. C: Effective radii measured in K-band imaging as a function of R eff we measure with the free Sérsic fit. The points are coloured according to the measured axis ratios in the K-band images. When the sources appear close to circular, axis ratios 0.75-1, in the K-band image, the measurements show a rough correlation (indicated by the dash-dotted curve) with a relative size ratio (K to 870 µm) of 2.2 ± 0.2. More elliptical K-band measurements show a larger scatter, which (by assuming disc morphologies) perhaps indicate that the K-band imaging is more affected by extinction when viewed edge-on or the geometry of the dust and stellar light is different and hence varies as function of orientation. We return to this point in § 3.5 when we model the structure of the dust in these galaxies and show that the dust morphology does not appear to be disc like. D: Distribution of circularised R eff for the n = 1 Sérsic fit, with a median R eff of 0. 10 ± 0. 0.04. E: The distribution of the axis ratios for the n = 1 fit overlaid with the axis ratio from the free Sérsic fit. The n = 1 fit distributions has a median of ∼ 0.64 ± 0.24, and the free fit distribution has a median of 0.63 ± 0.24. F: The Sérsic index distribution of the free fit has a median Sérsic index of n = 1.00 ± 0.12. The majority of the SMGs have n < 1.5, however, the distribution shows a tail out to n ∼ 5.
vidual n = 0.5 and n = 1 fits show a moderate preference (2.5σ) towards the n = 1 fits. This is supported by the composite profile, which is best fit with a n = 1 profile. The n = 1 fits results in χ 2 ∼ 1.5, compared to χ 2 ∼ 2.5 for a n = 0.5 fit. As we show in § 3.2 fitting a Sérsic model to each of the SMGs in our sample suggests a median Sérsic profile of n = 1.00 ± 0.12. This indicates that the majority (77%) of the SMGs in our sample have 870 µm continuum profiles that are consistent with Sérsic with n ∼ 1 (rather than Gaussian, n ∼ 0.5) light profiles. We therefore conclude that Sérsic models provide an appropriate description of the 870 µm brightness profiles of the SMGs in our sample. Hence to measure the spatial extent of the dust continuum emission, we adopt an n = 1 Sérsic and allow the effective radius as a free parameter and fit the uv-amplitude profile for each SMG. We derive a median effective radius for the 153 SMGs in our high-resolution/high-S/N ratio sample of 0. 10 ± 0. 04.
The typical effective radius we measure is comparable to estimates from previous studies. For example, Hodge et al. (2016) measure the spatial extent of the 870 µm dust continuum in 16 SMGs from ALMA observations with ∼ 0. 16 resolution and a sample with flux density range of S 870 = 3.4-9.0 mJy. They derive a median effective radius of 0. 15±0. 03, which is similar to our sample. Simpson et al. (2015a) also measure the 870 µm sizes of 30 from our ALMA Cycle 1 observations of the brightest SMGs in the AS2UDS pilot (there is no overlap in sources to the 153 SMGs analysed here). Using that ∼ 0. 3 resolution data, Simpson et al. (2015a) derive a median effective radius of 0. 13 ± 0. 02. While, as noted earlier, Ikarashi et al. (2017) analysed ALMA 1.1-mm observations of 65 SMGs from AS2UDS and derive an effective radius of 0. 13 ± 0. 06 using their lower resolution, 0. 7 data. Thus it appears that the extent of the dust continuum emission we measure for our sample is comparable to that estimated from earlier smaller-scale studies.

Sizes and shape measurements from the image plane
The azimuthally averaged amplitude measured from the visibility plane is well suited for deriving a characteristic radius for the emission from a galaxy, but provides a circularised average. Information about the Sérsic index and the axis ratio (b/a) of the emission can also be derived, in a computationally more tractable manner, from the image plane.
To measure the sizes from the image-plane maps of each SMG we fit a two dimensional Sérsic surface brightness profile, allowing the effective radius, axial ratio and Sérsic index to vary. We account for the beam by convolving each intrinsic model with the synthesised beam with a semi-major and semi-minor axes and position angle given by the beam parameters for each map. The fit returns measurements of the peak flux, the central position, the semi-major axes, the axis ratio (i.e. the ratio of the minor to major axes, b/a), the position angle, and for the case of free Sérsic fit, the Sérsic index (n) as well as uncertainties on all these parameters.
Before we present the results of the fitting, we first test the reliability of the deconvolved measurements, and calibrate their uncertainties. To do this, we generate a set of 1,000 simulated galaxies using casa which have a flux distribution similar to our sample. These simulated galaxies have semi-major axis between 0. 11 and 0. 24, and random inclination angles and Sérsic index. We use casa to simulate the observations of these galaxies with the same exposure time as our data, and hence these simulated maps have similar noise properties as our sample. We then fit these simulated galaxies with our code and return their best-fit parameters. On average we recover the effective radii to within 20% and the axis ratios are recovered within 12% of the value of the input parameters. The Sérsic index is the most difficult parameter to fit and recover at the signal to noise of our typical sources, with a typical error of 22%. To investigate the potential bias due to noise when measuring the shapes of round sources (which may cause an observed decline in the apparent numbers of round sources), we also test the code on Gaussian profiles with axis ratio b/a = 1 and find that the procedure with free fit returns the Sérsic index within 25% and the axis ratio within 14%. The axis ratio distribution peaks for b/a = 0.86 and has a standard deviation of 0.12 for a sample matched in signal-to-noise to our observations. We therefore conclude that the strong peak in axis ratio at b/a ∼ 0.65 is our observed distribution if not a result of this bias.
The measurement of the Sérsic can also be influenced by data sampling (e.g. Robotham et al. 2017), and so next we investigate the influence of the reconstructed map pixel sampling by testing the same procedure with different pixel scales. We create 1,000 model maps of SMGs at the same pixel scale as our data (0. 03 pixel −1 ), and also at three times smaller sampling (0. 01 pixel −1 ). We simulate observations of these maps with casa and then fit these these maps and infer their properties. This test shows that finer oversampling of the synthesised beam does not return more accurate or precise Sérsic indices ∆n/n = 0.34 ± 0.05. For low signal-tonoise profiles with input Sérsic indices n 1.25 the fittingprocedure on the over-sampled maps return systematically lower values of n, but with increasingly larger uncertainties.
To derive measurements of the effective radii of the dust continuum in our SMGs, we now perform two sets of fits. First, we perform a fit with Sérsic n as a free parameter. In these fits, we derive a median n = 1.00±0.12 and R eff = 0. 11± 0. 01. Since the Sérsic index is the least certain parameter, we then fix the Sérsic index to n = 1.0 and refit each SMG. For this n = 1.0 fit, the median effective radius for the sample is R eff = 0. 10 ± 0. 04 (see Fig. 4).
In Fig. 4 we compare the effective radius for the SMGs derived from the uv-fitting with that derived in the image plane. We first compare the effective radii derived from the fixed n = 1 Sérsic fit in both cases, deriving a median ratio of the image-plane to uv-plane of R eff (uv)/R eff (im)= 1.10±0.01. Although the two measurements are correlated, it is also clear from Fig. 4 that when fitting a Sérsic model in the image plane, for larger effective radii, the uv-derived effective radius is ∼ 30% larger than that derived in the image-plane, but with no strong trend with Sérsic index.
In Fig. 4 we compare the effective radii of the dust continuum in the image-plane for the free and fixed Sérsic models. The median ratio in effective radii of R eff (free)/R eff (n = 1) = 1.07±0.01. The scatter in this relation can be attributed to those SMGs with profile with a higher Sérsic (as indicated by the colour scaling of the points, which show that for higher Sérsic index, the free fit returns larger sizes). We also show the Sérsic index distribution derived from the free fit, which has a median of n = 1.00±0.12. This shows that ∼ 70% of the SMGs in our sample have Sérsic indices n = 0.7-2 (as also suggested by the uv-amplitude profiles).
Using the image-plane fits, we extract the distribution of axis ratios from the best-fit models from both the free and fixed Sérsic model fits, and show these in Figure 4. Both distributions are strongly peaked, with a median of b/a = 0.63 ± 0.02 for the fixed n fit and b/a = 0.64 ± 0.02 for the free n fits, and both distributions have the same standard deviations of σ b/a = 0.19. We confirm that there is no correlation between measured b/a and effective radius. Hodge et al. (2016) present axial ratios for their 16 ALESS SMGs observed at 870 µm with 0. 16 resolution, which span a narrow range in b/a = 0.3-0.7 and have a median of b/a = 0.55 ± 0.06 and σ b/a = 0.13. Their distribution is thus consistent with that from our larger sample, and both display a relative lack of "round" sources (compared to the naive expectation of randomly orientated circular, thin discs, Ryden 2004). We will return to a discussion of this distribution in §4.
In summary, both the 870 µm dust continuum in the uv-amplitude profiles and image-plane maps suggest that the majority of the SMGs in our sample are best fit by n ∼ 1 Sérsic models. Fitting Sérsic models to the uv and image plane returns consistent results, and suggests a median effective radii of 0. 10 ± 0. 04 (∼ 0.8 ± 0.3 kpc at the median redshift of our sample). Since the Sérsic fit in the image-plane maps also allows us to easily investigate the shape parameters for the SMGs, for the remaining analysis we will use this method, but will adopt the minimum uncertainties on any measurement of effective radius from the scatter determine from the correlation(s) in Fig. 4.

Comparison of the dust and rest-frame optical emission
In Fig. 1 we show the HST images of eight example SMGs in our sample which are also observed as part of the CAN-DELS survey (Grogin et al. 2011). It is clear from this figure that the rest-frame optical emission is much more extended than the 870 µm emission (see also Simpson et al. 2015a;Chen et al. 2015;Lang et al. 2019). To quantify this further, we perform a Sérsic fit to the K-band images of all of the SMGs in our sample using GALAPAGOS (Almaini et al. 2017, Maltby et al. in prep) which takes the PSF into account, and we compare the effective radii measured in the UDS K-band compared to that measured from 870 µm in Fig. 4. This shows that for galaxies with K-band axis ratio of b/a = 0.75-1 (i.e. close to circular) the K-band effective radii are a factor of 2.2 ± 0.2 larger than the effective radii measured at 870 µm, although with considerable scatter (especially for those with large apparent axial ratios). This implies that the stellar light distribution is typically twice that of dust (see also Lang et al. 2019). We will discuss the extended emission in the stars and dust further in §4. We note, however, that optical depth effects need to be considered, in cases of more detailed comparisons of sizes measured at different wavelengths, for example between short wavelength dust emission, and [CII] or low or high-J CO emission.

Size and shape evolution
Previous morphological studies (with sufficiently high signalto-noise detections) of SMGs at sub-/millimetre wavelengths have been limited either by moderate resolutions (Simpson et al. 2015b;Ikarashi et al. 2017) or modest sized samples (Hodge et al. 2016;Gullberg et al. 2018;Hodge et al. 2019). Our sample of 153 SMGs detected at S/N > 8 in 0. 18 resolution maps, allow for a statistical study of the morphology, sizes and axis ratios for a wide range in both 870 µm flux density (S 870 ∼ 2.7-11.5 mJy) and redshift (z ∼ 1-6). Our sample has a median redshift of z = 2.9 ± 0.1 and a quartile range in redshift of 2.5-3.5, so the corresponding rest-frame wavelength (λ rest ≥ 150âȂŞ-300 µm) of the dust emission from the galaxies should be generally optically thin (e.g. Simpson et al. 2017). However, at the highest redshift, the dust emission may be optically thick (e.g. at z 3.5-4, the rest-frame wavelength is λ < 200 µm), this would then produce a small bias such that galaxies whose orientation provides a larger apparent sky area (i.e. face-on for disc-like geometries) would have brighter S 870 . There is a possible hint of this in Figure 5 where we plot the axis ratios as a function of both 870 µm flux density and redshift (Dudzevičiūtė et al. 2019), and find weak positive correlations.
In the plot of effective radius (R eff ) versus 870 µm flux density, we identify a weak positive correlation. But this trend is marginal and so we conclude that the brighter SMGs are more luminous primarily due to their higher dust surface densities, rather than their larger sizes. The apparent trend suggests a doubling of dust mass surface density between SMGs with S 870 ∼ 4-11, which may imply a similar increase in gas density and a corresponding rise in mid-plane pressure in these systems, which are already thought to be extremely high (Swinbank et al. 2011(Swinbank et al. , 2015. More interestingly, we see a small variation of the effective radius of ∼ 10% (corresponding to ∼ 2.5 × 10 8 M kpc 2 ), by a weak decline with redshift (Dudzevičiūtė et al. 2019). We expect this behaviour to reflect both evolution in the physical size of the sources and also the influence of dust optical depth, dust temperature and source structure. For simple source geometries, the apparent sizes of sources are expected to decline with increasing L FIR /M d (Falcón-Barroso & Knapen 2013). Using the estimates of L FIR and M d for our sample from (Dudzevičiūtė et al. 2019), we expect an increase in median L FIR /M d of only ∼ 20% in the sample across z ∼ 2-4 (Falcón-Barroso & Knapen 2013), which would correspond to a ∼ 20% decline in apparent size at a fixed rest-frame wavelength. However, this will be countered by the effect of shifting to shorter (and hence optically thicker) rest-frame wavelengths as we observe higher redshift sources. Thus we expect the drop in observed 870 µm effective radius with redshift to be less than ≤ 20%, which would be consistent with the weak decline seen in Fig. 5.
Overall, we conclude there are a number of potentially competing effects which could influence the variation in apparent size of the SMGs with redshift, but none of these effects is strong and hence the absence of significant evolution in the effective radius with redshift most likely indicates that the intrinsic physical sizes of the SMGs do not evolve strongly with redshift. This would be in contrast with studies in the UV and optical of a variety of galaxy populations which have reduction of a factor of several in typical size with redshift for both quiescent galaxies and star-forming galaxies across z ∼ 0-4 (e.g. van der Wel et al. 2014;Shibuya et al. 2015;Kubo et al. 2018).

Modelling the axis ratio distribution
Our results above, as well as recent studies (e.g. Simpson et al. 2015b;Hodge et al. 2016), have shown that the 870 µm dust emission in SMGs follows an exponential surface brightness profile suggestive of a disc-like geometry. For a sample of circular exponential discs viewed at random viewing angles, the axis ratio distribution should be constant at high axial ratios, with a decline towards b/a = 0, the strength of which depends upon the relative thickness of the disc. However, as shown in Fig. 4 for our SMGs, the apparent axial ratio distribution is highly peaked, with a large fraction of the sample (62 ± 4%) having axis ratio in the range b/a = 0.5-0.8, and proportionally fewer at b/a < 0.5 and b/a > 0.8, where the latter account for 15 ± 8% of the sample. This "deficit" in the number of sources with near-circular shapes suggests that either the assumption of SMGs being circular discs with a uniform viewing angle distribution is incorrect, or that signal-to-noise effects cause us underestimate the numbers of high b/a sources.
We first confirm that the behaviour we see is not due to signal-to-noise effects. To ensure that our profiles are robust, we set a lower limit for the signal to noise of S/N > 8, since the fractional uncertainties on the measured sizes increases to over ∼ 35% at lower signal-to-noise ratios (Simpson et al. 2015b). The signal-to-noise ratio of our sample ranges from These data are observed 870 µm meaning that the rest-frame wavelengths vary across the redshift range z ∼ 1-6, and hence the effective optical depth is expected to vary at the observed wavelength. The influence of these optical depth effects, as well as potential evolution in the structure and physical properties of the SMGs (e.g. dust mass and far-infrared luminosity) mean it is hard to draw strong conclusions from this plot. Nevertheless, we suggest that the lack of any strong trends in effective radius with redshift most likely indicates that there is no strong evolution in the size of the SMGs with redshift. 8 to 29 with a median of 12. Dividing the sample in half at S/N 12 yields two distributions which both show a "deficit" of axis ratios at > 0.8, suggesting that the signalto-noise is not the cause of the observed "deficit" of high axis ratio sources. Hence we now explore three geometrical models for the sources to attempt to reproduce the observed axis ratio distribution (including the influence of potential selection effects): optically thick discs, optically thin discs and a model with a triaxial geometry for the sources.

Geometrical models
To assess the possible influence of selection effects on the axis ratio distribution we compare our observed axis ratio distribution with those predicted for sources which are modelled as optically thick or thin circular discs. We generate a simulated sample of circular discs, where the apparent axis ratio is only dependent on the (random) viewing angle. We follow the example of Ryden (2004), where the apparent axis ratio q is given by where, A, B and C are given by A = [1 − (2 − ) sin 2 φ] cos 2 θ + γ 2 sin 2 θ, B = 4 2 (2 − ) 2 cos 2 θ sin 2 φ cos 2 φ, Here, is the ellipticity of the source = 1 − b/a (where a and b are the intrinsic major and minor axes), and γ is the ratio between the third axis (c) and the major axis; γ = c/a. The two angles θ and φ are the two viewing angles. Only θ has an influence on the apparent axis ratio in the case of a circular disc, where a and b are equal. The resulting distribution of apparent axis ratios for circular discs is therefore dependent on the distribution of the viewing angle θ, the c/a ratio, and the flux distribution. We assume a flux distribution similar to our sample and calculate the apparent axis ratios for the possible combinations of the different parameters, adopting a uniform distribution for the viewing angle and for c/a. Model I. Optically thick disc -In the case of an optically thick disc, at an observed wavelength of 870 µm, the fraction of the emitted emission that is detected by the observer is given by the visible fraction of the surface area (i.e. the apparent axis ratio).
We attempt to fit the distribution of b/a with Eq. 1 Figure 6. Distribution of the axis ratios for n = 1 compared with the expected axis ratios from a) optically thick, b) optically thin disc models and c) from a triaxial model. For the axis ratio distribution of the optically thick and thin models to have a similar "deficit" of high axis ratios (face-on discs), we would have to assume an unphysical model where we are preferentially viewing the sources at a particular viewing angle e (∼ 40 • ). In contrast the triaxial model (with a uniform viewing angle distribution) can broadly reproduce the observed axis ratio distribution for our sample of SMGs. Right insert: In the best-fit triaxial model the intrinsic b/a and c/a distributions that result in the observed axis ratio distributions are Gaussian distributions with a peak at 0.68 ± 0.02 and a width of 0.12 ± 0.06 for b/a and a peak at 0.28 ± 0.01 with a width of 0.19 ± 0.01 for c/a.
for this optically thick model, but in all cases the best-fit is poor (see Fig. 6). The best-fit optically thick disc models with a uniform viewing angle distribution has a Gaussianshaped c/a distribution, peaking for c/a = 0.06 with standard deviation of 0.33. However, the best-fit is still poor, and a Kolmogorov-Smirnov (KS) test shows a negligible (0%) chance for these two distributions to be drawn from the same parent sample. While it would be possible to bias the viewing angle distribution and find a better fit, such non-uniform viewing angle distributions are unphysical in the absence of an identifiable selection bias as a cause, and so we discard this option. 2 Model II: Optically thin disc -For the optically thin case the emission only depends on the assigned source brightness. As for the optically thick case, we are unable to find acceptable fit to the observations for the optical thin model with a uniform viewing angle distribution. The closest model for the optically thin case, with a uniform viewing angle distribution, has a rise in axis ratio distribution beginning at lower axis ratios, and has a flatter and slower rise than in the optically thick case. The best-fit model with a uniform viewing angle distribution has a c/a distribution peaking at 0.09, with a standard deviation of 0.28. Again a KS test returns a negligible chance (0%) of the two distributions being drawn from the same parent sample. The only way to match the "deficit" at high axis ratio would then be to include a biased viewing angle distribution, which as noted earlier is not physically plausible.
Model III: Triaxial structure -The shape of our observed axis ratio distribution differs from those seen for late-type spiral galaxies (Ryden 2004), which lack the strong peak at b/a ∼ 0.6 and the associated deficit at high axis ratios which we observe. We have quantitatively confirmed this above and conclude that neither of the circular disc models is able to adequately fit the data without invoking assumptions of unphysical viewing angle biases. Thus we now explore a triaxial model, where a > b > c. In this model we determine the best fit b/a and c/a axis ratios distribution to fit the observed axis ratio distribution for an optically thin case. For the triaxial case both viewing angles (θ and φ) have an influence on the apparent axis ratio q. We assume uniform angle distributions for both θ and φ, and model the apparent axis ratio distribution for a range of model and width parameters for the b/a and c/a distributions. We assume that the b/a and c/a distributions follow Gaussian distributions given by where q int represents the intrinsic b/a and c/a axis ratios, q 0 is the mode value, and σ is the width of the distribution. We calculate axis ratios for two different ranges of mode and width parameters for b/a and c/a (see Table 1). As the distributions are expected to be continuous and smooth, we perform our search for the best-fit model by The profile shows that most of the emission is compact, but that 13 ± 1% of the emission is extended on larger scales (≥ 0. 5 or ≥ 10 kpc), and only detected in our maps by stacking. The RMS level of the image stack is illustrated by the dotted line. Right: The 870 µm continuum light profile obtained by stacking the 870 µm dust maps in the visibility and image planes, normalised and compared to the stacked stellar and molecular gas profile seen in HST and ALMA imaging (Calistro Rivera et al. 2018). While the 870 µm dust continuum emission has an excess in the nuclear region compared to the stellar and molecular gas profiles, the extended component follows that of the stellar and molecular gas, suggesting that they trace the same structural component in these systems.
running through the parameter space twice; first by using a large bin size to cover the full parameter range and select the combination of parameters that provide the closest match to the observed axis ratio distribution; this best-fit parameter combination is then used as the centre of the second run which uses a finer search grid. The parameters for the b/a and c/a distributions that result in an apparent axis ratio distribution best fit to our observed distribution (and see Fig. 6) are given in Table 1 (see also Fig. 6). A KS test shows that there is a 40% chance that the triaxial axis ratio distribution and our observed axis ratio distribution originate from the same parent sample. This is further supported by the Akaike Information Criterion (AIC), which takes into account the number of fitted parameters, and for the triaxial case is 10.8. The AIC values for the optical thick and thin cases with uniform viewing angle (discussed in §4.4) distributions are 33.8 and 15.2. The model resulting in the lowest AIC values yields the best-fit model, which is therefore the triaxial model. We conclude that the "deficit" in high axis ratios is most likely due to intrinsic triaxial morphologies, rather than the dust continuum emission of SMGs resembling randomly orientated circular disc galaxies.

Stacked emission profiles
Near-infrared HST imaging of ALMA SMGs shows that the rest-frame UV/optical emission is extended on ∼ 8-10 kpc scales (e.g. Fig. 1; see also Chen et al. (2016); Lang et al. (2019)). In comparison, as we showed above, the 870 µm dust continuum is much more compact, with an effective radius of just 0. 10 ± 0. 04 or ∼ 1 kpc. But is a more extended, lower surface brightness, emission component also present? Stach et al. (2019) showed that our 0. 18 ALMA resolution observations of the 870 µm emission from the SMGs recover ∼ 95% of the single-dish flux detected with SCUBA-2. This may indicate that a small fraction of the flux is resolved on larger scales, but we do not have the sensitivity to detect this extended emission on a case-by-case basis. Instead to search for this emission we can stack the 870 µm continuum maps of the SMGs. As with our size measurements, we perform this stacking in both the visibility and image planes to assess the reliability of our results.
First, we stack the SMGs in the visibility plane, which has the advantage of circumventing issues arising from inhomogeneous beam-sizes of the individual maps. We shift the phase centre of the ALMA primary beam to correspond to the position of the SMG (all of our targets are the sole detected SMGs in their maps) and employ the stacker library developed for use in casa to stack the visibilities (Lindroos et al. 2015). The resulting stacked visibilities are then imaged using casa. Since the SMGs in our sample have a range of flux and signal-to-noise ratio, we stack the data weighted by 1 / σ 2 . From the stacked visibilities weighted by 1 / σ 2 , we measure the flux as a function of radius, and show the resulting surface brightness profile as a function of radius in Fig. 7. This figure shows a resolved, high surface brightness central region, but we also clearly detect faint and extended emission on ≥ 0. 5 scales, with an integrated flux that is ∼ 10% of the total flux.
To assess the sensitivity of the derived properties of the extended component on details of the data processing, we also derive a stacked profile by combining the individual image-plane maps of the SMGs for comparison. We extract a 10 × 10 thumbnail centred on each SMG, and then average the thumbnails, weighted by 1 / σ 2 . We then again extract the surface brightness profile from the stacked map and overlay this on to the profile created from the uv-stack in Fig. 7. We apply the same procedure to the calibrator, and also overlay this in the same figure. The image-planeand uv-derived stacks are well matched, with both showing the same extended emission on ≥ 0. 5 scales. We stress that the profile of the (point-source) calibrator, which we scale to the same peak surface brightness, is much narrower than the compact emission and lacks the faint emission halo we see in the SMGs. To test the that this extended emission is not due to weak side-lobes in the individual maps, we use casa to simulate a source with an n = 1 light profile and a size and axis ratio equal to that of our sample, but 10 times brighter. This allows us to test if the extended emission we see in the stacked imaged is due to weak emission on large radii from a compact source. This test show that a compact source with R eff = 0. 11 observed at 0. 18 resolution does not show emission on ≥ 0. 5 scale. This imply that the extended emission we detect in the stacked image is indeed likely to be caused by a second component.
To characterise the surface brightness profile, we fit a two-component model, including an inner and an outer Sérsic profile each with n = 1. For the compact component we fix R eff to the value derived above, 0. 10 ± 0. 04, and for the extended component we adopt R eff ∼ 0. 5. The extended component accounts for 13 ± 1% of the total emission. This suggest that the SMGs generally comprise a centrally concentrated starburst which accounts for ∼ 90% of the total dust continuum flux density, with an extended star formation component on scales similar to that seen in the restframe UV/optical (as seen by HST ). The transition between the compact and extended components occur at 0. 15 ± 0. 01 (∼ 1.2 kpc at z ∼ 3), and the luminosity-weighted average effective radius of the two components, which provide the most appropriate "size" for the whole systems in 0. 15 ± 0. 05 corresponding to 1.2 ± 0.1 kpc. However, in terms of relative surface brightness -the extended component has a peak surface brightness (and hence implied dust mass surface density) which is around two orders of magnitude lower than the compact component.

Discs, spheroids or bars?
To provide a qualitative context for the structural properties of the SMGs in our sample, we compare in Fig. 8 the observed axis ratio and Sérsic index distributions from the 870 µm observations of the SMGs to those derived in the gband for the stellar light distributions of a morphologically classified sample of low-redshift galaxies from the GAMA survey (Kelvin et al. 2014). Consistent with the results from our modelling in the previous section, Fig. 8 shows that the apparent axis ratio distribution of disc galaxies is not a good match to that observed for the SMGs, with a 0% chance that they are drawn from the same parent population. Whereas the SMG distribution is better matched to that of the spheroids, with a KS test returning a 30% chance of the two samples sharing the same parent population.
However, Fig. 8 also shows that spheroids have light profiles which are better described by high Sérsic indexes. This is not the case for the SMGs, which have a median Sérsic index of n = 1.00 ± 0.12, with a distribution which differs significantly from the spheroids (a KS test returns a 0% probability that the two distributions are from the same parent sample). Whereas the Sérsic index distribution for late-type disks is a much closer match to that seen for the SMGs, with a strong peak at n ∼ 1, and a tail to higher n.
This combination of apparently triaxial structures with exponential surface brightness profiles resembles that seen in the central bars of barred spiral galaxies (e.g. Seigar & James 1998). This interpretation of the morphology of the dust emission in SMGs was first suggested by Hodge et al. (2019). They re-observed six of the ALESS SMGs from Hodge et al. (2016) with ALMA in deep integrations at 0. 07 resolution and identified complex structures which were unresolved in their earlier 0. 15 observations (comparable to the data analysed here). In their higher resolution and deeper maps, Hodge et al. (2019) find symmetric clump-like structures bracketing elongated nuclear emission. They interpret these morphologies as representing bars in galaxies where the "clump"-like structures are formed through orbit crowding or star-forming rings. Hodge et al. (2019) find a ratio of the diameters of the bar-to-ring structures of 1.9 ± 0.3 similar to that seen for these components in local barred galaxies. Our results on a larger statistical sample, while lacking the resolution and depth to directly detect these features, have structural properties consistent with the suggestion that much of the dust continuum emission from SMGs arises in bar-like structures in their central regions.

What is the extended dust component?
Our analysis of the stacked profiles of the SMGs in our sample indicates the presence of a spatially-extended component with a peak surface brightness which is roughly two orders of magnitude fainter than the compact component detected in the individual sources, and which contributes ∼ 10% to the total flux densities.
To investigate the relationship between the compact and extended dust continuum emission, we split our sample into four bins of star-formation rate (with equal numbers of SMGs in each bin) and stack the maps of these sources in the image plane. We show the resulting surface brightness profiles for the four independent subsamples in Fig. 9. Each of these profiles shows both a compact and an extended component. We fit these with the same double Sérsic model used above to derive the fraction of luminosity in the compact and extended components. We then plot this ratio as a function of star-formation rate in Fig. 10. This figure demonstrates that the luminosity density of the extended emission remains approximately constant, despite the luminosity density of the compact component increasing by a factor of 50% (over a range of a factor ∼ 6 in total starformation rate), suggesting that the star-formation surface density in the compact and extended regions are decoupled.
We also compare our measurements of the extended component with the median value of a sample of faint field galaxies in strongly lensed cluster in the Hubble Frontier Fields Survey (González-López et al. 2017;Laporte et al. 2017). These galaxies at z 1.0-2.9 are detected with Figure 8. A comparison of the distribution of axis ratios and Sérsic indices for 153 SMGs, to those measured at rest-frame optical wavelengths for morphologically-classified late-type discs and early-type spheroidal galaxies from the GAMA survey (Kelvin et al. 2014).
Left: Axis ratio distribution of our SMGs compared to that of the late-type disc and spheroids, with the distribution scaled to the peak of our SMG distribution. A KS test shows that our SMG distribution is most similar to that of the spheroidal galaxies, with a 30% chance that these are selected from the same parent sample (compared to 0% match to the late-type discs). This supports the validity of the triaxial model for the morphologies of the SMGs proposed in § 3.5.1. Right: Sérsic index distribution of our SMGs compared to that of late-type discs and spheroids, now showing a stronger similarity between the SMGs and the late-type disc galaxies, rather than the spheroids (which are rejected as a match based on KS test returning a 0% probability that they are selected from the same parent sample). These two plots thus suggest that the typical 870 µm dust emission of our SMGs has a triaxial morphology, but an exponential surface profile -these characteristics are seen for central bars in disc galaxies. Figure 9. Surface brightness as a function of radius for the image stacks split into bins of star-formation rates of 40-250, 250-350, 350-500 and 500-1500 M yr −1 . Both stacks in the image plane and visibility plane return similar profiles, which are composed of a compact and an extended component (we therefore only show the image stack). We fit double exponential profiles and find that the peak flux of the extended component is close to constant at a mean of 0.47 ± 0.03 mJy for the four stacks, while the compact component becomes brighter with increasing star-formation rate.
ALMA at 1.1 mm and represent more typical star-forming field galaxies with star-formation rates of 10-100 M yr −1 . Interestingly, the median value for these faint field galaxies show similar surface brightness and extent to the extended components seen in the SMGs, and suggests that the extended star formation we detect in SMGs has similar starformation surface density to that in "normal" star-forming galaxies in the field. In contrast, the intense star-formation surface density in the compact component is significantly higher.
We also wish to understand the relationship of the extended structure we have uncovered at 870 µm to the other baryonic components of these systems: the (unobscured) stellar emission and the cool molecular gas. To demonstrate how the observed near-infrared emission (rest-frame UV/optical at z ∼ 3) compares to the 870 µm emission, in Fig. 7 we overlay the surface brightness profile of the dust continuum emission and the near-infrared from HST (in this figure, we have scaled the surface brightness profiles to the same integral). We also overlay the gas emission as inferred from the molecular 12 CO(3-2) emission from four SMGs from ALESS (Calistro Rivera et al. 2018). This low-J CO transition is expected to arise from material in the interstellar medium which has low to moderate critical densities. This means its extent should trace the bulk of the underlying cool gas reservoir in these systems. As can be seen in Fig. 7, the extended emission seen at 870 µm from our stack seems to match the spatial extent of the sources as seen in the near-infrared and also in the low-J CO from the molecular gas. Since the extended 870 µm dust continuum component follows the same profile as the molecular gas and (less obscured) stellar emission, this suggests that the extended component traces a halo or outer disc in these galaxies, which are dominated by stars and with much lower star-formation surface densities and obscuration than the central starbursts.
Similar two-component profiles to that we see in our 870 µm stacking analysis have also been observed through stellar-mass surface-density profiles in high-resolution hydrodynamical simulations of merging high-redshift massive starburst discs with properties similar to those of SMGs (Hopkins et al. 2013). The merging galaxies in these simulations are initially disc-dominated, but form nuclear bulges (on kpc scales) dominated by in situ star formation fuelled by gas driven to the centre by strong torques (Mihos & Hernquist 1994;Hopkins et al. 2008). Using dynamical arguments Hopkins et al. (2009) suggest that because gas can dissipate energy, it can efficiently lose its angular momentum and rapidly fall into the centre in a merger event. This results in a concentrated starburst event seen in, for example, nearby merging ULIRGs and recent merger remnants (e.g. Scoville et al. 1986;Sargent et al. 1987Sargent et al. , 1989Kormendy & Sanders 1992;Rothberg & Joseph 2004). This process builds a clear bulge in the centre of the merger remnant.
The extended component or "envelope" in this simulation is dominated by stars formed before the merger event and gas at large radii with significant angular momentum. This gas does not lose its angular momentum in the merger and re-forms a disc as the remnant relaxes. The survival (or re-formation) of the disc is therefore dependent on how much gas loses its angular momentum (Hopkins et al. 2009). If all the gas in a merger event efficiently loses its angular momentum, it would all be consumed in the nuclear starburst, and no gas would be left to re-form a disc. However, high gas fractions have been shown to be inefficient at losing their angular momentum, leading to the fraction of gas available to fuel the central starburst scaling sub-linearly with the gas fraction, and so leaving gas to re-form a disc (Hopkins et al. 2013).
In the framework of the model developed by Hopkins et al. (2013) we can also investigate the physical nature of the SMGs and their triggering. Besides considering different models with feedback and effective equations of state, Hopkins et al. (2013) also consider the influence of prograde versus retrograde mergers of disc galaxies. The different feedback models and the relative angular momentum vectors of the discs have little influence on the stellar mass profile of the remnant. However, Hopkins et al. (2013) also show that the relative angular momentum vectors and orbit of the merging components has an influence on the time scales with which the merger remnant evolves. A prograde merger develops a morphology with a nuclear starburst and an envelope af- Figure 10. Fitted peak surface brightness as a function of the star-formation rate for the compact and extended components (in Fig. 9). We see that the peak surface brightness for the extended component is roughly constant, while the brightness of the compact component increases as a function of star-formation rate. This is illustrated by the linear fits. We compare the extended components to the surface brightnesses of faint dusty galaxies in the Hubble Frontier Fields, and find a comparable median surface brightness. Hence the extended component we find in the SMGs are comparable in dust surface density to "normal" star-forming field galaxies.
ter ∼ 1 Gyr, while the morphology of a retrograde merger after ∼ 1 Gyr still shows two separate disc galaxies. This suggest that if major-merger events are the cause of the nuclear starburst event that we observe for our SMGs, a large fraction will have to be remnants of, or late stage, prograde mergers. If they were retrograde then the median redshift of our sample of z ∼ 2.9 this means that the two merging disc galaxies would have to have appeared as highly star-forming systems at redshift 5. The redshift distribution of SMGs (Chapman et al. 2005;Weiß et al. 2013;Simpson et al. 2014;Dudzevičiūtė et al. 2019) has been shown to peak at z ∼ 2.5-3.5 and to have a tail out to redshifts z ∼ 7, meaning that starburst galaxies at z 5 do occur, but are rare in 870 µmselected samples. Hence we suggest that SMGs at z ∼ 2-3 are more likely to be late-stage mergers, rather than merger remnants.

SUMMARY AND CONCLUSIONS
We analyse the dust continuum morphologies and light profiles of 153 well-detected (S/N > 8) SMGs observed with ALMA at 0. 18 (∼ 1 kpc) resolution. We fit both Gaussians (in the visibility and image plane) and free and n = 1 Sérsic models (in the image plane), and measure the effective radii, axis ratio and Sérsic indices for the individual sources. We also stack (again in both the visibility and image planes) the 870 µm emission for the whole sample and selected subsets to trace fainter and more extended emission around these systems. Our main conclusions are: • The median effective radius for SMGs in our sample is 0. 10±0. 04. Accounting for the extended dust component we find, we derive a flux-weighted effective radius of 0. 15±0. 05 or 1.2 ± 0.4 kpc at z ∼ 3. This in consistent with estimates of the sub-millimetre sizes of SMGs from earlier smaller-scale studies. We show that there is a rough correlation between the 870 µm and observed K-band sizes of SMGs in our sample, with the 870 µm sizes being on average 2.2 ± 0.2 times smaller.
• The effective radii of SMGs in our sample show a very weak decline with increasing redshifts and a similarly marginal increase with S 870 . Using the physical properties of our sample from Dudzevičiūtė et al. (2019) and assuming a simple source geometry we expect that the typical apparent source size would decrease slightly at higher redshifts, although this evolution would be countered by optical depth effects. The weak decline in size we see with redshift suggests that the physical sizes of SMGs do not evolve strongly with redshift and the lack of variation in size with S 870 indicates that the more luminous systems are likely to exhibit higher pressures in their interstellar medium.
• We find that the apparent axis ratio distribution of the SMGs is best described by non-axisymmetric morphologies (triaxial) and the Sérsic index distribution has a median of n = 1.0 ± 0.1. By comparing these distributions with those of disc and spheroid galaxies, the axis ratio distribution of SMGs is most similar to those of spheroid galaxies, while their Sérsic index distribution is most similar to that of disc galaxies. This combination of exponential surface brightness profiles and triaxial structures are the characteristics of bars in galaxies. Higher resolution and deeper observations of six SMGs by Hodge et al. (2019) have identified potential bar and ring structures in those galaxies and we therefore suggest that the statistical properties of the SMGs in our sample point to bars being a ubiquitous feature of bright SMGs.
• We stack our SMGs in both the image and visibility planes and find that the continuum emission profiles are composed of not only the compact component we have directly detected, but also a much lower surface brightness, extended component. The extended component accounts for 13 ± 1% of the total emission and has a scale size of ∼ 0. 5 (∼ 4 kpc) Comparing with stacked CO(3-2) and HST imaging of samples of SMGs, we see that the extended component is comparable in size to the low-J CO molecular gas and stellar distributions. We conclude that it is likely that the extended component seen in the stacked 870 µm maps traces a surrounding disc or envelope around the central, compact far-infrared luminous starburst.
• By stacking the 870 µm maps in bins of star-formation rate we find that the size and luminosity of the extended component is roughly constant with total star-formation rate, while the compact component becomes brighter. This suggests that the star formation taking place in the compact component is broadly decoupled from the star formation taking place in the extended component.
We have studied a large sample of SMGs using moderate resolution ALMA data and find that the morphologies observed at 0. 18 resolution are best described by bars in galaxies. However, to confirm this requires deeper observations to detect the extended components in individual maps, higher resolution imaging (∼ 0. 08) to show that the 870 µm dust continuum trace bars structures, and ideally dynamical gas measurements of the extended component to determine if this has an order rotational motion as seen for disc or more chaotic as would be expected from a merger event.