© ESO, 2014

1. Introduction

Submillimeter (submm) and millimeter observations are efficient in detecting and studying dusty, star-forming galaxies, due to the effect of negative K-correction, which results in nearly constant observed brightness for galaxies with same infrared (IR) luminosity over a broad range of redshifts. However, most current deep submm surveys are limited to the brightest sources and submm-selected galaxies (SMGs; Blain et al. 2002), due to the limited spatial resolution and sensitivity of submm observations. SMGs are massive, highly dust obscured galaxies with extreme star formation rates (SFRs) of order 10³ M_⊙ yr^-1 (e.g., review by Blain et al. 2002; Casey et al. 2014), and are generally thought to represent the progenitors of local massive elliptical galaxies. While spectroscopic studies of SMGs originally gave a median redshift of z ~ 2.5 (Chapman et al. 2005), recent deep submm/mm continuum and radio observations show evidence for a significant population of higher redshift massive starbursts (SBs; e.g., Dannerbauer et al. 2004; Smolčić et al. 2012; Swinbank et al. 2014; Dowell et al. 2014), extending the redshift peak beyond z = 3. A substantial number of z> 4 SMGs have been identified to date (e.g., Dannerbauer et al. 2008; Daddi et al. 2009a,b; Capak et al. 2008, 2011; Schinnerer et al. 2008; Coppin et al. 2009; Knudsen et al. 2010; Riechers et al. 2010a, 2013; Smolčić et al. 2011; Walter et al. 2012; Combes et al. 2012; Vieira et al. 2013). The surface density of these galaxies is found to be significantly higher than that expected from theoretical models (e.g., Baugh et al. 2005; Hayward et al. 2013), suggesting that current models of galaxy formation underpredict the number of high-redshift starbursts.

Observations of the molecular gas in high-redshift galaxies reveal that while SMGs are highly gas-rich systems (Tacconi et al. 2008), the gas fractions of these systems are comparable to those of typical massive galaxies at similar epochs (~40–60%; Daddi et al. 2008, 2010), implying that SMGs have higher star formation efficiencies (SFEs; Daddi et al. 2010; Genzel et al. 2010). However, these results are complicated by the large uncertainties associated with the CO-to-H₂ conversion factor α_CO¹, which likely changes between normal disk galaxies and starbursts (see review by Bolatto et al. 2013; Carilli & Walter 2013). In the literature, an “ULIRG-like” value of α_CO = 0.8 (Downes & Solomon 1998; but see Papadopoulos et al. 2012 for a higher value of α_CO) is widely adopted for SMGs due to the lack of direct measurements, while a value of α_CO ~ 4 is favored for Milky Way and normal galaxies. This carries a significant uncertainty since high redshift SMGs may be dramatically different from local ULIRGs, given the more extended gas distribution and different physical conditions revealed in some SMGs (Riechers et al. 2011; Ivison et al. 2011; Carilli & Walter 2013; Scoville et al. 2014). Therefore, it is of significant importance to obtain a direct calibration of α_CO, since well-determined molecular gas masses are critical to study the variations in physical properties across the galaxy populations at high redshift.

Because of the extreme high specific star formation rates (sSFRs), some of the most luminous SMGs are placed as outliers above the main sequence (MS) of star formation, which is a tight correlation observed between the stellar mass and the SFR over a broad range of redshifts (e.g., Noeske et al. 2007; Elbaz et al. 2007; Daddi et al. 2007b; Rodighiero et al. 2010, and references therein). While galaxies on the MS are thought to form stars gradually with a long duty cycle and represent the bulk of the galaxy population, starbursts exhibit very intense and rapid star formation activity, likely driven by mergers (e.g., Daddi et al. 2007a,b; Tacconi et al. 2008, 2010; Elbaz et al. 2011; Rodighiero et al. 2011). Recent studies on the molecular gas of z> 3 Lyman break galaxies (LBGs) found these galaxies to be rather deficient in CO emission for their star formation activity (Magdis et al. 2012b; Tan et al. 2013). Similar results have also been reported for a luminous LBG at z = 6.595 called “Himiko”, for which the 1.2 mm dust continuum and [CII] 158 μm emission are much lower than predicted by local correlations and measured SFRs (Ouchi et al. 2013). It has been found that normal galaxies at z> 3 are increasingly metal poor, with metallicities dropping by about 0.6 dex as compared to local galaxies of similar stellar mass (Mannucci et al. 2010; Sommariva et al. 2012; Troncoso et al. 2014). This may suggest that metallicity effects could be a probable explanation for the deficit of CO emission, since the photodissociation of CO by far-UV radiation is enhanced at low metallicity (Leroy et al. 2011; Genzel et al. 2012; Narayanan et al. 2012; Bolatto et al. 2013). The decrease of CO emission in z> 3 normal galaxies for their IR luminosity is also predicted by simulations with a galaxy-formation model (Lagos et al. 2012), a result driven by the low metallicities in such objects. Similar detailed study of a local metal poor star-forming galaxy, I Zw 18, concluded that it would be much harder than hitherto anticipated to detect gas and dust in high-redshift galaxies like Himiko (several tens of days of integration with the complete ALMA; see Fisher et al. 2014), if assuming I Zw 18 is an analog of primitive galaxy population in the early Universe. Dust is expected and observed to be well mixed with gas in the interstellar medium (ISM), because it is composed of metals and regulates the gas phase abundances of the elements through accretion and destruction processes (Draine & Li 2007). The comparison of dust and gas properties of galaxies at different redshifts is thus crucial to explore the interplay between dust, gas, and metals in the ISM, and allow us to achieve a better understanding of galaxy evolution throughout cosmic time. For galaxies at z> 3, to date only very few luminous SMGs have been detected in both dust continuum and gas emission (e.g., Dannerbauer et al. 2008; Daddi et al. 2009a,b; Coppin et al. 2009, 2010; Walter et al. 2012; Riechers et al. 2013), while no dust continuum detection is available for indisputably normal galaxies at z> 3.5.

GN20 is one of the brightest SMGs in the GOODS-N field (Pope et al. 2006), of which the redshift (z = 4.055) was established by a serendipitous detection of its CO(4–3) emission (Daddi et al. 2009b). Two additional SMGs, GN20.2a and GN20.2b, were found to lie within ~25′′ of GN20 (projected physical separation ~170 kpc) and have redshifts of z ~ 4.055 ± 0.005. These two galaxies are separated by only a few arcseconds and hence are not spatially separated from each other in existing submm images (e.g., SCUBA 850 μm). Daddi et al. (2009b) also found 14 B-band dropouts (roughly z ~ 4) lying within 25′′ from GN20, which corresponds to an overdensity of 5.8σ in the GOODS-N field, suggesting a massive proto-cluster environment at z ~ 4.05, just 1.6 Gyr after the Big Bang. All three massive SMGs in the GN20 proto-cluster have been detected in CO emission, indicative of large amounts of molecular gas (Daddi et al. 2009b; Carilli et al. 2010, 2011; Hodge et al. 2012, 2013), feeding vigorous ongoing star formation (SFR ~ a few to ten times 100 M_⊙ yr^-1). The deep, high-resolution CO(2–1) observations reveal a clumpy, extended gas disk (14 ± 4 kpc in diameter) for GN20 (Hodge et al. 2012), and extended gas reservoirs (~5–8 kpc) for GN20.2a and GN20.2b (Hodge et al. 2013). Hodge et al. (2012, 2013) have attempted to constrain the estimate of CO-to-H₂ conversion factor by dynamical analysis, deriving α_CO of ~1–2 for these three galaxies. For the dust continuum emission, however, only GN20 has been reported in Herschel and (sub)mm detections (Daddi et al. 2009b; Magdis et al. 2011). Here we use PACS and PdBI mm data for GN20.2a and GN20.2b to study the dust properties of these two galaxies. By combining with the molecular gas properties and dynamical analysis, we aim to achieve a more comprehensive understanding of the nature of the massive SMGs in the GN20 proto-cluster environment. We further investigate the metallicity effects on molecular gas and dust emission by comparison of CO luminosity-to-dust mass ratio and dust-to-stellar mass ratio between normal galaxies and starbursts.

This paper is organized as follows. In Sect. 2 we present the new PdBI CO(4–3) observations and the reduction of data for the SMGs GN20, GN20.2a, and GN20.2b. Section 3 presents the results of CO and millimeter continuum observations, the methods used to compute SFR, dust mass, stellar mass, metallicity, dynamical mass, and CO-to-H₂ conversion factor. This section also describes the derived physical properties including sSFR, SFE, radio-IR correlation constraints, gas fraction, and gas depletion timescales. In Sect. 4 we discuss the nature of GN20, GN20.2a, and GN20.2b based on the physical properties of optical morphology, molecular gas, dust, and dynamical constraints. We further discuss the implications for the cosmic evolution of dust content in galaxies in Sect. 5. Finally, we summarize our results in Sect. 6. We adopt a cosmology with H₀ = 71 km s^-1 Mpc^-1, Ω_M = 0.3, Ω_Λ = 0.7, and a Chabrier (2003) initial mass function (IMF) throughout the paper.

Fig. 1 CO(4–3) spectra binned in steps of 25 km s-1 for GN20 (left), 50 km s-1 for GN20.2a (middle) and GN20.2b (right). The yellow color indicates the velocity ranges where line emission is detected. These regions have been used to derive the integrated fluxes. The red lines show best-fitting double Gaussian profiles to the spectra, and the blue dashed lines show the fits with a single Gaussian. Zero velocities correspond to the redshifts listed in Table 1.

2. Observations and data reductions

We used the IRAM PdBI to observe the CO(4–3) emission in the GN20 field in the AB, C and D configurations. The AB configuration observations were pointed close to a nearby z = 1.5 galaxy and have been reported in Daddi et al. (2009b), while the new C and D configuration observations, carried out in June 2009 and January–April 2013, were centered near a z = 4.058 LBG, which is located ~16′′ south of GN20 (Daddi et al. 2009b). Both GN20 and its companions, GN20.2a and GN20.2b, are within the primary beam of the PdBI observations. The primary beam attenuation (PBA) for GN20, GN20.2a, and GN20.2b are 1.29, 1.44, and 1.20 in C configuration observations, and 1.27, 1.18, and 1.06 in D configuration observations, respectively. All observations were tuned at a central frequency of 91.375 GHz and performed in dual polarization mode with six antennas available. For the observations performed in 2013, a 3.6 GHz WIDEX correlator was used. Further details on the observations can be found in Tan et al. (2013).

We reduced the data with the GILDAS software packages CLIC and MAPPING. After flagging bad and high phase noise data, and correcting for the PBA, the total effective on-source integration time for GN20, GN20.2a, and GN20.2b are 14.6 h, 11.9 h, and 15.1 h, respectively. We combined the C and D configuration data with published AB configuration observations to maximize the sensitivity. The spectra shown in Fig. 1 were extracted from the combined UV data by using a circular Gaussian model with a fixed full width at half maximum (FWHM) of 0.72′′, 0.53′′, 0.88′′ for GN20, GN20.2a, and GN20.2b (Carilli et al. 2010; Hodge et al. 2013), respectively. The CO(4–3) images integrating over the velocity range where emission was detected and the 3.3 mm maps averaging over the line-free channels for the three galaxies are shown in Fig. 2. In order to avoid the contamination from the side lobes of the bright GN20 galaxy, the emission at the position of GN20 was fitted in the UV data with a circular Gaussian source model and subtracted before creating the images of GN20.2a and GN20.2b shown in Fig. 2. The sensitivities of combined data set at the position of GN20, GN20.2a, and GN20.2b are 0.41, 0.63, and 0.52 mJy beam^-1 per 25 km s^-1 channel.

We make use of the 1.2 mm and 2.2 mm continuum observations of the GN20 field with the PdBI. The 2.2 mm continuum data of our three galaxies haven been reported in Carilli et al. (2010). More details on the 1.2 mm observations and reduction will be given in Riechers et al. (in prep.).

3. Results and analysis

3.1. CO(4–3) emission properties

Fig. 2 PdBI maps (from left to right) of CO(4–3) line, 3.3 mm, and 2.2 mm continuum emission for GN20 (top), GN20.2a (middle), and GN20.2b (bottom). The CO(4–3) images are averaged over the observed velocity range of CO emission (see Fig. 1) and were cleaned. The 2.2 mm and 3.3 mm continuum maps were created by averaging over the line-free emission channels and were not cleaned given the relatively low S/N. The crosses show the positions where we extract CO spectra. For GN20, the contours level of CO(4–3) are shown increase by a factor of 1.5 starting at ±2σ, with positive (negative) contours shown as solid (dashed) lines. For all the continuum maps and CO(4–3) maps of GN20.2a and GN20.2b, contours start at ±2σ in steps of 1σ. Beam sizes are displayed in the lower right corner of each integrated map.

The CO detections of GN20, GN20.2a, and GN20.2b at various transitions have been presented in some recent studies (Daddi et al. 2009b; Carilli et al. 2010, 2011; Hodge et al. 2012, 2013). Compared to these published work, our new deep CO(4–3) data allow us to improve the constraint on the CO redshift, flux density, and line width, which is a vital parameter for the determination of dynamic mass. In addition, we put useful constraints on the 3 mm continuum emission for these SMGs. Table 1 summarizes the observed and derived CO properties for these three galaxies.

To measure the velocity centroid and line width of CO emission of GN20, we fitted Gaussians to the observed spectra (Fig. 1), which is extracted at the fixed CO position from Daddi et al. (2009b), allowing for the presence of faint underlying continuum. The CO(4–3) spectrum of GN20 shows clear evidence of double-peaked structure. A single Gaussian fit to the spectrum shown in Fig. 1 yields a peak flux density of 2.13 ± 0.14 mJy and an FWHM of 630 ± 51 km s^-1, consistent with Daddi et al. (2009b). We also fitted the spectrum with a double Gaussian function by fixing the FWHM in each component to the same value. The peak flux densities for the two components are 2.33 ± 0.17 mJy and 2.02 ± 0.16 mJy, respectively. We find that the total FWHM (583 ± 36 km s^-1) derived from a double Gaussian fitting is better constrained than that from a single one. The redshift derived from a double Gaussian is 4.0553 ± 0.0002, consistent with the previous measurements (Daddi et al. 2009b; Carilli et al. 2011; Hodge et al. 2012). A circular Gaussian model fit with a fixed FWHM of 0.72′′ (we here adopt the FWHM of the CO(6–5) image measured by Carilli et al. 2010) to the CO(4–3) image of GN20 (Fig. 2) yields a velocity-integrated flux of I_CO = 1.68 ± 0.10 Jy km s^-1, in line with the measurements in Daddi et al. (2009b).

For GN20.2a, the fits with a single Gaussian and two Gaussian functions yield FWHM of 820 ± 237 km s^-1 and 763 ± 180 km s^-1, respectively. The spectrum is extracted at the fixed CO position from Daddi et al. (2009b). We find the CO spectrum of GN20.2a appears to be better described by a double Gaussian profile (Fig. 1), given the relatively low uncertainties of FWHM derived from double Gaussian fitting and the rapid decrease of the flux at the edges of the spectrum, though there is no clear indication of double-peaked emission in the spectrum. The favored fit with two Gaussian functions could in principle suggests either the existence of a rotating component, or kinematically distinct components undergoing a merger within the system. We will further discuss the possibility in the following section. The two components in the double Gaussian fitting have peak flux densities of 0.82 ± 0.21 mJy and 0.69 ± 0.22 mJy, respectively. A Gaussian fit to the CO spectrum gives a redshift of 4.0508 ± 0.0013. We used a circular Gaussian model with a fixed FWHM of 0.53′′ (derived from CO(2–1) image in Hodge et al. 2013) to fit the CO(4–3) map and derive a velocity-integrated flux of 0.65 ± 0.08 Jy km s^-1. All the measurements derived for GN20.2a are consistent with published results (Daddi et al. 2009b; Carilli et al. 2011; Hodge et al. 2013), but exhibit significantly improved sensitivity and signal-to-noise ratio (S/N).

The CO(4–3) spectrum of GN20.2b has been presented in Hodge et al. (2013). Compared to the large uncertainty in the line width fitted in Hodge et al. (2013), our deeper CO(4–3) data provide a much better constraint on the estimate of the line width. Similarly, we performed Gaussian fits to the spectrum (Fig. 1) extracted at the fixed position of 1.2 mm emission in GN20.2b (Riechers et al., in prep.) with a single Gaussian and double Gaussian functions, respectively. Both Gaussian fits give a similar FWHM with value of ~220 ± 40 km s^-1, which is much smaller than GN20 and GN20.2a, and also the SMGs at z ~ 2−4 with mean FWHM of 470 ± 80 km s^-1(Bothwell et al. 2013). The Gaussian fit to the spectrum has a peak flux density of 1.17 ± 0.19 mJy and a redshift of 4.0563 ± 0.0003. The velocity-integrated flux density derived from the CO(4–3) map with a circular Gaussian model (fixed FWHM of 0.88′′; see Hodge et al. 2013) is 0.27 ± 0.04 Jy km s^-1.

We calculate the CO(4–3) line luminosities (in K km s^-1 pc²) using the standard relation given by Solomon et al. (1997): $\begin{array}{ccc}{\mathit{L}}_{\mathrm{CO}}^{\mathrm{\prime }}\mathrm{=}\mathrm{3.25}\mathrm{\times }{\mathrm{10}}^{\mathrm{7}}{\mathit{S}}_{\mathrm{CO}}\mathrm{\Delta }\mathit{V}{\mathit{\nu }}_{\mathrm{obs}}^{-2}{\mathit{D}}_{\mathrm{L}}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{-3}& & \end{array}$(1)where S_COΔV is the velocity-integrated line flux in Jy km s^-1, ν_obs is the observed frequency in GHz and D_L is the luminosity distance in Mpc. With ${\mathit{L}}_{\mathrm{CO}\mathrm{(}\mathrm{2}\mathrm{-}\mathrm{1}\mathrm{)}}^{\mathrm{\prime }}$ measured by Carilli et al. (2011), we find CO(4–3)/CO(2–1) line brightness temperature ratios of r_{43 / 21} = 0.41, 0.36, and 0.44 for GN20, GN20.2a, and GN20.2b, respectively. These are consistent with the mean ratio (~0.48 ± 0.10) measured for SMGs at z ~ 2−4 (Bothwell et al. 2013).

3.2. The IR properties

3.2.1. Millimeter continuum emission

Figure 1 shows evidence for 3.3 mm continuum emission. Averaging the spectrum in the uv space for velocities outside of CO(4–3) line emission range, we find a bright source at the position of GN20, and 3.3 mm continuum detections with significance of 2.5σ and 3.3σ for GN20.2a and GN20.2b (Fig. 2), respectively. The flux densities are measured by fitting Gaussian functions to the maps using sizes estimated from the CO imaging and are summarized in Table 2. The 3.3 mm flux density of GN20 is consistent with the lower S/N measurement reported in Daddi et al. (2009b), while for GN20.2a and GN20.2b, they are presented here for the first time.

In addition, we detect 1.2 mm emission of GN20, GN20.2a, and GN20.2b in PdBI maps with significant S/N level (see Riechers et al., in prep.). Measured fluxes are in the range of 3.3–8.5 mJy (Table 2) and are used in the remainder of the paper. We also make use of 2.2 mm continuum data of these three galaxies, which have been reported in Carilli et al. (2010). Similar to the measure of CO(4–3) flux density, we used a circular Gaussian model with a fixed FWHM, same as the one used for CO(4–3) map, to fit the 1.2 and 2.2 mm images for each galaxy. For a consistency check, assuming an GN20-like spectral energy distribution (SED; Magdis et al. 2011), we extrapolate the flux density of 1.2 mm to 2.2 mm and 3.3 mm for GN20.2a and GN20.2b, and find these predictions match with the 2.2 mm and 3.3 mm flux densities derived from PdBI maps. The derived mm flux densities for these galaxies are summarized in Table 2.

3.2.2. Dust masses, dust temperatures and SFRs

Because they sample the rest-frame far-IR/submm bands at the Rayleigh-Jeans tail, the mm data have been demonstrated to have significant impact on the derivation of far-IR properties of high redshift galaxies (Magdis et al. 2012a), particularly crucial to reduce the uncertainties of dust mass estimates. Combining the Spitzer IRAC, MIPS photometry (Daddi et al. 2009b) with the Herschel PACS 100 μm, 160 μm and SPIRE 250 μm photometry (for GN20.2a and GN20.2b, the 250 μm data were not used for SED fitting because the 250 μm photometry is blended between these two galaxies) and following Magdis et al. (2012a), we derive the dust mass of the galaxies by fitting the IR SED using the Draine & Li (2007, hereafter DL07) dust models. The DL07 model describes the interstellar dust as a mixture of carbonaceous and amorphous silicate grains, with a size distribution mimicking the Milky Way extinction curve. The fraction of dust mass in the form of polycyclic aromatic hydrocarbon (PAH) grains (with <10³ carbon atoms) is parametrized by the PAH index, q_PAH. The PAH abundance has an important effect in shaping the galaxy SED at short wavelengths. According to the DL07 model, the spectrum of a galaxy can be described by a linear combination of one stellar component approximated by a blackbody with a color temperature of 5000 K and two dust components (see Magdis et al. 2012a). The majority (i.e., 1 − γ with 0.0 <γ< 0.3) of the dust grains are located in the diffuse ISM and heated by a radiation field contributed by many stars with a minimum intensity U_min (U is a dimensionless factor normalized to the local ISM), and the rest dust grains are exposed to a power-law radiation field ranging from U_min to U_max, associated with photodissociation regions (PDRs). Following the prescriptions of DL07, we fit the mid-IR to millimeter data points for each galaxy. For GN20, we use the updated 3.3 mm flux density for SED fitting and find a consistent dust mass with the estimates in Magdis et al. (2011, 2012a). In order to more accurately derive far-IR properties, we also add 1.4 GHz measurements (Table 2; Owen et al., in prep.) to the SED fitting, with the exception of GN20.2a, due to the presence of an AGN (Daddi et al. 2009b). The 1.4 GHz flux density of GN20.2b measured from our new deep radio map is revised to a factor of 2 lower than the value reported previously (Daddi et al. 2009b; Morrison et al. 2010). The best-fit models along with the observed SEDs for GN20, GN20.2a, and GN20.2b are shown in Fig. 3, and the derived parameters are summarized in Table 2.

The radiation field intensities ⟨ U ⟩ derived from the best-fit models are 27.2$\begin{array}{c}\mathrm{+}\mathrm{2.6}\\ -2.2\end{array}$, 26.4$\begin{array}{c}\mathrm{+}\mathrm{8.1}\\ -14.1\end{array}$, and 25.8$\begin{array}{c}\mathrm{+}\mathrm{6.3}\\ -12.5\end{array}$ for GN20, GN20.2a and GN20.2b, respectively. These can be represented as dust temperatures if assuming that the dust is heated only by radiation. Integrating over the best-fitting template from 8 to 1000 μm, we derive L_IR = (1.86 ± 0.09) × 10¹³L_⊙, (8.0 ± 0.7) × 10¹²L_⊙, and (6.9 ± 1.0) × 10¹²L_⊙ for GN20, GN20.2a, and GN20.2b, respectively. We find that the uncertainties of IR luminosity estimates are significantly reduced by adding 1.4 GHz measurements. To convert the observed IR luminosities into SFRs, we adopt a Kennicutt (1998) conversion of SFR[M_⊙ yr^-1 ] = L_IR[L_⊙]/10¹⁰, appropriate for a Chabrier (2003) IMF. The SFR derived in this way is used for analysis in the remainder of the paper.

3.2.3. The Radio-IR correlation

The IR-radio correlation is usually quantified by the so-called q-parameter (Helou et al. 1985): $\begin{array}{ccc}\mathit{q}\mathrm{=}\log \left(\frac{{\mathit{L}}_{\mathrm{IR}}}{\mathrm{3.75}\mathrm{\times }{\mathrm{10}}^{\mathrm{12}}\mathrm{W}}\right)\mathrm{-}\log \hspace{0.17em}\left(\frac{{\mathit{L}}_{\mathrm{1.4}\mathrm{GHz}}}{\mathrm{W}\mathrm{H}{\mathrm{z}}^{-1}}\right)\mathrm{·}& & \end{array}$(2)where L_{1.4 GHz} is the k-corrected radio luminosity assuming S_ν ∝ ν^α with α = − 0.8. Using the 1.4 GHz flux density measurements (Table 2; Owen et al. in prep.) and the IR luminosities derived above, we find values of q = 2.41 ± 0.07, 1.72 ± 0.05, and 2.60 ± 0.19 for GN20, GN20.2a, and GN20.2b, respectively. The q values of GN20 and GN20.2b are compatible with the local relation (q = 2.64 ± 0.02, scatter: 0.26; Bell 2003) and galaxies at z ~1–2 (Sargent et al. 2010), suggesting that dust heating in these two galaxies originates predominantly from star formation. The mid-IR observations of GN20 reveal a significant power-law component, likely related to the presence of an obscured AGN (Riechers et al. 2014b). However, the infrared SED fits to GN20 after correcting for the AGN fraction show that the AGN contribution to the total IR luminosity is only ≲15% (Riechers et al. 2014b). In contrast, GN20.2a falls below the range defined for star-forming galaxies (Yun et al. 2001; Bell 2003), the relatively low q value is suggestive of AGN activity boosting radio flux. In addition, we found that the q value of GN20.2a is similar to that of the Cloverleaf quasar at z = 2.55, which hosts a known radio-loud AGN (Beelen et al. 2006; Kayser et al. 1990).

Fig. 3 Observed mid-IR to millimeter SED of GN20, GN20.2a, and GN20.2b, overlaid with the best-fit Draine & Li (2007) models. The black line is the DL07 model consisting of one stellar component and two dust components: diffuse ISM and photodissociation region. Green dot-dashed and cyan dot-dashed lines show stellar component and “PDR” component, respectively. The best-fit parameters from DL07 models are listed within each panel. The arrows represent 5σ upper limits.

3.3. Multiwavelength counterparts

Figure 4 shows the CO(6–5) contours of GN20, GN20.2a, and GN20.2b with resolution of 0.9′′/6.1 kpc overlaid on the HST+ACS 850z-band, the HST+WFC3 F160W-band, the CHFT+WIRCam K-band, and the Spitzer+IRAC 3.6 μm images. The WFC3 images are publicly available from the CANDELS survey (Grogin et al. 2011; Koekemoer et al. 2011). For GN20.2b, the multiwavelength counterparts are found to be coincident with the CO source. We note that a nearby companion galaxy lies ~2′′ to the west of GN20.2b, which is also likely to be a B-band dropout based on the ACS images shown in Daddi et al. (2009b). For GN20.2a, we find a significant offset between CO position and optical counterpart, which has already been presented in Hodge et al. (2013) with higher resolution of CO(2–1) map. This is similar to the large offset found between CO and optical counterpart for GN20 (Carilli et al. 2010), revealing a substantial dust obscuration over a scale ~10 kpc. The spectroscopic redshift (z = 4.059 ± 0.007; Daddi et al. 2009) derived for the optical counterpart is found to be very close to the CO redshift of GN20.2a, makes it likely that the HST galaxy to the northeast of the radio position is related to the CO emission (see Fig. 4).

Fig. 4 PdBI CO(6–5)+continuum contours (Carilli et al. 2010) at resolution of 0.90′′ × 0.76′′ for GN20 (inset), GN20.2a and GN20.2b overlaid on the ACS 850z-band (top left), WFC3 F160W-band (top right), WIRCam Ks-band (bottom left), and IRAC 3.6 μm (bottom right) images. GN20 lies about 24′′ to the northeast of GN20.2a. Black and red crosses show the positions of the 1.4 GHz counterparts (Owen et al., in prep.) and optical counterparts of GN20, GN20.2a and GN20.2b, respectively. Contour levels are shown in steps of 1σ for GN20.2a and GN20.2b and 3σ for GN20 starting at ±3σ. The beam sizes of CO(6–5) are displayed in the left corner of each map.

3.4. Stellar mass and UV-based SFR estimates

The photometric information we use to measure the physical properties of our galaxy sample is drawn from the GOODS-N multiwavelength catalog presented in Pannella et al. (2014). For the sake of clarity we give here only a brief summary of the catalog properties while referring the reader to Pannella et al. (2014) for a more detailed description. The extended GOODS-N catalog is a Ks-band selected multiwavelength catalog spanning 20 passbands from GALEX NUV to IRAC 8 μm. PSF-matching corrections have been applied to account for the different angular resolution of the images. Aperture magnitudes (2′′ diameter) are used to sample the galaxy SED. Finally, all derived properties were extrapolated to “total” masses using the ratio between the total (FLUX_AUTO) and aperture flux in the K-band detection image. The catalog contains 53 675 objects over the WIRCAM Ks image field of 900 arcmin² and down to an AB magnitude of 24.5, the 5σ limiting magnitude of the image. The UBVRIzJHK photometry, plus the photometry in IRAC bands are shown for GN20, GN20.2a, and GN20.2b in Fig. 5.

The stellar masses were estimated by fitting the multiwavelength photometry of U- to IRAC 4.5 μm band to Bruzual & Charlot (2003) templates using FAST (Fitting and Assessment of Synthetic Templates; Kriek et al. 2009) through a χ² minimization. We fix the spectroscopic redshift as given by the CO line identification. The IRAC 3.6 μm measurements were not used in the SED fitting, to avoid the possible flux contamination from Hα emission, which happened to fall into this band. Given that all our galaxies are actively star forming, we adopt the templates with a constant star formation history (SFH), a large range of ages, a Chabrier (2003) IMF, and allowing the metallicity to vary within a range of Z = 0.2–2.5 solar. We used an extinction correction from Calzetti et al. (2000) with an optical extinction A_V = 0 − 6 mag. The typical uncertainty of stellar mass estimate is 0.2 dex. For comparison, we also assume an exponentially declining SFH for the modeling and find the stellar masses are on average increased by 0.05 dex for each galaxy, while the UV-based SFRs are lowered by ~0.5 dex. Therefore, the uncertainties of stellar mass estimates resulted from different assumptions of single-component SFH are relatively small for our sources. We also find that the stellar masses derived from the best-fitting model allowing a varied metallicity are not significantly different from the one assuming a solar metallicity. We adopt the modeling results derived by assuming a constant SFH for analysis in the remainder of the paper.

Fig. 5 Observed optical to mid-IR SED of GN20 (left), GN20.2a (middle), and GN20.2b (right). Red circles with error bars represent the multiwavelength photometry. The solid lines show the best-fitting templates from the libraries of Bruzual & Charlot (2003) at the derived CO redshift for each galaxy, assuming a constant SFH models with solar metallicity and using an extinction correction from Calzetti et al. (2000). We adopt the templates with a constant SFH to derive stellar mass by fitting the photometry of U- to IRAC 4.5 μm band.

For GN20.2b, the stellar masses derived from the SED fitting up to K-band and to IRAC 4.5 μm band are differing by 0.2 dex. Given that the IRAC photometry for GN20.2b is likely to be blended with the companion source (see Fig. 4), we adopt the average of these two estimates, M_⋆ = 1.1 × 10¹¹ M_⊙, and an uncertainty of about 0.2 dex. Looking at the IRAC image of GN20.2a, we find that the IRAC peak is coincident with the northeast optical counterpart and close to the CO peak, and thus we take into account the IRAC 4.5 μm photometry in the SED fitting for the stellar mass derivation, yielding M_⋆ = 3.8 × 10¹⁰ M_⊙. We stress that the stellar mass estimate of GN20.2a should be treated with caution, due to the possible contamination in the IRAC band emission from its nearby companion galaxy (see Fig. 4), and the possibility that the optical “counterpart” is a distinct unrelated galaxy from the CO-emitting dusty starburst galaxy. It is not easy to disentangle this issue given the current available data.

We estimated the UV-based SFR with the same SED fitting approach used to derive stellar masses. However, also to avoid possible uncertainties linked to blending in the IRAC bands, we only consider the photometry of U- to K-bands for the measure of UV-based SFR, as the observed K band corresponds to the rest-frame 4000 Å Balmer break. Comparing the SFR estimated from dust-corrected UV luminosity with the one converted from IR luminosity, we find that the IR-based SFR estimates are generally larger than the UV-based ones for these three SMGs. The ratios of SFR_IR/SFR_UV we derived are 13, 26, and 6 for GN20, GN20.2a, and GN20.2b, respectively. As discussed above, the UV-based SFRs would be lowered by a factor of ~3 if assuming an exponentially declining SFH, thus result in even higher SFR_IR/SFR_UV ratios. The extremely large value of SFR_IR/SFR_UV measured in GN20.2a is indicative of heavily obscured UV/optical emission at the positions of CO/radio emission, in agreement with the large offset observed between optical and CO peaks shown in Fig. 4. As discussed in Daddi et al. (2009b), both the complex morphology and large ratio of IR/UV-derived SFRs may suggest that GN20.2a is undergoing a major merger with a B-band dropout companion (~0.7′′ west to GN20.2a; see Fig. 4).

Fig. 6 Global UV to far-IR fits to our galaxies based on MAGPHYS. The black line is the best-fit SED and the blue line is the corresponding unattenuated stellar emission.

3.5. Physical parameters from global UV to far-IR SED fitting

In order to verify the stability of our results and test against possible systematic effects affecting our modeling, we obtained simultaneous estimates of the stellar masses, SFRs, IR luminosities and dust masses using MAGPHYS (da Cunha et al. 2008). We use an updated version of MAGPHYS where the parameter priors are appropriate for high-redshift galaxies (da Cunha et al., in prep.). This code includes a wide range of possible SFHs and dust properties, and fits the full SED consistently from the UV to the far-IR by requiring an energy balance between the UV/optical and far-IR output (Fig. 6).

We find in general excellent agreement with the results obtained from separate UV/optical and IR fitting described in previous sections. In all cases the SFRs are very close to what inferred from the IR luminosities applying a standard Kennicutt (1998) conversion. After correcting for the different κ values (κ is the dust absorption coefficient) adopted in MAGPHYS and in the DL07 model, the differences of dust mass estimates for these three sources are within 25–30%. For the stellar mass, except GN20.2b, the differences are also only 20%.

For GN20.2b, MAGPHYS provide two degenerate solutions for the stellar mass. One is consistent with our previous derivation (about 10¹¹ M_⊙) and a second one with much smaller values up to factors of 3–5.

3.6. The CO-to-H₂ conversion factor

Recently, a few attempts at the direct measurements of the CO luminosity to H₂ gas mass conversion factor for high redshift galaxies have been mainly based on dynamical modeling, gas-to-dust measurements, and radiative transfer modeling (e.g., Bolatto et al. 2013; Carilli & Walter 2013, and references therein). A pilot study performed by Magdis et al. (2011) showed that the gas-to-dust ratio method applied to GN20 leads to a conversion factor of α_CO< 1.0, consistent with the dynamical estimate measured by Carilli et al. (2010). Having derived the dust, stellar, and dynamical masses, we use two independent methods of dynamical modeling and gas-to-dust mass ratio to estimate α_CO. To avoid the uncertainties resulting from the extrapolation from high-order CO transitions, we adopt the low-J (J = 1,2) CO luminosities measured by Carilli et al. (2010, 2011). The CO luminosities are summarized in Table 3.

3.6.1. The dynamical mass method

Assuming that our SMGs are a mix of disc-like and virialized systems, we estimate the dynamical mass by taking the average of two different estimators: the isotropic virial estimator (e.g., Spitzer 1987; Pettini et al. 2001) $\begin{array}{ccc}{\mathit{M}}_{\mathrm{dyn}}\mathrm{(}\mathit{r}\mathit{<}{\mathit{r}}_{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{)}\mathrm{=}\frac{\mathrm{5}{\mathit{\sigma }}^{\mathrm{2}}{\mathit{r}}_{\mathrm{1}\mathit{/}\mathrm{2}}}{\mathit{G}}& & \end{array}$(3)where r_{1 / 2} is the half-light radius, σ = Δυ_FWHM/2.35 is the one-dimensional velocity dispersion (Δυ_FWHM is the line width), and G is the gravitational constant, and the rotating disk estimator (Neri et al. 2003), corrected for ⟨ sin²i ⟩ = 2 / 3 in mass: $\begin{array}{ccc}{\mathit{M}}_{\mathrm{dyn}}\mathrm{(}\mathit{r}\mathit{<}{\mathit{r}}_{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{)}\mathrm{=}\mathrm{6}\mathrm{\times }{\mathrm{10}}^{\mathrm{4}}\mathrm{\Delta }{\mathit{\upsilon }}_{\mathrm{FWHM}}^{\mathrm{2}}{\mathit{r}}_{\mathrm{1}\mathit{/}\mathrm{2}}\mathit{.}& & \end{array}$(4)This method has been widely applied to the dynamical mass estimate for SMGs (e.g., Tacconi et al. 2008; Bothwell et al. 2013; Hodge et al. 2013). Unlike the dynamical mass traced by HI emission line in local universe, the CO-based dynamical mass represents more compact star formation region (e.g., Solomon & Vanden Bout 2005; de Blok & Walter 2014). Adopting the FWHM of CO(4–3) line derived from our best-fitting double Gaussian models and the CO size measured from resolved CO(2–1) maps (Carilli et al. 2010; Hodge et al. 2013), the implied dynamical mass of GN20 within r_{1 / 2} ~ 4 kpc is (1.8 ± 0.2) × 10¹¹ M_⊙, of GN20.2a within r_{1 / 2} ~ 2.5 kpc is (1.9 ± 0.9) × 10¹¹ M_⊙, and of GN20.2b within r_{1 / 2} ~ 4 kpc is (2.6 ± 1.0) × 10¹⁰ M_⊙. The uncertainties were estimated based on the uncertainties on the FWHM measures. The dynamical mass estimates for GN20 and GN20.2a are consistent within uncertainties with the estimates by Carilli et al. (2010) and Hodge et al. (2012, 2013). For GN20.2b, our dynamical mass estimate is better constrained than the one measured by Hodge et al. (2013), as the uncertainty in its fitted line width is significantly reduced.

To estimate the gas mass from the dynamical modeling, we follow the method in Daddi et al. (2010)$\begin{array}{ccc}\mathit{M}\mathrm{(}\mathit{r}\mathit{<}{\mathit{r}}_{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{)}\mathrm{=}\mathrm{0.5}\mathrm{\times }\mathrm{(}{\mathit{M}}_{\mathit{\star }}\mathrm{+}{\mathit{M}}_{\mathrm{gas}}\mathrm{)}\mathrm{+}{\mathit{M}}_{\mathrm{dark}}\mathrm{(}\mathit{r}\mathit{<}{\mathit{r}}_{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{)}& & \end{array}$(5)where the dynamical mass within the half-light radius is composed of the half of the total amount of stellar mass and gas mass, and the amount of dark matter within r_{1 / 2}. Subtracting the stellar mass, including a 25% dark matter contribution, which is a typical value adopted for z ~ 1–2 normal galaxies (Daddi et al. 2010) and z ~ 2 SMGs (Tacconi et al. 2008), we find α_CO = 1.0 ± 0.2 and 2.8 ± 1.6 for GN20 and GN20.2a, respectively. The derived α_CO for GN20 is consistent with previous estimates (Carilli et al. 2010; Hodge et al. 2012). For GN20.2b, the dynamical mass is found to be much smaller than our preferred estimate of the stellar mass, suggesting that our method may possibly significantly underestimate the dynamical mass. If we adopt the rotating disk estimator and make an assumption for the extreme case that GN20.2b is composed of stellar mass and dark matter only, the maximum inclination angle derived for this galaxy would be about 20°, close to face-on. This is consistent with the narrow CO line width that we observed. On the other hand, the alternative lower stellar mass value suggested by MAGPHYS for this object is compatible with the inferred dynamical mass. Nevertheless, given all the uncertainties, no inference on dynamical mass and α_CO for GN20.2b can be directly obtained in this way.

3.6.2. The gas-to-dust ratio method

Some recent studies have revealed a tight correlation between the gas-to-dust ratio (δ_GDR) and the gas-phase oxygen abundance, with δ_GDR decreasing for more metal-rich galaxies (e.g., Leroy et al. 2011; Magdis et al. 2011). Under the assumption that the local observed δ_GDR − Z relation is valid at high redshift, Magdis et al. (2012a) derived indirect estimates of α_CO for a sample of star-forming galaxies at z ~ 2 and found that the variation tendency of α_CO is in line with previous studies (Leroy et al. 2011). Following the prescription in Magdis et al. (2012a), we attempt to apply this approach to the three galaxies in our sample to estimate α_CO.

Having obtained the dust masses for our sources, we need the measurements of their metallicities, for which we have to rely on indirect indicators. The fundamental metallicity relation (FMR) of Mannucci et al. (2010) that relates the metallicity to both SFR and M_⋆ is one of the methods that is commonly used for measuring metallicity. However, the FMR relation is valid only up to z ~ 2.5, and thus we also apply the mass-metallicity relation of present-day elliptical galaxies (e.g., Calura et al. 2009), assuming that the high redshift strongly star-forming galaxies are the progenitors of the local ellipticals, given that the large SFR (~690–1860 M_⊙ yr^-1) of our sources could be due to a final burst of star formation triggered by a major merger, and the galaxy will eventually evolve into a massive elliptical with unchanged mass and metallicity once there are no further star formation activities. The metallicity estimates based on this scenario and the FMR relation lead to 12 +log [O/H] = 8.8–9.2 for GN20 (agree with the values in Magdis et al. 2011), 8.6–9.0 for GN20.2a, and 8.9–9.2 for GN20.2b.

With the derived metallicity, we estimate the δ_GDR following the relation of Magdis et al. (2012a) which relates the δ_GDR to metallicity, i.e., $\begin{array}{ccc}\log {\mathit{\delta }}_{\mathrm{GDR}}\mathrm{=}\mathrm{(}\mathrm{10.54}\mathrm{\pm }\mathrm{1.0}\mathrm{)}\mathrm{-}\mathrm{(}\mathrm{0.99}\mathrm{\pm }\mathrm{0.12}\mathrm{)}\mathrm{\times }\mathrm{\left(}\mathrm{12}\mathrm{+}\mathrm{1}\mathrm{og}\mathrm{[}\mathrm{O}\mathit{/}\mathrm{H}\mathrm{]}\mathrm{\right)}\mathit{.}& & \end{array}$(6)The fit to this relation indicates a value for δ_GDR of ~77(27) for GN20, ~106(42) for GN20.2a, and ~53(27) for GN20.2b, assuming a metallicity estimated based on the FMR relation (mass-metallicity relation of present-day elliptical galaxies). With the gas-to-dust ratios, we can determine the M_gas and estimate α_CO from the equation ${\mathit{\alpha }}_{\mathrm{CO}}\mathrm{=}{\mathit{L}}_{\mathrm{CO}}^{\mathrm{\prime }}\mathit{/}{\mathit{M}}_{\mathrm{gas}}$. Here we have assumed that at high redshift M_H2 ≫ M_HI and therefore that M_gas ⋍ M_H2. This assumption is based on the observed high surface densities, above the characteristic threshold for HI saturation and thus most of the gas is in molecular form (e.g., Bigiel et al. 2008; Obreschkow & Rawlings 2009). The α_CO derived based on gas-to-dust method are 1.6 ± 0.7, 2.0 ± 1.1, and 2.8 ± 1.5 for GN20, GN20.2a, and GN20.2b, respectively. We note that our α_CO estimate for GN20 is larger than previously found by Magdis et al. (2011) with the same method, which is due to the revised estimates of dust and stellar mass used in our study. We caution that these α_CO estimates are systematically uncertain due to the large uncertainties of metallicity estimates.

3.6.3. Comparison of the CO-to-H₂ conversion factor

Comparing the α_CO estimates summarized in Table 3, we find that the α_CO determined based on the metallicity-dependent gas-to-dust ratio method are consistent within uncertainties with those measured from the dynamical modeling for each galaxy. The average α_CO are 1.3 ± 0.4, 2.4 ± 1.3, and 2.8 ± 1.5 for GN20, GN20.2a, and GN20.2b, respectively, which are found to be between the typical value (~0.8) determined for ULIRGs and the value (~4) appropriate for the Milky Way. However, these values might be well below the value appropriate for normal galaxies at those epochs, as the α_CO of z> 4 normal galaxies could be much higher than the Milky Way value (e.g., Genzel et al. 2012; Bolatto et al. 2013; Carilli & Walter 2013; Tan et al. 2013). We caution that the uncertainty of α_CO estimates derived from each method is significant. For the dynamical analysis, the main uncertainties are the crude estimates of dynamical mass based on marginally resolved imaging data and the stellar mass based on SED fitting, which is caused by the extreme dust obscuration and the potentially complex SFHs, while for the gas-to-dust method, the large uncertainties might result from the dust model assumed and the uncertainty on metallicity estimate. The consistency of α_CO estimates derived from two independent approaches further confirm the reliability of the gas-to-dust ratio method.

3.7. The specific star formation rates and star formation efficiency

With the stellar masses and SFR derived above, we estimate sSFR of 16.9 Gyr^-1, 20.9 Gyr^-1, and 6.3 Gyr^-1 for GN20, GN20.2a, and GN20.2b, respectively. Both GN20 and GN20.2a exhibit larger sSFRs with a factor of ~6 than the average of equally massive galaxies on the SFR-M_⋆ main-sequence at z ~ 4 (Sargent et al. 2014), while GN20.2b shows a sSFR-excess of 2.4 (which could be larger though, in case of the lower-mass MAGPHYS solution is adopted instead). The large sSFR-excess for GN20 and GN20.2a, similar to the z ~ 2−3 SMGs, have been revealed in Daddi et al. (2009b), suggesting that these galaxies are undergoing starburst events with short duty cycles, likely triggered by mergers.

The ratio between L_IR and ${\mathit{L}}_{\mathrm{CO}}^{\mathrm{\prime }}$ (i.e., SFE, in units of L_⊙ (K km s^-1 pc²)^-1) can be used as a measure of the efficiency with which molecular gas is converted into stars. Given that the total available reservoir of gas within a galaxy could be better traced by low-J CO emission, we use directly observed CO(1–0) luminosity or the one extrapolated from CO(2–1) observations for SFE derivation to avoid the uncertainties in the underlying gas excitation. We find ratios of ~90–238 for these three galaxies, comparable to those derived for other high-redshift SMGs (Ivison et al. 2011; Bothwell et al. 2013), and higher than in normal star-forming galaxies (~84 ± 12; Daddi et al. 2010), but slightly lower than in local ULIRGs (~250 ± 30; Solomon et al. 1997).

3.8. Gas fractions, depletion timescales, and the evolutionary status

Table 3 lists the molecular gas masses for our three galaxies, derived based on the average α_CO inferred in Sect. 3.6 and the measured ${\mathit{L}}_{\mathrm{CO}}^{\mathrm{\prime }}$. Combining the gas mass with the stellar mass estimates, we find molecular gas fraction of f_gas = M_gas/ (M_gas + M_⋆) = 0.66 ± 0.12, 0.85 ± 0.09, and 0.42 ± 0.18 for GN20, GN20.2a, and GN20.2b, respectively. These values are comparable to other SMGs at z = 2−4 (Bothwell et al. 2013) and high-z massive, gas-rich star-forming galaxies (e.g., Daddi et al. 2010; Tacconi et al. 2010). Bothwell et al. (2013) compared the gas fraction of SMGs to that of local LIRGs and found an increase of gas fraction up to z ~ 2 followed by a flattening toward higher redshift. This is similar to the redshift evolution trend of gas fraction revealed in normal star-forming galaxies (Magdis et al. 2012b; Tan et al. 2013), despite that SMGs have typically an order of magnitude higher SFRs than normal galaxies. The 2 star formation mode (2-SFM) predictions by Sargent et al. (2014) have also revealed that the average gas fraction of starbursts is only slightly lower than that of main sequence galaxies.

Dividing the stellar masses by the SFRs, we estimate the stellar mass building-up timescale of τ_build ~ 59 Myr, 48 Myr, and 159 Myr for GN20, GN20.2a, and GN20.2b, respectively, if assuming that these galaxies have sustained their current SFRs continuously. Similarly, assuming that the SFRs continues at the current rate and neglecting the effect of feedback, the minimum time for exhausting the molecular gas reservoir can be given by the gas depletion timescale, τ_gas = M(H₂)/SFR, which we find to be ~113 Myr, ~275 Myr, and ~116 Myr for GN20, GN20.2a, and GN20.2b, respectively. These timescales are relatively short than in massive, normal galaxies (~0.4–0.9 Gyr; e.g., Daddi et al. 2010; Tacconi et al. 2010), indicating a rapid star formation mode with intense burst for SMGs, probably undergoing major mergers or rapid cold accretion (e.g., Narayanan et al. 2010; Davé et al. 2010). In addition, the extremely high gas surface density revealed in GN20.2a may indicate a triggering mechanism of major merger, while the surface density of GN20.2b is found to be comparable to normal star-forming galaxies (Hodge et al. 2013). If we define starburst phase to be the fraction of stars have already been formed over the available gas reservoir, simply calculated as M_⋆/ (M_gas + M_⋆), taking these numbers at face value we would find that the relatively young stellar age and long gas-consumption timescale of GN20.2a place it at an early stage, ~15% ± 9% of the way through its starburst, while GN20 and GN20.2b have experienced ~35% ± 14% and ~60% ± 35% of their starburst phase, respectively. The CO excitation analysis has also revealed different merging states for these three galaxies, of which GN20.2b display lowest excitation (Carilli et al. 2010; Hodge et al. 2013).

4. What is the nature of GN20, GN20.2a and GN20.2b?

We have presented multiwavelength properties of GN20, GN20.2a and GN20.2b, including UV/optical, far-IR, and mm photometry and molecular gas content. Table 4 lists the ratios of sSFR/sSFR_MS and SFR_IR/SFR_UV for these three SMGs. The large sSFR-excess observed for GN20 and GN20.2a suggest that these two galaxies are starbursting outliers above the main sequence (Rodighiero et al. 2011; Sargent et al. 2014). Although GN20.2b might situate within the MS scatter with sSFR/sSFR_MS ~ 2.4, the large value of SFR_IR/SFR_UV indicates that a large fraction of UV emission from this galaxy is obscured, pointing also toward a population of starburst galaxy (Daddi et al. 2007b, 2010). For GN20 and GN20.2a, the large offsets between the CO positions and the optical counterparts show clear evidence of extreme obscuration, consistent with the heavy dust extinction observed in the UV emission. In addition, the high gas surface densities derived for GN20 and GN20.2a suggest that the star formation is dominated by the compact star-forming sites (Table 4; Carilli et al. 2010; Hodge et al. 2012, 2013). As discussed in Sect. 3.8, the gas-consumption timescales (see Table 4) of these three SMGs are found to be significantly shorter than those of normal galaxies, indicative of a more rapid star formation mode. The extreme SFRs (≳700 M_⊙ yr^-1) and high radiation field intensity (⟨ U ⟩ > 25) observed for our SMGs would also favor the scenario of intense starbursts, as suggested by Magdis et al. (2011) (see also Magnelli et al. 2012a). In addition, our three SMGs appear to be dynamically distinct from the normal disk galaxies, since these galaxies are clearly separated from disk galaxies in the velocity-size plane, which can be used as a tool to constrain the angular momentum properties of galaxies (Courteau 1997; Bouché et al. 2007). The low orbital angular momentum observed in our SMGs might be caused by a recent or ongoing merger.

Two main modes are typically considered for triggering star formation in high redshift SMGs: major mergers (Narayanan et al. 2010) and secular mode with smooth mass infall from the intergalactic medium and along the cosmic web (Dekel et al. 2009; Davé et al. 2010). However, it is still a matter of debate to identify which mode of star formation is responsible for the majority of the SMG population. High-resolution CO imaging of local ULIRGs and z ~ 2 SMGs have revealed large reservoirs of molecular gas concentrated in the galaxy nuclei with disturbed CO kinematics/morphologies, suggesting a merger-driven mechanism for the star formation activity (Downes & Solomon 1998; Tacconi et al. 2008; Engel et al. 2010). However, the resolved CO imaging of GN20 shows ordered rotation and an extended gas distribution (Carilli et al. 2010; Hodge et al. 2012), which could be suggestive of a rotating disk. On the other hand, hydrodynamical simulations of gas-rich mergers have shown that the merger remnant gas can cool quickly and produce an extended star-forming disk (Springel & Hernquist 2005; Robertson et al. 2006; Robertson & Bullock 2008).

GN20.2a has a smaller CO size than GN20 and GN20.2b, but exhibits a much larger line width, implying a much deeper gravitational potential well, as CO line emission traces the kinematics of the potential well where the molecular gas lies. Comparing the ${\mathit{L}}_{\mathrm{CO}}^{\mathrm{\prime }}$ between the one observed directly and the one predicted from the ${\mathit{L}}_{\mathrm{CO}}^{\mathrm{\prime }}$-FWHM relation discovered in Bothwell et al. (2013), we find the observed ${\mathit{L}}_{\mathrm{CO}}^{\mathrm{\prime }}$ is in good agreement with the prediction for all our three SMGs, indicating that the dynamics within the gas emission region is dominated by the baryons. Given that the stellar mass of GN20.2a is the smallest one among these three galaxies, we infer that a large reservoir of molecular gas dominates the dynamics of the region probed by our observations. This is consistent with the heavily dust-obscured UV emission observed in this galaxy. Similar to GN20, the physical properties of GN20.2a would point toward a major merger driven stage, likely approaching the final coalescence with substantial radio emission powered by an AGN, given the high gas surface density and extreme bright radio emission (Daddi et al. 2009b; Hodge et al. 2013). However, the complex optical morphology makes it difficult to understand the detailed physical process of star formation in this galaxy.

Comparing the physical properties of GN20.2b with those of GN20 and GN20.2a, we find that GN20.2b exhibits comparable gas surface density and perhaps even distance with respect to the MS of the SFR-M_⋆ plane with that observed in normal star-forming galaxy at z ~ 2. However, the optically thick UV emission and the high SFE observed in GN20.2b agree better with the case of a starburst. In combination with the analysis in Sect. 3.8, the relatively old stellar age and short gas-consumption timescale of GN20.2b are likely to suggest that this galaxy has already passed its peak of star formation activity, probably observed at the decaying stage of a major merger. In addition, we note that the line width (~220 ± 43 km s^-1) of GN20.2b is much narrower (by a factor of ~2.5–3.5) than those of GN20 and GN20.2a. Given the extreme obscuration of UV/optical emission in GN20.2a and GN20.2b, higher resolution and sensitivity CO imaging would be highly beneficial for the exploration of gas morphology and kinematics in these systems.

Despite the large uncertainties associated with the above analysis, our sample of three SMGs appears to be in different evolutionary stages: GN20 – a major merger undergoing final coalescence; GN20.2a – a major merger approaching final coalescence; and, GN20.2b – at the decaying stage of a merger-driven starburst. We consider that it could be reasonable to find starbursts at different stages in the same massive proto-cluster structure at z = 4.05, as it would be hard to expect absolute synchronization to much better than the 100 Myr timescale of a merger-driven starburst.

5. Implications for the cosmic evolution of dust content in galaxies

The aim of this section is to study the interplay between metals, gas, and dust grains, and the effect of metallicity evolution on observations of galaxies, by investigating the evolution of CO-luminosity-to-dust mass ratio and dust-to-stellar mass ratio across cosmic time for starbursts as well as for MS galaxies, thus placing the observations of the GN20 proto-cluster SMGs in the general context.

5.1. The evolution of L$\begin{array}{c}\mathrm{\prime }\\ \mathrm{CO}\end{array}$/M_dust

We start by looking at the cosmic evolution of the dust mass to CO luminosity ratios, looking for empirical trends and trying to understand possible expectations for the relative evolution of these two quantities versus redshift, the effect of metallicity, and searching for possible different behaviors of MS versus SB galaxies. Figure 7 shows the evolution of ${\mathit{L}}_{\mathrm{CO}}^{\mathrm{\prime }}\mathit{/}{\mathit{M}}_{\mathrm{dust}}$ with redshift for normal star-forming galaxies and starbursts. For MS galaxies, the sample available for our analysis consists of z ~ 0 local group galaxies from Leroy et al. (2011) and individually detected MS galaxies at z ~ 0.5 and z ~ 1.5 from Magdis et al. (2012a). We also show z ~ 1.4−3.1 lensed star-forming galaxies from Saintonge et al. (2013), although for lensed galaxies in general there is the possibility of a “lensing bias” affecting their selection (high surface brightness in the UV rest frame, and/or high intrinsic IR luminosity), so that it is not clear if they should be considered as typical MS galaxies, often times they could be starbursts (see also Tan et al. 2013). To reduce the possible uncertainty resulting from the assumption of excitation correction, for two of the UV-bright lensed sample galaxies, cB58 and Cosmic Eye, we adopt the CO(1–0) luminosity measured by Riechers et al. (2010b) instead of the one extrapolated from CO(3–2) shown in Saintonge et al. (2013). The dust masses of these galaxies are derived based on the DL07 model. For the starburst galaxies, besides our three z = 4.05 SMGs, we add additional objects located at z ~ 0−6.3 from the literature (Downes & Solomon 1998; Coppin et al. 2010; Gilli et al. 2014; Cox et al. 2011; Swinbank et al. 2010; Riechers et al. 2010a, 2011, 2013; Capak et al. 2011; Dwek et al. 2011; Walter et al. 2012). To be consistent with our targets, we also use DL07 models to estimate dust mass for the literature sources, except for ID 141 at z = 4.24, LESS J033229.3-275619 at z = 4.76, and HDF850.1 at z = 5.183, for which we adopt the dust mass estimates in the literature after a correction factor of 2, given that the dust mass estimates based on a single-temperature modified blackbody model are on average a factor of ~2 lower compared to those derived using DL07 model (Magdis et al. 2012a; Magnelli et al. 2012b). For the z = 6.3 SMG (Riechers et al. 2013), our DL07 estimate of M_dust is much higher than the published one, and we adopt a log-averaged value which is within 0.2 dex from both estimates.

Fig. 7 Redshift evolution of ${\mathit{L}}_{\mathrm{CO}}^{\mathrm{\prime }}\mathit{/}{\mathit{M}}_{\mathrm{dust}}$ for MS galaxies and starbursts. At z ~ 0 we use the average value derived from the samples of Leroy et al. (2011) for the spirals (gray circle). Black squares represent individually detected normal star-forming galaxies at z ~ 0.5–1.5 (Daddi et al. 2010; Magdis et al. 2012a), while the black triangles denote lensed galaxies at z ~ 1.4–3.1 (Saintonge et al. 2013). The best fit to the normal galaxies yields a relation of ${\mathit{L}}_{\mathrm{CO}}^{\mathrm{\prime }}\mathit{/}{\mathit{M}}_{\mathrm{dust}}\mathrm{\propto }\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{\mathrm{0.27}\mathrm{\pm }\mathrm{0.20}}$ (solid line). In addition to GN20, GN20.2a, and GN20.2b in our sample at z = 4.05, we also plot starbursts (grey and color stars) from the literature as follows: z ~ 0 ULIRGs (Downes & Solomon 1998); SMM J2135-0102 at z = 2.325 (Swinbank et al. 2010); HSLW-01 at z = 2.957 (Riechers et al. 2011); ID 141 at z = 4.243 (Cox et al. 2011); LESS J033229.4-275619 at z = 4.760 (Coppin et al. 2010; Gilli et al. 2014); HDF 850.1 at z = 5.183 (Walter et al. 2012); AzTEC-3 at z = 5.298 (Capak et al. 2011; Riechers et al. 2010a; Dwek et al. 2011); HFLS3 at z = 6.34 (Riechers et al. 2013). The best fit to the starburst galaxies gives a relation of ${\mathit{L}}_{\mathrm{CO}}^{\mathrm{\prime }}\mathit{/}{\mathit{M}}_{\mathrm{dust}}\mathrm{\propto }\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{\mathrm{-}\mathrm{0.62}\mathrm{\pm }\mathrm{0.23}}$ (dashed line).

The best fit to the MS galaxies yields a nearly flat ratio with redshift of ${\mathit{L}}_{\mathrm{CO}}^{\mathrm{\prime }}\mathit{/}{\mathit{M}}_{\mathrm{dust}}\mathrm{\propto }\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{\mathrm{0.27}\mathrm{\pm }\mathrm{0.20}}$ from z = 0 to 3.1, indicating that ${\mathit{L}}_{\mathrm{CO}}^{\mathrm{\prime }}\mathit{/}{\mathit{M}}_{\mathrm{dust}}$ remains roughly constant for normal galaxies. The best fit to the starburst galaxies gives a relation of ${\mathit{L}}_{\mathrm{CO}}^{\mathrm{\prime }}\mathit{/}{\mathit{M}}_{\mathrm{dust}}\mathrm{\propto }\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{\mathrm{-}\mathrm{0.62}\mathrm{\pm }\mathrm{0.23}}$. The GN20, GN20.2a and GN20.2b galaxies agree well with the general SB trend. Local ULIRGs have a mean ratio of ${\mathit{L}}_{\mathrm{CO}}^{\mathrm{\prime }}\mathit{/}{\mathit{M}}_{\mathrm{dust}}$ higher than normal galaxies by factors of 3–4, but the higher redshift starburst galaxies show consistent ${\mathit{L}}_{\mathrm{CO}}^{\mathrm{\prime }}\mathit{/}{\mathit{M}}_{\mathrm{dust}}$ with normal galaxies. These result do not change substantially, regardless of whether we include the lensed galaxies in the normal galaxy sample or in the starburst sample.

We now attempt a simple interpretation of these trends, starting from the MS galaxies. We recall that both ${\mathit{L}}_{\mathrm{CO}}^{\mathrm{\prime }}$ and M_dust are intimately related to the gas mass, with ${\mathit{L}}_{\mathrm{CO}}^{\mathrm{\prime }}\mathrm{=}{\mathit{M}}_{{\mathrm{H}}_{\mathrm{2}}}\mathit{/}{\mathit{\alpha }}_{\mathrm{CO}}$ and M_dust = M_gas/δ_GDR, respectively. Therefore, the ratio of ${\mathit{L}}_{\mathrm{CO}}^{\mathrm{\prime }}\mathit{/}{\mathit{M}}_{\mathrm{dust}}$ can be expressed as δ_GDR/α_CO, provided that most of the gas contained in the region where the CO and IR luminosity arise is molecular. As discussed in Sect. 3.6.2, H₂ is thought to dominate the gas mass of normal massive M_⋆> 10¹⁰ M_⊙ galaxies at low and high redshift, based on the large molecular gas fraction observed in high-z MS galaxies (Daddi et al. 2010; Tacconi et al. 2010) and the theoretical arguments (e.g., Obreschkow & Rawlings 2009; Bournaud et al. 2011; Lagos et al. 2011). It has been found that the gas-to-dust mass ratio is metallicity dependent with a relation of δ_GDR ∝ Z^-1 for normal galaxies (see Leroy et al. 2011; Sandstrom et al. 2013). Meanwhile, studies of CO-to-H₂ conversion factor have shown that α_CO also scales with metallicity on average as ∝Z^-1 for metallicities not much below 1 / 3 solar (see review by Bolatto et al. 2013; Genzel et al. 2012; Sargent et al. 2014). Based on these derivations, the ratio of ${\mathit{L}}_{\mathrm{CO}}^{\mathrm{\prime }}\mathit{/}{\mathit{M}}_{\mathrm{dust}}$ should not really depend strongly on metallicity (hence neither on stellar mass nor on redshift), for massive galaxies.

One might wonder if this expectation would change if the metallicity of massive galaxies in the high redshift Universe was substantially lower. Recent studies of dwarf galaxies have suggested that both α_CO and δ_GDR could rise more rapidly at metallicities <10 times smaller than solar (e.g., with α_CO ∝ Z^-2.4 in Schruba et al. 2012). Studies of gas-to-dust mass ratios of local galaxies over a large metallicity range reveal a broken power-law relation between δ_GDR and metallicity (Rémy-Ruyer et al. 2014; Fisher et al. 2014), with a steeper slope (∝Z^-3) for dwarf galaxies with metallicities lower than ~8.0. Also in this case, the ratio of ${\mathit{L}}_{\mathrm{CO}}^{\mathrm{\prime }}\mathit{/}{\mathit{M}}_{\mathrm{dust}}$ would remain roughly constant as normal star forming galaxies with higher metallicity, i.e., the ratio of ${\mathit{L}}_{\mathrm{CO}}^{\mathrm{\prime }}\mathit{/}{\mathit{M}}_{\mathrm{dust}}$ changes little with metallicity and thus cosmic time. Of course, the ratio could instead decrease substantially in any regime in which the molecular to total hydrogen fraction was substantially smaller than unity.

Focusing now on the SB galaxies, we already noted that the local ULIRGs show comparatively large ${\mathit{L}}_{\mathrm{CO}}^{\mathrm{\prime }}\mathit{/}{\mathit{M}}_{\mathrm{dust}}$ than that of normal galaxies of similar metallicity (see also Fig. 6 of Magdis et al. 2012a). It is interesting to try to interpret this empirical finding. As the ULIRGs are most likely H₂ dominated, this might be reconnected to either a larger δ_GDR and/or a lower α_CO than normal galaxies, at fixed metallicities (i.e., stellar masses). In Fig. 5 of Magdis et al. (2012a), the local ULIRGs are shown with on average lower α_CO than the one derived based on the α_CO − Z relation defined by normal galaxies, while the gas-to-dust mass ratios follow the local trend. Note that the metallicities derived in that study for ULIRGs are based on optical emission lines, and might be a biased estimate as they might be not very sensible to the most obscured regions.

As for the possible decline of the ${\mathit{L}}_{\mathrm{CO}}^{\mathrm{\prime }}\mathit{/}{\mathit{M}}_{\mathrm{dust}}$ ratio of SBs with redshifts, we suspect that it might be in part spurious. In fact, the local ULIRGs are more extreme starbursts than high-redshift ones, having offsets from the MS of order of 10–30, while for example GN20 and GN20.2a have only six times higher SFR than MS galaxies of similar masses. Matched samples in terms of the distance to the MS would be required to investigate the possible differential redshift evolution in the properties of SB galaxies.

All in all, these results suggest that the ratio of CO luminosity to dust mass remains approximately constant through cosmic time, albeit with larger ratio for starburst galaxies at z = 0. A decline in metallicity could affect both quantities strongly, but in the same way at first approximation. Hence, similarly to the scenarios briefly sketched in Tan et al. (2013) for the evolution of the CO luminosity at high redshift (as for example measured respect to the average CO luminosity at fixed SFR in z< 2 samples), in the case of rapidly declining metallicity in the distant Universe the dust mass of galaxies would be also substantially affected (and therefore their bolometric IR luminosities). In the next section we explore this issue in more detail, again trying to distinguish MS and SB galaxies.

5.2. The observed evolution of M_dust/M_⋆

To guide the understanding of what could regulate the evolution of dust masses to high redshift, and to clarify respect to which physical properties of the galaxies the dust masses should be compared to assess their behavior, it is useful to consider again the definition of δ_GDR, rewriting it in a more convenient form (e.g., Leroy et al. 2011; Magdis et al. 2012a): $\begin{array}{ccc}{\mathit{M}}_{\mathrm{dust}}\mathrm{~}\mathrm{0.5}\mathrm{\times }\mathit{Z}\mathrm{\times }{\mathit{M}}_{\mathrm{gas}}\mathit{.}& & \end{array}$(7)This uses a rough estimate of 50% for the fraction of metals in dust, as on average 1% of the gas is incorporated in dust in the local universe, while solar metallicity of Z_⊙ = 0.02 means that 2% of the gas is in metals. A similar value of the dust-to-metal ratio (~0.5) was found by studies of extinction and metal column densities for a sample of γ-ray burst afterglows and quasar foreground absorbers at z = 0.1−6.3 (Zafar & Watson 2013). Dividing by stellar mass we obtain: $\begin{array}{ccc}{\mathit{M}}_{\mathrm{dust}}\mathit{/}{\mathit{M}}_{\mathit{\star }}\mathrm{~}\mathrm{0.5}\mathrm{\times }\mathit{Z}\mathrm{\times }\mathrm{(}{\mathit{M}}_{\mathrm{gas}}\mathit{/}{\mathit{M}}_{\mathit{\star }}\mathrm{)}& & \end{array}$(8)which shows that the dust to stellar mass ratio depends on the metallicity and the gas fraction in the galaxy, and is expected to be only weakly dependent on stellar mass at fixed redshift for MS galaxies (scaling like ${\mathit{M}}_{\mathit{\star }}^{-0.35}$ based on the dependencies of Z and M_gas/M_⋆ discussed by Magdis et al. 2012a). Instead, at fixed stellar mass, the expected redshift evolution of M_gas/M_⋆ can be simply computed starting from the observed cosmic evolution of the sSFR in MS galaxies (see e.g., Fig. 19 in Sargent et al. 2014) and converting the SFR into gas masses using the Schmidt-Kennicutt law (S-K; $\mathit{SFR}\mathrm{\propto }{\mathit{M}}_{\mathrm{gas}}^{\mathrm{1.2}}$; see e.g., Daddi et al. 2010b; Genzel et al. 2010; Sargent et al. 2014). The remainder of the redshift evolution of the M_dust/M_⋆ ratio is thus contained in the evolution of metallicity (or, at fixed stellar mass, of the mass-metallicity relation).

Fig. 8 Evolution of Mdust/M⋆ as a function of redshift. The gray circle represents the averaged value of local normal galaxies from the sample of da Cunha et al. (2010a), while black solid circles represent the stacking results at z ~1 and z ~2. The light blue and blue solid circles represent “LBG-1” and “Himiko”, which are star-forming galaxies at z = 5.300 and z = 6.595, respectively. The black solid line shows the expected evolution of Mdust/M⋆ with redshift based on the relation observed in the local universe, Mdust ~ 0.5 × Mgas × Z. Predictions for the redshift-evolution of Mdust/M⋆ of MS galaxies are shown for the case of an universally valid FMR/broken FMR/evolving M⋆ − Z relation at z ≥ 3 (red/orange/green lines; the z> 2.5 extrapolations shown in colored dashed lines are observationally-unconstrained). The shaded region shows the range of Mdust/M⋆ ratio beyond z = 2.5 predicted by our models. All measurements for normal galaxies and predictions have been normalized to a common mass scale of M⋆ = 5 × 1010 M⊙. The star symbols indicate starbursts, same as in Fig. 7. The black dot-dashed line is the best fit to the starbursts with a relation of Mdust/M⋆ ∝ (1 + z)0.51 ± 0.35.

In Fig. 8 we show measurements of M_dust/M_⋆ for MS and SB galaxies at different redshifts, normalized to a stellar mass of 5 × 10¹⁰ M_⊙ using the scaling discussed above. We show average values for z ~ 0 spiral galaxies (da Cunha et al. 2010a) and stacked samples of z ~ 1 and z ~ 2 normal galaxies (Magdis et al. 2012a). The use of average samples is particularly useful as it should be representative of the typical behavior of MS galaxies, independently on object-to-object fluctuations. Recent studies of nearby galaxies of the Herschel Reference Survey derive a M_dust/M_⋆ of 3.5 × 10^-3 at z ~ 0 (Cortese et al. 2012; Ciesla et al. 2014), in good agreement with our estimate for SDSS galaxies (da Cunha et al. 2010a). For typical MS galaxies, the dust-to-stellar mass ratio rises by a factor of few from z = 0 to 2 (consistent with the findings of Scoville et al. 2014; Santini et al. 2014). Similar trends have also been revealed for large samples of Herschel-ATLAS galaxies at z< 0.5 (Dunne et al. 2011). Given that the lensed star-forming galaxies are intrinsically faint with low stellar mass (M_⋆ ≲ 10¹⁰ M_⊙) and probably not representative of massive galaxies at those epochs, we did not compare the dust-to-stellar ratio of these galaxies with MS galaxies, although some lensed galaxies at z ~ 1.5–3 are shown to have similar M_dust/M_⋆ ratio to local star-forming galaxies (Sklias et al. 2014). For the SB galaxies we plot the same sample as in Fig. 7 (Coppin et al. 2010; Gilli et al. 2014; Cox et al. 2011; Swinbank et al. 2010; Riechers et al. 2010a, 2011, 2013; Capak et al. 2011; Dwek et al. 2011; Walter et al. 2012). The dust to stellar mass ratio of SB galaxies is found to be fairly flat with redshift, with the best fitting trend being M_dust/M_⋆ ∝ (1 + z)^{0.51 ± 0.35} (see Fig. 8), indicative of still substantial metal enrichment at higher redshifts in these systems up to at least z ~ 6. The dispersion of the residuals from the fit is larger than what expected from the measurement uncertainties, indicating intrinsic dispersion of at least a factor of two in the properties of the SB population. Again, the GN20 proto-cluster galaxies behave similarly to the rest of the SBs.

No dust mass measurements exist for indisputably normal, MS galaxies at z> 3 (an attempt to study the evolution of M_dust/M_⋆ at 3 <z< 4 by stacking of normal MS galaxies will be presented in an upcoming paper by Béthermin et al., in prep.). However, the recent ALMA observations of Himiko, a star-forming Ly-α selected galaxy at z = 6.595 with SFR ~ 100 M_⊙ yr^-1 and M_⋆ ~ 1.5 × 10¹⁰ M_⊙, reveal a significant deficit of dust content and [CII] 158 μm emission (Ouchi et al. 2013), which is more than 30 times weaker than the one predicted by local correlations. With the 1.2 mm flux limit measured in Ouchi et al. (2013), we estimate a dust mass assuming a radiation field U to an appropriate value for a z = 6.6 normal galaxy by extrapolating the relation of ⟨ U ⟩ ∝ (1 + z)^1.15 derived in Magdis et al. (2012a). This upper limit in dust mass, coupled with the substantial stellar mass of a few 10¹⁰ M_⊙ in the system (Ouchi et al. 2013) implies a M_dust/M_⋆ ratio that is more than an order of magnitude lower than the trend defined by the SB galaxies at z> 5. Analogously, we estimate a dust mass of a z = 5.300 star-forming LBG called “LBG-1” based on the 1.0 mm flux limit observed with the ALMA (Riechers et al. 2014a) and measure the stellar mass by fitting stellar population models to the rest-frame ultraviolet–optical photometry (Capak et al. 2011). We find an upper limit of M_dust/M_⋆ of 6.5 × 10^-3 for this galaxy, significantly lower than the SB galaxies at the same epoch. These results are qualitatively similar to the recent observations toward two LBGs at z> 3 reporting a deficiency in CO emission, supporting the scenario of a rapid decline in metallicity (Tan et al. 2013).

5.3. A crucial ingredient: metallicity evolution in the distant Universe

Fig. 9 Upper panel: redshift evolution of the metallicity (normalized to solar) of a star-forming galaxy with M⋆ = 5 × 1010 M⊙ for the case of an universally valid FMR/ broken FMR/ evolving M⋆ − Z relation at z> 2.5 (dashed/solid/grey lines; the z> 3.3 extrapolations are observationally unconstrained). Lower panel: mass-metallicity relation between z = 0 and 5 for the three cases represented with same line styles as in the upper panel. Lines are color–coded according to the redshift.

We discuss in the following all these empirical findings. We start first with an attempt to interpret the behavior of MS galaxies, and to predict their properties at z> 3 where observations are very scarce. Following the derivation of Eq. (8) and the discussion above, together with the sSFR(z) evolution as inferred e.g., in Sargent et al. (2014), the crucially needed ingredient is the evolution of metallicity for a given mass range (e.g., within 10 < log M_⋆/M_⊙< 11). We thus briefly review here the current understanding for the evolution in metallicity of massive galaxies through cosmic time. There are clearly two regimes that can be emphasized: above or below z ~ 2.5–3. From z = 0 to 2.5 there is only a modest, factor of two or less decrease in the metallicity of massive galaxies (e.g., Erb et al. 2006; Tremonti et al. 2004; Zahid et al. 2013, many others), see Fig. 9-upper. Recent studies on metal content of galaxies have also shown that galaxies up to z ~2.5 follow the FMR defined locally with smooth variations (Mannucci et al. 2010). The factor of a few increase of M_dust/M_⋆ for normal galaxies up to z = 2 (Fig. 8) is thus understood following Eq. (8) in term of the rapid increase of gas fraction (e.g., Daddi et al. 2010; Magdis et al. 2012a; Tacconi et al. 2010, 2013), which strongly overpowers for the metallicity decrement.

However, while the exact evolution of average metallicity of massive galaxies at z> 3 is still in debate and therefore highly uncertain, there are observational claims suggesting that things might change quite drastically at these high redshifts. Evidence for rapid metallicity evolution at z> 3 has been proposed by studies of z ~ 3.5 LBGs (e.g., Maiolino et al. 2008; Mannucci et al. 2010; Sommariva et al. 2012; Troncoso et al. 2014), suggesting that the metallicity of massive M_⋆> 10¹⁰ M_⊙ galaxies might be decreasing fast, and the FMR no more valid over those redshift ranges (Fig. 9-upper). It is not clear if these measurements might be unrepresentative for typical z> 3 massive galaxies, as they might be biased down in reddening (hence in metallicity) due to the UV-selection of the spectroscopic samples (see Zahid et al. 2014). Moreover, these observational claims are in contrast with what is predicted by some theoretical (e.g., Finlator & Davé 2008) or heuristic (e.g., Lilly et al. 2013) models. Given the substantial uncertainties, for now we regard the range of scenarios as all plausible, and we explore in particular the implications within a reasonable range going from the case of stable metallicity of massive galaxies at z> 3 (equivalent to non evolution in the FMR) up to the case in which the claim for rapid metallicity decrease (e.g., Troncoso et al. 2014) is correct and representative for the whole population of z ~ 3.5 galaxies. The latter scenario can be computed either assuming a broken FMR framework, or parametrized in terms of the observed evolution of the mass-metallicity relation. For the case of a broken FMR, we interpolate linearly with redshift the metallicity drop of 0.6 dex reported in Mannucci et al. (2010) between z = 2.5 and z = 3.3 and continue this trend to extrapolate beyond z = 3.3. Similarly, we apply linear interpolation/extrapolation between z = 2.2 and z = 3 based on the M_⋆ − Z relation of Sommariva et al. (2012), Troncoso et al. (2014), and Zahid et al. (2013) for the case of evolving M_⋆ − Z relation. In Fig. 9-lower, we further plot the model predictions of the M_⋆-Z relation between z = 0 and z = 5 (color coded) for the above three different cases. For the scenarios involving the FMR, the SFR which was used to compute the metallicity at a given stellar mass is the SFR of the average MS galaxy. These plots illustrate how the mass-metallicity relation evolves with redshift under various hypothesis. For the case of an universally valid FMR, we find that the tracks at z ≥ 3 are almost identical (Fig. 9-lower), with no indication of evolution of the M_⋆ − Z relation above z ~ 3, which is due to the fact that a plateau or slow rise is observed in the evolution of sSFR at z> 3. For the case of accelerated metallicity evolution at z> 2.5, the gas metallicity decreases rapidly with redshift.

5.4. Interpreting the evolution of M_dust/M_⋆

Using again Eq. (8), we convert the different scenarios for the evolution of metallicity into related cosmic evolution of the M_dust/M_⋆ ratio for normal star-forming galaxies with M_⋆ = 5 × 10¹⁰ M_⊙ (Fig. 8). The expected evolution trend of the M_dust/M_⋆ ratio derived directly from Eq. (8) is shown in the black solid line. The metallicity is estimated based on the direct measurements of the mass-metallicity relation from optical surveys (Tremonti et al. 2004; Savaglio et al. 2005; Erb et al. 2006; Maiolino et al. 2008) at z< 3.3 and the extrapolation beyond higher redshift, while the ratio of M_gas/M_⋆ is computed from the observed cosmic evolution of the sSFR and of the integrated S-K relation (Sargent et al. 2014). The colored lines show the evolution trends predicted based on the above three scenarios for the evolution of metallicity at fixed stellar mass (M_⋆ = 5 × 10¹⁰ M_⊙; see Fig. 9), with dust mass converted from the gas mass computed similarly as above and assuming that the local δ_GDR-metallicity relation (Leroy et al. 2011; Magdis et al. 2012a) is valid at high redshift. We find the predictions based on these three realizations are basically undistinguishable up to z ~ 2.5, with an increase evolution of M_dust/M_⋆ with redshift. For galaxies beyond z ~ 2.5, however, the universally valid FMR would predict a mild increase of M_dust/M_⋆, while a rapid decline is expected for the case of accelerated evolution. The shaded region in Fig. 8 shows the predicted range of M_dust/M_⋆ ratio beyond z = 2.5, emphasizing how uncertain the situation is.

However, from Fig. 8 we find that the M_dust/M_⋆ ratios of LBG-1 and Himiko are significantly lower than the prediction from the universally valid FMR, consistent with the prediction assuming a rapid decline in metallicity at z> 3. Unless there are something special with LBG-1 and Himiko, these two galaxies observations appear to support scenarios with declining metallicities for massive galaxies, at least in the early Universe.

The scenarios with declining high-redshift metallicity would produce a broad peak of dust to stellar mass content in galaxies over at z> 1 and up to 3 (or maybe higher, should a decline happen at earlier redshifts) that would allow one to term the z ~ 2 range as the epoch of dusty galaxies, which would occur roughly at the same time as the peak in the SFRD. Recent studies of the IR to UV luminosity density ratio versus redshift depict a similar situation, with a broad z ~ 2 peak characterizing the ratio (Cucciati et al. 2012; Burgarella et al. 2013). There seems to be a very similar behavior in the dust-to-stellar mass ratio and the average dust attenuation in galaxies, at least to z ~ 2.5, which could somewhat support again the scenarios with declining metallicities at high redshifts.

We finally spend some words of caution noticing that, again following Eq. (8), a different evolution of gas-to-stellar mass ratio than what assumed here, could compensate in principle evolutionary trends in metallicity and influence the dust-to-stellar mass ratio evolution. Current studies show no direct evidence of further increase of gas fraction at z> 3 (e.g., Magdis et al. 2012b; Carilli & Walter 2013; Sargent et al. 2014; Tan et al. 2013), and the sSFRs are revealed not to be rapidly rising as well (e.g., Bouwens et al. 2012). However it might well be that the gas-to-stellar mass ratio is actually rapidly rising at z> 3 even if not seen by the sSFR (e.g., Dekel & Krumholz 2013). The latter case would imply a decrease of SFE and deviation from the S-K law.

We finish this discussion section with some consideration on the behavior of the SB galaxies. The local ULIRGs taken from the sample of da Cunha et al. (2010b) show substantially higher mean dust-to-stellar mass ratio than that of z ~ 0 disk galaxies, over a comparable stellar mass range, suggesting that they have more dust or more metals (or probably both as dust grains form from the available metals in the ISM) than normal star forming galaxies at fixed stellar mass. The dust-to-stellar mass ratios derived for the SBs in our sample ranges between 0.01 and 0.09 with mean value of 0.04, in agreement with the M_dust/M_⋆ ratios of z> 1 SMGs measured by Rowlands et al. (2014). Starburst galaxies appear to saturate their M_dust/M_⋆ ratios at a value of few × 10^-2, which is unlikely to be a coincidence. This corresponds to the solar to few times solar metallicities of massive elliptical galaxies, into which these objects might be evolving once the effects of the undergoing merger event has faded. If anything, the dust to stellar mass ratios for SBs might actually increase toward high redshifts, meaning that even in the early Universe SBs were substantially metal rich, pointing toward the need of rapid metal enrichment in the starburst phase. Studies of the origin of dust in galaxies using chemical evolution models found that, excepting asymptotic giant branch stars, supernovae and grain growth could be significant sources of dust to account for the high dust masses observed in star-bursting systems at high redshift (Rowlands et al. 2014; Michałowski et al. 2010a,b). Direct evidence for the supernova to be a primary source of dust in the early Universe has been found by ALMA observations of Supernova 1987A, of which a remarkably large dust mass is revealed to be concentrated at the center of the remnant (Indebetouw et al. 2014).

All of these suggest that the measurement of M_dust/M_⋆ could be a very powerful mean for distinguishing starbursts from normal galaxies at z> 4, during the phase of early metal enrichment in starburst galaxies.

6. Summary

We have presented results from deep IRAM PdBI CO(4–3) and 1.2–2.2–3.3 mm continuum observations of the GN20 proto-cluster at z = 4.05. The improved CO spectral profile of GN20, GN20.2a, and GN20.2b allow us to measure the line width more accurately and further constrain the dynamical mass. Combining with the ancillary multiwavelength photometry in the rest-frame UV, optical, and IR, we determine the stellar masses and dust properties of these three submm galaxies (SMGs). With the measured stellar masses, dynamical masses, CO luminosities, dust masses, and indirect metallicity estimates, we inferred the value of conversion factor α_CO in each of the three SMGs, using the dynamical modeling and the gas-to-dust ratio methods. Combining with literature data of normal galaxies and starbursts from local to high redshift, we discuss the effect of metallicity evolution on observations of dust and gas emission of galaxies across cosmic time. The main results and implications are summarized as follows:

• 1.

  All three SMGs are now detected in CO(4–3) with high S/Ns.The FWHM of the spectra derived from double Gaussian fits are583 ± 36 km s^-1,760 ± 180 km s^-1, and220 ± 43 km s^-1 for GN20, GN20.2a, andGN20.2b, respectively. With ${\mathit{L}}_{\mathrm{CO}\mathrm{(}\mathrm{4}\mathrm{-}\mathrm{3}\mathrm{)}}^{\mathrm{\prime }}$ derived from our study and ${\mathit{L}}_{\mathrm{CO}\mathrm{(}\mathrm{2}\mathrm{-}\mathrm{1}\mathrm{)}}^{\mathrm{\prime }}$measured by Carilli et al. (2011),we find CO(4–3)/CO(2–1) line ratio of ~0.4 for these three SMGs,which is consistent with the mean ratio (~0.48 ± 0.10) measured for SMGs at z ~ 2–4(Bothwell et al. 2013).

• 2.

  We report 3.3 mm continuum detections in GN20.2a and GN20.2b for the first time, and use continuum measurements at 1.2 and 2.2 mm (Carilli et al. 2010; Riechers et al., in prep.). The dust masses measured from the far-IR SED fitting are 2.1−5.2 × 10⁹ M_⊙. The IR-to-radio luminosity ratios for GN20 (q = 2.41 ± 0.05) and GN20.2b (q = 2.60 ± 0.13) are found to be comparable to the local value (~2.6), while GN20.2a shows a relatively low value (q = 1.72 ± 0.04), suggesting that this galaxy is radio-excess with large amounts of radio emission likely powered by an AGN (see also Daddi et al. 2009).

• 3.

  We find that the value of α_CO inferred from gas-to-dust ratio method is consistent with the one derived based on the dynamical mass for each galaxy. The α_CO derived for these three SMGs (~1.3–2.8 M_⊙ (K km s^-1 pc²)^-1) are found to be consistent with the typical value determined for local ULIRGs, but might be well below the value appropriate for normal galaxies at similar redshifts. The high gas fraction (~40%–80%) of these three SMGs at z=4.05 are found to be comparable to SMGs at z = 2–4 and high redshift normal galaxies.

• 4.

  Our study clearly distinguishes GN20, GN20.2a (and likely GN20.2b) as starbursts from normal star-forming galaxies by comparing observed physical properties between these galaxies. For GN20 and GN20.2a, the large sSFR-excess (sSFR/sSFR_MS ~ 6) compared to the normal galaxies at similar epochs place these SMGs as outliers above the main sequence. The extremely large value of SFR_IR/SFR_UV for GN20 and GN20.2a are consistent with the large offset between CO positions and optical counterparts (see Fig. 4), suggesting that the UV/optical emission is heavily obscured beyond optically thick dust. Although GN20.2b situates within the MS scatter with sSFR/sSFR_MS ~ 2.4, both the large value of SFR_IR/SFR_UV (~6) and short gas-consumption timescales (~116 Myr) suggest that this galaxy could also be a starburst.

• 5.

  We find that these three SMGs are likely to experience different evolutionary stages of starburst activity. Compared with GN20 and GN20.2a, GN20.2b displays relatively smaller sSFR-excess, older stellar age, and lower CO excitation. All these properties suggest that GN20.2b is probably observed at a decaying stage of a major merger, while GN20 and GN20.2a are likely to undergo and approach the final coalescence with intense starburst, respectively.

• 6.

  We compile a variety of archival data of normal galaxies and starbursts in order to investigate the effect of metallicity evolution on observations of galaxies across the cosmic time. For normal galaxies, the ratio of ${\mathit{L}}_{\mathrm{CO}}^{\mathrm{\prime }}\mathit{/}{\mathit{M}}_{\mathrm{dust}}$ is found to be almost constant from z = 0 to 3.1. And the same appear to be the case for high-z starbursts, implying that both the CO and dust emission could be affected in the same way by the metals in the ISM.

• 7.

  We calculate simple models for the dust emission of normal galaxies based on their properties and find a rapid decrease of dust emission for z> 3 normal galaxies at a given stellar mass if metallicities indeed decrease rapidly for these galaxies. The model predictions are well-matched to the measurements at z< 2.5, with a trend of increasing M_dust/M_⋆ with redshift. While no dust detection is available for indisputably normal galaxies beyond z = 2.5, the model based on the assumption of fast metallicity evolution predicts a sharp decline of dust emission for normal galaxies at z> 2.5. In contrast, the starburst galaxies show an increase of M_dust/M_⋆ at high redshift, providing evidence for rapid early metal enrichment in these systems. The different behaviors of normal galaxies lead to a significantly lower M_dust/M_⋆ compared to starbursts at z> 4, implying that the comparison of M_dust/M_⋆ ratio could also be used as a powerful tool for distinguishing starbursts from normal galaxies at z> 4.

Acknowledgments

This work was based on observations carried out with the IRAM PdBI, supported by INSU/CNRS (France), MPG (Germany), and IGN (Spain). The authors thank Raphaël Gobat for assistance with the SED fitting of LBG-1, and we thank the anonymous referee for useful comments and suggestions which greatly improved this manuscript. Q.T., E.D., and M.T.S. acknowledge the support of the ERC- StG UPGAL 240039 and ANR-08-JCJC-0008 grants. Q.T. was supported by CAS and CNRS, and partially supported by the NSFC grant #11173059, #11390373, and the CAS grant #XDB09000000.

References

All Tables

All Figures

Fig. 1 CO(4–3) spectra binned in steps of 25 km s-1 for GN20 (left), 50 km s-1 for GN20.2a (middle) and GN20.2b (right). The yellow color indicates the velocity ranges where line emission is detected. These regions have been used to derive the integrated fluxes. The red lines show best-fitting double Gaussian profiles to the spectra, and the blue dashed lines show the fits with a single Gaussian. Zero velocities correspond to the redshifts listed in Table 1.

In the text

Fig. 2 PdBI maps (from left to right) of CO(4–3) line, 3.3 mm, and 2.2 mm continuum emission for GN20 (top), GN20.2a (middle), and GN20.2b (bottom). The CO(4–3) images are averaged over the observed velocity range of CO emission (see Fig. 1) and were cleaned. The 2.2 mm and 3.3 mm continuum maps were created by averaging over the line-free emission channels and were not cleaned given the relatively low S/N. The crosses show the positions where we extract CO spectra. For GN20, the contours level of CO(4–3) are shown increase by a factor of 1.5 starting at ±2σ, with positive (negative) contours shown as solid (dashed) lines. For all the continuum maps and CO(4–3) maps of GN20.2a and GN20.2b, contours start at ±2σ in steps of 1σ. Beam sizes are displayed in the lower right corner of each integrated map.

In the text

Fig. 3 Observed mid-IR to millimeter SED of GN20, GN20.2a, and GN20.2b, overlaid with the best-fit Draine & Li (2007) models. The black line is the DL07 model consisting of one stellar component and two dust components: diffuse ISM and photodissociation region. Green dot-dashed and cyan dot-dashed lines show stellar component and “PDR” component, respectively. The best-fit parameters from DL07 models are listed within each panel. The arrows represent 5σ upper limits.

In the text

Fig. 4 PdBI CO(6–5)+continuum contours (Carilli et al. 2010) at resolution of 0.90′′ × 0.76′′ for GN20 (inset), GN20.2a and GN20.2b overlaid on the ACS 850z-band (top left), WFC3 F160W-band (top right), WIRCam Ks-band (bottom left), and IRAC 3.6 μm (bottom right) images. GN20 lies about 24′′ to the northeast of GN20.2a. Black and red crosses show the positions of the 1.4 GHz counterparts (Owen et al., in prep.) and optical counterparts of GN20, GN20.2a and GN20.2b, respectively. Contour levels are shown in steps of 1σ for GN20.2a and GN20.2b and 3σ for GN20 starting at ±3σ. The beam sizes of CO(6–5) are displayed in the left corner of each map.

In the text

Fig. 5 Observed optical to mid-IR SED of GN20 (left), GN20.2a (middle), and GN20.2b (right). Red circles with error bars represent the multiwavelength photometry. The solid lines show the best-fitting templates from the libraries of Bruzual & Charlot (2003) at the derived CO redshift for each galaxy, assuming a constant SFH models with solar metallicity and using an extinction correction from Calzetti et al. (2000). We adopt the templates with a constant SFH to derive stellar mass by fitting the photometry of U- to IRAC 4.5 μm band.

In the text

	Fig. 6 Global UV to far-IR fits to our galaxies based on MAGPHYS. The black line is the best-fit SED and the blue line is the corresponding unattenuated stellar emission.
In the text

Fig. 7 Redshift evolution of ${\mathit{L}}_{\mathrm{CO}}^{\mathrm{\prime }}\mathit{/}{\mathit{M}}_{\mathrm{dust}}$ for MS galaxies and starbursts. At z ~ 0 we use the average value derived from the samples of Leroy et al. (2011) for the spirals (gray circle). Black squares represent individually detected normal star-forming galaxies at z ~ 0.5–1.5 (Daddi et al. 2010; Magdis et al. 2012a), while the black triangles denote lensed galaxies at z ~ 1.4–3.1 (Saintonge et al. 2013). The best fit to the normal galaxies yields a relation of ${\mathit{L}}_{\mathrm{CO}}^{\mathrm{\prime }}\mathit{/}{\mathit{M}}_{\mathrm{dust}}\mathrm{\propto }\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{\mathrm{0.27}\mathrm{\pm }\mathrm{0.20}}$ (solid line). In addition to GN20, GN20.2a, and GN20.2b in our sample at z = 4.05, we also plot starbursts (grey and color stars) from the literature as follows: z ~ 0 ULIRGs (Downes & Solomon 1998); SMM J2135-0102 at z = 2.325 (Swinbank et al. 2010); HSLW-01 at z = 2.957 (Riechers et al. 2011); ID 141 at z = 4.243 (Cox et al. 2011); LESS J033229.4-275619 at z = 4.760 (Coppin et al. 2010; Gilli et al. 2014); HDF 850.1 at z = 5.183 (Walter et al. 2012); AzTEC-3 at z = 5.298 (Capak et al. 2011; Riechers et al. 2010a; Dwek et al. 2011); HFLS3 at z = 6.34 (Riechers et al. 2013). The best fit to the starburst galaxies gives a relation of ${\mathit{L}}_{\mathrm{CO}}^{\mathrm{\prime }}\mathit{/}{\mathit{M}}_{\mathrm{dust}}\mathrm{\propto }\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{\mathrm{-}\mathrm{0.62}\mathrm{\pm }\mathrm{0.23}}$ (dashed line).

In the text

Fig. 8 Evolution of Mdust/M⋆ as a function of redshift. The gray circle represents the averaged value of local normal galaxies from the sample of da Cunha et al. (2010a), while black solid circles represent the stacking results at z ~1 and z ~2. The light blue and blue solid circles represent “LBG-1” and “Himiko”, which are star-forming galaxies at z = 5.300 and z = 6.595, respectively. The black solid line shows the expected evolution of Mdust/M⋆ with redshift based on the relation observed in the local universe, Mdust ~ 0.5 × Mgas × Z. Predictions for the redshift-evolution of Mdust/M⋆ of MS galaxies are shown for the case of an universally valid FMR/broken FMR/evolving M⋆ − Z relation at z ≥ 3 (red/orange/green lines; the z> 2.5 extrapolations shown in colored dashed lines are observationally-unconstrained). The shaded region shows the range of Mdust/M⋆ ratio beyond z = 2.5 predicted by our models. All measurements for normal galaxies and predictions have been normalized to a common mass scale of M⋆ = 5 × 1010 M⊙. The star symbols indicate starbursts, same as in Fig. 7. The black dot-dashed line is the best fit to the starbursts with a relation of Mdust/M⋆ ∝ (1 + z)0.51 ± 0.35.

In the text

Fig. 9 Upper panel: redshift evolution of the metallicity (normalized to solar) of a star-forming galaxy with M⋆ = 5 × 1010 M⊙ for the case of an universally valid FMR/ broken FMR/ evolving M⋆ − Z relation at z> 2.5 (dashed/solid/grey lines; the z> 3.3 extrapolations are observationally unconstrained). Lower panel: mass-metallicity relation between z = 0 and 5 for the three cases represented with same line styles as in the upper panel. Lines are color–coded according to the redshift.

In the text