Resolving the ISM at the Peak of Cosmic Star Formation with ALMA: The Distribution of CO and Dust Continuum in z ∼ 2.5 Submillimeter Galaxies

We present ALMA Cycle 2 observations of the CO emission from four SMGs, ALESS49.1, ALESS57.1, ALESS67.1, and ALESS122.1. These sources were selected from the ALESS survey (Hodge et al. 2013b; Karim et al. 2013): an ALMA Cycle 0 follow-up program of 126 sources detected in the single-dish LABOCA Extended Chandra Deep Field South (ECDFS) Submm Survey (LESS, Weiß et al. 2009). These four sources were selected as they are spectroscopically confirmed at redshifts 2.1 < z < 2.9 as the result of an extensive spectroscopic follow-up campaign (Danielson et al. 2017). Three of the four sources in our sample were detected (ALESS57.1 and ALESS67.1) or marginally detected (ALESS122.1, through stacking) in the X-rays (Wang et al. 2013), using 4 Ms Chandra observations of the CDFS region (Xue et al. 2011) and 250 ks observations in the ECDFS (Lehmer et al. 2005). ALESS57.1 and ALESS122.1 are reported as AGN candidates (Wang et al. 2013; Danielson et al. 2017), while the origin of the X-ray emission in ALESS67.1 is most probably related to SF.

ALMA observations were taken in Cycle 2 as part of the project 2013.1.00470.S (PI: Hodge), with a total integration time of ∼2.5 hr in ALMA Band 3, covering the spectral range expected for the line emission of the CO(3–2) transition at these redshifts. The sources were observed on 2015 September 4, 6, and 20, using the 12 m array and under good phase stability/weather conditions, with PWV at zenith ranging between 1.56 and 3.03 mm. The antenna configuration consisted of 33, 36, and 35 antennas, respectively, achieving synthesized beam sizes that range between 034 and 067 major axis FWHM with the largest angular scales (LAS) between 19 and 29. The observations were calibrated based on Jupiter as the flux-calibrator, J0334–4008 as the band-pass calibrator, and J0334–401 as the phase calibrator.

New ALMA observations of ALESS49.1 (Wardlow et al. 2018) were taken during Cycle 4 on 2016 November 12, 16, and 20 as part of the project 2016.1.00754.S (PI: Wardlow). These observations were carried out using a total integration time of ∼2700 s and using the longest baselines of ∼650 m. With an angular resolution of 11, the additional data increase the signal-to-noise ratio (S/N) of our high-resolution data. We concatenated both Cycle 2 and Cycle 4 data sets achieving a combination of high-resolution imaging, high S/N and at the same time reducing concerns whether we might have resolved any significant flux from ALESS49.1, due to the short baselines covered. An analysis of the environment of ALESS49.1 is presented in Wardlow et al. (2018).

The ALMA data were calibrated following the ALMA pipeline and using the Common Astronomy Software Application package (CASA, McMullin et al. 2007). Manual flagging of a few corrupted time windows for ALESS57.1 was required after an inspection of the calibration output. The imaging process was carried out using CASA tasks (version 4.7.0). The UV data were Fourier transformed and deconvolved from the point-spread function using the CASA clean algorithm. After resampling of the visibilities at different spectral resolutions to optimize the S/N, we produced the final data cubes averaging the visibilities to channel widths of 16, 61, 82, and 77 MHz for ALESS49.1, 57.1, 67.1, 122.1, respectively. At the native resolution, the rms values achieved for the final data cubes range from 0.11 to 0.18 mJy beam⁻¹ (see Table 1).

Due to the high spatial resolution of the data, it is not trivial to estimate the masks on which to apply the clean task. We adopted an iterative cleaning technique (e.g., Chen et al. 2017), in order to optimize the mask size estimation to include possible extended low surface brightness emission. Iterative cleaning consists of drawing concentric circular mask regions at increasing radii and applying the clean task and line flux extraction within them. Plotting the resulting line fluxes against the corresponding circular region radii produces a curve of growth (see upper right panels of Figure 1). The expected behavior of the measured flux density in a curve of growth is to continuously increase as a function of radius, reaching a point of convergence at the maximum extent of the source. We used the masks inferred from the iterative cleaning method to extract the final line cubes. We explored the data for emission at different spatial scales and surface brightness, first at the native resolution using natural weighting, then using a Briggs weighting with a robustness parameter 0.5, to image the CO emission at lower S/N but slightly higher spatial resolution and finally, by tapering the visibilities to lower resolutions (1.5–2 × native beam size) to recover extended emission.

Zoom In Zoom Out Reset image size

We detect significant CO(3–2) line emission for the four sources in our sample even from the dirty data cubes, i.e., in the Fourier transformed visibilities prior to the beam deconvolution (cleaning process). The CO detections are coincident with the expected ALMA Cycle 0 positions, and at frequencies which confirm the previously reported spectroscopic redshifts (Table 1). We identify strong line detections for all the sources in the final line cubes constructed by applying the circular cleaning masks described in Section 2.1. Figure 1 shows the data as well as Gaussian profile fits to the extracted spectra. The estimated CO(3–2) line parameters are summarized in Table 2. We next briefly describe some important features of the line emission of each source.

The CO spectrum of ALESS49.1 is clearly best described by a double-peaked profile and was thus fitted using an asymmetrical two-Gaussian model. The red and blue component line widths are 280 ± 40 km s⁻¹ and 390 ± 30 km s⁻¹ (FWHM), respectively, and are separated by 300 ± 20 km s⁻¹.

The CO spectrum of ALESS57.1 is marginally spectroscopically resolved. Nevertheless, the line width is appropriately approximated by a Gaussian fit and the resulting values are summarized in Table 2.

The CO line of ALESS67.1 shows a very large line width (FWHM ∼ 800 km s⁻¹) and appears asymmetric although the asymmetry is not statistically significant, due to the low S/N in this source. We postulate that the CO line of ALESS67.1 is produced by two sources of emission. This is supported by the velocity averaged image (Figure 2), where two spatial components are recognizable independently in both CO (ALMA) and optical/near-infrared (HST) imaging: one main right component of around 1 arsec radius and a second extended emission in the east. We thus give preference to the description of this CO line as originating in two sources (possibly an early stage of a merger event) and focus on the emission extracted from the main component (1 arsec radius), which is shown as a black line overlaid on the total line emission in Figure 1.

Zoom In Zoom Out Reset image size

The CO spectrum of ALESS122.1 is the brightest in our sample, and its line profile is well described by a single Gaussian fit. Table 2 lists the CO(3–2) line properties and best-fit parameter values. A detailed description of the velocity field sampled by these CO lines will be provided in Section 3.2.1.

Continuum images were made for all the sources after excluding the channels contributing to the CO(3–2) line emission. Using natural weighting to achieve the highest sensitivity,¹⁷ we obtained images with σ_rms ∼ 6, 20, 18, 12 μJy beam⁻¹, as summarized in Table 1. The continuum emission is detected in only two of the four sources in our data set: ALESS49.1 and ALESS122.1. In both cases, the naturally weighted dust continuum emission at 3 mm (rest-frame ∼850 μm) is detected at 4-σ significance at sub-arsecond resolution (Table 2).

To test for consistency with the Cycle 0 continuum observations at 870 μm (Hodge et al. 2013b), we use an extrapolation following ${S}_{\nu }\propto {\nu }^{2+\beta }$, where S_ν is the measured flux density and β is the dust emissivity spectral index. Adopting β = 1.5 as a reasonable assumption for dusty galaxies (Casey et al. 2014; Swinbank et al. 2014), the values extrapolated from 870 μm are consistent with our measured S_3mm peak flux densities and upper limits, which are listed in Table 2.

High-resolution imaging reveals the detailed structure of the gas and dust continuum emission, which is key for placing constraints on the dynamical states of the sources. However, high-resolution interferometry is less sensitive to extended, low surface brightness emission, and thus care must be exercised in the estimation of the source sizes (e.g., Emonts et al. 2016; Dannerbauer et al. 2017; Ginolfi et al. 2017). Although the iterative cleaning method presented in Section 2.1 can provide a sense of the total source extent, it is prone to correlated noise effects and consequently may yield uncertain results. Moreover, iterative cleaning can only provide an (intrinsic) source size estimate when the extent of the source is greater than the synthesized beam. To determine the sizes independently of these possible beam-convolution effects intrinsic to the imaging process, we estimate the sizes directly from the UV data (lower-right panels of Figure 1).

We extract the UV data corresponding to the frequencies within the FWHM of the line from the continuum-subtracted cube, centered at the CO(3–2) observed frequency. For each source, we then average the visibilities at different UV distances and plot the resulting amplitudes against them, as shown in Figure 1, where the data points correspond to the mean amplitude of the visibilities within each bin and the errors are the corresponding standard deviation of the mean. We fit a Gaussian profile to the data, taking into account the uncertainties, and measure the half-light radius of the source, which corresponds to the radius within which half of the light of the galaxy is contained. We report the estimated half-light radii r_1/2 in Table 2 and adopt these values for our analysis. The half-light radii of our sources range between 04 and 08, which correspond to 2.6–6.9 kpc at the respective redshifts, implying total physical sizes (FWHM) of ∼10 kpc. In a study of the CO(3–2) emission of a similar sample of SMGs at z ∼ 2–3, Tacconi et al. (2008) found sizes ranging between 2 and 12 kpc. More recent high-resolution molecular gas studies have reported similar sizes for low-J (J ≤ 3) CO emission (e.g., Ivison et al. 2011; Riechers et al. 2011; Hodge et al. 2012). It is interesting to note that the available high-J CO observations for some of these sources (J = 6; Engel et al. 2010) reveal that the physical extent of the emission may decrease for higher J-transitions.

Figure 2 shows the CO(3–2) emission overlaid on the stellar emission as traced by the WFC3/IR (bands H₁₆₀, J₁₂₅, Y₁₀₅) and/or ACS imaging (I₈₁₄) from HST (e.g., Chen et al. 2015). We note that the astrometry in the HST images has been corrected based on the GAIA star positions, which are aligned to the Fifth Fundamental Catalog (FK5) (Gaia Collaboration et al. 2016). With the exception of ALESS122.1, the CO gas overlaps with the stellar distributions, although they have slight offsets (01–03, 1–3 kpc offsets). The sizes of the gas and stellar distribution in these sources are also roughly similar. Specifically, the ratio of the HST H₁₆₀ imaging (Chen et al. 2015) to CO half-light radii of these ALESS sources ranges between 1 and 1.5 (ALESS122.1 not included here as no H₁₆₀ is available). We will compare the extent of the different components in a statistical study for SMG populations in Section 4.2. ALESS122.1 displays the striking feature that the gas and stellar emission are completely offset. We must point out, however, that the stellar emission in ALESS122.1 is exclusively represented by the ACS I₈₁₄ band imaging (rest-frame 270 nm), in contrast to the other three sources which are covered by at least one WFC3/IR band. This particular case will be discussed in detail in Section 4.1.

Previous studies have shown that submillimeter continuum observations of similar galaxy populations reveal very compact rest-frame far-infrared-emitting regions, with median half-light radii of only ∼0.7–1.6 kpc (e.g., Ikarashi et al. 2015; Simpson et al. 2015; Hodge et al. 2016; Oteo et al. 2017). This is significantly more compact than the molecular gas emission in our sources. This difference in compactness is similarly observed between the dust continuum and radio emission of SMGs. Median radii of radio emitting regions originating from SMGs have been reported to be around 2.1 kpc on average (e.g., Chapman et al. 2004; Biggs & Ivison 2008; Miettinen et al. 2015, 2017). Motivated by these differences, we will explore in detail the relation between the radial profiles of the molecular gas, the dust continuum emission (at rest-frame 250 μm) and stellar emission in Section 4.2.

The masses of the CO(3–2) emitting gas in our sources can be estimated from the observed ¹²CO(J = 3–2) line luminosities L'_CO. These were calculated following Solomon & Vanden Bout (2005) as $L{{\prime} }_{\mathrm{CO}}=3.25\times {10}^{7}\,{S}_{\mathrm{CO}}{\rm{\Delta }}v\,{\nu }_{\mathrm{obs}}^{-2}\,{D}_{{\rm{L}}}^{2}{(1+z)}^{-3}\,{\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1}{\mathrm{pc}}^{2}$, where S_COΔv is the velocity integrated flux, in Jy km s⁻¹; ν_obs is the observed frequency, in GHz; and D_L is the luminosity distance in Mpc. The resulting values are listed in Table 2. Through an extrapolation from these CO masses we can calculate the total molecular content of the systems (dominated by H₂) by assuming a CO line luminosity to H₂ (+He) mass conversion factor, α_CO, and a brightness temperature ratio of CO(J = 3–2) to CO(J = 1–0). However, the combination of spatial and spectroscopic resolution of our data allows us to make an estimate of the total gas mass independently of the abovementioned assumptions, by estimating kinematic parameters (Section 3.2) and using further multiwavelength information to estimate the stellar mass contribution (Section 3.4). But to begin with we use the classical approach to estimate the molecular gas masses as a function of α_CO.

Using the measurements of ¹²CO(J = 1–0) by Huynh et al. (2017) for ALESS67.1 and ALESS122.1 (see Table 2), we calculate their brightness temperature ratio r_3/1 = L_CO(3–2)/L_CO(1–0), which yields r_3/1 = 1.01 ± 0.36 and r_3/1 = 0.77 ± 0.19, respectively. This is consistent with previous estimates for the SMG population (Ivison et al. 2011; Bothwell et al. 2013; Sharon et al. 2016). For ALESS49.1 and ALESS57.1 we estimate the ¹²CO(1–0) emission using the excitation ratio for SMGs derived by Sharon et al. (2016), r_3/1 = 0.78 ± 0.27, yielding L'_CO(1–0) = (0.5 ± 0.2) × 10¹¹ K km s⁻¹ pc², for both sources. The derived molecular gas masses as a function of α_CO (equivalent to α_CO = 1) are reported in Table 2.

In order to investigate the possible effects of the presence of an AGN on the molecular gas in our sample, we explore possible correlations between galaxy properties inferred from the CO measurements (listed in Table 2). We test for any bimodality in the distribution of properties such as the FWHM of the CO line, the gas-to-stellar mass fraction and SF efficiency (SFE = SFR/M_gas) corresponding to sources with and without an AGN (ALESS57.1 and ALESS122.1 are the AGN-host candidates in our sample, see Section 2.1). We find no clear bimodality correlated to the presence or absence of an AGN, although we note that this result is not conclusive given the very small size of our sample. However, we expect that the scales probed by our CO observations are not affected by AGN, since the CO transition probed (J = 3–2) is not expected to be a tracer of AGN excitation (e.g., Sharon et al. 2016).

The morphological and kinematic properties of our resolved CO images allow us to explore the nature of these star-forming systems, whether these are in an early stage of a merger with a chaotic structure, or if their velocity fields are described by an ordered rotating disk. It is important to note, however, that these scenarios are not mutually exclusive. Observing dynamics consistent with a rotating disk does not preclude the galaxy from being a merging system, as the gas is collisional and an ordered rotating disk can quickly reform after the final coalescence stage as has been shown observationally and theoretically (Robertson et al. 2006; Hopkins et al. 2009, 2013; Ueda et al. 2014).

While the few existing examples of CO detections at a similar resolution suggest a mixture of mergers and disk-like motion (e.g., Tacconi et al. 2008; Engel et al. 2010; Hodge et al. 2012), recent high-resolution continuum imaging of SMGs has shown that the dust continuum emission (at rest-frame 250 μm) is mostly described by compact disk profiles (e.g., Simpson et al. 2015; Hodge et al. 2016). Here we use our high-resolution CO observations to directly test the kinematics of these systems.

3.2.1. Position–Velocity (PV) Diagrams

To estimate the kinematic properties of our sources we first investigate their PV diagrams. We apply the interactive PV Diagram Creation task from the CASA viewer to the data cubes of our sources with spectral resolution corresponding to channel widths of ∼100 km s⁻¹. The lower panel of Figure 3 shows the kinematical major axis chosen for the extraction of the PV diagrams overplotted as blue arrows on the velocity averaged CO images. The upper panels show the resulting PV diagrams for the respective sources. We find velocity gradients in the PV diagrams of three of our four sources (with ALESS67.1 also showing some weak signs), suggesting that the bulk of emission in these sources is dominated by rotation.

Zoom In Zoom Out Reset image size

The velocity gradients of ALESS49.1, ALESS57.1 show double peaks at each side of the galactic center, although the detailed velocity structure is ambiguous given the modest S/N. In ALESS49.1, the double-peaked spatial structure gives rise to the double-peaked line profile observed in Figure 1. Since both sources have the CO centered on a single optical nucleus (as seen from their high-resolution HST imaging), the velocity structure suggests the presence of a disk structure or a disk-shaped merger remnant, rather than an early stage merger.

The CO emission in ALESS67.1 and its velocity structure appears more chaotic, which is probably accentuated by the combination of extended emission and the lower S/N achieved in this source. A detailed analysis of the velocity structure of this individual source has been conducted in a complementary study by Chen et al. (2017). They presented SINFONI  Hα observations of ALESS67.1, revealing that the extended CO emission (including the second component to the south-east) follows the bulk rotational motion of the rest-frame optical Hα emission line. Kinematic modeling of both lines concluded the bulk of the molecular gas in ALESS67.1 could be described as an ongoing merger, although without being able to reject a rotating disk given the errors. It is still uncertain whether ALESS67.1 is a multicomponent system in an early state of merging. This scenario has been discussed in Section 2.2 based on our CO and HST observations, and it has also been supported by the analysis of its kinematics based on ancillary data by Chen et al. (2017). For further dynamical mass calculations we will only consider the properties of the most luminous component of ALESS67.1 (Table 2) and assume that this component is well described by a rotating disk.

ALESS122.1 is the source with highest S/N in our sample and its velocity structure is consistent with disk rotation. Due to these reasons it was possible to analyze its velocity field quantitatively using a rotating disk model, as will be described in Section 3.2.2.

Although the velocity gradients may be consistent with large-scale disk rotation as a bulk motion, we emphasize that no conclusions can be drawn to exclude complex, disturbed gas motions on scales smaller than the resolution limit of our data. Based on this qualitative analysis of the PV diagrams in Figure 3, we will adopt the scenario of a rotating disk for our sources for the computation of their dynamical masses. The velocity line width (and thus an estimate of the velocity dispersion) and other morphological properties will be adopted from the values listed in Table 2.

3.2.2. Kinematic Modeling

To quantitatively describe the kinematics of the molecular gas in our sources, we model the data using the modeling package GalPak3D. GalPak3D is a Bayesian parametric Markov Chain Monte Carlo fitter for three-dimensional (3D) galaxy data that attempts to disentangle the galaxy kinematics from resolution effects (Bouché et al. 2015). Starting from a set of uniform priors on the disk parameters (source center, radius, inclination, velocity dispersion, etc.,) a 3D disk galaxy model is produced assuming an exponential radial flux profile. The fitting process consists then in comparing the observed data cube to mock data cubes modeled at each point of the parameter space. In each step the reduced ${\chi }_{\nu }^{2}$ value is then minimized to efficiently sample the probability density functions of the parameters. The reliability of the inferred kinematic parameters in GalPak3D goes approximately as $\tfrac{\delta p}{p}\propto \tfrac{{r}_{1/2}}{{r}_{\mathrm{beam}}}\times {\mathrm{SB}}_{1/2,\mathrm{obs}}$, i.e., the uncertainty in the inference of the parameter p is inversely proportional to the source size to resolution ratio ($\tfrac{{r}_{1/2}}{{r}_{\mathrm{beam}}}$) and surface brightness of the source (SB_1/2,obs). Given the dependence of these uncertainties on the quality of the data and due to the S/N limits in most of our sources, we were able to apply this method only to the source with the highest S/N and resolution, ALESS122.1. For the sources on which the dynamical modeling cannot be applied, we use the line profiles to draw estimates of kinematic properties such as rotational velocity (Section 3.3).

Figure 4 shows zeroth-moment (integrated intensity) and first-moment (intensity-weighted velocity) maps for ALESS122.1. The four images correspond to the moment maps of the observed data cube, the convolved model, deconvolved model and residual (model-substracted) cubes, where the latter three cubes are output products of GalPak3D. We estimated the noise per pixel in the zeroth-moment map as $\sqrt{N}\times {\sigma }_{\mathrm{chan}}$, where N is the number of channels used and σ_chan is the rms noise per channel. Using this value, for the computation of the first-moment maps of the four cubes, we masked out all pixels with S/N < 2. The bottom row of Figure 4 shows the first-moment maps (intensity-weighted velocity) of the corresponding cubes. The color scale has been chosen to represent the velocity width of CO line of the source.

Zoom In Zoom Out Reset image size

The resulting 3D model cube shows a modest agreement with the data, with a reduced ${\chi }_{\nu }^{2}$ of 1.87, which is calculated over the total area and frequency range covered by the observed and modeled data cubes, as represented in Figure 4. The slightly high, reduced ${\chi }_{\nu }^{2}$ value can be explained by the residual structure, which can be seen in the 0th and 1st moment maps of the residual cubes. The model parameters obtained for the source agree with comparable properties inferred directly from the line profile.

The bulk of the velocity field of ALESS122.1 is consistent with a rotation-dominated disk (mindful of the residual clumps), with an inclination of i ∼ 52°, a maximum rotational velocity of v_max ∼ 560 km s⁻¹, and a velocity dispersion of σ_v ∼ 130 km s⁻¹. The residual image in the right upper panel of Figure 4 reveals few residual clumps between 2 and 3σ significance, consistent with noise. We conclude that the bulk of the velocity field is to a first order consistent with disk rotation.

The kinematic properties of a galaxy, obtained e.g., through modeling its velocity field, can provide a reliable estimation of the mass enclosed within the region covered by the emitting medium. At high redshift, this is a complicated task, since observations of molecular gas are frequently poorly spatially resolved and the morphology of the mass distributions are thus usually unknown. Based on the kinematic study above and the size estimates from our analysis, we calculate the dynamical masses assuming the bulk of the emission in our sources can be well described by a rotating disk. The total dynamical mass within a radius r = 2r_1/2 is then given by:

$$\begin{eqnarray}&&{M}_{\mathrm{dyn}}(r\lt 2{r}_{1/2})=\displaystyle \frac{{({v}_{\max }/\sin (i))}^{2}\times 2{r}_{1/2}}{G},\end{eqnarray} \tag{ 1 }$$

where r_1/2 is the half-light radius estimated through UV fitting measured in kpc, i is the inclination of the galaxy (2), and G the gravitational constant (e.g., Solomon & Vanden Bout 2005; Erb et al. 2006; de Blok & Walter 2014). In cases where the velocity field cannot be modeled and v_max cannot be directly inferred, we can use the approximation that the velocity width is twice the maximum rotational velocity and thus we can replace v_max = 1/2 × ΔV_FWHM in Equation (1). This is a good assumption for a flat rotation curve with a steep inner rise and a filled exponential disk as shown by de Blok & Walter (2014). For ALESS122.1 we use the parameter values for v_max resulting from the kinematic modeling in Section 3.2.2, although both estimates provide consistent results within errors.

A large uncertainty in the estimation of the dynamical masses is contributed by the inclination angle of the disk. In previous dynamical mass studies, it has been customary to adopt a value of 〈sin ²(i)〉 = 2/3 (corresponding to i = 54.7 deg), which is the average value expected for a randomly oriented population of disk galaxies. However, given our high-resolution data, we have enough information to use a simple assumption motivated by the early observation that the apparent axis ratio is closely related to the inclination angle for a disk. We use the relation ${\cos }^{2}(i)\,={{((b/a)}^{2}-{q}_{0}^{2})(1-{q}_{0}^{2})}^{-1}$ (Hubble 1926), where (b/a) is the axis ratio and q₀ is the inherent thickness of the disk, assuming q₀ = 0.1 (Nedyalkov 1993). Obtaining the deconvolved axis ratios using the CASA routine imfit, we compute the inclination angles for ALESS49.1, ALESS57.1, ALESS67.1, and ALESS122.1 to be i ∼ 80° ± 30°, 40° ± 50°, 50° ± 20°, and 40° ± 20°, respectively. For ALESS122.1, where a more complex analysis was possible in Section 3.2, we will adopt the inclination angle estimated through the kinematic modeling, i ∼ 52° ± 2°, which is in agreement with the axis ratio approximation within the errors. Although this is an approximation and the uncertainties are large, this estimation is superior to assuming the single value of 〈sin ²(i)〉 = 2/3 for all systems, since the inclination in individual sources may anti-correlate with the line width. It is interesting to note the high inclination angle of ALESS49.1, which is also supported by the small CO size and double-peaked line profile. Although these profiles are not uncommon for rotating disks and have been found in at least 40% of SMGs (Tacconi et al. 2008; Bothwell et al. 2013), the asymmetry between the double-peaked profile suggests a non-uniform distribution of the gas in the "disk," produced possibly by minor instabilities or an unresolved merger of two gas disks.

The resulting dynamical masses for our sample range from 1 to 5 × 10¹¹ ${M}_{\odot }$ (Table 2), as calculated within 2× the half-light radii of the sources r_1/2. These values are in general agreement with previous measurements of the dynamical masses of two of these sources based on other tracers (CO(J = 1–0) and Hα, Chen et al. 2017; Huynh et al. 2017).¹⁸ This range corresponds to masses at the high end of the average values found for other SMG samples by Ivison et al. (2010) $(2.3\times {10}^{11}\,{M}_{\odot })$, Tacconi et al. (2008) (1.3 × 10¹¹ ${M}_{\odot }$), Engel et al. (2010) (∼1.9 × 10¹¹ ${M}_{\odot }$) and slightly below the values found for some extreme sources such as GN20 (5.4 ± 2.4 × 10¹¹ ${M}_{\odot }$, Hodge et al. 2012) and SMMJ131201 (9.5 ± 2.4 × 10¹¹ ${M}_{\odot }$, Engel et al. 2010). The uncertainties in these values are propagated from all of the parameters used in their calculation, including the inclination angles i (Table 2).

The dynamical mass estimates discussed in Section 3.3 can be related to the various mass components in galaxies following:

$$\begin{eqnarray}&&{M}_{\mathrm{dyn}}(r\leqslant 2{r}_{1/2})={M}_{\mathrm{baryons}}(r\leqslant 2{r}_{1/2})+{M}_{\mathrm{DM}}(r\leqslant 2{r}_{1/2}),\end{eqnarray} \tag{ 2 }$$

where M_baryons = M_gas + M_* and M_DM is the DM contribution. Thus, given assumptions on the DM content and stellar masses, the dynamical mass estimates can be used as an independent method to constrain the total gas masses in galaxies. Before we proceed to constrain the gas masses using this method, we summarize the unknown parameters intrinsic to this calculation. We note that given the scarcity of CO(J ≤ 3) gas data at high resolution, most of the existing literature on molecular gas in high-redshift galaxies usually assume fixed values for unknown parameters, estimating gas masses which are thus sensitive to these assumptions and lacking information on the inherent systematic uncertainties.

To begin with, the H₂ molecules that comprise the bulk of the molecular gas reservoir in galaxies are not directly observable. Measurements of molecular gas thus rely on CO observations, via a conversion from the ground state CO(1–0) luminosity which is parametrized with the factor α_CO. Indeed, high-redshift measurements are providing increasing evidence that the α_CO values in the early universe can be lower than in solar-metallicity galactic disks (α_CO ∼ 0.8–1.0,e.g., Tacconi et al. 2008; Hodge et al. 2012; Bothwell et al. 2013).

In addition to the uncertainties in the total gas mass estimation, stellar mass estimates are typically obtained via SED fitting, which relies on assumptions regarding the star formation history (SFH) of the stellar populations and dust attenuation, and can be uncertain especially for starburst systems (such as SMGs). In particular, it has been increasingly observed that SED fitting of starburst systems suffers from a degeneracy between the SFH, the mass-to-light ratio (M_*/L_H) and the age of the galaxy (Hainline et al. 2011; Michałowski et al. 2012; Simpson et al. 2014). This degeneracy means that in a galaxy with an instantaneous burst SFH the total stellar mass would be contained almost exclusively within young luminous stellar populations, resulting in a low (M_*/L_H) ∼ 0.05. Conversely, in a galaxy of intermediate age with a constant SFH, most of the stellar mass would be distributed in old faint populations whose emission is outshone by the luminous newly born stars, implying a higher mass-to-light ratio (M_*/L_H ∼ 0.3). As shown by recent resolved SED fitting studies (Sorba & Sawicki 2018), this effect could underestimate total stellar masses by factors of up to 5, especially in sources of high specific SF rates such as SMGs.

Rather than adopting a single (M_*/L_H) or α_CO, we can instead parametrize the baryonic mass as:

$$\begin{eqnarray}&&{M}_{\mathrm{baryons}}={L}_{H}\times \,{M}_{* }/{L}_{H}+{\alpha }_{\mathrm{CO}}\times {L}_{\mathrm{CO}},\end{eqnarray} \tag{ 3 }$$

where L_H is the rest-frame H-band luminosity corrected for dust obscuration (here estimated by da Cunha et al. 2015) and L_CO is the CO(J = 1–0) luminosity estimated in this work.

Finally, the dynamical mass of a galaxy also includes the mass of the non baryonic DM component, which is a source of uncertainty, as no independent measurement of this mass fraction is available for our SMGs. Unlike observations in local disk galaxies, which have reported DM fractions (f_DM) of ∼50% (e.g., Courteau & Dutton 2015), it has been claimed by recent spectroscopic surveys (e.g., Price et al. 2016; Wuyts et al. 2016; Genzel et al. 2017) that more compact star-forming disk galaxies at z > 2 appear to be heavily baryon dominated with f_DM ∼ 10%–20%, although these calculations involve the same unknown factors as those that are used for gas fractions. In contrast, recent simulations (e.g., Lovell et al. 2018) have reported that the DM fractions in disc-like galaxies could range up to 65% for galaxies with stellar masses similar to SMGs (∼10¹¹ M_⊙).

Combining Equations (2) and (3), the total mass can be expressed as:

$$\begin{eqnarray}&&{M}_{\mathrm{dyn}}=\displaystyle \frac{{L}_{H}\times ({M}_{* }/{L}_{H})+{L}_{\mathrm{CO}}\times {\alpha }_{\mathrm{CO}}}{(1-{f}_{\mathrm{DM}})},\end{eqnarray} \tag{ 4 }$$

where the conversion factor α_CO, the stellar mass-to-light ratio (M_*/L_H), and the DM fraction f_DM are the unknown parameters. Given the small number of independent data points available (one for each of our four sources), in order to reduce the parameter space we will keep the f_DM value fixed, and discuss the effect of low and high f_DM values to the other parameters distributions later on.

First, we explore the parameter space built up in Equation (4) by applying an MCMC technique, using the open-source algorithm emcee (Foreman-Mackey et al. 2013) for the sampling. Based on the likelihood of the measured dynamical masses M_dyn (estimated in Section 3.3 as a function of L_CO and L_H) given in this model, we sample the posterior probability density function (posterior PDF) for α_CO and (M_*/L_H). It is important to note that at this point the likelihood function accounts only for the uncertainties in the dynamical masses M_dyn. To account for the uncertainties in the L_CO and L_H variables as well, we use a Monte Carlo approach and repeat the inference exercise presented above for 100 iterations, each time using values randomly chosen from within a normal distribution that corresponds to the mean and standard deviation of our L_CO and L_H measurements.¹⁹ Finally, we average the posterior PDFs gathered in these 100 iterations and obtain a final probability distribution for the parameters that takes into account the uncertainties in the measurements of all observables (M_dyn, L_CO and L_H).

Figure 5 shows the one- and two-dimensional final posterior PDFs of the (M_*/L_H) and α_CO parameters adopting three different DM fractions f_DM = 0%, 15%, and 30%. The covariance between the parameters can be recognized in the two-dimensional PDF (Figure 5) as an elongated ellipse in the central contour level plot. Higher DM fractions decrease the value of α_CO and consequently M_gas. This trend is similarly followed by the (M_*/L_H) parameter. For a DM fraction of 15% we find a median value of ${\alpha }_{\mathrm{CO}}={1.1}_{-0.7}^{+0.8}[{M}_{\odot }/{({\rm{K}}\mathrm{km}{{\rm{s}}}^{-1}{\mathrm{pc}}^{2})}^{-1}]$, while for a larger DM fraction of 30%, these values decrease to ${\alpha }_{\mathrm{CO}}\,={0.9}_{-0.6}^{+0.7}[{M}_{\odot }/{({\rm{K}}\mathrm{km}{{\rm{s}}}^{-1}{\mathrm{pc}}^{2})}^{-1}]$. For f_DM = 0 we obtain an upper limit to ${\alpha }_{\mathrm{CO}}\leqslant {1.4}_{-0.9}^{+0.9}[{M}_{\odot }/{({\rm{K}}\mathrm{km}{{\rm{s}}}^{-1}{{pc}}^{2})}^{-1}]$. For the mass-to-light ratio parameter we find $({M}_{* }/{L}_{{\rm{H}}})={0.22}_{-0.15}^{+0.16}\,{M}_{\odot }/{L}_{\odot }$, $({M}_{* }/{L}_{{\rm{H}}})={0.19}_{-0.13}^{+0.15}\,{M}_{\odot }/{L}_{\odot }$, and $({M}_{* }/{L}_{{\rm{H}}})={0.16}_{-0.10}^{+0.12}\,{M}_{\odot }/{L}_{\odot }$ for DM fractions of 0%, 15%, and 30%. Using these values, the gas reservoirs in our sources have an average gas mass of ∼0.9 × 10¹¹ ${M}_{\odot }$ (DM contribution of 15%) and ∼0.7 × 10¹¹ ${M}_{\odot }$ (DM contribution of 30%), corresponding to average gas-to-total-mass fractions of 0.26% and 0.33%, respectively. These sources show similar gas fractions than those reported for typical star-forming galaxies (≤50%, Tacconi et al. 2017). For comparison, Bothwell et al. (2013) presented a study of a representative sample of SMGs at z ∼ 2 from predominantly CO(3–2) observations, finding mean gas masses of (3.2 ± 2.1) × 10¹⁰ ${M}_{\odot }$ assuming α_CO = 1. Since the resulting PDF for the M_*/L_H and age parameters for our galaxies is broad and poorly constrained, the galaxy stellar population properties inferred through this method (∼0.1–2.0 Gyr) easily agree with those estimated through SED fitting by da Cunha et al. (2015), which are also associated with large errors (∼70%). They find the age of these galaxies range between ∼0.04 and 1.22 Gyr, corresponding to M_*/L_H between 0.16 and 0.2.

Zoom In Zoom Out Reset image size

As discussed above, the values recovered for the parameters are subject to large uncertainties due to our small sample size. However, we note that one goal of this exercise is to shed light upon the degeneracies among the key parameters and large uncertainties inherent to gas mass calculations. Most importantly, given the increasing availability of dynamical mass measurements for high-redshift galaxies with ALMA, we consider that this approach has potential for constraining the molecular gas properties with more accuracy, especially since larger sample sizes will allow the exploration of a more detailed parameter space in the future.

An interesting case among our samples is the apparently uncorrelated distribution of dust continuum, gas, and stars in ALESS122.1 (Figure 6). The gas and stellar emission extend across regions of similar size, separated by ∼05 (∼5 kpc at z = 2) as measured from their centroids, but are tightly aligned next to each other. There is almost no overlap of gas and stellar emission, even when considering the combined astrometric uncertainties²⁰ (<01). Moreover, while the gas and stars have half-light radii r_1/2 ∼ 4 kpc and are thus extended on scales of ∼10 kpc, the detected dust continuum emission (3 mm observed-frame, corresponding to ∼850 μm rest-frame) is confined to a central region of ∼5 kpc, with the emission peak spatially coincident with that of the molecular gas emission. Similar physical offsets between the dust continuum, gas, and stellar emission, and specifically between ALMA and HST, have been previously found in a number of high-redshift sources (e.g., Riechers et al. 2010; Chen et al. 2015; Hodge et al. 2015; Elbaz et al. 2017; Fujimoto et al. 2017; Simpson et al. 2017). For example, based on low-resolution 870 μm continuum data (∼15), Chen et al. (2015) reported a statistical offset with the existing stellar component traced by HST as large as Δp = 04 (∼4 kpc at z = 2) in 66% of the ALESS SMGs. Similarly, Hodge et al. (2015) found in the z = 4 galaxy, GN20, that the FIR and CO emission are offset by 06 (4 kpc) from the peak of the rest-frame UV emission as traced by the HST/WFC3 F105W image.

Zoom In Zoom Out Reset image size

We consider three scenarios for explaining the mismatch of the distributions of stellar and gas/dust continuum emission in ALESS122.1. The first and most plausible scenario assumes that the optical component extends to the regions of gas and dust continuum emission, while suffering extreme extinction in these regions. Such high-extinction scenarios have been suggested previously in starbursts with extreme SFRs and in cases of partial or total absence of optical counterparts (e.g., Hodge et al. 2012; Walter et al. 2012; Simpson et al. 2015). This is not uncommon in submm selected samples, where optical detection rates are only ∼70%–80% even with deep data (e.g., Simpson et al. 2014).

The second scenario implies that the observed offset reflects a physical misalignment between the gas/dust continuum distributions and that of the stellar mass. This argument has been tested by Chen et al. (2015) for the full ALESS sample. They tested this hypothesis by comparing the positional offsets found in low-redshift (z < 2) to those in high redshift (z > 2) subsamples, taking into account that at lower redshift the optical photometry would probe wavelengths less affected by obscuration. They found no statistical difference in the measured positional scatter between the low- and high-redshift samples, implying that obscuration is not the dominant factor for the population as a whole.

The third scenario implies that ALESS122.1 is an ongoing, merging system comprised of two components, one of which is a heavily obscured SMG that has no detected optical counterpart, and the other is an optically bright galaxy in the HST-ACS F814W imaging (which corresponds to rest-frame ∼3000 Å). High-resolution optical-MIR counterparts at other wavelengths are needed in order to confirm the nature of this source and explain the uncorrelated distributions of its physical components.

Regardless of the underlying scenario, the significant anticorrelation between the dust continuum and optical imaging disprove some assumptions that are intrinsic to multiwavelength studies such as energy balance in SED fitting (this has been previously discussed by e.g., Hainline et al. 2011; Chen et al. 2015, 2017; Hodge et al. 2016; Simpson et al. 2017).

To gain a general understanding of the distributions of the molecular gas, dust continuum, and stellar emission in SMGs and how they relate to each other, next we conduct a stacking analysis using the CO emission measured in our sources and the available multiwavelength data at similar resolution.

We first produced the stacked radial profiles of the CO emission of the four ALESS sources presented here. To achieve this, we create 10 × 10'' thumbnails of our sources, align their centers and scale them to have the same peak flux. We then produce the radial profiles of the individual sources by radially averaging the emission centered at the source peak, and finally, we average the individual radial profiles in our sample (Figure 7).²¹ To account for the effect of outlier structures of the individual sources, we estimate the uncertainties of the radial profiles through bootstrapping, e.g., by creating 1000 different average profiles, where each iteration used a randomly resampled combination of sources.

Zoom In Zoom Out Reset image size

Using the same technique, we also stack high-resolution dust continuum emission (016 FWHM at rest-frame 250 μm) of 16 luminous ALESS sources studied by Hodge et al. (2016). This sample of sources is representative of the ALESS survey as a whole (ALESS median z = 2.5 ± 0.2 and L_IR = (3.0 ± 0.3) × 10¹² L_⊙, Simpson et al. 2014; Swinbank et al. 2014) as it has a median redshift of z = 2.6 ± 0.5 and infrared luminosity of ${L}_{\mathrm{IR}}=(3.6\pm 0.9)\times {10}^{12}\,{L}_{\odot }$. Most importantly, these properties are also comparable to those of the sample used for the CO stacked profiles (median redshift of z = 2.5 ± 0.4 and infrared luminosity of ${L}_{\mathrm{IR}}=(5.6\pm 2.2)\,\times {10}^{12}\,{L}_{\odot }$), although both samples have only one source in common (ALESS67.1).

Finally, we stack high-resolution HST imaging (017 in the H₁₆₀ band) available for the same 16 ALESS sources following the procedure described above. The stacked surface brightness profiles for the gas, dust continuum, and stellar components are shown in Figure 7 as green, red, and blue data points, respectively.

To quantify the extent and distribution of the molecular, submm and optical/near-infrared emission, first we assume that the observed emission can be described by an exponential profile (Sérsic profile with index n = 1) convolved with a Gaussian profile, which corrects for the effect of the shape of the beam. Based on this we fit a one-dimensional exponential to the CO and dust continuum profiles and find that they are well described by such a model. We note that for fitting the dust profile we discarded the data points at radius of >7 kpc, since these are evidently affected by artifacts produced by the side lobes of the beam (shown as the shaded area in Figure 7).²² The radial profile of the stellar emission results in a bad quality fit when using an exponential, i.e., when fixing the Sérsic index n = 1. To improve this, we leave the Sérsic index as a free parameter and achieve a better quality fit with n = 0.8 ± 0.2, which is consistent within the errors with the median value found for the ALESS sample as whole (n = 1.2 ± 0.3, Chen et al. 2015). The best-fit models consist then of exponential profiles with radii of r = 1.7 ± 0.1, 3.8 ± 0.1, and 4.0 ± 2.0 kpc for the dust continuum, CO, and stellar emission, respectively.

These results reveal that both the stellar (H₁₆₀) and molecular gas (CO(J = 3–2)) emission have similar sizes, but are >2 × more extended than the dust continuum emission, measured at rest-frame 250 μm. In addition, we note that the estimated size of the CO(3–2) emission should be considered a lower limit for the extent of the total molecular-mass distribution, as the ground state transition CO(1–0) has been observed to be even more extended than the higher J-transitions (Engel et al. 2010).

A crucial consequence of the different scales probed by the molecular gas, dust, and stellar tracers in these galaxies (as shown in Section 4.2) is the uncertainty in using conversion factors between the luminosities of individual tracers. The compactness of the dust continuum emission in SMGs has been discussed in the literature (Ikarashi et al. 2015; Simpson et al. 2015; Barro et al. 2016; Hodge et al. 2016; Simpson et al. 2017; Tadaki et al. 2017a), as ALMA continuum observations at high angular resolution (013–04) have revealed median radii of ∼0.7–1.5 kpc.

Apart from CO and optical imaging, high-resolution observations of other SF tracers such as synchrotron and free–free emission have allowed a comparison of relative sizes of dust continuum and gas emission regions. Indeed, radio interferometric observations of SMGs with Very Large Array at 1.4 GHz (Chapman et al. 2004; Biggs & Ivison 2008), at 3 GHz (Miettinen et al. 2015, 2017), and at 6 GHz have shown that dust continuum sizes appear to be 1.4–4.4 times smaller than the ratio.

What is the origin of the significantly smaller sizes (>2) of the dust continuum emission in SMGs compared to the other physical components? There are three parameters that could affect the emission of dust at larger radii: the dust-to-gas ratio, the dust temperature, and the optical depth of the dust and gas. The small extent of the dust continuum emission at a rest-frame 250 μm is likely to trace the heating of dust by SF or AGN, as it is confined to central compact regions, as expected for highly star-forming systems. A decrease in optical depth at larger radii would similarly affect the observed submm continuum emission, as would a radially dependent dust-to-gas ratio. Studies in the nearby universe, however, (e.g., Sandstrom et al. 2013; Groves et al. 2015) suggest that no large variations in the latter are expected across the galaxy, making this scenario less plausible.

The different spacial extent for the rest-frame 250 μm dust continuum emission and the CO line emission can be explained in the context of self-consistent radiative transfer of the CO(3–2) and dust continuum emission. Weiß et al. (2007) presented such a model, where the line and dust continuum emission are linked through a constant dust-to-gas mass ratio throughout the galaxy, specifically, by connecting the gas column density ${N}_{{{\rm{H}}}_{2}}$ to the optical depth of the dust, ${\tau }_{\mathrm{dust}}(r)\sim {N}_{\mathrm{dust}}\sim {N}_{{{\rm{H}}}_{2}}$. The different dust continuum and CO sizes can then be explained by assuming a radially decreasing gas column density and a radially decreasing temperature distribution T(r). Following the radiative transfer equation, and neglecting the background temperature, we have

$$\begin{eqnarray}&&S(r)\propto T(r)\cdot (1-{e}^{-\tau }),\end{eqnarray} \tag{ 5 }$$

where T(r) implies T(r) = T_dust(r) for the dust and T(r) = T_ex(r) for the molecular gas. We can expect that the long wavelength dust continuum emission, which is optically thin (τ ≪ 0), will be sensitive to both the column density and the temperature gradients (S_cont(r) ∝ T_dust(r) · τ_dust(r)), whereas the optically thick CO line emission (τ ≫ 1) will to first order only respond to the temperature gradient (${S}_{\mathrm{CO}}(r)\propto {T}_{\mathrm{ex}}(r)\sim {T}_{\mathrm{kin}}(r)$ for the low-J CO transitions, which are close to thermalization. This would have the effect that the dust continuum emission becomes fainter more rapidly, as a function of radius, than the gas emission.

Based on these arguments, we have used an extended version of the dust and line radiative transfer models of Weiß et al. (2007) to demonstrate this effect for the radial distribution of the stacked CO(3–2) and rest-frame 250 μm dust continuum data shown in Figure 7. In the model we use an exponentially decreasing column density distribution (Sérsic index of n = 1) and a linear temperature gradient from the center to the outer parts of the disk, and we solve the dust and line radiative transfer in several radial bins (more details on the model will be presented in A. Weiß et al. 2018, in preparation). Figure 7 shows the fit of an exponentially decreasing gas component, taking into account the different beam sizes for the dust continuum and CO line observations. Given the uncertainties of the data, the model can reproduce the different radial behavior of the dust continuum and the CO(3–2) line emission. From this fit we derive an intrinsic exponential scale length of 3.4 kpc (04) for the underlying H₂ distribution, intermediate between the intrinsic sizes we have estimated earlier based on the deconvolved intensity profiles for the dust and the CO(3–2) line emission.

In this section we have proposed that assuming a unique intrinsic distribution of dust and gas, temperature, and optical depth gradients alone could give rise to the apparent size differences observed between the CO and dust continuum emission. If these parameters are the origin of the compactness of dust continuum emission observed in the high-redshift universe (here probed at a rest-frame ∼250 μm), caution must be exercised when extrapolating properties of high-resolution continuum observations to conclusions on the molecular phase of the ISM. Specifically, these results advise against combining unresolved molecular gas observations with radius estimates from resolved dust continuum observations when calculating key parameters such as dynamical masses. In order to empirically test the assumptions placed into this model on the temperature and optical depth gradients, high-resolution observations of the dust continuum at different frequencies would be necessary, as a gradient in observed sizes of the emission regions is a direct tracer of the dust optical depth across the SED.