Starburst to Quiescent from HST/ALMA: Stars and Dust Unveil Minor Mergers in Submillimeter Galaxies at z ∼ 4.5

We selected a sample of six of the highest-redshift unlensed SMGs from Toft et al. (2014; see Table 1), which are part of the extensive (sub)millimeter interferometric and optical/millimeter spectroscopic follow-up campaigns in the COSMOS field (Scoville et al. 2007; Younger et al. 2007, 2008; Capak et al. 2008, 2011; Schinnerer et al. 2008; Riechers et al. 2010, 2014a; Smolčić et al. 2011, 2015; Yun et al. 2015). All of our sample sources had been spectroscopically confirmed to be at 4.3 ≲ z ≲ 4.8, except AzTEC5 at a slightly lower (photometric) redshift (see Table 3). We refer the reader to Smolčić et al. (2015) for a detailed description of the selection of each source.

HST WFC3/IR observations of AzTEC1, J1000+0234, and Vd-17871 were taken in the F125W and F160W bands at a two-orbit depth on each filter (program 13294; PI: A. Karim). For AK03 and AzTEC5, WFC3 F125W and F160W imaging were taken from the CANDELS survey (Grogin et al. 2011; Koekemoer et al. 2011). AzTEC/C159 was not in our HST program due to its faintness at near-IR wavelengths. Additionally, we included COSMOS HST ACS/WFC F814W images (Koekemoer et al. 2007) available for the full sample. At the redshift of the sources, these three bands probe the UV continuum regime in the range ∼140–300 nm (175–345 nm for AzTEC5).

In order to process the HST observations from our program, we made use of the DrizzlePac 2.0 package (Gonzaga et al. 2012). First, we assured a good alignment between the four dithered frames on each band using the TweakReg task. Next, we combined the frames with AstroDrizzle, employing the same parameters as used in the CANDELS reduction procedure: final_scale = 0.06 and final_pixfrac = 0.8 (Koekemoer et al. 2011).

For the purpose of this work, it is important that all three bands are properly aligned, sharing a common World Coordinate System (WCS) frame with accurate absolute astrometry. In order to guarantee the absolute astrometric accuracy, we chose the COSMOS ACS F814W image as the reference frame. The fundamental astrometric frame for COSMOS uses the CFHT Megacam i-band image (Capak et al. 2007). The latter is tied to the USNO-B1.0 system (Monet et al. 2003), which is also tied to the VLA 1.4 GHz image (Schinnerer et al. 2004), ensuring an absolute astrometric accuracy of 005–01 or better, corresponding to ∼1–1.5 pixel for our pixel scale. To align the F125W and F160W images to the F814W WCS, we used TweakReg along with SExtractor (Bertin & Arnouts 1996) catalogs of the three bands, with the F814W catalog and frame as references. Once the three bands shared the same WCS frame, we propagated the WCS solution back to the original flt.fits frames using the TweakBack task and then ran AstroDrizzle once again to produce the final drizzled images. In the case of AK03 and AzTEC5, where the F125W and F160W data came from CANDELS, this alignment procedure is not necessary, since the images are already matched to the COSMOS WCS. The final drizzled images in the three bands were resampled to a common grid and a pixel scale of 006 pixel⁻¹ using SWarp (Bertin 2010).

Our galaxies were observed in ALMA's Cycle 2 as part of Cycle 1 (program 2012.1.00978.S; PI: A. Karim). We used the ALMA band 7 and tuned the correlator such that a single spectral window (SpW) would cover the [C ii] line emission of our galaxies, while three adjacent SpWs with a total bandwidth of 5.7 GHz would be used for continuum detection. These continuum SpWs are those analyzed in our study, while the [C ii] line data cubes are presented in A. Karim et al. (2018, in preparation).

All observations were taken in 2014 June using 34 12 m antennas in configuration C34-4 with a maximum baseline of ∼650 m. For all galaxies, J1058+0133 and J1008+0621 were used as bandpass and phase calibrators, respectively. In contrast, the flux calibrator is not the same for all galaxies, varying among Titan, J1058+0133, Ceres, and Pallas. Calibration was performed with the Common Astronomy Software Applications (CASA; version 4.2.2) using the scripts provided by the ALMA project. Calibrated visibilities were systematically inspected, and additional flaggings were added to the original calibration scripts. Flux calibrations were validated by checking the flux density accuracies of our phase and bandpass calibrators. Continuum images were created by combining the three adjacent continuum SpWs with the CASA task CLEAN in multifrequency synthesis imaging mode and using a standard Briggs weighting scheme with a robust parameter of −1.0. The effective observing frequencies, synthesized beams, and resulting noise of these continuum images are listed in Table 2.

Each galaxy yields a significant continuum detection signal-to-noise ratio (S/N) > 10 at the phase center of our images. Their fluxes were measured via 2D Gaussian fits using the python package PyBDSF and are given in Table 2. These fluxes are consistent with the measured 850 μm fluxes from the S2COSMOS/SCUBA2 survey (J. M. Simpson et al. 2018, in preparation). This suggests that there is no extended emission that is resolved out in the higher-resolution ALMA observations.

In terms of the WCS, we do not expect a significant offset in the ALMA absolute astrometry with respect to the COSMOS WCS. The main source of uncertainty for the relative astrometry between ALMA and HST is the uncertainty in the HST absolute astrometry with respect to the COSMOS WCS, which is $\lt 0\buildrel{\prime\prime}\over{.} 1$, as shown in the previous section. Schreiber et al. (2017) tested the relative astrometry between an ALMA single pointing and an HST image tied to the COSMOS WCS. Following Schreiber et al. (2017), at our S/N and resolution, the combined pointing accuracy between our ALMA and HST images is <012, corresponding to <2 pixels for our pixel scale.

The HST data span three different bands from two different instruments; consequently, the spatial resolution is different. It is essential to compare the same physical regions when obtaining resolved color information. We therefore degraded the ACS F814W and WFC3 F125W images to the resolution of the WFC3 F160W data (018 FWHM), which has the broadest point-spread function (PSF). First, we created a stacked PSF in the different bands, selecting stars that were not saturated and did not show irregularities on their light profiles. Second, we derived the kernels to match the ACS F814W and WFC3 F125W PSFs to the PSF in the WFC3 F160W image using the task PSFMATCH in IRAF. We applied a cosine bell function tapered in frequency space to avoid introducing artifacts in the resulting kernel from the highest frequencies. To get the best size for the convolution box, we iterated over different values. Finally, we implemented the kernel on the ACS F814W and WFC3 F125W images. The matched PSF FWHMs in the different bands deviate by less than 2%.

ALMA continuum images also show different spatial resolution compared to that in the PSF-matched HST images (median synthesized beam size of 030 × 027 versus 018 FWHM, respectively). It is important to perform the measurements in the same physical regions when comparing HST and ALMA photometry as well, such as to derive rest-frame FIR/UV ratios. When this is required, we used HST images matched to the resolution of the ALMA continuum images constructed following the procedure explained above. In this case, the kernel was computed from the WFC3 F160W PSF and the ALMA clean beam and then applied to the PSF-matched HST images. The matched PSF FWHMs in the HST and ALMA images deviate by less than 2%.

We applied a smoothing technique to the HST images to enhance the S/N and improve our ability to detect low surface brightness features and color gradients between neighboring pixels.

The code employed for this purpose was ADAPTSMOOTH (Zibetti 2009), which smooths the images in an adaptive fashion, meaning that at any pixel, only the minimum smoothing length to reach the S/N requested is applied. In this way, the images retain the original resolution in regions where the S/N is high, and only low-S/N regions are smoothed.

We required a minimum S/N = 5 and a maximum smoothing length of two neighboring pixels in the code. The former holds true for uncorrelated noise, which is not the case for the drizzled HST images analyzed here. In our images, the chosen value of 5 corresponds to S/N ∼ 2 when taking into account the noise correlation and pixelation effects in the code. The chosen smoothing length prevents cross-talking between pixels, also reduced by calculating the median of the pixel distribution inside the smoothing radius as opposed to the mean. Such a smoothing length was chosen to match the resolution in the HST data, so the smoothing technique does not smear out the images.

We generated a smoothing mask for each band, which is a mask of the required smoothing length to reach the requested minimum S/N for each pixel. When applying a mask to the images, the pixels that do not reach the minimum S/N level are blanked out by the code. If a pixel reached the minimum S/N in at least two bands, we replaced the smoothing length in the mask by the maximum value of them. This guarantees that the same physical regions are probed in different bands, at the same time maintaining the signal if a pixel is above the minimum S/N in only one band.

A series of additional multiwavelength-imaging data sets in the optical/IR were employed in this work: the Subaru Hyper-Suprime-Cam (HSC) from the HSC Subaru Strategic Program (SSP) team and the University of Hawaii (UH) joint data set in the g, r, i, z, and y bands (Tanaka et al. 2017), with spatial resolutions (seeing FWHM) of 092, 057, 063, 064, and 081, respectively; the UltraVISTA DR3 survey (McCracken et al. 2012) covering the near-IR J, H, and K_s bands, which have resolutions of 08, 07, and 07, respectively; and the Spitzer Large Area Survey with Hyper-Suprime-Cam (SPLASH; P. L. Capak et al. 2018, in preparation) mid-IR Spitzer/IRAC 3.6 and 4.5 μm, with PSF FWHMs of 166 and 172, respectively.

In this section, we provide a detailed discussion of the individual systems and their subcomponents detected in the HST data (see Figure 1 and Table 3).

AK03: This system has two main UV components separated by ∼1'' (AK03-N and AK03-S), with the F125W image suggesting a bridge connecting the two at an integrated S/N =2.4. The spectroscopic confirmation refers to AK03-N, but AK03-S has a comparable photometric redshift (Smolčić et al. 2015; see Section 5.3). All of these may be considered evidence for a merger. The dust continuum emission is associated with AK03-S and shows two very compact emission peaks (unresolved at the current resolution), whereas AK03-N remains undetected. Therefore, the bulk of the star formation is associated with AK03-S.

AzTEC1: The source shows a compact UV component (AzTEC1-S) and a very faint companion source ∼2'' toward the north, which is detected at 2 < S/N < 3 in all three HST bands (AzTEC1-N). Despite the low S/N of this companion feature being detected in all three bands, the probability of being spurious is ∼10⁻⁵. More importantly, it is detected at S/N > 3 in the HSC r, i, and z bands and Spitzer/IRAC data, confirming that it is a real source. We derived a photometric redshift consistent with lying at the same redshift as AzTEC1-S (Yun et al. 2015) within the uncertainties (see Section 5.3). The rest-frame FIR emission is also compact and centered on AzTEC1-S.

AzTEC5: For this system, three main UV components are detected in all three HST bands (AzTEC5-2, AzTEC5-3, and AzTEC5-4), and a fourth component is detected only in F814W (AzTEC5-1). AzTEC5 is the only source in our sample that lacks spectroscopic confirmation, but photometric redshift estimation indicates a plausible solution for all four components at the same redshift (see Section 5.3). The irregular rest-frame UV morphology of AzTEC5-2 and AzTEC5-4, with emission connecting both in F160W, is suggestive of an ongoing merger. The rest-frame FIR has three emission peaks: two bright peaks associated with AzTEC5-1 and AzTEC5-2 and a fainter peak between them, which is not detected in any HST bands. Besides, the FIR peaks related to AzTEC5-1 and AzTEC5-2 are aligned with the position of the two peaks in the IRAC images, suggesting that the bulk of the stellar mass is associated with these two components, which are probably merging.

AzTEC/C159: As mentioned in Section 2.2, this source was excluded from the HST program and remains undetected in the F814W band image, so we do not have any constraints on its UV morphology. The rest-frame FIR emission is compact and associated with detections in the IRAC bands.

J1000+0234: This system has three main UV components. J1000+0234-N and J1000+0234-S are spectroscopically confirmed at the same redshift (Capak et al. 2008; Schinnerer et al. 2008; A. Karim et al. 2018, in preparation), and J1000+0234-X is a foreground source at z_spec = 1.41 (Capak et al. 2008). An additional companion is detected west of J1000+0234-S in all three HST bands, but the HSC images show diffuse features rather than a concentrated source, consistent with Capak et al. (2008). The north and south components show a connection between them in all three HST bands, suggesting a merger. The rest-frame FIR emission is compact and associated with J1000+0234-N.

Vd-17871: This system has two main UV components ∼15 apart (Vd-17871-N and Vd-17871-S), both with elongated morphologies. Both the north and the south components are spectroscopically confirmed at the same redshift (Smolčić et al. 2015; A. Karim et al. 2018, in preparation). The compact rest-frame FIR emission is associated with the north component.

Having disentangled different stellar components at rest-frame UV wavelengths, we performed photometry in the lower-resolution data sets mentioned in Section 2.6, aiming at fitting the resulting spectral energy distributions (SEDs) to constrain stellar masses for every major stellar component (see Table 3), corresponding to those circled in Figure 1. In the case of Spitzer/IRAC, with a significantly lower resolution, the components appear blended, so it is particularly important to know the number of them to properly deblend the fluxes.

From the g to the K_s bands, the sources are resolved into the stellar components defined from the rest-frame UV HST data, appearing unresolved themselves but separated enough so potential blending is not a concern. To estimate the fluxes in these bands, we carried out aperture photometry. The size of the apertures varied for each component and source, being the same across bands, and correspond to those plotted in Figure 1. We chose the apertures in the K_s band to be as large as possible, enclosing the component we wanted to study without overlapping with a neighboring component aperture. We performed aperture corrections for every band. In order to do so, we traced the growth curve of a PSF in the different bands and applied a correction factor to the fluxes accounting for the missing flux outside the aperture. We performed aperture corrections on each band instead of measuring in PSF-matched data to take advantage of the resolution, important for this kind of multiple-component system, that otherwise would be degraded to the lowest-resolution band. The uncertainties in the magnitudes were derived from empty aperture measurements. To assure a good SED fit, we only use detections above 3σ (upper limits are included in Figure 2).

Zoom In Zoom Out Reset image size

For the blended Spitzer/IRAC 3.6 and 4.5 μm images, we employed the magnitudes from a PSF model using the two-dimensional surface brightness distribution fitting algorithm GALFIT (Peng et al. 2002). We required at least a 5σ detection to perform the fit, which was the case for all sources in the 3.6 and 4.5 μm bands. The number of PSFs was set to the number of stellar components the source has as defined from the HST data, and the PSFs centroids were placed at the positions of K_s-band centroids used as priors, allowing a shift in both the X and Y axes that turns out to be <1 pixel from the initial positions (the IRAC image pixel scale is 06 pixel⁻¹). The uncertainties in the photometry due to the deblending were calculated by performing a number of realizations varying the centroid coordinates randomly within 1 pixel of the best-fit centroid and fixing those coordinates for each realization. Additionally, we checked for detections in the IRAC 5.8 and 8.0 μm bands from the S-COSMOS survey (Sanders et al. 2007), but the sources are not detected at the required 5σ level (upper limits are included in Figure 2).

We fitted the resulting 13-band SEDs (g to 4.5 μm, including the three HST bands) using LePHARE (Arnouts et al. 1999; Ilbert et al. 2006). We adopted Bruzual & Charlot (2003) stellar population synthesis models with emission lines to account for contamination from Hα, which at the redshift probed in this work is redshifted into the IRAC 3.6 μm band. A Chabrier (2003) IMF, exponentially declining star formation histories (SFHs), and a Calzetti et al. (2000) dust law were assumed. We explored a large parameter grid in terms of SFH e-folding times (0.1–30 Gyr), extinction (0 < A_V < 5), stellar age (1 Myr–age of the universe at the source redshift), and metallicity (Z = 0.004, 0.008, and 0.02, i.e., solar). The redshift was fixed to the spectroscopic redshift, if available, or to the photometric redshift if not (see Table 3). In Figure 2, we show the derived SEDs, with the fitted models being in good agreement with the data.

At the redshift of the galaxies, the three HST bands trace the rest-frame UV continuum. This makes it possible to directly determine their spatially resolved UV slopes (β).

In Figure 3, we present β maps constructed by fitting a linear slope to pixels that have S/N > 2 detections in at least two smoothed images (see Section 2.5). The 1σ uncertainty maps (inserts) were constructed by computing β-values in ∼10,000 realizations of the data, varying in each realization the measured pixel flux values within their uncertainties. Note that the pixel size is 006, but the PSF FWHM is 018. Consequently, spatially independent regions are those separated by at least 3 pixels. Since the UV slope maps were obtained using at least two detections in the HST bands, we see more clearly the presence of faint companions toward the north in AK03, AzTEC1, AzTEC5, and J1000+0234, as mentioned in Section 3.

Zoom In Zoom Out Reset image size

In general, the objects present blue UV slopes, but the values are not homogeneous over the extent of the galaxies. The color gradients could be caused by structure in the distribution of dust, stellar age, or metallicity. The relative importance of these cannot be disentangled from the available data, but we expect a patchy dust distribution to be the dominant cause. However, as most of the extent of the rest-frame UV emission is not detected in our ALMA observations, revealing the underlying dust structure in emission would require deeper observations.

The rest-frame FIR dust emission is in all cases associated with the reddest components. These components show evidence of gradients in their UV slopes. AK03-S is redder toward the northeast and bluer at the southwest. J1000+0234-N has an extended redder feature at the northeast. Vd-17871-N is slightly redder toward the southwest and bluer toward the northeast. In AzTEC1-N, the red-to-blue gradient goes along the north–south axis. These color gradients may be due to a star formation gradient with higher dust content toward the redder areas. Another possibility could be close mergers between red and blue galaxies.

In AK03-S, two close FIR peaks are detected. At the current resolution and sensitivity and without dynamical information, we cannot determine whether these are part of a larger dynamical structure like a clumpy disk or remnants of a past interaction/merger. Note that in AzTEC5, we were unable to constrain the resolved UV slope of AzTEC5-1, since it is only detected in F814W, suggesting an extremely high extinction with strong rest-frame FIR emission but also a very blue rest-frame UV component.

The bluer components in all five systems remain undetected in the ALMA continuum. This indicates less dusty star formation.

Spatially integrated values for the UV slopes (see Table 3) show a median and median absolute deviation of β =−0.59 ± 0.57 for the components associated with the ALMA continuum emission (namely, AK03-S, AzTEC1-S, AzTEC5-2, J1000+0234-N, and Vd-17871-N and excluding the AzTEC5-1 and AzTEC/C159 upper limits). The rest of the components are bluer, with β = −1.73 ± 0.54, consistent with estimates of Lyman break galaxies (LBGs) at similar redshift (e.g., Bouwens et al. 2009; Castellano et al. 2012; Finkelstein et al. 2012). By performing the photometry over larger apertures, enclosing all the components per source, we derive β = −0.91 ± 0.85, which is in between the derived values for the red and blue components.

Having identified which UV components are associated with the dust continuum emission, we can relate the SFR in the IR (SFR_IR), tracing the obscured star formation, with that in the UV (SFR_UV), probing the unobscured star formation. The former was obtained from the FIR SEDs presented in Smolčić et al. (2015; Toft et al. 2014 for AzTEC5) covering 100 μm–1.1 mm, updated with new 850 μm fluxes from the S2COSMOS/SCUBA2 survey (J. M. Simpson et al. 2018, in preparation). The procedure is as follows. The FIR SED is modeled using the Draine & Li (2007) dust model (DL07; e.g., Magdis et al. 2012, 2017; Berta et al. 2016), L_IR is calculated by integrating the best fit to the SED in the range 8–1000 μm, and the SFR_IR is obtained using the L_IR-to-SFR_IR conversion from Kennicutt (1998) for a Chabrier IMF. The SFR_UV was calculated employing the Salim et al. (2007) prescription that relates L_UV to SFR_UV for a Chabrier IMF. Note that SFR_UV derived this way corresponds to the observed value, i.e., not corrected from extinction. The total SFR can be accounted for by adding both IR and UV estimates (SFR_IR+UV). Not surprisingly, the star formation is dominated by SFR_IR, with SFR_UV only contributing at the level of 2%–20% to the total SFR (SFR_IR+UV), in agreement with other previous works comparing obscured and unobscured star formation in starburst galaxies (e.g., Puglisi et al. 2017) and galaxies with similar stellar mass (e.g., Whitaker et al. 2017).

In relation to the L_IR and SFR_IR estimates, it is important to consider whether an important fraction of the IR emission could be related to active galactic nucleus (AGN) activity. As reported in Smolčić et al. (2015), none of the sources are detected in the X-ray catalog in the COSMOS field (Chandra COSMOS Legacy Survey; Civano et al. 2016; Marchesi et al. 2016). In terms of the radio emission, Smolčić et al. (2015) studied the IR–radio correlation of the sample, which shows a discrepancy when compared with low-redshift star-forming galaxies due to a mild radio excess. This excess would be in line with studies showing an evolving IR–radio ratio depending on the age of the starburst. In any case, while many SMGs host AGNs, their L_IR is dominated by star formation, with the AGN contribution being <33% (e.g., Pope et al. 2008; Riechers et al. 2014b). This translates into a maximum overestimation in the SFR_IR of 33%, below the SFR_IR+UV sample scatter.

We measured the sizes of the rest-frame FIR dust continuum emission by modeling the ALMA continuum images using GALFIT. Sérsic and PSF profiles were fitted to compare both resolved and unresolved modeling of the objects. The only object that was better fitted by a point source than a Sérsic model (and thus unresolved) is AK03. For this galaxy, we derived an upper limit on the size from the PSF.

For the rest of the galaxies, we fitted models with the Sérsic index fixed to n = 0.5, 1, and 4, corresponding to a Gaussian, exponential disk, and de Vaucouleurs profiles, respectively, and also leaving the index free. The size of the emitting regions was obtained through the effective radius of the models (r_e). We cannot constrain which Sérsic index better explains the data at the current resolution and S/N. From higher-resolution observations, Hodge et al. (2016) found a median Sérsic index of n = 0.9 ± 0.2 for a sample of 15 SMGs and concluded that the dust emission follows an exponential disk profile. Motivated by this, we fixed n = 1 to report the rest-frame FIR sizes for our sample in Table 4. We also performed fits varying the axis ratio (b/a) and found that no particular value with b/a ≥ 0.3 fitted the data better than others, so we fixed it to the circular value b/a = 1. We take into account the possible systematic errors associated with the assumed Sérsic index and axis ratio in the listed effective radii errors. These were computed by adding in quadrature the statistical uncertainty from GALFIT for the circular disk model and the difference between this model and the full range of models with varying n and b/a. Therefore, the uncertainties conservatively account for the inability of the data to robustly constrain the detailed shape of the surface brightness profiles. We note that the ALMA continuum fluxes are consistent with the 850 μm fluxes from the S2COSMOS/SCUBA2 survey (J. M. Simpson et al. 2018, in preparation); thus, there is no evidence for resolved-out or missing flux that could affect the size estimates.

Finally, we cross-checked the results, analyzing the data directly in the (u, v) plane employing UVMULTIFIT (Martí-Vidal et al. 2014) following the procedure described in Fujimoto et al. (2017). In this case, for a direct comparison with the GALFIT image-plane fits, we also assumed a circular disk model (to obtain secure results, we omit AK03 and AzTEC5 for this comparison, as they show two and three components, respectively, in our ALMA continuum images). We find that these estimates are in agreement with the results derived in the reconstructed images using GALFIT (see Table 4). In the following, we use the estimates derived from GALFIT for further calculations.

The median and median absolute deviation of the size estimate for our sample are then r_e = 0.70 ± 0.29 kpc at ∼870 μm, which corresponds to ∼160 μm rest frame at z = 4.5 (excluding the AK03 upper limits and only considering the brightest peak in AzTEC5, associated with AzTEC5-2). This result is in good agreement with that of Ikarashi et al. (2015), who found similar compact sizes of r_e = 0.67_−0.14^+0.13 kpc for a sample of 13 1.1 mm selected SMGs at a comparable redshift 3 < z < 6. Oteo et al. (2017) presented an average value of r_e = 0.91 ± 0.26 kpc (converting the reported FWHM into a circularized effective radius) in a sample of 44 DSFGs at z ∼ 4–6 observed at ∼870 μm and selected as Herschel 500 μm risers (SED rise from 250 to 500 μm). On the other hand, the typical sizes derived for SMGs at a median redshift of z ∼ 2.5 were reported to be r_e = 1.8 ± 0.2 kpc from Hodge et al. (2016) and r_e = 1.2 ±0.1 kpc from Simpson et al. (2015), both targeting ∼870 μm. This suggests that SMGs may be more compact at z > 3 than at z < 3 (e.g., Fujimoto et al. 2017; Oteo et al. 2017). Other individual sources at z > 4 also point toward very compact dust continuum emission (e.g., Riechers et al. 2013, 2014a; Díaz-Santos et al. 2016) and pairs of compact interacting starburst galaxies detected in gas and dust continuum, which suggests a gas-rich major merger (e.g., Oteo et al. 2016; Riechers et al. 2017). Spilker et al. (2016) found no evidence for a difference in the size distribution of lensed DSFGs compared to unlensed samples from a sample of 47 DSFGs at z = 1.9–5.7. Our results are also similar to the compact morphologies of local ultra-luminous infrared galaxies (ULIRGs; r_e = 0.5 kpc; Lutz et al. 2016) at 70 μm rest frame. We note that caution should be exercised when comparing samples tracing different rest-frame FIR wavelengths and based on different selection methods. Another caveat for a fair comparison is the stellar mass, since more massive galaxies are typically larger (e.g., van der Wel et al. 2014).

From the SFR_IR obtained for these sources (see Table 3) and their rest-frame FIR sizes, we calculated the SFR surface density (Σ_SFR = 0.5SFR/πr_e,circ²; see Table 4). Ranging from Σ_SFR = 150 to 1300 M_⊙ yr⁻¹ kpc⁻² (excluding AK03 lower limits), the most extreme cases are AzTEC1, AzTEC5-1, and Vd-17871-N, but the last two are poorly constrained due to the large uncertainty on their sizes. At such extreme values, they are candidates for Eddington-limited starbursts (Σ_SFR ∼ 1000 M_⊙ yr⁻¹ kpc⁻²; Andrews & Thompson 2011; Simpson et al. 2015).

The dust masses derived for this sample are very high at ∼10⁹ M_⊙ (see Table 3). Dust masses are a free parameter in the DL07 model employed, controlling the normalization of the SED. In terms of the dust opacity, DL07 assumes optically thin dust (τ < <1) at all wavelengths (e.g., Magdis et al. 2012, 2017; Berta et al. 2016). Very high dust masses combined with the small sizes derived for the dust-emitting regions implies very high dust mass surface densities (Σ_dust = 0.5M_dust/πr_e,circ²) with values in the range Σ_dust = 0.33–5.8 × 10⁹ M_⊙ kpc⁻² (see Table 4) and, consequently, very high extinction.

We calculated the expected extinction assuming that the dust is distributed in a sheet with uniform density. We inferred the mean extinction from the dust mass surface density–to–extinction ratio (Σ_dust/A_V). To calculate Σ_dust/A_V, we assumed a gas-to-dust mass ratio (GDR) appropriate for SMGs of GDR = 90 (Swinbank et al. 2014) and the gas surface number density–to–extinction ratio N_H/A_V = 2.2 × 10²¹ cm⁻² mag⁻¹ (Watson 2011). Therefore, Σ_dust/A_V = (N_H/A_V) · m_H/GDR = 2.44 M_⊙ pc⁻² mag⁻¹. With this number, the mean extinction is $\langle {A}_{V}\rangle ={{\rm{\Sigma }}}_{\mathrm{dust}}/2.44$. The values for our sample are extreme $\langle {A}_{V}\rangle =130\mbox{--}2400\,\mathrm{mag}$, even when the numbers are halved to account for the dust behind the sources (see also Simpson et al. 2017).

Comparing Figures 1 and 3, we see that while the dust emission is always associated with the reddest (likely most dust-extincted) component, in most cases it is not centered on the reddest part of that component (with the possible exception of Vd-17871). This suggests that the extinction seen in emission and absorption is disconnected, consistent with the expected extreme A_V, which implies that no emission can escape at any wavelength.

The fact that we do see blue UV emission at the peak of the dust emission suggests that a fraction of the light is able to escape due to a clumpy dust distribution and/or that the dust and stars are seen in different projections; e.g., the stars responsible for the UV emission could be in front of the dusty starbursts.

In any case, it is clear that the rest-frame UV and FIR emissions are spatially disconnected and originate from a different physical region. This implies that the dust as seen in absorption from the UV slope inhomogeneities in Figure 3 is not tracing the dust seen in emission from the ALMA continuum.

The IR-to-UV luminosity ratio, commonly referred to as IR excess (IRX = L_IR/L_UV), is known to correlate with the UV continuum slope (β). This so-called Meurer relation (Meurer et al. 1999, hereafter M99) is well established for normal star-forming galaxies (e.g., Overzier et al. 2011; Takeuchi et al. 2012; Casey et al. 2014b). Its origin is thought to be that galaxies get redder as the dust absorbs the rest-frame UV emission and reradiates it at IR wavelengths. For galaxies on the relation, the amount of dust absorption can thus be directly inferred from the UV slope. Therefore, in the absence of FIR data, the relation can be used to obtain a total extinction-corrected SFR from UV data (e.g., Bouwens et al. 2009). Furthermore, this relation physically motivates energy balance codes that require that dust extinction inferred from rest-frame UV–optical SED fits must match the observed emission measured at IR wavelengths (e.g., da Cunha et al. 2008).

Spatially unresolved observations have shown that DSFGs do not follow the M99 relation (e.g., Buat et al. 2005; Howell et al. 2010; Casey et al. 2014b). An excess of dust and UV/FIR decoupling have both been suggested as a possible origin of the offsets by Howell et al. (2010), who showed that the deviation from the nominal M99 relation (ΔIRX) increases with L_IR but does not correlate with L_UV. Following this argument, the authors postulated that a concentration parameter might correlate with ΔIRX as an indicator of the decoupled UV/FIR. Casey et al. (2014b) reinforced these results, showing that the deviation from the M99 relation increases with L_IR above a threshold of $\mathrm{log}({L}_{\mathrm{IR}}/{L}_{\odot })\gt 11.0$. Faisst et al. (2017a) proposed that the blue colors of sources with high IRX values could be due to holes in the dust cover, tidally stripped young stars, or faint blue satellite galaxies. In addition, simulations propose recent star formation in the outskirts and low optical depths in UV-bright regions as plausible explanations of the offset (Safarzadeh et al. 2017; Narayanan et al. 2018). Simple models placing a dust screen in front of a starburst have been studied to provide a detailed explanation of all the possible effects that might lead to a deviation in the IRX–β plane (Popping et al. 2017b).

The sample studied here has IR luminosities in the range $\mathrm{log}({L}_{\mathrm{IR}}/{L}_{\odot })=12.1$–13.4, above the mentioned threshold $\mathrm{log}({L}_{\mathrm{IR}}/{L}_{\odot })\gt 11.0$, and the spatially resolved rest-frame UV/FIR data make it possible to study the origin of the DSFG offsets in the IRX–β plane (see Figure 4).

Zoom In Zoom Out Reset image size

To confirm that the galaxies in this sample are representative of previous DSFG studies in spatially unresolved data, we first derived UV and IR luminosities in large apertures enclosing all the components of each source. In Figure 4, these measurements are plotted as large open symbols, confirming that the sample does not follow the M99 relation and is located in the same region as previous spatially unresolved measurements for DSFGs (e.g., Casey et al. 2014a at z > 2). Second, we take advantage of the spatial resolution to pinpoint the origin of the FIR emission and recalculate the UV luminosity in smaller apertures defined by the 3σ contour in the ALMA images (ALMA apertures). In this case, both the HST and ALMA images were PSF-matched as described in Section 2.4.

In Figure 4, we plot the sample of DSFGs at z > 2 from Casey et al. (2014b) for comparison. Note that this study employed similar methods to obtain L_IR, L_UV, and the UV slopes as we did: L_IR by integrating over the wavelength range 8–1000 μm and using a single temperature modified graybody plus mid-IR power law, which properly accounts for the warm dust contribution as the DL07 dust model; L_UV by interpolating the observed photometry to rest-frame 1600 Å; and the UV slopes by fitting a power law to the photometry, which is equivalent to our linear fit in magnitude space. Additionally, we include other IRX–β relations from the literature: the original M99 and follow-up corrections (e.g., Overzier et al. 2011; Takeuchi et al. 2012), although the methodology they followed to obtain the quantities shown in the IRX–β diagram slightly differs from that of Casey et al. (2014b) and us.

All of the galaxies have higher IRX in the ALMA aperture than in the large aperture. This is expected from their smaller extent in the rest-frame FIR compared to that in the rest-frame UV, which effectively lowers the L_UV contribution to the L_IR/L_UV ratio. Furthermore, three of the galaxies have redder UV slopes and two have similar UV slopes in the ALMA apertures compared with the large apertures. Again, this can be understood as a result of removing the contribution from the extended irregular UV features and companion satellite galaxies that appear bluer than the dust-emitting region detected in ALMA. These results agree with the model proposed in Faisst et al. (2017a) to explain the blue colors of DSFGs with high IRX.

On the other hand, even after accounting for the correction that implies going from the large to the ALMA aperture, our sample does not follow the M99 relation and lies 1.75 dex (median) above it. However, while the rest-frame FIR dust continuum emission is associated with the reddest component in the mergers, it is generally not centered on the reddest part of the component, and the component is too blue to be consistent with a physical connection between the dust seen in emission and absorption, suggesting that the UV and FIR emissions of DSFGs are spatially disconnected.

This provides morphological and geometrical evidence for the origin of the DSFG offsets from the M99 relation (see also Chen et al. 2017) being consistent with the extreme extinction expected from the compact and intense dust emission for this sample (see Section 4.3), implying that UV emission should be expected not to escape the starbursts.

A possible scenario for the origin of the UV and FIR emissions could be a patchy dust distribution causing some of the UV to be completely extincted and some to leak relatively unextincted, in a similar way as proposed by the holes in the dust cover by Faisst et al. (2017a).

The UV/FIR lack of spatial coincidence has important implications for energy balance codes, as noted by Hodge et al. (2016), where the detected stellar light will have no information about the obscured starburst (Simpson et al. 2015).

Therefore, the results here support that IRX and β are unrelated for such FIR-bright sources and that extinction-correction prescriptions based on the nominal IRX–β relation are inappropiate for DSFGs.

In Section 4.1, we interpreted the UV slope differences over the source extent as variations in the dust content not detected in emission in the ALMA observations. It is possible that this regime of star formation is compatible with the M99 relation. In order to check this, we calculated the expected L_IR below the 3σ dust continuum detection limit over the components detected in the rest-frame UV for each source by rescaling their FIR SEDs (Toft et al. 2014; Smolčić et al. 2015). The resulting upper limits lie above the M99 relation for all cases, not being useful for putting constraints about whether these galaxies follow M99 or lie above or below it, a subject of main focus in current studies (e.g., Capak et al. 2015; Barisic et al. 2017; Faisst et al. 2017a; Fudamoto et al. 2017).

Major mergers between gas-rich galaxies are often assumed to be the triggering mechanism for starburst galaxies, as local universe IR-luminous galaxies are exclusively associated with major mergers with L_IR > 10^11.5 L_⊙ (e.g., Sanders & Mirabel 1996). The multiplicity of close, approximately equally bright galaxies in the HST images studied here would naively support a similar triggering mechanism at z > 4. However, as the images trace the rest-frame UV, a stellar mass analysis of the individual merging components is needed to test this picture.

In Table 3, we list the stellar masses of the stellar components of each system derived from the SED fits described in Section 3.2. Also listed is the stellar mass ratio relative to the most massive component in the system (M_*,prim).

The median stellar mass of the most massive component is $\mathrm{log}({M}_{* }/{M}_{\odot })=10.49\pm 0.32$ (where the uncertainty is the median absolute deviation). For the remaining less-massive components, the median is $\mathrm{log}({M}_{* }/{M}_{\odot })=9.56\pm 0.10$. A stellar mass ratio of 1:3–4 is often adopted to distinguish between major and minor mergers (e.g., Conselice et al. 2003; Tacconi et al. 2008; Kaviraj et al. 2009; Lotz et al. 2011; Man et al. 2016). Adopting this definition, AzTEC5 is formally classified as a major merger, with a stellar mass ratio for the two most massive components (AzTEC5-1 and AzTEC5-2) of M_*/M_*,prim = 3.0. Vd-17871 could be classified as a major or minor merger, depending on the exact distinction ratio (M_*/M_*,prim = 3.5). The rest of the systems are consistent with undergoing at least one minor merger (including AzTEC5, which might undergo minor merging with AzTEC5-3 and AzTEC5-4). Furthermore, it is important to note that regardless of the precise distinction between major and minor mergers, the components detected in dust continuum with ALMA are undergoing starbursts with SFRs that overwhelm those of the companions; therefore, the stellar mass ratios are expected to decrease. Taking this into account, all systems could be classified as minor mergers (except AzTEC5-1 and AzTEC5-2, both starbursting systems).

In addition to the components that were bright enough to estimate stellar masses, AK03, AzTEC5, and J1000+0234 present additional low-S/N companions detected in one or more of the HST images (marked with arrows in Figure 1), which may be additional minor merger components if they are at the same redshift. The residuals in the modeling of the Spitzer/IRAC images do not show significant detections at their positions; thus, they must be less massive than the detected companions. In fact, the HST images display 2 < S/N < 3 potential additional low-mass components in the case of F814W particularly, as expected if they are small, blue star-forming galaxies. If their redshifts are confirmed, it will be further evidence for the starbursts in z ∼ 4.5 SMGs being triggered by multiple minor mergers, a picture consistent with living in overdense environments (e.g., Blain et al. 2004; Smolčić et al. 2017). Indeed, Smolčić et al. (2017) showed evidence that AzTEC1, AzTEC5, J1000+0234, and Vd-17871 have statistically significant small-scale overdensities.

Note, however, that these results do not rule out that major mergers played a role in triggering these starbursts if they have already coalesced or are so close that they are not resolved in the HST and ALMA data. Indeed, the multiple FIR peaks in AK03 and AzTEC5 and the color gradients observed in the most massive components of the systems (most prominently in AzTEC1, J1000+0234, and Vd-17871) are consistent with such a picture.

Previous estimates of the stellar mass of the galaxies in this sample, derived using MAGPHYS (da Cunha et al. 2008), led to a median value of $\mathrm{log}({M}_{* }/{M}_{\odot })=10.92\pm 0.13$ (Toft et al. 2014; Smolčić et al. 2015). This is ∼0.4 dex higher than our derived median value for the most massive component. Adding up all the components per source, the median total stellar mass would be slightly higher, $\mathrm{log}({M}_{* }/{M}_{\odot })\,=10.63\pm 0.11$, but still ∼0.3 dex lower than the previous estimates for this sample. Recent results from Miettinen et al. (2017), also employing MAGPHYS, are also systematically higher by at least 0.3 dex for the sources in common with our sample (AzTEC1, AzTEC5, and J1000+0234). Such systematic discrepancies are consistent with the expected overestimation of MAGPHYS-derived stellar masses and slight underestimation of exponentially declining models employed here, according to Michałowski et al. (2014) SMG stellar mass studies from simulated data sets.

We also compared our stellar mass estimates with those listed in the 3D-HST survey catalog (Brammer et al. 2012; Skelton et al. 2014; Momcheva et al. 2016) for the sources covered in the CANDELS fields (e.g., AK03 and AzTEC5) and the COSMOS2015 catalog (Laigle et al. 2016). In general, the catalogs successfully extract the majority of the components for these complex objects and list photometric redshifts consistent with the available spectroscopic redshifts (see Table 3). However, for a subset, we found significant discrepancies in the derived stellar masses. The discrepancy might be due to the different approach in the photometry measurements. While we measured fluxes in apertures carefully chosen to minimize the effect of blending and applied aperture corrections, COSMOS2015 employs automated PSF-matched photometry, which can be more contaminated by blending of close objects.

Furthermore, J1000+0234-N is not in the COSMOS2015 catalog, and the bulk of its stellar mass is associated with J1000+0234-S (likely due to a mismatch between the Spitzer/IRAC and optical/near-IR data). AzTEC/C159 is also missing from the catalog, due to its extreme faintness in the optical/near-IR. Similarly, there is no entry corresponding to the location of AzTEC5-1 in either 3D-HST or COSMOS2015. The absence and misidentifications of massive and optically faint sources could affect the photometry and, thus, the stellar mass estimates. It could also affect the stellar mass functions at high redshifts (e.g., Davidzon et al. 2017).

J1000+0234 is also present in the recent work by Brisbin et al. (2017), and the assigned shorter-wavelength counterpart to the ALMA detection is also J1000+0234-S, since J1000+0234-N remains undetected. This indicates that significant offsets between sub-mm/radio sources and UV/optical/near-IR counterparts could indeed be due to the presence of multiple blended, and perhaps merging, components if the depth and resolution of the data are not enough to detect all those components (provided a good relative astrometry between the different instruments).

Compared with previous estimates of the average stellar masses of SMGs, our results are in line with studies indicating that most SMGs have M_* < 10¹¹ M_⊙ (e.g., Hainline et al. 2011; Wardlow et al. 2011; Casey et al. 2013; Simpson et al. 2014). Other studies also report higher values of M_* > 10¹¹ M_⊙ for z ∼ 4.5 sources (Michałowski et al. 2010, 2012, 2017). The median stellar mass of the satellite galaxies is consistent with estimates for faint LBGs at similar redshifts (Magdis et al. 2010).

Stellar masses of highly obscured starburst galaxies are notoriously difficult to estimate. In this work, we took advantage of high-resolution HST imaging to identify the positions of multiple stellar components in the systems, which in turn was used to deblend the rest-frame optical Spitzer/IRAC fluxes that are tracing the stellar mass available for these high-redshift systems. However, our stellar mass estimates are potentially subject to a number of additional systematic uncertainties.

One caveat is that some of the components lack spectroscopic confirmation. That is the case for AK03-S, AzTEC1-N, and all components of AzTEC5. When possible, we assumed that these components were at the same redshift as their spectroscopically confirmed companions. For AK03-S, Smolčić et al. (2015) found z_phot = 4.40 ± 0.10 or 4.65 ± 0.10, depending on the template used. Therefore, the two components are likely at the same redshift. AzTEC1-N is a very faint component with S/N < 3 in all HST bands, but it is detected above this threshold in the HSC r, i, and z bands and the IRAC bands, where the residuals from AzTEC1-S fitting showed that there is indeed a secondary component toward the north. We derived a photometric redshift consistent with being at the same redshift as AzTEC1-S within the uncertainties. Its probability distribution peaks at 3.77_−0.22^+0.32 (where the uncertainties are the 1σ percentiles of the maximum-likelihood distribution), being not null in the redshift range 3.2 < z < 4.7. In the case of AzTEC5, none of the components have spectroscopic redshifts, but the 3D-HST survey catalog (Brammer et al. 2012; Skelton et al. 2014; Momcheva et al. 2016) lists z_phot = 3.63_−0.15^+0.14 for AzTEC5-2, z_phot = 4.02_−0.08^+0.08 for AzTEC5-3, and z_phot =3.66_−0.43^+0.40 for AzTEC5-4. Therefore, it seems plausible that all components in AzTEC5 lie at the same redshift within the uncertainties.

Another caveat in the stellar mass estimates comes from the assumptions made in the SED fits. Michałowski et al. (2014) studied the importance of the assumed SFHs (see also Hainline et al. 2011) over several SED fitting codes, concluding that the exponentially declining SFHs used here are able to recover the stellar masses of their simulated SMGs with slight underestimation and significant scatter. Regardless of the model employed, the derived photometry and the color of the sources already indicate that there is a component more massive than the other. The most massive components have higher IRAC fluxes, and they are also redder than their fainter IRAC companions.

Given the extreme dust mass surface densities derived for this sample (see Table 4), if the stars formed in situ in the starburst that created the dust, it is possible that some stellar mass is so obscured that it is not detectable even by IRAC and, thus, not accounted for in the SED fitting. Higher spatial resolution rest-frame FIR continuum observations would be needed to disentangle the underlying structure of the dust-emitting region and measure its degree of homogeneity or clumpiness. This could reveal how much of the stellar light is completely obscured beneath the dust and the implied systematic error in the derived stellar masses. To estimate how big this effect could be, using the empirical dust–to–stellar mass ratio (DTS) for local ULIRGs in Calura et al. (2017) $\mathrm{logDTS}=-2.83$, the median stellar mass of this sample would increase to $\mathrm{log}({M}_{* }/{M}_{\odot })\sim 11.6$. However, assuming the ratio from simulations in Popping et al. (2017c), $\mathrm{logDTS}\sim -1.8$, the effect would not be that significant, increasing to $\mathrm{log}({M}_{* }/{M}_{\odot })\sim 10.9$.

Over the last decade, several studies have uncovered a tight correlation between the SFR and the stellar mass of star-forming galaxies, the so-called main sequence (MS) of star formation (e.g., Daddi et al. 2007; Elbaz et al. 2007; Noeske et al. 2007). Strong outliers to the MS are present at all redshifts, and this is often used as a formal definition of starburst galaxies. These systems exhibit elevated specific star formation rates (sSFRs) compared with typical MS galaxies. For the components with ALMA detection, from the total SFR_IR+UV and stellar masses, we obtain sSFR = 2.5–100 Gyr⁻¹. Considering the MS as defined in Schreiber et al. (2015), the distance to the MS is in the range sSFR/sSFR_MS = 0.5–22, calculated at the redshift of each source. Consequently, all of the sources studied here would formally fall into the starburst regime, with AK03 on the MS but also consistent with the starburst region, given its large SFR uncertainty (see Figure 5). If an important fraction of the stellar mass is undetectable hidden beneath the dust, the objects will move toward smaller distances to the MS, as represented by the bottom arrows in Figure 5.

Zoom In Zoom Out Reset image size