# Illuminating the dark side of cosmic star formation two billion years after the Big Bang

Margherita Talia, Andrea Cimatti, Marika Giulietti, Gianni Zamorani, Matthieu Bethermin, Andreas Faisst, Olivier Le Fèvre, Vernesa Smolčić

► To cite this version:

Margherita Talia, Andrea Cimatti, Marika Giulietti, Gianni Zamorani, Matthieu Bethermin, et al. Illuminating the dark side of cosmic star formation two billion years after the Big Bang. The Astrophysical Journal, 2021, 909 (1), pp.23. 10.3847/1538-4357/abd6e3. hal-03022738

HAL Id: hal-03022738

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

Submitted on 2 Aug 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Illuminating the Dark Side of Cosmic Star Formation Two Billion Years after the Big Bang

Margherita Talia¹,² ID, Andrea Cimatti³,⁴ ID, Marika Giulietti³,⁵ ID, Gianni Zamorani² ID, Matthieu Bethermin⁶ ID, Andreas Faisst⁷ ID,
Olivier Le Fèvre⁶ ID, and Vernesa Smolčić⁸ ID

¹ University of Bologna—Department of Physics and Astronomy “Augusto Righi” (DIFA), Via Gobetti 93/2, I-40129 Bologna, Italy; margherita.talia2@unibo.it

² INAF—Osservatorio di Astrofisica e Scienza dello Spazio, Via Gobetti 93/3, I-40129 Bologna, Italy

³ University of Bologna—Department of Physics and Astronomy, Via Gobetti 93/2, I-40129 Bologna, Italy

⁴ INAF—Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, I-50125 Firenze, Italy

⁵ SISSA, Via Bonomea 265, I-34136 Trieste, Italy

⁶ Aix Marseille Univ, CNRS, LAM, Laboratoire d'Astrophysique de Marseille, Marseille, France

⁷ IPAC, M/C 314-6, California Institute of Technology, 1200 East California Boulevard, Pasadena, CA 91125, USA

⁸ Department of Physics, Faculty of Science, University of Zagreb Bijenička cesta 32, 10000 Zagreb, Croatia
Received 2020 August 7; revised 2020 December 22; accepted 2020 December 24; published 2021 March 2

## Abstract

How and when did galaxies form and assemble their stars and stellar mass? The answer to these questions, so crucial to astrophysics and cosmology, requires the full reconstruction of the so-called cosmic star formation rate density (SFRD), i.e., the evolution of the average star formation rate per unit volume of the universe. While the SFRD has been reliably traced back to 10–11 billion years ago, its evolution is still poorly constrained at earlier cosmic epochs, and its estimate is mainly based on galaxies luminous in the ultraviolet and with low obscuration by dust. This limited knowledge is largely due to the lack of an unbiased census of all types of star-forming galaxies in the early universe. We present a new approach to finding dust-obscured star-forming galaxies based on their emission at radio wavelengths coupled with the lack of optical counterparts. Here, we present a sample of 197 galaxies selected with this method. These systems were missed by previous surveys at optical and near-infrared wavelengths, and 22 of them are at very high redshift (i.e., z > 4.5). The contribution of these elusive systems to the SFRD is substantial and can be as high as 40% of the previously known SFRD based on UV-luminous galaxies. The mere existence of such heavily obscured galaxies in the first two billion years after the Big Bang opens new avenues to investigate the early phases of galaxy formation and evolution, and to understand the links between these systems and the massive galaxies that ceased their star formation at later cosmic times.

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

## 1. Introduction

How efficiently did gas transform into stars as a function of cosmic time? The answer to this key question requires us to reconstruct the cosmic SFRD to the highest possible redshifts. However, despite the major progress achieved in the last decades in understanding galaxy evolution (Madau & Dickinson 2014), several key questions remain open. The integration of the SFRD (z) over redshift, making appropriate corrections for stellar evolution processes, yields the current stellar mass density $\rho_*(z)$. The results obtained so far (Madau & Dickinson 2014; Oesch et al. 2018) show a rather consistent picture up to $z \approx 3$. Several independent results indicate that the SFRD rapidly increases from $z=0$ to $z \approx 1$, and flattens around $z \approx 2$ (the so-called “cosmic noon”).

However, major uncertainties remain at $z > 3$ (Casey et al. 2014; Magnelli et al. 2019). At these high redshifts, it is still unclear whether the SFRD rapidly declines or remains rather flat (Gruppioni et al. 2013; Rowan-Robinson et al. 2016; Novak et al. 2017; Gruppioni et al. 2020). The poor knowledge of the SFRD evolution at $z > 3$ is primarily due to two main limitations. First, the vast majority of SFRD estimates at these redshifts comes from the observation of Lyman-break galaxies (LBGs) in the rest-frame UV. This makes the results strongly dependent on the adopted dust extinction correction. Moreover, LBGs may not be fully representative of the whole population of the star-forming galaxies (SFGs) existing at these redshifts.

Second, the available surveys in the IR (Spitzer, Herschel; Lutz et al. 2011; Magnelli et al. 2013) and submillimeter/millimeter bands have either a limited sensitivity to galaxies at $z > 3$ and/or are often plagued by large beam sizes, which imply significant source blending. Indeed, since the first deep surveys in the submillimeter regime (850–870 µm) uncovered the presence of submillimeter galaxies at mJy flux levels (SMGs; Smail et al. 1997; Hughes et al. 1998), the coarse angular resolution of single-dish telescopes and the faintness, or sometimes complete lack, of the optical/near-infrared (OPT/NIR) counterparts posed serious challenges to their identification and characterization (e.g., Dannerbauer et al. 2004; Frayer et al. 2004). In this respect, the intense starburst HDF 850.1 is a notable case because of the 15 yr gap between its first discovery (Hughes et al. 1998) and the secure spectroscopic identification (Walter et al. 2012). The ALMA follow-up of increasingly larger samples of SMGs first identified with single-dish observations (e.g., Simpson et al. 2014, 2017; Dudzevičiūtė et al. 2020) proved fundamental to establishing the physical properties of these bright ($S_{870\mu m} > 1–2$ mJy) galaxies, which are typically located around $z \sim 2.5–3.0$ (e.g., Simpson et al. 2014, 2017, 2020; Dudzevičiūtė et al. 2020), are characterized by very high far-infrared (FIR) luminosities and SFRs (several hundreds of $M_\odot$ yr$^{-1}$; Swinbank et al. 2014), and constitute a significant fraction of the SFRD at cosmic noon. ALMA deep fields provide a complementary approach, typically reaching fainter flux limits of $S_{870\mu m} \sim 0.1–1$ mJy,
but they are currently too small to map sufficiently large volumes (e.g., Aravena et al. 2016; Walter et al. 2016; Dunlop et al. 2017; Franco et al. 2018; Hatsukade et al. 2018). Different works report that a fraction of the ALMA-identified SMGs (10%–20%; e.g., Franco et al. 2018; Dudzevičiūtė et al. 2020; Gruppioni et al. 2020) lack an OPT/NIR counterpart. A favored explanation is that these UV- or HST-dark galaxies, as they are often called, are extremely dust-obscured and/or lie at a much higher redshift than the bulk of the SMGs population (>4). For example, Wang et al. (2019) recently reported the results from the ALMA follow-up of a population of optically dark galaxies (H-dropouts), and confirmed a fraction of them to be massive dusty galaxies at high redshift. They concluded that this population constitutes a significant fraction of the SFRD at high redshift (>3), but also left open the question of whether the fraction of currently missed SFGs might be even higher. A handful of extreme SFGs heavily obscured by dust and missed in OPT/NIR surveys has indeed been identified out to very high redshifts (z ~ 5–6; Dannerbauer et al. 2008; Wang et al. 2009; Walter et al. 2012; Riechers et al. 2013, 2017; Marrone et al. 2018; Riechers et al. 2020). However, it is unclear if these seemingly very rare objects are representative of a more vast and elusive population, and whether their high luminosities are partly amplified by gravitational lensing due to intervening matter along the line of sight (Bakx et al. 2020; Ciesla et al. 2020). To answer this question, several efforts are being pursued to uncover dusty systems at very high redshifts (>3–4) with a combination of space- and ground-based data in the FIR to millimeter spectral range (Ivison et al. 2016; Casey et al. 2018; Duivenvoorden et al. 2018; Magnelli et al. 2019; Béthermin et al. 2020; Faisst et al. 2020a; Le Fèvre et al. 2020).

In this work, we present the results of a search for dusty UV-
dark galaxies at $z \gtrsim 3$ selected at radio frequencies, taking
advantage of the depth and area of the VLA–COSMOS 3 GHz
Large Project (Smolčić et al. 2017). SFGs display radio
emission due to a combination of synchrotron radiation emitted
by electrons accelerated in supernova remnants and free–free
continuum from H II regions. As a consequence, the radio
luminosity of SFGs is an indicator of the SFR, provided that
the contamination from AGN emission is negligible. As radio
photons are immune to dust extinction, it is therefore possible
to exploit deep radio surveys to search for dust-obscured,
highly star-forming systems at high redshifts (e.g., Chapman
et al. 2001, 2002, 2004). Moreover, the angular resolution of
interferometric data allows us to localize the OPT/NIR
counterparts much more reliably than for sources selected in
the FIR/millimeter where the beam size and the consequent
blending of different sources can be a severe limitation. The
selection in the radio can be advantageous with respect to FIR/
submillimeter surveys also because it is independent of dust
temperature. It has been suggested that the dust temperature
($T_{\text{dust}}$) correlates with the infrared luminosity, specific star
formation rate, and redshift (Béthermin et al. 2015; Faisst et al.
2017; Schreiber et al. 2018a). For instance, the average $T_{\text{dust}}$
almost doubles from $z \approx 0$ to $z \approx 4$. This would imply that FIR-
to-millimeter surveys could be affected by selection biases,
which depend on $T_{\text{dust}}$, and hence the wavelength at which the
graybody emission peaks; instead, radio emission is free from
these effects. The nature of the $T_{\text{dust}}$–redshift relation is indeed
still debated (e.g., Dudzevičiūtė et al. 2020; Faisst et al. 2020b),
with recent works suggesting that the observed trend could be

mostly due to the relatively high median SFR of the current
sample of dusty SFGs at $z > 5$ (Riechers et al. 2020).

We give magnitudes in the AB photometric system, assume a Chabrier (2003) initial mass function, and use the following cosmological parameters: $\Omega_M = 0.3$, $\Omega_\Lambda = 0.7$, and $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$.

## 2. Radio Selection of Dust-obscured Systems

Our study relies on the data set collected with the VLA–COSMOS 3 GHz Large Project (Smolčić et al. 2017), a survey based on 384-hour observations of the COSMOS field (2 deg²) with the Karl G. Jansky Very Large Array (VLA) interferometer at 3 GHz ($\lambda = 10$ cm). This survey already proved to be deep enough to allow the identification of SFGs at $z > 3$ (Novak et al. 2017), but our aim is to go beyond the results obtained so far and explore the existence of dusty galaxies that were missed in previous studies.

We built our sample of radio-selected, UV-dark galaxies
following these steps:

1. Starting from the VLA source catalog (Smolčić et al. 2017), we selected a parent sample of 8850 objects with flux density at 3 GHz (10 cm) $S_{3\,\text{GHz}} > 12.65 \, \mu\text{Jy beam}^{-1}$, corresponding to $5.5\sigma$, which is the threshold over which the estimated fraction of spurious sources over the entire catalog is only 0.4% (Smolčić et al. 2017).

2. We subsequently excluded all sources that were located inside masked or bad areas of the UltraVISTA footprint (Laigle et al. 2016). The total effective area used in our analysis is $1.38 \, \text{deg}^2$, inside which we counted 5982 sources above the S/N cut.

3. We excluded all multi-component radio sources as identified by Smolčić et al. (2017), because the radio emission in these sources is likely associated with AGN activity, rather than to star formation.

4. Then, we cross-correlated our catalog of radio sources with the COSMOS2015 photometric catalog (Laigle et al. 2016) within a search radius of 0.8" (Smolčić et al. 2017). We identified 476 VLA sources without a COSMOS2015 counterpart.

5. Finally, we restricted our selection to isolated sources, in order to avoid the risk of contamination of the photometry of the galaxies in our sample by nearby, physically unrelated objects, since source blending would make photometry and identification of multiwavelength counterparts uncertain. In particular, we visually inspected the so-called $\chi^2$-image used as the detection map to build the COSMOS2015 catalog (Laigle et al. 2016), which is a map obtained by coadding the $z^+$, Y, J, H, and Ks images, and we excluded from our final sample all sources whose radio $3\sigma$ isophote intersects the NIR $3\sigma$ isophote of nearby objects in the $\chi^2$-image. We did not find strong evidence of gravitational lensing effects.

The final sample analyzed in this work counts 197 galaxies.

## 3. Stacking Analysis and Possible AGN Contribution

Our first approach was to statistically analyze the properties of
these galaxies, because of their extreme faintness at OPT/NIR
wavelengths (down to limiting magnitudes AB = 24.0–24.7 in
the NIR bands). We built the median SED of our total sample of
197 galaxies (Figure 1) by stacking the COSMOS images in
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Top: median SED of the total sample of the 197 radio sources analyzed in this work (red and salmon points). Error bars are plotted for each data point, although in some cases they are smaller than the points. The black line is the best fit from *MAGPHYS*. Detections and 3σ upper limits are derived from stacked images up to 24 µm, while median fluxes from the *Super-deblended* catalog (Jin et al. 2018) are used in the FIR regime. The inset shows the likelihood distribution of the photometric redshift for the total sample of 197 galaxies. We also show the median OPT-to-MIR SED of the secondary subsample of 99 galaxies (black and gray symbols). Bottom: normalized to the peak flux OPT-to-MIR stacked images of the total sample.</figcaption></figure>

each photometric filter, from the optical to 24 µm, at the radio positions. In particular, in the optical regime we used the same maps as in Laigle et al. (2016; see their Table 1 for a summary), in the NIR bands the UltraVISTA DR3 maps, and in the MIR (mid-infrared) regime the SPLASH/IRAC, and MIPS maps (see the next section). In the FIR regime we did not stack directly Herschel and SCUBA maps. Instead, we computed the median of the flux in each photometric band from the *Super-deblended* catalog (Jin et al. 2018) and we employed survival analysis (Feigelson & Nelson 1985) to properly account for the upper limits for the undetected sources.

Quite remarkably, even in the stacked images, which have an effective depth ~14 times higher than the original maps, there are only marginal (slightly more than 2σ) detections in the r+, i+, z++, and J bands, at the current depths, while significant

fluxes emerge only at longer wavelengths, starting from the H band.

We fitted the broadband photometry in the entire wavelength range up to radio with the *MAGPHYS* code (da Cunha et al. 2008; Battisti et al. 2019). *MAGPHYS* is a physically motivated model package that self-consistently models the OPT-to-radio emission of a galaxy. In particular, the emission from stellar populations in galaxies is computed using the Bruzual & Charlot (2003) models, with delayed exponentially declining star formation histories, together with random bursts superimposed onto the continuous model, with a range of ages and exponential timescale values. The effects of dust attenuation are included as prescribed by Charlot & Fall (2000). The models are uniformly distributed in metallicity between 0.02 and 2 times solar. The emission from dust accounts both for the
<table><caption>**Table 1**
Summary of the Median Physical Properties of the Total Sample of 197 Galaxies Analyzed in This Work</caption><thead><tr><th colspan="3">From SED fitting to the median stacked photometry<sup>a</sup></th></tr></thead><tbody><tr><td>z<sub>phot</sub></td><td></td><td>3.1 ± 0.2</td></tr><tr><td>M<sub>*</sub></td><td>M<sub>☉</sub></td><td>(1.7 ± 0.3) × 10<sup>11</sup></td></tr><tr><td>L<sub>IR</sub><sup>b</sup></td><td>L<sub>☉</sub></td><td>(2.3 ± 0.5) × 10<sup>12</sup></td></tr><tr><td>A<sub>V</sub></td><td>mag</td><td>4.2 ± 0.3</td></tr><tr><td>T<sub>dust</sub></td><td>K</td><td>31.6 ± 6.1</td></tr></tbody><thead><tr><th colspan="3">Median value from catalog</th></tr></thead><tbody><tr><td>S<sub>3 GHz</sub></td><td>μJy</td><td>18.2 ± 0.5</td></tr></tbody><thead><tr><th colspan="3">Derived quantities<sup>a</sup></th></tr></thead><tbody><tr><td>L<sub>2–10 keV</sub><sup>c</sup></td><td>erg s<sup>-1</sup></td><td>(2.3 ± 0.7) × 10<sup>42</sup></td></tr><tr><td>L<sub>1.4 GHz</sub><sup>d</sup></td><td>erg s<sup>-1</sup> Hz<sup>-1</sup></td><td>(1.7 ± 0.1) × 10<sup>31</sup></td></tr><tr><td>q<sub>TIR</sub></td><td></td><td>2.13 ± 0.13</td></tr><tr><td>SFR<sub>IR</sub><sup>e</sup></td><td>M<sub>☉</sub> yr<sup>-1</sup></td><td>309 ± 70</td></tr><tr><td>SFR<sub>rad</sub><sup>f</sup></td><td>M<sub>☉</sub> yr<sup>-1</sup></td><td>328 ± 8</td></tr><tr><td>SFR<sub>X-ray</sub><sup>e</sup></td><td>M<sub>☉</sub> yr<sup>-1</sup></td><td>398 ± 115</td></tr></tbody></table>

**Notes.**

<sup>a</sup> Uncertainties on SED fitting-derived properties are the average of 16th and 84th percentiles of the parameter distributions. Uncertainties on derived properties come from the propagation of errors on fluxes.

<sup>b</sup> $L_{3-1100\mu m}$.

<sup>c</sup> Derived from X-ray stacking and assuming (Civano et al. 2016) a power-law photon index $\Gamma = 1.8$.

<sup>d</sup> Derived from $S_{3 GHz}$ assuming (Novak et al. 2017) a spectral index $\alpha = -0.7$.

<sup>e</sup> Derived following Kennicutt & Evans (2012), scaled to a Chabrier (2003) IMF.

<sup>f</sup> Derived following Novak et al. (2017) (see their Equation (13)).

dust emission originating from the stellar birth clouds and for the dust emission originating from the ambient (i.e., diffuse) ISM (da Cunha et al. 2008). Radio emission is also included as prescribed by da Cunha et al. (2015). The optical and infrared libraries are linked together via energy balance. The broad wavelength range used in the analysis helps in mitigating the degeneracy between stellar age, dust, and photometric redshift (da Cunha et al. 2015).

We derived a photometric redshift for the median SED of the total sample of 197 galaxies of $z_{med} = 3.1 \pm 0.2$, which is slightly higher than the typical median redshift of the general population of ALMA-identified SMGs ($z \sim 2.6 \pm 0.1$; e.g., Simpson et al. 2017), but comparable to the median redshift of K-dark SMGs in the AS2UDS survey ($z = 3.0 \pm 0.1$ Dudzevičiūtė et al. 2020).<fnref n="9" /> We summarize the other median properties of the sample, as derived from the stacked SED, in Table 1.

The uncertainties quoted in Figure 1 and Table 1 take into account the fact that, at our radio flux threshold, ~0.4% of spurious sources are expected (Smolčić et al. 2017). Under the assumption that all spurious sources inside the sampled area were picked up by our selection, we estimated that ~10 out of the 197 sources that contributed to the stacked SED, might be spurious. We performed 100 realizations of the stacking analysis by substituting 10 random sources and we added quadratically the standard deviation of the distributions of the output properties to the MAGPHYS errors from the main fit.

The median $L_{IR} (2.3 \times 10^{12} L_{☉})$ is in the so-called ultraluminous infrared galaxy (ULIRG) regime. The effective $T_{dust}^{10}$ is $31.6 \pm 6.1$ K. We also derived it by using a different approach, i.e., by fitting a modified blackbody with $\beta = 1.5$ to the median stacked FIR photometry. We obtain $T_{dust} = 33 \pm 4$ K, which is perfectly in line with the MAGPHYS estimate. Our median $T_{dust}$ is slightly colder, but broadly consistent, with that expected from the redshift evolution in main sequence galaxies from the literature ($T_{dust} = 42 \pm 3$ K at $z = 3.1$)<fnref n="11" /> (Schreiber et al. 2018a).

The inferred high dust extinction ($A_V \sim 4.2$ mag) confirms the strong obscuration of these galaxies. It is also interesting to note that the median values of SFR<sub>IR</sub> and $M_*$ lie within 0.3 dex of the star formation main sequence at $z \sim 3$ (Rodighiero et al. 2011; Speagle et al. 2014; Talia et al. 2015; Tasca et al. 2015).

The IR-based and radio-based SFR estimates are in very good agreement. The $q_{TIR}$ parameter (IR-radio correlation) is consistent with the results ($q_{TIR} = 2.28$ at $z = 3.1$) obtained for a sample of radio-selected SFGs (Novak et al. 2017), suggesting the lack of strong AGN activity in the radio band (Delhaize et al. 2017).

In order to investigate further the possible presence of hidden AGN activity in our sample, we also performed X-ray stacking. We used the publicly available CSTACK<fnref n="12" /> tool to stack Chandra soft ([0.5–2] keV) X-ray images from the Chandra-LEGACY survey (Civano et al. 2016) at the radio position of the objects in our sample. We excluded from the stack one source that has a counterpart within $l''$ in the Chandra-LEGACY point-source catalog. The stacked count rate detection, at $3.4\sigma$ significance, was converted into a stacked 2–10 keV luminosity by assuming a power-law spectrum with a slope (Luo et al. 2017) $\Gamma$ of 1.8. The X-ray based SFR is somewhat higher than, but consistent with, the IR- and radio-based estimates (Table 1), suggesting that, on average, star formation alone is enough to produce the observed $L_{2–10 keV}$. We point out that our estimate of the SFR<sub>X-ray</sub> is only a lower limit, because we did not consider the effects of intrinsic absorption. For a column density of $N_H \sim 10^{22} \text{cm}^{-2}$, typical of massive ($10^{10}-10^{11} M_{\odot}$) SFGs (Buchner et al. 2017), the SFR<sub>X-ray</sub> would be a factor of 1.3 higher, not significantly altering our conclusions.

Summarizing, the stacking analysis gave us important information at the statistical level on the nature of our sample of UV-dark galaxies: in particular it highlights that the bulk of the population is constituted by extremely dust-obscured galaxies at $z \sim 3$, which were unaccounted for, up to now, by surveys based on selections at OPT/NIR observed wavelengths.

# 4. Analysis of Individual Sources

## 4.1. The Multiwavelength Catalog

We tried to go beyond this statistical analysis and to investigate more in detail the nature of these galaxies on an individual basis. To achieve this scope, we exploited a multiwavelength catalog assembled from the deepest maps and catalogs available in the

<footnote n="9">⁹ Dudzevičiūtė et al. (2020) report that 17% of their sample of bright SMGs is K-dark. We point out that the limiting K-band magnitude in the AS2UDS field is slightly deeper than in the UltraVISTA DR4, namely $AB = 25.1 (5\sigma)$.</footnote>

<footnote n="10">¹⁰ $T_{dust}$ derived by MAGPHYS is an average luminosity-weighted value obtained by the fit of the multiple dust components.</footnote>

<footnote n="11">¹¹ We quote the light-weighted average $T_{dust}^{\text{light}}$, obtained by applying Equation (6) from Schreiber et al. (2018a). This is comparable to the temperature one would measure by using a modified blackbody model with a single temperature and an emissivity of $\beta = 1.5$.</footnote>

<footnote n="12">¹² CSTACK was created by Takamitsu Miyaji and is available at http://lambic.astroen.unam.mx/cstack/.</footnote>
COSMOS field. We searched for counterparts to each radio source of our sample in the following catalogs and/or data sets.

1. The latest public data release catalog (DR4) of the UltraVISTA survey (McCracken et al. 2012), where the Ks-band data reach a limiting $5\sigma$ magnitude $AB = 24.5 - 24.9$, fainter than the $AB = 24.0 - 24.7$ limits of the DR2 images used for the COSMOS2015 catalog (Laigle et al. 2016).

2. IRAC (channels 1–4) fluxes from the v2.0 mosaics of the Spitzer Large Area Survey with Hyper-Suprime-Cam (SPLASH<fnref n="13" />; Steinhardt et al. 2014), extracted inside an aperture with radius of $2''$.9. In particular, fluxes were extracted blindly with the code SExTRACTOR (Bertin & Arnouts 1996) and then matched with the radio positions using a search radius (Smolčić et al. 2017) of $1''$.7. We estimated the fraction of possible false associations of our radio sources with IRAC counterparts in the following way. We constructed a sample of ~170,000 mock sources in empty regions of the radio map, but close enough (within $60''$) to our 197 sources in order to sample regions with similar IRAC depth. We excluded mock sources with a counterpart in the COSMOS15 catalog within $1''$.7, in order to mimic the fifth step of our selection (Section 2) aimed at excluding radio galaxies potentially contaminated by nearby bright sources.<fnref n="14" /> Finally, we matched the mock sources to the IRAC catalogs in the same way as the real radio data and we found a 0.7% percentage of spurious associations.

3. MIPS 24 µm fluxes from the latest release (G03) of the SCOSMOS data (Sanders et al. 2007).

4. The Super-deblended catalog (Jin et al. 2018), which includes photometry from Spitzer, Herschel, SCUBA, AzTEC, and MAMBO using as priors the radio positions at 3 GHz.

5. A3COSMOS catalog (Liu et al. 2019), which includes photometry at submillimeter wavelengths from the ALMA archive. In particular, we searched for counterparts to our sources in version v.20180801 of the catalog within a radius of $0''$.8.

6. 1.4 GHz catalog from the VLA–COSMOS survey (Schinnerer et al. 2010).

We find that only ~10% of the UV-dark radio galaxies in our total sample have an observed $S_{850\mu m} \ge 4$ mJy ($>3\sigma$), which is the usual flux threshold to define bright SMGs in single-dish surveys. We point out that the flux densities at 850 µm for our sample come mainly from super-deblended photometry of SCUBA-2 sources (we find only 14 ALMA counterparts in the A3COSMOS catalog). An ALMA follow-up of all our UV-dark galaxies would be fundamental to assess the actual overlap of our radio selection with the general population of SMGs.

Thanks to our extended data set we could build the photometric NIR(MIR)-to-radio SED of 98 individual sources for which we could collect at least one FIR and one NIR-to-MIR photometric point (our primary subsample). We point out that all galaxies with an IRAC counterpart in this sample (all but one, with an estimated $z_{phot} \sim 0.5$), have radio-IRAC separations below $1''$.3. From our simulations we estimated that

within such a radius the percentage of potentially spurious associations between the radio and IRAC bands is 0.1%.

On average we had significant detections in at least ~6–7 filters from NIR to FIR per galaxy, not counting the radio band at 3 GHz, while for only ~3% of the primary subsample we could collect only two photometric points (plus upper limits), which allowed us, however, to derive an estimate of their redshift and $L_{IR}$.

We stress that the definition of dark galaxies is not absolute. Specifically, in our case it is related to the lack of counterparts for our radio-selected galaxies in the COSMOS15 catalog. By taking advantage of deeper data that became available following the COSMOS15 release, together with the availability of precise positions from the radio data, we were nonetheless able to assign a faint NIR counterpart to a fraction of our galaxies. However, since the counterpart identification of the VLA–COSMOS 3 GHz general sample was based on the COSMOS15 catalog (Novak et al. 2017; Smolčić et al. 2017), not on the deeper UltraVISTA maps, our sources were not considered in previous works for the determination of the cosmic SFRD.

We found that 24 radio sources do not have a counterpart in any NIR-to-FIR bands (we only consider detections at a $3\sigma$ significance level). For the remaining 75 sources, though being detected in at least one NIR-to-FIR band, we could not collect enough photometric points to model the observed SED. These latter two groups of galaxies constitute our secondary subsample.

## 4.2. Physical Properties and Very-high-redshift Candidates

For 98 out of 197 galaxies we could collect at least one FIR point and one NIR-to-MIR photometric point and partially reconstruct the NIR-to-radio individual SED. We call this group of 98 sources the primary subsample, and the remaining 99 sources the secondary subsample. We used again the MAGPHYS code to perform SED fitting and derive the photometric redshifts and physical properties of the individual galaxies belonging to the primary subsample. In the optical regime we assumed the upper limits from Laigle et al. (2016), while in the NIR-to-MIR bands we derived the upper limits for each undetected source directly from the maps. Uncertainties on the photometric redshifts were derived from the 16th and 84th percentiles of the probability distribution produced by MAGPHYS. This accounts for uncertainty in the photometry as well as on the model galaxy templates.

In Figure 2 we compare the redshift distribution of our primary sources with a complementary radio-selected sample with OPT/NIR counterparts (Novak et al. 2017), taken from the same VLA–COSMOS 3 GHz Large Project. The derived redshift distribution and median physical properties (see Table 2) confirm that the bulk of our population of UV-dark radio sources is indeed consisting of dusty SFGs at $z \sim 3$, but we highlight a tail of 22 newly identified very-high-z candidates ($z > 4.5$). We also compared our redshift estimates with those available in the Jin et al. (2018) catalog and found a mild agreement. We stress that the estimates of photometric redshifts in Jin et al. (2018) were computed using FIR data alone, while in our case we used the entire UV-to-radio SED, including upper limits.

A spectroscopic redshift, from $\text{Ly}\alpha$ and submillimeter emission lines, was available in the literature only for one source, COSMOSVLA3_181 (a.k.a. AzTEC-C159 Smolčić et al. 2015; Gómez-Guijarro et al. 2018), and it was consistent with our photometric estimate: $z_{spec} = 4.57$ versus our $z_{phot} = 5.11 \pm 0.4$.

<footnote n="13">¹³ http://splash.caltech.edu/</footnote>

<footnote n="14">¹⁴ The average separation between the not-isolated UV-dark radio galaxies and the COSMOS15 possible contaminant is $1''$.7.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Redshift distributions (normalized to their maximum) of our primary subsample (red), compared to the complementary sample of radio-selected galaxies with OPT/NIR counterpart (Novak et al. 2017; gray).</figcaption></figure>

<table><caption>**Table 2**
Median<sup>a</sup> of the Distribution of Physical Properties of the Individual Galaxies Belonging to the Primary Subsample</caption><thead><tr><th>Parameter</th><th>Units</th><th>Median</th></tr></thead><tbody><tr><td>z<sub>phot</sub></td><td></td><td>3.1 ± 0.1</td></tr><tr><td>M<sub>*</sub></td><td>M<sub>☉</sub></td><td>(2.0 ± 0.2) × 10<sup>11</sup></td></tr><tr><td>L<sub>IR</sub><sup>b</sup></td><td>L<sub>☉</sub></td><td>(3.2 ± 0.4) × 10<sup>12</sup></td></tr><tr><td>A<sub>V</sub></td><td>mag</td><td>4.2 ± 0.1</td></tr><tr><td>T<sub>dust</sub></td><td>K</td><td>39.5 ± 0.7</td></tr><tr><td>S<sub>3 GHz</sub></td><td>μJy</td><td>18.4 ± 0.7</td></tr><tr><td>L<sub>1.4 GHz</sub><sup>c</sup></td><td>erg s<sup>-1</sup> Hz<sup>-1</sup></td><td>(1.7 ± 0.2) × 10<sup>31</sup></td></tr></tbody></table>

**Notes.**

<sup>a</sup> For each parameter we quote the median value of the distribution. The associated errors have been estimated using the median absolute deviation (MAD), defined (Hoaglin et al. 1983) as MAD = 1.482 × median(|xᵢ - median(xᵢ)|), divided by $\sqrt{N}$, where N is the number of i objects.

<sup>b</sup> L<sub>3-1100 µm</sub>.

<sup>c</sup> Derived from S<sub>3 GHz</sub> assuming (Novak et al. 2017) a spectral index α = -0.7.

Although this latter comparison is reassuring about the robustness of our photometric redshifts, a spectroscopic follow-up of our entire population of radio-selected UV-dark galaxies would be crucial to the full understanding of their properties.

Regarding the secondary subsample, where no photometric information is available at FIR wavelengths for individual galaxies, the constraints on the physical properties are necessarily weaker. In Figure 3 we compare the mean stacked FIR SEDs of the three samples (total, primary, and secondary) cited in this work. The points were derived by mean stacking the images of our sources in five Herschel bands (from 100 to 500 µm), SCUBA 850 µm, and AzTEC 1.1mm, following the procedure by Béthermin et al. (2015). We did not apply any correction to account for possible contamination of the stacked flux by clustered neighbors (see Appendix A of Béthermin et al. 2015), because galaxies in our sample were selected to be isolated (Section 2). The almost identical position of the FIR peak in the primary and secondary stacked samples suggests that the redshift distributions are also likely similar, while the flux in the FIR bands of the secondary sample is on average

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Mean far-infrared SED of the total sample of radio-selected galaxies (red points and dashed curve), the primary subsample (gold squares), and the secondary subsample (silver diamonds).</figcaption></figure>

<table><caption>**Table 3**
SFRD<sup>a</sup></caption><thead><tr><th>z<sub>phot</sub> bin</th><th>SFRD (M<sub>☉</sub> yr<sup>-1</sup> Mpc<sup>-3</sup>)</th></tr></thead><tbody><tr><td>[0.0 – 2.0]</td><td>(6.8 ± 1.3) × 10<sup>-4</sup></td></tr><tr><td>[2.0 – 3.0]</td><td>(2.8 ± 0.5) × 10<sup>-3</sup></td></tr><tr><td>[3.0 – 4.5]</td><td>(7.1 ± 1.7) × 10<sup>-3</sup></td></tr><tr><td>[>4.5]</td><td>(5.2 ± 1.3) × 10<sup>-3</sup></td></tr></tbody></table>

**Note.**

<sup>a</sup> The tabulated values, added to their uncertainties, correspond to the upper boundary in Figure 4. The SFRDs plotted as lower boundary are a factor of ~1.2 lower.

~35% lower than the primary sample. Also, it is reasonable to deduce that their median stellar mass is about 50% lower, based on the NIR-to-IRAC fluxes (see the gray points in Figure 1).

A detailed analysis of the physical properties of our sample of UV-dark radio-selected galaxies will be presented in a forthcoming paper (M. Giulietti et al. 2021, in preparation).

# 5. Star Formation Rate Density

## 5.1. Correction for Incompleteness

To compute the density of sources and subsequently the SFRD at different cosmic times, we employed the 1/V<sub>max</sub> method (Schmidt 1968). We chose four redshift bins, the highest-redshift one being z<sub>phot</sub> > 4.5 (see Table 3). For each galaxy, we computed the maximum observable volume 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 the sum adds together comoving volume spherical shells in small redshift steps of Δz = 0.005 between z<sub>min</sub>, which is the lower boundary of the considered redshift bin, and z<sub>max</sub>, which is the minimum between the maximum redshift at which a source would still be included in the sample given the limiting flux of the survey (assumed to be constant over the entire field) and the upper
boundary of the redshift bin. The C_A and C_I constants take into account the incompleteness due to our selection criteria, and are further described below. Finally, to derive the SFRD as a function of cosmic time we added up the SFRIR of all galaxies in a given redshift bin, weighted by their individual V_max. SFRD values and their errors, quoted in Table 3, were derived via bootstrap analysis, taking into account the uncertainty on the individual redshifts. In particular, we have first replaced a random number of galaxies in the primary sample, we have assigned to each source a redshift by randomly sampling its likelihood distribution with an inverse transform sampling, and finally we have computed the SFRD. The SFRD values in each redshift bin were defined as the median of the distribution of the values over 400 realizations, while the errors were derived from the 16th and 84th percentiles of the distribution.

The C_A parameter is a correction factor that takes into account the observed area:

$$C_A = \frac{A_{\text{eff}}}{41253 \text{ deg}^2}, \qquad (2)$$

where A_{eff} = 1.38 deg^2 is the effective unflagged area covered by UltraVISTA observations (see Section 2) and 41253 deg^2 is the total sky area. The C_I parameter (correction for incompleteness) is defined 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 (3)$$

where N_primary = 98 (i.e., the primary subsample), N_secondary = 99, N_tot = 476 (i.e., the total number of radio sources with no counterpart in the COSMOS2015 catalog), and f is the ratio between the mean SFRIR of the primary and secondary subsamples. In order to derive the upper and lower boundaries of the SFRD of our UV-dark sample, we made the following considerations about our selection procedure. Given the purely geometrical nature of our fourth selection step (see Section 2), i.e., selection of isolated sources in order to ensure uncontaminated photometry, there is no reason to assume that the 197 galaxies analyzed in this work were drawn from a different distribution than the total 476 radio sources with no counterpart in the COSMOS2015 catalog. On the other hand, the 99 (out of 197) isolated galaxies for which we do not have the individual SED (i.e., the secondary subsample) might be indeed less star-forming, or at even higher redshift, than the primary subsample, and might provide a lower contribution to the SFRD. We considered three scenarios to derive the C_I parameter.

1. Case 1: the primary subsample is fairly representative of the entire radio-selected UV-dark population. Under this assumption, f = 1 in Equation (3). The SFRD values corresponding to this scenario are reported in Table 3 and are plotted (added to their uncertainties) as the upper boundary in Figure 4

2. Case 2: the galaxies in the secondary subsample do not contribute at all to the SFRD. Under this assumption, f = 0 in Equation (3). As demonstrated in the previous section, this extreme scenario is not realistic, since we did detect a signal in Herschel stacked images of the secondary subsample. Therefore we discarded this hypothesis and we did not report it at all in Figure 4.

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Cosmic star formation rate density (SFRD) history. Light-red shaded regions indicate our confidence interval, depending on the assumptions on the incompleteness correction (see Methods). The black dashed line indicates the estimate by Madau & Dickinson (2014), scaled to a Chabrier (2003) IMF. At z > 3 the SFRD estimate is mainly based on LBGs (Bouwens et al. 2012). As a comparison, we also show the samples of H-dropouts by Wang et al. (2019), radio-selected galaxies with optical counterpart by Novak et al. (2017), serendipitous ALMA Band-7 (860 or 1000 µm) continuum detections from ALPINE (Gruppioni et al. 2020), LBGs in the ASPECS field (including contribution from ULIRG-type galaxies; Bouwens et al. 2020), and HST-dark galaxies from ALPINE (Gruppioni et al. 2020).</figcaption></figure>

3. Case 3: we assume that the redshift distributions of the primary and secondary subsamples are similar, as hinted by the relative positions of the FIR peaks (see Section 4.2), and we fix f to the ratio between the stacked FIR fluxes of the two subsamples (Figure 3). In this scenario f = 0.65 in Equation (3) and the resulting SFRDs are a factor of ~1.2 lower than the upper limit set by case 1. The SFRD values corresponding to the case 3 scenario, minus their uncertainties, are plotted as the lower boundaries in Figure 4.

## 5.2. Results

After accounting for our selection function, we estimated the number density n of all UV-dark radio sources at z > 3 to be (1.3 ± 0.3) × 10⁻⁵ Mpc⁻³. If we only consider the highest-redshift bin (z > 4.5), n = (0.5 ± 0.1) × 10⁻⁵ Mpc⁻³, which is comparable to other estimates in the literature derived from smaller samples of dark galaxies at z > 5 (Gruppioni et al. 2020; Riechers et al. 2020). Existing semianalytical models (Henriques et al. 2015) and hydrodynamical simulations (Snyder et al. 2017; Pillepich et al. 2018) in the literature do not predict the early formation of such a large number of massive, dusty galaxies, and underestimate their number density by one to two orders of magnitude (see also Wang et al. 2019), with respect to our findings. Interestingly, we notice that the number density of our UV-dark radio-selected galaxies at high-redshift is comparable to that of massive quiescent galaxies at 3 < z < 4 (~2 × 10⁻⁵ Mpc⁻³) (Straatman et al. 2014; Schreiber et al. 2018b). UV-selected galaxies at z > 4 are presumed not to be abundant and star-forming enough to produce the earliest known massive quiescent galaxies (Straatman et al. 2014), suggesting that most of the star formation in the progenitors of quiescent galaxies at high-z was obscured by dust. However, it is difficult to establish a conclusive connection between high-z quiescent galaxies and
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Observed (points as in Figure 1) and best-fit (black line) SEDs of the highest-redshift galaxies from our primary subsample (z ≥ 4.5). Error bars are plotted for each data point, although in some cases they are smaller than the points. For each galaxy we quote the SFR<sub>IR</sub> and show the likelihood distribution of the photometric redshift in the inset (continued in the next page).</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 5.** (Continued.)</figcaption></figure>

the various samples of dusty starburst galaxies identified up to $z \sim 6$ (Marrone et al. 2018; e.g., the so-called 500 µm or 850 µm risers: Cox et al. 2011; Riechers et al. 2013; Ivison et al. 2016; Riechers et al. 2017; Ma et al. 2019), because of their low and uncertain number density, which is one order of magnitude lower than that of our UV-dark sources. The number density and high levels of star formation of our sample of UV-dark sources could imply an evolutionary link between this population and that of high-$z$ quiescent galaxies, where a significant fraction of the latter ones could be the descendants of the former ones at higher redshifts.

Wang et al. (2019) reported, for their sample of H-dropout-selected dusty galaxies, a space density ($\sim 2 \times 10^{-5}$ Mpc$^{-3}$ at $z>3$) similar to that of our high-redshift UV-dark galaxies, in the same redshift range. However, comparing the two populations, we find that our radio-selected galaxies produce stars at a rate that is on average three times higher. In fact, the median $L_{IR}$ of our galaxies at $z>3$ is $(6.3 \pm 0.5) \times 10^{12}$, which is about three times higher than the median infrared luminosity of $(2.2 \pm 0.4) \times 10^{12}$ reported by Wang et al. for their sample. This means that their contribution to the cosmic SFRD is also higher. In particular, it reaches $7.0 \times 10^{-3} M_{\odot} \text{yr}^{-1} \text{Mpc}^{-3}$ in the redshift range $3.0 < z < 4.5$, and $5.4 \times 10^{-3} M_{\odot} \text{yr}^{-1} \text{Mpc}^{-3}$ at $4.5 < z < 7.7$, with uncertainties of the order of 20% (Figure 4 and Table 3), therefore doubling the contribution of the

population of H-dropouts previously cited at the same cosmic epoch (Wang et al. 2019). A similar result was reported for a sample of HST-dark galaxies identified at submillimeter wavelengths in the ALPINE fields Gruppioni et al. (2020).

We compared our own selection criteria also with those by Wang et al. (2019), again to estimate the possible overlap between the two samples (see the Appendix), and we concluded that at least $\sim 82\%$ of the UV-dark radio-selected galaxies are not consistent with the H-dropout selection by Wang et al. (2019) and therefore constitute a different galaxy population. It is worth noting that we find a similar percentage (83%) when focusing only on the common redshift range between the two samples ($z > 3$).

We compared the trend toward high redshift of the cosmic SFRD of our radio-selected UV-dark galaxies with the estimate based on galaxies identified in the rest-frame UV regime (Madau & Dickinson 2014; Bouwens et al. 2020).$^{15}$ At $3.0 < z < 4.5$, the contribution of UV-dark radio galaxies to the SFRD corresponds to $\sim 10-25\%$ of the SFRD of UV-bright galaxies. This fraction rises to $\sim 25\%–40\%$ at $z > 4.5$, where we identified 22 very-high-$z$ candidates, whose SEDs are illustrated in Figure 5 along with their redshift likelihood

<footnote>$^{15}$ Bouwens et al. (2020) include the SFRD contribution from ULIRG-type galaxies in the ASPECS volume, derived from various submillimeter surveys.</footnote>
distributions. We stress that the actual contribution of radio-
selected UV-dark galaxies to the cosmic SFRD could be even
higher, because our estimates are derived without any
extrapolation to lower fluxes (hence SFRs) than our selection,
which would have been extremely uncertain. In fact, the
correction for incompleteness, which we apply only takes into
account how representative the primary sample is with respect
to our total sample of radio-selected UV-dark galaxies.

The cosmic SFRD of UV-dark radio-selected galaxies is
flatter than that of UV-bright galaxies and similar to the
complementary sample of radio-selected galaxies (Novak et al.
2017) with OPT/NIR counterparts and to that of serendip-
tously detected submillimeter galaxies in the ALPINE fields
from Gruppioni et al. (2020; see also Loiacono et al. 2021).
The estimates of the SFRD based on far-IR/submillimeter and
radio data, from different works (see also Rowan-Robinson
et al. 2016), show an almost flat trend at $z > 2$, suggesting a
significant contribution of dust-obscured activity that cannot be
recovered by the dust extinction corrected UV data. Our results
quantify which fraction of the missing amount of star formation
activity could be explained by UV-dark galaxies.

**6. Summary**

Our findings have several implications. We demonstrated
that starting from deep radio surveys, it is possible to identify a
population of massive (median $M_{\text{star}} \sim 1.7 \times 10^{11} M_{\odot}$,
extremely dust-obscured ($A_V \sim 4$) SFGs at $z \sim 3$, which are
invisible in current optical surveys and near the detection limit
of the deepest available NIR surveys. The radio selection
turned out to be particularly effective in identifying candidates
at very high redshift ($z > 4.5$), whose number density is
underpredicted by current simulations and which provide a
significant contribution to the cosmic SFRD. The comparison
of different selections of UV-/HST-dark galaxies from the
literature showed only a partial overlap between the various
samples, suggesting a possible diversity of galaxy populations
under the common dark label and the need of a multi-
wavelength approach to the search of such objects. These
results highlight the limits of our current understanding of the
processes that govern galaxy formation. To enhance our
knowledge of the dust-obscured part of the high-redshift
universe, observational efforts, especially spectroscopic follow-
ups with current and upcoming facilities like ALMA and the
James Webb Space Telescope, should be focused on UV-dark
radio-selected galaxies, since the confirmation of their red-
shifts, along with information on chemical abundances, stellar
masses and dust properties would prove essential to our
understanding of massive galaxy assembly.

This paper is dedicated to the memory of Olivier Le Fèvre.
We thank the anonymous referee for useful suggestions to
improve the paper. M.T. thanks Francesca Pozzi for useful
discussions on dust temperature. M.T., A.C., and M.G.
acknowledge the support from grant PRIN MIUR 2017
20173ML3WW_001. We acknowledge the use of Python
(v.2.7) 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 Ultra-
VISTA consortium.

**Appendix**

**Comparison with the Wang et al. Sample**

In the Wang et al. (2019) sample of H-dropouts in the COSMOS field, 4 out of 18 galaxies have a detected radio counterpart at 3 GHz. Two of them are below our chosen detection threshold (i.e., S/N > 5.5). The remaining two are compatible with our own selection, although they are not included in the sample of 197 isolated galaxies analyzed in this work. Therefore ~11% of the Wang et al. sample is overlapping with our sample. The Wang et al. selection is based on the H and IRAC2 (4.5 µm) bands. In particular, their galaxies are H-dropouts with a 5σ limiting magnitude of H < 27.4–27.8 (in the COSMOS field), with a counterpart at 4.5 µm, IRAC2 > 24.0. With respect to the Wang et al. selection criteria, our sample of 197 radio-selected galaxies can be divided into three groups.

1. *Group 1*: galaxies that are detected in the H-band (43% of the sample). These galaxies are not consistent with the Wang et al. criteria, since they are not H-dropouts.

2. Group 2: galaxies that are not detected either in the H-band or in the IRAC2 band (20% of the sample). These galaxies are excluded by the Wang et al. criteria, since they do not have an IRAC2 counterpart. We notice here that the IRAC limits are similar in the two works.

3. *Group 3*: galaxies that are not detected in the H-band, but that have a counterpart in the IRAC2 band (37% of the sample). These galaxies could be potentially overlapping with Wang et al. criteria, because the 5σ limiting magnitude in the H-band of the UltraVISTA DR4 maps is $\bar{H} = 25.2-24.1$.

In order to quantify the actual overlap between the radio-
selected galaxies in *Group 3* and the Wang et al. sample we
performed two tests. First, we stacked the galaxies in *Group 3*
in the H band. The measured median flux is $H = 26.57 \pm 0.23$,
which is brighter than the limiting magnitude of the Wang et al.
selection, meaning that at least 50% of *Group 3* galaxies are not
consistent with the Wang et al. selection. As a second test we
examined the best-fit magnitudes of *Group 3* galaxies in the
primary sample: 70% of these galaxies have $H < 27.4$ mag,
confirming the result from the first test. By summing up all
*Group 1* and 2 and 50% of *Group 3* galaxies we conclude that
at least ~82% of the UV-dark radio-selected galaxies are not
consistent with the H-dropout selection by Wang et al. (2019).
We find a similar percentage (83%) when focusing only on the
common redshift range between the two samples ($z > 3$).

**ORCID iDs**

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

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

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

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

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

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

Olivier Le Fèvre https://orcid.org/0000-0001-5891-2596

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

**References**

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

Bakx, T. J. L. C., Eales, S., & Amvrosiadis, A. 2020, MNRAS, 493, 4276

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

Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393
Béthermin, M., Daddi, E., Magdis, G., et al. 2015, A&A, 573, A113

Béthermin, M., Fudamoto, Y., Ginolfi, M., et al. 2020, A&A, 643, A2

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

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

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

Buchner, J., Schulze, S., & Bauer, F. E. 2017, MNRAS, 464, 4545

Casey, C. M., Hodge, J., Zavala, J. A., et al. 2018, ApJ, 862, 78

Casey, C. M., Scoville, N. Z., Sanders, D. B., et al. 2014, ApJ, 796, 95

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

Ciesla, L., Béthermin, M., Daddi, E., et al. 2020, A&A, 635, A27

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

Cox, P., Krips, M., Neri, R., et al. 2011, ApJ, 740, 63

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

Dannerbauer, H., Walter, F., & Morrison, G. 2008, ApJL, 673, L127

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

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

Faisst, A. L., Fudamoto, Y., Oesch, P. A., et al. 2020b, MNRAS, 498, 4192

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

Faisst, A. L., Schaerer, D., Lemaux, B. C., et al. 2020a, ApJS, 247, 61

Feigelson, E. D., & Nelson, P. I. 1985, ApJ, 293, 192

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

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

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

Gruppioni, C., Pozzi, F., Rodighiero, G., et al. 2013, MNRAS, 432, 23

Hatsukade, B., Kohno, K., Yamaguchi, Y., et al. 2018, PASJ, 70, 105

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

Hoaglin, D. C., Mosteller, F., & Tukey, J. W. 1983, Understanding Robust and Exploratory Data Analysis (New York: Wiley)

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

Ivison, R. J., Lewis, A. J. R., Weiss, A., et al. 2016, ApJ, 832, 78

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

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

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

Le Fèvre, O., Béthermin, M., Faisst, A., et al. 2020, A&A, 643, A1

Liu, D., Schinnerer, E., Groves, B., et al. 2019, ApJ, 887, 235

Loiacono, F., Decarli, R., Gruppioni, C., et al. 2021, A&A, 646, A76

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

Ma, J., Cooray, A., Nayyeri, H., et al. 2019, ApJS, 244, 30

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

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

Magnelli, B., Karim, A., Staguhn, J., et al. 2019, ApJ, 877, 45

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

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

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

Oesch, P. A., Bouwens, R. J., Illingworth, G. D., Labbé, I., & Stefanon, M. 2018, ApJ, 855, 105

Pillepich, A., Springel, V., Nelson, D., et al. 2018, MNRAS, 473, 4077

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Speagle, J. S., Steinhardt, C. L., Capak, P. L., & Silverman, J. D. 2014, ApJS, 214, 15

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

Swinbank, A. M., Simpson, J. M., Smail, I., et al. 2014, MNRAS, 438, 1267

Talia, M., Cimatti, A., Pozzetti, L., et al. 2015, A&A, 582, A80

Tasca, L. A. M., Le Fèvre, O., Hathi, N. P., et al. 2015, A&A, 581, A54

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

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

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

Wang, W.-H., Barger, A. J., & Cowie, L. L. 2009, ApJ, 690, 319