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

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 ρ_*(z). The results obtained so far (Madau & Dickinson 2014; Oesch et al. 2018) show a rather consistent picture up to z ≈ 3. Several independent results indicate that the SFRD rapidly increases from z = 0 to z ≈ 1, and flattens around z ≈ 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μm > 1–2 mJy) galaxies, which are typically located around z ∼ 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_⊙ yr⁻¹; 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μm ∼ 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 ≳ 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_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_dust almost doubles from z ≈ 0 to z ≈ 4. This would imply that FIR-to-millimeter surveys could be affected by selection biases, which depend on T_dust, and hence the wavelength at which the graybody emission peaks; instead, radio emission is free from these effects. The nature of the T_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: Ω_M = 0.3, Ω_Λ = 0.7, and H₀ = 70 km s⁻¹ Mpc⁻¹.

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 (λ = 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 GHz} > 12.65 μ Jy beam⁻¹, corresponding to 5.5σ, 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 deg², 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 08 (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 χ²-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σ isophote intersects the NIR 3σ isophote of nearby objects in the χ²-image. We did not find strong evidence of gravitational lensing effects.

The final sample analyzed in this work counts 197 galaxies.

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 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.

Zoom In Zoom Out Reset image size

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 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 ± 0.2, which is slightly higher than the typical median redshift of the general population of ALMA-identified SMGs (z ∼ 2.6 ± 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 ± 0.1 Dudzevičiūtė et al. 2020). ⁹ 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 × 10¹² L_⊙) is in the so-called ultraluminous infrared galaxy (ULIRG) regime. The effective T_dust ¹⁰ is 31.6 ± 6.1 K. We also derived it by using a different approach, i.e., by fitting a modified blackbody with β = 1.5 to the median stacked FIR photometry. We obtain T_dust = 33 ±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 ± 3 K at z = 3.1) ¹¹ (Schreiber et al. 2018a).

The inferred high dust extinction (A_V ∼ 4.2 mag) confirms the strong obscuration of these galaxies. It is also interesting to note that the median values of SFR_IR and ${{ \mathcal M }}_{\star }$ lie within 0.3 dex of the star formation main sequence at z∼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 ¹² 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 ${1}^{{\prime \prime} }$ in the Chandra-LEGACY point-source catalog. The stacked count rate detection, at 3.4σ significance, was converted into a stacked 2–10 keV luminosity by assuming a power-law spectrum with a slope (Luo et al. 2017) Γ 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_X-ray is only a lower limit, because we did not consider the effects of intrinsic absorption. For a column density of N_H ∼ 10²² cm⁻², typical of massive (10¹⁰–10¹¹ M_⊙) SFGs (Buchner et al. 2017), the SFR_X-ray 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 ∼ 3, which were unaccounted for, up to now, by surveys based on selections at OPT/NIR observed wavelengths.

4.1. The Multiwavelength Catalog

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 ∼ 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.

Zoom In Zoom Out Reset image size

Table 2. Median ^a of the Distribution of Physical Properties of the Individual Galaxies Belonging to the Primary Subsample

Parameter	Units	Median
zphot		3.1 ± 0.1
M⋆	M⊙	(2.0 ± 0.2) × 1011
LIR b	L⊙	(3.2 ± 0.4) × 1012
AV	mag	4.2 ± 0.1
Tdust	K	39.5 ± 0.7
S3 GHz	μJy	18.4 ± 0.7
L1.4 GHz c	erg s−1 Hz−1	(1.7 ± 0.2) × 1031

Notes.

^a 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_i − median(x_i )∣), divided by $\sqrt{N}$, where N is the number of i objects. ^b L_{3−1100 μm.} ^c Derived from S_{3 GHz} assuming (Novak et al. 2017) a spectral index α = −0.7.

Download table as:  ASCIITypeset image

A spectroscopic redshift, from Lyα 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 ± 0.4. 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 ∼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).

Zoom In Zoom Out Reset image size

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.1. Correction for Incompleteness

To compute the density of sources and subsequently the SFRD at different cosmic times, we employed the $1/{V}_{\max }$ method (Schmidt 1968). We chose four redshift bins, the highest-redshift one being z_phot > 4.5 (see Table 3). For each galaxy, we computed the maximum observable volume as:

$$\begin{eqnarray}&&{V}_{\max }=\displaystyle \sum _{z={z}_{\min }}^{{z}_{\max }}[V(z+{\rm{\Delta }}z)-V(z)]\times \displaystyle \frac{{C}_{A}}{{C}_{I}},\end{eqnarray} \tag{ 1 }$$

where the sum adds together comoving volume spherical shells in small redshift steps of Δz = 0.005 between z_min, which is the lower boundary of the considered redshift bin, and z_max, 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 SFR_IR 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.

Table 3. SFRD ^a

zphot bin	SFRD (M⊙ yr−1 Mpc−3)
[0.0 − 2.0]	(6.8 ± 1.3) × 10−4
[2.0 − 3.0]	(2.8 ± 0.5) × 10−3
[3.0 − 4.5]	(7.1 ± 1.7) × 10−3
[ >4.5]	(5.2 ± 1.3) × 10−3

Note.

^a 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.

Download table as:  ASCIITypeset image

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

$$\begin{eqnarray}&&{C}_{A}=\displaystyle \frac{{A}_{\mathrm{eff}}}{41253\,{\deg }^{2}},\end{eqnarray} \tag{ 2 }$$

where A_eff = 1.38 deg² is the effective unflagged area covered by UltraVISTA observations (see Section 2) and 41253 deg² is the total sky area. The C_I parameter (correction for incompleteness) is defined as:

$$\begin{eqnarray}&&{C}_{I}=\displaystyle \frac{({N}_{\mathrm{primary}}+f\times {N}_{\mathrm{secondary}})}{{N}_{\mathrm{primary}}}\times \displaystyle \frac{{N}_{\mathrm{tot}}}{{N}_{\mathrm{primary}}+{N}_{\mathrm{secondary}}},\end{eqnarray} \tag{ 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 SFR_IR 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.
• 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.

Zoom In Zoom Out Reset image size

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 the various samples of dusty starburst galaxies identified up to z ∼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 (∼2 × 10⁻⁵ Mpc⁻³ 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 ± 0.5) × 10¹², which is about three times higher than the median infrared luminosity of (2.2 ± 0.4) × 10¹² 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 × 10⁻³ M_⊙ yr⁻¹ Mpc⁻³ in the redshift range 3.0 < z < 4.5, and 5.4 × 10⁻³ M_⊙ yr⁻¹ Mpc⁻³ 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 ∼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). ¹⁵ At 3.0 < z < 4.5, the contribution of UV-dark radio galaxies to the SFRD corresponds to ∼10-25% of the SFRD of UV-bright galaxies. This fraction rises to ∼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 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 serendipitously 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.

Our findings have several implications. We demonstrated that starting from deep radio surveys, it is possible to identify a population of massive (median M_star ∼ 1.7 × 10¹¹ M_⊙), extremely dust-obscured (A_V ∼ 4) SFGs at z ∼ 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 multiwavelength 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 redshifts, 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 UltraVISTA consortium.

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 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 ± 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).