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

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

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

On the one hand, Rowan-Robinson et al. (2016), using IR data from the Herschel Telescope, found that, at z > 3, the contribution of the dusty galaxies to the SFRD could be about ten times that due to UV-bright sources. On the other hand, studying high-redshift galaxies in the infrared and submillimeter bands has always been technically tricky since observations are deeply affected by confusion problems (Lutz et al. 2011; Magnelli et al. 2013). Moreover, the faintness and the lack of counterparts in the UV/optical/NIR bands made identifying dusty galaxies at z > 3 even more challenging. The first deep surveys in the (sub)millimetric bands opened the path to studies of these galaxies at high redshift, taking advantage of the negative k-correction in this part of the spectrum (Smail et al. 1997; Hughes et al. 1998). Smail et al. (1997) observed for the first time dusty submillimeter galaxies with the Submillimeter Common-User Bolometric Array (SCUBA) instrument of the James Clerk Maxwell Telescope. These are the well-known submillimeter galaxies, i.e., galaxies with S_850μm > 3 mJy (Hayward et al. 2011), and many of them were found to be at high redshift. The South Pole Telescope has made a notable contribution to this field, allowing the identification of ∼40 candidate dusty star-forming galaxies (DSFGs) at 5 < z < 6 (Strandet et al. 2016), together with the Atacama Large Millimeter/submillimeter Array (ALMA). ALMA enabled the detection of many of these objects by tracing dust in the continuum (e.g., Scoville et al. 2014, 2016, 2017; Brisbin et al. 2017; González-López et al. 2017; Franco et al. 2018; Dudzevičiūtė et al. 2020; Smail et al. 2021).

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

Primordial dusty galaxies are now ubiquitous in theoretical models of galaxy formation as possible progenitors of massive red and dead galaxies (Cimatti et al. 2004, 2008; Cimatti 2008; Casey et al. 2014; Mancuso et al. 2016a, 2016b; Lapi et al. 2018). Moreover, it has recently been suggested that these galaxies could be the most likely hosts of high-redshift gravitational waves that the Laser Interferometer Space Antenna and the Einstein Telescope will study in the near future (Boco et al. 2020). Our work fits into this context, aiming to provide a deeper understanding of the role of obscured galaxies in galaxy formation and evolution. We propose an alternative approach to the traditional ones used to select obscured galaxies (e.g., (sub)millimeter and FIR selections), along the path of Talia et al. (2021, hereafter T21). They selected radio sources in the Cosmic Evolution Survey (COSMOS) field having no (or extremely faint) NIR counterparts (see Section 2). We call these sources radio-selected NIR-dark (RS-NIRdark) galaxies (see also Chapman et al. 2001, 2002, 2004; Enia et al. 2022). The value of ${K}_{{s}_{{\rm{lim}}}}$ adopted is ∼24.7, which is the K_s limiting flux density at 3σ of the COSMOS2015 catalog (Laigle et al. 2016).

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

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

Throughout this paper we adopt a Chabrier (2003) initial mass function (IMF) and we assume a cosmology with matter density Ω_M = 0.3, dark energy density Ω_Λ = 0.7, and Hubble parameter H₀ = 70 km s⁻¹ Mpc⁻¹. Magnitudes are given in the AB system.

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

Zoom In Zoom Out Reset image size

This work inherits the first steps of the selection process from T21. They selected 8850 sources detected with a signal-to-noise ratio S/N > 5.5 from the VLA-COSMOS 3 GHz (Smolčić et al. 2017). They choose this S/N threshold to minimize the fraction of spurious sources. Subsequently, T21 excluded all the sources within the bad areas of the VISTA Ultra-Deep (UltraVISTA) field (Laigle et al. 2016), reducing the effective area to 1.38 deg², and discarded all the multicomponent objects. They cross-correlated the resulting catalog of 5982 sources with the COSMOS2015 catalog (Laigle et al. 2016) within a search radius of 08 (Smolčić et al. 2017), finding 476 sources without a counterpart. These are the RS-NIRdark (according to the ${K}_{{s}_{{\rm{lim}}}}$ of COSMOS2015, i.e., ∼24.7 mag) galaxies, which constitute the starting point of this work. As an additional step, T21 narrowed down their selection to the 197 sources whose 3σ isophotes in the 3 GHz and χ² maps, ¹³ respectively, do not intersect. In our work we relax this constraint and include all the sources whose 5σ 3 GHz isophote does not intersect the 5σ NIR isophote—a total of 272 objects. Thus, T21 selection is stricter than ours and includes only the most isolated galaxies, while we extend the analysis to more contaminated sources, which requires the use of more refined photometric techniques. The analysis of the remaining 204 sources, which are heavily contaminated by nearby NIR-bright galaxies, will be presented in a forthcoming paper (F. Gentile et al., in preparation).

We split the 272 objects into four subsamples based on visual inspection of the NIR band maps from the second UltraVISTA Data Release (DR2, McCracken et al. 2012) and InfraRed Array Camera (IRAC) and Multiband Imaging Photometer for Spitzer (MIPS) maps from Spitzer ¹⁴ :

• 1.  
  Type 0: Isolated galaxies, with no counterparts up to MIPS 24 μm.
• 2.  
  Type 1: Isolated galaxies in both the χ² and IRAC-1 maps, having a counterpart in one or more bands.
• 3.  
  Type 2: Galaxies with a counterpart in one or more bands which are also contaminated by nearby sources. Objects belonging to this class might show a distinct brightness peak, with respect to the contaminant, in all the bands, or they might present a distinct peak in the optical and NIR bands but be unresolved in the MIR images.
• 4.  
  Type 3: Like Type 0, but close to another radio source with a counterpart in one or more bands.

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

Zoom In Zoom Out Reset image size

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

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

• 1.  
  Primary sample: sources having at least one detection above 3σ in one of the FIR/(sub)millimeter bands. This accounts for 145 sources.
• 2.  
  Secondary sample: sources with no detections above 3σ in any FIR bands. This accounts for 127 sources.

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

3.1. Median Photometry

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

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

In the FIR regime we compute the median fluxes by performing survival analysis (Isobe & Feigelson 1986), which takes into account the presence of upper limits, on the Herschel counterparts from the aforementioned Super-Deblended catalog by Jin et al. (2018). The choice of not performing the image stacking for these bands is due to the coarse resolution of the FIR observations and the ensuing difficult extraction of the sources from the maps.

3.1.1. Type 0

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

Zoom In Zoom Out Reset image size

3.2. Median SED Fitting

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

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

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

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

The resulting best-fit SED is shown in Figure 4. The median properties of the sample are consistent with those of ultraluminous infrared galaxies (ULIRGs), since the median infrared luminosity is as high as L_IR = (8.48 ± 0.05) × 10¹² L_⊙, corresponding to an SFR ∼ 100 M_⊙ yr⁻¹. This value is broadly consistent with, though slightly higher than, that by T21 (L_IR = (2.3 ± 0.5) × 10¹² L_⊙). The median photometric redshift ¹⁷ is z_med = 3.3 ± 0.2, which is also consistent with T21 (z_med = 3.1 ± 0.2) and suggests the presence of a relevant number of high-redshift (z > 4.5) galaxies in our sample.

Zoom In Zoom Out Reset image size

In order to investigate the possible presence of obscured AGNs in our sample, we extend our analysis toward the high-energy side of the spectrum. Thanks to the CSTACK online tool, we stacked X-ray data in the soft band ([0.5–2] keV) from the Chandra-LEGACY survey (Civano et al. 2016). After the exclusion of X-ray point sources, assuming N_H = 10²² cm⁻² and a photon index Γ = 1.8, the resulting rest-frame luminosity is L_{2–10 keV} ∼ 10⁴² erg⁻¹, which suggests the absence of powerful AGNs in our sample. In fact, this X-ray emission can be attributed to the star formation process alone, since it corresponds to an SFR ∼ 100 M_⊙ yr⁻¹, according to the Kennicutt & Evans (2012) relation, in perfect agreement with the SFR computed through L_IR.

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

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

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

4.1. Optical and NIR

4.2. MIR

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

4.2.1. The Deblending Procedure

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

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

We use the Python package c (Lang et al. 2016), which is particularly efficient in the deblending of sources similar to ours, ¹⁹ as already shown by Enia et al. (2022; see also Weaver et al. 2022).

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

Zoom In Zoom Out Reset image size

Table 1. Description of the Input Parameters for the Deblending Procedure with Tractor, and the Specific Values Used in Our Work

Input Parameter	Description	This Work
Map	Blended map	IRAC-1 27 × 27 pixels2 map
Positions	Positions of contaminants	The COSMOS2015 position of the
	and target	contaminant + the VLA-COSMOS 3 GHz
		position of the target
Fluxes	Best-guess fluxes for the target	SE blended flux/N;
	and the contaminants	N = number of contaminants + 1
		(= 2 if 1 contaminant)
Point-spread function	FWHM of the IRAC band considered	1.7"
Per-pixel image noise	FWHM of the map	0.1
dlnP	Precision of the linearized least squared method for the best fit	10−3
Profile	Gaussian with σ = FWHM	Gaussian with σ = 1.7

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

Download table as:  ASCIITypeset image

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

We compared the residual maps obtained after subtracting the blended source+contaminant detection using SE with those obtained after subtracting the deblended source and contaminant using Tractor: they are practically identical, meaning that the two methods recover the same total flux densities in the source+contaminant region. This constitutes additional proof of the reliability of our results. With this method, we successfully deblended 46 of the 52 Type-2 primary sources.

4.3. FIR and (sub)millimeter

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

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

4.3.1. 500 μm Risers

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

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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

4.4. Radio

4.5. X-Ray

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

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

Zoom In Zoom Out Reset image size

In the following, we focus on the sources with an estimated z_phot < 5, given the uncertainties in the dust properties at higher redshifts.

The dust luminosity L_dust has a median value over the sample of 2.8 × 10¹² L_⊙ and it varies between 2.9 × 10¹⁰ L_⊙ and 1.3 × 10¹³ L_⊙. Thus, we conclude that most of the objects in our catalog are ULIRGs. We also find six sources with L_dust > 10¹³ L_⊙, i.e., hyperluminous infrared galaxies. These very high values for the infrared luminosity might indicate that many of these sources are going through a starburst phase. In fact, using the Kennicutt relation (Kennicutt & Evans 2012) we derive a median SFR of ∼4.2 × 10² M_⊙ yr⁻¹.

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

Furthermore, one of the most interesting parameters computed by MAGPHYS is the fraction f_μ of dust luminosity due to the diffused component of the ISM, compared to the one due to the birth clouds (1 − f_μ ). The younger the galaxy is, the higher the contribution of the birth clouds, since they should host large star-forming regions with newborn stars. The parameter f_μ strongly correlates with the SED shape for dusty galaxies: larger values of f_μ correspond to lower dust temperatures and SFR (Martis et al. 2019). We can decompose the FIR luminosity into a hot (T ∼ 60–100 K) dust component and a warm one (T ∼ 20 K) and the ratio of the luminosities of these two is an efficient indicator of the dominant power in the galaxy (Kirkpatrick et al. 2012). The hot dust arises from the birth clouds or possible heating from an AGN; the warm dust is associated with the diffuse ISM and is mostly composed of larger dust grains (Kirkpatrick et al. 2015).

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

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

6.1. Computation of the SFRD

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

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

$$\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 Δz = 0.005, z_min is the lower boundary of the redshift bin, and z_max is the minimum value between the redshift upper boundary and the highest redshift at which we can detect the source, given the sensitivity of the survey.

The C_A and C_I constants take into account the incompleteness due to our selection criteria. We recall that the 272 galaxies analyzed in this work were selected from an initial sample of 476 RS-NIRdark sources only on a geometrical basis; therefore there is reason to think that the selected sources follow a distribution different than that of the whole sample. We further divided the 272 galaxies into primary and secondary samples, depending on their photometric coverage in the FIR/(sub)millimeter regime, and derived individual redshifts and physical properties only for the former. The galaxies in the secondary sample might indeed be less star-forming, or at even higher redshift, than the primary subsample, and might provide a lower contribution to the SFRD. Following T21, we define the C_I parameter (correction for incompleteness) as

$$\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{ 2 }$$

where N_primary = 145 (i.e., the primary subsample), N_secondary = 127, N_tot = 476, and f is the ratio between the mean SFR_IR of the primary and secondary subsamples.

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

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

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

6.2. Results

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

Zoom In Zoom Out Reset image size

Table 2. Values for the SFRD Found Using the Correction of Case 1 (Column 2) and Case 3 (Column 3) Described in Section 6.1

Bin	SFRD (f = 0.56)	SFRD (f = 1)
	(M⊙ Mpc−3 yr−1)	(M⊙ Mpc−3 yr−1)
0.0 < zphot < 2.0	(0.6 ± 0.1) × 10−3	(0.8 ± 0.1) × 10−3
2.0 < zphot < 3.0	(2.4 ± 0.4) × 10−3	(3.2 ± 0.5) × 10−3
3.0 < zphot < 5.0	(5.0 ± 0.7) × 10−3	(6.2 ± 0.6) × 10−3

Download table as:  ASCIITypeset image

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

In particular, T21 find an SFRD equal to (7.1 ± 1.7) × 10⁻³ M_⊙ yr⁻¹ Mpc⁻³ in the third redshift bin. Gruppioni et al. (2020) computed the SFRD up to z ∼ 6 for 56 sources detected in the ALMA Band 7 in the ALPINE fields. They find an SFRD of (1.5 ± 0.9) × 10⁻² M_⊙ yr⁻¹ Mpc⁻³ between z = 3.5 and z = 4.5 and of (0.9 ± 0.7) × 10⁻² M_⊙ yr⁻¹ Mpc⁻³ between z = 4.5 and z = 6.0. Enia et al. (2022) analyzed 17 H-dark galaxies in the GOODS-N field using the same procedure adopted in this paper and in T21 and find the SFRD at z ∼ 3 to be ∼4.5 × 10⁻³ M_⊙ yr⁻¹ Mpc⁻³, consistent with what we find in this work at the same redshift.

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

Comparing the SFRD from our sample of RS-NIRdark galaxies to the ones computed through UV-selected galaxies (Madau & Dickinson 2014; Bouwens et al. 2020), we find that the contribution to the SFRD of our galaxies at 3 < z < 5 is ∼10%–25% of the SFRD of UV-bright ones at the same redshifts, e.g., Bouwens et al. (2020), who studied the SFRD in a sample of 1362 UV-selected Lyman-break galaxies in the ALMA Spectroscopic Survey in the Hubble Ultra-Deep Field (ASPECS) between redshifts 1.5 and 10. They also include the contribution to SFRD from ULIRGs. The best-fit trend of the SFRD, computed by Madau & Dickinson (2014), is also based mainly on UV-selected sources. The SFRD trend based on UV data is steeper than that found through radio data, e.g., by Novak et al. (2017), who compute the SFRD for the full VLA-COSMOS 3 GHz Large Project (Smolčić et al. 2017) sample, complementary to the sample presented in this work. This means that UV-based studies miss an important part of the SFRD at high redshift due to dusty galaxies. Moreover, as emerges from this analysis, the SFRD related to the more dusty objects selected in the radio band, like our RS-NIRdarks, seems important at high redshift. Furthermore, we might still be missing an important fraction of obscured galaxies because of technological limitations that will be overcome thanks to new generation of telescopes, such as JWST and the Square Kilometer Array.

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

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

• 1.  
  These objects are likely highly infrared luminous (median L_dust ∼ 4 × 10¹² L_⊙), dust-obscured galaxies at a median photometric redshift of z ∼ 3.2: they are possibly young high-redshift DSFGs, as also suggested by the median value of f_μ ∼ 0.56.
• 2.  
  76 of the 272 selected sources have a photometric redshift > 3.
• 3.  
  Through the SPIRE color–color diagram, we have further proof, independent of the SED fitting procedure, that there are high-redshift candidates in our sample.
• 4.  
  The contribution of RS-NIRdark galaxies to the SFRD at 3 < z < 5 is at least ∼40% of that due to UV-selected sources, which means that the role of the dark galaxies might be crucial to understanding the still widely unknown evolution of SFRD at z > 3.

Our results remark on the importance of the obscured galaxies in galaxy formation and evolution and the need for a better understanding of their role in the broader cosmological framework. This work is a step toward a less biased view of the nature of these kinds of objects. The inclusion of more contaminated sources than in previous studies (Talia et al. 2021) allowed us to start building an improved photometric catalog and to get more information from our data on the physical properties of these galaxies.

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

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

The biggest difference in the photometric depth between our catalog and the COSMOS2020 is in the MIR bands, where IRAC-1 and IRAC-2 maps (Moneti et al. 2022) are slightly deeper than the ones we used, as summarized in Table 3. We checked the MIR fluxes of the sources that we have in common with the COSMOS2020 catalog and verified that using the deeper IRAC maps would not change our results.

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

Zoom In Zoom Out Reset image size

In this appendix, we show the results when excluding from the SED fitting of the sources with A3COSMOS counterparts the 350 and 500 μm data points, which could be less reliable than the others because of their coarse resolution, especially with respect to the ALMA data available for this particular subsample. We show an example in Figure 11 and we compare the redshifts and reduced χ² computed in the two ways in Figure 12. On average, the quality of the SED fitting stays almost the same: the median reduced χ² for the whole photometry approach is 2.69, while for the reduced photometry it amounts to 2.87. Also, the values of the redshifts stay almost the same: the median redshift in the former case is 3.5, while in the latter it is 3.3. Moreover, the median relative error is about 0.09 for both approaches.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size