# Illuminating the Dark Side of Cosmic Star Formation.
## II. A Second Date with RS-NIRdark Galaxies in COSMOS

Meriem Behiri, Margherita Talia, Andrea Cimatti, Andrea Lapi, Marcella Massardi, Andrea Enia, Cristian Vignali, Matthieu Bethermin, Andreas Faisst, Fabrizio Gentile, et al.

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

To cite this version:

Meriem Behiri, Margherita Talia, Andrea Cimatti, Andrea Lapi, Marcella Massardi, et al. Illuminating the Dark Side of Cosmic Star Formation. II. A Second Date with RS-NIRdark Galaxies in COSMOS. The Astrophysical Journal, 2023, 957, 10.3847/1538-4357/acf616. insu-04271569

HAL Id: insu-04271569

https://insu.hal.science/insu-04271569

Submitted on 9 Nov 2023

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Illuminating the Dark Side of Cosmic Star Formation. II. A Second Date with RS-NIRdark Galaxies in COSMOS

Meriem Behiri<fnref n="1" />,<fnref n="2" />,<fnref n="3" /> ID, Margherita Talia<fnref n="2" />,<fnref n="4" /> ID, Andrea Cimatti<fnref n="2" />,<fnref n="5" /> ID, Andrea Lapi<fnref n="1" />,<fnref n="3" />,<fnref n="6" />,<fnref n="7" /> ID, Marcella Massardi<fnref n="1" />,<fnref n="3" />,<fnref n="8" /> ID,
Andrea Enia<fnref n="2" />,<fnref n="4" /> ID, Cristian Vignali<fnref n="2" />,<fnref n="4" /> ID, Matthieu Bethermin<fnref n="9" />,<fnref n="10" /> ID, Andreas Faisst<fnref n="11" /> ID, Fabrizio Gentile<fnref n="2" />,<fnref n="4" /> ID, Marika Giulietti<fnref n="1" />,<fnref n="4" /> ID,
Carlotta Gruppioni<fnref n="4" /> ID, Francesca Pozzi<fnref n="2" />,<fnref n="4" /> ID, Vernesa Smolčić<fnref n="12" /> ID, and Gianni Zamorani<fnref n="4" /> ID

Received 2022 October 27; revised 2023 August 17; accepted 2023 August 31; published 2023 October 30

## Abstract

About 12 billion years ago, the Universe was first experiencing light again after the dark ages, and galaxies filled the environment with stars, metals, and dust. How efficient was this process? How fast did these primordial galaxies form stars and dust? We can answer these questions by tracing the star formation rate density (SFRD) back to its widely unknown high-redshift tail, traditionally observed in the near-infrared (NIR), optical, and UV bands. Thus, objects with a large amount of dust were missing. We aim to fill this knowledge gap by studying radio-selected NIR-dark (RS-NIRdark) sources, i.e., sources not having a counterpart at UV-to-NIR wavelengths. We widen the sample of Talia et al. from 197 to 272 objects in the Cosmic Evolution Survey (COSMOS) field, including also photometrically contaminated sources, which were previously excluded. Another important step forward consists in the visual inspection of each source in the bands from $u^*$ to MIPS 24 µm. According to their “environment” in the different bands, we are able to highlight different cases of study and calibrate an appropriate photometric procedure for the objects affected by confusion issues. We estimate that the contribution of RS-NIRdark sources to the cosmic SFRD at $3 < z < 5$ is ~10%–25% of that based on UV-selected galaxies.

*Unified Astronomy Thesaurus concepts: Extragalactic radio sources (508); Galaxy formation (595); Galaxy evolution (594); High-redshift galaxies (734); Star formation (1569)*

## 1. Introduction

Understanding the cosmic star formation history is fundamental to improving our knowledge of galaxy formation and evolution, the production of heavy elements, and the reionization process through cosmic time. In the 1990s Lilly et al. (1996) and Madau et al. (1996) built the Lilly–Madau plot with Hubble Space Telescope (HST) data to study the evolution of the star formation rate (SFR) per unit of cosmic volume (i.e., the star formation rate density, SFRD) as a function of redshift. Thanks to the critical technological improvements, scientists have added more observational points to this plot, improving the constraints on the SFRD and pushing this study toward higher redshifts. Madau & Dickinson (2014) show an evolution of the SFRD $\propto (1+z)^{1.7}$ up to $z \sim 2$, the so-called cosmicnoon, while between $z \sim 3$ and $z \sim 8$ the SFRD decreases $\propto (1+z)^{-2.9}$. The evolution of the SFRD through cosmic time is now well constrained until $z \sim 3$. However, the possible decline of the SFRD after $z \sim 3$ is still an open question (Casey et al. 2014; Magnelli et al. 2015; Zavala et al. 2021), and some

works argue that it might flatten (Gruppioni et al. 2013, 2020; Rowan-Robinson et al. 2016; Novak et al. 2017).

Many high-redshift studies of the SFRD rely mainly on rest-frame UV/optically bright sources, mostly from Lyman-break galaxies (LBGs). The lack of dusty galaxies in the census introduces biases in the estimations of the SFRD since the LBGs represent only a fraction of the whole population of galaxies, and the computation of the SFR for these sources depends strongly on the dust corrections adopted. We have recently witnessed a surge of candidates at very high z discovered by the James Webb Space Telescope (JWST) (e.g., Adams et al. 2023; Atek et al. 2023; Castellano et al. 2022; Finkelstein et al. 2022; Naidu et al. 2022; Yan et al. 2023); however, these measurements are mostly in the UV rest frame, therefore a dust-unbiased view of the high-z Universe is still fundamental to reconstruct the SFRD.

On the one hand, Rowan-Robinson et al. (2016), using IR data from the Herschel Telescope, found that, at $z > 3$, the contribution of the dusty galaxies to the SFRD could be about ten times that due to UV-bright sources. On the other hand, studying high-redshift galaxies in the infrared and submillimeter bands has always been technically tricky since observations are deeply affected by confusion problems (Lutz et al. 2011; Magnelli et al. 2013). Moreover, the faintness and the lack of

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹SISSA, Via Bonomea 265, I-34136 Trieste, Italy; mbehiri@sissa.it</footnote>

<footnote n="2">²University of Bologna, Department of Physics and Astronomy “Augusto Righi” (DIFA), Via Gobetti 93/2, I-40129 Bologna, Italy</footnote>

<footnote n="3">³INAF/IRA, Istituto di Radioastronomia, Via Gobetti 101, I-40129, Bologna, Italy</footnote>

<footnote n="4">⁴INAF—Osservatorio di Astrofisica e Scienza dello Spazio, Via Gobetti 93/3, I-40129, Bologna, Italy</footnote>

<footnote n="5">⁵INAF—Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, I-50125, Firenze, Italy</footnote>

<footnote n="6">⁶INFN—Sezione di Trieste, via Valerio 2, Trieste I-34127, Italy</footnote>

<footnote n="7">⁷IFPU—Institute for fundamental physics of the Universe, Via Beirut 2, I-34014 Trieste, Italy</footnote>

<footnote n="8">⁸INAF, Istituto di Radioastronomia—Italian ARC, Via Piero Gobetti 101, I-40129 Bologna, Italy</footnote>

<footnote n="9">⁹Aix Marseille Univ, CNRS, LAM, Laboratoire d'Astrophysique de Marseille, Marseille, France</footnote>

<footnote n="10">¹⁰Université de Strasbourg, CNRS, Observatoire astronomique de Strasbourg, UMR 7550, F-67000 Strasbourg, France</footnote>

<footnote n="11">¹¹IPAC, M/C 314-6, California Institute of Technology, 1200 East California Boulevard, Pasadena, CA 91125, USA</footnote>

<footnote n="12">¹²Department of Physics, Faculty of Science, University of Zagreb, Bijenićka cesta 32, 10000 Zagreb, Croatia</footnote>
counterparts in the UV/optical/NIR bands made identifying dusty galaxies at $z > 3$ even more challenging. The first deep surveys in the (sub)millimetric bands opened the path to studies of these galaxies at high redshift, taking advantage of the negative *k*-correction in this part of the spectrum (Smail et al. 1997; Hughes et al. 1998). Smail et al. (1997) observed for the first time dusty submillimeter galaxies with the Submillimeter Common-User Bolometric Array (SCUBA) instrument of the James Clerk Maxwell Telescope. These are the well-known submillimeter galaxies, i.e., galaxies with $S_{850\mu m} > 3$ mJy (Hayward et al. 2011), and many of them were found to be at high redshift. The South Pole Telescope has made a notable contribution to this field, allowing the identification of ~40 candidate dusty star-forming galaxies (DSFGs) at $5 < z < 6$ (Strandet et al. 2016), together with the Atacama Large Millimeter/submillimeter Array (ALMA). ALMA enabled the detection of many of these objects by tracing dust in the continuum (e.g., Scoville et al. 2014, 2016, 2017; Brisbin et al. 2017; González-López et al. 2017; Franco et al. 2018; Dudzevičiūtė et al. 2020; Smail et al. 2021).

In fact, prior to the arrival of ALMA data, a thorough comprehension of this kind of object was not conceivable (e.g., Dannerbauer et al. 2004; Frayer et al. 2004). For instance, thanks to the high spatial and spectral resolution of ALMA, Walter et al. (2012) spectroscopically confirmed the redshift of the starburst galaxy HDF 850.1, detected with SCUBA by Hughes et al. (1998). As more data have become available, the idea that most high-redshift dust-obscured galaxies were rare, incredibly star-forming or active galactic nuclei (AGNs) has been revisited, since normal galaxies started to emerge from the dust (e.g., Simpson et al. 2014, 2017, 2020; Franco et al. 2018; Wang et al. 2019), and they could be more common than expected. These studies confirmed the hypothesis that optical/NIR studies missed a significant fraction of objects in the census of high-redshift galaxies. Nevertheless, (sub)millimeter and far-infrared (FIR) surveys also have some drawbacks. These bands sample the dust graybody: since the dust temperature might correlate with infrared luminosity and redshift (Béthermin et al. 2015; Faisst et al. 2017; Schreiber et al. 2018b), this could generate a selection bias. Furthermore, ALMA deep fields are still not wide enough (~4.5–10 arcmin², Aravena et al. 2016; Walter et al. 2016; Dunlop et al. 2017; Franco et al. 2018; Hatsukade et al. 2020) to provide a complete view of these objects. In the end, since every selection method affects the results differently, combining different approaches and exploring new strategies is vital.

Primordial dusty galaxies are now ubiquitous in theoretical models of galaxy formation as possible progenitors of massive red and dead galaxies (Cimatti et al. 2004, 2008; Cimatti 2008; Casey et al. 2014; Mancuso et al. 2016a, 2016b; Lapi et al. 2018). Moreover, it has recently been suggested that these galaxies could be the most likely hosts of high-redshift gravitational waves that the Laser Interferometer Space Antenna and the Einstein Telescope will study in the near future (Boco et al. 2020). Our work fits into this context, aiming to provide a deeper understanding of the role of obscured galaxies in galaxy formation and evolution. We propose an alternative approach to the traditional ones used to select obscured galaxies (e.g., (sub) millimeter and FIR selections), along the path of Talia et al. (2021, hereafter T21). They selected radio sources in the Cosmic Evolution Survey (COSMOS) field having no (or extremely faint) NIR counterparts (see Section 2). We call these sources

radio-selected NIR-dark (RS-NIRdark) galaxies (see also Chapman et al. 2001, 2002, 2004; Enia et al. 2022). The value of $K_{s,\text{sim}}$ adopted is ~24.7, which is the $K_s$ limiting flux density at $3\sigma$ of the COSMOS2015 catalog (Laigle et al. 2016).

Radio-selected sources have many advantages: radio interferometers reach high sensitivity and resolution, plus dust temperature bias does not affect the selection. The most significant bias possibly affecting radio selection is the possible presence of AGNs, but we can deal with this problem thanks to multiwavelength ancillary data. Galaxies in our study are, by the selection, not (or barely) detected in optical and NIR bands; therefore, we need mid-infrared (MIR), FIR, and millimetric data in order to reconstruct their spectral energy distributions (SEDs). Thus, the availability of multiwavelength catalogs is fundamental to studying these sources. Another drawback of the radio selection is the effect of the positive *k*-correction, which makes going toward higher redshifts more challenging. For these reasons, we selected our sources from the VLA-COSMOS 3 GHz Large Project (Smolčić et al. 2017) catalog, a deep radio survey in the COSMOS field that benefits from many ancillary data.

This paper is organized as follows. In Section 2 we explain the sample selection process and the features of the different (sub)samples. In Section 3 we describe the median properties of the sample. In Section 4 we present the study of the individual photometry of our objects, with a particular focus on the new deblending procedure and on the presence of 500 µm risers. In Section 5 we discuss the results from the SED fitting of the individual sources. In Section 6 we derive the SFRD and the number density for our sample.

Throughout this paper we adopt a Chabrier (2003) initial mass function (IMF) and we assume a cosmology with matter density $\Omega_M = 0.3$, dark energy density $\Omega_\Lambda = 0.7$, and Hubble parameter $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$. Magnitudes are given in the AB system.

## 2. Sample Selection

We selected the galaxies studied in this work from the VLA-COSMOS 3 GHz Large Project (Smolčić et al. 2017) catalog, based on a survey covering 2.6 deg² in the COSMOS field in 348 hr of observations with the Karl G. Jansky Very Large Array (VLA) interferometer at 3 GHz. Thanks to its limiting flux density of $12.6 \ \mu\text{Jy beam}^{-1}$ at $5.5\sigma$, this data set is one of the deepest ever obtained, and suitable for estimating the contribution to the cosmic SFRD of galaxies at least up to $z \sim 5$, as shown in Novak et al. (2017) and T21. We summarize the selection process in Figure 1.

This work inherits the first steps of the selection process from T21. They selected 8850 sources detected with a signal-to-noise ratio S/N > 5.5 from the VLA-COSMOS 3 GHz (Smolčić et al. 2017). They choose this S/N threshold to minimize the fraction of spurious sources. Subsequently, T21 excluded all the sources within the bad areas of the VISTA Ultra-Deep (UltraVISTA) field (Laigle et al. 2016), reducing the effective area to 1.38 deg², and discarded all the multicomponent objects. They cross-correlated the resulting catalog of 5982 sources with the COSMOS2015 catalog (Laigle et al. 2016) within a search radius of 0.8 (Smolčić et al. 2017), finding 476 sources without a counterpart. These are the RS-NIRdark (according to the $K_{s,\text{sim}}$ of COSMOS2015, i.e., ~24.7 mag) galaxies, which constitute the starting point of this work. As an additional step, T21 narrowed down their selection
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Flow chart representing the selection process of the sample. The dark blue part of the chart indicates those steps made by T21, while the purple contours are used to highlight the ones made in this work.</figcaption></figure>

to the 197 sources whose $3\sigma$ isophotes in the 3 GHz and $\chi^2$ maps,<fnref n="13" /> respectively, do not intersect. In our work we relax this constraint and include all the sources whose $5\sigma$ 3 GHz isophote does not intersect the $5\sigma$ NIR isophote—a total of 272 objects. Thus, T21 selection is stricter than ours and includes only the most isolated galaxies, while we extend the analysis to more contaminated sources, which requires the use of more refined photometric techniques. The analysis of the remaining 204 sources, which are heavily contaminated by nearby NIR-bright galaxies, will be presented in a forthcoming paper (F. Gentile et al., in preparation).

We split the 272 objects into four subsamples based on visual inspection of the NIR band maps from the second UltraVISTA Data Release (DR2, McCracken et al. 2012) and

InfraRed Array Camera (IRAC) and Multiband Imaging Photometer for Spitzer (MIPS) maps from Spitzer<fnref n="14" />:

1. *Type 0:* Isolated galaxies, with no counterparts up to MIPS 24 µm.

2. *Type 1:* Isolated galaxies in both the $\chi^2$ and IRAC-1 maps, having a counterpart in one or more bands.

3. *Type 2:* Galaxies with a counterpart in one or more bands which are also contaminated by nearby sources. Objects belonging to this class might show a distinct brightness peak, with respect to the contaminant, in all the bands, or they might present a distinct peak in the optical and NIR bands but be unresolved in the MIR images.

4. *Type 3:* Like Type 0, but close to another radio source with a counterpart in one or more bands.

We display examples and explanations of the different subsamples in Figure 2.

<footnote n="13">¹³ The NIR detection map used to build the COSMOS2015 catalog is a $\chi^2$ composite image of the $z^+$, Y, J, H, and Ks maps, where each pixel has a certain probability of belonging or not to the background. Knowing the probability distribution of the background pixels (a $\chi^2$ distribution), it is possible to use the data to obtain the distribution of the pixels dominated by the source emission and to determine the best threshold to identify the background pixels (Szalay et al. 1999).</footnote>

<footnote n="14">¹⁴ We used the UltraVISTA DR2 maps for the selection for consistency with T21. The IRAC maps are the v2.0 mosaics of the Spitzer Large Area Survey with Hyper-Suprime-Cam (SPLASH; Steinhardt et al. 2014): https://splash.caltech.edu/.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Galaxies belonging to each of the four subsamples. We present examples for each group and show their χ² and IRAC-1 cutouts, on which we superimposed the radio (solid lines) and the χ² and IRAC-1 (dashed lines) contours at 3σ, 4σ, and 5σ (crimson). From the upper left: (a) Type 0: galaxies with no counterparts up to MIPS 24 µm; (b) Type 1: isolated galaxies in all bands and with a counterpart in the MIR filters; (c) Type 2: galaxies having a counterpart and a nearby contaminant in one or more bands; (d) Type 3: like Type 0, but close to another radio source with a counterpart in one or more bands.</figcaption></figure>

We stress that T21 selected the most isolated galaxies based on inspection of the radio and NIR maps only. This means that a fraction of T21 sources in fact presented some, albeit low, degree of contamination by a nearby source in the MIR bands. In this work, thanks to the visual inspection of the maps at all wavelengths, we have not only enlarged the sample but also improved the overall photometry with respect to the previous study.

We matched our catalog also to the Herschel Super-Deblended Catalog by Jin et al. (2018) for the FIR bands and, where available, the Automated Mining of the ALMA Archive in the COSMOS Field (A3COSMOS; Liu et al. 2019) for the (sub)millimeter bands. Finally, we introduce an additional classification based on the FIR/(sub)millimeter coverage of our sources:

1. *Primary sample:* sources having at least one detection above 3σ in one of the FIR/(sub)millimeter bands. This accounts for 145 sources.

2. *Secondary sample:* sources with no detections above 3σ in any FIR bands. This accounts for 127 sources.

Since we want to estimate the contribution to the cosmic SFRD of the RS-NIRdark sources, we only consider the former class since it leads to better constraints on redshift and the star formation rate, which we compute from the IR luminosity.

# 3. Median Properties of the Sample

## 3.1. Median Photometry

We start by considering the average properties of the sample as a whole through statistical analysis, by performing a median

stacking. In this procedure, we exclude Type-3 sources, because these might be jets belonging to the nearby radio source rather than galaxies themselves. After the individual analysis of the sources (Section 5), we have considered the possibility of performing a tomographic stacking, dividing the sample into redshift bins. Nevertheless, the number of sources in the high-redshift bins would be tiny, and thus the median stacking would not have the desired effect.

We build the median maps from $u^*$ to MIR by stacking individual maps at the VLA 3 GHz coordinates of the radio source. In particular, we used: Y, J, H, and $K_s$ maps from UltraVISTA DR4,<fnref n="15" /> IRAC maps from SPLASH<fnref n="16" />, and MIPS maps from the Spitzer COSMOS (SCOSMOS) program. In the optical regime we used the same maps as in Laigle et al. (2016, see their Table 1 for a summary). From the stacked maps, we extract the median fluxes using Source Extractor (SExtractor; Bertin & Arnouts 1996) for the bands from H to MIPS 24 µm, while from $u^*$ to J, where the stacked flux is at best marginally detected, we choose the Aperture Photometry Tool (Laher 2012) to execute this analysis.

In the FIR regime we compute the median fluxes by performing survival analysis (Isobe & Feigelson 1986), which takes into account the presence of upper limits, on the Herschel counterparts from the aforementioned Super-Deblended catalog by Jin et al. (2018). The choice of not performing the image

<footnote n="15">¹⁵ This is the same UltraVISTA data release as used in the recent COSMOS2020 catalog (Weaver et al. 2022).</footnote>

<footnote n="16">¹⁶ The IRAC maps used in the COSMOS2020 catalog are instead those from Moneti et al. (2022): the difference is minimal and only concerns channels 1 and 2, as shown in Figure 3 of Weaver et al. (2022).</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Median stacked images of the Type-0 subsample in the *u*\* to MIPS 24 µm bands. The red circles indicate the bands where we find a stacked flux with S/N > 2.</figcaption></figure>

stacking for these bands is due to the coarse resolution of the FIR observations and the ensuing difficult extraction of the sources from the maps.

### 3.1.1. Type 0

We perform the same procedure also for the Type-0 subsample alone and we show the stacked images in Figure 3: although these sources do not have any counterpart up to the MIPS 24 µm map, the median stacking reveals in *H*, *Ks*, IRAC-1, IRAC-2, and MIPS 24 µm an emission at S/N > 2. Many of the individual sources belonging to this category have extremely bright radio fluxes, up to 5.2 ± 0.3 µJy, and show peculiar radio morphologies; furthermore, 12 out of 34 of these objects belong to the primary sample, i.e., they have a counterpart in at least one of the FIR/(sub) millimeter bands. These sources could be heavily obscured galaxies and high-redshift candidates. However, SED fitting is inconclusive for these sources, as would be expected given the very poor photometric coverage currently available. Only with more data might we be able to retrieve some information on this subsample, which is particularly suited for deeper follow-up observations with ALMA and JWST.

### 3.2. Median SED Fitting

We derive the median properties of the Types 0+1+2 sources by performing SED fitting on the median stacked photometric data, from *u*\* to 1.4 GHz, described in the previous section. To this purpose, we use the MAGPHYS+PHOTOZ code (da Cunha et al. 2015; Battisti et al. 2019).

MAGPHYS is a SED fitting code based on energy balance. It derives the physical properties of a galaxy through the comparison of the observed photometry to a large number of models: for each of the possible parameters, the code builds a likelihood distribution to choose the one that best fits the observations. MAGPHYS can model a large variety of galaxies, as shown by Battisti et al. (2019). The photo-z version of MAGPHYS, MAGPHYS+PHOTOZ, allows us to compute the redshift along with the other properties of the galaxy. This feature brings numerous advantages: fitting all the properties simultaneously means having more robust uncertainty parameters. One of the main reasons for choosing MAGPHYS in this work is that the code uses all the photometric points up to the radio band, which is particularly suited to a sample like ours where the coverage in the optical and NIR bands is very scarce.

For the stellar component, MAGPHYS uses the simple stellar populations from Bruzual & Charlot (2003), assuming a Chabrier (2003) IMF, with delayed exponentially declining
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Top: median SED of Types 0 + 1 + 2 galaxies. Crimson and empty teal points and arrows represent, respectively, this work and T21 detections and $1\sigma$ upper limits. The best-fit model derived for our stack is the dark gray line for this work, while for T21 it is represented by the light gray dashed line. The SED fitting is performed on the median stacked photometry from the UV to the mid-IR bands and on the median fluxes from the Super-Deblended catalog (Jin et al. 2018) in the FIR regime, as explained in the text. The table shows some of the properties from the SED fitting, i.e., photometric redshift, stellar mass, dust luminosity, and extinction, for both T21 and this work. The errors are computed as the semi-interval between the 16th and the 84th percentiles of the probability distribution functions of each parameter. The quoted error on the redshift includes also a component from bootstrap analysis on the stacked photometry. Bottom: median stacking of the Types 0 + 1 + 2 subsamples in the $u^*$ to MIPS 24 µm bands.</figcaption></figure>

star formation histories combined with random bursts. Dust attenuation is taken into account as described by Charlot & Fall (2000), with the addition of the 2175 Å feature (Battisti et al. 2020), while dust emission modeling takes into account both the components due to stellar birth clouds and the one due to the ambient interstellar medium (ISM) (da Cunha et al. 2008). The radio component model relies on the hypothesis that the IR-radio correlation parameter is $q_{\text{TIR}} = 2.34$. For the metallicity, MAGPHYS uses a prior uniformly distributed between 0.2 and 2.0 solar metallicities. For further details on the code, see Battisti et al. (2020).

From the SED fitting procedure on the stacked photometry of the Types 0 + 1 + 2 subsample, we obtain the median properties of the galaxies in our catalog.

The resulting best-fit SED is shown in Figure 4. The median properties of the sample are consistent with those of ultraluminous infrared galaxies (ULIRGs), since the median

infrared luminosity is as high as $L_{\text{IR}} = (8.48 \pm 0.05) \times 10^{12} L_{\odot}$, corresponding to an SFR $\sim 100 M_{\odot} \text{yr}^{-1}$. This value is broadly consistent with, though slightly higher than, that by T21 ($L_{\text{IR}} = (2.3 \pm 0.5) \times 10^{12} L_{\odot}$). The median photometric redshift$^{17}$ is $z_{\text{med}} = 3.3 \pm 0.2$, which is also consistent with T21 ($z_{\text{med}} = 3.1 \pm 0.2$) and suggests the presence of a relevant number of high-redshift ($z > 4.5$) galaxies in our sample.

In order to investigate the possible presence of obscured AGNs in our sample, we extend our analysis toward the high-energy side of the spectrum. Thanks to the CSTACK online tool, we stacked X-ray data in the soft band ([0.5–2] keV) from the Chandra-LEGACY survey (Civano et al. 2016). After the exclusion of X-ray point sources, assuming $N_{\text{H}} = 10^{22} \text{cm}^{-2}$ and a photon index $\Gamma = 1.8$, the resulting rest-frame luminosity

<footnote>$^{17}$ The error on the redshift is the quadratic sum of the output from MAGPHYS and the result of bootstrap analysis on the stacked photometry.</footnote>
is $L_{2-10\,\text{keV}} \sim 10^{42} \, \text{erg}^{-1}$, which suggests the absence of powerful AGNs in our sample. In fact, this X-ray emission can be attributed to the star formation process alone, since it corresponds to an SFR $\sim 100 \, \text{M}_\odot \, \text{yr}^{-1}$, according to the Kennicutt & Evans (2012) relation, in perfect agreement with the SFR computed through $L_{\text{IR}}$.

The small differences in the median results with respect to T21 are mainly due both to the widening of the sample, which enabled us to have more constraining upper limits, and to the fact that the new sources are overall more disturbed than the ones studied in the previous work.

To test the representativeness of the results obtained with the stacking analysis, we also compared the best-fit SED obtained from the median stacked photometry of the primary sample alone to the median of the best-fit SEDs of the individual galaxies described in the next section. From this test, it emerges that the two composite SEDs are almost coincident, the main difference being that in the former the FIR peak is wider. This is due to the redshift distribution of the sample sources, which artificially enlarges the bell in the Herschel region of the SED. This deviation, however, does not have any dramatic effect on the redshift nor on the recovered median properties.

# 4. Individual Photometry

We perform the analysis of the individual sources as well, where the visual inspection allows us to choose the best technique to process the four different types. This is a step forward with respect to T21: in this work, in fact, we present brand new photometry for the sources with confusion issues in the IRAC bands that takes into account the possible contamination by nearby sources.

## 4.1. Optical and NIR

By selection, our sources are either not detected or extremely faint in the optical bands, and, as such, in these bands the flux densities are given as upper limits. In particular, we set the threshold for considering whether the source is detected or not in a given band at $S/N = 3$.

In particular we use those sources reported by Laigle et al. (2016).

As for the NIR bands, we took advantage of the UltraVISTA DR4 catalog,<fnref n="18" /> which is deeper than the UltraVISTA DR2 on which the COSMOS2015, and hence our selection, is based.

We can distinguish two different cases. First, radio sources that have a counterpart with $S/N > 3$ in the deeper images (128 objects). In this case we took the flux density and associated errors from the existing catalog. For the marginally detected sources we set the error on the flux as the $1\sigma$ upper limit. We note that most of these sources therefore also have a counterpart in COSMOS2020 (Weaver et al. 2022), the most recent version of the multiwavelength photometric catalog in the COSMOS field, which is also based on the UltraVISTA DR4. We briefly discuss this point in Appendix A. We stress that, given that our selection was based on the lack of a counterpart in the COSMOS2015 catalog, we decided not to treat these sources separately.

Second, sources with no counterpart in the deeper images. In this case, we measure the $1\sigma$ upper limit directly on the maps

within a fixed aperture of $2''$ using the Python PHOTUTILS package (Bradley et al. 2020).

## 4.2. MIR

We extract the flux densities in the MIR bands from the IRAC-1, -2, -3, and -4 and MIPS 24 µm maps from the v2.0 mosaics of SPLASH using SEExtractor (SE; Bertin & Arnouts 1996) for all objects, except those belonging to Type-2 class (see Section 2), whose MIR photometry is affected by blending problems due to the poor spatial resolution. For these sources we perform the photometry with specific methods, described in Section 4.2.1. For the objects not suffering from blending, we first make a blind run on the MIR maps with SE on all the sources. This procedure generates a catalog that we match with our sources using a $1''/7$ radius (Smolčić et al. 2017) to find the MIR flux densities. If the MIR counterpart has $S/N < 3$, we place an upper limit equal to the error on the extracted flux itself. If the source has no counterpart in the SE catalog, we measure the upper limits with PHOTUTILS with an aperture size of $3''/6$.

### 4.2.1. The Deblending Procedure

Following a visual inspection, we isolate a subsample of 95 galaxies (Type 2; see Section 2) that suffer from contamination by one or more adjacent sources in the IRAC bands.

For these galaxies we did not obtain satisfying results using SE, even when forcing the aperture photometry at the radio position in double mode and then correcting for the aperture loss. Thus, we decide for a different approach to extract the fluxes.

We use the Python package `c` (Lang et al. 2016), which is particularly efficient in the deblending of sources similar to ours,<fnref n="19" /> as already shown by Enia et al. (2022; see also Weaver et al. 2022).

Starting from defined astronomical sources with certain specified properties (i.e., position, flux, and so on), the code gives estimates in the pixel space of what should be found in the observed image. The input map is the IRAC map centered on the radio coordinates of the RS-NIRdark source. Therefore, the IRAC map should show the MIR counterpart of this source. However, as we said before, sometimes this counterpart is contaminated by the flux from one or more other sources. The contaminants are usually sources belonging to the COSMOS2015 catalog. Thus, we give as prior positions to TRACTOR the radio and COSMOS2015 coordinates and an initial guess on the flux value of all the sources involved (target and contaminants). We set this initial guess as a fraction of the total blended flux extracted with SE. We summarize the input parameters in Table 1, while in Figure 5 we show the results on an IRAC-1 map of a galaxy in our sample.

TRACTOR produces a mask containing the positions and the fluxes given as input. We let TRACTOR be free to move inside the map to find the brightness peaks. Through a Markov Chain Monte Carlo process, the code builds a best-fit model and creates a catalog containing the flux of the deblended sources.

We compared the residual maps obtained after subtracting the blended source+contaminant detection using SE with those obtained after subtracting the deblended source and contaminant using TRACTOR: they are practically identical, meaning

<footnote n="18">¹⁸ https://ultravista.org</footnote>

<footnote n="19">¹⁹ In IRAC-4 and MIPS the deblending procedure failed because of their too coarse resolution, so we used the blended flux by SE in these bands.</footnote>
<table><caption>**Table 1**

Description of the Input Parameters for the Deblending Procedure with TRACTOR, and the Specific Values Used in Our Work</caption><thead><tr><th>Input Parameter</th><th>Description</th><th>This Work</th></tr></thead><tbody><tr><td>Map</td><td>Blended map</td><td>IRAC-1 27 × 27 pixels<sup>2</sup> map</td></tr><tr><td>Positions</td><td>Positions of contaminants and target</td><td>The COSMOS2015 position of the contaminant + the VLA-COSMOS 3 GHz position of the target</td></tr><tr><td>Fluxes</td><td>Best-guess fluxes for the target and the contaminants</td><td>SE blended flux/N;<br/>N = number of contaminants + 1 (= 2 if 1 contaminant)<br/>1.7"</td></tr><tr><td>Point-spread function</td><td>FWHM of the IRAC band considered</td><td>0.1</td></tr><tr><td>Per-pixel image noise</td><td>FWHM of the map</td><td>10<sup>-3</sup></td></tr><tr><td>dlnP</td><td>Precision of the linearized least squared method for the best fit</td><td>Gaussian with σ = 1.7</td></tr><tr><td>Profile</td><td>Gaussian with σ = FWHM</td><td></td></tr></tbody></table>

Note. In Figure 5 we show the input map, the output model, and the residuals relative to the example described in this table.

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Example of the deblending procedure on an IRAC-1 map. The first image is the input map, with the input positions (white circles). The second image is the model produced by TRACTOR, where the white crosses indicate the best-fit positions. The last image is the residual map.</figcaption></figure>

that the two methods recover the same total flux densities in the source+contaminant region. This constitutes additional proof of the reliability of our results. With this method, we successfully deblended 46 of the 52 Type-2 primary sources.

### 4.3. FIR and (sub)millimeter

We extend our photometric range to the FIR thanks to the Herschel Super-Deblended Catalog of Jin et al. (2018). We find a counterpart at the $3\sigma$ significance level in at least one Herschel band for 139 galaxies out of 272. For marginally detected sources we set as $1\sigma$ upper limits the error on the flux, when available. We also perform a match within a radius of 0.8" between our data and the A3COSMOS catalog,<fnref n="20" /> which contains photometric data from the ALMA archive. We find counterparts at $S/N > 3$ in at least one ALMA band for 34 galaxies. Given the extremely heterogeneous nature of the observations used to build the catalog, it is impossible to define upper limits. Therefore, when no counterpart is present in the A3COSMOS catalog in the ALMA bands, we exclude the ALMA point from the photometric SED that we give as input to MAGPHYS.

We recall that the 145 galaxies with at least one $3\sigma$ detection in the FIR/(sub)millimeter range constitute our primary sample (see Section 2), which we use as the basis of our SFRD computation. We perform individual SED fitting analysis for the primary sample only, for which we expect more reliable results than for the whole sample because of the broader photometric coverage from MIR to FIR/(sub)millimeter.

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Histogram showing the percentage of galaxies in our primary sample in a given $S_{500}/S_{350}$ bin. The dashed line indicates the median value.</figcaption></figure>

#### 4.3.1. 500 µm Risers

Given that most of the radio-selected dust-obscured galaxies have detections in the FIR bands (145 out of 272 in our sample), it is possible to exploit the trend of the flux between 250 and 500 µm to infer a lower limit for the redshift estimate. In fact, it has been shown that galaxies with $S_{250} < S_{350} < S_{500}$ are likely to sit at $z > 4$ (Donevski et al. 2018; Duivenvoorden et al. 2018; Greenslade et al. 2020). We find 22 of these so-called 500 µm risers in our sample (see Figure 6). The SPIRE diagram in Figure 7 displays the colors of the 500 µm risers in our sample, which are compatible with other samples of DSFGs at $z > 4$ from the literature. In our whole catalog, only one galaxy is above the threshold of 30 mJy defined by Donevski et al. (2018), with $S_{500} \sim 34$ mJy. All the others are well below this limit, suggesting a negligible impact of lensing in this analysis, as suggested by Lapi et al. (2014), Negrello et al. (2017), and Donevski et al. (2018).

According to the models (Blain et al. 2004), at $S_{500} \sim 100$ mJy the number count of nonlensed galaxies drops, due to the intrinsically steep luminosity function at $z > 1.5$; therefore, this flux limit is a good threshold for identifying magnified sources. The 500 µm risers in our sample have relatively low fluxes compared to those studied in the works mentioned above. This feature might be due not only to the lack of lensed sources but also to the intrinsic nature of these galaxies, which are likely normal dust-obscured SFGs at $z > 4$, rather than extreme objects.

### 4.4. Radio

Our analysis is based on radio detections at 3 GHz from COSMOS VLA 3 GHz Large Project, therefore all our sources

<footnote n="20">²⁰ Version v.20200310.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 7.** SPIRE color–color diagram of our 500 µm risers, overlaid with redshift tracks of Arp 220 (Rangwala et al. 2011) and Cosmic Eyelash (Swinbank et al. 2010) and other samples of DSFGs as a reference (Oteo et al. 2017). The color scale indicates the photometric redshifts as computed by MAGPHYS.</figcaption></figure>

have at least one radio photometric point. In order to improve our knowledge about the radio band, we also use 1.4 GHz data included in the Herschel Super-Deblended Catalog, taken from Schinnerer et al. (2010).

### 4.5. X-Ray

Finally, we match our sample with the Chandra catalog of Civano et al. (2016): 11 sources have a counterpart within 5". However, the low resolution in these bands makes it challenging to understand whether the X-ray emission belongs to our radio target or a nearby COSMOS2015 source. After an accurate analysis, we conclude that for five of these sources, the X-ray is more likely related to the nearby NIR-bright object. For another object the association remains ambiguous, while for the remaining six the X-ray emission belongs most likely to our targets: two show evidence of being heavily obscured, close to the Compton Thick regime, one is an AGN with $L_X \sim 10^{45}$ erg$^{-1}$, and three are obscured AGNs ($n_H \gtrsim 10^{23}$ cm$^{-2}$). We exclude these galaxies from the computation of the SFRD and will defer their study to future works.

## 5. Results from the Individual SED Fitting

From the SED fitting of the individual sources in the primary sample, we identify some of the main properties of the objects in our sample and estimate the possible presence of obscured AGNs.

First of all, we show the redshift distribution in Figure 8. The typical error on the photometric redshift computed via SED fitting is ~10%. The usage of the photometry (and upper limits) over the full available wavelength range improves the reliability of our results for the redshift estimation. As a test, we perform an SED fitting run using only the FIR photometry, as done in other works (e.g., Jin et al. 2018): we find that the mean relative uncertainty is about seven times higher.

In the following, we focus on the sources with an estimated $z_{phot} < 5$, given the uncertainties in the dust properties at higher redshifts.

The dust luminosity $L_{dust}$ has a median value over the sample of $2.8 \times 10^{12} L_{\odot}$ and it varies between $2.9 \times 10^{10} L_{\odot}$ and $1.3 \times 10^{13} L_{\odot}$. Thus, we conclude that most of the objects in our catalog are ULIRGs. We also find six sources with $L_{dust} > 10^{13} L_{\odot}$, i.e., hyperluminous infrared galaxies. These very high values for the infrared luminosity might indicate that many of these sources are going through a starburst phase. In fact, using the Kennicutt relation (Kennicutt & Evans 2012) we derive a median SFR of $\sim 4.2 \times 10^2 M_{\odot} yr^{-1}$.

From the radio flux at 3 GHz, we derive the luminosity at 1.4 GHz for all primary sources, assuming a spectral index $\alpha = -0.7$. This is a fair assumption, since it corresponds to the peak of the spectral index distribution computed for the subsample of sources that have a counterpart at 1.4 GHz (see also Novak et al. 2017). Through the relation between infrared luminosity $L_{dust}$ and radio luminosity $L_{1.4}$, we are then able to explore the possible presence of AGNs in our sample by looking at the $q_{TIR}$ parameter. Its median value $q_{TIR} \sim 2.27$ is consistent with the results from T21 and Novak et al. (2017; see also Algera et al. 2020), suggesting no evidence for strong radio AGN emission in our sample (Delhaize et al. 2017).

Furthermore, one of the most interesting parameters computed by MAGPHYS is the fraction $f_{\mu}$ of dust luminosity due to the diffused component of the ISM, compared to the one due to the birth clouds ($1 - f_{\mu}$). The younger the galaxy is, the higher the contribution of the birth clouds, since they should host large star-forming regions with newborn stars. The parameter $f_{\mu}$ strongly correlates with the SED shape for dusty galaxies: larger values of $f_{\mu}$ correspond to lower dust temperatures and SFR (Martis et al. 2019). We can decompose the FIR luminosity into a hot ($T \sim 60–100$ K) dust component and a warm one ($T \sim 20$ K) and the ratio of the luminosities of these two is an efficient indicator of the dominant power in the galaxy (Kirkpatrick et al. 2012). The hot dust arises from the birth clouds or possible heating from an AGN; the warm dust is associated with the diffuse ISM and is mostly composed of larger dust grains (Kirkpatrick et al. 2015).
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Distribution of the photometric redshifts of the individual sources in this work (crimson) compared with the results of T21 (teal blue) for their similarly defined primary samples.</figcaption></figure>

In high-redshift dusty galaxies, we expect a lower $f_μ$ than in the local ones since a large fraction of large grains did not have time to form at early times, and these objects are usually efficient star-formers, hosting a large number of star-forming regions. On the other hand, for galaxies without an AGN contribution, the warmer component should not be dominant. We find a median value of $f_μ \sim 0.56$, in agreement with those reported for high-redshift dusty galaxies by Kirkpatrick et al. (2015) and Martis et al. (2019). 21% of our galaxies have $f_μ < 0.3$: this could be a hint of an extreme starburst phase or of possible AGN activity.

In conclusion, this is further confirmation that the galaxies in our primary sample are mostly young dusty star-forming galaxies. We also stress that for 101 out of 145 objects in our primary sample, we have photometric points in the Rayleigh-Jeans regime of the dust graybody, i.e., in the (sub)millimeter bands, and hence a direct constraint on $f_μ$. However, further observations with interferometers like ALMA and NOEMA are necessary to improve our knowledge of this topic.

# 6. Star Formation Rate Density

## 6.1. Computation of the SFRD

To compute the SFRD for the primary sample, we split it into three redshift bins: $0 < z < 2$, $2 < z < 3$, $3 < z < 5$. The third bin contains about 85% of the objects, and they are also the ones we are more interested in for our final purpose.

We use the $1/V_{max}$ method of Schmidt (1968). We find the maximum observable volume for each galaxy as

$$V_{\text{max}} = \sum_{z=z_{\text{min}}}^{z_{\text{max}}} [V(z + \Delta z) - V(z)] \times \frac{C_A}{C_I}, \quad (1)$$

where $\Delta z = 0.005$, $z_{min}$ is the lower boundary of the redshift bin, and $z_{max}$ is the minimum value between the redshift upper boundary and the highest redshift at which we can detect the source, given the sensitivity of the survey.

The $C_A$ and $C_I$ constants take into account the incompleteness due to our selection criteria. We recall that the 272 galaxies analyzed in this work were selected from an initial sample of 476 RS-NIRdark sources only on a geometrical basis; therefore there

is reason to think that the selected sources follow a distribution different than that of the whole sample. We further divided the 272 galaxies into primary and secondary samples, depending on their photometric coverage in the FIR/(sub)millimeter regime, and derived individual redshifts and physical properties only for the former. The galaxies in the secondary sample might indeed be less star-forming, or at even higher redshift, than the primary subsample, and might provide a lower contribution to the SFRD. Following T21, we define the $C_I$ parameter (correction for incompleteness) as

$$C_I = \frac{(N_{\text{primary}} + f \times N_{\text{secondary}})}{N_{\text{primary}}} \times \frac{N_{\text{tot}}}{N_{\text{primary}} + N_{\text{secondary}}} \quad (2)$$

where $N_{primary} = 145$ (i.e., the primary subsample), $N_{secondary} = 127$, $N_{tot} = 476$, and $f$ is the ratio between the mean SFR$_IR$ of the primary and secondary subsamples.

To derive the $C_I$ parameter we considered the three scenarios outlined in T21. In Case 1, $f=1$ gives an upper boundary to the SFRD under the assumption that the primary sample is representative of the whole catalog. In Case 3 we fix $f$ to the ratio between the mean FIR fluxes of the two subsamples, derived by stacking the images of our sources in five Herschel bands (from 100 to 500 µm), SCUBA 850 µm, and AzTEC 1.1 mm, following Béthermin et al. (2015). In this work we find $f=0.56$, and we consider the resulting values of the SFRD as a lower boundary.

The other correction factor, $C_A$, considers the effective area (see Equation (2) in T21).

To find the SFRD in each redshift bin, we weight the SFR of every galaxy by its $V_{max}$ and sum together all those belonging to the same redshift range. We performed a bootstrap analysis to find the values of the SFRD and the related errors, taking into account the likelihood distribution of the redshift given by the SED fitting code. We extracted random values of the redshift for each galaxy in the primary sample from this likelihood distribution. Thus, we inferred the value of the infrared luminosity for each of these redshifts. In this way, we obtain a distribution of infrared luminosities and redshifts for every object in the sample, and we can use them to compute a distribution of SFRDs. In each bin, the final value of the SFRD is the median of the distribution. We compute the errors from the 16th and 84th percentile.

## 6.2. Results

The SFRD as a function of redshift is plotted in Figure 9 and the values in each bin are reported in Table 2.

The values at $3 < z < 5$ are consistent with those found by T21 for a subsample of our galaxies, by Enia et al. (2022) for a sample of radio-selected H-dark galaxies in the GOODS-N field, and by Gruppioni et al. (2020) for a sample of (sub) millimeter selected NIR-dark galaxies in the ALPINE fields, as reported in Figure 9.

In particular, T21 find an SFRD equal to $(7.1 \pm 1.7) \times 10^{-3} M_{\odot} \text{yr}^{-1} \text{Mpc}^{-3}$ in the third redshift bin. Gruppioni et al. (2020) computed the SFRD up to $z \sim 6$ for 56 sources detected in the ALMA Band 7 in the ALPINE fields. They find an SFRD of $(1.5 \pm 0.9) \times 10^{-2} M_{\odot} \text{yr}^{-1} \text{Mpc}^{-3}$ between $z = 3.5$ and $z = 4.5$ and of $(0.9 \pm 0.7) \times 10^{-2} M_{\odot} \text{yr}^{-1} \text{Mpc}^{-3}$ between $z = 4.5$ and $z = 6.0$. Enia et al. (2022) analyzed 17 H-dark galaxies in the GOODS-N field using the same procedure adopted in this paper and in T21 and find the SFRD at $z \sim 3$ to
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Star formation rate density as a function of redshift. The crimson patterned and the violet rectangles show the confidence intervals for the results of this work and T21, where the lower boundary is the value obtained for Case 3 and the upper boundary the value obtained for Case 1. We indicate with green, pink, and turquoise the values found by Gruppioni et al. (2020) HST-dark, Gruppioni et al. (2020) HST-bright ALMA Band-7 ALPINE detections, and Yamaguchi et al. (2019) respectively. The violet arrow is the result of Enia et al. (2022) for H-dark galaxies. Red crosses are the results of Novak et al. (2017) for VLA-COSMOS 3 GHz Large Project sources with an optical/NIR counterpart, while light green diamonds are the results of Williams et al. (2019), purple diamonds are from Wang et al. (2019), blue diamonds from Fudamoto et al. (2021), and gray circles from Barrufet et al. (2023). The best-fit curve of Madau & Dickinson (2014), scaled to a Chabrier (2003) IMF, is displayed as a dashed black line.</figcaption></figure>

<table><caption>**Table 2**

Values for the SFRD Found Using the Correction of Case 1 (Column 2) and Case 3 (Column 3) Described in Section 6.1</caption><thead><tr><th>Bin</th><th>SFRD (f = 0.56)<br>(M<sub>☉</sub> Mpc<sup>-3</sup> yr<sup>-1</sup>)</th><th>SFRD (f = 1)<br>(M<sub>☉</sub> Mpc<sup>-3</sup> yr<sup>-1</sup>)</th></tr></thead><tbody><tr><td>0.0 &lt; z<sub>phot</sub> &lt; 2.0</td><td>(0.6 ± 0.1) × 10<sup>-3</sup></td><td>(0.8 ± 0.1) × 10<sup>-3</sup></td></tr><tr><td>2.0 &lt; z<sub>phot</sub> &lt; 3.0</td><td>(2.4 ± 0.4) × 10<sup>-3</sup></td><td>(3.2 ± 0.5) × 10<sup>-3</sup></td></tr><tr><td>3.0 &lt; z<sub>phot</sub> &lt; 5.0</td><td>(5.0 ± 0.7) × 10<sup>-3</sup></td><td>(6.2 ± 0.6) × 10<sup>-3</sup></td></tr></tbody></table>

be ~4.5 × 10⁻³ M☉ yr⁻¹ Mpc⁻³, consistent with what we find in this work at the same redshift.

We stress that our estimate of the RS-NIRdark contribution to the SFRD is only a lower limit (similarly to the cases presented in T21; Gruppioni et al. 2020; Enia et al. 2022) because we are not extrapolating to lower fluxes than our detection limit.

Comparing the SFRD from our sample of RS-NIRdark galaxies to the ones computed through UV-selected galaxies (Madau & Dickinson 2014; Bouwens et al. 2020), we find that the contribution to the SFRD of our galaxies at 3 < z < 5 is ~10%–25% of the SFRD of UV-bright ones at the same redshifts, e.g., Bouwens et al. (2020), who studied the SFRD in a sample of 1362 UV-selected Lyman-break galaxies in the ALMA Spectroscopic Survey in the Hubble Ultra-Deep Field (ASPECS) between redshifts 1.5 and 10. They also include the contribution to SFRD from ULIRGs. The best-fit trend of the

SFRD, computed by Madau & Dickinson (2014), is also based mainly on UV-selected sources. The SFRD trend based on UV data is steeper than that found through radio data, e.g., by Novak et al. (2017), who compute the SFRD for the full VLA-COSMOS 3 GHz Large Project (Smolčić et al. 2017) sample, complementary to the sample presented in this work. This means that UV-based studies miss an important part of the SFRD at high redshift due to dusty galaxies. Moreover, as emerges from this analysis, the SFRD related to the more dusty objects selected in the radio band, like our RS-NIRdarks, seems important at high redshift. Furthermore, we might still be missing an important fraction of obscured galaxies because of technological limitations that will be overcome thanks to new generation of telescopes, such as JWST and the Square Kilometer Array.

The value of the number density at 3 < z < 5, (6.7 ± 0.9) × 10⁻⁶ Mpc⁻³, is overall compatible with those found by T21, Riechers et al. (2020), Enia et al. (2022), and Wang et al. (2019). Semianalytical models (Henriques et al. 2015) and hydrodynamical simulations (Snyder et al. 2017) do not predict the presence of infrared ultraluminous dust-obscured galaxies at high redshift, but our results are in agreement with the number density computed by Straatman et al. (2014), Schreiber et al. (2018a), and Girelli et al. (2019) for massive quiescent galaxies at 3 < z < 4. This result supports the hypothesis of these galaxies being the most likely missing progenitors of massive quiescent galaxies, according to the in situ scenario (Lapi et al. 2018). Future observations will be fundamental to truly understanding the role of these objects in the galaxy formation and evolution scenarios.

## 7. Summary

In this work, we estimate the contribution to the SFRD made by RS-NIRdark galaxies selected from the VLA-COSMOS 3 GHz Large Project. We find that:

1. These objects are likely highly infrared luminous (median $L_{dust} \sim 4 \times 10^{12} L_{☉}$), dust-obscured galaxies at a median photometric redshift of $z \sim 3.2$: they are possibly young high-redshift DSFGs, as also suggested by the median value of $f_{μ} \sim 0.56$.

2. 76 of the 272 selected sources have a photometric redshift > 3.

3. Through the SPIRE color–color diagram, we have further proof, independent of the SED fitting procedure, that there are high-redshift candidates in our sample.

4. The contribution of RS-NIRdark galaxies to the SFRD at 3 < z < 5 is at least ~40% of that due to UV-selected sources, which means that the role of the dark galaxies might be crucial to understanding the still widely unknown evolution of SFRD at $z > 3$.

Our results remark on the importance of the obscured galaxies in galaxy formation and evolution and the need for a better understanding of their role in the broader cosmological framework. This work is a step toward a less biased view of the nature of these kinds of objects. The inclusion of more contaminated sources than in previous studies (Talia et al. 2021) allowed us to start building an improved photometric catalog and to get more information from our data on the physical properties of these galaxies.
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Redshift distribution of the 76 sources in the primary sample with a counterpart in the COSMOS2020 catalog: in green we show the redshifts computed in this work; in red the redshifts reported in the COSMOS2020 catalog.</figcaption></figure>

## Acknowledgments

We acknowledge the use of Python (v.3.8) libraries in the analysis. This work is partly based on data products from observations made with ESO Telescopes at the La Silla Paranal Observatory under ESO program ID 179.A-2005 and on data products produced by CALET and the Cambridge Astronomy Survey Unit on behalf of the UltraVISTA consortium. M.T., A.C., A.L., M.M., A.E., F.G., C.G., and F.P. acknowledge the support from grant PRIN MIUR 2017 20173ML3WW_001. M.B. acknowledges A. Battisti for adding the A3COSMOS filters to the SED fitting code MAGPHYS+PHOTOZ. M.B. acknowledges A. Traina and N. Borghi for all the concrete support given from the very beginning of this work and F. Gabrielli for useful inspiring discussions.

## Appendix A
### Comparison with COSMOS2020

In this work we based our selection of RS-NIRdark galaxies on the lack of a counterpart in the NIR-based COSMOS2015 catalog (Laigle et al. 2016). However, as previously stated, in our analysis we then took advantage of the UltraVISTA DR4 maps, deeper than those used to build the COSMOS2015 catalog. During the writing of this paper, an updated version (COSMOS2020) of the multiwavelength catalog in the COSMOS field has been published (Weaver et al. 2022), which is based on UltraVISTA DR4 maps in the NIR bands. We matched our primary sample to the COSMOS2020 within a 0.′′8 radius and we found 76 counterparts. In our work we did not use the results from COSMOS2020 for these sources, for consistency with the rest of our analysis: in fact the rest of our primary sample is not present in the COSMOS2020 catalog. Also, we decided not to treat these sources separately.

The biggest difference in the photometric depth between our catalog and the COSMOS2020 is in the MIR bands, where IRAC-1 and IRAC-2 maps (Moneti et al. 2022) are slightly

<table><caption>**Table 3**
Comparison between the Depth of the IRAC Maps Used in This Work and Those Used in COSMOS2020</caption><thead><tr><th></th><th>IRAC-1</th><th>IRAC-2</th><th>IRAC-3</th><th>IRAC-4</th></tr></thead><tbody><tr><td>This work; COSMOS2015</td><td>25.5</td><td>25.5</td><td>23.0</td><td>22.9</td></tr><tr><td>COSMOS2020</td><td>25.7</td><td>25.6</td><td>22.6</td><td>22.5</td></tr></tbody></table>

deeper than the ones we used, as summarized in Table 3. We checked the MIR fluxes of the sources that we have in common with the COSMOS2020 catalog and verified that using the deeper IRAC maps would not change our results.

We find also that, for the sources in common, the redshift distributions from our own analysis and from the COSMOS2020 catalog (EAZY version, Brammer et al. 2008) are in broad agreement (Figure 10). We stress that in our work we took advantage of the full photometric range up to radio frequencies to derive the redshifts.

## Appendix B
### SED Fitting without 350 and 500 µm

In this appendix, we show the results when excluding from the SED fitting of the sources with A3COSMOS counterparts the 350 and 500 µm data points, which could be less reliable than the others because of their coarse resolution, especially with respect to the ALMA data available for this particular subsample. We show an example in Figure 11 and we compare the redshifts and reduced $\chi^2$ computed in the two ways in Figure 12. On average, the quality of the SED fitting stays almost the same: the median reduced $\chi^2$ for the whole photometry approach is 2.69, while for the reduced photometry it amounts to 2.87. Also, the values of the redshifts stay almost the same: the median redshift in the former case is 3.5, while in the latter it is 3.3. Moreover, the median relative error is about 0.09 for both approaches.
<figure><img src="image_14.png" /><figcaption>**Figure 11.** The best-fit SED and the probability distribution functions of the derived parameters of an example galaxy in our sample performed using the whole available photometry (left) and without the 350 and 500 µm (right).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Left: comparison between the photometric redshifts computed with and without the 350 and 500 µm photometric data. Right: distribution of the reduced χ² of the SED fitting for the sources in the A3COSMOS sample obtained using the whole photometry (crimson) and without the 350 and 500 µm.</figcaption></figure>

## ORCID iDs

Meriem Behiri https://orcid.org/0000-0002-6444-8547

Margherita Talia https://orcid.org/0000-0003-4352-2063

Andrea Cimatti https://orcid.org/0000-0002-4409-5633

Andrea Lapi https://orcid.org/0000-0002-4882-1735

Marcella Massardi https://orcid.org/0000-0002-0375-8330

Andrea Enia https://orcid.org/0000-0002-0200-2857

Cristian Vignali https://orcid.org/0000-0002-8853-9611

Matthieu Bethermin https://orcid.org/0000-0002-3915-2015

Andreas Faisst https://orcid.org/0000-0002-9382-9832

Fabrizio Gentile https://orcid.org/0000-0002-8008-9871

Marika Giulietti https://orcid.org/0000-0002-1847-4496

Carlotta Gruppioni https://orcid.org/0000-0002-5836-4056

Francesca Pozzi https://orcid.org/0000-0002-7412-647X

Vernesa Smolčić https://orcid.org/0000-0002-3893-8614

Gianni Zamorani https://orcid.org/0000-0002-2318-301X

## References

Adams, N. J., Conselice, C. J., Ferreira, L., et al. 2023, MNRAS, 518, 4755

Algera, H. S. B., Smail, I., Dudzevičiūtė, U., et al. 2020, ApJ, 903, 138

Aravena, M., Decarli, R., Walter, F., et al. 2016, ApJ, 833, 68

Atek, H., Shuntov, M., Furtak, L. J., et al. 2023, MNRAS, 519, 1201

Barrufet, L., Oesch, P. A., Weibel, A., et al. 2023, MNRAS, 522, 449

Battisti, A. J., da Cunha, E., Grasha, K., et al. 2019, ApJ, 882, 61

Battisti, A. J., da Cunha, E., Shivaei, I., & D., C. 2020, ApJ, 888, 108

Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393

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

Blain, A., Alexander, D., Armus, L., et al. 2004, Deep Infrared Imaging Around Ultraluminous High-redshift Galaxies, 3473

Boco, L., Lapi, A., & Danese, L. 2020, ApJ, 891, 94

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

Bradley, L., Sipőcz, B., Robitaille, T., et al. 2020, astropy/photutils: v1.0.0, Zenodo, doi:10.5281/zenodo.4044744

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

Brisbin, D., Miettinen, O., Aravena, M., et al. 2017, A&A, 608, A15

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

Casey, C. M., Narayanan, D., & Cooray, A. 2014, PhR, 541, 45

Castellano, M., Fontana, A., Treu, T., et al. 2022, ApJL, 938, L15

Chabrier, G. 2003, PASP, 115, 763

Chapman, S. C., Lewis, G. F., Scott, D., Borys, C., & Richards, E. 2002, ApJ, 570, 557

Chapman, S. C., Richards, E. A., Lewis, G. F., Wilson, G., & Barger, A. J. 2001, ApJL, 548, L147

Chapman, S. C., Smail, I., Blain, A. W., & Ivison, R. J. 2004, ApJ, 614, 671

Charlot, S., & Fall, S. M. 2000, ApJ, 539, 718

Cimatti, A. 2008, Natur, 451, 253

Cimatti, A., Cassata, P., Pozzetti, L., et al. 2008, A&A, 482, 21

Cimatti, A., Daddi, E., Renzini, A., et al. 2004, Natur, 430, 184

Civano, F., Marchesi, S., Comastri, A., et al. 2016, ApJ, 819, 62

da Cunha, E., Charlot, S., & Elbaz, D. 2008, MNRAS, 388, 1595

da Cunha, E., Walter, F., Smail, I. R., et al. 2015, ApJ, 806, 110

Dannerbauer, H., Lehnert, M. D., Lutz, D., et al. 2004, ApJ, 606, 664

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

Donevski, D., Buat, V., Boone, F., et al. 2018, A&A, 614, A33

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

Duivenvoorden, S., Oliver, S., Scudder, J. M., et al. 2018, MNRAS, 477, 1099

Dunlop, J. S., McLure, R. J., Biggs, A. D., et al. 2017, MNRAS, 466, 861

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

Faisst, A. L., Capak, P. L., Yan, L., et al. 2017, ApJ, 847, 21

Finkelstein, S. L., Bagley, M. B., Arrabal Haro, P., et al. 2022, ApJL, 940, L55

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

Frayer, D. T., Chapman, S. C., Yan, L., et al. 2004, ApJS, 154, 137

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

Girelli, G., Bolzonella, M., & Cimatti, A. 2019, A&A, 632, A80
González-López, J., Bauer, F. E., Romero-Cañizales, C., et al. 2017, A&A, 597, A41

Greenslade, J., Clements, D. L., Petitpas, G., et al. 2020, MNRAS, 496, 2315

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

Gruppioni, C., Pozzi, F., Rodighiero, G., et al. 2013, MNRAS, 436, 2875

Hatsukade, B., Kohno, K., Yamaguchi, Y., et al. 2020, in IAU Symp. 352, Uncovering Early Galaxy Evolution in the ALMA and JWST Era, ed. E. da Cunha (Cambridge: Cambridge Univ. Press), 239

Hayward, C. C., Kereš, D., Jonsson, P., et al. 2011, ApJ, 743, 159

Henriques, B. M. B., White, S. D. M., Thomas, P. A., et al. 2015, MNRAS, 451, 2663

Hughes, D. H., Serjeant, S., Dunlop, J., et al. 1998, Natur, 394, 241

Isobe, T., & Feigelson, E. D. 1986, BICDS, 31, 209

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

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

Kirkpatrick, A., Pope, A., Alexander, D. M., et al. 2012, ApJ, 759, 139

Kirkpatrick, A., Pope, A., Sajina, A., et al. 2015, ApJ, 814, 9

Laher, R., 2012 APT: Aperture Photometry Tool, Astrophysics Source Code Library, ascl:1208.003

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

Lang, D., Hogg, D. W., & Mykytyn, D., 2016 The Tractor: Probabilistic Astronomical Source Detection and Measurement, Astrophysics Source Code Library, ascl:1604.008

Lapi, A., Pantoni, L., Zanisi, L., et al. 2018, ApJ, 857, 22

Lapi, A., Raimundo, S., Aversa, R., et al. 2014, ApJ, 782, 69

Lilly, S. J., Le Fevre, O., Hammer, F., et al. 1996, in IAU Symp. 171, New Light on Galaxy Evolution, ed. R. Bender & R. L. Davies (Dordrecht: Springer), 209

Liu, D., Lang, P., Magnelli, B., et al. 2019, ApJS, 244, 40

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

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

Madau, P., Ferguson, H. C., Dickinson, M. E., et al. 1996, MNRAS, 283, 1388

Magnelli, B., Ivison, R. J., Lutz, D., et al. 2015, A&A, 573, A45

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

Mancuso, C., Lapi, A., Shi, J., et al. 2016a, ApJ, 833, 152

Mancuso, C., Lapi, A., Shi, J., et al. 2016b, ApJ, 823, 128

Martis, N. S., Marchesini, D. M., Muzzin, A., et al. 2019, ApJ, 882, 65

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

Moneti, A., McCracken, H. J., Shuntov, M., et al. 2022, A&A, 658, A126

Naidu, R. P., Oesch, P. A., Setton, D. J., et al. 2022, arXiv:2208.02794

Negrello, M., Amber, S., Amvrosiadis, A., et al. 2017, MNRAS, 465, 3558

Novak, M., Smolčić, V., Delhaize, J., et al. 2017, A&A, 602, A5

Oteo, I., Zhang, Z.-Y., Yang, C., et al. 2017, ApJ, 850, 170

Rangwala, N., Maloney, P. R., Glenn, J., et al. 2011, ApJ, 743, 94

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

Rowan-Robinson, M., Oliver, S., Wang, L., et al. 2016, MNRAS, 461, 1100

Schinnerer, E., Sargent, M. T., Bondi, M., et al. 2010, ApJS, 188, 384

Schmidt, M. 1968, ApJ, 151, 393

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

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

Scoville, N., Aussel, H., Sheth, K., et al. 2014, ApJ, 783, 84

Scoville, N., Lee, N., Vanden Bout, P., et al. 2017, ApJ, 837, 150

Scoville, N., Sheth, K., Aussel, H., et al. 2016, ApJ, 820, 83

Simpson, J. M., Smail, I., Dudzevičiūtė, U., et al. 2020, MNRAS, 495, 3409

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

Simpson, J. M., Swinbank, A. M., Smail, I., et al. 2014, ApJ, 788, 125

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

Smail, I., Ivison, R. J., & Blain, A. W. 1997, ApJL, 490, L5

Smolčić, V., Delvecchio, I., Zamorani, G., et al. 2017, A&A, 602, A2

Snyder, G. F., Lotz, J. M., Rodriguez-Gomez, V., et al. 2017, MNRAS, 468, 207

Steinhardt, C. L., Speagle, J. S., Capak, P., et al. 2014, ApJL, 791, L25

Straatman, C. M. S., Labbé, I., Spitler, L. R., et al. 2014, ApJL, 783, L14

Strandet, M. L., Weiss, A., Vieira, J. D., et al. 2016, ApJ, 822, 80

Swinbank, A. M., Smail, I., Longmore, S., et al. 2010, Natur, 464, 733

Szalay, A. S., Connolly, A. J., & Szokoly, G. P. 1999, AJ, 117, 68

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

Walter, F., Banados, E., Venemans, B., et al. 2016, AAS Meeting, 227, 243.39

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

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

Weaver, J. R., Kauffmann, O. B., Ilbert, O., et al. 2022, ApJS, 258, 11

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

Yan, H., Ma, Z., Ling, C., et al. 2023, ApJL, 942, L9

Zavala, J. A., Casey, C. M., Manning, S. M., et al. 2021, ApJ, 909, 165