© S. Jin et al. 2022

1. Introduction

The interstellar medium (ISM) is key to diagnosing the evolutionary state of galaxies in the early Universe. In particular, the available reservoir of molecular hydrogen (M_H2) dictates the ongoing star- and subsequent dust formation. However, measuring the molecular gas mass M_gas of galaxies, both in the local and the high-z Universe is not a trivial task, as H₂ lacks dipole moment and thus is not detectable under the most typical, most representative ISM conditions. As gas and dust are well mixed (e.g., Bohlin et al. 1978; Boulanger et al. 1996), an efficient method to infer M_gas indirectly is to use far-infrared(FIR)/submillimeter(submm) dust continuum observations to derive the (observationally cheaper) dust mass (M_dust) of a galaxy that is subsequently converted to M_gas through the metallicity-dependent dust-to-gas-mass ratio technique (e.g., Leroy et al. 2011; Eales et al. 2012; Magdis et al. 2012, 2017; Santini et al. 2014; Genzel et al. 2015; Berta et al. 2016; Schinnerer et al. 2016; Tacconi et al. 2018, 2020; Liu et al. 2019b; Kokorev et al. 2021). Similarly, single band measurements of the dust emission flux on the Rayleigh-Jeans side of the SED have also been calibrated against CO observations and under a series of assumptions (e.g., dust temperature, α_CO conversion factor) can be used to infer gas mass estimates (e.g., Groves et al. 2015; Scoville et al. 2016, 2017a; Schinnerer et al. 2016; Liu et al. 2019b).

However, the derivation of M_gas from M_dust or from a single band Rayleigh-Jeans (RJ) dust continuum flux density are usually based on the assumption that the observed FIR/submm emission is optically thin. If instead the dust opacity of a galaxy extends to λ > 100 μm, its M_dust based on optically thin models can be severely overestimated (a factor of 2–20 e.g., Galliano et al. 2011; Jin et al. 2019; Cortzen et al. 2020; Fanciullo et al. 2020). Subsequently, the inferred M_gas would also be problematic. While there is mounting evidence for optically thick dust emission in both local and high-z dusty star-forming galaxies (e.g., Blain et al. 2003; Draine & Li 2007; Conley et al. 2011; Casey 2012; Spilker et al. 2016; Hodge et al. 2016; Scoville et al. 2017b; Simpson et al. 2017; Riechers et al. 2014, 2017, 2020; Cortzen et al. 2020), the large degeneracy between the solutions of the thick and thin dust model SEDs prevents a definitive conclusion.

Besides the dust SEDs, submm lines provide more powerful diagnostics for ISM properties, revealing the presence of optically thick dust emission at FIR wavelengths in the starbursting cores of submillimeter galaxies (SMGs) in the local Universe. For example, Papadopoulos et al. (2010) developed a method to identify the optically thick dust using ratios between optically thick and optically thin lines in local galaxies (e.g., HCN(1-0)/CO(1-0) and CO(6-5)/CO(3-2)) and proposed that dust optical depths can significantly suppress high-J CO and C⁺ emission. Based on James Clerk Maxwell Telescope (JCMT) observations, they found that the very faint CO(6-5) lines in Arp220 and other ultraluminous infrared galaxies (ULIRGs) can be attributed to high dust optical depths at FIR wavelengths immersing those lines in a strong dust continuum. Furthermore, using high-resolution ALMA CO(1-0) mapping, Scoville et al. (2017b) confirmed that the dust in the two nuclei of Arp 220 is optically thick well into the FIR, and in one of them even at 2.6 mm with a high dust brightness temperature of 147 K.

Such high optical depths are also expected in the ISM of high-z star-forming galaxies that were more compact and rich in gas and dust. For example, Jin et al. (2019) discovered a sample of “apparently” extremely cold dusty galaxies at z = 4–6 using ALMA line scans. Their dust spectral energy distributions (SEDs) peak at longer rest-frame wavelengths than most literature sources at z > 4, and their dust temperatures (24–42 K) as inferred with optically thin models are between 1.5 and 2.0 times colder than that of main sequence (MS) galaxies at these redshifts. Given the massive dust content and compact morphology of these galaxies, the “apparent” cold dust temperatures are best understood as evidence for large optical depths that effectively attenuate the dust continuum emission in the FIR and shift the peak of the observed SED to longer wavelengths. In the same direction, Cortzen et al. (2020) used the two neutral carbon transitions in GN20, a star-bursting galaxy at z = 4.05, and found that the excitation temperature as inferred by the [CI](2-1)/(1-0) ratio is significantly higher than the abnormally cold dust temperature derived from SEDs under the assumption of an optically thin FIR dust emission. On the other hand, the excitation temperature is fully consistent with the dust temperature of a general opacity model with λ_eff ≈ 170 μm. These findings suggest that high FIR optical depths could be a widespread feature of starbursts across comic time. Follow-up observations of seemingly “cold”, high-z star-bursting galaxies with compact-dust-emitting sizes are necessary to push this investigation forward.

In this paper, we present an in-depth investigation of the ISM properties of six such galaxies based on new NOEMA and ALMA observations. These six galaxies are used to extend the original sample of Jin et al. (2019) to a total of ten, which is used in the remainder of this paper. In Sect. 2 we describe the sample selection and data processing. In Sect. 3 we discuss redshift confirmation, SED modeling, and size measurements. We present evidence of optically thick dust in Sect. 4. We discuss the implication of these results in Sect. 5 and provide our summary and conclusions in Sect. 6. We adopt cosmology H₀ = 73, Ω_M = 0.27, Λ₀ = 0.73; and a Chabrier initial mass function (IMF; Chabrier 2003).

2. Sample, observations, and data reduction

Six new dusty star-forming galaxies were selected for follow-up observations in this study, after being culled for potentially lying at very high redshifts (z > 6) based on their SEDs. We show their multi-wavelength images in Fig. 1, including deep ALMA and NOEMA continuum data. The NOEMA and ALMA observations are summarized in Table 1. Five are in the COSMOS field and one is in the GOODS-North field.

Fig. 1. Multi-wavelength images with 15″ × 15″ sizes. The instrument used and wavelength are shown in green text in each cutout. ID5340, 6309, and 7549 have no HST F160W data available and thus are shown as blank panels.

Among the five COMSOS sources, ID3117, 5340, 6309, and 7549 were selected from the COSMOS super-deblended catalog (Jin et al. 2018). They were originally detected in the VLA 3 GHz catalog (Smolčić et al. 2017) and are optically dark, that is, with no detection in HST I or H band, and are therefore not included in the COSMOS2015 catalog (Laigle et al. 2016). In the COSMOS super-deblended catalog (Jin et al. 2018), these four galaxies are detected with significant FIR emission (S/N_FIR + mm> 5) in Herschel, SCUBA2, AzTEC, and/or MAMBO (Geach et al. 2017; Cowie et al. 2017; Aretxaga et al. 2011; Bertoldi et al. 2007). Fitting of their dust SEDs returns very high photometric redshifts of z_phot, FIR > 6 using the z = 6.3 HFLS3 template (Riechers et al. 2013). We carried out follow-up observations of the four galaxies with NOEMA 3mm line scans (ID: S18DI, PI: S. Jin) in the winter 2018 semester.

ID9316 is a bright submm galaxy close to a galaxy group RO0959 at z = 3.1 (Daddi et al. 2022), and was proposed to be at z > 7 by Oteo and collaborators (see ALMA proposal ID: 2018.1.00874.S, PI: I. Oteo) and was observed with ALMA standard line scans in band 3. Meanwhile, this galaxy has also snapshot imaging in ALMA band 6 (ID: 2013.1.00118.S, PI: M. Aravena) and band 7 (ID: 2016.1.00478.S, PI: O. Miettinen).

ID12646 was selected from the FIR/(sub)mm super-deblended catalog in the GOODS-North field (Liu et al. 2018) as the source with the highest photometric redshift z_phot, FIR ∼ 6. During 2017 to 2019, Liu et al. carried out follow-up analyses of this source with NOEMA 3mm (IDs: W17EK, W19DQ) and 1mm observations (ID: W18EY).

2.1. NOEMA

The four galaxies, ID3117, 5340, 6309, and 7549, were observed with NOEMA 3mm line scans in the winter of 2018 (ID: S18DI, PI: S. Jin), continuously covering the frequency range 74.7–106 GHz. The observations were conducted in track-sharing mode with array configurations C and D. Each source was observed with two to three tracks in each tuning, reaching rms sensitivities of 0.13–0.16 mJy per 500 km s⁻¹ with synthesized beams as listed in Table 1. We re-observed the two sources ID5340 and 6309 with one NOEMA tuning setup in January and February 2021 (ID: W20DM, PIs: E. Daddi & S. Jin) in order to obtain solid [CI](1-0) detections. ID7549 is also well-detected at 2mm continuum with NOEMA (ID: W20DM, PI: S. Jin).

ID12646 was firstly followed up with NOEMA 3mm line scans in the winter of 2017 (IDs: W17EK, PI: D. Liu); however, only one line was detected at 89.7 GHz. As the line might have been CO(6-5) at z = 6.7, a NOEMA 1mm observation (ID: W18EY, PI: D. Liu) was carried out to search for the [CII]158 um line, but no line was detected. In the winter of 2019, this galaxy was observed with a NOEMA 3mm setup (ID: W19DQ, PI: D. Liu).

The NOEMA data were reduced and calibrated using IRAM GILDAS software packages. We then produced uv tables to perform further analysis with IRAM GILDAS working on the uv space (visibility) data. The resulting synthesized beam at 3mm is in most cases rather elongated (Table 1) and no source is expected to be resolved given the resolution from the compact NOEMA configurations used. We therefore extracted the spectra by fitting a point source model in the uv space at the known spatial positions from the high-resolution ALMA imaging. The NOEMA spectra are shown in Fig. 2 and the line measurements are shown in Table 2.

Fig. 2. NOEMA and ALMA (ID9316) 3mm spectra. Detected lines are highlighted in green with their names labeled. The red line shows a power-law fit to the continuum that increases with frequency ν3.7.

2.2. ALMA

The above-mentioned four NOEMA sources in COSMOS have public imaging data in the ALMA archive at 345 GHz (ID: 2016.1.00463.S, PI: Y. Matsuda), in which two of them are also observed at 210 GHz (ID: 2016.1.00279.S; PI: I. Oteo). The four galaxies are solidly detected in the dust continuum and their fluxes were accurately measured (Table 3). No serendipitous line was found in the ALMA data cubes.

ID9316 was observed with five ALMA setups at 3mm, continuously covering frequencies from 89.5 to 103.5 GHz. We processed each setup by reproducing the observatory calibration with their custom-made script based on the Common Astronomy Software Application package (CASA; McMullin et al. 2007). We then exported the data into uvfits format to generate uv tables to perform further analysis with the IRAM GILDAS working on the uv space (visibility). The spectrum of each setup is extracted in the uv space on the position of the continuum peak using a Gaussian kernel determined from the continuum and kept fixed in size and position in each spectral channel. The final spectra, as shown in Fig. 2, are produced by stacking the spectra of the five setups. We measure a redshift of z = 4.072 via detecting multiple lines, including CO(4-3), [CI](1-0), and CO(5-4). The 1mm and 870 μm ALMA data of ID9316 were well processed as part of the A³COSMOS catalog (Liu et al. 2019a). We therefore directly adopt their measurements.

As a sanity check, we compared the ALMA 870 μm photometry for all COSMOS sources with the SCUBA2 850 μm photometry from the super-deblended catalog, finding that they are in excellent agreement (see a similar comparison, with consistent results, in Jin et al. 2018, Figs. 18 and 19). The source ID3117 has ACA (Atacama Compact Array) 630 μm observations (2019.1.01832.S, PI: J. Zavala), and its 630 μm flux from our SED fitting also agrees well with the ACA photometry. Together, these checks confirm the expected high fidelity of the flux measurements from the super-deblended technique.

We emphasize that we analyzed both NOEMA and ALMA data in uv space using the same algorithms that were applied in Jin et al. (2019), Puglisi et al. (2019), and Valentino et al. (2018), Valentino et al. (2020a,b). FIR/(sub)mm photometry is measured using the same Super-deblending technique (Jin et al. 2018; Liu et al. 2018). The consistencies therefore allow us to directly compare results in this work to these previous studies and to merge the samples.

3. Results

3.1. Line detection and redshift identification

Adopting the line-search algorithm from Jin et al. (2019), we blindly searched for lines by sweeping through the full 3 mm continuum-subtracted spectrum. We present their NOEMA and ALMA 3mm spectra in Fig. 2 with detected transitions highlighted in green, and the line measurements are listed in Table 2. Every galaxy is detected in at least two transitions, where the brightest lines are detected with 6–14σ and secondarily bright ones are detected at 4.5–6.0σ. We thus confirm redshifts of all sources to be at z = 3.545–4.147. ID3117 has the lowest redshift in this sample, of namely z = 3.545, which is confirmed by CO(3-2) and CO(4-3) emission lines (S/N = 9, 13). ID5340, 6309, and 7549 are confirmed at z = 3.815, 4.032, and 4.095 with CO(4-3) and [CI](1-0) emission lines, respectively. ID12646 is detected with three lines: CO(5-4), CO(4-3), and [CI](1-0) at z = 4.147. ID9316 is confirmed at z = 4.074 with ALMA, detecting CO(4-3), CO(5-4), and [CI](1-0) emission with S/N = 30, 23, and 19, respectively. An absorption of H₂O(1₁₀-1₀₁) is also detected with S/N = 6.2, which appears two times stronger than the continuum at the same frequency, similar to the strong H₂O absorption found in HFLS3 by Riechers et al. (2022). Notably, five of the six sources are detected with [CI](1-0) transition at 4.5–19σ. Given that most CI detections are in the local Universe (e.g., Jiao et al. 2017, 2019, 2021), the largest high-z [CI] sample is at z ∼ 1.2 from (Valentino et al. 2018, 2020b), and several [CI](2-1) detections at z ∼ 3 were presented by Yang et al. (2017). This [CI] sample is comparable in size to the SPT ones (Vieira et al. 2013; Bothwell et al. 2017), thus constituting one of the first sizable [CI] samples at z ≳ 4.

Remarkably, the redshifts of the five NOEMA sources are all lower than the original photometric redshifts. In the same way as the four galaxies found by Jin et al. (2019), this sample was selected as photometric z > 6 candidates but eventually confirmed at z ∼ 4. We note that the redshift of ID9316 is consistent with our own z_phot, IR = 3.9 ± 0.8 from the super-deblended catalog (Jin et al. 2018), but we keep this galaxy in the sample for completeness, distinguishing it from the remaining galaxies where appropriate. The full sample of ten galaxies is studied here as an enlarged version of the sample studied by Jin et al. (2019).

3.2. Cold dust SEDs and CMB impact

As mentioned above, the sources are confirmed at lower redshifts than the early expectations. To estimate the characteristic dust temperatures, we performed a detailed SED analysis by fitting their deblended FIR photometry and ALMA/NOEMA continuum fluxes. Following the identical SED recipes adopted in Jin et al. (2019), we fitted the FIR/(sub)mm photometry of this sample in three ways, using: (1) optically thin MBB (modified black body) models from Magdis et al. (2012); (2) optically thin MBB models from Magdis et al. (2012) accounting for the CMB impact on the dust continuum following da Cunha et al. (2013); and (3) optically thick MBB models from Magdis et al. (2012) also including CMB following da Cunha et al. (2013). The results are listed in Table 4 and optically thick SEDs are presented in Fig 3. The radio data points are not directly used for the fittings, but we extrapolated the expected radio emission based on the measured IR luminosities using the mass-stratified IR–radio correlation from Delvecchio et al. (2021). The predicted radio fluxes are in good agreement with the observations for most cases in this sample, showing little evidence for radio AGN emission (Fig. 3).

Fig. 3. SEDs of galaxies in this study. Photometry is from the super-deblended catalogs (Jin et al. 2018; Liu et al. 2018) and ALMA and NOEMA measurements. The dust SEDs are fitted with optically thick MBB models accounting for CMB impact on the continuum (da Cunha et al. 2013). The black (red) curve shows the best fit to the observed SED accounting (not accounting) for CMB. The blue curve marks the stellar component. We note that the radio data are not included in the fitting; we extrapolated a radio component to the dust SEDs using the IR–radio correlation from Delvecchio et al. (2021).

The five NOEMA sources show clearly steep Rayleigh-Jeans slopes β > 2.5 − 2.7 for optically thin MBB models; see the ${\beta }_{\mathrm{thin}}^{\mathrm{noCMB}}$ listed in Table 4. These Rayleigh-Jeans slopes are similar and even steeper than the findings in Jin et al. (2019). Benefiting from high signal-to-noise ratios (S/Ns) from the ALMA and NOEMA photometry, the constraints on the Rayleigh-Jeans slopes are more robust than those in Jin et al. (2019). However, no indisputable evidence has yet been produced in the literature for galaxies with steep Rayleigh-Jeans slope β > 2. A recent ALMA large survey revealed that the ALESS SMG sample can be characterised by β ≃ 1.9 ± 0.4 (da Cunha et al. 2021), which is consistent with the emissivity index of the dust in the Milky Way and other local and high-z galaxies, indicating little evolution in dust grain properties between high-z SMGs and local dusty galaxies. Therefore, the steep Rayleigh-Jeans slopes β > 2 appear not to be genuine in these galaxies, similarly to what is discussed in Jin et al. (2019). On the other hand, the CMB temperature is warmer at higher redshifts, reaching 14 K at z = 4. This reduces the observable dust continuum of cold dusty galaxies notably at the longest wavelengths (da Cunha et al. 2013; Zhang et al. 2016); e.g., 3mm. This suggests that the CMB has an impact on the Rayleigh-Jeans slopes and accounting for it can explain their abnormally steep observed values (Jin et al. 2019).

After accounting for the CMB effect on the continuum using da Cunha et al. (2013) models, the Rayleigh-Jeans slopes reduce to lower best-fitting values β = 2.0–2.5. We note that the best-fitting Rayleigh-Jeans slopes of ID3117, 5340, and 6309 remain at β > 2 in optically thin models even though the CMB effect has been accounted for. Their slopes can be further reduced to β ∼ 2 if we fit their SEDs using an optically thick dust template (plus the CMB effect).

In Fig. 4, we compare the dust temperatures from optically thin models to literature data. In the five NOEMA sources, their dust temperatures from optically thin models are 24–26 K, which is two times colder than the average T_d of main-sequence (MS) galaxies at the same redshift (Schreiber et al. 2018), and similar to the extreme cold sample in Jin et al. (2019). This work therefore confirms the existence of a new population of apparently cold dusty galaxies at z ≈ 4, further strengthening the conclusion drawn by Jin et al. (2019). However, their dust temperatures would be two times warmer if the dust is optically thick in FIR (Table 4), thus becoming comparably warm as typical dusty star-forming galaxies at high-z (Fig. 4-right). We present a diagnosis of the dust opacity of this sample in Sect. 4.

Fig. 4. Dust temperature Td versus redshift for galaxies in the literature and this work. Left panel: dust temperatures in optically thin models. Right panel: dust temperatures in optically thick models for galaxies in Jin et al. (2019) and this work, and in optically thin Td for remaining data points. This figure shows a population of cold dusty galaxies at z = 4–6 in optically thin models, but they might be comparably warm, as in main sequence galaxies, if they have optically thick dust.

3.3. Size measurements

The ancillary 870 μm and 1mm data from ALMA provide high-resolution imaging for our sample. We measure the sizes of the dust continuum in the uv space using 870 μm and/or 1mm observations, following the same method applied in Puglisi et al. (2019) and Valentino et al. (2020b). The size measurements are listed in Table 6. The six sources are well resolved with FWHM sizes in the range of 1.7–3.4 kpc. In Fig. 5, we compare sizes against a control sample of MS galaxies at z ∼ 1.2 (Valentino et al. 2020b; Puglisi et al. 2019) and from the GOODS-ALMA sample (Franco et al. 2020; Gómez-Guijarro et al. 2022a,b). Our sample has a dust continuum FWHM which is 2.4 ± 1.3 times smaller than that of comparable FIR luminous galaxies at z ∼ 1–2, while it has comparable FWHM sizes with the GOODS-ALMA sample that was selected from high-resolution ALMA 1mm images (resolution = 0.6″, Franco et al. 2020). Interestingly, an anti-correlation is seen in the figure, that is,

$$\begin{array}{c}\hfill log(FWHM\mathrm{size}/\mathrm{kpc})=-0.38\times \log (L_{\mathrm{IR}}/{10}^{10}{L}_{\odot })+1.42,\end{array}$$(1)

Fig. 5. Size of dust continuum vs. IR luminosity for this sample and data in the literature. Grey diamonds represent MS galaxies at z ∼ 1.2 in Valentino et al. (2020b) and blue squares mark the GOODS-ALMA sample (Franco et al. 2020; Gómez-Guijarro et al. 2022b). Green crosses show the medians of four local galaxy samples (Spirals: Bolatto et al. 2017; LIRGs: Bellocchi et al. 2022; ULIRGs: Pereira-Santaella et al. 2021); the dispersion of each sample is indicated by an error bar. The solid line shows the best fit to the high-z data points, and the dashed lines mark the 1σ scatter. Both local and high-z galaxies show an anti-correlation between dust continuum size and infrared luminosity.

with 1σ scatter of 0.24 dex, indicating galaxies with higher IR luminosity tend to have a more compact dust morphology. This correlation appears to be intrinsic rather than due to selection effects: as the local spirals, LIRGs, and ULIRGs (green crosses in Fig. 5) follow a similar anti-correlation, where high-z dusty galaxies have either larger sizes than the local ones, higher luminosities, or a combination there of. This is consistent with the findings of Tacconi et al. (2006), who showed that SMGs at z = 2–3 have larger sizes than local ULIRGs, suggesting that high-z SMGs are scaled-up versions of local ULIRGs at the maximum starburst limit.

This correlation is in opposition with the positive correlation reported in Fujimoto et al. (2017). We note that the sample of Fujimoto et al. (2017) was selected from the ALMA archive, which mixes all types of galaxies observed with variable configurations and sensitivities.

3.4. Near-infrared and IRAC photometry and stellar mass

We cross-matched this sample to the deepest wide-field near-infrared(NIR)/IRAC catalogs from COSMOS2020 (Weaver et al. 2022). Two sources, ID6309 and 7549, are found in both The Farmer and Classic versions. Because of the typically greater precision of profile-fitting, we chose to adopt photometry from The Farmer for these two sources. ID5340 is not included in the catalog derived using The Farmer because it is in a masked region where photometry was not performed. However, it is found in the Classic catalog but the IRAC channel 2, 3, and 4 measurements are missing, and so for this source we use Classic photometry (VISTA and IRAC channel 1), and IRAC channel 2–4 photometry from the S-COSMOS catalog (Sanders et al. 2007). For the remaining sources in COSMOS (including the four in Jin et al. 2019), their UltraVISTA and IRAC photometry are measured by performing new measurements with the The Farmer (Weaver et al. 2022), adopting point-like models centered on their ALMA positions. Models of known sources neighboring these positions from COSMOS2020 were re-fitted to improve deblending, which is consistent with the methods of the original COSMOS2020 catalog from The Farmer. The IRAC photometry of the GOODS-N source ID12646 is from the GOODS-Spitzer project (PI: M. Dickinson). Stellar masses, attenuation, and luminosity-weighted ages are summarized in Table 5, which are estimated using the MICHI2 code (Liu et al. 2021) by fitting the NIR+IRAC photometry with star-formation templates from BC03 (Bruzual & Charlot 2003), allowing for dust attenuation (Calzetti et al. 2000) and multiple ages. The E(B − V) values indicate high attenuation, that is, 2 < A_V < 3, if adopting the Milky Way reddening parameter R_V = 3.1 and 3 < A_V < 4.6 if adopting R_V = 5 for dense dust, which is in line with the optically thick dust suggested by FIR data (Sect. 4). The age is not reliably constrained because of the high attenuation.

Note that we are fitting stellar and dust SEDs separately rather than using an energy balance between optical and FIR. This is because dust and stellar components are often found to be spatially unassociated in high-z dusty galaxies, e.g., Franco et al. (2018), Hodge et al. (2019). Literature work has shown that in this type of galaxy the star formation from optical emission is just a tiny fraction of the total SFR (Puglisi et al. 2017; Calabrò et al. 2018). Therefore, there is no reason to assume that the dust and stellar are well mixed.

4. Evidence for optically thick dust in the FIR

4.1. Cold dust temperature against high IR surface brightness T_d–Σ_IR

Figure 6 shows that the new sample has cold dust temperatures and high IR surface densities that are more extreme than those of the other dusty galaxy samples. By performing a bootstrapped Kolmogorov-Smirnov test in T_d–Σ_IR space, we found a low probability of similarity (median = 2.4 × 10⁻³, semi-interquartile = 3.5 × 10⁻³) between this sample (including Jin et al. 2019) and that of Simpson et al. (2017), indicating that these galaxies belong to a new population and are significantly different from general high-z submm galaxies in the literature. Formally, these measurements are inconsistent with the Stefan-Boltzmann law even when accounting for their uncertainties and error bars. The Stefan-Boltzmann (S-B) law applies to black body emission, and acts as an upper envelope to the measurable SFR density in real galaxies. It has been shown empirically to confine the range of observations of previously known samples of distant dusty galaxies; see for example Simpson et al. (2017). However, the Stefan-Boltzmann law is only valid for black body emission that is optically thick at all wavelengths, and it is unclear whether or not grey bodies like galaxies could break the Stefan-Boltzmann envelope.

Fig. 6. IR luminosity surface density vs. dust temperature for dusty galaxies in the literature and this work. Left and right panels show dust temperatures from optically thin and thick models, respectively. Blue squares mark data points from Simpson et al. (2017); red diamonds show the ALESS sample in Hodge et al. (2019). The AS2UDS sources (Stach et al. 2019) are shown in open circles, which are calculated using results from Gullberg et al. (2019) and Dudzevičiūtė (2020). The Stefan-Boltzmann (S-B) law is shown as a solid curve, representing the expected limit for a black body. The dot-dashed curve marks the 10% systematic bias of Td, thin, where above the curve Td would be underestimated more than 10% by thin models. Dashed curves show the expected relations for general models with different λeff as labeled. This figure provides a diagnostic method identifying optically thick dust.

Here we demonstrate that the Stefan-Boltzmann law is a strict upper limit that should not be violated in any case. The Stefan-Boltzmann law states that the emission of a black body depends only on its size and temperature T. The emission is optically thick at all wavelengths, and the emitted energy is E = AσT⁴. For the case of a grey body, it has E′=AϵσT⁴ where ϵ is the emissivity (ϵ = 1 − e^−τ, τ = (ν/ν₀)^β), and E′< E because ϵ < 1. We integrated the formula of grey bodies with different effective wavelengths λ_eff from 1 μm up to 1E7 μm, that is, from optically thin to black body. The resulting Stefan-Boltzmann law for grey body can be described as log L_IR = a + b × log T_d in cgs units, and the coefficients are listed in Table 7 for different λ_eff.

When a galaxy is thick at all wavelengths, it emits like a black body (L_IR ∝ T⁴) rather than a grey body (L_IR ∝ T^4 + β in the thin limit).

In Fig. 6, we show the Σ_IR − T_d for different galaxy samples and the normalized grey-body models with different λ_eff (dashed curves). The thicker dust models tend to approach higher IR surface density with colder T_d, but they are limited at the Stefan-Boltzmann law (${\mathrm{\Sigma }}_{\text{IR}}=5.5\times {10}^5{T}_{\text{d},\text{thin}}^4$) where λ_eff is infinite. The SFR densities derived from global size measurements are lower limits to the actual Σ_IR. It is therefore noticeable that most of the galaxies in our new sample explicitly violate the Stefan-Boltzmann envelop. In order to quantify the bias introduced by the misuse of thin models, we fit thick dust templates using thin models and compare the intrinsic temperatures to the measured ones T_d, thin. The left panel of Fig. 6 shows the 10% systematic bias of T_d, thin as a dot-dashed curve; for any data points consistent with or above the curve, their T_d will be underestimated by at least 10% by thin models. The expression of the dot-dash curve is ${\mathrm{\Sigma }}_{\mathrm{IR}}=1.5\times {10}^5{T}_{\mathrm{d},\mathrm{thin}}^{4.21}$, providing a useful limit for thin models; above this curve the systematic bias of T_d, thin becomes significant and optically thick dust is needed. On the other hand, using T_d in optically thick models (Fig. 6-right) resolves the tension, as the galaxies show a distribution that appears comparable with what was seen for other regular dusty galaxy samples; for example, the SCUBA2 sources (Simpson et al. 2017) and the ALESS sample (Hodge et al. 2019). Therefore, the Stefan-Boltzmann law disproves the extremely cold T_d from the thin model and indicates optically thick dust.

In the right panel of Fig. 6 we also compare with the AS2UDS sample (Stach et al. 2019) with quantities calculated by Gullberg et al. (2019) and Dudzevičiūtė (2020), finding that the AS2UDS sample has more scatter violating the Stefan-Boltzmann limit. We note that the SFR and T_d of the AS2UDS are derived using a different method from this work, namely MAGPHYS optical+FIR+radio SED fitting with optical-FIR energy balance; this could lead to systematic errors.

4.2. High star formation efficiency from [CI](1-0)

Dust temperature can be expressed in terms of the average intensity of the radiation field ⟨U⟩, according to the relation: ⟨U⟩=(T_d/18.9)^6.04 (Magdis et al. 2012). For the parameter ⟨U⟩, the following relation applies ⟨U⟩∝SFE/Z (i.e., star formation efficiency SFE over metallicity Z, Magdis et al. 2012). The ⟨U⟩ measurements for our cold sample are approximately ten times lower than those of MS galaxies at z = 4 from Magdis et al. (2017), or about 3.3 times lower than those of MS galaxies at z ∼ 2 (see Fig. 4-right in Jin et al. 2019). Hence, in Jin et al. (2019) we discussed several interpretations for the cold dust temperatures and low ⟨U⟩ values: (1) the galaxies have optically thick dust in FIR, which are hot in the cores but appear cold for observers; (2) the galaxies have very low SFEs (i.e., long depletion timescales) and/or (3) the ISM is cooling efficiently due to very high metallicities.

Given that five of the six sources in this work have been detected with [CI](1-0) and all FIR SEDs are well constrained, we can directly measure their SFEs to disentangle the above SFE − Z degeneracy. We estimate their molecular gas masses by adopting the empirical calibration from Valentino et al. (2018), assuming a CI excitation temperature T_exc, CI = T_d, thick. SFRs are converted from IR luminosity from the FIR SED fitting, that is, SFR = L_IR × 10⁻¹⁰. In Table 6, we list gas depletion times τ_dep = M_mol, CI/SFR, the inverse of SFEs. We note that the IR luminosity is invariant upon optical depth, and the [CI](1-0) emission is almost always optically thin (see specific discussion in Jin et al. 2019 for these cold galaxies), and therefore the measurements of gas depletion times and SFEs are robust regardless of dust optical depth. We show molecular gas masses and SFRs in the top panel of Fig. 7, and compare them to scaling relations for typically star-forming and starbursting galaxies at z ∼ 2 (Sargent et al. 2014). The two sources with the highest [CI](1-0) luminosities, ID6309 and 7549, are in line with the MS scaling. The remaining galaxies fall between the MS and starburst scaling relations, showing SFEs higher than MS galaxies by a factor of 1.8–2.5. Even if adopting the canonical CI-to-gas-mass conversion from Weiß et al. (2005), which scales up the M_gas, CI by a factor of two, the SFEs would still be in line with the MS scaling. Hence, there is no evidence for anomalously low SFE.

Fig. 7. Top: SFR vs. molecular gas mass derived from CI. Solid and dashed lines mark the relations for MS and starbursting galaxies at z ∼ 2 in Sargent et al. (2014). Bottom: SFR vs stellar mass. Black line shows the MS at z = 4 from Schreiber et al. (2015), and the blue dashed line marks three times the MS. Our sample galaxies show high SFEs with respect to MS galaxies and appear to be starbursts on the SFR–M* diagram.

Metallicity might be responsible for the anomalously cold dust. To explain these ⟨U⟩ measurements together with the observed SFEs, the metallicity Z would need to be between six to eight times higher than that seen in MS galaxies at z ∼ 2. This appears unrealistic for galaxies at z ∼ 4. However, in principle, a much higher metallicity could also affect the [CI](1-0) luminosity-to-gas-mass conversion factor, for example, CI/H₂ abundance scales linearly with metallicity Z (Heintz & Watson 2020). This would reduce the requirement for excess metallicity over MS galaxies to factors of two to three, but it would still be problematic for these gas-rich galaxies. Nevertheless, we note that dynamical constraints discussed in the following section rule out an anomalous [CI](1-0) luminosity-to-gas-mass conversion factor. Therefore, the low SFE/high-Z explanation for cold dust is strongly disfavored, and the optically thick dust becomes the only rational explanation left.

We note that recent studies questioned the molecular gas mass estimates from CI. For example, Izumi et al. (2020) found that [CI](1-0) luminosity is enhanced by a factor of ten by type 1 AGN in a local galaxy. Similar cases are also found in other local galaxies (Liu et al. 2022). Such extreme cases occur only in the very central region (∼300 pc) around luminous AGN where there is non-local thermodynamic equilibrium (non-LTE). However, such bias would be diluted in unresolved galaxies or starbursting cores, significantly reducing the impact on global molecular gas estimates. This is also supported by a study of CI+CO in a local starburst galaxy by Salak et al. (2019), who found [CI](1-0)/CO(2–1) is uniform within statistical errors in the central 1 kpc, and suggested that [CI](1-0) luminosity can be used as a CO-equivalent tracer of molecular gas mass. More recently, Papadopoulos et al. (2022) studied line ratios of the two CI transitions of 106 galaxies and found their average excitation conditions to be strongly subthermal. This non-LTE excitation of the CI lines can add uncertainty in CI-based molecular gas mass estimates, but is only effective when the [CI](2-1) line is used, while the [CI](1-0)-based gas mass is substantially less affected (Papadopoulos et al. 2022). Therefore, [CI](1-0) is still a good gas tracer for galaxies in this work, as supported by the dynamical estimates in the following sections.

In the bottom panel of Fig. 7, we compare this sample to the MS at z = 4 (Schreiber et al. 2015). Intriguingly, all these cold galaxies have high sSFRs (2–30 × MS), and therefore clearly belong predominantly to the starburst population.

4.3. Gas masses from optically thin dust models are overestimated

Jin et al. (2019) argued that the dust mass of the cold galaxy ID85001674 is overestimated when derived using optically thin dust models, as it leads to an abnormally high gas mass that cannot be reconciled with the gas mass derived from [CI](1-0) emission and from dynamics. On the other hand, the dust mass would be three times smaller if the dust were optically thick in the FIR and the appropriate thick models were used, bringing the gas mass into agreement with the gas mass derived from [CI](1-0) and from dynamics. Here the same test can be performed for a larger sample thanks to the larger numbers of [CI](1-0) detections and resolved submm sizes (e.g., Jiao et al. 2017, 2019, 2021). We adopt the same approach for this sample, and estimate gas mass in three ways:

(1) The gas mass M_gas, dust is derived from dust mass assuming a gas-to-dust ratio, G/D, depending on metallicity. The dust masses are derived from SED fitting (Table 4), for both the thick and thin dust cases. For G/D, the metallicities are derived using the FMR relation (Magdis et al. 2012). We added a 20% uncertainty to the derived gas mass estimates to account for variations in the assumed metallicity (see Jin et al. 2019 for examples of the effect).

(2) The gas mass M_gas, dyn is derived from the dynamical masses of the galaxies, which are based on CO line widths and submm ALMA sizes, following Daddi et al. (2010) and Coogan et al. (2018); see Table 6. The dynamical masses are converted to total gas mass using the relation: M_gas, dyn = M_dyn, 2Re × 0.8 − M_*. The factor 0.8 is used to account for dark matter (Daddi et al. 2010). The inclinations are measured from the sizes of major and minor axes, which are obtained by fitting elliptical models to the ALMA emission of the galaxies¹. The statistical dispersion in the average is propagated to the dynamical mass error.

(3) The gas mass M_gas, CI is derived from [CI](1-0) emission. To account for the small dependence on CI excitation temperature, we assume that CI excitation temperature equals to the dust temperature estimated from thick models (Cortzen et al. 2020). The molecular gas mass is derived using Eq. (2) in Valentino et al. (2018), as listed in Table 6. We note that high-J CO transitions are impractical for constraining molecular gas mass due to uncertain excitation (Liu et al. 2015; Daddi et al. 2015) and can be significantly suppressed by thick dust (Papadopoulos et al. 2010). The [CI](1-0) transition is optically thin in all reasonable conditions, and is therefore an optimal tracer of molecular gas.

In Fig. 8, we compare the three gas mass derivations. The gas masses from optically thin dust, M_gas, dust, are explicitly (3 to 5 times) higher than those derived from dynamical mass M_gas, dyn for three galaxies, and 1.4–4 times higher than gas mass from [CI](1-0) M_gas, CI for five galaxies. In the left panel of Fig. 8, only the less extreme galaxy ID9316 falls well on the identity line of M_gas, dust − M_gas, dyn, while the other source ID12646 appears on the identity line but with large uncertainty on M_gas, dyn. In terms of M_gas, dust − M_gas, CI comparison, ID12646 is the only with M_gas, dust agreeing well with M_gas, CI. Both comparisons indicate that the dust masses from optically thin dust models are overestimated by factors of 2–3 on average. On the other hand, optically thick models predict dust masses that are indeed lower than the thin ones by just about the same factors of 2–3 (Fig. 8-right) and produce estimates that are fully consistent on average with both CI and dynamical derivations.

Fig. 8. Gas mass derived from dust mass vs. gas mass for this sample and the z = 3.62 source in Jin et al. (2019). Left panels: comparison of gas mass from optically thin dust mass to gas mass from dynamical and [CI] estimations (Mgas, dyn, Mgas, CI), respectively. Right panel: gas mass derived from optically thick dust to Mgas, dyn (red dots) and Mgas, CI (blue squares). Solid line shows the 1:1 identity relation and the dashed line is three times the identity line. Optically thin dust models significantly and systematically overestimate the dust mass of galaxies in this work, indicating that they have optically thick dust.

We note that the dynamical masses are calculated using high-J CO(4-3) velocity widths, while sizes are from 1mm dust continuum; that is, we are assuming that high-J CO gas and dust emission are equivalent. Although in some bright SMGs the low-J CO(1-0) sizes can be much larger than the high-J ones (e.g., Ivison et al. 2011; Riechers et al. 2011), the dust continuum sizes in typical distant star-forming galaxies are found to be consistent with the high-J CO sizes (Puglisi et al. 2019), and the high-J CO accurately reflects the star forming gas (Daddi et al. 2015; Liu et al. 2015) that emits from the same region as the dust. We were only able to verify this equivalence for the ALMA 3mm-detected ID9316 in our sample, and find that its CO(4-3) size (0.39″ ± 0.07″) is in agreement with the dust continuum size (Table 6).

4.4. High dust opacity at FIR

Figure 9 shows direct estimates for the dust opacity at 100 μm rest-frame, that is, τ_100 μm = κρR_e, where for κ we are adopting the value from Jones et al. (2013) and ρ is the dust volume density assuming spherical geometry. Similarly to what was discussed in terms of the Stefan-Boltzmann law, because the galaxy might actually be composed of clouds, the opacity derived in this way is an upper limit. For our sample, we show τ_100 μm derived for both optically thin and thick cases. We directly adopted the dust masses (optically thin) and sizes from Valentino et al. (2020a), as a comparison sample. These were measured in a way that is consistent with the method used in this work. As shown in Fig. 9, a correlation is seen between τ_100 μm and Σ_SFR following τ_100 μm = 3 × 10⁻²Σ_SFR. This correlation is in fact probably driven by the correlation between dust mass and IR luminosity.

Fig. 9. Dust opacity τ100 μm vs. SFR surface density ΣSFR. The dashed line shows a fit to the data points under the optically thin assumption, i.e., τ100 μm = 3.1 × 10−2ΣSFR. In both the optically thin and thick models, our sample shows high dust opacity at FIR wavelengths.

Remarkably, the τ_{100 μm, thin} of our galaxies are all greater than 1, even for the optically thick case, which already suggests that the FIR emission at 100 μm is generally optically thick, and sometimes very much so. In the optically thick case, as shown by the blue squares, the opacity is systematically lower simply because the dust masses are lower, but it is still high (0.5 < τ_100 μm < 3.0 with a median of 1.1). This further demonstrates the need for optically thick dust in these galaxies.

5. Discussion

5.1. Caution for gas mass estimates based on dust continuum

Nowadays, as an alternative approach to CO observations, the Rayleigh-Jeans dust continuum emission, is widely used to estimate dust and molecular ISM masses at high redshift (Eales et al. 2012; Magdis et al. 2012; Scoville et al. 2016, 2017a). However, dust emission is presumed to be optically thin in literature studies, and calibrated relations are based on this hypothesis.

This assumption will work for typical galaxy populations, but would become problematic for dusty galaxies with compact optically thick dust similar to those of our sample. As already discussed in Sect. 4, their M_dust, M_gas, and SFEs can be misjudged by the simple Rayleigh-Jeans flux-to-gas conversion if using optically thin fits. The dust masses in this work are overestimated by 2.2 ± 0.5 using optically thin fits, and simulations indicate that such overestimation can even be up to a factor of 20 when accounting for variable emissivity (Fanciullo et al. 2020). Therefore, conventional gas conversion based on dust continuum should be treated with caution particularly for the most luminous systems.

Moreover, the Rayleigh-Jeans slopes of cold dusty galaxies are more substantially impacted by the CMB at high-z, which further complicates the dust and gas mass estimates. According to the SED results with and without CMB (the red and black curves in Fig. 3), the 3mm continua can be lower than the intrinsic fluxes by a factor of 1.3–1.6. This indicates that models that do not account for the impact of CMB would bias the dust and gas masses in apparently cold galaxies, this time towards underestimates. Several literature studies are using 870 μm and 1mm continuum observations to constrain gas masses of galaxies at intermediate redshifts for convenience. It is true that the CMB impact is much smaller at these wavelengths, but the dust mass measurements are less accurate closer to the peak of the SED. Going to longer wavelengths like 2mm or 3mm can provide potentially more reliable derivations, but this only works when properly accounting for the CMB.

5.2. A procedure to identify optically thick sources

Based on our results, we suggest using the following procedure in order to check whether optically thick models should be considered for a given source or sample.

(1) Calculate the dust opacity (upper limit) τ_100 μm and SFR surface density Σ_SFR (Fig. 9). The threshold above which one should seriously consider optically thick dust is τ_100 μm ≳ 1 and/or Σ_SFR ≳ 20 M_⊙ yr⁻¹ kpc⁻².

(2) If a line width measurement is available (e.g., from CO), check the dynamical estimate of gas masses and compare them with gas masses derived from dust using optically thin models. If the latter have substantially larger estimates than the former, this points towards optically thick IR dust. If CI is available, do the same with gas estimates from CI. We note that CI gas estimation also depends on CI/H2 abundance.

(3) If a dust temperature T_d can be measured from the SED, compute its location in a Σ_IR − T_d diagram (Fig. 6) and compare this to the Stefan-Boltzmann limit. Optically thick dust is needed if ${\mathrm{\Sigma }}_{\mathrm{IR}}\gtrsim 1.5\times {10}^5{T}_{\mathrm{d},\mathrm{thin}}^{4.21}$(10% bias limit) and is indispensable for galaxies violating the Stefan-Boltzmann law. Also, one can compare to the T_d versus redshift trend for normal galaxies.

(4) If a [CI] excitation temperature is available, a significantly higher value than T_d, thin from optically thin models will also suggest optically thick dust (e.g., Cortzen et al. 2020).

5.3. The abundance of optically thick and compact dusty galaxies

Given the relevant properties of the IR optically thick and compact dusty galaxy population discussed in this work, an important question to tackle is how common these galaxies are. Fully answering this question would require much more extensive follow-up spectroscopy on significantly larger samples. Meanwhile, we can attempt some estimate using the larger photometric samples.

One interesting aspect is their fraction among the general population of FIR-bright galaxies, for example at z ∼ 4. We can estimate this using the super-deblended catalogs in COSMOS and GOODS-North (Jin et al. 2018; Liu et al. 2018). Adopting the equivalent IR luminosity limit that characterises this sample L_IR > 4 × 10¹² L_⊙, there are 154 (COSMOS) and 19 (GOODS-N) galaxies found with NIR/optical photometric redshift 3.5 < z < 4.5, respectively. We assume conservatively that all these NIR-detected sources are optically thin (and thus at the correct redshift). This places a strict limit on the fraction of thick dust cases ≳4% (> 7/161) in COSMOS and ≳10% (> 2/21) in GOODS-North, which are roughly consistent given the small-number statistics.

This is a strict lower limit, particularly because our spectroscopic sample of optically thick galaxies at z ∼ 4 is not complete in absolute terms. The majority (78%) of the z > 6 candidates that we followed up were actually found to be z ∼ 4. Jin et al. (2018) presented a sample of 85 FIR-detected candidates at z > 4–7. A similar fraction of those z_phot, FIR > 6 candidates would be expected at lower redshift, with their z_phot, FIR being overestimated because of thick dust. Counting only those z > 6 candidates with L_IR > 6 × 10¹² L_⊙ would be consistent with our luminosity limit L_IR > 4 × 10¹² L_⊙ if they were instead at z ∼ 4; we can add up to 30 sources to the sample of optically thick galaxies at 3.5 < z < 4.5. This would imply a more realistic fraction of 20% IR-luminous sources at z ∼ 4 that are optically thick in their dust emission. This is similar to the thick fraction that we estimated from the Simpson et al. (2017) sample, applying the proposed Stefan-Boltzmann limit in Fig. 6.

In terms of comoving density, we have to take into account that the completeness of the super-deblended catalogs is differentially affected by blending, meaning that effectively deep limits are only reachable in those regions where crowding is less severe than average. Using an approximate reduction of 2 − 3 in the effective volume, our source counts imply a number density of order 0.5 × 10⁻⁵ Mpc⁻³. Their contribution to the cosmic SFRD would be of order 2.5 × 10⁻³ M_⊙ yr⁻¹ Mpc⁻³. This is comparable to the contribution of optically dark galaxies (Wang et al. 2019), which are also z ∼ 4 candidates; although these remain largely unconfirmed with spectroscopy.

6. Conclusions

Using NOEMA and ALMA observations, we spectroscopically confirmed redshifts for six dusty galaxies at z ∼ 4, solidly detecting [CI](1-0) and CO transitions. We investigated the nature of these galaxies together with the similar sample of four analogous objects from Jin et al. (2019). Our conclusions are as follows:

1. These galaxies show cold dust temperatures, compact morphology, and abnormally steep Rayleigh-Jeans slopes. This work confirms the presence of a galaxy population at z ∼ 4 with seemingly cold temperatures.

2. In all “cold” galaxies, we detect the impact of the CMB on their mm observables. The CMB reduces their 3mm continuum observed fluxes and results in steep Rayleigh-Jeans slopes.

3. This sample shows compact dust morphology. Combining with literature data, we find a negative correlation between dust continuum size and IR luminosity.

4. We provide multiple pieces of evidence for optically thick dust in these galaxies: (1) abnormally cold dust temperature from optically thin models with respect to high IR surface brightness, violating the Stefan-Boltzmann law; (2) they have higher SFEs than MS galaxies at z ∼ 2, which disfavors the mechanism of low-efficiency star formation with rapid metal enrichment, and they will evolve to massive quiescent galaxies if there is no external gas supply; (3) the high gas mass from the thin dust model shows explicit discrepancy with the gas mass derived from CO gas dynamics and [CI](1-0) emission, indicating that dust mass is overestimated with thin dust models; (4) high optical opacity at FIR after accounting for optically thick dust mass. These taken together confirm the presence of optically thick dust in this sample, paving the way for selection of large samples of optically thick dust in the early Universe.

5. We roughly estimated the abundance of galaxies with optically thick dust, which consists of a substantial population (at least 10%) of dusty star-forming galaxies with L_IR > 4 × 10¹² L_⊙ at z ∼ 4.

Why these dusty star-forming galaxies are so compact (and therefore optically thick in FIR) is still an open question. The high dust compactness in lower redshift samples has been ascribed to mergers (e.g., as in Puglisi et al. 2021), gas inflow, or auto-regulation from star formation compression (Gómez-Guijarro et al. 2022b). Further follow-up studies are needed to clarify the nature of compact, dusty star-forming galaxies through the highest redshifts.

Acknowledgments

This work is based on observations carried out under projects number S18DI, W20DM, W17EK, W18EY and W19DQ with the IRAM Interferometer NOEMA. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). This paper makes use of the following ALMA data: ADS/JAO.ALMA 2013.1.00118.S, 2016.1.00279.S, 2016.1.00463.S, 2016.1.00478.S and 2018.1.00874.S. ALMA is a partnership of ESO (representing its member states), NSF (USA), and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ. The Cosmic Dawn Center (DAWN) is funded by the Danish National Research Foundation under grant No. 140. G.E.M., S.J. and D.B.S. acknowledge the Villum Fonden research grants 37440, 13160 and the Cosmic Dawn Center of Excellence funded by the Danish National Research Foundation under the grant No. 140. S.J. is supported by the European Union’s Horizon 2021 research and innovation program under the Marie Sklodowska-Curie grant agreement No 101060888. J.R.W. acknowledges support from the European Research Council (ERC) Consolidator Grant funding scheme (project ConTExt, grant No. 648179). F.V. acknowledges support from the Carlsberg Foundation Research Grant CF18-0388 “Galaxies: Rise and Death”. Y.G. acknowledges the National Key R&D Program of China No. 2017YFA0402704 and National Natural Science Foundation of China grant Nos. 11861131007 and 12033004. Q.S.G. is supported by the National Key Research and Development Program of China (No. 2017YFA0402703), and by the National Natural Science Foundation of China (No. 11733002, 12121003, 12192220 and 12192222), and by the science research grants from the China Manned Space Project with No. CMS-CSST-2021-A05.

References

All Tables

All Figures

	Fig. 1. Multi-wavelength images with 15″ × 15″ sizes. The instrument used and wavelength are shown in green text in each cutout. ID5340, 6309, and 7549 have no HST F160W data available and thus are shown as blank panels.
In the text

	Fig. 2. NOEMA and ALMA (ID9316) 3mm spectra. Detected lines are highlighted in green with their names labeled. The red line shows a power-law fit to the continuum that increases with frequency ν3.7.
In the text

Fig. 3. SEDs of galaxies in this study. Photometry is from the super-deblended catalogs (Jin et al. 2018; Liu et al. 2018) and ALMA and NOEMA measurements. The dust SEDs are fitted with optically thick MBB models accounting for CMB impact on the continuum (da Cunha et al. 2013). The black (red) curve shows the best fit to the observed SED accounting (not accounting) for CMB. The blue curve marks the stellar component. We note that the radio data are not included in the fitting; we extrapolated a radio component to the dust SEDs using the IR–radio correlation from Delvecchio et al. (2021).

In the text

Fig. 4. Dust temperature Td versus redshift for galaxies in the literature and this work. Left panel: dust temperatures in optically thin models. Right panel: dust temperatures in optically thick models for galaxies in Jin et al. (2019) and this work, and in optically thin Td for remaining data points. This figure shows a population of cold dusty galaxies at z = 4–6 in optically thin models, but they might be comparably warm, as in main sequence galaxies, if they have optically thick dust.

In the text

Fig. 5. Size of dust continuum vs. IR luminosity for this sample and data in the literature. Grey diamonds represent MS galaxies at z ∼ 1.2 in Valentino et al. (2020b) and blue squares mark the GOODS-ALMA sample (Franco et al. 2020; Gómez-Guijarro et al. 2022b). Green crosses show the medians of four local galaxy samples (Spirals: Bolatto et al. 2017; LIRGs: Bellocchi et al. 2022; ULIRGs: Pereira-Santaella et al. 2021); the dispersion of each sample is indicated by an error bar. The solid line shows the best fit to the high-z data points, and the dashed lines mark the 1σ scatter. Both local and high-z galaxies show an anti-correlation between dust continuum size and infrared luminosity.

In the text

Fig. 6. IR luminosity surface density vs. dust temperature for dusty galaxies in the literature and this work. Left and right panels show dust temperatures from optically thin and thick models, respectively. Blue squares mark data points from Simpson et al. (2017); red diamonds show the ALESS sample in Hodge et al. (2019). The AS2UDS sources (Stach et al. 2019) are shown in open circles, which are calculated using results from Gullberg et al. (2019) and Dudzevičiūtė (2020). The Stefan-Boltzmann (S-B) law is shown as a solid curve, representing the expected limit for a black body. The dot-dashed curve marks the 10% systematic bias of Td, thin, where above the curve Td would be underestimated more than 10% by thin models. Dashed curves show the expected relations for general models with different λeff as labeled. This figure provides a diagnostic method identifying optically thick dust.

In the text

Fig. 7. Top: SFR vs. molecular gas mass derived from CI. Solid and dashed lines mark the relations for MS and starbursting galaxies at z ∼ 2 in Sargent et al. (2014). Bottom: SFR vs stellar mass. Black line shows the MS at z = 4 from Schreiber et al. (2015), and the blue dashed line marks three times the MS. Our sample galaxies show high SFEs with respect to MS galaxies and appear to be starbursts on the SFR–M* diagram.

In the text

Fig. 8. Gas mass derived from dust mass vs. gas mass for this sample and the z = 3.62 source in Jin et al. (2019). Left panels: comparison of gas mass from optically thin dust mass to gas mass from dynamical and [CI] estimations (Mgas, dyn, Mgas, CI), respectively. Right panel: gas mass derived from optically thick dust to Mgas, dyn (red dots) and Mgas, CI (blue squares). Solid line shows the 1:1 identity relation and the dashed line is three times the identity line. Optically thin dust models significantly and systematically overestimate the dust mass of galaxies in this work, indicating that they have optically thick dust.

In the text

	Fig. 9. Dust opacity τ100 μm vs. SFR surface density ΣSFR. The dashed line shows a fit to the data points under the optically thin assumption, i.e., τ100 μm = 3.1 × 10−2ΣSFR. In both the optically thin and thick models, our sample shows high dust opacity at FIR wavelengths.
In the text