# Discovery of four apparently cold dusty galaxies at z=3.62-5.85 in the COSMOS field: direct evidence of CMB impact on high-redshift galaxy observables

S. Jin, E. Daddi, G.E. Magdis, D. Liu, E. Schinnerer, P.P. Papadopoulos, Q.

Gu, Y. Gao, A. Calabrò

► To cite this version:

S. Jin, E. Daddi, G.E. Magdis, D. Liu, E. Schinnerer, et al. Discovery of four apparently cold dusty galaxies at z=3.62-5.85 in the COSMOS field: direct evidence of CMB impact on high-redshift galaxy observables. The Astrophysical Journal, 2019, 887 (2), pp.144. 10.3847/1538-4357/ab55d6 . hal-02279083

HAL Id: hal-02279083

https://hal.science/hal-02279083v1

Submitted on 13 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
<figure><img src="image_2.png" /></figure>

# Discovery of Four Apparently Cold Dusty Galaxies at $z = 3.62–5.85$ in the COSMOS Field: Direct Evidence of Cosmic Microwave Background Impact on High-redshift Galaxy Observables

S. Jin<sup>1,2,3,4,5</sup>①, E. Daddi<sup>2</sup>①, G. E. Magdis<sup>6,7,8</sup>①, D. Liu<sup>9</sup>①, E. Schinnerer<sup>9</sup>①, P. P. Papadopoulos<sup>10,11</sup>, Q. Gu<sup>1,5</sup>②, Y. Gao<sup>12</sup>③, and A. Calabrò<sup>2</sup>④

<sup>1</sup> School of Astronomy and Space Science, Nanjing University, Nanjing 210093, People's Republic of China; shuowen.jin@gmail.com, qsgu@nju.edu.cn

<sup>2</sup> CEA, IRFU, DAP, AIM, Université Paris-Saclay, Université de Paris, CNRS, F-91191 Gif-sur-Yvette, France

<sup>3</sup> Instituto de Astrofísica de Canarias (IAC), E-38205 La Laguna, Tenerife, Spain

<sup>4</sup> Universidad de La Laguna, Dpto. Astrofísica, E-38206 La Laguna, Tenerife, Spain

<sup>5</sup> Key Laboratory of Modern Astronomy and Astrophysics, Nanjing University, Ministry of Education, People's Republic of China

<sup>6</sup> Cosmic Dawn Center at the Niels Bohr Institute, University of Copenhagen and DTU-Space, Technical University of Denmark, Denmark

<sup>7</sup> Niels Bohr Institute, University of Copenhagen, DK-2100 Copenhagen, Denmark

<sup>8</sup> Institute for Astronomy, Astrophysics, Space Applications and Remote Sensing, National Observatory of Athens, 15236, Athens, Greece

<sup>9</sup> Max Planck Institute for Astronomy, Konigstuhl 17, D-69117 Heidelberg, Germany

<sup>10</sup> Department of Physics, Section of Astrophysics, Astronomy and Mechanics, Aristotle University of Thessaloniki, Thessaloniki, Macedonia, 54124, Greece

<sup>11</sup> Research Center for Astronomy, Academy of Athens, Soranou Efesiou 4, GR-115 27 Athens, Greece

<sup>12</sup> Purple Mountain Observatory/Key Laboratory for Radio Astronomy, Chinese Academy of Sciences, 8 Yuanhua Road, Nanjing 210034, People's Republic of China

Received 2019 May 31; revised 2019 November 5; accepted 2019 November 7; published 2019 December 17

## Abstract

We report Atacama Large Millimeter Array (ALMA) observations of four high-redshift dusty star-forming galaxy candidates selected from far-infrared (FIR)/submillimeter observations in the COSMOS field. We securely detect all galaxies in the continuum and spectroscopically confirm them at $z = 3.62–5.85$ using ALMA 3 mm line scans, detecting multiple CO and/or [C I] transitions. This includes the most distant dusty galaxy currently known in the COSMOS field, ID85001929 at $z = 5.847$. These redshifts are lower than we had expected, as these galaxies have substantially colder dust temperatures (i.e., their spectral energy distributions peak at longer rest-frame wavelengths) than most literature sources at $z > 4$. The observed cold dust temperatures are best understood as evidence for optically thick dust continuum in the FIR, rather than the result of low star formation efficiency with rapid metal enrichment. We provide direct evidence that, given their cold spectral energy distributions, cosmic microwave background (CMB) plays a significant role in biasing their observed Rayleigh-Jeans (RJ) slopes to unlikely steep values and, possibly, reducing their CO fluxes by a factor of two. We recover standard RJ slopes when the CMB contribution is taken into account. High-resolution ALMA imaging shows compact morphology and evidence for mergers. This work reveals a population of cold dusty star-forming galaxies that were underrepresented in current surveys and are even colder than typical main-sequence galaxies at the same redshift. High FIR dust optical depth might be a widespread feature of compact starbursts at any redshift.

*Unified Astronomy Thesaurus concepts: Dust continuum emission (412); Far infrared astronomy (529); Starburst galaxies (1570); Luminous infrared galaxies (946); Submillimeter astronomy (1647); High-redshift galaxies (734); Cosmic microwave background radiation (322); Interstellar medium (847); Galaxy evolution (594); Galaxy formation (595)*

## 1. Introduction

Dust temperature is an important parameter of the interstellar medium (ISM) of galaxies, encoded in the shape of their infrared spectral energy distribution (SED; Draine & Li 2007). It is strongly correlated with the mean intensity of the radiation field $\langle U \rangle$, which in turn can be used to infer the metallicity-weighted star formation efficiency (i.e., SFE/Z; Magdis et al. 2012).

It is now well established that for main-sequence (MS) galaxies,<sup>13</sup> which dominate cosmic star formation at all observed redshifts (Rodighiero et al. 2011; Schreiber et al. 2015), the dust temperature is rising with redshift up to at least $z \sim 4$. From the analysis of infrared SEDs at mid-infrared to millimeter wavelengths for individual MS galaxies and for stacked ensembles, Magdis et al. (2012) were first to measure

the evolution of dust temperature with redshift, finding that MS galaxies have more intense radiation fields and thus warmer temperatures as the redshift increases from $z = 0$ to 2. By stacking far-infrared and submillimeter data from the Herschel Space Observatory, Magnelli et al. (2014) studied the evolution of the dust temperature of MS galaxies up to $z \sim 2$, largely confirming these results. Later on, Béthermin et al. (2015) found that the mean intensity of the radiation field $\langle U \rangle$ increases with increasing redshift up to $z = 4$ in MS galaxies, following the trend $(1+z)^{1.8}$ (somewhat faster than in Magdis et al. 2012) and consistent with models that account for the decrease in the gas metallicity with redshift. Using both individual detections and stacks of Herschel and Atacama Large Millimeter Array (ALMA) imaging, Schreiber et al. (2018) also found a consistent trend of increasing dust temperature with redshifts up to $z \sim 4$ for MS galaxies. The typical scatter across the MS at fixed redshift is about 0.2–0.3 dex in terms of the intensity of the radiation field $\langle U \rangle$ and 10%–13% in terms of $T_{\text{dust}}$ (Magdis et al. 2012; Schreiber et al. 2018). Similar trends of

<footnote n="13"><sup>13</sup> Here and in the rest of the paper we are discussing massive galaxies, above $10^{10} M_{\odot}$.</footnote>
dust temperature have been reported by simulations (Liang et al. 2019).

On the other hand, starburst galaxies (SBs) have a different trend of dust temperature with redshift with respect to MS galaxies and show an almost constant $\langle U \rangle \sim 30$ from $z=0$ to $z \sim 2$ and perhaps further (Béthermin et al. 2015; Schreiber et al. 2018), when fitted with optically thin dust models like MS galaxies. These SBs are galaxies situated substantially above the MS and are typically (U)LIRGs in the local universe. At moderate redshifts (up to $z \sim 1–2$), SB galaxies are observed to have hotter dust and higher star formation efficiency (SFE) with respect to MS galaxies. Ensemble average dust temperatures for SBs are not currently measured with good accuracy at $z > 2$. However, quite intriguingly, the extrapolation of the trends found by Béthermin et al. (2015) suggests that dust in typical SBs might become colder than MS galaxies starting somewhere at $z > 2.5$, and possibly all the way to higher redshifts if the dust temperature–redshift trends discussed above persist in the earlier universe. It is currently unclear what could cause such a reversal from the local universe, where ultraluminous infrared galaxies (ULIRGs) are much warmer than spirals and display much higher star formation efficiencies (SFEs), as a result of the merger process (e.g., Renaud et al. 2018, 2019).

A case of a fairly cold starburst at high redshift has already been observed, individually, in the GN20 galaxy ($z = 4.05$; Daddi et al. 2009; Tan et al. 2014). Its dust temperature has been estimated as $T_{\text{dust}} = 33$ K (Magdis et al. 2012), colder than the average dust temperature of MS galaxies at the same redshift. Still, most superluminous dusty systems at $z > 4$ are typically found to show fairly hot dust temperature ($T_{\text{dust}} \sim 40–70$ K; e.g., Riechers et al. 2013; Smolčić et al. 2015; Riechers et al. 2017; Pavesi et al. 2018). While dusty galaxies colder than GN20 have not been seen yet at $z > 4$, it is unclear whether the GN20 case is just a curiosity or the tip of the iceberg for a cold galaxy population missing in current studies. Thus, investigating galaxy SEDs on larger samples at $z > 4$ is fundamental to answering this question. Although hundreds of square degrees have been mapped at far-infrared (FIR)/(sub) millimeter wavelengths to sufficient depths that should allow detection of dusty star-forming galaxies (DSFGs) with a roughly fixed star formation rate (SFR) threshold up to $z \sim 5–10$ (particularly from Herschel SPIRE and submillimeter/millimeter ground-based detectors; Zavala et al. 2018), only a handful of sources have been spectroscopically confirmed to lie at $z > 5$ (Capak et al. 2011; Walter et al. 2012; Vieira et al. 2013; Smolčić et al. 2015; Riechers et al. 2017; Pavesi et al. 2018), and only three at $z > 6$ (Riechers et al. 2013; Fudamoto et al. 2017; Strandet et al. 2017; Zavala et al. 2018). These samples were preselected by their red colors in FIR bands (Riechers et al. 2013, 2017) and/or strong detections at (sub)millimeter wavelengths (e.g., Walter et al. 2012; Strandet et al. 2017). The SED shapes and FIR colors of galaxies are subject to the well-known degeneracy between dust temperature and redshift, so that the same FIR SED can be reproduced either by relatively hot dust continuum emission at high redshift or by intrinsically colder dust at lower redshift. Nevertheless, the spectroscopically confirmed samples at $z > 4–5$ are all displaying fairly hot dust content, and cold dusty samples similar to GN20 have not been reported yet.

The persistent sparsity of these very high-$z$ submillimeter galaxy (SMG) samples is likely due not only to the intrinsic

rarity of massive dusty galaxies in the early universe, going along with the decrease in the SFR density to early times (Lilly et al. 1996; Madau et al. 1998; Madau & Dickinson 2014; Liu et al. 2018), but also to incompleteness, i.e., missing detections of more typical objects at faint fluxes in heavily blended FIR/ (sub)millimeter images. To counter these problems, a number of new-generation FIR catalogs have been built (Béthermin et al. 2010, 2012, 2015; Roseboom et al. 2010; Elbaz et al. 2011; Lee et al. 2013; Safarzadeh et al. 2015; Hurley et al. 2017), which should allow in principle selection of DSFGs down to lower luminosities and up to highest redshifts. Notably, we have been developing new “super-deblended” catalogs (Jin et al. 2018; Liu et al. 2018) that provide state-of-the-art FIR photometry with well-behaved quasi-Gaussian uncertainties while limiting as much as possible the effects of blending from the poor IR point-spread functions (PSFs) of current ground-based facilities. Our super-deblended catalogs detect a substantial number of $z > 4$ FIR-detected galaxy candidates for which only a small fraction of spectroscopic identifications are as yet available. Given the less extreme luminosities of these candidates, it is more difficult to identify their redshifts spectroscopically. Detection of rest-frame UV/optical lines is time-consuming (Chapman et al. 2005; Capak et al. 2011) and infeasible for DSFGs that are invisible in *Hubble Space Telescope (HST)* imaging (e.g., Walter et al. 2012; Riechers et al. 2013). IR/(sub)millimeter fluxes and FIR lines are more easily accessible in these objects than optical features (partially by construction), so that (sub)millimeter line scans (e.g., CO, [C II]) are more efficient for confirming redshifts for IR-selected star-forming samples at $z > 5–7$ (Walter et al. 2012; Riechers et al. 2013; Vieira et al. 2013; Fudamoto et al. 2017; Strandet et al. 2017). This method is becoming more prevalent in the ALMA era, although the lower-luminosity galaxies still require comparably longer efforts for succeeding.

For cold sources in particular, there are extra difficulties stemming not only from the intrinsically fainter luminosities but also from the effects caused by the cosmic microwave background (CMB) at high redshift. The CMB, increasing in temperature proportionally to $(1+z)$, provides an additional source of heating to the ISM, thus increasing the CO luminosities for high rotational number transitions (e.g., Silk & Spaans 1997). However, taking into account the higher background against which the respective lines must be detected, it actually suppresses any such boost (Papadopoulos et al. 2000). By modeling CO line emission in starburst galaxies, Combes et al. (1999) found that the hotter CMB at high-$z$ does not help the detection of CO lines and the net effect is negative. For Milky Way (MW) type galaxies, Obreschkow et al. (2009) found that the CO lines will be dramatically suppressed by the weak contrast against the CMB at $z > 5$ and might become totally undetectable at $z \sim 7$. Da Cunha et al. (2013) quantified the CMB effects on (sub)millimeter observations of high-redshift galaxies at $z > 5$ and concluded that the inferred dust and molecular gas masses can be severely underestimated for cold systems if the impact of the CMB is not properly taken into account. Zhang et al. (2016) studied the CMB effect on the observed structural and dynamical characteristics of galaxies. They found that elevated CMB at high-$z$ can dramatically shrink the emergent continuum and low-$J$ CO brightness distributions of the cold molecular gas, thus erasing spatial and spectral contrasts between their
brightness distributions and the CMB. This can strongly affect the measurements of dust and H₂ gas distribution scale lengths and velocity fields, and thus the enclosed dynamical masses in galaxies.

It is clear that the potentially dramatic effects of the CMB, and the resulting biases against the coldest sources, will need to be considered carefully when designing and conducting observations. At the same time, studying cold dusty galaxies at higher z, discovering more examples of them, will be providing the chance to obtain observational confirmation of the predicted CMB effects, while at the same time helping to clarify the cosmic evolution of dust temperature in galaxies with less biases. This is crucial for constraining the physics of star formation and galaxy formation in the very early universe.

In this paper, we report the discovery of four remarkably cold dusty galaxies at $z = 3.62–5.85$. These galaxies were selected in the FIR/(sub)millimeter super-deblended catalog in the COSMOS field (Jin et al. 2018). We present here their spectroscopic observations with ALMA Band 3 scans, which confirm their redshifts and allow a detailed investigation of their physical properties. The paper is organized as follows. We discuss the target selection and ALMA observations in Section 2. The analysis of ALMA spectra, line identifications, and redshift estimates is discussed in Section 3. Section 4 focuses on the derivation of physical parameters for our target galaxies, with an evaluation of the impact of the CMB on observables. We discuss the implication of these results in Section 5 and provide our summary and conclusions in Section 6. We use standard cosmology (73, 0.23, 0.73) and a Chabrier initial mass function.

We emphasize here that throughout the paper we discuss the coldness of the SED as an objective, directly measured property, reflected in the wavelength location of the peak of the SED, in the rest frame. This is thus independent of any physical, model-dependent estimate of the actual effective dust temperature. In fact, a hypothesis we will address in the discussion of this paper is that many of these cold sources are only seemingly cold owing to optical depth, while their physical dust temperatures could be much higher. Regardless of the true dust temperature, the CMB effects depend only on the shape of the SED.

## 2. Sample Selection, ALMA Observations, and Data Reduction

In the FIR/(sub)millimeter super-deblended catalog in the COSMOS field, Jin et al. (2018) selected a sample of 85 $z > 4$ high-redshift candidates with significant IR detection ($S/N_{FIR+mm} > 5$ combined) from Herschel, SCUBA2, AzTEC, and/or MAMBO images (Bertoldi et al. 2007; Aretxaga et al. 2011; Cowie et al. 2017; Geach et al. 2017). We singled out four galaxies among these candidates for further ALMA follow-up, based on their highest photometric redshifts $z_{phot, FIR} > 6$ that were derived using the $z = 6.3$ HFLS3 template (Riechers et al. 2013). Among these four galaxies, ID20010161 was originally detected in the Very Large Array (VLA) 3 GHz catalog (Smolčić et al. 2017), while the others were found in residual images from the SCUBA2 850 µm maps from Geach et al. (2017), after subtracting all known sources in our super-deblending process (see details in Section 4.5 of Jin et al. 2018). The four sources have no detection in HST or in UltraVISTA images (McCracken et al. 2012; Laigle et al. 2016), while we found red IRAC counterparts in the Spitzer Large Area Survey with Hyper-Suprime-Cam (SPLASH, P. Capak et al. 2019, in preparation) imaging data sets, well aligned with their VLA/SCUBA2 positions with offsets <3". The subsequent ALMA observations with a subarcsecond beam, described in the next sections, leave no further doubts on the correctness of the IRAC counterpart identification.

We followed up the four sources with ALMA Band 3 spectroscopic scans in Cycle 5 (Project ID: 2017.1.00373.S, PI: S. Jin). The Band 3 observations were performed by combining three tunings covering 84–108 GHz with resolution of 16 MHz, which is a more efficient way than performing formal scans albeit at the price of small gaps (850 km s⁻¹) within the spectral range. The observations were carried out with the array configuration C43-6 giving a synthesized beam of ~0"×85. The four galaxies were observed with track sharing, and each source was observed for 1.3 hr of on-source time in each tuning, reaching typical rms sensitivities of 0.1 mJy per 500 km s⁻¹. We also used ALMA Band 6 imaging (Project ID: 2016.1.00279.S; PI: I. Oteo) available in the archive for two of our targets. The Band 6 observations are taken in a single tuning at 230 GHz with 56 s of integration time per source. One of our sources has also 345 GHz imaging publicly available from the program 2016.1.00463.S (PI: Y. Matsuda). The sources are very strongly detected in the continuum but show no evidence for lines at 230 or 345 GHz.

We processed the 3 mm spectra by reproducing the observatory calibration with their custom-made script based on the Common Astronomy Software Application package (CASA; McMullin et al. 2007). We converted the data into uvfits format to perform further analysis with the IRAM GILDAS tool working on the *u-v* space (visibility) data. The source positions are determined by the centroid of collapsed 3 mm (and 1 mm when available) data cubes, which mostly show highly significant continuum detections, and we used these positions for extracting spectra by fitting source models in the *u-v* space. We iteratively searched for evidence of resolved emission either in the continuum, in the emission-line candidates, or both. Only one source was found to be resolved at more than the $3\sigma$ level, as discussed below. We further extracted spectra for this source with Gaussian models, and we use point-source models for the rest. The original spectra with noise are extracted with a resolution of 16 MHz per channel. We verified that the noise estimates per channel are reliable from global $\chi^2$ statistics. We present the binned and noise-normalized spectra in Figure 1.

We also verified the FIR/(sub)millimeter photometry of the sources that resulted from SCUBA2 residual maps by obtaining new super-deblended photometry using their better determined positions from ALMA. The difference with the original photometry in Jin et al. (2018) is entirely negligible.

## 3. Redshift Confirmation

By sweeping through their full 3 mm (continuum subtracted) spectrum, we blindly search for the single emission-line feature with the highest significance for each galaxy, following the same line-searching algorithm applied in Daddi et al. (2015), Coogan et al. (2018), and Puglisi et al. (2019). This algorithm returns both the signal-to-noise ratio (S/N) and the optimum velocity range of the line, by determining the pair of starting channel and ending channel, across which the line flux is integrated, that corresponds to the lowest chance probability for
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Top panels: ALMA 3 mm spectra at a resolution of 32 MHz per channel. In order to flatten out noise variation effects (e.g., at sidebands edges), the fluxes are divided by the noise at each channel and multiplied by the average noise of the entire spectrum. Detected lines are highlighted in green, and their identifications are labeled. Gray shaded areas show frequency gaps. The blue dashed line indicates the $1\sigma$ rms noise for each bin, which is also labeled in the upper right corner, computed for a 32 MHz channel, corresponding to $100 \text{ km s}^{-1}$ at the average frequency of the spectra. The red dashed line shows a power-law fit to the continuum that increases with frequency $\nu^{3.7}$. Bottom panels: stacked line spectra after continuum subtraction. Green areas mark line ranges, the same as in the top panels, as determined by the primary line. Gray areas show the location of gaps in velocity space.</figcaption></figure>

the line (hence highest detection significance). The flux density $I_{1\text{st}}$ (Table 1) is integrated over this velocity range, which basically corresponds to the full width at zero intensity

(FWZI) of the line, which is highlighted in green in Figure 1 and shown in Table 1. Simulations and double-Gaussian fits carried out in Coogan et al. (2018) show that typically this
<table><caption>**Table 1**
Line Detections</caption><thead><tr><th>ID</th><th>z<sub>spec</sub></th><th>S/N<sub>1st</sub></th><th>I<sub>1st</sub><br>(Jy km s<sup>-1</sup>)</th><th>N<sup>EFF</sup><sub>1st</sub></th><th>P<sub>1st</sub></th><th>FWZI<br>(km s<sup>-1</sup>)</th><th>S/N<sub>2nd</sub></th><th>I<sub>2nd</sub><br>(Jy km s<sup>-1</sup>)</th><th>N<sub>2nd</sub></th><th>P<sub>2nd</sub></th><th>P<sub>c,l</sub></th></tr></thead><tbody><tr><td>85001929</td><td>5.847</td><td>5.22 (CO 6–5)</td><td>0.35 ± 0.07</td><td>3925</td><td>7.0e-4</td><td>885</td><td>3.27 (CO 5–4)</td><td>0.22 ± 0.07</td><td>12</td><td>0.013</td><td>0.9e-5</td></tr><tr><td>20010161</td><td>5.051<sup>a</sup></td><td>5.97 (CO 5–4)</td><td>0.33 ± 0.06</td><td>4188</td><td>1.2e-5</td><td>850</td><td>...</td><td>...</td><td>...</td><td>...</td><td>1.2e-5</td></tr><tr><td>85000922</td><td>4.440</td><td>4.68 (CO 5–4)</td><td>0.29 ± 0.06</td><td>5162</td><td>0.014</td><td>608</td><td>3.55 (CO 4–3)</td><td>0.17 ± 0.05</td><td>13</td><td>0.005</td><td>0.7e-4</td></tr><tr><td>85001674</td><td>3.623</td><td>6.33 (CO 4–3)</td><td>0.42 ± 0.07</td><td>4294</td><td>1.3e-6</td><td>747</td><td>4.44 ([C I] 1–0)</td><td>0.25 ± 0.06</td><td>13</td><td>1.4e-4</td><td>1.8e-10</td></tr></tbody></table>

**Notes.** S/N<sub>1st</sub>: significance of the first line; I<sub>1st</sub>: flux density of the first line; N<sup>EFF</sup><sub>1st</sub>: number of effective trials following Equation (2); P<sub>1st</sub>: chance probability of the first line; FWZI: full width at zero intensity of the lines; S/N<sub>2nd</sub>: significance of the second line; I<sub>2nd</sub>: flux density of the second line; N<sub>2nd</sub>: number of trials for the second line; P<sub>2nd</sub>: chance probability of the second line; P<sub>c,l</sub>: chance probability of two lines.

FWZI is very close to the actual FWHM of the line. As shown in Figure 1, the most significant lines in the four spectra are detected at 5σ–6σ with full velocity widths of 608–885 km s<sup>-1</sup> (Table 1). These line widths are large but comparable to those found in the literature for SMGs (e.g., Daddi et al. 2009; Riechers et al. 2013, 2017; Marrone et al. 2018).

### 3.1. Chance Probability of Emission Lines

We run extensive simulations to determine how to compute the chance probability of finding emission lines (with the same line-searching algorithm) as produced by noise fluctuations in a given spectral scan, depending on the S/N of the feature and its line width (broader lines are more unlikely to happen by chance because fewer independent realizations are available within the observed range; see also González-López et al. 2019), and also based on the full velocity range spanned by the data, the velocity sampling, and the range of acceptable line widths (from 100 to 1000 km s<sup>-1</sup> in this case). We assume well-behaved Gaussian noise in the data, as verified to hold for our ALMA spectra. We define the probability that a line detection is due to noise fluctuations as

$$P_{\text{line}} = 1 - P_{0}^{\text{EFF}}, \qquad (1)$$

where P<sub>0</sub> is the probability to have a line with S/N up to that observed in a single trial (approaching 1 for high-σ events), and the exponent N<sup>EFF</sup><sub>trials</sub> is the effective number of trials, which, based on our simulations, can be well approximated by

$$N_{\text{trials}}^{\text{EFF}} \sim 10 \frac{N_{\text{total,ch}}}{N_{\text{line,ch}}} N_{\text{line,ch}}^{0.58} \log \frac{N_{\text{line,ch}}^{\text{max}}}{N_{\text{line,ch}}^{\text{min}}}, \qquad (2)$$

where N<sub>total,ch</sub> and N<sub>line,ch</sub> are the number of channels in the entire spectrum scanned and in the recovered line, respectively, while N<sup>max</sup><sub>line,ch</sub> and N<sup>min</sup><sub>line,ch</sub> define the allowed velocity range of lines to be searched, expressed as number of channels. Notice that a naive treatment would account only for the first term (ratio of full velocity in spectrum to line velocity), while in reality the number of effective realizations is considerably larger. In fact, less than a full line-velocity shift is required to imply a new realization, and the number of effective realizations depends also on the range of velocity widths searched, as the line width is not known a priori. We emphasize that this recipe, and in particular Equation (2), holds only in the assumption that there is no active spatial search for the line position, as in our case, where the positions for the spectral extractions are all strongly constrained by the continuum

detections. When one is blindly searching for emission lines, having also to determine their spatial position, the number of effective trials would be much larger, as Equation (2) would have to be multiplied by the number of effective spatial positions searched.<sup>14</sup>

### 3.2. Measuring Redshifts for Our Sample

We report the spurious probability for the strongest line in each spectrum in Table 1 (P<sub>1st</sub>). They are below 0.1% in all cases except for the 4.7σ line of ID85000922, which has a 1.4% probability of being spurious, by itself.

Then, we searched for additional matching lines based on the redshift solutions determined by the first ones, using the same velocity range as determined by the primary line, assuming that they are CO transitions. We excluded, in fact, solutions in which the primary line was [C I] or H<sub>2</sub>O, given that stronger accompanying CO lines would have been detected in the spectra in those cases. Meaningful second lines were detected in three galaxies with significance of 3.3σ–4.4σ. We compute chance probabilities to find matching second lines (P<sub>2nd</sub>) based on the significance of each line and the number of search trials performed (N<sub>2nd</sub>). The combined probability for the two-line match is in all cases ~10<sup>-4</sup> or lower; hence, we consider these three galaxies as reliable spectroscopic confirmations. The redshifts (and line identifications) are determined as CO (5–4)/CO (6–5) at z = 5.85, CO (4–3)/CO (5–4) at z = 4.44, and CO (4–3)/[C I] (1–0) at z = 3.62 (Figure 1). In the bottom panels of Figure 1, we present zoom-ins to the weighted averages of lines detected for each galaxy after continuum subtraction. Given that the CO (4–3) line of ID85001674 partially falls in a frequency gap, we ultimately calculate its redshift and velocity width from the [C I] (1–0) line. We note that the confirming/second line for ID85001929 is located at the edge of the spectral range and only partially covered by our data. Still, the available signal within the velocity range predetermined by the first line is strong enough that its spurious probability is at 0.013. We thus maintain that this z = 5.85 identification is secure.<sup>15</sup> Therefore, ID85001929 at z = 5.85 is

<footnote><sup>a</sup> This redshift is less secure, as it is based on single-line identification (see text for details).</footnote>

<footnote n="14"><sup>14</sup> This is typically assumed to be the area within the primary beam divided by the area covered by half the synthesized beam, although dedicated simulations would be required to confirm this accurately.</footnote>

<footnote n="15"><sup>15</sup> Interestingly, some excess flux at 1.1 mm is found in its FIR SED (Figure 3), which hints at the presence of a bright [C II] 158 µm line at a consistent z ~ 5.9 redshift that is boosting the AzTEC and MAMBO photometry. However, to make an impact, this would need to be at the level of many dozens of Jy km s<sup>-1</sup>, which would imply L<sub>[C II]</sub>-to-L<sub>IR</sub> ratios at the level of local galaxies >3 × 10<sup>-3</sup>. Direct observations are required to confirm this tentative evidence.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** ALMA 1 mm (top panels) and 3 mm (middle panels) images, with contours in steps of 2σ. The size and orientation of the beam are indicated in the lower left corner. Bottom panels: color images of the four galaxies (blue: UltraVISTA $K_s$; green: IRAC 3.6 µm; red: IRAC 4.5 µm). Cyan contours show ALMA 3 mm continuum starting from 3σ in steps of 2σ.</figcaption></figure>

now the most distant known SMG in the COSMOS field, at a higher redshift than the recently found dusty galaxy CRLE at $z = 5.67$ (Pavesi et al. 2018).

The ID20010161 has a single line detection at 95.2 GHz with $5.97\sigma$ significance and a spurious probability of $\sim 10^{-5}$ (Table 1); hence, it is highly secure. However, the redshift identification is less secure, as there are multiple solutions given the lack of a second line. We disfavor the single-line identification as CO (7–6) at $z = 7.47$ given the lack of [C I] (2–1) (the formal signal within the expected velocity range gives $-2.7\sigma$ at 95.5 GHz) and H$_2$O (2$_{11}$–2$_{02}$) (0.07$\sigma$ at 88.8 GHz) lines, since the H$_2$O line is typically 0.4–1.1 of the high-J CO flux in high-$z$ lensed ULIRGs (Omont et al. 2013; Yang et al. 2016, 2019). We also exclude CO (6–5) at $z = 6.26$ because of the lack of detection of H$_2$O (2$_{11}$–2$_{02}$) ($-1.075\sigma$ at 103.6 GHz). The solution with CO (4–3) at $z = 3.84$ is also less likely owing to the lack of any [C I] (1–0) signal ($-0.7\sigma$ at 101.7 GHz), although it is not unconceivable that [C I] (1–0) might be intrinsically faint for this source, as it can reach $\sim 25\%$ of the CO (4–3) flux (e.g., Walter et al. 2011; Valentino et al. 2018). This redshift is less secure given that no confirming lines are present in its spectrum. However, we consider the most likely identification of this line as CO (5–4)

at $z = 5.05$, for which we would not expect other significant lines in the observed range. In the following we will use this redshift identification for this source. If the redshift should turn out to be instead $z = 3.84$, its dust temperature would be even lower and the conclusions derived from these sources even further strengthened.

These redshift identifications consistently confirm that no strong lines are expected in the 230 GHz spectral ranges of the Oteo data or in the 345 GHz spectral range of the Matsuda data, for either of the two galaxies for which 230 GHz or 345 GHz observations are available.

# 4. Results

## 4.1. Multiwavelength Imaging and Morphology

In Figure 2, we present ALMA (clean) images for this sample. We solidly detected continuum emission with peak significance of 9σ–14σ at 3 mm for the four sources, 13σ–19σ at 1 mm, and 29σ at 870 µm for those in which these are available (Table 2). We measured their sizes by fitting models in u-v space in the ALMA 1/0.8 mm and 3 mm collapsed data sets (also adding lines when available, but these provide only a modest contribution, being at substantially lower S/Ns). We
<table><caption>**Table 2**
ALMA Continuum Measurements</caption><thead><tr><th>ID</th><th>S<sub>3</sub> mm (mJy)</th><th>S<sub>1.3</sub> mm (mJy)</th><th>S<sub>870</sub> µm (mJy)</th></tr></thead><tbody><tr><td>85001929</td><td>0.088 ± 0.007</td><td>1.52 ± 0.12</td><td>...</td></tr><tr><td>20010161</td><td>0.115 ± 0.008</td><td>...</td><td>...</td></tr><tr><td>85000922</td><td>0.068 ± 0.007</td><td>...</td><td>...</td></tr><tr><td>85001674</td><td>0.124 ± 0.009</td><td>2.59 ± 0.14</td><td>7.48 ± 0.26</td></tr></tbody></table>

combine the different tracers for a single galaxy in the u-v space, following the procedure presented in Puglisi et al. (2019). Only the lowest-redshift galaxy, ID85001674 at $z = 3.62$, is resolved with an FWHM size of $0.′′42 ± 0.′′04$ (3.1 ± 0.3 kpc), implying an SFR surface density ($\Sigma_{SFR}$) of $22_{-9}^{+13} M_{\odot} yr^{-1} kpc^2$. The remaining sources are unresolved; we show their $2\sigma$ size upper limits (and lower limits of SFR surface density $\Sigma_{SFR}$) in Table 3. ID85000922 at $z = 4.4$ has a size limit of $<2.8$ kpc, and the other two galaxies at $z > 5$ appear to be more compact with FWHM sizes $<1.9$ kpc. The measured sizes of our ALMA sources are significantly smaller than the 3 GHz sizes reported by Miettinen et al. (2017) for a sample of 115 known SMGs in the COSMOS field, with a median size of $4.6 ± 0.4$ kpc. In contrast, our measured sizes are close to the ALMA sizes of high-z SMGs from Ikarashi et al. (2015) and Simpson et al. (2015), who report median sizes of $1.6 ± 0.14$ kpc at 1.1 mm and $2.4 ± 0.2$ kpc at 870 µm, respectively. Our ALMA sizes also agree with the 10 GHz selected sample in Murphy et al. (2017), which has a median size of $1.20 ± 0.28$ kpc.

The high compactness of this sample is unlikely to be caused by AGN activity, because the emission from any dusty torus would be negligible at FIR/(sub)millimeter wavelength, and it cannot dominate the 3 mm emission (rest frame 440–650 µm for this sample). In any case, we do not see any evidence for an AGN component in their cold SEDs (see Section 5.2), and none are X-ray detected from Jin et al. (2018). Therefore, we speculate that the compact morphology in this very highly star-forming sample is a consequence of ongoing mergers, similarly to local ULIRGs (Soifer et al. 2000; Juneau et al. 2009) and high-z compact starbursts (Puglisi et al. 2017; Marrone et al. 2018; Calabrò et al. 2019).

In the bottom panels of Figure 2, we show $K_s$+SPLASH color images for this sample. All sources have a red counterpart in SPLASH images quite consistent with the ALMA imaging; however, only two of them are entirely consistent in position and extension with the ALMA detections. The SPLASH counterpart of the $z = 5.85$ galaxy appears to be more extended than its ALMA morphology, the IRAC peak being located about $0.′′8$ to the south, where some faint extension is seen at 1 mm in the ALMA imaging as well (Figure 2, top left panel). ID85000922 at $z = 4.44$ has a very faint SPLASH IRAC 3.6 µm counterpart at its ALMA position, while the blue source with $\sim 2.′′0$ offset appears to be unrelated, having a well-constrained $z_{phot} = 0.86$ (Laigle et al. 2016).

The ID20010161 and ID85001674 galaxies are detected in the 3 GHz image with $5.0\sigma$ and $3.6\sigma$ (Daddi et al. 2017; Smolčić et al. 2017), showing infrared-to-1.4 GHz radio luminosity ratio $q_{IR} = 2.3$ and $2.4$, which agree well with the evolution of $q_{IR}(z) = (2.88 ± 0.03)(1 + z)^{0.19±0.01}$ for star-forming galaxies in Delhaize et al. (2017), while the other sources are either marginally or not detected in the radio ($<2.8\sigma$). No signal is observed in the UltraVISTA $K_s$ band (or

any shorter wavelength) for our sources, which thus could be qualified as $K_s$-band dropouts. No photometric redshift is, of course, available for these sources in near-IR or optically selected catalogs.

## 4.2. Cold Dust against the CMB

We collected super-deblended FIR/(sub)millimeter photometry together with the ALMA 0.8/1 mm and 3 mm continuum measurements for deriving their full dust SEDs. We start by fitting these SEDs using modified blackbody (MBB) models from Casey et al. (2012) with free $T_{dust}$ and $\beta$, not accounting for the effect of the CMB. We find that these galaxies are fitted with very cold dust temperatures $T_{dust} = 20–41$ K and abnormally steep Rayleigh–Jeans (RJ) slopes $\beta = 2.4–3.7$, which appear to be well determined with relatively small errors, particularly for the two objects with highly accurate 1 mm and 3 mm ALMA continuum measurements in the RJ tail. Such steep slopes had not been reported in literature previously. However, galaxies are always observed against the CMB, and the CMB effect on the continuum would be nonnegligible on dusty SEDs, especially for our systems that appear to be fairly cold for their high $z$ (da Cunha et al. 2013).

We thus use alternative optically thin MBB SEDs (Magdis et al. 2012) accounting for the CMB effect on dust continuum following the prescriptions by da Cunha et al. (2013). MBB models are not ideal to describe the Wien part of galaxies SEDs but are appropriate here, given the photometric sampling available for our galaxies. We convert the luminosity of MBB models $L_{MBB}$ to $L_{IR,8–1000 \mu m}$ by multiplying by a constant of 1.35, a median value for starburst samples (Magdis et al. 2012) based on comparison with Draine & Li (2007). In Figure 3, we show the best fit to the observed photometry with solid curves, while dashed curves indicate the intrinsic SEDs before accounting for the CMB effect on the continuum (green curve). We can clearly see that the CMB depresses the dust continuum at longer wavelengths, severely impacting its level at 3 mm (see Section 4.4), thus making the observed RJ slope steeper. In Table 3, we list derived parameters for fits with and without accounting for the CMB. In both cases we find consistently cold dust temperature $T_{dust} = 20–42$ K, which are also consistent with those derived by fitting Casey et al. (2012) models. The consistency of $T_{dust,obs}$ and $T_{dust,thin+CMB}$ shows that accounting for the CMB effect does not alter $T_{dust}$ for this sample, implying that the effects of the $T_{dust-\beta}$ degeneracy are negligible for our galaxies, within the uncertainties. This is because the CMB is affecting the observed continuum at long enough wavelengths, while it has barely any effect on the peak of the dust SED where $T_{dust}$ is measured. We verified that the impact of the CMB on the derived $T_{dust}$ is negligible out to $z = 5$; even for a galaxy with intrinsic $T_{dust} = 18$ K, the marginal overestimate of $T_{dust}$ would be less than 1 K in that case. Strikingly, the SEDs corrected for CMB effect yield $\beta \sim 2–2.2$, consistent with typical $\beta$ values found in the literature for star-forming galaxies, within the uncertainties (of about 0.2 in $\beta$).

Have we thus really directly detected the effect of the CMB on galaxy SEDs, for the first time (as far as we know)? This appears to be the case, as we argue that the $\beta$ values of order $2.4–3$ are not plausible and the CMB is indeed required to obtain a physically meaningful interpretation of the SEDs of
<table><caption>Table 3
CMB Impact on Observables</caption><thead><tr><th>ID</th><th>L<sub>IR</sub><br>(10<sup>12</sup> L<sub>☉</sub>)</th><th>β<sup>noCMB</sup><sub>thin</sub></th><th>β<sup>CMB</sup><sub>thin</sub></th><th>β<sup>CMB</sup><sub>thick</sub></th><th>T<sup>noCMB</sup><sub>dust,thin</sub><br>(K)</th><th>T<sup>CMB</sup><sub>dust,thin</sub><br>(K)</th><th>T<sup>CMB</sup><sub>dust,thick</sub><br>(K)</th><th>M<sup>noCMB</sup><sub>dust,thin</sub><br>(10<sup>8</sup> M<sub>☉</sub>)</th><th>M<sup>CMB</sup><sub>dust,thin</sub><br>(10<sup>8</sup> M<sub>☉</sub>)</th><th>M<sup>CMB</sup><sub>dust,thick</sub><br>(10<sup>8</sup> M<sub>☉</sub>)</th><th>L<sup>'obs</sup><sub>CO(5-4)</sub><br>(a)</th><th>L<sup>'CMB</sup><sub>CO(5-4)</sub><br>(a)</th></tr></thead><tbody><tr><td>85001929</td><td>12.8 ± 1.8</td><td>2.5 ± 0.2</td><td>2.2 ± 0.2</td><td>2.1 ± 0.2</td><td>40 ± 3</td><td>42 ± 3</td><td>61 ± 8</td><td>3.5 ± 0.6</td><td>4.3 ± 1.1</td><td>2.2 ± 0.9</td><td>2.0 ± 0.4<sup>b</sup></td><td>2.9 ± 5.5<sup>b</sup></td></tr><tr><td>20010161</td><td>6.2 ± 1.0</td><td>2.4 ± 0.2</td><td>2.0 ± 0.2</td><td>1.9 ± 0.2</td><td>32 ± 3</td><td>35 ± 4</td><td>40 ± 6</td><td>8.7 ± 2.5</td><td>10.4 ± 2.9</td><td>7.7 ± 2.1</td><td>1.3 ± 0.2</td><td>1.8 ± 0.3</td></tr><tr><td>85000922</td><td>5.1 ± 1.0</td><td>3.7 ± 0.6</td><td>2.3</td><td>2.1 ± 0.4</td><td>20 ± 4</td><td>30 ± 4</td><td>42 ± 6</td><td>18.3 ± 9.8</td><td>10.2 ± 1.3</td><td>7.0 ± 2.3</td><td>0.9 ± 0.2</td><td>1.3 ± 0.3</td></tr><tr><td>85001674</td><td>3.3 ± 0.8</td><td>2.7 ± 0.2</td><td>2.2 ± 0.2</td><td>2.1 ± 0.2</td><td>23 ± 2</td><td>24 ± 2</td><td>41 ± 5</td><td>29.2 ± 8.9</td><td>33.6 ± 8.1</td><td>13.0 ± 3.7</td><td>0.8 ± 0.1<sup>b</sup></td><td>1.2 ± 0.2<sup>b</sup></td></tr></tbody></table>

Notes. L<sub>IR</sub>: IR luminosity at 8–1000 µm corrected for CMB effect (dashed curve in Figure 3); β<sup>noCMB</sup><sub>thin</sub>, T<sup>noCMB</sup><sub>dust,thin</sub>, and M<sup>noCMB</sup><sub>dust,thin</sub>: RJ slope, dust temperature, and dust mass from optically thin MBB fits without accounting for CMB; β<sup>CMB</sup><sub>thin</sub>, T<sup>CMB</sup><sub>dust,thin</sub>, and M<sup>CMB</sup><sub>dust,thin</sub>: RJ slope, dust temperature, and dust mass from optically thin MBB fits accounting for CMB effect; β<sup>CMB</sup><sub>thick</sub>, T<sup>CMB</sup><sub>dust,thick</sub>, and M<sup>CMB</sup><sub>dust,thick</sub>: RJ slope, dust temperature, and dust mass from optically thick MBB fits accounting for CMB effect; L<sup>'CMB</sup><sub>CO(5-4)</sub>: CO (5–4) luminosities corrected for CMB effect assuming T<sub>exc</sub> = T<sup>CMB</sup><sub>dust,thin</sub> in LTE condition.

<footnote>a Units of 10¹⁰ K km s⁻¹ pc².</footnote>

<footnote>b The observed CO (5–4) luminosities for these two sources are extrapolated as 1.2 × L'<sub>CO(4–3)</sub> and 0.8 × L'<sub>CO(6–5)</sub>, respectively.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Dust SEDs of the four galaxies in this work. Photometries are taken from the super-deblended catalog (Jin et al. 2018), ALMA (magenta) 3 mm, 1 mm, and/or 870 µm observations. Arrows mark the 3σ upper limits. The SEDs are fitted by MBB models accounting for CMB impacts on the continuum (da Cunha et al. 2013). The solid (dashed) curve shows the best fit to the observed SED accounting (not accounting) for CMB. A radio component has been added assuming an evolving $q_{IR}$ from Delhaize et al. (2017). The green curves indicate the flux removed from the continuum SEDs by the effects of the CMB. The most impacted fluxes are the ones at 3 mm.</figcaption></figure>

our targets. In fact, large Herschel surveys have been used to constrain the slopes of local galaxies, finding $\beta \approx 1.5$ and up to a maximum of 2 (Boselli et al. 2012; Bianchi 2013). None of the local galaxies with an apparent best fitting $\beta > 2$ are at more than $1\sigma$ from $\beta = 2$ (Rémy-Ruyer et al. 2013). Also, although our high-redshift sample might be somewhat unusual because of particularly cold dust for its redshift (as discussed in more detail later in this paper), there is no evidence for steeper $\beta$ slope in cold galaxies in literature. For example, our MW ($T_{\text{dust}} = 19$ K) is colder than this sample but has no steep $\beta$ slope, where the analysis of the full MW dust SED yields $\beta \approx 1.8 \pm 0.2$ (Planck Collaboration et al. 2011). Without accounting for the CMB, the best-constrained galaxies in our sample (not considering ID85000922 given its poor SED$^{16}$) have $\beta > 2$ at significance ranging from $2\sigma$ to $3.5\sigma$. Taken together, their average SED is steeper than $\beta = 2$ at about the $5\sigma$ level. Therefore, we conclude that the steep observed $\beta$ slopes are unambiguous evidence that the CMB is having a measurable effect on our target galaxies, which arises from a severe reduction of the dust continuum observed at 3 mm, or beyond about 500 µm in the rest frame.

### 4.3. Cold Galaxies in the Distant Universe

In general, the impact of CMB on galaxy SEDs can be expressed as a function of the contrast between dust and CMB temperatures, with the highest impact when these temperatures are close. We recall that what counts is the apparent temperature, as encoded in the SED shape and rest-frame peak rather than the intrinsic temperature (which could differ in the case of high optical depths). It is therefore relevant to evaluate how the $T_{\text{dust}}$ for our galaxies relates to the cosmic evolution of temperatures in galaxies and the CMB. In the left panel of Figure 4, we compare dust temperature $T_{\text{dust}}$ for our sample to various objects taken from the literature (Magdis et al. 2011; Riechers et al. 2013, 2014, 2017; Ciesla et al. 2014; Béthermin et al. 2015; Pavesi et al. 2018; Schreiber et al. 2018), and to the evolving CMB temperature, versus redshift. We also fitted the available photometry for the $z = 6.9$ SPT0311-58 galaxy (Strandet et al. 2017; Marrone et al. 2018) and obtained $T_{\text{dust}} = 41 \pm 2$ K using the same MBB models including the CMB effect. In the right panel of Figure 4, we present the comparison in terms of the intensity of the radiation field $\langle U \rangle$: all values are consistently and homogeneously measured using the DL07 (Draine & Li 2007) models, eliminating any systematics in the $T_{\text{dust}}-z$ diagram coming from different assumptions of $\beta$ and fitting techniques of the MBB models. In terms of both $T_{\text{dust}}$ and $\langle U \rangle$, this sample has significant cold dust content with respect to both $z < 4$ MS galaxies and $z > 5$ SMGs.

<footnote n="16">16 Note that the galaxy ID85000922 ($z = 4.44$) also has an abnormally steep $\beta_{\text{obs}} = 3.7$ but with very large uncertainty. When including the CMB, we still fit its SED with a fixed $\beta = 2.3$, which seems required for producing both a good fit to the data and a conservatively cold dust temperature $T_{\text{dust}} = 30 \pm 4$ K. However, it should be noted that its dusty SED is not well constrained given the weak signal in the Herschel bands, and we suspect that the SCUBA2 850 µm photometry is boosted by noise.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Optically thin dust temperature $T_{\text{dust}}$ and intensity of the radiation field $\langle U \rangle$ vs. redshift for galaxies in the literature and the DSFGs studied in this work. Left: data in this work are marked as red circles, where the filled circles show the MBB results and the open ones show results derived by $\langle U \rangle$ scaled to DL07 models. Dust temperatures of MS samples (gray and black filled squares) are converted from $\langle U \rangle$ in Béthermin et al. (2015) and Magdis et al. (2012, 2017) by $\langle U \rangle = (T_{\text{dust}}/18.9)^{6.04}$. The solid and dashed magenta lines represent the evolutionary trends expected for a broken and universal FMR, respectively. Right: analog to Figure 5 in Magdis et al. (2017), overlaid by this work (red circles) and GN20 (blue cross). The filled and open circles correspond to $\langle U \rangle$ from MBB SED fitting and $\langle U \rangle$ derived from dust temperature according to $\langle U \rangle = (T_{\text{dust}}/18.9)^{6.04}$, respectively.</figcaption></figure>

The $z = 5.85$ galaxy has the highest $T_{\text{dust}}$ and $\langle U \rangle$ in this sample, though its $T_{\text{dust}}$ is still lower than the extrapolation of the $T_{\text{dust}}-z$ correlation for MS galaxies, and its $\langle U \rangle$ is lower than the extrapolation from a broken fundamental metallicity relation (FMR) as in Béthermin et al. (2015). The $z = 5.05$ source shows a dust temperature similar to the $z = 4$ starburst GN20 (Daddi et al. 2009). The $z = 4.44$ and $z = 3.62$ ones have the coldest $T_{\text{dust}}$ in this sample, which are just moderately above the CMB temperature by factors of 1.5–2.0. Not counting this study, among previously known galaxies, only GN20 is robustly colder than MS galaxies, while CRLE is lower than the MS at only ~1.2σ significance. The cold dust in SPT0311-58 at $z = 6.9$ had not been reported previously. A significant population of cold galaxies is now clearly detected.

Please note that the increase of $T_{\text{dust}}$ with redshift within our sample of four targets is simply a selection effect: we had originally selected our targets to be at $z_{\text{phot}} > 6$ based on photometric redshifts derived using an HFLS3 template. As a result, the lower the redshift with respect to $z = 6$, the intrinsically colder the temperature necessarily is (and the larger the implied dust masses). We show the $T_{\text{dust}}-z$ degeneracy trends for each galaxy in our sample in Figure 4 (and 7): any redshift trend that our four targets might define is spurious as linked to this degeneracy and affected by selection. If we had used a cold dust template like GN20, we would have obtained $z_{\text{phot}} = 4.0–5.5$ for our targets, which appear to be closer to the observations.

### 4.4. CMB Impacts Observables

We review and summarize here how CMB impacts various observables for this sample, as a necessary step before proceeding to infer their global properties for understanding their nature, and also as a reference for designing and conducting future observations of high-$z$ dusty galaxies.

**Continuum:** The CMB acts as a nonnegligible background for cold systems against which the line emission and

continuum emission are measured (Combes et al. 1999; Papadopoulos et al. 2000; da Cunha et al. 2013). According to Equation (18) in da Cunha et al. (2013), given the observed dust temperatures in our galaxies, the ratios of intrinsic/observed continuum are in the range of 1.32–1.45 at 3 mm and 1.06–1.09 at 1 mm. These ratios are monotonically increasing toward lower redshifts for our sample, an effect due to the fact that lower-redshift objects are colder and closer in relative terms to the CMB temperature at each redshift. The induced underestimate on the intrinsic dust continuum is more severe at 3 mm, and this has to be carefully accounted for when designing and conducting observations for high-redshift galaxies. The IR luminosities of the galaxies are also somewhat affected by the CMB, with a small underestimate by a factor of ≈1.05.

**CO emission:** A higher CMB temperature enhances the population of the high-$J$ CO levels and also corresponds to a higher background against which the lines must be detected (Combes et al. 1999; Papadopoulos et al. 2000; Obreschkow et al. 2009). Ultimately for fundamental thermodynamic reasons, $T_{\text{CMB}} \le T_{\text{excitation}} \le T_{\text{kinetic}}$ for any collisionally excited transition. Then, when the kinetic temperature of a cold gas reservoir gets closer to the CMB temperature, its low-$J$ CO line brightness diminishes by contrast to the background (Papadopoulos et al. 2000; da Cunha et al. 2013). Given the limited range of transitions observed for our galaxies and the degeneracy between temperature and density in the analysis of CO spectral line energy distributions (SLEDs; e.g., Weiß et al. 2005; Dannerbauer et al. 2009; Daddi et al. 2015), it is not possible to properly estimate the expected effect for our galaxies. Nevertheless, one could obtain some guidance by assuming, for example, that the kinetic temperature of the gas in the CO transitions we observed is close to the dust temperature in the same galaxies, as expected for local thermodynamic equilibrium (LTE). This would imply observed CO fluxes at the level of 70% of the intrinsic ones according to
<figure><img src="image_7.png" /><figcaption>**Figure 5.** $L'_{CO}-L_{IR}$ diagram with literature data (Carilli et al. 2013; Daddi et al. 2015; Liu et al. 2015a) and data in this work. Blue filled data points with error bars show observed values in this work, while the red ones show intrinsic values that are corrected for the CMB effect, assuming LTE conditions. Open circles with error bars are used for the $z = 3.623$ and $z = 5.847$ galaxies to emphasize that their $L'_{CO(5-4)}$ were extrapolated as $1.2P \times L'_{CO(4-3)}$ and $0.8 \times L'_{CO(6-5)}$, respectively. Note that a moderate line flux boosting (20%; see also Coogan et al. 2018, Figure 5) in this work has been accounted for in the error bar.</figcaption></figure>

Equation (32) in da Cunha et al. (2013). This seems to be reasonably consistent with observations. In Figure 5 we present the CO (5-4) luminosities of this sample in the $L'_{CO}-L_{IR}$ diagram. The average ratio of the observed CO (5-4) luminosity (blue data points) to $L_{IR}$ is 0.7 times lower than the best fit and 0.5 times lower than the linear fit in Liu et al. (2015a). Thus, the CMB effect provides a reasonable explanation for the seemingly low CO fluxes observed in our sample, although we cannot fully demonstrate this, lacking a measurement of the gas kinetic temperature. This is nevertheless interesting information for planning follow-up searches for CO emission in high-redshift galaxies, as reduction factors of up to ×2 in fluxes, due to the CMB, might not be uncommon to be encountered, particularly for galaxies with low CO excitation temperatures.

**Size measurement:** The CMB might also affect the 3 mm continuum sizes of our sources, in the case in which $T_{\text{dust}}$ gradients are present, similarly to local spiral galaxies in which the outskirts are colder. In that case, the outskirt surface brightness would be impacted more strongly by the CMB than the hotter centers, resulting in apparently smaller half-light radii. Taking as examples the temperature gradients found in some typical local galaxies M31, M33, and NGC 628, Zhang et al. (2016) computed that this effect could be at the level of 10%-20% at $z = 4-5$ for 3 mm bands, while a negligible size bias is expected at ~1 mm. Given that this effect is in any case not very strong (even for surface brightness it would be at most 40% for the local examples in Zhang et al. 2016), and that our sources are very compact and possibly merging driven starbursts where we consider it unlikely that they might have strong $T_{\text{dust}}$ gradients, we conclude that the CMB should not impact strongly the size measurements for this sample (the situation could be different for extended cold reservoirs in high-$z$ disks). We can verify this directly for the galaxy

ID85001674 at $z = 3.623$. This object is the coldest and closest in relative terms of the CMB temperature and shows a combined size $0.42 \pm 0.04$ from 3 mm+1 mm+0.8 mm observations. The sizes measured from 1 mm/0.8 mm data alone are $0.32 \pm 0.07$ and $0.38 \pm 0.04$, both about $1\sigma-1.5\sigma$ lower than the combined size. This is contrary to what would be expected if CMB were affecting sizes, showing that indeed any CMB-induced size bias should be quite small.

**Line widths and dynamical masses:** Together with sizes, velocity widths are used to estimate dynamical masses of galaxies, as we report in Table 3 ($M_{\text{dyn,2Re}}$), following the relations given in Daddi et al. (2010) and Coogan et al. (2018). The CMB could also impact the observed CO line width $V_{\text{FWHM}}$, in the presence of gas kinetic temperature gradients, as described above for the sizes. The magnitude of this effect has not been quantified so far, but we could expect that it should be roughly comparable to the one on sizes, i.e., probably small. The CO spectra of the galaxies in this work are still relatively low in S/N, and thus we assume that any CMB effect on line widths would likely be in the noise for our sources. In any case, the net effect of the CMB on derivation of dynamical mass from observables could imply that they are somewhat underestimated given $M_{\text{dyn,2Re}} \propto V_{\text{FWHM}}^2 \times R_e$.

Being much less affected by CMB and with higher intrinsic luminosities, [C II] 158 µm observations would be particularly valuable for this and other high-redshift samples. They would provide robust measurements on both size and dynamics, with minimal effects from the CMB. However, see discussion later in Section 5.2 for other reasons why [C II] 158 µm might also be affected in these galaxies.

## 5. Discussion

We have presented evidence for the existence of a population of fairly cold DSFGs in the distant universe, with SED shapes seemingly much colder (i.e., peaking at longer rest-frame wavelengths) than what is expected for average MS galaxies at the same redshift. We also showed that such cold galaxies are substantially affected by the CMB, so that observed millimeter fluxes are suppressed by factors up to 1.5, and CO fluxes might be reduced by factors up to 2, or perhaps more for lower-$J$ transitions. This implies that there are observational biases against such cold galaxies, and their real number density might be higher than what has been recovered so far. They might be the majority among SB galaxies in the distant universe, if the trend determined by Béthermin et al. (2015) for SBs from $z = 0-2$ continues to higher redshifts. This raises obvious questions: Why are these galaxies so cold? And, particularly, why can SB galaxies become systematically colder than MS galaxies, when at $z = 0$ the opposite situation is observed? In order to address these questions, we took a closer look at the overall physical properties that can be derived for these galaxies, given the observations available. In physical terms, the $\langle U \rangle$ values reflect the ratio of SFE to metallicity. Hence, these questions can be rephrased in equivalent terms: are these galaxies low-efficiency star formers, even lower efficiency than typical MS galaxies? Or particularly metal-rich? Should we consider them to be typical disk-like star formers at their redshifts, or might they be merger-affected typical starbursts?
<figure><img src="image_8.png" /><figcaption>**Figure 6.** IR luminosity surface density vs. dust temperature for DSFGs. Blue squares show optically thin data of SMGs in Simpson et al. (2017); red diamonds show the ALESS sample in Hodge et al. (2019). The Stefan-Boltzmann law is shown as a dashed line, which is only valid for blackbody dust clouds that are optically thick at all FIR wavelengths.</figcaption></figure>

## 5.1. Physical Properties

**SFR surface density:** We measure the SFR surface densities by converting the $L_{IR}$ estimates into SFRs (Kennicutt 1998) and dividing by the area $\Sigma_{SFR} = SFR/(2\pi R_e^2)$. We compared the SFR surface density to data from the literature, particularly taking samples shown in an analogous figure from Simpson et al. (2017). In Figure 6, our sample shows a trend of increasing luminosity surface density with dust temperature, consistent with trends in the lower-$z$ sample of Simpson et al. (2017) and Hodge et al. (2019). Meanwhile, our sample shows a higher-luminosity surface density than any lower-$z$ sample. Our values are at the limit of the Stefan-Boltzmann law, consistent with a single emitting, homogeneous starburst core. Our sample shows comparable SFR surface density with respect to local (U)LIRGs (Liu et al. 2015b), which are 2–3 orders of magnitude higher $\Sigma_{SFR}$ than normal spirals. The coldest galaxy ($z=3.623$) in this sample appears to be even colder than local IR-detected spirals from Liu et al. (2015b), while it has higher $\Sigma_{SFR}$ by 2–3 orders of magnitude. The high SFR surface densities suggest that these galaxies are more efficient starbursts than low-efficiency disks. Their extreme compactnesses indeed suggest mergers. In fact, using the evolutionary trend found by van der Wel et al. (2014), we would expect a disk-like galaxy at $z=5$ with $M_* \sim 5 \times 10^{10} M_\odot$ to have an FWHM size of $\sim 4.5$ kpc, substantially larger than the sizes we have measured for our objects (Table 4).

**Dust and stellar masses:** MBB dust masses are listed in Table 4. We calculated dust mass from DL07 templates, where we use the intrinsic SED (i.e., corrected for the CMB effect) of each galaxy as recovered from the MBB, take six data points (from 80 to 3000 µm in rest frame), and fit those data points with DL07. The MBB dust masses are consistent within 20% with the mass scaled to DL07; both result in a massive dust content for this sample. To constrain their stellar masses, we collected IRAC photometry in literature (Sanders et al. 2007) for the $z=5.05$ source and obtained SPLASH photometry for the others by PSF fitting the SPLASH images on ALMA positions. We fit stellar SEDs using a 200 Myr template

allowing for dust attenuation, following Liu et al. (2018). We list stellar masses in Table 4. Despite the IRAC detections well in the optical rest frame, we expect uncertainties at least at the level of $\times 2$ for our galaxies, given the lack of detections of age-dependent features in their optical SEDs. In some cases the dust-to-stellar-mass ratios can reach 10% levels (Figure 7), much higher than what is found for typical MS galaxies and consistent with merger-driven systems like bright SMGs (see also Tan et al. 2014; Rujopakarn et al. 2019).

**Starburstiness:** Defined by the ratio of the specific SFR for a given galaxy to the same specific SFR of an average MS galaxy at that redshift. At redshifts $4 < z < 6$ the MS value is not very securely determined, so this classification is uncertain. However, taking the extrapolation of the trends defined in Schreiber et al. (2015), we found that the $z=4.44$ galaxy has a starburstiness SFR/SFMSR $> 4-60$ and appears thus to be a secure starburst galaxy despite its large uncertainty on the stellar mass. The remaining galaxies have fairly high masses and starburstiness of 1.5–2.7, which are all above the MS extrapolation at $z>3.6-6$ but not by large factors. In fact, they could be considered MS galaxies in the definition by Rodighiero et al. 2011. However, as shown by several works (e.g., Elbaz et al. 2018; Puglisi et al. 2019), a substantial number of compact, probably merger-driven galaxies are observed even within the MS boundary at high masses, already at intermediate redshift $1 < z < 3$.

**Gas masses:** Inferring gas masses would be the most direct way to establish whether these objects are high- or low-efficiency star formers. We constrain the total gas mass (molecular plus any atomic gas) of our sample using three different methods: (1) based on the dynamical mass (three of which are only upper limits), after subtracting stellar mass using the Equation (3) in Daddi et al. (2010), assuming 20% of dark matter; (2) based on dust mass scaled to DL07, using FMR, FMR MS, elliptical, and solar-metallicity-based G/D conversion factors, respectively; (3) for the $z=3.62$ galaxy with a [C I] (1–0) detection, given the tight correlations between [C I] and low-$J$ CO in local and high-$z$ galaxies (Jiao et al. 2017, 2019; Valentino et al. 2018), we estimated its gas mass using the relation in Valentino et al. (2018) assuming an excitation temperature of $T_{ex}=30$ K and show the resulting $M_{gas,[C I]}$ in Table 4. We cannot estimate gas masses from CO, as the lowest-$J$ transition we detect is CO (4–3), whose conversion to CO (1–0) is plagued by too large and uncertain excitation corrections. In Figure 8, we compared the three different derivations of gas masses in the Schmidt-Kennicutt plane. Comparing to the relations for SBs and MS galaxies taken from Sargent et al. (2014), one can see that most values tend to be toward the high SFEs typical of SBs, with quite some large range and large uncertainties as exemplified by the scatter in the derived gas masses from different methods and assumptions. However, it is interesting to consider the $z=3.62$ source, which is the coldest in absolute terms and closest to the CMB temperature and has also alternative derivations of gas mass from the [C I] luminosity and the dynamical mass. The two latter estimates agree well between them but are lower than the estimate based on $M_{dust}$ by factors of 4–5, suggesting that $M_{dust}$ are probably overestimated.

## 5.2. Optically Thick Dust Emission in the FIR?

The behavior displayed by the galaxy ID85001674, as discussed in the previous subsection, suggests an alternative
<table><caption>**Table 4**

Physical Quantities</caption><thead><tr><th>ID</th><th>FWHM Size<br>(arcsec/kpc)</th><th>M<sub>*</sub><br>(10<sup>10</sup> M<sub>☉</sub>)</th><th>M<sub>dyn, 2 Re</sub><br>(10<sup>10</sup> M<sub>☉</sub>)</th><th>Σ<sub>SFR</sub><br>(M<sub>☉</sub> yr<sup>-1</sup> kpc<sup>2</sup>)</th><th>M<sub>gas,dyn</sub><br>(10<sup>10</sup> M<sub>☉</sub>)</th><th>M<sub>gas, C I</sub><br>(10<sup>10</sup> M<sub>☉</sub>)</th><th>M<sub>gas,FMR</sub><br>(10<sup>10</sup> M<sub>☉</sub>)</th></tr></thead><tbody><tr><td>85001929</td><td>&lt;0.33 (1.90)</td><td>11.6 ± 5.8</td><td>&lt;15.9</td><td>&gt;196</td><td>&lt;8.2</td><td>...</td><td>5.7 ± 1.3</td></tr><tr><td>20010161</td><td>&lt;0.25 (1.55)</td><td>9.0 ± 5.2</td><td>&lt;12.0</td><td>&gt;135</td><td>&lt;6.8</td><td>...</td><td>12.5 ± 2.8</td></tr><tr><td>85000922</td><td>&lt;0.43 (2.83)</td><td>1.8 ± 1.6</td><td>&lt;19.4</td><td>&gt;31</td><td>&lt;16.9</td><td>...</td><td>25.1 ± 8.6</td></tr><tr><td>85001674</td><td>0.42 ± 0.04 (3.10 ± 0.30)</td><td>3.9 ± 2.0</td><td>17.6 ± 2.6</td><td>22<sup>+13</sup><sub>-9</sub></td><td>10.2 ± 4.2</td><td>9.7 ± 2.2</td><td>48.3 ± 12.1</td></tr></tbody></table>

Note. $M_*$: stellar mass; $M_{dyn, 2 Re}$: dynamic mass within FWHM size; SFR surface density $\Sigma_{SFR} \equiv SFR/(2\pi R_e^2)$, where $R_e$ is half of the FWHM size.

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Dust-to-stellar-mass ratio vs. redshift. Data are taken from Tan et al. (2014) and overlaid by data in this work.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Gas mass vs. SFR for this work. Gas mass is calculated by multiple methods as labeled (see text for details).</figcaption></figure>

explanation for the nature of these cold galaxies, which could solve a number of problems with them, particularly the coldest ones. In fact, the dust-to-stellar-mass ratios of the coldest galaxies are very large, reaching 10%, which seems very high (Figures 7 and 8). Given that the mass in dust can be taken as a proxy of the mass in metals, this would imply that the metallicity in these objects could be truly extreme. Similarly, also the gas-to-stellar-mass fraction would be seemingly ×10

or higher, as is only expected for primordial galaxies, which would be in contradiction with the high metallicities (Figure 8). All these values are derived from $M_{dust}$ estimates, and in all cases assuming models with optically thin dust emission. These measurements are not supported by a dynamical estimate of the gas masses, nor by derivation from [C I] (1–0), a transition known to be in most cases optically thin given the low fractional abundance of carbon in the neutral form (Valentino et al. 2018). It could well be that the dust continuum emission is instead optically thick in these systems, even at the peak of the emission and possibly somewhat beyond, similar to the optically thick SB nuclei found in Arp 220 (Papadopoulos et al. 2010; Wilson et al. 2014; Scoville et al. 2017). This is supported by the extreme compactness of our sources and their very high SFR surface densities. A simple calculation for ID85001674 assuming a spherical geometry with constant density of dust $\rho$ within the measured radius $R_e$, using the measured dust mass from the thin models, and using $\kappa$ from the Jones et al. (2013) opacity models would imply an optical depth ($\tau = \kappa \rho R_e$) at 100 µm rest (the observed SED peak) toward the center of about 1. Of course, there are many implicit assumptions and unknowns in this calculation,<fnref n="17" /> but certainly the case for optically thick dust until the observed SED peak is a plausible one for these galaxies (and for the others, which are likely even more compact).

If the dust emission is optically thick around 100 µm (and possibly beyond), the suppression of the Wien emission would make the observed SEDs seem cold, while the intrinsic temperatures would be much higher (e.g., Scoville & Kwan 1976; Condon et al. 1991; Conley et al. 2011; Scoville 2012). Conversely, together with the underestimate in $T_{dust}$, we would be overestimating $M_{dust}$ by some factor that would depend ultimately on the optical depth (i.e., from the ratio between the intrinsic $T_{dust}$ and the real $T_{dust}$). In all cases $L_{IR}$ (hence SFR) derivations would be unaffected.

The standard treatment of optically thick emission assumes a constant $T_{dust}$ at all wavelengths (e.g., Casey et al. 2012) and is unlikely to be realistic, given that radial gradients of $T_{dust}$ are expected (e.g., Scoville 2012). Explicitly accounting for such effects is difficult, and we hope to present calculations of more physical thick models in future works. Nevertheless, we have fitted our targets using the more standard approach of optically thick sources with constant temperature and accounting for the effects of the CMB. Results are reported in Table 3, along with the other estimates. The empirical evidence discussed above for ID85001674 requires that we might be overestimating dust masses by factors of up to 4–5 for our sample. These thick

<footnote n="17">¹⁷ See also the similar calculation for SMG attenuation in SMG by Simpson et al. (2017), suggesting hundreds of magnitudes of attenuation of the optical rest-frame light.</footnote>
models indeed already account for a ×3 reduction of $M_{\text{dust}}$ for ID85001674 and ≈×2 for the other galaxies. True corrections might be even larger.

If these sources are really optically thick, then the ultimate reason why CMB is so strongly affecting these galaxies is not because their physical dust temperatures are cold but because they seem cold. Similarly, this would also explain why these galaxies seem to be colder than MS galaxies, when they might not be in reality. $T_{\text{dust}}$ values in Table 3 derived for the thick fits are already comparable to or larger than MS galaxies at the same redshifts. In this case their CO kinetic temperatures would also be higher, rendering CMB effects on their line fluxes negligible: the reduction of CO fluxes might be instead due to high optical depths (see also Papadopoulos et al. 2010).

A strong reduction of $M_{\text{dust}}$ in SBs would, finally, also solve the problematic observations that their dust-to-stellar-mass ratios are ~5× larger than those of MS galaxies (Béthermin et al. 2015), while dynamical estimates suggest that the gas-to-stellar-mass ratios are near unity (Silverman et al. 2015, 2018a, 2018b). The case should be thus seriously considered that many of the merging driven SBs are actually optically thick in the FIR, a case previously proposed for local ULIRGs by Papadopoulos et al. (2010).

More observations would be needed to confirm the optically thick case. One interesting path would be comparing to [C I]-derived excitation temperatures, given that [C I] lines should remain thin. Another possibility would be to use [C II] 158 µm as a direct tracer of the gas mass (Zanella et al. 2018) and from there more realistic $M_{\text{dust}}$ estimates via assumptions of G/D ratios: while [C II] could be in principle heavily affected by optical depth effects, it seems to provide fair gas mass estimates in local ULIRGs (perhaps by lucky coincidence). More observations are also needed to establish whether these cold (i.e., with long peak rest wavelengths) SEDs are prevalent in the distant universe, and if the optically thick case could provide a general explanation for the odd inversion that seems to take place at high redshifts, with SB galaxies becoming colder than MS galaxies.

## 6. Summary and Conclusions

We have used ALMA to obtain spectroscopic redshifts and investigate the properties of four galaxies detected in the FIR in the *super-deblended* catalog by Jin et al. (2018) and selected to be at $z > 6$ based on photometric redshift derived using a well-studied $z = 6.3$ galaxy, HFLS3. Our findings can be summarized as follows:

1. We securely confirm the redshifts to be at $3.6 < z < 5.85$ based on multiple CO/C I transitions, with the exception of one object that has only one secure line detected, for which we derive a less secure $z = 5.05$.

2. Our sample contains the most distant spectroscopically confirmed IR-detected galaxy in COSMOS, ID85001929 at $z = 5.85$.

3. When not accounting for CMB effects, the galaxies display highly unusual steep slopes in the RJ regimes. The slopes are back to normal when accounting for the CMB. We think that this could be the first direct evidence of the impact of the CMB on galaxy observables at high redshifts.

4. Our galaxies are anomalously cold (i.e., peak at long rest wavelengths), colder than what typical MS galaxies are expected to be at these redshifts.

5. While in general cold dust emission is observed in galaxies with low star formation efficiency, which is typically found in quiescent star-forming disks, our galaxies appear instead to be efficient star formers with compact sizes and high SFR surface densities, more typical of galaxy mergers.

6. We investigate possible reasons for the coldness of the observed FIR SEDs. An interesting possibility is that the galaxies have optically thick dust emission up until the peak of their emission at about 100 µm, and possibly beyond. We propose a scheme of high dust optical depths “hiding” a warm compact dust mass associated with star-forming, merger-driven “cores,” rather than fitting the associated IR SED with large cold dust reservoirs.

7. A large fraction of high (and even low) redshift SBs being optically thick in the FIR would provide an appealing solution to a number of odd observations for this population, including claims of extreme gas fractions and their seemingly redshift-independent $T_{\text{dust}}$ evolution. Future work is required to demonstrate whether this is actually the case.

We are grateful to Simone Bianchi and Zhi-Yu Zhang for discussions and to the referee for constructive criticism. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2017.1.00373.S, #2016.1.00463.S, and #2016.1.00279.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. S.J. and Q.G. acknowledge support from the National Key Research and Development Program of China (No. 2017YFA0402703) and the National Natural Science Foundation of China (Key Project No. 11733002). S.J. acknowledges financial support from the Spanish Ministry of Science, Innovation and Universities (MICIU) under grant AYA2017-84061-P, co-financed by FEDER (European Regional Development Funds). G.E.M. acknowledges support from the Villum Fonden research grant 13160 “Gas to stars, stars to dust: tracing star formation across cosmic time,” the Cosmic Dawn Center of Excellence funded by the Danish National Research Foundation, and the ERC Consolidator Grant funding scheme (project ConTExt, grant No. 648179) D.L. and E.S. acknowledge funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement No. 694343). Y.G. is partially supported by National Key Basic Research and Development Program of China (grant No. 2017YFA0402704), NSFC grant Nos. 11861131007 and 11420101002, and Chinese Academy of Sciences Key Research Program of Frontier Sciences (grant No. QYZDJ-SSW-SLH008). This paper makes use of the GILDAS software developed by IRAM and available at http://www.iram.fr/IRAMFR/GILDAS.

## ORCID iDs

S. Jin https://orcid.org/0000-0002-8412-7951

E. Daddi https://orcid.org/0000-0002-3331-9590

G. E. Magdis https://orcid.org/0000-0002-4872-2294
D. Liu https://orcid.org/0000-0001-9773-7479

E. Schinnerer https://orcid.org/0000-0002-3933-7677

Q. Gu https://orcid.org/0000-0002-3890-3729

Y. Gao https://orcid.org/0000-0003-0007-2197

A. Calabrò https://orcid.org/0000-0003-2536-1614

References

Aretxaga, I., Wilson, G. W., Aguilar, E., et al. 2011, MNRAS, 415, 3831

Bertoldi, F., Carilli, C., Aravena, M., et al. 2007, ApJS, 172, 132

Béthermin, M., Daddi, E., Magdis, G., et al. 2015, A&A, 573, A113

Béthermin, M., Dole, H., Cousin, M., & Bavouzet, N. 2010, A&A, 516, A43

Béthermin, M., Le Floc'h, E., Ilbert, O., et al. 2012, A&A, 542, A58

Bianchi, S. 2013, A&A, 552, A89

Boselli, A., Ciesla, L., Cortese, L., et al. 2012, A&A, 540, A54

Calabrò, D., Daddi, E., Puglisi, A., et al. 2019, A&A, 623, A64

Capak, P. L., Riechers, D., Scoville, N. Z., et al. 2011, Natur, 470, 233

Carilli, C. L., Riechers, D., Walter, F., et al. 2013, ApJ, 763, 120

Casey, C. M., Berta, S., Béthermin, M., et al. 2012, ApJ, 761, 140

Chapman, S. C., Blain, A. W., Smail, I., & Ivison, R. J. 2005, ApJ, 622, 772

Ciesla, L., Boquien, M., Boselli, A., et al. 2014, A&A, 565, A128

Combes, F., Maoli, R., & Omont, A. 1999, A&A, 345, 369

Condon, J. J., Huang, Z.-P., Yin, Q. F., & Thuan, T. X. 1991, ApJ, 378, 65

Conley, A., Cooray, A., Vieira, J. D., et al. 2011, ApJL, 732, L35

Coogan, R. T., Daddi, E., Sargent, M. T., et al. 2018, MNRAS, 479, 703

Cowie, L. L., Barger, A. J., Hsu, L.-Y., et al. 2017, ApJ, 837, 139

da Cunha, E., Groves, B., Walter, F., et al. 2013, ApJ, 766, 13

Daddi, E., Bournaud, F., Walter, F., et al. 2010, ApJ, 713, 686

Daddi, E., Dannerbauer, H., Liu, D., et al. 2015, A&A, 577, A46

Daddi, E., Dannerbauer, H., Stern, D., et al. 2009, ApJ, 694, 1517

Daddi, E., Jin, S., Strazzullo, V., et al. 2017, ApJL, 846, L31

Dannerbauer, H., Daddi, E., Riechers, D. A., et al. 2009, ApJL, 698, L178

Delhaize, J., Smolčić, V., Delvecchio, I., et al. 2017, A&A, 602, A4

Draine, B. T., & Li, A. 2007, ApJ, 657, 810

Elbaz, D., Dickinson, M., Hwang, H. S., et al. 2011, A&A, 533, A119

Elbaz, D., Leiton, R., Nagar, N., et al. 2018, A&A, 616, A110

Fudamoto, Y., Ivison, R. J., Oteo, I., et al. 2017, MNRAS, 472, 2028

Geach, J. E., Dunlop, J. S., Halpern, M., et al. 2017, MNRAS, 465, 1789

González-López, J., Decarli, R., Pavesi, R., et al. 2019, ApJ, 882, 139

Hodge, J. A., Smail, I., Walter, F., et al. 2019, ApJ, 876, 130

Hurley, P. D., Oliver, S., Betancourt, M., et al. 2017, MNRAS, 464, 885

Ikarashi, S., Ivison, R. J., Caputi, K. I., et al. 2015, ApJ, 810, 133

Jiao, Q., Zhao, Y., Lu, N., et al. 2019, ApJ, 880, 133

Jiao, Q., Zhao, Y., Zhu, M., et al. 2017, ApJL, 840, L18

Jin, S., Daddi, E., Liu, D., et al. 2018, ApJ, 864, 56

Jones, A. P., Fanciullo, L., Köhler, M., et al. 2013, A&A, 558, A62

Juneau, S., Narayanan, D. T., Moustakas, J., et al. 2009, ApJ, 707, 1217

Kennicutt, R. C., Jr. 1998, ApJ, 498, 541

Laigle, C., McCracken, H. J., Ilbert, O., et al. 2016, ApJS, 224, 24

Lee, N., Sanders, D. B., Casey, C. M., et al. 2013, ApJ, 778, 131

Liang, L., Feldmann, R., Kereš, D., et al. 2019, MNRAS, 489, 1397

Lilly, S. J., Le Fevre, O., Hammer, F., & Crampton, D. 1996, ApJL, 460, L1

Liu, D., Daddi, E., Dickinson, M., et al. 2018, ApJ, 853, 172

Liu, D., Gao, Y., Isaak, K., et al. 2015a, ApJL, 810, L14

Liu, L., Gao, Y., & Greve, T. R. 2015b, ApJ, 805, 31

Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415

Madau, P., Pozzetti, L., & Dickinson, M. 1998, ApJ, 498, 106

Magdis, G. E., Daddi, E., Béthermin, M., et al. 2012, ApJ, 760, 6

Magdis, G. E., Elbaz, D., Dickinson, M., et al. 2011, A&A, 534, A15

Magdis, G. E., Rigopoulou, D., Daddi, E., et al. 2017, A&A, 603, A93

Magnelli, B., Lutz, D., Saintonge, A., et al. 2014, A&A, 561, A86

Marrone, D. P., Spilker, J. S., Hayward, C. C., et al. 2018, Natur, 553, 51

McCracken, H. J., Milvang-Jensen, B., Dunlop, J., et al. 2012, A&A, 544, A156

McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in ASP Conf. Ser. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell (San Francisco, CA: ASP), 127

Miettinen, O., Novak, M., Smolčić, V., et al. 2017, A&A, 602, A54

Murphy, E. J., Momjian, E., Condon, J. J., et al. 2017, ApJ, 839, 35

Obreschkow, D., Heywood, I., Klöckner, H.-R., & Rawlings, S. 2009, ApJ, 702, 1321

Omont, A., Yang, C., Cox, P., et al. 2013, A&A, 551, A115

Papadopoulos, P. P., Isaak, K., & van der Werf, P. 2010, ApJ, 711, 757

Papadopoulos, P. P., Röttgering, H. J. A., van der Werf, P. P., et al. 2000, ApJ, 528, 626

Pavesi, R., Riechers, D. A., Sharon, C. E., et al. 2018, ApJ, 861, 43

Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2011, A&A, 536, A16

Puglisi, A., Daddi, E., Liu, D., et al. 2019, ApJL, 877, L23

Puglisi, A., Daddi, E., Renzini, A., et al. 2017, ApJL, 838, L18

Rémy-Ruyer, A., Madden, S. C., Galliano, F., et al. 2013, A&A, 557, A95

Renaud, F., Athanassoula, E., Amram, P., et al. 2018, MNRAS, 473, 585

Renaud, F., Bournaud, F., Agertz, O., et al. 2019, A&A, 625, 65

Riechers, D. A., Bradford, C. M., Clements, D. L., et al. 2013, Natur, 496, 329

Riechers, D. A., Carilli, C. L., Capak, P. L., et al. 2014, ApJ, 796, 84

Riechers, D. A., Leung, T. K. D., Ivison, R. J., et al. 2017, ApJ, 850, 1

Rodighiero, G., Daddi, E., Baronchelli, I., et al. 2011, ApJL, 739, L40

Roseboom, I. G., Oliver, S. J., Kunz, M., et al. 2010, MNRAS, 409, 48

Rujopakarn, W., Daddi, E., Rieke, G. H., et al. 2019, ApJ, 882, 107

Safarzadeh, M., Ferguson, H. C., Lu, Y., Inami, H., & Somerville, R. S. 2015, ApJ, 798, 91

Sanders, D. B., Salvato, M., Aussel, H., et al. 2007, ApJS, 172, 86

Sargent, M. T., Daddi, E., Béthermin, M., et al. 2014, ApJ, 793, 19

Schreiber, C., Elbaz, D., Pannella, M., et al. 2018, A&A, 609, A30

Schreiber, C., Pannella, M., Elbaz, D., et al. 2015, A&A, 575, A74

Scoville, N., Murchikova, L., Walter, F., et al. 2017, ApJ, 836, 66

Scoville, N. Z. 2012, in Evolution of Star Formation and Gas, ed. J. Falcón-Barroso & J. H. Knapen (Cambridge: Cambridge Univ. Press), 491

Scoville, N. Z., & Kwan, J. 1976, ApJ, 206, 718

Silk, J., & Spaans, M. 1997, ApJL, 488, L79

Silverman, J. D., Daddi, E., Rodighiero, G., et al. 2015, ApJL, 812, L23

Silverman, J. D., Daddi, E., Rujopakarn, W., et al. 2018a, ApJ, 868, 75

Silverman, J. D., Rujopakarn, W., Daddi, E., et al. 2018b, ApJ, 867, 92

Simpson, J. M., Smail, I., Swinbank, A. M., et al. 2015, ApJ, 799, 81

Simpson, J. M., Smail, I., Swinbank, A. M., et al. 2015, ApJ, 839, 58

Smolčić, V., Karim, A., Miettinen, O., et al. 2015, A&A, 576, A127

Smolčić, V., Novak, M., Bondi, M., et al. 2017, A&A, 602, A1

Soifer, B. T., Neugebauer, G., Matthews, K., et al. 2000, AJ, 119, 509

Strandet, M. L., Weiss, A., De Breuck, C., et al. 2017, ApJL, 842, L15

Tan, Q., Daddi, E., Magdis, G., et al. 2014, A&A, 569, A98

Valentino, F., Magdis, G. E., Daddi, E., et al. 2018, ApJ, 869, 27

van der Wel, A., Franx, M., van Dokkum, P. G., et al. 2014, ApJ, 788, 28

Vieira, J. D., Marrone, D. P., Chapman, S. C., et al. 2013, Natur, 495, 344

Walter, F., Decarli, R., Carilli, C., et al. 2012, Natur, 486, 233

Walter, F., Weiß, A., Downes, D., Decarli, R., & Henkel, C. 2011, ApJ, 730, 18

Weiß, A., Walter, F., & Scoville, N. Z. 2005, A&A, 438, 533

Wilson, C. D., Rangwala, N., Glenn, J., et al. 2014, ApJL, 789, L36

Yang, C., Gavazzi, R., Beelen, A., et al. 2019, A&A, 624, 138

Yang, C., Omont, A., Beelen, A., et al. 2016, A&A, 595, A80

Zanella, A., Daddi, E., Magdis, G., et al. 2018, MNRAS, 481, 1976

Zavala, J. A., Montaña, A., Hughes, D. H., et al. 2018, NatAs, 2, 56

Zhang, Z.-Y., Papadopoulos, P. P., Ivison, R. J., et al. 2016, RSOS, 3, 160025