# The hidden side of cosmic star formation at z > 3

M.-Y. Xiao, D. Elbaz, C. Gómez-Guijarro, L. Leroy, L.-J. Bing, E. Daddi, B.

Magnelli, M. Franco, L. Zhou, M. Dickinson, et al.

► To cite this version:

M.-Y. Xiao, D. Elbaz, C. Gómez-Guijarro, L. Leroy, L.-J. Bing, et al. The hidden side of cosmic star formation at z > 3. Astronomy & Astrophysics - A&A, 2023, 672, pp.A18. 10.1051/0004-6361/202245100. hal-04289564

HAL Id: hal-04289564

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

Submitted on 16 Nov 2023

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>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# The hidden side of cosmic star formation at $z > 3$

## Bridging optically dark and Lyman-break galaxies with GOODS-ALMA

M.-Y. Xiao<sup>1,2,3</sup>, D. Elbaz<sup>2</sup>, C. Gómez-Guijarro<sup>2</sup>, L. Leroy<sup>2</sup>, L.-J. Bing<sup>4</sup>, E. Daddi<sup>2</sup>, B. Magnelli<sup>2</sup>, M. Franco<sup>5</sup>, L. Zhou<sup>1,3</sup>, M. Dickinson<sup>6</sup>, T. Wang<sup>1,3</sup>, W. Rujopakarn<sup>7,8,9</sup>, G. E. Magdis<sup>10,11,12</sup>, E. Treister<sup>13</sup>, H. Inami<sup>14</sup>, R. Demarco<sup>15</sup>, M. T. Sargent<sup>16,17</sup>, X. Shu<sup>18</sup>, J. S. Kartaltepe<sup>19</sup>, D. M. Alexander<sup>20</sup>, M. Béthermin<sup>4</sup>, F. Bournaud<sup>2</sup>, L. Ciesla<sup>4</sup>, H. C. Ferguson<sup>21</sup>, S. L. Finkelstein<sup>5</sup>, M. Giavalisco<sup>22</sup>, Q.-S. Gu<sup>1,3</sup>, D. Iono<sup>23,24</sup>, S. Juneau<sup>6</sup>, G. Lagache<sup>4</sup>, R. Leiton<sup>15</sup>, H. Messias<sup>25,26</sup>, K. Motohara<sup>27</sup>, J. Mullaney<sup>28</sup>, N. Nagar<sup>15,29</sup>, M. Pannella<sup>30,31</sup>, C. Papovich<sup>32,33</sup>, A. Pope<sup>25</sup>, C. Schreiber<sup>34</sup>, and J. Silverman<sup>9</sup>

*(Affiliations can be found after the references)*

Received 30 September 2022 / Accepted 1 February 2023

## ABSTRACT

Our current understanding of the cosmic star formation history at $z > 3$ is primarily based on UV-selected galaxies (Lyman-break galaxies, i.e., LBGs). Recent studies of H-dropouts (HST-dark galaxies) have revealed that we may be missing a large proportion of star formation that is taking place in massive galaxies at $z > 3$. In this work, we extend the H-dropout criterion to lower masses to select optically dark or faint galaxies (OFGs) at high redshifts in order to complete the census between LBGs and H-dropouts. Our criterion ($H > 26.5$ mag & [4.5] < 25 mag) combined with a de-blending technique is designed to select not only extremely dust-obscured massive galaxies but also normal star-forming galaxies (typically $E(B-V) > 0.4$) with lower stellar masses at high redshifts. In addition, with this criterion, our sample is not contaminated by massive passive or old galaxies. In total, we identified 27 OFGs at $z_{phot} > 3$ (with a median of $z_{med} = 4.1$) in the GOODS-ALMA field, covering a wide distribution of stellar masses with $log(M_{\star}/M_{\odot}) = 9.4–11.1$ (with a median of $log(M_{\star med}/M_{\odot}) = 10.3$). We find that up to 75% of the OFGs with $log(M_{\star}/M_{\odot}) = 9.5–10.5$ were neglected by previous LBGs and H-dropout selection techniques. After performing an optical-to-millimeter stacking analysis of the OFGs, we find that rather than being limited to a rare population of extreme starbursts, these OFGs represent a normal population of dusty star-forming galaxies at $z > 3$. The OFGs exhibit shorter gas depletion timescales, slightly lower gas fractions, and lower dust temperatures than the scaling relation of typical star-forming galaxies. Additionally, the total star formation rate ($SFR_{tot} = SFR_{IR} + SFR_{UV}$) of the stacked OFGs is much higher than the $SFR_{UV}^{corr}$ ($SFR_{UV}$ corrected for dust extinction), with an average $SFR_{tot}/SFR_{UV}^{corr} = 8 \pm 1$, which lies above (~0.3 dex) the 16–84th percentile range of typical star-forming galaxies at $3 \le z \le 6$. All of the above suggests the presence of hidden dust regions in the OFGs that absorb all UV photons, which cannot be reproduced with dust extinction corrections. The effective radius of the average dust size measured by a circular Gaussian model fit in the uv plane is $R_{e(1.13 mm)} = 1.01 \pm 0.05$ kpc. After excluding the five LBGs in the OFG sample, we investigated their contributions to the cosmic star formation rate density (SFRD). We found that the SFRD at $z > 3$ contributed by massive OFGs ($log(M_{\star}/M_{\odot}) > 10.3$) is at least two orders of magnitude higher than the one contributed by equivalently massive LBGs. Finally, we calculated the combined contribution of OFGs and LBGs to the cosmic SFRD at $z = 4–5$ to be $4 \times 10^{-2} M_{\odot} yr^{-1}$ Mpc$^{-3}$, which is about 0.15 dex (43%) higher than the SFRD derived from UV-selected samples alone at the same redshift. This value could be even larger, as our calculations were performed in a very conservative way.

**Key words.** galaxies: high-redshift – galaxies: evolution – galaxies: star formation – galaxies: photometry – submillimeter: galaxies

## 1. Introduction

Our current knowledge of the first two billion years of cosmic star formation history is based mainly on (i) UV-selected galaxies, such as Lyman-break galaxies (LBGs; e.g., Giavalisco et al. 2004; Bouwens et al. 2012a, 2015, 2020; Oesch et al. 2014, 2015, 2018; Madau & Dickinson 2014), which are known to be biased against massive galaxies; and (ii) the most massive and extremely dusty starburst galaxies (e.g., Walter et al. 2012; Marrone et al. 2018), which are limited to a rare population and are not representative of the most common galaxies typically on the star-formation main sequence (SFMS; e.g., Elbaz et al. 2007, 2011; Noeske et al. 2007; Magdis et al. 2010; Whitaker et al. 2012, 2014; Speagle et al. 2014; Schreiber et al. 2015; Lee et al. 2015; Leslie et al. 2020). Recent Atacama Large Millimeter and Submillimeter Array (ALMA) and *Spitzer* observations have identified a more abundant and less extreme population of obscured galaxies at $z > 3$ (e.g., H-dropouts in Wang et al. 2019;

HST-dark galaxies in Zhou et al. 2020, optically dark or faint galaxies in Gómez-Guijarro et al. 2022a), revealing that a significant population of high-$z$ optically dark or faint galaxies have been missed, and they may dominate the massive end of the stellar mass function. The contribution of these optically dark or faint galaxies to the cosmic star formation rate density (SFRD) at $z > 3$ could be substantial, corresponding up to ~10–25% of the SFRD from LBGs, or even up to ~40%, depending on the methodology (e.g., Wang et al. 2019; Williams et al. 2019; Gruppioni et al. 2020; Fudamoto et al. 2021; Talia et al. 2021; Enia et al. 2022; Shu et al. 2022; Barrufet et al. 2022).

Optically dark or faint galaxies have generally been completely undetected or tentatively detected with very low significance even in the deepest HST/WFC3 images (typical $5\sigma$ depth of $H > 27$ mag), but brighter at longer wavelengths such as *Spitzer*/IRAC 3.6 and 4.5 µm, (e.g., Franco et al. 2018; Yamaguchi et al. 2019; Zhou et al. 2020; Smail et al. 2021; Gómez-Guijarro et al. 2022a). In GOODS-ALMA 1.0,
Franco et al. (2020a) reported six optically dark galaxies (i.e., HST-dark galaxies) out of 35 galaxies detected above $3.5\sigma$ at 1.13 mm. With the ALMA spectroscopic follow-up, Zhou et al. (2020) further analyzed these six optically dark galaxies in detail and found that four (~70%) could be associated with a $z \sim 3.5$ overdensity (corresponding to OFG1, 2, 25, and 27 in the southwest region of Fig. 1). Afterward, in the deeper GOODS-ALMA 2.0, Gómez-Guijarro et al. (2022a) updated the sample with 13 optically dark or faint galaxies (including six in the GOODS-ALMA 1.0), among a total of 88 sources detected above $3.5\sigma$ at 1.13 mm. So far, we do not have a unified and clear definition of optically dark or faint galaxies. The six optically dark galaxies in GOODS-ALMA 1.0 have no optical counterparts in the deepest H-band based on the CANDELS catalog down to $H = 28.16$ AB ($5\sigma$ limiting depth in CANDELS-deep field). However, two of them show H-band magnitudes of approximately 25 mag and 27 mag following a de-blending process (Zhou et al. 2020). The remaining seven sources were classified as optically dark or faint galaxies because they are currently undetected or very faint in the H-band of the deepest fields and other shorter wavelength bands (Gómez-Guijarro et al. 2022a).

Therefore, the purpose of our work is to first make a clear definition of the selection of optically dark or faint galaxies. Furthermore, by systematically studying optically dark or faint galaxies in the GOODS-ALMA field, we aim to obtain a more complete picture of the cosmic star formation history in the $z > 3$ Universe. In this work, our sample includes not only sources detected by ALMA 1.13 mm, but also those that are currently undetected (i.e., no millimeter counterparts in the GOODS-ALMA 2.0 catalog) to obtain a somewhat complete sample of optically dark or faint galaxies. By stacking their optical to millimeter emission, we can, however, investigate the differences between the optically dark or faint galaxies detected by ALMA 1.13 mm and those that remain undetected.

This paper is organized as follows. In Sect. 2, we describe the GOODS-ALMA survey and the multiwavelength data used. In Sect. 3, we present our selection criterion for optically dark or faint galaxies at $z > 3$. In Sect. 4, we study the properties of individual sources in our sample, such as the redshift, stellar mass, star formation rate (SFR), molecular gas mass, and dust temperature. In Sect. 5, we present and discuss the properties of optically dark or faint galaxies mainly based on our optical to millimeter stacking analysis. In Sect. 6, we calculate the cosmic SFRD contributed by optically dark or faint galaxies and discuss the level of the incompleteness of our understanding of dust-obscured star formation in the $z > 3$ Universe. Finally, we summarize our main conclusions in Sect. 7.

Throughout this paper, we adopt a Chabrier initial mass function (IMF; Chabrier 2003) to estimate SFR and stellar mass. We assume cosmological parameters of $H_0 = 70$ km s$^{-1}$ Mpc$^{-1}$, $\Omega_M = 0.3$, and $\Omega_\Lambda = 0.7$. When necessary, data from the literature have been converted with a conversion factor of $M_\star$ (Salpeter 1955, IMF) $= 1.7 \times M_\star$ (Chabrier 2003, IMF). All magnitudes are in the AB system (Oke & Gunn 1983), such that $m_{AB} = 23.9-2.5 \times \log(S_v [\mu Jy])$.

## 2. Data

### 2.1. GOODS-ALMA survey

GOODS-ALMA is an ALMA 1.13 mm survey in the deepest part of the Great Observatories Origins Deep Survey South field (GOODS-South; Dickinson et al. 2003; Giavalisco et al. 2004). It covers a continuous area of $72.42 \text{ arcmin}^2$ (effective

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Sky distributions of the optically dark or faint galaxies (OFGs) in the GOODS-ALMA 2.0 map at 1.13 mm. North is up and east is to the left.</figcaption></figure>

area of ~69 arcmin$^2$ if the shallower areas at the edges are trimmed off) with ALMA band 6 receivers, centered at $\alpha = 3^h 32^m 30.0^\circ$, $\delta = -27^\circ 48' 00''$ (J2000). The observations were obtained from Cycle 3 and Cycle 5, with two different array configurations to include both small and large spatial scales. The ALMA Cycle 3 observations (high-resolution dataset; Project ID: 2015.1.00543.S; PI: D. Elbaz) were conducted between August and September 2016 in the C40-5 array configuration with a total on-source exposure time of approximately 14.06 h, providing a high-resolution image with the longest baseline of 1808 m. The ALMA Cycle 5 observations (low-resolution dataset; Project ID: 2017.1.00755.S; PI: D. Elbaz) were performed between July 2018 and March 2019 with the C43-2 array configuration with a total on-source exposure time of 14.39 h, providing a low-resolution image with the longest baseline of 360.5 m.

The calibration was processed using the Common Astronomy Software Application package (CASA; McMullin et al. 2007) with the standard pipeline. We systematically inspected calibrated visibilities and added a few additional flags to the original calibration scripts. The calibrated visibilities were then time- and frequency-averaged over 120 s and 8 channels, respectively, to reduce the computational time for subsequent continuum imaging. Given the excellent coverage of the uv plane and the absence of very bright sources, we used the task TCLEAN in CASA version 5.6.1–8 to produce a dirty map with 0.05" pixels and a natural weighting scheme to avoid potential biases from the CLEAN algorithm. The resulting high- and low-resolution 1.13 mm continuum maps have similar root mean square (rms) sensitivities, that is, $\sigma \approx 89.0$ and $95.2 \mu Jy beam^{-1}$, with spatial resolutions of full width at half maximum (FWHM) $\approx 0.251 \times 0.232$ and $1.330 \times 0.935$, respectively. To improve the sensitivity, we concatenated these two data configurations in the uv plane with visibility weights proportional to $1:1$. The combined map achieves an rms sensitivity of $\sigma \approx 68.4 \mu Jy beam^{-1}$ with a spatial resolution of $0.447 \times 0.418$ (see Fig. 1). For more details on the same data reduction, we refer to Franco et al. (2018) for the high-resolution dataset (GOODS-ALMA 1.0) and
Gómez-Guijarro et al. (2022a) for the low-resolution dataset and the combined dataset (GOODS-ALMA 2.0).

## 2.2. Multiwavelength images

Here we list the multiwavelength data we used for ultra-violet (UV) to mid-infrared (MIR) and MIR to millimeter (mm) spectral energy distribution (SED) fitting, as well as those used for the stacking analysis (see Sects. 4.1, 4.2, and 5.2): (i) X-ray data: *Chandra* 7 Ms (0.5–7.0 keV, 0.5–2.0 keV, and 2–7 keV bands) images in the *Chandra* Deep Field-South (CDF-S) field (Luo et al. 2017); (ii) UV, optical (OPT), and near-infrared (NIR) data: HST/ACS (F435W, F606W, F775W, F814W, F850LP) and HST/WFC3 (F105W, F125W, F140W, F160W) images from the *Hubble* Legacy Fields Program (HLF v2.0; Whitaker et al. 2019), VLT/VIMOS (U, R; Nonino et al. 2009) images, VLT/ISAAC (H, J, Ks; Retzlaff et al. 2010) images, VLT/HAWK-I (Ks; Fontana et al. 2014) images, *Magellan/FourStar* (Hs, H1, J1, J2, J3, Ks) images from ZFOURGE (Straatman et al. 2016), and CFHT/WIRCAM (J, Ks; Hsieh et al. 2012) images; (iii) MIR data: deepest *Spitzer/IRAC* images from the GREATS program (3.6 µm, 4.5 µm, 5.8 µm, 8 µm; Stefanon et al. 2021), which were obtained by combining programs of the IUDF (PI: I. Labbé), IGOODS (PI: P. Oesch), GOODS (PI: M. Dickinson), ERS (PI: G. Fazio), S-CANDELS (PI: G. Fazio), SEDS (PI: G. Fazio), UDF2 (PI: R. Bouwens), and GREATS (PI: I. Labbé); (iv) far-infrared (FIR) data: *Spitzer/MIPS* (24 µm) images from GOODS (PI: M. Dickinson). We used *Herschel/PACS* images (100 µm, 160 µm; Magnelli et al. 2013) combined from the PEP (Lutz et al. 2011) and *GOODS-Herschel* (Elbaz et al. 2011) programs and *Herschel/SPIRE* (250 µm, 350 µm, 500 µm; Elbaz et al. 2011) images; (v) millimeter data: 1.13 mm map of GOODS-ALMA 2.0, which is a combination of high- and low-resolution 1.13 mm continuum maps (see details in Sect. 2.1); and (vi) radio data: radio image at 3 GHz (10 cm) from the VLA (PI: W. Rujopakarn, private communication), which covers the entire *GOODS-ALMA* field (Rujopakarn et al. 2016, and in prep.). In Table 1, we summarize the multiwavelength dataset used in this work from UV to mm.

## 2.3. The multiwavelength catalogs

We list below the main catalogs we used in this work. Briefly, we combined catalogs from the HST (HLF v2.0; Whitaker et al. 2019) and IRAC to select our sample, and from *Herschel* and ALMA (GOODS-ALMA 2.0; Gómez-Guijarro et al. 2022a) for the MIR-to-mm SED fitting. In addition, we also used the *GOODS-ALMA* 2.0 catalog to classify our sources detected at ALMA 1.13 mm and the X-ray catalog (CDF-S 7 Ms; Luo et al. 2017) to help identify candidate X-ray active galactic nuclei (AGN). We note that here we have corrected for systematic and local astrometric offsets in different catalogs following Franco et al. (2020a) and Whitaker et al. (2019) to ensure a consistent astrometric reference frame.

1. The HST H-band catalog: The HLF v2.0 in the *GOODS-South* region (HLF-GOODS-S) uses a deep noise-equalized combination of four HST bands (F850LP, F125W, F140W, F160W) for detection (Whitaker et al. 2019). The catalog includes all UV, optical, and NIR data (13 filters in total) taken by HST over 18 years across the field. The 5σ point-source depth in the H-band is approximately 27.0–29.8 mag. We note that the

<table><caption>**Table 1.** Broad and intermediate bands (UV to mm) in this work.</caption><thead><tr><th>Telescope/Camera</th><th>Filter</th><th>λ<sub>c</sub> (μm)</th><th>Ref.</th></tr></thead><tbody><tr><td rowspan="2">VLT/VIMOS</td><td>U</td><td>0.3759</td><td rowspan="2">(a)</td></tr><tr><td>R</td><td>0.6481</td></tr><tr><td rowspan="5">HST/ACS</td><td>F435W</td><td>0.4347</td><td rowspan="5">(b)</td></tr><tr><td>F606W</td><td>0.6033</td></tr><tr><td>F775W</td><td>0.7730</td></tr><tr><td>F814W</td><td>0.8143</td></tr><tr><td>F850LP</td><td>0.9085</td></tr><tr><td rowspan="4">HST/WFC3</td><td>F105W</td><td>1.0644</td><td rowspan="4">(b)</td></tr><tr><td>F125W</td><td>1.2561</td></tr><tr><td>F140W</td><td>1.4064</td></tr><tr><td>F160W</td><td>1.5463</td></tr><tr><td rowspan="2">CFHT/WIRCAM</td><td>J</td><td>1.2554</td><td rowspan="2">(c)</td></tr><tr><td>K<sub>s</sub></td><td>2.1630</td></tr><tr><td rowspan="3">VLT/ISAAC</td><td>J</td><td>1.2423</td><td rowspan="3">(d)</td></tr><tr><td>H</td><td>1.6560</td></tr><tr><td>K<sub>s</sub></td><td>2.1709</td></tr><tr><td>VLT/HAWK-I</td><td>K<sub>s</sub></td><td>2.1586</td><td>(e)</td></tr><tr><td rowspan="5">Magellan/FourStar</td><td>J<sub>1</sub></td><td>1.0552</td><td rowspan="5">(f)</td></tr><tr><td>J<sub>2</sub></td><td>1.1472</td></tr><tr><td>J<sub>3</sub></td><td>1.2819</td></tr><tr><td>H<sub>s</sub></td><td>1.5564</td></tr><tr><td>H<sub>1</sub></td><td>1.7038</td></tr><tr><td rowspan="4">Spitzer/IRAC</td><td>K<sub>s</sub></td><td>2.1599</td><td rowspan="4">(g)</td></tr><tr><td>CH1</td><td>3.5763</td></tr><tr><td>CH2</td><td>4.5289</td></tr><tr><td>CH3</td><td>5.7875</td></tr><tr><td></td><td>CH4</td><td>8.0449</td><td></td></tr><tr><td>Spitzer/MIPS</td><td>MIPS</td><td>24</td><td>(h)</td></tr><tr><td rowspan="3">Herschel/PACS</td><td>blue</td><td>70</td><td rowspan="3">(i)</td></tr><tr><td>green</td><td>100</td></tr><tr><td>red</td><td>160</td></tr><tr><td rowspan="4">Herschel/SPIRE</td><td>PSW</td><td>250</td><td></td></tr><tr><td>PMW</td><td>350</td><td></td></tr><tr><td>PLW</td><td>500</td><td></td></tr><tr><td>band 6</td><td>1130</td><td>(j)</td></tr></tbody></table>

**Notes.** These bands are used for the UV-to-MIR and MIR-to-mm SED fitting. From left to right: Telescope and/or instrument, the filter name, the central wavelength of the filters, and the reference for the survey, including the images and catalogs we used. (a)Nonino et al. (2009). (b)HLF program (Whitaker et al. 2019). (c)Hsieh et al. (2012). (d)Retzlaff et al. (2010). (e)Fontana et al. (2014). (f)ZFOURGE program (Straatman et al. 2016). (g)Image: the GREATS program (Stefanon et al. 2021); catalog: this work (see Appendix A). (h)PI: M. Dickinson; catalog: Magnelli et al. (2011). (i)Image: Magnelli et al. (2013) and Elbaz et al. (2011); catalog: T. Wang (priv. comm.) and Elbaz et al. (2011). (j)Image: this work (see Sect. 2.1); catalog (and also image): Gómez-Guijarro et al. (2022a).

GOODS-South field includes the HUDF, which is much deeper than other parts of the HST data, but covers only a small portion of the field.

2. IRAC catalog: We built the IRAC catalog in the *GOODS-ALMA* field using the deepest IRAC 3.6 and 4.5 µm images from the GREATS program (Stefanon et al. 2021). The source detection was performed using *Source Extractor* (SE version 2.25.0; Bertin & Arnouts 1996) on the background-subtracted 3.6 and 4.5 µm images. The total sources in the catalogs are 125 338 for 3.6 µm and 154 234 for 4.5 µm in the *GOODS-S* field (~150 arcmin²).
To ensure the purity of detections, we then cross-matched these two catalogs with a radius of 1.0" (~0.5 FWHM). Furthermore, for sources in the GOODS-S-ALMA 2.0 catalog (Gómez-Guijarro et al. 2022a) that were detected at least in one IRAC band, we also considered them to be real sources and kept them in the final catalog. Finally, we end up with 71 899 sources in the IRAC catalog, of which 5127 are in the GOODS-S-ALMA region (see Appendix A for more details).

3. *Spitzer/MIPS* 24 µm and *Herschel* catalog: We mainly used the catalog of Wang et al. (in prep.), which was built based on a state-of-the-art de-blending method, similar to that used in the 'super-deblended' catalogs in GOODS-North (Liu et al. 2018) and COSMOS (Jin et al. 2018). Briefly, Wang et al. used *Spitzer/MIPS* 24 µm detections as priors for source extraction on the PACS and SPIRE images. They performed super-deblending in one band each time, from shorter to longer wavelengths, and predicted fluxes at longer wavelengths based on the redshift and photometry information given by the shorter wavelengths. Then, the faint priors at longer wavelengths were removed, which helped break the blending degeneracies.

For one source (OFG27 in our work) still affected by blending with problematic photometry from MIR to FIR, we used the catalog of Elbaz et al. (2011).

4. **GOODS-S-ALMA 2.0 catalog:** This catalog contains 88 sources detected by ALMA at 1.13 mm in the GOODS-S-ALMA 2.0 field (Gómez-Guijarro et al. 2022a). These include 50% of sources detected above $5\sigma$ with a purity of 100% and 50% detected within the 3.5 and $5\sigma$ range aided by priors. The median redshift and stellar mass of the 88 sources are $z_{med} = 2.46$ and $\log(M_{\star\text{med}}/M_{\odot}) = 10.56$ (Chabrier 2003 IMF), respectively.

5. X-ray catalog: We used the CDF-S 7 Ms catalog (Luo et al. 2017). It contains 1008 sources in the main source catalog, observed in three X-ray bands (0.5–7.0 keV, 0.5–2.0 keV, and 2–7 keV) and 47 lower-significance sources in a supplementary catalog. This catalog includes the candidate X-ray AGNs identified by Luo et al. (2017), which we used in this work to search for X-ray AGN counterparts of our sources.

### 3. Selection of optically dark or faint galaxies at $z>3$

We aim to obtain a more complete picture of the cosmic star formation history, that is, to bridge the extreme population of optically dark galaxies (e.g., Wang et al. 2019) with the most common population of lower-mass, less-attenuated galaxies, such as those selected using the LBG selection technique (e.g., Bouwens et al. 2012a, 2015, 2020). To reach this goal, we chose to select our sample with a less strict cut than the one used in Wang et al. (2019) for the H-band and 4.5 µm magnitudes, which would allow for our sample to encompass lower-mass and less-attenuated galaxies, while still including extremely dust-obscured galaxies. Here, we call them optically dark or faint galaxies (hereafter, OFGs) and we select them with the criteria defined below (see Sect. 3.1).

#### 3.1. A selection criterion for OFGs uncontaminated by massive passive galaxies

To start, we used a combination of H-band and IRAC 4.5 µm photometry, where the H-band measures rest-frame UV light ($\lambda < 4000$ Å) for galaxies at $z > 3$, while the 4.5 µm band probes

the rest-frame J-band. Galaxies with red color in H-[4.5] are either quiescent or passive galaxies or dusty star-forming galaxies with significant dust extinction (UV attenuation).

For optically dark galaxies at high redshifts, a common selection approach in the literature is to target H-dropouts, which are defined to be undetected in the H-band (i.e., absent in the H-band catalog) and bright in the IRAC band (e.g., [4.5] < 24 mag; Wang et al. 2019; see the blue triangular region in Fig. 2), and/or extremely red in color (e.g., H-[4.5] > 4 in Shu et al. 2022). These methods can help to select extremely dust-obscured massive galaxies. However, the detection of an H-dropout obviously depends on the depth of the H-band. For example, in the HLF-GOODS-S field, the $5\sigma$ depth in the H-band ranges from 27.0 mag to 29.8 mag (Whitaker et al. 2019). To avoid this imprecise selection and to extend the selection in view of bridging the heavily obscured star-forming galaxies with more common galaxies, we defined OFGs based on the following characteristics: 1) $H > 26.5$ mag; 2) [4.5] < 25 mag (see the red triangular region in Fig. 2). Instead of only selecting galaxies undetected in the H-band and/or extremely red in H-[4.5], we used the criterion of $H > 26.5$ mag to select not only optically dark sources but also optically faint galaxies with less dust obscuration. The $H > 26.5$ mag cut also helps to distinguish massive passive galaxies with stellar masses $\log(M_{\star}/M_{\odot}) > 10$ (the grey region in Fig. 2; see more details afterward) from our selected OFGs (the red triangular region in Fig. 2). In Fig. 2, the faintest modeled passive galaxy in the H-band ($H = 26.0$ mag) has a similar magnitude to the brightest OFGs (OFG25 and OFG26: $H = 26.3 \pm 0.2$ mag), considering the $1\sigma$ uncertainty in the flux measurements. The criterion of [4.5] < 25 mag can help to select not only massive galaxies, but also galaxies with intermediate stellar masses.

To verify the reliability of our selection criteria, we investigated the evolutionary tracks of theoretical galaxy templates at $z = 2–6$ in the color-magnitude diagram (H-[4.5] vs. [4.5]; Fig. 2). The templates are based on the BC03 models (Bruzual & Charlot 2003), including an instantaneous burst model (i.e., simple stellar population: SSP) formed at $z = 10$ and a non-evolving constant star formation (CSF) model with an age of 300 Myr with varying degrees of reddening. The templates have stellar masses in the range of $\log(M_{\star}/M_{\odot}) = 10–11$, with the Calzetti et al. (2000) attenuation law and solar metallicity. The SSP model corresponds to passive or old galaxies with an age of 1 Gyr at $z = 3.5$, while the CSF models represent star-forming galaxies with different dust obscurations.

Our selection criteria for OFGs are shown in Fig. 2 as the red triangular region, which encompasses high redshift galaxies with lower stellar mass and less dust attenuation that were excluded by previous H-dropout selection criteria such as those of Wang et al. (2019; the blue triangular region in Fig. 2). For instance, the selected OFGs include those with $E(B-V) = 0.4$, $\log(M_{\star}/M_{\odot}) = 10$, and $z = 4–5$, as well as those with $E(B-V) = 0.6$, $\log(M_{\star}/M_{\odot}) = 10$, and $z = 3–4$. Similarly, extremely dust-obscured massive galaxies, such as those with $E(B-V) = 0.8$, $\log(M_{\star}/M_{\odot}) = 11$, and $z > 4$ can also be selected by our criteria. We note that although a few OFGs ($E(B-V) = 0.8$ and $\log(M_{\star}/M_{\odot}) = 10$) at $z = 2.5–3$ were selected, the total OFGs dominate at $z > 3$. Overall, in our selection, the majority of OFGs have $E(B-V) > 0.4$ and are at $z > 3$. In addition, with these criteria, our sample is not contaminated by massive passive or old galaxies ($\log(M_{\star}/M_{\odot}) > 10$; the grey region in Fig. 2). Therefore, the selection of OFGs at high redshifts with $H > 26.5$ mag and [4.5] < 25 mag is a reliable approach. In summary, the selected OFGs contain not only extremely
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Selection criteria for OFGs at $z \ge 3$. Left: color-magnitude diagram for the selection of OFGs. This diagram shows the evolutionary tracks of theoretical galaxy templates at $z = 2-6$, with stellar mass $\log(M_*/M_\odot) = 10-11$ and solar metallicity. The numbers in the circles indicate the redshift. The solid and dashed evolutionary tracks correspond to stellar masses with $\log(M_*/M_\odot) = 11$ and $\log(M_*/M_\odot) = 10$, respectively. These templates are based on the BC03 models (Bruzual & Charlot 2003), including an instantaneous burst (i.e., simple stellar population; SSP) model formed at $z = 10$ (grey area) and a non-evolving constant star formation (CSF) model with an age of 300 Myr with varying degrees of reddening (blue, green, olive, and red areas). The SSP model corresponds to a passive or old galaxy with an age of 1 Gyr at $z = 3.5$. The blue triangular region shows the selection criteria for H-dropouts in Wang et al. (2019), i.e., no H-band flux above $5\sigma$ ($H > 27$ mag; blue dashed line) and $[4.5] < 24$ mag (blue solid line). The red triangular region shows the selection criteria for OFGs in this work. Right: different theoretical galaxy templates in HST/F160W (H-band) and IRAC 4.5 µm filters. The templates include a passive galaxy with an age of 1 Gyr at $z = 3.5$ and three star-forming galaxies with different degrees of reddening at $z = 4$. At a similar [4.5], the three star-forming galaxies with different dust extinction values exhibit different H-band magnitudes.</figcaption></figure>

dust-obscured massive galaxies at $z > 4$, but also lower-mass and less-attenuated (typically $E(B-V) > 0.4$) galaxies at $z > 3$, without contamination from massive passive galaxies.

## 3.2. Sample selection

We selected candidate OFGs by cross-matching our IRAC catalog at 3.6 and 4.5 µm (see Appendix A for details on the catalog construction) with the HLF catalog (Whitaker et al. 2019) in the GOODS-ALMA field. Candidates were required to have none or very faint HST counterparts ($H > 26.5$ mag) within a 0.6" radius around the IRAC positions. The small radius of 0.6" (roughly one-third of the point spread function size of IRAC 4.5 µm) was taken to select as many individual candidates as possible, while avoiding excessive contamination of our final sample by fake candidates. After cross-matching, we had 88 candidates.

## 3.3. Photometry

We visually inspected all the candidates and noted that blending was common in the IRAC images due to their relatively worse spatial resolution (e.g., ~2" at 4.5 µm). To obtain photometric values in different bands without contamination from neighboring galaxies, we simultaneously de-blended sources in the multi-wavelength images (from UV to 8 µm) by applying the de-blending code<fnref n="1" /> described in Schreiber et al. (2018b). The de-blending method is briefly summarized in the following three steps. Step 1: for each OF candidate, to save computational time, we first cut the stacked HST image (four bands

of HST/F105W, F125W, F140W, and F160W) into a 10" × 10" area around the IRAC position of the candidate. Then, we detected all sources in the clipped stacked HST image.

For some optically dark galaxies undetected even in the stacked HST image, we used positions of their $K_s$-band, 4.5 µm, or ALMA counterparts and modeled them as point sources. Step 2: following Schreiber et al. (2018b), we fitted all the sources detected in the clipped stacked HST image simultaneously with a single Sérsic profile to obtain a best-fit deconvolved model (intrinsic light profile) for each source. Step 3: these models were then convolved with a point spread function (PSF) for each image at all wavelengths (up to 8 µm). We then used the positions and the PSF-convolved models as priors for all objects to fit the multi-wavelength images. The uncertainties of the fluxes were calculated by Monte Carlo simulations. We find that 60% (53/88) of the sources in our candidate sample needed to be de-blended at 4.5 µm.

## 3.4. Incompleteness correction and final sample

To identify high-$z$ OFGs, we used the simple selection techniques discussed in Sect. 3.1: (1) $H > 26.5$ mag and (2) $[4.5] < 25$ mag. Considering the $1\sigma$ uncertainty of the flux measurements, we finally identified 26 individual OFGs in total in the GOODS-ALMA field (see Figs. 1 and 3).

Here, we discuss the corrections for the incompleteness of our sample selection approach. Considering the criterion of no/very faint HST counterparts within the search radius of 0.6" at the IRAC position, we may have missed some target sources simply due to random bright HST sources falling within the

<footnote n="1">¹ qdeblend: https://github.com/cschreib/qdeblend</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 3. Color-magnitude diagram color-coded by photometric redshift. Our criteria (in red) for selecting OFGs are: $H > 26.5$ mag & [4.5] < 25 mag. We note that we include sources outside the wedge whose $1\sigma$ photometric uncertainties overlap the wedge, so we have points outside the red triangle. The arrows show our H-dropouts ($S/N < 3$), with the typical depth of $H = 27$ mag ($5\sigma$) in the shallowest region of the HLF survey as their lower limits. The blue shaded area describes the distribution of all the IRAC detected sources (see Appendix A for more details) in the GOODS-ALMA field. The blue and red triangular regions are the same as in Fig. 2.</figcaption></figure>

radius. Following the same method as in Lilly et al. (1999); also in Wang et al. 2019), at a given position, the probability of finding one or more random galaxies within a given radius is defined as:

$$p(n \ge 1) = 1 - \exp(-N \times \pi \times \text{radius}^2), \quad (1)$$

where N represents the surface density of bright HST sources ($H < 26.5$ mag) in our case. In the GOODS-ALMA field, N is 0.05 galaxies arcsec$^{-2}$, so the derived $p(n \ge 1) = 0.056$. It suggests that using this selection approach, we may have missed around 5% of OFGs that we wrongly associated with a counterpart due to projection effects. That is, one source could have been missed due to the serendipitous presence of a bright H detection within our 0.6" radius search circle.

After comparing our sample with 13 ALMA-detected OFGs in the GOODS-ALMA 2.0 catalog, we confirmed that the missing galaxy is OFG27 (A2GS7 in Gómez-Guijarro et al. 2022a). A random bright source (ID55970 in the HLF catalog) with $H = 24.5$ mag is located at a distance of 0.33" (<0.6" searching radius) from OFG27. This bright source has been confirmed not to be the H-band counterpart of OFG27 (AGS17 in Zhou et al. 2020). Therefore, we included OFG27 in our catalog (see Table 2) to correct the incompleteness of our selection approach. Also, we included OFG27 in the analysis of the main discussion presented in this paper.

In addition, an extra IRAC 4.5 µm dropout candidate was detected only in the longer wavelength images: JCMT/SCUBA-2 850 µm and ALMA 870 µm and 1.13 mm and 1.2 mm (OFG28; see Table 2). Including this one, we have 28 OFGs in our final catalog (Table 2). Considering that OFG28 is an IRAC 4.5 µm dropout candidate with [4.5] ≫ 25 mag, in the following analysis, we focus only on the first 27 OFGs, which meet our criteria:

$$H > 26.5 \text{ mag} \text{ & [4.5] < 25 mag.}$$

### 3.5. Lyman-break galaxies and H-dropouts in final sample

To compare our sample with LBGs (e.g., Bouwens et al. 2012a) and H-dropouts (e.g., Wang et al. 2019), we need to determine how many OFGs are LBGs or H-dropouts in our sample.

#### 3.5.1. Lyman-break galaxies

LBGs are a UV-selected population of high-z star-forming galaxies. To understand how many galaxies in our OFG catalog are missed by this UV-selected approach and to further know their contribution to the cosmic SFRD, we identify LBGs from our OFG catalog by employing the Lyman-break color criteria used in Bouwens et al. (2020; also see similar methods in Bouwens et al. 2012a, 2015). The redshift range of our OFGs is $z = 3-7$ (see Sect. 4.1, Fig. 4, and Table 2). The Lyman-break color criteria are as follows:

$$z \sim 3 : (U_{336} - B_{435} > 1) \land (B_{435} - V_{606} < 1.2) \land (i_{775} - Y_{105} < 0.7) \land (\chi^2_{UV_{225},UV_{275}} < 2),$$

$$z \sim 4 : (B_{435} - V_{606} > 1) \land (i_{775} - J_{125} < 1) \land (B_{435} - V_{606} > 1.6(i_{775} - J_{125}) + 1),$$

$$z \sim 5 : (V_{606} - i_{775} > 1.2) \land (z_{850} - H_{160} < 1.3) \land (V_{606} - i_{775} > 0.8(z_{850} - H_{160}) + 1.2),$$
<table><caption>**Table 2.** Derived properties of the OFGs.</caption><thead><tr><th>ID<br>(1)</th><th>RA<br>(deg)<br>(2)</th><th>Dec<br>(deg)<br>(3)</th><th>H<br>(mag)<br>(4)</th><th>[4.5]<br>(mag)<br>(5)</th><th>S<sub>1.13 mm</sub><br>(mJy)<br>(6)</th><th>z<br>(7)</th><th>log(M<sub>*</sub>)<br>(8)</th><th>log(L<sub>IR</sub>)<br>(9)</th><th>SFR<br>(M<sub>☉</sub> yr<sup>-1</sup>)<br>(10)</th><th>Other ID<br>(11)</th></tr></thead><tbody><tr><td>OFG1</td><td>53.087184<sup>(†)</sup></td><td>-27.840242<sup>(†)</sup></td><td>(...)</td><td>22.10<sup>+0.11</sup><sub>-0.10</sub></td><td>0.85 ± 0.11</td><td>3.47<sup>o</sup></td><td>10.79<sup>+0.11</sup><sub>-0.16</sub></td><td>12.41 ± 0.05</td><td>384 ± 48</td><td>AGS24, A2GS29</td></tr><tr><td>OFG2</td><td>53.108810<sup>(†)</sup></td><td>-27.869037<sup>(†)</sup></td><td>(...)</td><td>23.72<sup>+0.11</sup><sub>-0.10</sub></td><td>1.24 ± 0.10</td><td>3.47<sup>o</sup></td><td>10.76<sup>+0.11</sup><sub>-0.17</sub></td><td>12.39 ± 0.03</td><td>365 ± 30</td><td>AGS11, A2GS15</td></tr><tr><td>OFG3</td><td>53.102536</td><td>-27.806531</td><td>(...)</td><td>23.75<sup>+0.11</sup><sub>-0.10</sub></td><td>(...)</td><td>7.04<sup>+1.01</sup><sub>-1.19</sub></td><td>10.91<sup>+0.22</sup><sub>-0.14</sub></td><td>(...)</td><td>(...)</td><td>(...)</td></tr><tr><td>OFG4</td><td>53.206064<sup>(†)</sup></td><td>-27.819142<sup>(†)</sup></td><td>(...)</td><td>25.16<sup>+0.11</sup><sub>-0.10</sub></td><td>0.70 ± 0.13</td><td>6.13<sup>+1.02</sup><sub>-1.36</sub></td><td>10.60<sup>+0.16</sup><sub>-0.36</sub></td><td>12.56 ± 0.07</td><td>537 ± 100</td><td>A2GS38</td></tr><tr><td>OFG5</td><td>53.119150<sup>(†)</sup></td><td>-27.814066<sup>(†)</sup></td><td>(...)</td><td>23.34<sup>+0.11</sup><sub>-0.10</sub></td><td>0.43 ± 0.13</td><td>3.81<sup>+0.74</sup><sub>-0.94</sub></td><td>10.17<sup>+0.34</sup><sub>-0.27</sub></td><td>12.26 ± 0.09</td><td>268 ± 61</td><td>A2GS87, GDS44539</td></tr><tr><td>OFG6</td><td>53.197493</td><td>-27.813789</td><td>(...)</td><td>23.23<sup>+0.11</sup><sub>-0.10</sub></td><td>(...)</td><td>4.10<sup>+0.31</sup><sub>-0.32</sub></td><td>10.93<sup>+0.06</sup><sub>-0.06</sub></td><td>(...)</td><td>(...)</td><td>(...)</td></tr><tr><td>OFG7</td><td>53.183697<sup>(†)</sup></td><td>-27.836500<sup>(†)</sup></td><td>(...)</td><td>24.55<sup>+0.11</sup><sub>-0.10</sub></td><td>1.23 ± 0.12</td><td>4.58<sup>+0.42</sup><sub>-0.42</sub></td><td>10.34<sup>+0.08</sup><sub>-0.09</sub></td><td>12.86 ± 0.04</td><td>1070 ± 101</td><td>AGS25, A2GS17</td></tr><tr><td>OFG8</td><td>53.210736</td><td>-27.813706</td><td>27.64<sup>+0.31</sup><sub>-0.24</sub></td><td>24.93<sup>+0.60</sup><sub>-0.40</sub></td><td>(...)</td><td>4.14<sup>+0.41</sup><sub>-0.44</sub></td><td>9.92<sup>+0.18</sup><sub>-0.27</sub></td><td>(...)</td><td>(...)</td><td>(...)</td></tr><tr><td>OFG9</td><td>53.191204</td><td>-27.835791</td><td>27.61<sup>+0.22</sup><sub>-0.19</sub></td><td>24.82<sup>+0.11</sup><sub>-0.10</sub></td><td>(...)</td><td>3.99<sup>+0.18</sup><sub>-0.17</sub></td><td>9.59<sup>+0.08</sup><sub>-0.08</sub></td><td>(...)</td><td>(...)</td><td>(...)</td></tr><tr><td>OFG10</td><td>53.191348</td><td>-27.737554</td><td>27.54<sup>+0.23</sup><sub>-0.19</sub></td><td>24.25<sup>+0.11</sup><sub>-0.10</sub></td><td>(...)</td><td>5.04<sup>+0.34</sup><sub>-0.31</sub></td><td>10.37<sup>+0.06</sup><sub>-0.08</sub></td><td>(...)</td><td>(...)</td><td>(...)</td></tr><tr><td>OFG11</td><td>53.196569<sup>(†)</sup></td><td>-27.757065<sup>(†)</sup></td><td>27.44<sup>+0.35</sup><sub>-0.27</sub></td><td>23.59<sup>+0.11</sup><sub>-0.10</sub></td><td>0.62 ± 0.12</td><td>3.60<sup>+0.27</sup><sub>-0.27</sub></td><td>10.41<sup>+0.06</sup><sub>-0.07</sub></td><td>12.41 ± 0.09</td><td>384 ± 90</td><td>A2GS40, GDS48885</td></tr><tr><td>OFG12</td><td>53.154787</td><td>-27.806529</td><td>27.36<sup>+0.17</sup><sub>-0.15</sub></td><td>24.43<sup>+0.11</sup><sub>-0.10</sub></td><td>(...)</td><td>5.55<sup>+0.29</sup><sub>-0.28</sub></td><td>10.31<sup>+0.05</sup><sub>-0.06</sub></td><td>(...)</td><td>(...)</td><td>(...)</td></tr><tr><td>OFG13</td><td>53.047834</td><td>-27.829186</td><td>27.34<sup>+0.24</sup><sub>-0.20</sub></td><td>25.07<sup>+0.11</sup><sub>-0.10</sub></td><td>(...)</td><td>3.56<sup>+0.09</sup><sub>-0.09</sub></td><td>9.50<sup>+0.05</sup><sub>-0.03</sub></td><td>(...)</td><td>(...)</td><td>(...)</td></tr><tr><td>OFG14</td><td>53.105489</td><td>-27.830711</td><td>27.34<sup>+0.12</sup><sub>-0.11</sub></td><td>24.38<sup>+0.11</sup><sub>-0.10</sub></td><td>(...)</td><td>3.39<sup>+0.92</sup><sub>-0.57</sub></td><td>9.86<sup>+0.23</sup><sub>-0.22</sub></td><td>(...)</td><td>(...)</td><td>(...)</td></tr><tr><td>OFG15</td><td>53.132675</td><td>-27.765496</td><td>27.08<sup>+0.16</sup><sub>-0.14</sub></td><td>24.23<sup>+0.11</sup><sub>-0.10</sub></td><td>(...)</td><td>3.192<sup>sp</sup></td><td>9.96<sup>+0.04</sup><sub>-0.05</sub></td><td>(...)</td><td>(...)</td><td>(...)</td></tr><tr><td>OFG16</td><td>53.080379</td><td>-27.869420</td><td>26.93<sup>+0.29</sup><sub>-0.23</sub></td><td>23.58<sup>+0.11</sup><sub>-0.10</sub></td><td>(...)</td><td>3.69<sup>+0.27</sup><sub>-0.25</sub></td><td>10.38<sup>+0.08</sup><sub>-0.08</sub></td><td>(...)</td><td>(...)</td><td>(...)</td></tr><tr><td>OFG17</td><td>53.062276</td><td>-27.875036</td><td>26.84<sup>+0.15</sup><sub>-0.13</sub></td><td>23.93<sup>+0.11</sup><sub>-0.10</sub></td><td>(...)</td><td>4.23<sup>+0.23</sup><sub>-0.21</sub></td><td>10.29<sup>+0.07</sup><sub>-0.05</sub></td><td>(...)</td><td>(...)</td><td>(...)</td></tr><tr><td>OFG18</td><td>53.188278<sup>(†)</sup></td><td>-27.801928<sup>(†)</sup></td><td>26.82<sup>+0.14</sup><sub>-0.12</sub></td><td>24.88<sup>+0.11</sup><sub>-0.10</sub></td><td>0.37 ± 0.11</td><td>3.81<sup>+0.11</sup><sub>-0.11</sub></td><td>9.44<sup>+0.03</sup><sub>-0.11</sub></td><td>12.28 ± 0.11</td><td>286 ± 86</td><td>A2GS47</td></tr><tr><td>OFG19</td><td>53.162978<sup>(†)</sup></td><td>-27.841940<sup>(†)</sup></td><td>26.76<sup>+0.20</sup><sub>-0.17</sub></td><td>24.14<sup>+0.11</sup><sub>-0.10</sub></td><td>0.40 ± 0.12</td><td>4.09<sup>+0.35</sup><sub>-0.32</sub></td><td>10.30<sup>+0.06</sup><sub>-0.09</sub></td><td>12.26 ± 0.09</td><td>271 ± 59</td><td>A2GS82</td></tr><tr><td>OFG20</td><td>53.064807<sup>(†)</sup></td><td>-27.862613<sup>(†)</sup></td><td>26.64<sup>+0.27</sup><sub>-0.22</sub></td><td>22.85<sup>+0.11</sup><sub>-0.10</sub></td><td>0.54 ± 0.10</td><td>4.74<sup>+0.42</sup><sub>-0.50</sub></td><td>10.88<sup>+0.16</sup><sub>-0.26</sub></td><td>12.79 ± 0.04</td><td>913 ± 86</td><td>A2GS57</td></tr><tr><td>OFG21</td><td>53.060144</td><td>-27.793838</td><td>26.59<sup>+0.12</sup><sub>-0.12</sub></td><td>24.68<sup>+0.11</sup><sub>-0.10</sub></td><td>(...)</td><td>5.89<sup>+0.36</sup><sub>-0.38</sub></td><td>9.75<sup>+0.32</sup><sub>-0.13</sub></td><td>(...)</td><td>(...)</td><td>(...)</td></tr><tr><td>OFG22</td><td>53.043745</td><td>-27.804347</td><td>26.51<sup>+0.19</sup><sub>-0.16</sub></td><td>24.56<sup>+0.12</sup><sub>-0.11</sub></td><td>(...)</td><td>4.56<sup>+0.26</sup><sub>-0.24</sub></td><td>9.97<sup>+0.07</sup><sub>-0.07</sub></td><td>(...)</td><td>(...)</td><td>(...)</td></tr><tr><td>OFG23</td><td>53.081890</td><td>-27.828815</td><td>26.45<sup>+0.11</sup><sub>-0.10</sub></td><td>23.83<sup>+0.11</sup><sub>-0.10</sub></td><td>(...)</td><td>3.88<sup>+0.18</sup><sub>-0.18</sub></td><td>10.25<sup>+0.04</sup><sub>-0.04</sub></td><td>(...)</td><td>(...)</td><td>(...)</td></tr><tr><td>OFG24</td><td>53.109771</td><td>-27.807466</td><td>26.40<sup>+0.11</sup><sub>-0.10</sub></td><td>25.07<sup>+0.11</sup><sub>-0.10</sub></td><td>(...)</td><td>6.27<sup>+0.16</sup><sub>-0.16</sub></td><td>9.87<sup>+0.10</sup><sub>-0.08</sub></td><td>(...)</td><td>(...)</td><td>(...)</td></tr><tr><td>OFG25</td><td>53.074868<sup>(†)</sup></td><td>-27.875889<sup>(†)</sup></td><td>26.31<sup>+0.30</sup><sub>-0.24</sub></td><td>23.10<sup>+0.11</sup><sub>-0.10</sub></td><td>1.67 ± 0.10</td><td>3.47<sup>o</sup></td><td>9.99<sup>+0.46</sup><sub>-0.11</sub></td><td>12.92 ± 0.03</td><td>1227 ± 74</td><td>AGS15, A2GS10</td></tr><tr><td>OFG26</td><td>53.207252</td><td>-27.791408</td><td>26.29<sup>+0.22</sup><sub>-0.19</sub></td><td>23.63<sup>+0.11</sup><sub>-0.10</sub></td><td>(...)</td><td>4.16<sup>+0.25</sup><sub>-0.25</sub></td><td>10.33<sup>+0.08</sup><sub>-0.09</sub></td><td>(...)</td><td>(...)</td><td>(...)</td></tr><tr><td>OFG27<sup>(*)</sup></td><td>53.079416<sup>(†)</sup></td><td>-27.870820<sup>(†)</sup></td><td>(...)</td><td>21.53<sup>+0.11</sup><sub>-0.10</sub></td><td>2.05 ± 0.12</td><td>3.467<sup>sp</sup></td><td>11.11<sup>+0.15</sup><sub>-0.19</sub></td><td>13.08 ± 0.02</td><td>1795 ± 90</td><td>AGS17, A2GS7</td></tr><tr><td>OFG28<sup>(**)</sup></td><td>53.120402<sup>(†)</sup></td><td>-27.742111<sup>(†)</sup></td><td>(...)</td><td>(...)</td><td>0.95 ± 0.12</td><td>(...)</td><td>(...)</td><td>(...)</td><td>(...)</td><td>A2GS33</td></tr></tbody></table>

**Notes.** (1) Source ID; (2)(3) Right ascension and declination (J2000) of sources. Coordinates detected in the ALMA 1.13 mm image are marked with a “(†)” exponent; (4)(5) H-band and IRAC 4.5 µm AB magnitudes. These magnitudes are given for the best-fitting Sérsic profile during de-blending procedure (see Sect. 3.3); (6) ALMA 1.13 mm flux density: obtained from the GOODS-ALMA 2.0 catalog (Gómez-Guijarro et al. 2022a); (7) Photometric redshifts: determined with the EAZY code (see Sect. 4.1; spectroscopic redshifts expressed in three decimal places and flagged with a “sp” exponent). The spectroscopic redshifts of OFG15 and OFG27 are from Herenz et al. (2017) and Zhou et al. (2020), respectively. OFG1, 2, 25, 27 were discovered in an overdensity region with the redshift peak to be 3.47 (flagged with a “o” exponent; Zhou et al. 2020); (8) Stellar masses: determined with the FAST++ (see Sect. 4.1); (9) Infrared luminosities: derived from CIGALE for three sources (OFG2, OFG20, and OFG27) with a Herschel counterpart or from the IR template library (Schreiber et al. 2018c) for the galaxies without a Herschel counterpart (see Sect. 4.2); (10) SFR = SFR<sub>IR</sub> + SFR<sub>UV</sub> (see Sect. 4.3); (11) Source IDs in other work: AGS (GOODS-ALMA 1.0 catalog; Franco et al. 2018, 2020a); A2GS (GOODS-ALMA 2.0 catalog; Gómez-Guijarro et al. 2022a); GDS (H-dropouts catalog; Wang et al. 2019). (*) The missing galaxy in our selection approach is due to incompleteness of the search radius of 0.6". We add it to our catalog to correct this incompleteness. (**) OFG28 is a candidate IRAC 4.5 µm dropout, which is only detected in longer wavelength images, e.g., 850 µm from the JCMT/SCUBA-2 and 870 µm from the ALMA (ID68 in Cowie et al. 2018), 1.13 mm from the GOODS-ALMA 2.0 100% pure source catalog (A2GS33 in Gómez-Guijarro et al. 2022a), and 1.2 mm from the ALMA (ID20 in Yamaguchi et al. 2019). We add it here to refine our OFG catalog in the GOODS-ALMA field. We note that OFG28 is not used in our analysis. We also note that OFG9, OFG13, OFG18, OFG22, and OFG24 are identified as LBGs (see Sect. 3.5.1). In this catalog, eight sources (OFG1-OFG7, and OFG27) are H-dropouts, which are not detected in the H-band (<3σ).

$$z \sim 6: (i_{775} - z_{850} > 1.0) \wedge (Y_{105} - H_{160} < 1.0) \wedge \\ (i_{775} - z_{850} > 0.77(Y_{105} - H_{160}) + 1.0),$$

$$z \sim 7: (z_{850} - Y_{105} > 0.7) \wedge (J_{125} - H_{160} < 0.45) \wedge \\ (z_{850} - Y_{105} > 0.8(J_{125} - H_{160}) + 0.7),$$

here, ∧ and ∨ represent the logical AND and OR symbols, respectively. The $\chi^2_{UV_{225,UV_{275}}}=\Sigma_i \text{SGN}(f_i)(f_i/\sigma_i)^2$, where $f_i(\sigma_i)$ is the flux (uncertainty) in the i-band and SGN($f_i$) is equal to 1 if $f_i > 0$ and -1 if $f_i < 0$. As in Bouwens et al. (2015), we use a $1\sigma$ upper limit as the flux in the dropout band in the case of non-detection. The selected sources are required to be

undetected (<2σ) in all bands blueward of the Lyman break and detected (>3σ) in all of the above bands redward of the break. We note that we do not include $U_{336}$, $UV_{225}$, and $UV_{275}$ in our work because our OFGs are undetected in these bands. Even if we only use the criteria of $z \sim 3$: ($B_{435} - V_{606} < 1.2$) ∧ ($i_{775} - Y_{105} < 0.7$), no galaxy at $z \sim 3$ in our catalog is classified as LBG. The color criteria for Lyman-break of Bouwens et al. (2020) are slightly different from those of Bouwens et al. (2012a, 2015), thus we also used the criteria of Bouwens et al. (2012a, 2015) to select LBGs. All three methods identified the same 5 LBGs in our OFG catalog: OFG9, OFG13, and OFG18 at $z \sim 4$; OFG22 at $z \sim 5$; and OFG24 at $z \sim 6$.
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Photometric redshift and stellar mass distributions for our total 27 OFGs (grey). Among them, five OFGs are identified as LBGs (blue; see Sect. 3.5.1) and seven OFGs meet the criteria of H-dropouts ($H > 27$ mag & [4.5] < 24 mag) in Wang et al. (2019; red; see Sect. 3.5.2). *Top*: photometric redshifts of our OFGs (including two sources with spectroscopic redshifts). *Bottom*: LBGs cover the low stellar mass end and H-dropouts cover the high stellar mass end. Our OFGs have a wide range of stellar masses with $log(M_*/M_\odot) = 9.4–11.1$. The median redshift and stellar mass of the OFGs (dashed lines) are $z_{med} = 4.1$ and $log(M_{*med}/M_\odot) = 10.3$, respectively.</figcaption></figure>

### 3.5.2. H-dropouts

Galaxies bright in IRAC but not detected in the H-band are commonly referred to as H-dropouts. However, this definition can be confusing since different fields have been observed at different depths in the H-band. Here, we used the deepest H-band image to date in the GOODS-South field (HLF; Whitaker et al. 2019), with a $5\sigma$ point-source depth of approximately 27.0–29.8 mag, to identify H-dropouts and to extend the sample to our more general definition of OFGs. In our sample, eight galaxies (OFG1-OFG7 and OFG27) are classified as H-dropouts, that is, there is no detection above $3\sigma$ in the H-band (see Table 2). We note that in Wang et al. (2019), H-dropouts include all sources with no H-band flux above $5\sigma$, that is, $H > 27$ mag – instead of $3\sigma$ here – and [4.5] < 24 mag. If we apply the same criterion, we find seven OFGs that meet this definition: OFG1, OFG2, OFG3, OFG5, OFG6, OFG11, and OFG27 (see Fig. 3).

In the GOODS-ALMA 2.0 catalog, Gómez-Guijarro et al. (2022a) reported 13 OFGs out of 88 galaxies detected above $3.5\sigma$ at 1.13 mm in the GOODS-ALMA field. Among them,

12 OFGs are included in our sample (see Table 2 with a “†” exponent). The remaining one (A2GS2 in Gómez-Guijarro et al. 2022a or AGS4 in Franco et al. 2018, 2020a; Zhou et al. 2020) does not meet our criterion of $H > 26.5$ mag, since after applying our de-blending procedure (Sect. 3.3), we measured an H-band magnitude of $H = 24.76$ mag. This value is consistent with the findings of Zhou et al. (2020), who measured $H = 25.23$ mag.

## 4. Properties of individual galaxies

In this section, we focus on a set of properties of individual galaxies: redshift, stellar mass, infrared luminosity, star formation rate, gas mass, dust mass, and dust temperature. The methodologies used to derive these properties are also used for the stacked samples, as described in Sect. 5.

### 4.1. Redshifts and stellar masses

We fit the SED from the UV to MIR (rest-frame UV to NIR) to measure the photometric redshifts ($z_{phot}$) and stellar masses ($M_*$). The photometric redshifts were determined with the code EAzY<fnref n="2" /> (Brammer et al. 2008). Then we fixed the redshift and derived $M_*$ with the code FAST++<fnref n="3" />. The setup is described below.

Photometric redshifts were obtained with the galaxy template set “eazy_v1.3”, which includes, in particular, a dusty starburst model to account for extremely dusty galaxies. We did not apply the redshift prior based on K-band magnitudes, as this prior is based on models that do not reproduce high-redshift mass functions (see discussion in Schreiber et al. 2018a).

Two sources, OFG15 and OFG27, have spectroscopic redshifts ($z_{spec}$) confirmed by one-line detections. The galaxy OFG15 has a $z_{spec}$ measured by the Lyman $\alpha$ line (~7$\sigma$) from the MUSE-Wide survey (ID139013229 in Herenz et al. 2017; Urrutia et al. 2019). The galaxy OFG27 has a $z_{spec}$ identified from the CO(6-5) line detection (~10$\sigma$) with ALMA (AGS17 in Zhou et al. 2020). In addition, OFG1, 2, 25, and 27 were discovered in an overdensity region with a peak redshift of $z_{overdensity} = 3.47$ (AGS24, 11, 15, 17 in Zhou et al. 2020, where they were studied in detail). For the above five galaxies, we used their $z_{spec}$ or $z_{overdensity}$ in the following analysis.

Stellar masses were then derived using the code FAST++, assuming a delayed, exponentially declining star formation history (SFH), with Bruzual & Charlot (2003) stellar population models and a Calzetti et al. (2000) dust attenuation law. The parameters used in FAST++ are shown in Table 3.

In Fig. 4, we show the distributions of derived redshifts and stellar masses of our OFGs. The redshift distribution confirms that our OFGs exhibit redshifts of $z_{phot} > 3$, which are consistent with the theoretical galaxy templates (see Fig. 2). The median redshift of the distribution is $z_{med} = 4.1$. Compared to the LBGs covering the low stellar mass end and the H-dropouts in Wang et al. (2019) covering the high stellar mass end, our sample presents a broad distribution of stellar masses with $log(M_*/M_\odot) = 9.4–11.1$ and a median value of $log(M_{*med}/M_\odot) = 10.3$. The individual redshift and stellar mass values are listed in Table 2, and the individual SEDs are presented in Figs. C.1 and C.2.

We investigate in Fig. 5, the proportions of LBGs, H-dropouts, and remaining OFGs (after removing LBGs and

<footnote n="2">² EAZY: https://github.com/gbrammer/eazy-photoz</footnote>

<footnote n="3">³ FAST++: https://github.com/cschreib/fastpp, an updated version of the SED fitting code FAST (Kriek et al. 2009).</footnote>
<table><caption>**Table 3.** Input parameters used in the UV to 8 µm SED fitting procedures with FAST++ to derive stellar masses.</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td></td><td>Delayed SFH</td></tr><tr><td>Age [log(yr<sup>-1</sup>)]</td><td>6.0–10.2, step 0.1</td></tr><tr><td>τ [log(yr<sup>-1</sup>)]</td><td>6.5–11, step 0.1</td></tr><tr><td>Metallicity</td><td>0.02 (solar)</td></tr><tr><td colspan="2">Dust attenuation: Calzetti et al. (2000)</td></tr><tr><td>A<sub>v</sub></td><td>0–6, step 0.02</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Number fractions of galaxies in our OFG sample that are identified as LBGs (blue shaded region), H-dropouts (red shaded region), and the remaining OFGs (orange filled region) that are missed by the first two selection techniques. The OFGs undetected by the LBGs and H-dropouts criteria represent approximately 75% at stellar masses between $log(M_{\star}/M_{\odot}) = 9.5–10.5$.</figcaption></figure>

<table><caption>**Table 4.** Input parameters used in the 24 µm to mm SED fitting procedures with CIGALE.</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td colspan="2">Dust emission: Draine et al. (2014)</td></tr><tr><td>q<sub>PAH</sub></td><td>0.47, 1.12, 1.77, 2.50, 3.19, 3.90, 4.58, 5.26, 5.95, 6.63, 7.32</td></tr><tr><td>U<sub>min</sub></td><td>0.1, 0.5, 1, 5, 12, 15, 20, 25, 30, 35, 40, 50</td></tr><tr><td>α</td><td>1, 1.5, 2, 2.5, 3.0</td></tr><tr><td>γ</td><td>0.0001, 0.001, 0.01, 0.1, 0.5, 1</td></tr></tbody></table>

H-dropouts) in our sample at different stellar masses. At stellar masses of $log(M_{\star}/M_{\odot}) = 9.5–10.5$, the fraction of OFGs is about three times the sum of LBGs and H-dropouts. In other words, up to 75% of the galaxies with a stellar mass of $log(M_{\star}/M_{\odot}) = 9.5–10.5$ at $z > 3$ are missed by the previous LBG and H-dropout selection techniques.

## 4.2. Infrared SED fitting

The infrared luminosities (8–1000 µm; $L_{IR}$) and dust mass ($M_{dust}$) are derived from the FIR SED-fitting. We used two

methodologies in the fit, depending on whether the galaxies have a Herschel counterpart or not (see Figs. C.1 and C.2 for the entire sample).

For galaxies with a Herschel counterpart (3/27), we performed the FIR SED-fitting with CIGALE<fnref n="4" /> (Code Investigating GALaxy Emission; Burgarella et al. 2005; Noll et al. 2009; Boquien et al. 2019). We fit data from 24 µm up to millimeter wavelengths from the catalogs of Wang et al. (in prep.), Elbaz et al. (2011), and GOODS-ALMA v2.0 1.13 mm (Gómez-Guijarro et al. 2022a). The dust infrared emission model is the one of Draine et al. (2014). The parameters used in CIGALE are shown in Table 4.

For galaxies without a Herschel counterpart but with ALMA detections at 1.13 mm (8/27; Gómez-Guijarro et al. 2022a), applying the dust emission model (Draine et al. 2014) of CIGALE would fit a single point in the FIR with four parameters, which would leave us with much less meaningful results. As a compromise, we used the IR template library<fnref n="5" /> described in Schreiber et al. (2018c). In brief, it consists of two ingredients: i) dust continuum created by big dust grains (silicate + amorphous carbon grains) and ii) mid-infrared features contributed by polycyclic aromatic hydrocarbon (PAH) molecules. To form a full dust spectrum, the relative contribution of big grains and PAHs are characterized by the mid-to-total infrared color ($IR8 = L_{IR}/L_8$), which is related to the redshift and starburstiness of a galaxy. The starburstiness is defined by $R_{SB} \equiv \Delta MS \equiv SFR/SFR_{MS}$ (Elbaz et al. 2011), where $SFR_{MS}$ is the average SFR of the main sequence galaxies presented in Schreiber et al. (2015).

We fit the data iteratively with two templates: the star formation main sequence (MS) template and the starburst (SB) template, with fixed $R_{SB}$ values of 1 and 5, respectively. For each template, given the known redshift and fixed $R_{SB}$, we calculated the dust temperature ($T_{dust}$) and IR8 using Eqs. (18) and (19) from Schreiber et al. (2018c). The templates we used were normalized to $M_{dust} = 1 M_{\odot}$. After re-normalizing the SED to the ALMA flux density at 1.13 mm, we obtained the total $M_{dust}$ and total $L_{IR}$ by integrating the SED in the 8–1000 µm rest-frame range. Then we computed $R_{SB}$ using the output $L_{IR}$ (Kennicutt & Evans 2012; the contribution of UV to the SFR is negligible, as shown later in Sect. 5.3.1). We computed two $R_{SB}$ values for each galaxy derived from both templates (MS and SB). If both $R_{SB}$ values are less (greater) than 3, we consider a galaxy a MS (SB) galaxy. The best-fit SED was then generated with the typical value $R_{SB} = 1$ (MS) or $R_{SB} = 5$ (SB). Otherwise, that is, if both templates do not agree with each other for $R_{SB} < 3$ or $>3$, we kept the two SEDs given by both templates as upper and lower limits and used the average template as the best SED. This approach is similar to the one used by Gómez-Guijarro et al. (2022b) but slightly more conservative. Compared to the $L_{IR}$ in Gómez-Guijarro et al. (2022b) for the same sources, the results are generally consistent, with a median relative difference of $(L_{IR}^{This work} - L_{IR}^{G22})/L_{IR}^{G22} = 0.04 \pm 0.18$. The relatively large dispersion was expected because the IR template fit is based on only one observed data point of ALMA 1.13 mm with large uncertainty in $z_{phot}$.

As a consistency test, for the three sources detected by Herschel and ALMA, we performed a SED fit using the IR template library normalized only to the ALMA point as if they had no Herschel values. We then derived $L_{IR}$. The ratios of the IR luminosities, $L_{IR}^{CIGALE}/L_{IR}^{IR template}$ (where $L_{IR}^{IR template}$ is derived

<footnote n="4">⁴ CIGALE: https://cigale.lam.fr</footnote>

<footnote n="5">⁵ S17 library: http://cschreib.github.io/s17-irlib/</footnote>
with only one photometric point), are 0.49, 2.95, and 1.95 for OFG2, OFG20, and OFG27, respectively. The sample is obviously statistically limited but we do not find a systematic offset when using only one photometric point.

One caveat for the $M_{\text{dust}}$ estimates from CIGALE or the IR template library is that they are based on different dust models. Compared to the more standard dust models of Draine et al. (2014; an updated version of Draine & Li 2007) that we adopted in CIGALE, the one used in the IR template library of Schreiber et al. (2018c) assumes that the carbonated grains are amorphous carbon grains rather than graphites. Schreiber et al. (2018c) stated that different dust grain species from the IR template library have different emissivities, systematically lowering the derived $M_{\text{dust}}$ by a factor of about two. Therefore, to have comparable $M_{\text{dust}}$ for galaxies with and without a Herschel counterpart, we have corrected the differences in $M_{\text{dust}}$ obtained using the IR template library in Table 5 and also in the following sections.

### 4.3. SFRs

The total SFR was measured from the contributions of dust-obscured star formation (SFRIR) and unobscured star formation (SFRUV). The SFRIR was calculated based on the total infrared luminosity ($L_{\text{IR}}$), derived from integrating the best-fitted SED between 8 and 1000 µm in the rest frame, following Kennicutt & Evans (2012). The SFRUV was derived from the luminosity emitted in the UV ($L_{\text{UV}}$), which was not corrected for dust attenuation, following Daddi et al. (2004) (scaled to a Chabrier 2003 IMF). We calculated the total SFR:

$$ \begin{aligned} \text{SFR}_{\text{tot}} [M_{\odot} \text{yr}^{-1}] &= \text{SFR}_{\text{IR}} + \text{SFR}_{\text{UV}} \\ &= 1.49 \times 10^{-10} L_{\text{IR}} + 1.27 \times 10^{-10} L_{\text{UV}}, \end{aligned} \quad (2) $$

both $L_{\text{IR}}$ and $L_{\text{UV}}$ in units of $L_{\odot}$, and

$$ L_{\text{IR}} [L_{\odot}] = 4\pi D_L^2 \int_{8 \mu\text{m}}^{1000 \mu\text{m}} F_{\nu}(\lambda) \times \frac{c}{\lambda^2} d\lambda, \quad (3) $$

$$ L_{\text{UV}} [L_{\odot}] = 4\pi D_L^2 \frac{\nu_{1500}}{(1+z)} \frac{10^{-0.4(48.6+m_{1500})}}{3.826 \times 10^{33}}, \quad (4) $$

where $D_L$ is the luminosity distance (cm), $\nu_{1500}$ is the frequency (Hz) corresponding to the rest-frame wavelength 1500 Å, and $m_{1500}$ is the AB magnitude at the rest-frame 1500 Å. Here, the value for $m_{1500}$ was derived from the best-fitting templates using EAZY, with a top-hat filter centered at 1500 Å and a width of 350 Å.

For individual OFGs with a Herschel and/or ALMA counterpart (and therefore SFRIR), we present their SFR (SFRtot = SFRUV+SFRIR) in Table 2. Further discussion will be provided later in Sect. 5.3.1 for the stacked samples.

### 4.4. Molecular gas mass

The gas mass, $M_{\text{gas}}$, can be determined from $M_{\text{dust}}$ by employing the gas-to-dust ratio ($\delta_{\text{GDR}}$) with a metallicity dependency (e.g., Magdis et al. 2012):

$$ M_{\text{gas}} = M_{\text{H}_2} + M_{\text{H}_i} = \delta_{\text{GDR}} M_{\text{dust}}, \quad (5) $$

$$ \log(\delta_{\text{GDR}}) = (10.54 \pm 1.0) - (0.99 \pm 0.12) \times (12 + \log(\text{O/H})). \quad (6) $$

<table><caption>Table 5. Dust and gas properties of the OFGs.</caption><thead><tr><th>ID<br>(1)</th><th>log(<i>M</i><sub>dust</sub>)<br>log(<i>M</i><sub>☉</sub>)<br>(2)</th><th><i>T</i><sub>dust</sub><br>(K)<br>(3)</th><th>log(<i>M</i><sub>gas</sub>)<br>log(<i>M</i><sub>☉</sub>)<br>(4)</th><th><i>f</i><sub>gas</sub><br>(5)</th><th>τ<sub>dep</sub><br>(Myr)<br>(6)</th></tr></thead><tbody><tr><td>OFG1</td><td>8.49<sup>+0.05</sup><sub>-0.06</sub></td><td>(...)</td><td>10.64<sup>+0.05</sup><sub>-0.06</sub></td><td>0.42<sup>+0.08</sup><sub>-0.08</sub></td><td>115 ± 20</td></tr><tr><td>OFG2</td><td>8.80<sup>+0.10</sup><sub>-0.13</sub></td><td>37.0<sup>+1.7</sup><sub>-1.6</sub></td><td>10.96<sup>+0.10</sup><sub>-0.13</sub></td><td>0.61<sup>+0.09</sup><sub>-0.10</sub></td><td>249 ± 67</td></tr><tr><td>OFG4</td><td>8.04<sup>+0.07</sup><sub>-0.09</sub></td><td>(...)</td><td>10.39<sup>+0.07</sup><sub>-0.09</sub></td><td>0.38<sup>+0.11</sup><sub>-0.14</sub></td><td>46 ± 12</td></tr><tr><td>OFG5</td><td>8.09<sup>+0.06</sup><sub>-0.07</sub></td><td>(...)</td><td>10.49<sup>+0.06</sup><sub>-0.07</sub></td><td>0.67<sup>+0.26</sup><sub>-0.11</sub></td><td>114 ± 31</td></tr><tr><td>OFG7</td><td>8.37<sup>+0.04</sup><sub>-0.04</sub></td><td>(...)</td><td>10.75<sup>+0.04</sup><sub>-0.04</sub></td><td>0.72<sup>+0.05</sup><sub>-0.04</sub></td><td>52 ± 7</td></tr><tr><td>OFG11</td><td>8.28<sup>+0.06</sup><sub>-0.07</sub></td><td>(...)</td><td>10.57<sup>+0.06</sup><sub>-0.07</sub></td><td>0.59<sup>+0.05</sup><sub>-0.05</sub></td><td>97 ± 27</td></tr><tr><td>OFG18</td><td>7.96<sup>+0.11</sup><sub>-0.16</sub></td><td>(...)</td><td>10.71<sup>+0.11</sup><sub>-0.16</sub></td><td>0.95<sup>+0.02</sup><sub>-0.02</sub></td><td>180 ± 77</td></tr><tr><td>OFG19</td><td>8.02<sup>+0.06</sup><sub>-0.06</sub></td><td>(...)</td><td>10.39<sup>+0.06</sup><sub>-0.06</sub></td><td>0.55<sup>+0.05</sup><sub>-0.06</sub></td><td>90 ± 23</td></tr><tr><td>OFG20</td><td>8.40<sup>+0.04</sup><sub>-0.05</sub></td><td>68.5<sup>+4.2</sup><sub>-3.7</sub></td><td>10.59<sup>+0.04</sup><sub>-0.05</sub></td><td>0.34<sup>+0.10</sup><sub>-0.10</sub></td><td>43 ± 6</td></tr><tr><td>OFG25</td><td>8.66<sup>+0.03</sup><sub>-0.03</sub></td><td>(...)</td><td>11.11<sup>+0.03</sup><sub>-0.03</sub></td><td>0.93<sup>+0.12</sup><sub>-0.02</sub></td><td>105 ± 9</td></tr><tr><td>OFG27</td><td>9.01<sup>+0.04</sup><sub>-0.05</sub></td><td>51.4<sup>+2.2</sup><sub>-2.1</sub></td><td>11.07<sup>+0.04</sup><sub>-0.05</sub></td><td>0.48<sup>+0.11</sup><sub>-0.09</sub></td><td>66 ± 8</td></tr></tbody></table>

**Notes.** (1) Source ID; (2) Dust mass obtained from CIGALE for the galaxies with a Herschel counterpart or from the IR template library (Schreiber et al. 2018c) for the galaxies without a Herschel counterpart but with an ALMA counterpart (see Sect. 4.2). Since the dust emissivity used in Schreiber et al. (2018c) is different from the one of Draine et al. (2014; used in CIGALE), resulting in a systematic twice lower $M_{\text{dust}}$, here we have corrected the $M_{\text{dust}}$ derived from the IR template library by multiplying by two. (3) Dust temperature obtained from a single temperature MBB model for the galaxies with a Herschel counterpart (see Sect. 4.5); (4) Gas mass obtained from the metallicity-dependent gas-to-dust mass ratio technique (see Sect. 4.4); (5) Gas fraction: $f_{\text{gas}} = M_{\text{gas}}/(M_{\star}+M_{\text{gas}})$; (6) Gas depletion time: $\tau_{\text{dep}} = M_{\text{gas}}/(SFR_{\text{IR}}+SFR_{\text{UV}})$, which is the inverse of the star formation efficiency ($SFE = 1/\tau_{\text{dep}}$).

The metallicity was determined from the redshift-dependent mass-metallicity relation (MZR; Genzel et al. 2015):

$$ 12 + \log(\text{O/H}) = a - 0.087[\log(M_{\star}) - b]^2, \quad (7) $$

where $a = 8.74$ and $b = 10.4 + 4.46 \times \log(1+z) - 1.78 \times \log(1+z)^2$. We adopted an uncertainty of 0.2 dex in the metallicities (Magdis et al. 2012).

With the estimates of $M_{\text{gas}}$, we can calculate the gas fraction ($f_{\text{gas}}$) and gas depletion time ($\tau_{\text{dep}}$) as $f_{\text{gas}} = M_{\text{gas}}/(M_{\star} + M_{\text{gas}})$ and $\tau_{\text{dep}} = M_{\text{gas}}/(SFR_{\text{IR}} + SFR_{\text{UV}})$. The $\tau_{\text{dep}}$ is the inverse of the star formation efficiency ($SFE = 1/\tau_{\text{dep}}$). The $M_{\text{gas}}$, $f_{\text{gas}}$, and $\tau_{\text{dep}}$ for the individual sources are presented in Table 5. We underline that only three OFGs have a Herschel counterpart, and the rest have $M_{\text{dust}}$ only based on ALMA 1.13 mm. Thus, there is a large uncertainty in the values of $M_{\text{dust}}$ and, consequently, $M_{\text{gas}}$, $f_{\text{gas}}$, and $\tau_{\text{dep}}$. Therefore, this paper does not go deeper into the $M_{\text{dust}}$ and gas properties of individual galaxies. Instead, for the study of gas properties of the OFGs, we focus on the stacked sample, described in Sect. 5.4.

### 4.5. Dust temperatures

For the comparison with previous studies, we measured the effective dust temperatures ($T_{\text{dust}}$) by fitting single-temperature modified black-body (MBB) models to the FIR to mm photometry of the individual galaxies with a Herschel counterpart, following:

$$ S_{\nu} \propto \frac{\nu^{3+\beta}}{\exp\left(\frac{\hbar\nu}{k_{\text{B}}T_{\text{dust}}}\right) - 1}, \quad (8) $$

under the assumption of optically thin dust, where $S_{\nu}$ is the flux density, $k_{\text{B}}$ is Boltzmann constant, $\hbar$ is the Planck constant,
and $\beta$ is the dust emissivity index. We assumed $\beta = 1.5$, a typical value for dusty star-forming galaxies (e.g., Hildebrand 1983; Kovács et al. 2006; Gordon et al. 2010). We note that changing $\beta$ does not have a significant effect on $T_{\text{dust}}$, as $\beta$ is affecting the slope of the Rayleigh-Jeans (RJ) tail of the dust emission at the rest-frame $\lambda \ge 200$ µm, while the peak of the dust SED is what determines $T_{\text{dust}}$ (e.g., Casey 2012; Jin et al. 2019).

Following the criteria used in Hwang et al. (2010); also used in Franco et al. 2020b; Gómez-Guijarro et al. 2022b), we only fit the observed data points at $\lambda \ge 0.55\lambda_{\text{peak}}$ to avoid contamination from small dust grains, polycyclic aromatic hydrocarbon (PAH) molecules, and/or AGNs in the MIR, where $\lambda_{\text{peak}}$ is the peak of IR SEDs of the CIGALE best fit. The fitted galaxies should satisfy the following conditions: (i) at least one data point at $0.55\lambda_{\text{peak}} \le \lambda < \lambda_{\text{peak}}$; and (ii) at least one data point at $\lambda > \lambda_{\text{peak}}$ (to exclude the synchrotron contribution from radio data). In our case, we finally fitted the photometry from Herschel/SPIRE bands (250 µm, 350 µm, and 500 µm) and ALMA 1.13 mm. We note that we did not consider the CMB effect in the MBB fit because of the lack of data points in the RJ-tail where the CMB plays an important role (e.g., Jin et al. 2019).

The MBB fit was performed using a Markov Chain Monte Carlo (MCMC) approach with 12000 iterations using the Python package PyMC3⁶. The derived $T_{\text{dust}}$ for individual galaxies with a Herschel counterpart are listed in Table 5. For individual galaxies (three galaxies with Herschel counterparts were fitted with MBB), their $T_{\text{dust}}$ exhibit a large dispersion, that is, $T_{\text{dust}} = 37-69$ K, and only one ($T_{\text{dust}} = 51 \pm 2$ K) is in agreement with the expected value from the redshift evolution in MS galaxies from the literature (Schreiber et al. 2018c). The fitting results are also shown in Figs. C.1 and C.2. Considering that only three OFGs have $T_{\text{dust}}$ values, we do not go deeper into the discussion of individual galaxies. Instead, we further discuss the stacked OFG samples in Sect. 5.4.

# 5. Properties of the stacked OFGs

## 5.1. Stacking analysis

We performed a stacking analysis to study the global properties of our sample. Considering that OFGs are very faint in the optical/NIR ($H > 26.5$ mag) and only 11/27 have Herschel and/or ALMA counterparts, performing a stacking analysis helps improving the accuracy of the median photometric redshift and SFR measurements. To build the SED of our sample, we generated a median and mean stacked image in each filter, from the optical to 1.13 mm. Specifically, we used images from the HST/ACS (F435W, F606W, F775W, F814W, F850LP), HST/WFC3 (F105W, F125W, F140W, F160W), ZFOURGE Ks-band, Spitzer/IRAC (3.6, 4.5, 5.8, and 8 µm), Spitzer/MIPS (24 µm), Herschel/PACS&SPIRE (100, 160, 250, 350 and 500 µm), and ALMA 1.13 mm maps.

The photometry was obtained mostly using aperture photometry techniques, except for the Herschel bands, where appropriate aperture corrections were applied to account for flux losses outside the aperture. This procedure is very similar to that used previously in the deep surveys, which we summarize here. In the HST/ACS and HST/WFC3 bands, fluxes were extracted on the PSF-matched images (to the F160W) using the same aperture of 0.7"-diameter as in Whitaker et al. (2019), which maximizes the signal-to-noise ratio (S/N) of the resulting aperture photometry. In the Ks-band, we used a 1.2" diameter circular aperture

to measure flux on the ZFOURGE Ks-band image following Straatman et al. (2016), whose PSF was matched to a Moffat profile with $FWHM = 0.9"$. In the IRAC bands, fluxes were extracted separately without PSF matching due to the broader PSFs. We adopted a 2.2" diameter aperture to maximize the S/N of the resulting aperture photometry. In the MIPS 24 µm band, we used a large aperture of 6" in diameter corresponding to its full width at half maximum. At 1.13 mm, we used a diameter of 1.6" to measure the flux, which is the optimal trade-off between total flux and SNR.

Uncertainties on the photometry were derived from the Monte Carlo simulations. For each band, we carried out the same stacking analysis as above, but at random positions, and measured the flux value on the stacked image. This was repeated 1000 times. We then calculated the 16th and 84th percentiles of the distribution of values as flux uncertainties.

For the Herschel/PACS and SPIRE bands, we used the PSF fitting with a free background to fit the stacked image following Schreiber et al. (2015). The uncertainties were obtained using the following methods: 1) a bootstrap approach; specifically, as an example, we generated a sample of 27 sources from 27 OFGs, allowing the same galaxy to be picked repeatedly, and measured the stacked flux. This procedure was repeated 100 times, and we calculated their standard deviation as the flux uncertainty, and 2) a Monte Carlo simulation approach, which is the same as that used for lower wavelength images and 1.13 mm images. Here we adopted a $0.9 \times FWHM$ diameter circular aperture. We note that the results given from the bootstrap approach include the uncertainties from i) the PSF fitting, ii) the clustering bias effect, and iii) background fluctuation. Thus, the derived values of uncertainties from bootstrap are larger than those from the Monte Carlo simulation. We conservatively take the former values as our uncertainties.

## 5.2. Fitting of the stacked SEDs

We obtained the stacked full-wavelength SEDs in the same way as for the individual galaxies (see Sects. 4.1, 4.2, Tables 3, 4). In brief, first, we fitted the broad photometry at OPT to MIR with the EAZY code to obtain photometric redshifts. Then, we independently performed the OPT to MIR SED fitting with FAST++ and the MIR to mm SED fitting with CIGALE, respectively, at the previously obtained redshifts of the stacked sources. This approach helps to 1) disentangle the degeneracy between redshift and other parameters, such as stellar age and dust temperature; and 2) break the energy balance principle (the total energy emitted in the MIR and FIR is determined by the attenuation of observed starlight in the UV and optical) used in CIGALE. For dusty star-forming galaxies, especially for H-dropouts and Ks-dropouts with strong dust obscuration, there could exist regions with strong UV extinction due to strong dust obscuration, which may not participate in the UV to optical part, but emit FIR light (e.g., Simpson et al. 2015; Gómez-Guijarro et al. 2018; Elbaz et al. 2018). Assuming an energy balance with a fixed redshift will lead to an underestimation of the $L_{\text{IR}}$, hence, the SFR.

The best-fit SED is shown in Fig. 6. The median redshift for the total sample is $z_{\text{med,stacked}} = 4.5 \pm 0.2$, which is consistent with $z_{\text{med}} = 4.1$ derived from the median value of individual OFGs with a wide distribution (see Fig. 4, top). In addition, the median stacked SED peak (and the mean stacked SED peak; see Fig. 7) is between 350 and 500 µm, also in agreement with being at $z \sim 4$. Thus, these agreements confirm that the bulk population of OFGs consists of dusty star-forming galaxies at $z \sim 4-5$. Remarkably, most fluxes in the stacked images are above the $3\sigma$

<footnote n="6">⁶ PyMC3 is available at: https://docs.pymc.io/en/v3/</footnote>
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Median stacked SED and images of the total sample of 27 OFGs in this work. **Top:** best-fit SED of the total sample (black line). The measured fluxes (red points) are derived from the stacked images. Error bars ($1\sigma$) and upper limits ($3\sigma$) are obtained from the Monte Carlo simulation (except Herschel) and bootstrap approach (Herschel; see Sect. 5.2). We also show the best-fit SED for 22 non-LBGs (grey line). These 22 non-LBGs will be used to calculate the cosmic SFRD. The inset shows the likelihood distributions of the photometric redshift of our samples (total sample in red, 22 non-LBGs in grey), based on the UV to MIR SED fitting from EAzY, which is normalized to the peak value. The redshift obtained from the maximized likelihood is $z \sim 4.5$ for the total 27 OFGs and $z \sim 4.2$ for the 22 non-LBGs. **Bottom:** stacked images of the total sample with peak fluxes normalized. Each panel is $6'' \times 6''$ except for the MIPS 24 µm, which is $24'' \times 24''$.</figcaption></figure>

confidence level, especially in the H, Ks, and IRAC bands, helping to establish the position of the Balmer and 4000 Å breaks very well, hence determining a robust redshift. The median properties derived from the stacked SED for the total sample are summarized in Table 6.

In addition, to further investigate the characteristics of different subpopulations of our OFGs, we performed median and mean stacked SED fitting for four sub-samples of OFGs. The four sub-samples and the purpose of our investigation are: 1) OFGs that are LBGs: in our sample, five OFGs are classified as LBGs. Given that the traditional approach to estimating the cosmic SFRD at $z > 3$ is mainly based on the LBGs

(e.g., Madau & Dickinson 2014; Bouwens et al. 2020), studying the differences in the properties of LBGs and OFGs (after removing 5 LBGs) in our sample can help us understand the importance of OFGs in the cosmic SFRD; 2) OFGs that are not LBGs: 22 OFGs are not classified as LBGs. To quantify the level of underestimation of the cosmic SFRD at $z > 3$ (Madau & Dickinson 2014) based mainly on LBGs, we used the mean SFR of this sub-sample to calculate the SFRD to avoid contamination by LBGs (describe later in Sect. 6.1); 3) Massive OFGs that are not LBGs: there are 16 OFGs not classified as LBGs with $log(M_*/M_\odot) > 10.3$. To compare with the results of Wang et al. (2019) on the SFRD, here we
<table><caption>**Table 6.** Median physical properties of the total sample of 27 OFGs.</caption><thead><tr><th colspan="3">Derived from SED fitting with median stacked photometry <sup>(*)</sup></th></tr></thead><tbody><tr><td>z<sub>phot</sub> <sup>(a)</sup></td><td></td><td>4.5 ± 0.2</td></tr><tr><td>M<sub>⋆</sub> <sup>(b)</sup></td><td>M<sub>☉</sub></td><td>(2.8<sup>+0.2</sup><sub>-0.1</sub>) × 10<sup>10</sup></td></tr><tr><td>L<sub>IR</sub> <sup>(c)</sup></td><td>L<sub>☉</sub></td><td>(1.6 ± 0.3) × 10<sup>12</sup></td></tr><tr><td>A<sub>V</sub> <sup>(b)</sup></td><td>mag</td><td>0.9<sup>+0.3</sup><sub>-0.0</sub></td></tr><tr><td>M<sub>dust</sub> <sup>(c)</sup></td><td>M<sub>☉</sub></td><td>(1.2 ± 0.2) × 10<sup>8</sup></td></tr><tr><td>T<sub>dust</sub> <sup>(d)</sup></td><td>K</td><td>45.5<sup>+2.1</sup><sub>-2.1</sub></td></tr><tr><th colspan="3">Median stacked photometry</th></tr><tr><td>H</td><td>mag</td><td>27.4 ± 0.1</td></tr><tr><td>[4.5]</td><td>mag</td><td>23.92 ± 0.04</td></tr><tr><td>S<sub>1.13 mm</sub></td><td>μJy</td><td>334 ± 24</td></tr><tr><td>S<sub>3 GHz</sub></td><td>μJy</td><td>4.1 ± 0.7</td></tr><tr><th colspan="3">Derived quantities <sup>(**)</sup></th></tr><tr><td>L<sub>1.4 GHz</sub> <sup>(e)</sup></td><td>erg s<sup>-1</sup> Hz<sup>-1</sup></td><td>(9.53 ± 1.53) × 10<sup>30</sup></td></tr><tr><td>q<sub>TIR</sub> <sup>(f)</sup></td><td></td><td>2.23 ± 0.03</td></tr><tr><td>SFR<sub>rad,med</sub> <sup>(g)</sup></td><td>M<sub>☉</sub> yr<sup>-1</sup></td><td>287.59 ± 53.86</td></tr><tr><td>SFR<sub>IR,med</sub> <sup>(h)</sup></td><td>M<sub>☉</sub> yr<sup>-1</sup></td><td>235.33<sup>+47.37</sup><sub>-49.77</sub></td></tr><tr><td>SFR<sub>UV,med</sub> <sup>(i)</sup></td><td>M<sub>☉</sub> yr<sup>-1</sup></td><td>0.33 ± 0.02</td></tr><tr><td>ΔMS <sup>(j)</sup></td><td></td><td>1.45<sup>+0.29</sup><sub>-0.31</sub></td></tr><tr><td>M<sub>gas</sub> <sup>(k)</sup></td><td>M<sub>☉</sub></td><td>(2.6 ± 0.4) × 10<sup>10</sup></td></tr><tr><td>f<sub>gas</sub> <sup>(l)</sup></td><td></td><td>0.48 ± 0.05</td></tr><tr><td>τ<sub>dep</sub> <sup>(m)</sup></td><td>Myr</td><td>110<sup>+29</sup><sub>-30</sub></td></tr></tbody></table>

**Notes.** <sup>(*)</sup>Uncertainties are the 16–84th percentile ranges of the probability distribution function given by the SED fitting. <sup>(**)</sup>Uncertainties on derived quantities were calculated from the propagation of the errors in the parameter values. <sup>(a)</sup>Photometric redshift, determined with the code EAZY. <sup>(b)</sup>M<sub>⋆</sub> and A<sub>V</sub>, derived from the UV to MIR SED fitting with the code FAST+. <sup>(c)</sup>Given by IR SED fitting with CIGALE. <sup>(d)</sup>Measured by MBB model fit (see Sect. 4.5). <sup>(e)</sup>Derived from S<sub>3 GHz</sub> assuming a radio spectral index α = −0.75 (see Eq. (10)). <sup>(f)</sup>Calculated from the IR-radio correlation (see Eq. (12)). <sup>(g)</sup>Calculated following Delhaize et al. (2017; see Eq. (9)), which was simply estimated from the radio emission without correction for AGN. <sup>(h)</sup>Derived following Kennicutt & Evans (2012; see Eq. 2). <sup>(i)</sup>Derived following Daddi et al. (2004), scaled to a Chabrier (2003) IMF (see Eq. (2)). <sup>(j)</sup>Distance to the SFMS: ΔMS = SFR/SFR<sub>MS</sub>, where SFR<sub>MS</sub> is the average SFR of MS galaxies at fixed stellar mass and redshift (Schreiber et al. 2015, see Fig. 9 and Sect. 5.4). <sup>(k)</sup>M<sub>gas</sub>, computed based on gas-to-dust ratio. <sup>(l)</sup>Gas fraction: f<sub>gas</sub> = M<sub>gas</sub>/(M<sub>⋆</sub> + M<sub>gas</sub>). <sup>(m)</sup>Gas depletion time: τ<sub>dep</sub> = M<sub>gas</sub>/(SFR<sub>IR</sub> + SFR<sub>UV</sub>), which is the inverse of the star formation efficiency (SFE = 1/τ<sub>dep</sub>).

used the same stellar mass cut for this sub-sample; and 4) ALMA-detected OFGs: 11 OFGs in our sample are detected by ALMA at 1.13 mm (>3.5σ; Gómez-Guijarro et al. 2022a). Since using ALMA detections to select OF sources is a very efficient method (e.g., Franco et al. 2020a; Zhou et al. 2020; Gómez-Guijarro et al. 2022a), we study the properties of the ALMA-detected OFGs and compare them with other sub-samples to understand whether there is a selection bias using this approach.

The best-fit mean SEDs of the total sample and the four sub-samples are shown in Fig. 7. For the OFGs that are LBGs (Fig. 7b), there is a 3.3σ detection at 24 µm, no detection in all the Herschel bands, and a 1.2σ detection at 1.13 mm. To successfully perform the FIR SED fitting, we fit the fluxes at 24 µm and 1.13 mm. We then compared the best-fit model with the 3σ upper limits in the Herschel bands (red arrows in Fig. 7b). The best-fit model is below the red arrows, showing a good consistency with the Herschel no detections. The mean properties derived from

the mean stacked SEDs for the sub-samples are summarized in Table 7.

For the mean stacked SED fit, one hypothesis here is that all OFGs have a similar SED shape. This is because the stacked SED we used to derive the SFR is a flux-weighted average in each band and if the brightest galaxy has a different SED shape, then the fitting results will be biased towards the properties of the brightest galaxy. For example, it has been shown that the dust temperature derived from the mean stacked SED is biased by 1.5 K to a higher temperature than the true value since the starburst galaxies in the sample are warmer and brighter (Schreiber et al. 2018c). However, we do not yet know the true IR SED shapes of most OFGs in our sample because of the lack of the Herschel detections. It is also unclear whether the brightest OFGs have a different SED shape compared to the remaining OFGs, therefore causing a bias. Hence, we cannot correct this potential bias here. Instead, we performed a median stacked SED fitting as a comparison. Although the median stacked SED fitting exhibits a lower confidence level (because it is less influenced by the brightest sources) compared to the mean one, their properties are more robust against outliers and are representative of the vast majority of galaxies in the sample. On the other hand, and most importantly, we need the SFR derived from the mean stacked SED fitting to calculate the cosmic SFRD (described later in Sect. 6.1).

## 5.3. SFRs and AGN

### 5.3.1. SFRs

We obtained the SFR<sub>tot</sub> of the stacked OFG (sub)samples using the same method as for individual galaxies (see Sect. 4.3), following Eq. (2). With the 3 GHz VLA observations in the GOODS-South, we can also calculate the radio-based SFR (SFR<sub>rad</sub>; assuming a Chabrier 2003 IMF) following Delhaize et al. (2017):

$$ \text{SFR}_{\text{rad}} [\text{M}_{\odot} \text{yr}^{-1}] = 10^{-24} \times 10^{q_{\text{TIR}}(z)} L_{1.4 \text{ GHz}} [\text{W Hz}^{-1}], \quad (9) $$

where L<sub>1.4 GHz</sub> is the rest-frame 1.4 GHz luminosity converted from the 3 GHz flux density (S<sub>3 GHz</sub> at observed-frame; W m<sup>-2</sup> Hz<sup>-1</sup>) using:

$$ L_{1.4 \text{ GHz}} [\text{W Hz}^{-1}] = \frac{4\pi D_L^2}{(1+z)^{\alpha+1}} \left( \frac{1.4 \text{ [GHz]}}{3 \text{ [GHz]}} \right)^{\alpha} S_{3 \text{ GHz}}, \quad (10) $$

here, the radio spectral index α<sub>1.4 GHz</sub><sup>3 GHz</sup> is assumed to be α = −0.75. The q<sub>TIR</sub>(z) in Eq. (9) is the IR-to-radio luminosity ratio, which was recently found to evolve primarily with the stellar mass and depend secondarily on the redshift (Delvecchio et al. 2021):

$$ q_{\text{TIR}}(M_{\star}, z) = (2.646 \pm 0.024) \times (1 + z)^{-0.023 \pm 0.008} \\ - (0.148 \pm 0.013) \times (\log M_{\star}/M_{\odot} - 10). \quad (11) $$

The derived IR-based and radio-based SFR values (SFR<sub>IR</sub> and SFR<sub>rad</sub>) are in good agreement (except for the sub-sample of OFGs that are LBGs), as shown in Tables 6 and 7. For the sub-sample of OFGs that are LBGs, the SFR<sub>rad</sub> is about four times higher than the SFR<sub>IR</sub>, although with large uncertainty, hinting at the existence of radio AGNs in the sub-sample of OFGs that are LBGs. The median SFRs for our total 27 OFGs are given in Table 6, while the mean SFRs for our OF sub-samples are summarized in Table 7. For our entire sample, the median contribution from SFR<sub>UV</sub> to SFR<sub>tot</sub> is only 0.1%, which is negligible.
<figure><img src="image_8.png" /><figcaption>Fig. 7. Mean stacked SEDs for the total sample and four sub-samples of our OFGs. Panel a: total sample of 27 OFGs; (b)–(e): sub-samples of OFGs that are LBGs, OFGs that are not LBGs, massive OFGs (log(M*/M☉) > 10.3) that are not LBGs, and OFGs with ALMA detections, respectively. We use the sub-samples in panel c and panel d to calculate the SFRD. Top: best-fit SED of the stacked OFGs (black line), which is composed of the uncorrected dust-attenuated stellar component (dark blue line) and the IR dust contribution (orange solid line). The corrected UV emission is shown as a faint blue line. We also plot the best fit of a MBB (orange dashed line) with a fixed dust emissivity index β = 1.5. We note that there is no significant AGN contribution in our sample, except for the five OFGs that are LBGs displaying the presence of radio AGNs (see Sect. 5.3.2). The measured fluxes (red points) are derived from the stacked images. Error bars (1σ) and upper limits (3σ) are obtained from the Monte Carlo simulation (except Herschel) and bootstrap approach (Herschel; see Sect. 5.2). The inset shows the likelihood distribution of the photometric redshift based on the UV to MIR SED fitting from EAZY (see Sect. 4.1), which is normalized to the peak value. Bottom: residuals calculated from (observation – model)/observation. The SEDs for individual OFGs are given in the appendix (see Figs. C.1 and C.2).</figcaption></figure>

### 5.3.2. AGN

As our selection criterion has been designed to avoid selecting passive galaxies, the OFGs in our sample are mainly dusty star-forming galaxies. Studying the presence of AGN in our sample can help us understand the co-evolution between AGN and star formation activity in the early Universe. It is also crucial for ensuring that our calculations of the SFR and, eventually, the cosmic SFRD are correct (uncontaminated by AGN). Here, we examine our sample for AGN contributions using three different methods, that is, studying their IR, radio, and X-ray excesses.

First, we fit the stacked IR SED with an additional AGN template using CIGALE (Fritz et al. 2006, see Sect. 5.2 for the infrared SED fitting). The contribution of the IR-bright AGN to the total IR luminosity ($f_{AGN}$) can lead to an overestimation of the dust IR emission and thus of the total SFR. We found that the SED fitting yields a $f_{AGN} < 0.01$, indicating the absence of IR-bright AGN in our sample.

Secondly, the $q_{TIR}$ is defined as the IR-to-radio luminosity ratio (e.g., Helou et al. 1985; Yun et al. 2001):

$$q_{\text{TIR}} \equiv \log \left( \frac{L_{\text{IR}} [\text{W}]}{3.75 \times 10^{12} [\text{Hz}]} \right) - \log \left( L_{1.4 \text{ GHz}} [\text{W Hz}^{-1}] \right). \quad (12)$$

The derived $q_{TIR}$ from the IR-radio correlation (Eq. (12)) is presented in Tables 6 and 7. Except for the sub-sample of OFGs that are LBGs, $q_{TIR}$ values of the OFGs are consistent with those in Delvecchio et al. (2021) for star-forming galaxies at the same

redshift and stellar mass. These agreements suggest that for the OFGs that are not LBGs, there is a lack of strong AGN activity in the radio band. On the other hand, the sub-sample of OFGs that are LBGs present a mean $q_{TIR} = 1.8$, much smaller than the typical $q_{TIR}$ of 2.6 for star-forming galaxies at the same redshift and stellar mass, and would thus be classified as radio AGNs (Fig. 12 in Delvecchio et al. 2021).

We also searched for X-ray-bright AGN in the CDF-S 7 Ms catalog (Luo et al. 2017). Among 1008 sources in the main catalog and 47 lower-significance sources in the supplementary catalog, we did not find any X-ray counterpart for the individual OFGs within a 0.6" radius. None of the sources in our catalog exhibit a total X-ray luminosity integrated over the entire 0.5–7 keV range larger than $L_X = 10^{42.5} \text{ erg s}^{-1}$ (AGN definition in Luo et al. 2017). Hence, we find no evidence for any bright X-ray AGN in our catalog. We also performed mean and median stacking for the 27 OFGs in 0.5–7 keV images and did not find any significant detections (≪3σ) in either of the stacked images.

In addition, we considered the MIR-AGN selection criterion developed by Donley et al. (2012) to diagnose the presence of a power-law AGN based on IRAC colors. However, this criterion does not apply to our high-z OFGs. At $z > 3$, the IRAC bands mainly collect emission from stars below 2 µm in the rest frame, outside the typical domain where power-law AGNs contribute. In summary, we do not find evidence for significant contamination by AGNs in our OFG sample, except for the five OFGs that are LBGs displaying the presence of radio AGNs.
<table><caption>**Table 7.** Mean physical properties of sub-samples of OFGs.</caption><thead><tr><th>Parameter</th><th>Unit</th><th>OFGs that are LBGs<br>(5 Obj.)</th><th>OFGs that are not LBGs<br>(22 Obj.)</th><th>Massive OFGs that are not LBGs<br>(16 Obj.)</th><th>ALMA-detected OFGs<br>(11 Obj.)</th></tr></thead><tbody><tr><td colspan="6"><strong>Derived from SED fitting with mean stacked photometry</strong></td></tr><tr><td>z<sub>phot</sub></td><td></td><td>4.0 ± 0.3</td><td>4.3 ± 0.2</td><td>4.4 ± 0.2</td><td>4.2 ± 0.2</td></tr><tr><td>M<sub>★</sub></td><td>M<sub>☉</sub></td><td>(4.9<sup>+0.7</sup><sub>-3.8</sub>) × 10<sup>9</sup></td><td>(4.7 ± 0.3) × 10<sup>10</sup></td><td>(6.6<sup>+0.6</sup><sub>-0.4</sub>) × 10<sup>10</sup></td><td>(5.9<sup>+0.4</sup><sub>-0.3</sub>) × 10<sup>10</sup></td></tr><tr><td>L<sub>IR</sub></td><td>L<sub>☉</sub></td><td>(2.9 ± 2.1) × 10<sup>11</sup></td><td>(2.6<sup>+0.6</sup><sub>-0.5</sub>) × 10<sup>12</sup></td><td>(3.7 ± 0.6) × 10<sup>12</sup></td><td>(4.3<sup>+0.6</sup><sub>-0.7</sub>) × 10<sup>12</sup></td></tr><tr><td>A<sub>V</sub></td><td>mag</td><td>0.7<sup>+0.4</sup><sub>-0.3</sub></td><td>1.4<sup>+0.1</sup><sub>-0.3</sub></td><td>1.2<sup>+0.3</sup><sub>-0.1</sub></td><td>1.7<sup>+0.1</sup><sub>-0.2</sub></td></tr><tr><td>M<sub>dust</sub></td><td>M<sub>☉</sub></td><td>(0.3 ± 1.3) × 10<sup>8</sup></td><td>(2.9 ± 0.7) × 10<sup>8</sup></td><td>(2.8 ± 0.6) × 10<sup>8</sup></td><td>(5.3 ± 1.1) × 10<sup>8</sup></td></tr><tr><td>T<sub>dust</sub></td><td>K</td><td>(...)</td><td>42.3<sup>+1.2</sup><sub>-1.2</sub></td><td>45.5<sup>+1.4</sup><sub>-1.5</sub></td><td>41.5<sup>+1.2</sup><sub>-1.2</sub></td></tr><tr><td colspan="6"><strong>Mean stacked photometry</strong></td></tr><tr><td>H</td><td>mag</td><td>26.95 ± 0.17</td><td>27.09 ± 0.08</td><td>26.99 ± 0.08</td><td>26.85 ± 0.10</td></tr><tr><td>[4.5]</td><td>mag</td><td>24.74 ± 0.31</td><td>23.34 ± 0.03</td><td>23.33 ± 0.02</td><td>22.99 ± 0.03</td></tr><tr><td>S<sub>1.13 mm</sub></td><td>μJy</td><td>65 ± 56</td><td>650 ± 26</td><td>660 ± 30</td><td>1153 ± 40</td></tr><tr><td>S<sub>3 GHz</sub></td><td>μJy</td><td>2.7 ± 0.5</td><td>6.5 ± 0.5</td><td>6.9 ± 0.5</td><td>10.6 ± 0.8</td></tr><tr><td colspan="6"><strong>Derived quantities</strong></td></tr><tr><td>L<sub>1.4 GHz</sub></td><td>erg s<sup>-1</sup> Hz<sup>-1</sup></td><td>(4.81 ± 0.92) × 10<sup>30</sup></td><td>(1.40 ± 0.11) × 10<sup>31</sup></td><td>(1.54 ± 0.11) × 10<sup>31</sup></td><td>(2.18 ± 0.17) × 10<sup>31</sup></td></tr><tr><td>q<sub>TIR</sub></td><td></td><td>1.80<sup>+0.15</sup><sub>-0.43</sub></td><td>2.28 ± 0.03</td><td>2.39 ± 0.02</td><td>2.30 ± 0.03</td></tr><tr><td>SFR<sub>rad,avg</sub></td><td>M<sub>☉</sub> yr<sup>-1</sup></td><td>189.95 ± 58.40</td><td>392.17 ± 49.02</td><td>408.91 ± 49.54</td><td>590.15 ± 74.07</td></tr><tr><td>SFR<sub>IR,avg</sub></td><td>M<sub>☉</sub> yr<sup>-1</sup></td><td>44.52 ± 33.08</td><td>387.93<sup>+82.03</sup><sub>-72.08</sub></td><td>556.32<sup>+84.94</sup><sub>-85.59</sub></td><td>635.15<sup>+96.61</sup><sub>-104.83</sub></td></tr><tr><td>SFR<sub>UV,avg</sub></td><td>M<sub>☉</sub> yr<sup>-1</sup></td><td>0.63 ± 0.07</td><td>0.37 ± 0.03</td><td>0.45 ± 0.03</td><td>0.59 ± 0.04</td></tr><tr><td>ΔMS</td><td></td><td>1.81 ± 1.32</td><td>1.46<sup>+0.31</sup><sub>-0.27</sub></td><td>1.45 ± 0.22</td><td>1.96<sup>+0.30</sup><sub>-0.32</sub></td></tr><tr><td>M<sub>gas</sub></td><td>M<sub>☉</sub></td><td>(1.3 ± 5.7) × 10<sup>10</sup></td><td>(5.0 ± 1.2) × 10<sup>10</sup></td><td>(4.4 ± 0.9) × 10<sup>10</sup></td><td>(8.4 ± 1.8) × 10<sup>10</sup></td></tr><tr><td>f<sub>gas</sub></td><td></td><td>0.73 ± 0.87</td><td>0.52 ± 0.06</td><td>0.40 ± 0.06</td><td>0.59 ± 0.05</td></tr><tr><td>τ<sub>dep</sub></td><td>Myr</td><td>291 ± 1277</td><td>130<sup>+42</sup><sub>-40</sub></td><td>79 ± 21</td><td>133<sup>+34</sup><sub>-35</sub></td></tr><tr><td>R<sub>e(1.13 mm)</sub></td><td>kpc</td><td>(...)</td><td>1.09 ± 0.05</td><td>1.05 ± 0.06</td><td>0.80 ± 0.03</td></tr><tr><td>Σ<sub>SFR</sub><sup>(a)</sup></td><td>M<sub>☉</sub> yr<sup>-1</sup> kpc<sup>-2</sup></td><td>(...)</td><td>52 ± 11</td><td>80 ± 15</td><td>158 ± 28</td></tr><tr><td colspan="6"><strong>Cosmic SFRD<sup>(a)</sup></strong></td></tr><tr><td>V<sup>(a)</sup></td><td>Mpc<sup>3</sup></td><td>(...)</td><td>7.4 × 10<sup>5</sup></td><td>7.4 × 10<sup>5</sup></td><td>7.4 × 10<sup>5</sup></td></tr><tr><td>SFRD</td><td>M<sub>☉</sub> yr<sup>-1</sup>Mpc<sup>-3</sup></td><td>(...)</td><td>(1.2 ± 0.2) × 10<sup>-2</sup></td><td>(1.2 ± 0.2) × 10<sup>-2</sup></td><td>(0.9 ± 0.2) × 10<sup>-2</sup></td></tr></tbody></table>

Notes. Same as Table 6 but for the four sub-samples of OFGs. (a)SFR surface density: $\Sigma_{SFR} = 0.5SFR_{tot}/(\pi R_{e(1.13 mm)}^2)$, where $SFR_{tot} = SFR_{IR,avg} + SFR_{UV,avg}$. (b)The SFRD is discussed with details in Sect. 6. (c)Survey volume, calculated using Eq. (16) with a broad redshift range of $z = 3.2-7.0$.

## 5.4. The main sequence of star-forming galaxies and the properties of gas and dust for the stacked OFGs

In this section, we investigate the properties of the stacked samples derived from the SED fitting. We examine their locations in the star-formation main sequence, their gas depletion timescales, gas fractions, and dust temperatures in the framework of the scaling relations for galaxy evolution, and their dust sizes.

In Fig. 8, we place the stacked total OFG sample and the four sub-samples in the SFR-$M_★$ plane, showing the locations compared to the SFMS. In the SFR-$M_★$ plane, it is well known that the SFMS as a whole moves to higher SFRs with increasing redshift (e.g., Elbaz et al. 2007, 2011; Noeske et al. 2007; Magdis et al. 2010; Whitaker et al. 2012, 2014; Speagle et al. 2014; Schreiber et al. 2015; Lee et al. 2015; Leslie et al. 2020). We adopted a fixed $z = 4$ for the SFMS (Schreiber et al. 2015) as a comparison since it is close to $z_{med} = 4.1$ from the median value of the individual OFGs and $z_{med,stacked} = 4.5$ from the median stacked total OFGs. This figure shows that all the (sub-)samples are located within the SFMS region (~0.5 dex) at $z = 4$, and most of them lie within the $1\sigma$ scatter of the SFMS ($0.5 < \Delta\text{MS} < 2$, i.e., ~0.3 dex), consistent with being normal star-forming galaxies at the same redshift. It suggests that unlike studies limited to a rare population of extreme starburst galaxies (e.g., Riechers et al. 2013; Strandet et al. 2017; Marrone et al. 2018; Dudzevičiūtė et al. 2020; Riechers et al. 2020), our OFGs represent a less extreme population of dusty star-forming galaxies at $z > 3$.

Furthermore, we study the gas and dust properties of the OFGs by focusing on their gas depletion timescales, gas frac-

tions, dust temperatures, and dust sizes. In Fig. 9, we show the normalized $\tau_{\text{dep}}$, $f_{\text{gas}}$, and $T_{\text{dust}}$ by scaling them to the observed relation (scl; which is the median of the MS) of $\tau_{\text{dep,scl}}(z, M_★, \Delta\text{MS})$ and $f_{\text{gas,scl}}(z, M_★, \Delta\text{MS})$ from Tacconi et al. (2018) and of $T_{\text{dust,scl}}(z, \Delta\text{MS})$ from Schreiber et al. (2018c) as a function of $\Delta\text{MS}$ and $M_★$. The $\Delta\text{MS}$ is the SFR of each stacked (sub)sample normalized by the SFR of the SFMS ($\Delta\text{MS} = \text{SFR}/\text{SFR}_{\text{MS}}$; Schreiber et al. 2015) at its own redshift and stellar mass. The $\tau_{\text{dep,scl}}$, $f_{\text{gas,scl}}$, and $T_{\text{dust,scl}}$ are calculated for each data point at fixed redshift, stellar mass, and $\Delta\text{MS}$. We also present dust continuum sizes at 1.13 mm ($R_{e(1.13 mm)}$) of the mean stacked OFG (sub)samples as a function of $\Delta\text{MS}$ and $M_★$. Here, the half-light radius $R_{e(1.13 mm)}$ was measured in the uv plane by fitting a circular Gaussian (task uvmodelit) after performing uv plane stacking according to the method described by Gómez-Guijarro et al. (2022a). We note that we did not scale $R_{e(1.13 mm)}$ to the observed relations (e.g., van der Wel et al. 2014) because the redshifts ($z > 3$) of our OFGs exceed the limits of these relations.

In Fig. 9, there is no global offset between mean and median results of normalized $\tau_{\text{dep}}$, $f_{\text{gas}}$, and $T_{\text{dust}}$ for the stacked (sub)samples. It indicates no significant differences in the SED shapes of the brightest OFGs compared to the remaining ones, which would otherwise cause a strong bias (as discussed in Sect. 5.2) and further show a global offset even for the different stacked samples. The mean SFR is larger than the median SFR (see Fig. 8), which is expected since the former is influenced by the brightest sources in the flux-weighted average in each band (as discussed in Sect. 5.2).
<figure><img src="image_9.png" /><figcaption>Fig. 8. Locations of OFGs compared to the SFMS in the SFR-$M_{\star}$ plane. The SFMS at $z=4$ (Schreiber et al. 2015), $1\sigma$ scatter ($0.5 < \Delta\text{MS} < 2$, i.e., ~0.3 dex), and $\pm 3 \times \Delta\text{MS}$ region ($0.33 < \Delta\text{MS} < 3$, i.e., ~0.5 dex) are highlighted with a grey dashed line, a grey shaded area, and a light grey shaded area, respectively. $\Delta\text{MS} > 3$ is commonly used to separate MS and SB galaxies. Squares and stars respectively represent the median and mean stacking results for our total sample (black) and four sub-samples. The four sub-samples are the OFGs that are LBGs (blue), the OFGs that are not LBGs(green), the massive OFGs that are not LBGs(red), and the OFGs with ALMA detections (orange). When necessary, data from the literature have been converted to a Chabrier (2003) IMF.</figcaption></figure>

In the first and second columns of Fig. 9, the stacked (sub)samples show $\tau_{\text{dep}}$ values below the scatter of the scaling relation, while $f_{\text{gas}}$ is at the lower boundary of the scaling relation. That is to say, the OFGs have shorter $\tau_{\text{dep}}$ and slightly lower $f_{\text{gas}}$ values compared to normal star-forming galaxies. This indicates that galaxies with stronger dust obscuration tend to have lower gas fractions and shorter gas depletion times. Their gas is consumed more rapidly, hence, they form their stars with a high efficiency, which sets them in the so-called class of starbursts in the main sequence (Elbaz et al. 2018; Gómez-Guijarro et al. 2022b).

Among all the stacked (sub)samples, the ALMA-detected OFGs have the longest gas depletion timescale and the highest gas fraction. We believe this is due to a selection effect, as galaxies with higher dust content are more easily detected by ALMA at 1.13 mm. The $M_{\text{gas}}$ was derived from $M_{\text{dust}}$ in our study by employing the gas-to-dust ratio (see Sect. 4.4). Thus, the ALMA-detected galaxies tend to have higher $M_{\text{gas}}$ and, consequently, higher values of $\tau_{\text{dep}}$ and $f_{\text{gas}}$ as well. Furthermore, the SFR is positively correlated with $M_{\text{dust}}$ for star-forming galaxies at fixed $T_{\text{dust}}$ (e.g., Genzel et al. 2015; Orellana et al. 2017; Donevski et al. 2020). It explains why they show a higher SFR in the SFR-$M_{\star}$ plane compared to the total stacked OFGs (in Fig. 8). Notably, it raises the caveat that the approach of selecting only ALMA-detected galaxies in studies of OFGs will end up biasing the sample toward larger SFRs, longer $\tau_{\text{dep}}$, and larger $f_{\text{gas}}$.

In addition, the massive OFGs (excluding LBGs) present the lowest gas fraction and the shortest gas depletion timescale of all stacked (sub)samples in Fig. 9. Yet, we did not see any significant difference in the $\Delta\text{MS}$ of the massive OFGs compared with

the other stacked (sub)samples. This suggests that these galaxies are observed just before becoming passive.

The median $T_{\text{dust}} = 46 \pm 2$ K for the stacked total OFGs (see Table 6) is consistent with the scaling relation of $T_{\text{dust}}(z, \Delta\text{MS})$ (Schreiber et al. 2018c, black squares in third column of Fig. 9). However, surprisingly, most of the stacked (sub-)samples show slightly colder dust temperatures compared to the scaling relation. In particular, the ALMA-detected OFGs have the most abundant dust but show the lowest $T_{\text{dust}}$, indicating that the dust is colder in the more obscured sources ($A_{\text{V}} = 1.7_{-0.2}^{+0.1}$ in Table 7). The mean $T_{\text{dust}} = 42 \pm 1$ K of the ALMA-detected OFGs (see Table 7) is consistent with the $T_{\text{dust}} = 40 \pm 2$ K<fnref n="7" /> of the median stacked ALMA-detected massive H-dropouts at $z=4$ (Wang et al. 2019). The median $T_{\text{dust}}$ of the ALMA-detected OFGs is much lower, with $T_{\text{dust}} = 33 \pm 3$ K. The low $T_{\text{dust}}$ of the ALMA-detected OFGs cannot be explained by current studies (see, e.g., Magnelli et al. 2014; Schreiber et al. 2018c), which suggest that an increasing $T_{\text{dust}}$ is correlated to an enhanced specific star formation rate. Furthermore, this is contrary to the findings of Sommovigo et al. (2022), for instance, where the authors conclude that dust is warmer in obscured sources because a larger obscuration leads to more efficient dust heating. However, quite intriguingly, cases of cold dusty star-forming galaxies at high redshifts have already been reported in the literature, such as GN20 at $z=4.05$ with $T_{\text{dust}} = 33$ K (Magdis et al. 2012; Cortzen et al. 2020) and four ALMA-detected sources at 3 mm at $z \sim 5$ (Jin et al. 2019). A possible reason for the cold dust temperature is that the dust emission in the FIR of the dust-obscured sources may be optically thick rather than optically thin, where a warm and compact dust core is hidden (Jin et al. 2019, 2022). Indeed, the compact dust core is shown in the last column of Fig. 9. Among (sub-)samples of our OFGs, the ALMA-detected OFGs with the highest dust obscuration (largest $A_{\text{V}}$) present the most compact dust core with a half-light radius $R_{e(1.13\text{mm})} = 0.80 \pm 0.03$ kpc. The SFR surface density ($\Sigma_{\text{SFR}}$) of the ALMA-detected OFGs is about two to three times higher than the others. This would imply that the measured dust temperature underestimates the actual dust temperature, that would be higher after correcting for the attenuation in the shorter FIR bands. Making this correction is out of the scope of this paper due to the limited information that we have on those galaxies.

### 5.5. The hidden side of the dust region

An important check is to test whether UV continuum emission alone (after correcting for dust extinction) provides a robust estimate of the total SFR, especially for those highly dust-obscured galaxies. We again used the stacked total sample and four sub-samples of OFGs.

We derived the $SFR_{\text{UV}}$ corrected for dust extinction (i.e., $SFR_{\text{UV}}^{\text{corr}}$) using the Calzetti et al. (2000) reddening law and assuming a constant star formation history from the UV to MIR SED fitting with the code FAST++. Specifically, similar to Sect. 4.3, the $SFR_{\text{UV}}^{\text{corr}}$ was obtained from the $L_{\text{UV}}^{\text{corr}}$ following Daddi et al. (2004), which was calculated based on the AB magnitude at the rest-frame 1500 Å (see Eq. (4)). The intrinsic flux was derived with

$$f_{\text{int}}(\lambda) = f_{\text{obs}}(\lambda) 10^{0.4 A_{\lambda}}, \qquad (13)$$

<footnote n="7">⁷ $T_{\text{dust}}$ is derived using the IR template library (Schreiber et al. 2018c). To compare with our results, it has been scaled to the light-weighted dust temperature by applying Eq. (6) in Schreiber et al. (2018c).</footnote>
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Normalized $\tau_{\text{dep}}$ (left), $f_{\text{gas}}$ (mid-left), and $T_{\text{dust}}$ (mid-right) by the scaling relation of $\tau_{\text{dep,scl}}(z, M_{\star}, \Delta\text{MS})$ and $f_{\text{gas,scl}}(z, M_{\star}, \Delta\text{MS})$ from Tacconi et al. (2018) and of $T_{\text{dust,scl}}(z, \Delta\text{MS})$ from Schreiber et al. (2018c) as well as the half-light radius $R_e$ (right) at ALMA 1.13 mm as a function of $\Delta\text{MS}$ (first row) and $M_{\star}$ (second row). Here, $\Delta\text{MS}$ is the distance to the SFMS (Schreiber et al. 2015), $\Delta\text{MS} = \text{SFR}/\text{SFR}_{\text{MS}}$, at its own stellar mass and redshift for each sample. The blue dashed line and shaded area show the scaling relation and scatter. For each data point, the $\tau_{\text{dep,scl}}$, $f_{\text{gas,scl}}$, and $T_{\text{dust,scl}}$ are calculated at a fixed redshift, stellar mass, and $\Delta\text{MS}$. Squares, stars, and shaded grey areas are the same as in Fig. 8.</figcaption></figure>

where $f_{\text{int}}$ and $f_{\text{obs}}$ are the intrinsic and observed fluxes, respectively. The extinction $A_{\lambda}$ is related to the reddening curve $k(\lambda)$:

$$A_{\lambda} = \frac{k(\lambda) A_{V}}{R_{V}}. \qquad (14)$$

From the Calzetti et al. (2000) reddening law, we have $R_V = 4.05$ and

$$K_{\lambda} = 2.659 \left( -2.156 + \frac{1.509}{\lambda} - \frac{0.198}{\lambda^2} + \frac{0.011}{\lambda^3} \right) + R_V, \quad (15)$$

with $0.12 \mu\text{m} \le \lambda \le 0.63 \mu\text{m}$.

We then compared the derived $SFR_{UV}^{\text{corr}}$ with the total SFR ($SFR_{\text{tot}} = SFR_{\text{IR}} + SFR_{\text{UV}}$) and calculated their ratios. The results are presented in Table 8. We find that all the stacked (sub-)samples have $SFR_{\text{tot}}$ much larger than $SFR_{UV}^{\text{corr}}$. In addition, they have similar $SFR_{\text{tot}}/SFR_{UV}^{\text{corr}}$ ratios within uncertainties, with $SFR_{\text{tot}}/SFR_{UV}^{\text{corr}} = 8 \pm 1$ for the mean stacked total OFGs. It suggests: (i) the existence of a hidden dust region in the OFGs (even for the LBGs) that absorbs all the UV photons, which cannot be reproduced with a dust extinction correction, indicating that the dust emission in these OFGs might be optically thick; and (ii) that it is fundamental to include IR/mm band observations when studying extremely dusty star-forming galaxies. Otherwise, the total SFR and, therefore, the cosmic SFRD will be strongly underestimated.

Furthermore, we compared the stacked OFG (sub)samples with the star-forming galaxies from the ZFOURGE catalog (Straatman et al. 2016). These star-forming galaxies are from the CDF-S field with a Herschel detection, and are split into two redshift bins ($0 \le z \le 3$ and $3 \le z \le 6$). We calculated their $SFR_{\text{tot}}/SFR_{UV}^{\text{corr}}$ ratios using the same method as for the stacked OFG (sub)samples. As shown in Fig. 10, the $SFR_{\text{tot}}/SFR_{UV}^{\text{corr}}$ ratio increases with increasing starburstiness, indicating the presence

<table><caption>**Table 8.** UV-corrected SFR vs. total SFR for the mean stacked OFGs.</caption><thead><tr><th>Mean stacked OFGs</th><th>SFR<sub>UV</sub><sup>corr</sup><br>(M<sub>☉</sub> yr<sup>-1</sup>)</th><th>SFR<sub>tot</sub><br>(M<sub>☉</sub> yr<sup>-1</sup>)</th><th>SFR<sub>tot</sub>/SFR<sub>UV</sub><sup>corr</sup></th></tr><tr><th>(1)</th><th>(2)</th><th>(3)</th><th>(4)</th></tr></thead><tbody><tr><td>Total OFGs</td><td>46<sup>+3</sup><sub>-4</sub></td><td>348 ± 59</td><td>8 ± 1</td></tr><tr><td>OFGs that are LBGs</td><td>5<sup>+1</sup><sub>-1</sub></td><td>45 ± 33</td><td>9 ± 7</td></tr><tr><td>OFGs that are not LBGs</td><td>57<sup>+6</sup><sub>-4</sub></td><td>388 ± 82</td><td>7 ± 2</td></tr><tr><td>Massive OFGs that are not LBGs</td><td>88<sup>+9</sup><sub>-8</sub></td><td>557 ± 86</td><td>6 ± 1</td></tr><tr><td>ALMA-detected OFGs</td><td>83<sup>+10</sup><sub>-7</sub></td><td>636 ± 105</td><td>8 ± 2</td></tr></tbody></table>

**Notes.** (1) Mean stacked total sample and four sub-samples of OFGs; (2) $SFR_{UV}$ corrected for dust extinction (see Sect. 5.5); (3) $SFR_{tot} = SFR_{IR} + SFR_{UV}$; (4) Ratio of total SFR to $SFR_{UV}^{\text{corr}}$.

of more hidden dust regions in galaxies that are likely to be optically thick. It suggests that using the UV emission alone to determine the total SFR of starburst galaxies, even after dust attenuation correction, could result in strong underestimates, consistent with the findings of Elbaz et al. (2018) and Puglisi et al. (2017). We further find that the strong underestimations appear at both redshift bins, suggesting that this may be a general phenomenon for starburst galaxies, regardless of the redshift. In addition, for MS galaxies ($\Delta\text{MS} \sim 1$) with $0 \le z \le 3$, their $SFR_{tot}$ and $SFR_{UV}^{\text{corr}}$ are very similar, showing that both SFR estimators agree with each other for typical MS galaxies at low redshifts. However, for MS galaxies with $3 \le z \le 6$, their $SFR_{tot}$ is about twice (~0.3 dex) larger than the $SFR_{UV}^{\text{corr}}$. Generally, the median $SFR_{tot}/SFR_{UV}^{\text{corr}}$ ratio is about twice higher for star-forming galaxies with $3 \le z \le 6$ than those with $0 \le z \le 3$. We note that in Fig. 10, we did not perform a stellar-mass cut for the star-forming galaxies from the ZFOURGE catalog due to the small number of galaxies. As a test, we selected galaxies with $\log(M_{\star}/M_{\odot}) > 9$ and obtained similar results but with a larger dispersion of the galaxy distribution because of their
<figure><img src="image_11.png" /><figcaption>Fig. 10. Ratio of the total SFR (SFR<sub>tot</sub> = SFR<sub>IR</sub> + SFR<sub>UV</sub>) and the SFR<sub>UV</sub> corrected for the dust extinction (SFR<sub>UV</sub><sup>corr</sup>) as a function of the starburstiness (ΔMS = SFR<sub>tot</sub>/SFR<sub>MS</sub>). Stars represent the mean stacked OFG (sub)samples as in Fig. 8. The solid and dashed lines show the sliding median and 16–84th percentile range of star-forming galaxies from the CDF-S field with a Herschel detection at 0 ≤ z ≤ 3 (blue) and 3 ≤ z ≤ 6 (orange) from the ZFOURGE catalog (Straatman et al. 2016).</figcaption></figure>

small number. A more detailed study of the stellar masses, and the SFR<sub>tot</sub>/SFR<sub>UV</sub><sup>corr</sup> ratio is beyond the scope of this paper.

In Fig. 10, the stacked OFG (sub)samples, with z<sub>med,stacked</sub> = 4.5, lie above the 16–84th percentile range of the star-forming galaxies at 3 ≤ z ≤ 6. This is consistent with the fact that these dusty star-forming galaxies are the more extreme cases (more dust-obscured), with lower dust temperatures compared to typical star-forming galaxies at similar redshifts and with relatively compact dust sizes (R<sub>e</sub>(1.13 mm) = 1.01 ± 0.05 kpc; black star in Fig. 9).

# 6. Cosmic star formation rate density

## 6.1. Star formation rate density

In this section, we calculate the SFRD in a very conservative way. In our calculation, the derived value would be a lower limit for the OFGs.

As our OFGs were discovered randomly within a blind GOODS-ALMA survey area, we can simply calculate their SFRD by using their total SFR divided by the survey volume. The survey volume is the volume between the shells defined by the redshift range of the sources and within a solid angle:

$$V = \frac{\Omega}{4\pi} (V_{z1} - V_{z0}) = \frac{\Omega}{4\pi} \left[ \frac{4\pi}{3} (d_{z1}^3 - d_{z0}^3) \right] = \frac{\Omega}{3} (d_{z1}^3 - d_{z0}^3), \quad (16)$$

where the solid angle Ω corresponding to the effective area of 72.42 arcmin<sup>2</sup>, in units of steradian (6.1 × 10<sup>-6</sup> sr), of the GOODS-ALMA survey; d<sub>z0</sub> and d<sub>z1</sub> are the comoving distances at given redshifts of z<sub>0</sub> and z<sub>1</sub>. To be conservative, here we use a broad redshift range of z = 3.2–7.0 from the whole OFG sample (see Fig. 4) for our total 22 OFGs that are not LBGs, instead of using the redshift and 1σ confidence interval given by the stacked SED.

The cosmic star-formation-rate density for our total OFG sample is:

$$\psi [M_{\odot} \text{yr}^{-1} \text{Mpc}^{-3}] = \frac{\text{SFR}_{\text{mean}} \times 22}{V}, \quad (17)$$

where the SFR<sub>mean</sub> is the UV+IR-based SFR given from the mean stacked SED fitting of the 22 OFGs that are not LBGs (Table 7). Assuming that the intrinsic infrared SED of our OFGs is the same as the SED derived from mean stacking, the SFRD of our 22 OFGs (excluding LBGs) reaches approximately 1.2 × 10<sup>-2</sup> M<sub>☉</sub> yr<sup>-1</sup> Mpc<sup>-3</sup> with the mean redshift of z<sub>mean</sub> = 4.3 (and z<sub>median</sub> = 4.1). The survey volume, V, is 7.4 × 10<sup>5</sup> Mpc<sup>3</sup>. The number density is n ~ 3 × 10<sup>-5</sup> Mpc<sup>-3</sup>, which is slightly higher than the one of massive H-dropouts.

## 6.2. Incompleteness of our understanding of the cosmic star-formation history

The reference cosmic star-formation history at z ≥ 4 (Madau & Dickinson 2014) is based on measurements of the LBGs corrected for dust attenuation in the rest-frame UV (Bouwens et al. 2012a,b). The study of OFGs (excluding LBGs) can help us quantify the incompleteness of our understanding of the cosmic star-formation history at high redshifts.

The left panel of Fig. 11 shows the contribution of OFGs (z<sub>mean</sub> = 4.3) to the cosmic SFRD. At z = 4–5, the contribution of OFGs (red-filled circles) reaches about 52% of the reference SFRD from LBGs (Madau & Dickinson 2014, black solid line with grey open triangles). If we combine the contributions of OFGs (excluding LBGs) and LBGs (Madau & Dickinson 2014) as total SFRD, the OFGs contribute about 34% and the LBGs contribute about 66% of the total SFRD. It shows that the LBGs dominate the total SFRD, while the OFGs (excluding LBGs) make up about a third of the total SFRD.

In addition, we investigate the contribution of massive OFGs to the cosmic SFRD. Out of a total of 22 OFGs (excluding LBGs), 16 massive OFGs have log(M<sub>⋆</sub>/M<sub>☉</sub>) > 10.3, the same stellar mass cut as the so-called ALMA-detected massive H-dropouts (Wang et al. 2019). The SFRD of massive OFGs is approximately 1.2 × 10<sup>-2</sup> M<sub>☉</sub> yr<sup>-1</sup> Mpc<sup>-3</sup> with z<sub>mean</sub> = 4.4 (and z<sub>median</sub> = 4.1), shown as the blue-shaded area in the left panel of Fig. 11 (see also Table 7). If we limit the redshift range to z = 4–5, the SFRD of massive OFGs is approximately 1.6 × 10<sup>-2</sup> M<sub>☉</sub> yr<sup>-1</sup> Mpc<sup>-3</sup>, as illustrated by the blue-filled circles. This value is consistent with the SFRD of the total 22 OFGs within errors (red-filled circles in Fig. 11). Since the massive and total OFGs have similar SFRD values, it indicates that among the OFGs, the massive ones dominant the SFRD compared to the remaining less massive ones.

We further compare the massive OFGs with equivalently massive LBGs and H-dropouts. Firstly, comparing the massive OFGs with equivalently massive LBGs, the SFRD of massive OFGs is at least two orders of magnitude higher. Thus, we conclude that for massive galaxies (log(M<sub>⋆</sub>/M<sub>☉</sub>) > 10.3), the cosmic SFRD is dominated by massive OFGs (99%) rather than massive LBGs (1%) at high redshifts. Secondly, the SFRD of massive OFGs is more than four times higher than that of massive H-dropouts detected by ALMA (Wang et al. 2019). We recall that the selection criteria for H-dropouts are H > 27 mag & [4.5] < 24 mag (Wang et al. 2019), whereas our OFGs are selected with H > 26.5 mag & [4.5] < 25 mag. It implies that the optically faint galaxies that contribute significantly to the SFRD of massive galaxies have been neglected in previous
<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Cosmic star formation history of the Universe. **Left:** contribution of our OFGs to the cosmic SFRD. The black line is the cosmic SFRD, $\psi$, as a function of redshift from Madau & Dickinson (2014), which is based on LBGs at $z \ge 4$ with a dust correction applied (grey open triangles; Bouwens et al. 2012a,b). The 22 OFGs that are not LBGs are shaded in red and the 16 massive OFGs that are not LBGs ($\log(M_*/M_\odot) > 10.3$) are shaded in blue, where the sources at $z = 4-5$ are shown in faint red and blue filled circles, respectively. Grey filled circles are ALMA-detected massive H-dropouts with $\log(M_*/M_\odot) > 10.3$ (Wang et al. 2019). Grey filled triangles are massive LBGs with $\log(M_*/M_\odot) > 10.3$ whose SFRD is based on the dust-corrected UV (Wang et al. 2019). We also show the combined contribution of OFGs and LBGs to the cosmic SFRD. The red (blue) open circle indicates the combined contribution of OFGs (massive OFGs) and all LBGs (Madau & Dickinson 2014) at $z = 4-5$. **Right:** summary of the contribution of optically dark or faint galaxies to the cosmic SFRD from the literature. These optically dark or faint galaxies were selected in different methods, including our OFGs (red point; this work), ALMA-detected massive H-dropouts (grey filled circles; Wang et al. 2019), JWST-selected HST-dark galaxies from CEERS (blue area; Barrufet et al. 2022), ALMA-detected optically dark sources from REBELS (brown point and arrow; Fudamoto et al. 2021), an ALMA-detected NIR-dark galaxy in COSMOS (light grey point; Williams et al. 2019), ALMA-detected HST+near-IR dark galaxies from ALPINE (orange area; Gruppioni et al. 2020), ALMA-detected K-dropouts from ASAGAO (red area; Yamaguchi et al. 2019), radio 3 GHz-selected UV-dark galaxies in COSMOS (green area; Talia et al. 2021), and radio 1.4 GHz-selected HST-dark galaxies in GOODS-N (grey arrows; Enia et al. 2022). All data from the literature have been converted to a Chabrier (2003) IMF.</figcaption></figure>

studies of LBGs as well as H-dropouts. Therefore, we emphasize the importance of considering moderately obscured, non-LBG sources in complementing the high-$z$ SFRD measurements, rather than focusing only on those extremely dusty star-forming galaxies.

Table 7 presents the SFRD of the total OFGs, the massive OFGs, and the ALMA-detected OFGs, with LBGs excluded. Table 7 shows that these three (sub)-samples have the same SFRD values within errors. A comparison of the total OFGs and the massive OFGs has already been discussed above and shown in the left panel of Fig. 11. For the ALMA-detected OFGs, their SFRD value is lower than the total OFGs but still within errors. This implies that the investigation of the ALMA-detected sources alone might be sufficient to account for the contribution of OFGs in the SFRD. However, limited by the number of sources, this result needs to be further confirmed with a larger OFG sample in future studies. Besides, we would like to remind that the dust and gas properties of the ALMA-detected OFGs are not representative of all OFGs, as described in Sect. 5.4.

In our work, we find that the OFGs contribute ~52% to the reference SFRD at $z = 4-5$. Many works have also studied optically dark or faint galaxies selected or originally detected with different methodologies. The right panel of Fig. 11 shows a summary of the cosmic SFRD of optically dark or faint galaxies from the literature. Specifically, in Wang et al. (2019), the contribution of the ALMA-detected massive H-dropouts accounts for ~10% of the reference SFRD from LBGs at $z \sim 4-6$. At $z \sim 5$, this value is broadly consistent with that from the

ALMA-detected NIR-dark galaxy in COSMOS (Williams et al. 2019) and the ALMA-detected HST+near-IR dark galaxies from ALPINE (Gruppioni et al. 2020). Similar results have also been found for the ALMA-detected K-dropouts from ASAGAO (~10% at $z \sim 3-5$; Yamaguchi et al. 2019) and the JWST-selected HST-dark galaxies from CEERS (~10% at $z \sim 4-6$, but increasing to ~36% at $z \sim 7$; Barrufet et al. 2022). At $z \sim 7$, Fudamoto et al. (2021) identified two optically dark galaxies in two separate REBELS pointings that contribute ~10–25% to the reference SFRD. In radio, the 3 GHz-selected UV-dark galaxies from VLA-COSMOS (Talia et al. 2021) contribute 10–25% at $z \sim 3-4.5$ and 25–40% at $z > 4.5$ to the reference SFRD. More recently, Enia et al. (2022) found a 7–58% contribution from 1.4 GHz radio-selected HST-dark galaxies in GOODS-N at $z \sim 3$. Despite the broad range of contributions from these galaxies, the high values all agree that the contribution of optically dark or faint galaxies to the cosmic SFRD cannot be neglected, suggesting that highly dust-obscured star formation is relatively common in the $z > 3$ Universe.

Finally, we calculated the combined contribution of the OFGs and LBGs (values given from Madau & Dickinson 2014) to the cosmic SFRD. As shown by the red open circle in the left panel of Fig. 11, the total SFRD is $4.0 \times 10^{-2} M_\odot \text{ yr}^{-1} \text{ Mpc}^{-3}$ at $z = 4-5$, about 0.15 dex higher (43%) than the SFRD at the same redshift from Madau & Dickinson (2014). As we mention previously, the calculation of the SFRD values was very conservative and the true total SFRD of the OFGs and LBGs could be greater.
7. Conclusions

This work aims to obtain a more complete picture of the cosmic star formation history in the z > 3 Universe, that is, to bridge the extreme population of optically dark or faint galaxies (or H-dropouts) with the most common population of lower mass, less attenuated galaxies, such as LBGs. We use a more permis- sive criterion (H > 26.5 mag & [4.5] < 25 mag) to select opti- cally dark or faint galaxies (i.e., OFGs) at high redshifts, which avoids limiting the sample to the most extreme cases. This cri- terion selects extremely dust-obscured massive galaxies that are normal star-forming galaxies, with dust obscuration typically at E(B − V ) > 0.4, with lower stellar masses at high redshifts than H-dropouts. In addition, our selection method has the capacity to select OFGs without contamination from passive or old galaxies. In the GOODS-ALMA region, we have a total of 28 OFGs (includ- ing a candidate IRAC 4.5 µ m dropout). We analyzed the properties of individual and stacked OFGs, respectively. We calculated their SFRD and quantified the incompleteness of our understanding of the cosmic star-formation history in the z > 3 Universe.

Here are the main results of this work:

1. After performing SED analyses with the code EAZY and FAST++, we find that the OFGs cover a redshift of zphot > 3, as indicated by the theoretical galaxy templates. The median redshift of individual OFGs is zmed = 4.1, with a wide distri- bution (z = 3.2 − 7.0), consistent with zmed,stacked = 4.5 ± 0.2 derived from the median stacked SED. The OFGs have a broad stellar mass distribution with log(M*/M⊙) = 9.4 − 11.1, with a median of log(M*med/M⊙) = 10.3.

2. We investigate the proportions of LBGs, H-dropouts, and remaining OFGs (after removing LBGs and H-dropouts) in our sample. Remarkably, at stellar masses of log(M*/M⊙) = 9.5–10.5, the fraction of remaining OFGs is about three times the sum of LBGs and H-dropouts. In other words, up to 75% of the OFGs with log(M*/M⊙) = 9.5–10.5 at z > 3 are neglected by the previous LBGs and H-dropout selection techniques.

3. All stacked OFGs are located within the SFMS region (0.33 < ΔMS < 3, i.e., ~0.5 dex), which is consistent with being normal star-forming galaxies at the same redshift. It suggests that rather than being limited to a rare population of extreme starburst galaxies, our OFGs represent a normal population of dusty star-forming galaxies at z > 3.

4. The gas properties of the OFGs imply that the OFGs have shorter τ$_{dep}$ and slightly lower f$_{gas}$ values compared to the scaling relation followed by typical main sequence galax- ies. Their gas is consumed more rapidly, hence they form their stars with a high efficiency, setting them in the so- called class of starbursts in the main sequence. In addition, the massive OFGs have the shortest τ$_{dep}$ and lowest f$_{gas}$ of all stacked (sub)samples, indicating that they are in the process of becoming passive. Finally, we point out that the approach of selecting only ALMA-detected galaxies in stud- ies of OFGs will end up biasing the sample toward larger SFRs, longer τ$_{dep}$ , and larger f$_{gas}$ .

5. Studying dust temperatures in the OFGs, we find that, sur- prisingly, most of the stacked (sub)samples show colder T$_{dust}$ than the scaling relation. In particular, the ALMA-detected OFGs have the most abundant dust but show the lowest T$_{dust}$, indicating that the dust is colder in more obscured sources. A possible reason for the cold dust temperature is that the dust emission in the FIR of the dust-obscured sources may be optically thick rather than optically thin, where a warm and compact dust core is hidden.

6. In the comparison of SFRcorr UV and SFRtot, we find that (i) all the stacked (sub-)samples have SFRtot much larger than SFRcorr UV . There could be a hidden dust region of OFGs (even for LBGs as well) that absorbs all the UV photons, which cannot be reproduced with a dust extinction correction; and (ii) it is fundamental to include IR/mm band observations when studying extremely dusty star-forming galaxies; other- wise, the total SFR and, therefore, the cosmic SFRD will be underestimated.

7. After excluding five LBGs in the OFG sample, we study the hidden side of cosmic SFRD at high redshift (z ≥ 3) and find that: (i) among all galaxies, the total SFRD is dominated by LBGs, followed by OFGs. At z = 4–5, the contribution of OFGs reaches about 52% of the SFRD (Madau & Dickinson 2014), which is calculated mainly based on the LBGs (Bouwens et al. 2012a,b); (ii) for the OFGs, the massive and total OFGs have similar SFRD values, indicating that the massive OFGs make a major contribution to the SFRD compared to the remaining low-mass OFGs; (iii) for mas- sive galaxies, the SFRD is dominated by massive OFGs rather than massive LBGs. The SFRD contributed by mas- sive OFGs is at least two orders of magnitude higher than the one contributed by massive LBGs; (iv) the contribu- tion of massive OFGs to the SFRD is more than four times higher than that of H-dropouts (Wang et al. 2019). It implies that optically faint galaxies also contribute significantly to the SFRD, which has been neglected in previous studies of LBGs and H-dropouts; and (v) the ALMA-detected OFGs and the total OFGs have similar SFRD values, but with dif- ferent gas and dust properties (as we mentioned in point 4), such as τ$_{dep}$ , f$_{gas}$ , T$_{dust}$ , and R$_{e}$ (1.13 mm), which require atten- tion when studying the OFGs selected by ALMA detection alone.

8. Finally, we calculate the combined contribution of the OFGs and LBGs to the cosmic SFRD at z = 4–5, which is 4.0 × 10$^{-2}$ M$_{⊙}$ yr$^{-1}$ Mpc$^{-3}$, about 0.15 dex higher (43%) than the SFRD derived from LBGs alone (Madau & Dickinson 2014) at the same redshift. This value could be even larger as our calculation was very conservative.

Acknowledgements. We are very grateful to the anonymous referee for instruc- tive comments, which helped to improve the overall quality and strengthened the analyses of this work. We thank Alain Omont, Boris Sindhu Kalita, Dan- gning Hu, Shuowen Jin, Pascal Oesch, and Bing Luo for valuable discussions and suggestions that improved this paper. This paper makes use of the fol- lowing ALMA data: ADS/JAO.ALMA #2015.1.00543.S and ADS/JAO.ALMA #2017.1.00755.S. ALMA is a partnership of ESO (representing its member states), NSF (USA), and NINS (Japan), together with NRC (Canada), NSC, ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Repub- lic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ. This work is supported by the National Key Research and Develop- ment Program of China (No. 2017YFA0402703), and by the National Natural Science Foundation of China (No. 11733002, 12121003, 12192220, 12192222, 12173017, and 12141301). This work is also supported by the Programme National Cosmology et Galaxies (PNCG) of CNRS/INSU with INP and IN2P3, co-funded by CEA and CNES. M.-Y. X. acknowledges the support by China Scholarship Council (CSC). M.F. acknowledges NSF grant AST-2009577 and NASA JWST GO Program 1727. G.E.M. acknowledges the Villum Fonden research grants 13160 and 37440 and the Cosmic Dawn Center of Excel- lence funded by the Danish National Research Foundation under the grant No. 140. R.D. gratefully acknowledges support by the ANID BASAL projects ACE210002 and FB210003.

References

Barrufet, L., Oesch, P. A., Weibel, A., et al. 2022, MNRAS, submitted, [arXiv:2207.14733]

Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393
Boquien, M., Burgarella, D., Roehlly, Y., et al. 2019, A&A, 622, A103

Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2012a, ApJ, 754, 83

Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2012b, ApJ, 752, L5

Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2015, ApJ, 803, 34

Bouwens, R., González-López, J., Aravena, M., et al. 2020, ApJ, 902, 112

Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, ApJ, 686, 1503

Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000

Burgarella, D., Buat, V., & Iglesias-Páramo, J. 2005, MNRAS, 360, 1413

Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682

Casey, C. M. 2012, MNRAS, 425, 3094

Chabrier, G. 2003, PASP, 115, 763

Cortzen, I., Magdis, G. E., Valentino, F., et al. 2020, A&A, 634, L14

Cowie, L. L., González-López, J., Barger, A. J., et al. 2018, ApJ, 865, 106

Daddi, E., Cimatti, A., Renzini, A., et al. 2004, ApJ, 617, 746

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

Delvecchio, I., Daddi, E., Sargent, M. T., et al. 2021, A&A, 647, A123

Dickinson, M., Giavalisco, M., & GOODS Team 2003, The Mass of Galaxies at Low and High Redshift, 324

Donevski, D., Lapi, A., Małek, K., et al. 2020, A&A, 644, A144

Donley, J. L., Koekemoer, A. M., Brusa, M., et al. 2012, ApJ, 748, 142

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

Draine, B. T., Aniano, G., Krause, O., et al. 2014, ApJ, 780, 172

Dudzevičiūtė, U., Smail, I., Swinbank, A. M., et al. 2020, MNRAS, 494, 3828

Elbaz, D., Daddi, E., Le Borgne, D., et al. 2007, A&A, 468, 33

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

Enia, A., Talia, M., Pozzi, F., et al. 2022, ApJ, 927, 204

Fontana, A., Dunlop, J. S., Paris, D., et al. 2014, A&A, 570, A11

Franco, M., Elbaz, D., Béthermin, M., et al. 2018, A&A, 620, A152

Franco, M., Elbaz, D., Zhou, L., et al. 2020a, A&A, 643, A53

Franco, M., Elbaz, D., Zhou, L., et al. 2020b, A&A, 643, A30

Fritz, J., Franceschini, A., & Hatziminaoglou, E. 2006, MNRAS, 366, 767

Fudamoto, Y., Oesch, P. A., Schouws, S., et al. 2021, Nature, 597, 489

Genzel, R., Tacconi, L. J., Lutz, D., et al. 2015, ApJ, 830, 20

Giavalisco, M., Ferguson, H. C., Koekemoer, A. M., et al. 2004, ApJ, 600, L93

Gómez-Guijarro, C., Toft, S., Karim, A., et al. 2018, ApJ, 856, 121

Gómez-Guijarro, C., Elbaz, D., Xiao, M., et al. 2022a, A&A, 658, A43

Gómez-Guijarro, C., Elbaz, D., Xiao, M., et al. 2022b, A&A, 659, A196

Gordon, K. D., Galliano, F., Hony, S., et al. 2010, A&A, 518, L89

Gruppioni, C., Béthermin, M., Loiacono, F., et al. 2020, A&A, 643, A8

Helou, G., Soifer, B. T., & Rowan-Robinson, M. 1985, ApJ, 298, L7

Herenz, E. C., Urrutia, T., Wisotzki, L., et al. 2017, A&A, 606, A12

Hildebrand, R. H. 1983, QJRAS, 24, 267

Hsieh, B.-C., Wang, W.-H., Hsieh, C.-C., et al. 2012, ApJS, 203, 23

Hwang, H. S., Elbaz, D., Magdis, G., et al. 2010, MNRAS, 409, 75

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

Jin, S., Daddi, E., Magdis, G. E., et al. 2019, ApJ, 887, 144

Jin, S., Daddi, E., Magdis, G. E., et al. 2022, A&A, 665, A3

Kennicutt, R. C., & Evans, N. J. 2012, ARA&A, 50, 531

Kovács, A., Chapman, S. C., Dowell, C. D., et al. 2006, ApJ, 650, 592

Kriek, M., van Dokkum, P. G., Labbé, I., et al. 2009, ApJ, 700, 221

Lee, N., Sanders, D. B., Casey, C. M., et al. 2015, ApJ, 801, 80

Leslie, S. K., Schinnerer, E., Liu, D., et al. 2020, ApJ, 899, 58

Lilly, S. J., Eales, S. A., Gear, W. K. P., et al. 1999, ApJ, 518, 641

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

Luo, B., Brandt, W. N., Xue, Y. Q., et al. 2017, ApJS, 228, 2

Lutz, D., Poglitsch, A., Altieri, B., et al. 2011, A&A, 532, A90

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

Magdis, G. E., Rigopoulou, D., Huang, J.-S., et al. 2010, MNRAS, 401, 1521

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

Magnelli, B., Elbaz, D., Chary, R. R., et al. 2011, A&A, 528, A35

Magnelli, B., Popesso, P., Berta, S., et al. 2013, A&A, 553, A132

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, Nature, 553, 51

McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, ASP Conf. Ser., 376, 127

Noeske, K. G., Weiner, B. J., Faber, S. M., et al. 2007, ApJ, 660, L43

Noll, S., Burgarella, D., Giovannoli, E., et al. 2009, A&A, 507, 1793

Nonino, M., Dickinson, M., Rosati, P., et al. 2009, ApJS, 183, 244

Oesch, P. A., Bouwens, R. J., Illingworth, G. D., et al. 2014, ApJ, 786, 108

Oesch, P. A., Bouwens, R. J., Illingworth, G. D., et al. 2015, ApJ, 808, 104

Oesch, P. A., Bouwens, R. J., Illingworth, G. D., et al. 2018, ApJ, 855, 105

Oke, J. B., & Gunn, J. E. 1983, ApJ, 266, 713

Orellana, G., Nagar, N. M., Elbaz, D., et al. 2017, A&A, 602, A68

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

Retzlaff, J., Rosati, P., Dickinson, M., et al. 2010, A&A, 511, A50

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

Riechers, D. A., Hodge, J. A., Pavesi, R., et al. 2020, ApJ, 895, 81

Rujopakarn, W., Dunlop, J. S., Rieke, G. H., et al. 2016, ApJ, 833, 12

Salpeter, E. E. 1955, ApJ, 121, 161

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

Schreiber, C., Labbé, I., Glazebrook, K., et al. 2018a, A&A, 611, A22

Schreiber, C., Glazebrook, K., Nanayakkara, T., et al. 2018b, A&A, 618, A85

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

Shu, X., Yang, L., Liu, D., et al. 2022, ApJ, 926, 155

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

Smail, I., Dudzevičiūtė, U., Stach, S. M., et al. 2021, MNRAS, 502, 3426

Sommovigo, L., Ferrara, A., Pallottini, A., et al. 2022, MNRAS, 513, 3122

Speagle, J. S., Steinhardt, C. L., Capak, P. L., et al. 2014, ApJS, 214, 15

Stefanon, M., Labbé, I., Oesch, P. A., et al. 2021, ApJS, accepted [arXiv:2110.06226]

Straatman, C. M. S., Spitler, L. R., Quadri, R. F., et al. 2016, ApJ, 830, 51

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

Tacconi, L. J., Genzel, R., Saintonge, A., et al. 2018, ApJ, 853, 179

Talia, M., Cimatti, A., Giulietti, M., et al. 2021, ApJ, 909, 23

Urrutia, T., Wisotzki, L., Kerutt, J., et al. 2019, A&A, 624, A141

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

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

Wang, T., Schreiber, C., Elbaz, D., et al. 2019, Nature, 572, 211

Whitaker, K. E., van Dokkum, P. G., Brammer, G., et al. 2012, ApJ, 754, L29

Whitaker, K. E., Franx, M., Leja, J., et al. 2014, ApJ, 795, 104

Whitaker, K. E., Ashas, M., Illingworth, G., et al. 2019, ApJS, 244, 16

Williams, C. C., Labbe, I., Spilker, J., et al. 2019, ApJ, 884, 154

Yamaguchi, Y., Kohno, K., Hatsukade, B., et al. 2019, ApJ, 878, 73

Yun, M. S., Reddy, N. A., & Condon, J. J. 2001, ApJ, 554, 803

Zhou, L., Elbaz, D., Franco, M., et al. 2020, A&A, 642, A155

<footnote n="1">¹ School of Astronomy and Space Science, Nanjing University, Nanjing 210093, PR China e-mail: my.xiao@sjtu.edu.cn</footnote>

<footnote n="2">² AIM, CEA, CNRS, Université Paris-Saclay, Université Paris Diderot, Sorbonne Paris Cité, 91191 Gif-sur-Yvette, France</footnote>

<footnote n="3">³ Key Laboratory of Modern Astronomy and Astrophysics (Nanjing University), Ministry of Education, Nanjing 210093, PR China</footnote>

<footnote n="4">⁴ Aix- Marseille Université, CNRS, LAM, Laboratoire d'Astrophysique de Marseille, 13013 Marseille, France</footnote>

<footnote n="5">⁵ Department of Astronomy, The University of Texas at Austin, 2515 Speedway Blvd Stop C1400, Austin, TX 78712, USA</footnote>

<footnote n="6">⁶ Community Science and Data Center/NSF's NOIRLab, 950 N. Cherry Ave., Tucson, AZ 85719, USA</footnote>

<footnote n="7">⁷ Department of Physics, Faculty of Science, Chulalongkorn University, 254 Phayathai Road, Pathumwan, Bangkok 10330, Thailand</footnote>

<footnote n="8">⁸ National Astronomical Research Institute of Thailand (Public Organization), Don Kaeo, Mae Rim, Chiang Mai 50180, Thailand</footnote>

<footnote n="9">⁹ Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan</footnote>

<footnote n="10">¹⁰ Cosmic Dawn Center (DAWN), Jagtvej 128, 2200 Copenhagen N, Denmark</footnote>

<footnote n="11">¹¹ DTU-Space, Technical University of Denmark, Elektrovej 327, 2800 Kgs. Lyngby, Denmark</footnote>

<footnote n="12">¹² Niels Bohr Institute, University of Copenhagen, Jagtvej 128, 2200 Copenhagen N, Denmark</footnote>

<footnote n="13">¹³ Instituto de Astrofísica, Facultad de Física, Pontificia Universidad Católica de Chile, Casilla 306, Santiago 22, Chile</footnote>

<footnote n="14">¹⁴ Hiroshima Astrophysical Science Center, Hiroshima University, 1-3-1 Kagamiyama, Higashi-Hiroshima, Hiroshima 739-8526, Japan</footnote>

<footnote n="15">¹⁵ Departamento de Astronomía, Facultad de Ciencias Físicas y Matemáticas, Universidad de Concepción, Concepción, Chile</footnote>

<footnote n="16">¹⁶ Astronomy Centre, Department of Physics and Astronomy, University of Sussex, Brighton BN1 9QH, UK</footnote>

<footnote n="17">¹⁷ International Space Science Institute (ISSI), Hallerstrasse 6, CH-3012 Bern, Switzerland</footnote>

<footnote n="18">¹⁸ Department of Physics, Anhui Normal University, Wuhu, Anhui 241000, PR China</footnote>

<footnote n="19">¹⁹ Laboratory for Multiwavelength Astrophysics, School of Physics and Astronomy, Rochester Institute of Technology, 84 Lomb Memorial Drive, Rochester, NY 14623, USA</footnote>
20 Centre for Extragalactic Astronomy, Department of Physics, Durham University, Durham DH1 3LE, UK

21 Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA

22 Astronomy Department, University of Massachusetts, Amherst, MA 01003, USA

23 National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan

24 Department of Astronomical Science, SOKENDAI (The Graduate University for Advanced Studies), Mitaka, Tokyo 181-8588, Japan

25 Joint ALMA Observatory, Alonso de Córdova 3107, Vitacura 763-0355, Santiago, Chile

26 European Southern Observatory, Alonso de Córdova 3107, Vitacura, Casilla 19001, 19 Santiago, Chile

27 Institute of Astronomy, Graduate School of Science, The University of Tokyo, 2-21-1 Osawa, Mitaka, Tokyo 181-0015, Japan

28 Department of Physics and Astronomy, University of Sheffield, Sheffield S3 7RH, UK

29 Department of Astronomy, Universidad de Concepción, Casilla 160-C, Concepción, Chile

30 Astronomy Unit, Department of Physics, University of Trieste, via Tiepolo 11, 34131 Trieste, Italy

31 Fakultät für Physik der Ludwig-Maximilians-Universität, 81679 München, Germany

32 Department of Physics and Astronomy, Texas A&M University, College Station, TX 77843-4242, USA

33 George P. and Cynthia Woods Mitchell Institute for Fundamental Physics and Astronomy, Texas A&M University, College Station, TX 77843-4242, USA

34 Astrophysics, Department of Physics, Keble Road, Oxford OX1 3RH, UK
# Appendix A: IRAC catalog compilation

We constructed the IRAC catalog using the deepest IRAC 3.6 and 4.5 µm images from the GREATS program (Stefanon et al. 2021). Source detection was performed using Source Extractor (SE version 2.25.0; Bertin & Arnouts 1996) on the background-subtracted 3.6 and 4.5 µm images, respectively. This is to help us exclude as many false detections as possible and get pure sources after cross-matching two catalogs from 3.6 and 4.5 µm images. Each detected source was required to have a minimum area of 4 pixels, with each pixel achieving the threshold of 0.25σ. The deblending parameters were optimized by setting `DEBLEND_THRESHOLD = 64` and `DEBLEND_MINCONT = 0.0001`, which were a compromise between deblending neighboring sources and minimizing a potential division of a larger object into multiple components. The sky subtraction was performed with SE, using a bicubic interpolation of the background with an adopted mesh size of 512 pixels and a median filter size of 5 pixels. A Mexhat filter was used to smooth the images before detection, to help detect faint and extended objects. The Mexhat filter can very well de-blend faint sources around the bright source, but it also extracts artifacts around the bright source caused by diffraction, which needs to be cleaned. We set `CLEAN_PARAM = 1` to clean spurious detections with SE. Then, we visually inspected images to further remove artifacts around the bright sources. For the 3.6 and 4.5 µm images, we obtained their source catalogs respectively. The total sources in the catalogs are 125,338 for 3.6 µm and 154,234 for 4.5 µm in the entire GOODS-S field.

To ensure the purity of detections, we then cross-matched two catalogs from 3.6 and 4.5 µm images with a radius of 1.0" (~0.5×FWHM). For sources that are simultaneously in the GOODS-ALMA 2.0 catalog (Gómez-Guijarro et al. 2022a) and detected in at least one IRAC band, we considered them as real sources and kept them in the final catalog. Finally, we have 5,127 pure sources detected by both 3.6 and 4.5 µm bands and/or ALMA 1.13 mm in the GOODS-ALMA field.

To estimate the completeness of our detection strategy, we employed a Monte Carlo approach where we simulated 10,000 artificial model sources with random magnitudes between 14.5 mag and 29 mag. They are point sources with the same PSF profiles as 3.6 and 4.5 µm images, respectively. The simulated sources were allowed to fall at random positions on the real image, including on top of other sources, to account for the impact of source blending in the real image. We injected ten sources each time on the IRAC 3.6µm and 4.5µm images to avoid excessively artificial source confusion caused by bright fake sources. After each injection, we performed the same blind source detection procedure using SE.

In Fig. A.1 (left), we show the completeness as a function of the input flux density ($S_{\text{in}}$) for the simulated sources in 3.6 and 4.5 µm images. The survey reaches a 100% completeness for all simulated sources for flux densities $S_{\text{in}} \ge 20$ µJy at both 3.6 and 4.5 µm wavelength. To understand the incompleteness of our OFG sample, we need to know the behavior of completeness as a function of output flux density ($S_{\text{out}}$) shown in Fig. A.1 (right). Here, $S_{\text{out}}$ is the parameter being measured directly.

<figure><img src="image_13.png" /><figcaption>Fig. A.1. Completeness as a function of the input flux density ($S_{\text{in}}$; top) and the output flux density ($S_{\text{out}}$; bottom). The dashed lines represent the 50% and 80% completeness limits as a reference.</figcaption></figure>
Appendix B: Multiwavelength Postage-stamp

<figure><img src="image_14.png" /><figcaption>**Fig. B.1.** Postage-stamps (4" × 4") of our OFGs. From left to right: HST/WFC3 (F160W), ZFOURGE ($K_s$), *Spitzer*/IRAC (3.6 µm and 4.5 µm), and ALMA band 6 (1.13 mm). The white circles mark our target galaxies with 0.6" in radius. The black circles indicate the positions of the ALMA detections. North is up and east is to the left.</figcaption></figure>
Appendix C: SEDs for the individual OFGs

<figure><img src="image_15.png" /><figcaption>**Fig. C.1.** UV to millimeter SEDs for 27 OFGs. The symbol convention and lines are the same as those in Figs. 7. The data from the UV to MIR (16 µm) bands are fitted with the FAST++ code (see §4.1). From 24 µm up to millimeter wavelengths, we fitted galaxies with a *Herschel* counterpart using the CIGALE code (orange solid line) and galaxies without a *Herschel* counterpart but with the ALMA 1.13 mm detection using the dust template libraries (green line; Schreiber et al. 2018c, see §4.2).</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Fig. C.2. continued.</figcaption></figure>