Illuminating the Dark Side of Cosmic Star Formation. III. Building the Largest Homogeneous Sample of Radio-selected Dusty Star-forming Galaxies in COSMOS with PhoEBO

In the last decades, an increasing scientific interest has been growing in the elusive population of dark (i.e., lacking an optical/near-IR, hereafter NIR, counterpart) dusty star-forming galaxies (DSFGs). Although extremely promising for their likely contribution to the cosmic star formation rate density (SFRD) and for their possible role in the evolution of the first massive and passive galaxies around z ∼ 3, the difficulty in selecting statistically significant samples of dark DSFGs is limiting their scientific potentialities. This work presents the first panchromatic study of a sample of 263 radio-selected NIR-dark (RS-NIRdark) galaxies discovered in the COSMOS field following the procedure by Talia et al. These sources are selected as radio-bright galaxies (S 3 GHz > 12.65 μJy) with no counterpart in the NIR-selected COSMOS2020 catalog (Ks ≳ 25.5 mag). For these sources, we build a new photometric catalog including accurate photometry from the optical to the radio obtained with a new deblending pipeline (Photometry Extractor for Blended Objects, or PhoEBO). We employ this catalog to estimate the photo-zs and the physical properties of the galaxies through an spectral energy distribution-fitting procedure performed with two different codes (Magphys and Cigale). Finally, we estimate the active galactic nucleus contamination in our sample by performing a series of complementary tests. The high values of the median extinction (A v ∼ 4) and star formation rate (SFR ∼ 500 M ⊙ yr−1) confirm the likely DSFG nature of the RS-NIRdark galaxies. The median photo-z (z ∼ 3) and the presence of a significant tail of high-z candidates (z > 4.5) suggest that these sources are important contributors to the cosmic SFRD and the evolutionary path of galaxies at high redshifts.


INTRODUCTION
Our picture of the Universe is essentially the offspring of the instruments we use to observe it.Until a few years ago, the Hubble Space Telescope (HST) was the most powerful facility available to study the high-z Universe.The data from this facility have put several critical constraints on important quantities as a function of cosmic time (see, e.g., the review by Madau & Dickinson 2014 and references therein).One of the main results concerned the reconstruction of the cosmic Star Formation Rate Density (SFRD; i.e. the amount of stellar mass formed in the Universe per each year and Mpc 3 ) back to the epoch of reionization (z ∼ 8).The vast majority of these studies agreed in describing a SFRD constantly growing ∝ (1 + z) 1.7 from z ∼ 0 up to z ∼ 2.5 (the socalled cosmic noon) and -then -decreasing ∝ (1+z) −2.9 up to z ∼ 8 (Madau & Dickinson 2014).Given the spectral coverage of HST, however, our census of galaxies in the first few billion years of cosmic history was forcedly limited to optically-bright galaxies.Given the redshift ranges explored, these sources were mainly rest-frame UV /optically-bright galaxies.
Even in that time, however, an increasing number of evidences suggested that our picture of galaxy evolution based on these UV-bright sources was quite far from complete.Above all, the discovery of a significant population of massive (M ⋆ > 10 11 M ⊙ ) and passive (sSFR< 10 −11 yr −1 ) galaxies already in place at z ∼ 3 (i.e. when the Universe was just 2 Gyr old; see some examples in Straatman et al. 2014 andSchreiber et al. 2018b) was -probably -the most outstanding.The number density of these sources (n ∼ 2 × 10 −5 Mpc −3 ) resulted in being more than two orders of magnitude higher than that computed for UV-bright galaxies at z > 4.5 (see e.g. the discussion in Straatman et al. 2014;Schreiber et al. 2018b;Valentino et al. 2020).This suggested that, when not accounting for UV-dark star forming galaxies, we are missing a significant fraction of the progenitors of high-z passive galaxies.
Observations carried out with facilities able to explore other regions of the electromagnetic spectrum such as the sub-millimeter array camera SCUBA on the James Clerk Maxwell Telescope (e.g.Smail et al. 1997;Hughes et al. 1998) or the Herschel Space Observatory (e.g.Gruppioni et al. 2013;Burgarella et al. 2013) confirmed that selections based only on the optical regime actually miss a significant fraction of high-z galaxies (see e.g.Blain & Phillips 2002;Weiß et al. 2013).Among the main missing population of galaxies, the so-called "Dusty Star-Forming Galaxies" (DSFGs; see e.g.Lagache et al. 2005;Casey et al. 2014 for reviews) are worth a deeper discussion.These sources are characterized by high values of the Star Formation Rate and high stellar extinction due to the presence of a significant amount of dust.Hence, these sources are extremely faint or even undetected in the optical/NIR regime observed by HST.Despite their elusiveness, several studies targeting these sources (e.g.Wang et al. 2019;Gruppioni et al. 2020;Talia et al. 2021;Zavala et al. 2021;Enia et al. 2022;van der Vlugt et al. 2022) are unanimous in finding that the contribution of these galaxies to the cosmic SFRD in the high-z regime is definitely significant (up to 40% the contribution of optically bright galaxies at z > 4.5; see e.g.Talia et al. 2021;Enia et al. 2022;Behiri et al. 2023).The inclusion of these new sources in the cosmic census of the high-z galaxies could even change the behavior of the SFRD at z > 3 (favoring a flatter behavior before the cosmic noon; see e.g.Gruppioni et al. 2020) and solve the puzzle of the missing progenitors of massive galaxies (e.g.Toft et al. 2014).
It is crucial to underline how all the aforementioned studies on DSFGs are generally based on limited samples of galaxies.The main reason for this is the extreme faintness of these sources in the optical/NIR regime, making their identification feasible only at longer wavelengths.For this reason, these galaxies are generally selected in the IR or (sub)mm regimes (e.g.Smail et al. 2021;Dudzevičiūtė et al. 2021;Manning et al. 2022;Giulietti et al. 2022), taking advantage of the bright thermal emission of the dust and the intensely negative k -correction that makes these sources bright even at high redshifts.
IR/(sub)mm selections, however, are affected by -at least -two significant issues.The first one is the low sensitivity of single-dish (sub)mm telescopes.This, combined with their large beam sizes, makes it hard to identify the multi-wavelength counterparts of each galaxy.In principle, this problem could be overcome by employing state-of-the-art facilities such as the Atacama Large (sub)Millimeter Array (ALMA) with higher sensitivities and smaller beams.However, the deep surveys conducted with these telescopes (e.g.Franco et al. 2018;Casey et al. 2021) are currently mapping small volumes, hence are affected by low statistics and are prone to cosmic variance.The second issue concerning IR/(sub)mm selections resides in the still-debated properties of the dust in the high-z regime.For instance, a correlation between the dust temperature and the redshift (e.g. the one suggested by Béthermin et al. 2015, Faisst et al. 2017and Schreiber et al. 2018a) could produce a significant selection bias.
A solution to both the aforementioned issues can reside in the radio selection of DSFGs (e.g.Talia et al. 2021;Enia et al. 2022;Behiri et al. 2023).In "normal" (i.e.non-active) galaxies, radio photons can be generated by free-free emission in HII regions and synchrotron emission from relativistic electrons accelerated in supernovae remnants.In addition, they are not affected by the presence of dust.For these reasons, once excluded the nuclear origin of these photons, they could represent an excellent unbiased tracer of star formation.Therefore, as noted by Talia et al. (2021), the selection of radio-bright galaxies with a faint or undetected counterpart in the optical/NIR regime could provide a significant sample of likely DSFGs.
This paper expands the studies by Talia et al. (2021) and Behiri et al. (2023), which aimed to assemble the widest homogeneously-selected sample of candidate DS-FGs available in the current literature.We take advantage of the multi-wavelength coverage in the Cosmic Evolution Survey (COSMOS; Scoville et al. 2007) field and the high sensitivity of the VLA-3GHz COSMOS Large Project (Smolčić et al. 2017) to collect a sample of 323 Radio-Selected NIRdark Galaxies (RS-NIRdark galaxies in the following).In this paper, we present the photometric catalog of these sources with new accurate photometry extracted from the optical to the radio regimes.We also employ a SED-fitting procedure to estimate photo-z s and physical properties of our galaxies.
The paper is structured as follows.In Section 2, we introduce the sample selection.In Section 3, we describe the maps employed to extract the photometry in the optical-to-MIR regime and the additional photometry retrieved at longer wavelengths.In Section 4, we present the new PhoEBO pipeline (Photometry Extractor for Blended Objects) used to extract the photometry and its validation on simulated data.Section 5 is focused on the SED-fitting procedure employed to assess the photo-z s and the physical properties of the galaxies in our sample.Furthermore, in Section 6 we analyze in detail the possible AGN contamination in our sample.In Section 7, we discuss what the results presented in this paper tell us about the nature of the RS-NIRdark galaxies and compare them with previous analogous studies in the current literature.Finally, we draw our conclusions in Section 8. Throughout this paper, we use AB magnitudes (Oke & Gunn 1983), employ a Chabrier (2003) Initial Mass Function (IMF) and assume a concordance cosmology ΛCDM with H 0 = 70 km s −1 Mpc −1 and (Ω tot , Ω Λ , Ω m )=(1,0.7,0.3).

SAMPLE SELECTION
To assemble a large sample of candidate DSFGs with complete photometry from the optical to the radio, we focus on the galaxies in the COSMOS field.We collect our sources by performing a selection analogous to that employed by Talia et al. (2021) and summarized here: 1. We start from the VLA-COSMOS 3GHz Large Project catalog (Smolčić et al. 2017).We select 8850 radio sources with an S/N> 5.5 (i.e brighter than 12.65 µJy beam −1 ) over the full 2 deg 2 coverage of the survey.This cut allows us to assemble a sample of radio-bright galaxies, limiting the likely contamination of fake sources to 0.4% (Smolčić et al. 2017).
2. To limit the expected presence of galaxies hosting AGN (common in radio-selected catalogs; see e.g.Bonzini et al. 2013 andNovak et al. 2018), we remove from the sample all the sources flagged as "multi-component" in the initial catalog.We underline that this flag is the result of a visual inspection of the sources in the full catalog by Smolčić et al. (2017).Therefore -at this stage -we cannot exclude that some multi-component radio sources are still present in our sample (see the discussion in Section 6).
3. To take advantage of the multi-wavelength coverage of the COSMOS field, we exclude from the sample all the radio sources lying outside the Ul-traVISTA survey footprint (Laigle et al. 2016).This further limits the sample to 5982 galaxies and the effective area mapped by this study to 1.38 deg 2 .
4. Finally, we cross-match the resulting sample with the two versions of the COSMOS2020 catalog (the "Classic" and "Farmer"; see Weaver et al. 2022), removing all the sources with a match within 0.7".The final sample is composed of 323 galaxies.
The last step aims to include in the catalog only sources without a significant counterpart in the optical/NIR bands.However, it is crucial to underline how the source detection in the COSMOS2020 catalog is performed on the χ 2 -image produced through the SWarp software (Bertin et al. 2002) by combining the maps in the i and z band from the Subaru telescope and the Y, J, H and Ks bands from the VISTA telescope (see Weaver et al. 2022).The flux in the generic pixel I of this image is computed as a weighted average of the fluxes f i in the different photometric bands, with the weights w i provided by the uncertainty maps: Given this definition, we cannot exclude that some of the sources in our catalog could have a significant counterpart in a few individual NIR bands.Moreover, through a visual inspection of the χ 2 -map employed by Weaver et al. (2022) for detecting the sources in the COSMOS2020, we notice that 60 galaxies of our sample have a detection in that image with a S/N higher than 1.5σ (i.e. the threshold used in Weaver et al. 2022).Some of these sources (∼ 10%) are inside regions masked in the COSMOS2020 that were employed in the COS-MOS2015, while most of the lasting sources are close to a bright companion, therefore it is likely that they were not properly deblended in the COSMOS2020.Since these galaxies are not included in the COSMOS2020, they should be part of our sample of RS-NIRdark galaxies.However, due to the different photometry, we expect their properties to differ from the rest of the sample.Therefore, we do not consider these sources in all the statistical analyses performed in this paper aimed to characterize our population of galaxies: for all these studies, we will focus on the lasting 263 galaxies.Finally, we underline that the S/N cut applied to the radio catalog excludes from the sample all the faintest radiosources.Although this step could potentially cause the lack of significantly high-z galaxies (see e.g.Casey et al. 2019), it is necessary to exclude a high contamination by fake sources (this rate would increase up to 25% just including the sources with 5 < S/N < 5.5, see Smolčić et al. 2017).
A significant difference between this study and the previous papers of this series (Talia et al. 2021 andBehiri et al. 2023) resides in the last step of the selection: while the previous works selected a sample of 476 RS-NIRdark galaxies without a counterpart in the COS-MOS2015 catalog (Laigle et al. 2016), we employ the updated (and deeper) COSMOS2020 catalog (Weaver et al. 2022).This step allows us to exclude from the sample 153 galaxies un-detected in the COSMOS2015 but revealed by COSMOS2020.
More in detail, 16 sources were only included in the "Farmer" catalog, 21 only in the "Classic", and 116 in both the catalogs.Finally, improving the works by Talia et al. (2021) and Behiri et al. (2023), we do not employ any additional selection aimed to avoid sources with a bright contaminant in the vicinity, but we analyze the entire sample, performing a more accurate photometry extraction in presence of contamination or source blending.A more comprehensive comparison with Talia et al. (2021) and Behiri et al. (2023) can be found in Appendix A.

Analyzed maps
The analysis of DSFGs requires the most comprehensive wavelength coverage to account for all the physical processes taking place in these complex objects (i.e.stellar emission, dust obscuration, thermal emission, and non-thermal processes).For this purpose, we extract accurate photometry in the optical-to-MIR regimes by analyzing the following maps: Further details on the maps employed in the photometry extraction can be found in Table 1.

Additional photometry
To analyze the dust emission in the FIR/(sub)mm and the non-thermal processes emitting at radio frequencies, we retrieve additional photometry for our sources by cross-matching our sample with other catalogs in the current literature: • FIR: We cross-match our catalog with the 2020 version of the SuperDeblended catalog by Jin et al. (2018).Since the procedure followed in building that catalog employs as positional prior the sources in the VLA-COSMOS 3GHz Large Project, it is possible to retrieve FIR photometry obtained with Spitzer (24 µm; Le Floc'h et al.  et al. 2007) for all the galaxies in our sample.Further details on the employed maps and on their depth can be found in Jin et al. (2018).
• (sub)mm: We cross-match our sample with the catalog of the Automated Mining of the ALMA Archive in the COSMOS Field (A3COSMOS) survey (v.20200310;Liu et al. 2019).This catalog contains (sub)mm continuum fluxes for all the sources in the COSMOS field observed with ALMA.Since the coverage of the A3COSMOS survey is not uniform, only a tiny percentage of our RS-NIRdark sources has a counterpart in this catalog (32 galaxies within a matching radius of 1", ∼ 10%; Liu et al. 2019), with central bandwidth and depth strongly variable for different sources.
• Radio: Finally, we retrieve additional radio photometry by cross-matching our catalog with the public catalog by Schinnerer et al. (2010)  Extracting accurate photometry for the sources in our catalog is a challenging task.The main limitation is represented by the possible presence of bright contaminants close to our (extremely faint or even undetected) galaxies in the optical/NIR bands.This issue becomes highly significant in the IRAC channels, where the large (up to FWHM= 2 ′′ in ch3 and ch4) and irregular Point Spread Function (PSF) makes it almost impossible to blindly deblend multiple sources without priors on their positions and shapes.
Basic deblending algorithms such as that implemented in sExtractor (Bertin & Arnouts 1996) need a mini- mum contrast between the different blended components inside the same cluster of pixel to recognize the presence of multiple sources.This level of contrast is not reached in the IRAC bands, making these algorithms unfit to analyze our galaxies.
More complex softwares such as Tractor (Lang et al. 2016) and Farmer (Weaver et al. 2022), relying on profile-fitting, can overcome this problem by extracting prior information on the position and shape of the different objects through a highresolution image in which all the blended components are present and distinguishable at the same time (e.g. the χ 2 -map employed by Weaver et al. 2022).As before, these techniques cannot be applied to our RS-NIRdark galaxies, since the only bands in which our galaxies are robustly detected -together with the contaminants -are the IRAC channels, where it is generally impossible to distinguish the different sources.Therefore, in order to build the photometric catalog, we developed a new deblending pipeline called PhoEBO (Photometry Extractor for Blended Objects; Gentile et al. 2023)1 relying on a slightly modified implementation of the method employed in several previous studies (e.g.Labbé et al. 2006;Endsley et al. 2021;Whitler et al. 2022).The whole procedure is summarized in Figure 1 and follows these steps: 1. We start from two high-resolution "detection images" and a low-resolution "target image" to be deblended.One detection image must contain information on the RS-NIRdark galaxy we want to analyze, the other on the contaminant sources.In our case, we choose the 3GHz and the Ks maps, respectively.We underline that we minimize the risk of biases by employing two maps with a comparable PSF FWHM (0.75" and 0.78", respectively).
2. We employ the basic deblending algorithm included in the Sep library (Bertin & Arnouts 1996;Barbary 2016) to identify the different sources in the two detection images and to produce the segmentation maps (i.e. a set of images assigning each pixel to one of the galaxies in the field).To include all the galaxies in the maps, we set a low detection threshold (2σ) and a minimum of 5 contiguous pixels in each detected source.
3. We combine the two segmentation maps and use their product as a binary mask to isolate the sources present in the detection images and produce a different image for each of them.
4. We convolve each image produced in the previous step with a matching kernel to homogenize the PSF with that of the target image.The matching kernel is produced with the "photutils" library (Bradley et al. 2020) using the ratio of Fourier transforms (Gordon et al. 2008;Aniano et al. 2011).The PSF of the detection images is modeled as a 2D Gaussian with a FWHM equal to the PSF FWHM of the Ks band reported in McCracken et al. (2012) and Weaver et al. (2022).
The PSF of the optical/NIR bands are obtained in the same way (i.e. with Gaussians with fixed widths, see Table 1), while the PSFs of the four IRAC channels are downloaded from the IRSA archive 2 5. We normalize all the resulting images and co-add them into a "model image".This one is resampled to match the pixel size of the target image (if needed) and then fitted to the target image by multiplying each normalized source by a freeparameter α i .The fit is performed with the Scipy library (Virtanen et al. 2020), aiming to minimize the χ 2 between the model and the target image .
6.We multiply all the components of the model image for the relative α i obtained through the fitting procedure.Then, we subtract all the resulting images of the contaminants from the original target image.In doing so, we get a residual image containing only the source present in our sample.
7. Finally, we perform aperture photometry with Photutils on this residual image by employing a fixed diameter of 2 arcsec.The local background is computed in an annulus between 1.5 and 2 times the radius of the aperture and subtracted from the extracted counts.The counts are then converted to AB magnitudes and to micro-Jansky through the photometric zeropoints employed in Weaver et al. (2022), corrected for the systematic offset reported in Weaver et al. (2022).This last correction is needed to account for the systematic mismatch between the spec-z s and the photo-z s computed in the COSMOS field (Laigle et al. 2016;Weaver et al. 2022).
The whole procedure described above allows us to estimate the fluxes in all the NIR and IRAC bands reported 2 https://irsa.ipac.caltech.edu/data/SPITZER/docs/irac/ in Section 3.1.As prescribed by the IRAC Instrument Handbook, we correct the fluxes measured through aperture photometry by the factors reported in the IRAC Instrument Handbook3 .For the optical bands -with a smaller PSF FWHM than the radio map at 3GHz and the Ks map -we subtract the contaminants with PhoEBO by using the HSC-i band as a detection image.Finally, we estimate the photometric uncertainties as: where the σ are obtained from the weight maps and the sum is extended to all the pixels in the aperture employed in the estimation of the flux.As prescribed by Weaver et al. (2022), we correct the photometric uncertainties in the optical/NIR bands with the multiplicative factors reported in Table 1.This step is required to account for the expected under-estimation of the uncertainties through Equation 2 due to the presence of correlated noise in the analyzed maps (e.g.Leauthaud et al. 2007).Since the weight maps are not affected by the subtraction of the contaminants, we underline that the procedure followed here to estimate the photometric uncertainties is totally consistent with that employed by Weaver et al. (2022) for the optical/NIR bands.Regarding the IRAC bands, we can compare our uncertainties with those by Weaver et al. (2022) by running PhoEBO on the 153 RS-NIRdark galaxies selected by Talia et al. (2021) and excluded from this study being revealed in the COSMOS2020 (see Section 2).We obtain that the median ratio between our uncertainties and those included in the "Classic" COSMOS2020 is in the order of ∼ 3.This result can be explained by the likely underestimation of the IRAC uncertainties found by Weaver et al. (2022) in the COSMOS2020 catalog.

Validation of PhoEBO
Thanks to the small difference in the PSFs, the photometry extraction performed by PhoEBO in the optical/NIR bands does not differ significantly from that performed by "classic" algorithms such as sExtractor (Bertin & Arnouts 1996).Nevertheless, the extraction in the IRAC channels is quite different.Therefore, in order to obtain reliable results on the physical properties extracted from the photometric catalog, we need to validate the performances of the pipeline in the MIR regime.We validate the results of the PhoEBO pipeline performing extensive simulations of blended galaxies in the four IRAC channels.The simulation procedure recalls the philosophy discussed in Section 4.1: 1. We simulate two noise maps through a random Gaussian generator, requiring an rms compatible with the 3GHz images observed in the VLA-COSMOS survey at 3 GHz and with the Ks band images observed in the UltraVISTA survey (conservatively, we employ the sensitivity reached in the deep stripes).
2. We simulate a radio source and a NIR-bright contaminant.Both the galaxies are simulated as 2D Gaussians on the noise maps generated in the previous step.Since the vast majority of the RS-NIRdark galaxies are unresolved at 3 GHz, we choose a FWHM=0.7"for the radio sources (i.e. the width of the synthesized beam in the VLA-COSMOS survey).On the contrary, to account for the presence of partially resolved contaminants, the FWHMs of the Gaussians in the Ks band are uniformly sampled in the range [0.7,1.4]arcsec (i.e. between one and two PSF FWHM).The normalizations of the two Gaussians are chosen to obtain a S/N uniformly sampled in the range [5.5,8] and [3,8], for the radio and Ks images, respectively.The radio sources are placed in the center of each image, while the contaminants are allowed to space in the range [0.7,1.4]arcsecfrom the center.This range allows us to study the accuracy of the algorithm as a function of the blending between the two sources.The lower limit on the distance is chosen as 0.7", since we recall that -due to the selection described in Section 2 -all the NIR sources with a separation lower than this threshold are considered NIR counterparts of the radio signal.
3. We convolve each of these Gaussians with a set of matching kernels computed as prescribed in Section 4.1 to obtain the sources as observed in each of the four IRAC channels.5. Finally, we co-add the two images and add Gaussian noise with the same rms expected in the four IRAC channels.
To explore the whole parameter space of the randomly sampled quantities, we simulate ∼ 10 3 images.Then, we run the PhoEBO pipeline on the IRAC-like images, employing as detection images the high-resolution data in the radio and Ks bands.In order to assess the performances of the pipeline, we compare the fluxes reported by PhoEBO with those obtained by performing a standard aperture photometry on the isolated IRAC-like images (i.e.those obtained in point 3 of the simulation procedure, before adding the contaminant).
Computing the ∆mag on the whole dataset, we obtain a median(∆mag)∼0.03 in all the IRAC channels, without any significant difference between the different bands.Similarly, we obtain a std(∆mag)∼0.15 in all the channels.Moreover, the employment of simulations allows us to estimate the ∆mag as a function of the different parameters employed during the simulation procedure (Figure 2).We find interesting -albeit expected -trends with the IRAC flux and with the angular separation between the high-resolution images, with lower accuracy achieved on more blended sources and IRACfainter RS-NIRdark galaxies.

PHYSICAL PROPERTIES FROM SED-FITTING
To determine the nature of the RS-NIRdark galaxies, we need to assess their photo-z s and physical properties.We do this through an SED-fitting procedure.To test the robustness of our results against the characteristics of different codes, we base our analysis on two algorithms: Magphys+photo-z (da Cunha et al. 2008;Battisti et al. 2019) and Cigale (Boquien et al. 2019).

SED-fitting with Magphys
Magphys (da Cunha et al. 2008) is a physicallymotivated SED-fitting code based on the energy balance between stellar attenuation and thermal dust emission.The software estimates the physical properties of a galaxy by comparing its optical-to-radio broad-band photometry with more than a million templates, including stellar emission, dust attenuation, thermal dust emission and non-thermal radio emission.The stellar emission is considered by combining the Single Stellar Populations (SSPs) by Bruzual & Charlot (2003) with an exponentially declining Star Formation History (SFH) with random bursts of star formation superimposed on the continuum.The dust attenuation is included as prescribed by Charlot & Fall (2000), with the addition of a 2175A feature accounting for the young stars born in denser clouds.A three-component model accounts for thermal dust emission.These components include the "hot dust" (i.e.warmed up by young stars in the birth clouds), the "cold dust" (i.e., present in the diffuse interstellar medium), and the characteristic emission by the Polycyclic Aromatic Hydrocarbures (PAHs).The first two components are modeled with a modified gray body, while the PAH emission is modeled with an

SED-fitting with Cigale
To avoid possible biases arising from the use of a single SED-fitting code, we involve the software Cigale (Boquien et al. 2019) in the analysis.This code is based on a similar energy-balance principle as Magphys but allows a larger customization of the libraries employed in building the templates.In this paper, we start by using a setup as close as possible to Magphys to achieve consistent results.A detailed investigation of the different modules included in Cigale will be discussed in a forthcoming paper.The SSPs, SFH, dust attenuation and radio emission are the same as discussed in Section 5.1.A significant difference between the two codes resides in the treatment of the thermal dust emission.Cigale does not have a single model including both the gray-body thermal emission and the PAH emission as Magphys.

SED-fitting results and comparison between the codes
The two codes are applied to the photometric catalog presented in Section 4. Before running the codes, we correct the photometry in the catalog for the galactic extinction.We employ the dust maps by Lenz et al. (2017) and the extinction law by Fitzpatrick & Massa (2007).This correction is performed through the python libraries DustMaps (Green 2018) and Extinction (Barbary 2016).To account for the possible biases in the photometry extraction (see Section 4.2), to account for the known under-estimation of the photometric uncertainties by the Photutils library em-ployed in PhoEBO (see e.g Leauthaud et al. 2007;Laigle et al. 2016;Weaver et al. 2022) and to allow the SED-fitting codes to explore a wider region in the photometry space, we add in quadrature 0.15 mag to the photometric uncertainties included in the catalog for the optical/NIR/MIR bands, following Laigle et al. (2016) and Weaver et al. (2022).The output of the SED-fitting codes Magphys and Cigale are summarized in Figure 4, while the median values of the photoz s and the main physical properties for both codes are reported in Table 2.In generating the histograms and in computing the medians, we exclude from the sample 57 sources that could host an AGN (see the discussion in Section 6), in order to avoid possible biases coming from the SED-fitting performed with templates not including nuclear activity.Similarly, we do not include in our analysis 10 sources with no other robust detections than in the radio (the so-called "type 0" of Behiri et al. 2023, see that study for the possible models of these sources) for which the properties estimated through SED-fitting would be highly unreliable.The good convergence of the SED-fitting procedure is ensured by the median reduced χ 2 < 2.5 obtained by both codes on the analyzed sample.Since for our galaxies the constraints on stellar population in the optical/NIR regimes are limited -in the majority of the cases -to upper limits, we compute the SFR starting from the infrared luminosity through The uncertainties on the median properties are estimated as the Median Absolute Deviation: M AD = 1.482 × median(|xi − median(xi)|) (Hoaglin et al. 1983) divided by √ N , where N is the number of galaxies in the sample.
For each quantity, we also report the dispersion computed as half the symmetrized interval between the 16th and the 84th percentiles.
the relation by Kennicutt & Evans (2012) rescaled to a Chabrier (2003) IMF.This quantity is more robustly estimated, since ∼ 60% of the galaxies in the sample have a detection at FIR/(sub)mm wavelengths.However, we must underline how this quantity is more unconstrained for the other ∼40% of the sources in the sample without these detections.Nevertheless, since the two codes employed in the SED-fitting procedure rely on the energy balance and include the radio fluxes, some constraints on the infrared luminosity can also come by the modeled dust attenuation and by the radio luminosity (da Cunha et al. 2015;Battisti et al. 2019;Boquien et al. 2019).Finally, we notice that no significant discrepancy is visible when comparing the SFR estimated through the L IR and the L 1.4GHz (see e.g.Kennicutt & Evans 2012;Novak et al. 2017) for the galaxies with FIR/(sub)mm detections and those without.This result can be partly explained by the fact that Magphys assumes a constant radio-infrared correlation with a q TIR centered on 2.34 and with a 1σ dispersion of 0.25 (da Cunha et al. 2015;Battisti et al. 2019).This value is broadly compatible with the median q T IR = 2.1 ± 0.3 computed for the galaxies in our sample (see Section 6.3).
As it can be seen, the results of the two software for what concerns the photo-z s, and most of the physical properties are broadly compatible.The slight difference between the outputs for these quantities can be explained by some minor differences between the codes (e.g. the treatment of the upper limits, representing the majority of the constraints in the optical/NIR regime for our galaxies; see e.g.Battisti et al. 2019 andBoquien et al. 2019).Major differences hold for the dust mass and temperature.As discussed in Sections 5.1 and 5.2, Magphys and Cigale report the luminosity-(T L Dust ) and mass-weighted (T M Dust ) dust temperature, respectively.Since these quantities weight differently the hot and cold components of the dust (see e.g.Liang et al. 2019;Sommovigo et al. 2020), it is not surprising the difference in the two distributions reported in Figure 4.Moreover, this difference in the two codes can also explain the discrepancy in the two estimates of the dust mass, being this quantity computed starting from the dust temperature.In the following, for consistency with the previous papers of the series (Talia et al. 2021;Behiri et al. 2023), we will consider the results obtained by Magphys.A final note concerns the accuracy of our photo-z s.Unfortunately, due to the elusive nature of the sources in our sample, it is not easy to retrieve spectroscopic redshifts for our RS-NIRdark galaxies from the current literature.One spectroscopic redshift can be found in Jin et al. (2022), where they targeted one of our galaxies with a spectral scan performed with the ALMA and NOEMA interferometers.For our galaxy RSN-436 (ID 3117 in Jin et al. 2022), they obtained z spec = 3.545, which is in good agreement with our estimates (z M = 3.3 ± 0.3 and z C = 3.4 ± 0.2.To increase the spectroscopic coverage of our sample, we already have planned approved spectroscopic follow-ups with (sub)mm facilities such as ALMA and NOEMA (Gentile et al., subm.).

AGN CONTRIBUTION
While in Section 1 we discussed the possible biases affecting the FIR and (sub)mm selection of DSFGs, this section focuses on the main bias possibly affecting our radio selection: the presence of strong Active Galactic Nuclei (AGN).Since both star formation and radio activity can generate radio emission, we cannot exclude that some of our galaxies are hosting an AGN.For instance, Bonzini et al. (2013); Novak et al. (2017) and Enia et al. (2022), analyzing samples of radio-bright galaxies selected at 1.4 GHz, reported significant fractions of AGN in their catalogs, spanning from 10% (for the faintest radio fluxes) to 100% (for the radio-brightest sources).Fortunately, our selection procedure focuses on galaxies without a significant optical/NIR counterpart.This step allows us to remove from the sample all the brightest AGN.This selection, however, does not allow us to remove from the sample the obscured AGN, where a significant amount of dust absorbs the optical/NIR emission (see e.g. the review by Hickox & Alexander 2018 and references therein).The presence of these sources in our sample can be unveiled by searching for AGN tracers at longer wavelengths (namely the IR and radio regimes), or in the X-ray, taking advantage of the broad wavelength coverage of our sample.Estimating the fraction of obscured AGN and properly accounting for their presence is crucial for two main reasons.Firstly, the presence of an extra IR component due to a dusty torus surrounding the AGN and/or a radio excess due to nuclear activity can bias the SFR when obtained from the IR/radio luminosity through the relation by Kennicutt & Evans (2012).Similarly, the employment of galaxy-only templates could bias all the other physical properties (and the photo-z s) estimated through SEDfitting.Secondly, the contamination affecting the selection procedure has to be considered in the determination of the statistical properties of the sample (e.g. the luminosity function and the contribution to the SFRD).

Visual inspection
The first selection of galaxies hosting AGN in our sample is performed through a visual inspection of the 3GHz radio maps.This procedure aims to select all the galaxies with radio morphologies generally associated with AGN (e.g.radio blobs and relativistic jets).This visual inspection allows us to mark as likely AGN 17 galaxies.Some examples of these sources are reported in Figure 5.We underline that -although this selection is 100% pure (since non-active galaxies do not have these peculiar radio morphologies) -it is rather far from complete since most of the galaxies in the sample are unresolved at the 0.7" resolution of the 3GHz maps by Smolčić et al. 2017.Moreover, 9 RS-NIRdark galaxies are also part of the catalog by Vardoulaki et al. (2021) containing AGN in the COSMOS field selected through visual inspection thanks to their morphology in the same radio maps employed for our selection.Cross-matching our sample with this catalog, we find out that all the 9 sources are classified by AGN by both the selections.This results strengthen the reliability of our classification.

X-ray stacking
Another test to unveil the presence of AGN in our sample can be performed in the X-ray regime.We conduct two complementary tests to characterize the sample of RS-NIRdark galaxies as a whole and the individual sources.
• Sample characterization: The first test is based on the median X-ray flux computed on all the galaxies in our sample and on the possibility of explaining it by only invoking the star formation.
We follow this procedure: 1. We extract the X-ray flux for each RS-NIRdark galaxy through aperture photometry on the event file in the range [0.5,7] keV produced by the Chandra telescope in the C-COSMOS (Elvis et al. 2009) and COSMOS legacy (Civano et al. 2016) surveys.
We perform the extraction through the dmextract function of the CIAO library (Fruscione et al. 2006), employing circular apertures centered on the radio position of each source.We choose the radius of each aperture by assuming that each source is not resolved (a reasonable hypothesis for possible AGN) and employing the radius used in Elvis et al. (2009) and Civano et al. (2016) to extract the flux of the X-ray source closest to the considered galaxy.More in detail, we employ the median radius (considering all the observations in which the source fell) corresponding to an area encompassing 90% of the PSF at the X-ray closest source.The local background is computed and subtracted from each galaxy through an annulus with an outer radius 1.5 times the circular aperture.
2. Once the net counts are obtained, these are converted to X-ray luminosities through the PIMMS4 software by assuming a power-law model with a slope of Γ = 1.8 (i.e. that expected for possible AGN) and a galactic n H = 1.7 × 10 20 cm −2 (Civano et al. 2016).3. Finally, we convert the luminosities in fluxes by assuming the photometric redshifts computed by Magphys (Section 5) and these in SFR through the empiric relation by Ranalli et al. (2003).
Considering the median SFR obtained by these X-ray fluxes, we obtain a value of log(SF R X ) = (2.24 ± 0.02)M ⊙ yr −1 , slightly lower than that obtained through the FIR flux (log(SF R IR ) = (2.38±0.02)M⊙ yr −1 .This result confirms the lack of strong un-obscured AGN activity in our sample.
• Source characterization: If the analysis performed in the previous paragraph ensures that the bulk of the sample of RS-NIRdark is composed of SFGs, it does not give insights on the presence of single AGN within the total sample.To address this point, we perform an additional analysis.We cross-match the two catalogs of X-ray sources in the COSMOS field by Elvis et al. (2009) and Civano et al. (2016) with our sample by employing as matching radii the positional uncertainty included in these catalogs (for the X-ray sources) and 0.7" (for the NIRdark galaxies).Given the relatively shallow depth of the X-ray coverage in the COSMOS field, all these sources have an X-ray luminosity higher than 10 42 erg s −1 , therefore, we can assume that their X-ray flux is largely due to the presence of an AGN (see e.g.Hickox & Alexander 2018).Once the list of RS-NIRdark galaxies with a possible X-ray counterpart is obtained, we visually inspect the NIR and radio maps at 3GHz.We obtain two main cases: 1.The positional uncertainty on the X-ray source includes only the RS-NIRdark galaxy.In this case, we can safely assume that the X-ray signal is produced by the galaxy in our sample, suggesting the presence of an AGN in that galaxy.The 3 sources in this class are marked in the photometric catalog with an appropriate flag.
2. The positional uncertainty of the X-ray source includes both an RS-NIRdark galaxy and a NIR-bright galaxy.In this case, we cannot unambiguously associate the X-ray signal to one of these galaxies.The 15 NIR-dark galaxies in this class are marked in the catalog with a different flag accounting for the possibility of hosting an AGN.

q TIR analysis
A standard method to identify AGN relies on the socalled infrared-radio correlation.It is well established that the radio luminosity measured at 1.4 GHz and the infrared luminosity measured in the range [8;1000]µm are tightly (σ ∼ 0.16 dex; e.g.Molnár et al. 2021) correlated in star-forming galaxies (see e.g. de Jong et al. 1985;Helou et al. 1985).This correlation is generally measured through the q TIR parameter defined as (see e.g.Helou et al. 1985;Yun et al. 2001): and mainly arises because of the connection between the star formation, the radio synchrotron emission in star forming regions and the thermal FIR emission of dust in the same regions.Several authors studied the possible evolution of the q TIR parameter with cosmic time or possible correlations with other physical quantities.In this work, we explore two of the main studies on this point.Delhaize et al. (2017) found a possible evolution with the redshift through the relation On the contrary, more recently, Delvecchio et al. (2021) suggested a q TIR almost constant with the redshift, but strongly dependent on the stellar mass of the hosting galaxy through the relation where A = (1 + z) and B = (log(M/M ⊙ ) − 10).
We compute the radio luminosity at 1.4 GHz for the galaxies with a counterpart in the 1.4 GHz catalog by Schinnerer et al. (2010) through the relation by employing the photometric redshifts estimated by Magphys and the spectral slope α computed through the flux densities reported in the 3 GHz and 1.4 GHz catalogs.For the galaxies without a counterpart in the 1.4 GHz sample, we evaluate the 1.4 GHz flux starting from the 3GHz flux density through the relation and assuming a spectral slope of α = −0.7 (generally considered for star-forming galaxies; see e.g.Novak et al. 2017).
Once we applied Equations 4 and 5 to our sample (Figure 6), marking as likely AGN the sources distant more than 3σ from the relations (see e.g.Delvecchio et al. 2017;Enia et al. 2022), we obtain that 40 sources are classified as AGN following Delhaize et al. (2017) and 57 following Delvecchio et al. (2021).We include the results from both the tests in the final catalog.However, given the high uncertainties affecting our photo-z s and our stellar masses, in the following we consider as likely AGNs only the 37 sources classified as AGN according to both the relations.A final interesting remark concerns the overall distribution of the q TIR as a function of the photometric redshift.As visible in Figure 6, the bulk of the population of the RS-NIRdark galaxies have a median q TIR compatible with that expected at the median redshift of the sample from the relation by Delhaize et al. (2017) (2.1 ± 0.3 and 2.21 ± 0.05, respectively).However, the distribution of these values at low redshift appear to significantly differ from that expected from the relation.This discrepancy can be explained with the different selection operated in this study with respect to Delhaize et al. (2017).In particular, our selection of NIR-dark sources is expected to produce a sample of more extremely-obscured sources at low-z to account for the extinction of redder rest-frame bands.Therefore, we do not expect in this regime the properties of our (incomplete) sample to completely resemble those of the total population of radio-selected galaxies analysed in Delhaize et al. (2017).

SED decomposition
A further sign of AGN hosted in our galaxies can be found at MIR wavelengths.It arises from the presence of a dusty torus surrounding the supermassive black hole and heated up from the high-energy radiation coming from there.Trying to model the SED with a template not accounting for this additional component generally produces a best-fitting SED under-estimating all the MIR fluxes (see e.g.Hickox & Alexander 2018).This problem can be solved by adding a torus component to the templates fitted to the galaxy photometry.We perform this test with Cigale (Magphys does not allow the addition of an AGN component).The overall setup is the same as discussed in the Section 5.2, but adding the dusty torus component as modeled by Fritz et al. (2006).The model's parameters are the same as employed in Donevski et al. (2020) in modeling a sample of likely DSFGs.The main parameter describing the effect of the AGN on the modeled SED is the AGN fraction f AGN , defined as the ratio between the torus luminosity and the dust luminosity in the range [5,40]µm.Defining as likely AGN all the galaxies with a f AGN > 10%, we mark 11 sources among the entire sample of RS-NIRdark galaxies.We underline that the median value of f AGN computed on the whole sample is compatible with zero, suggesting an overall small contribution of AGN in the RS-NIRdark galaxies.However, we underline that -due to the limited coverage of the MIR regime in our sample -also this estimation should be considered as a lower limit on the actual AGN contribution in our galaxies.

Final remarks on AGN contamination
The sample of likely AGN reported by the different methods have some overlap, as shown in Figure 7. Considering all the galaxies marked as possible AGN by at least one method, we obtain a sample of 64 sources (∼ 23% of the full sample).Excluding the sources with an uncertain X-ray counterpart (i.e.those at point 2 in Section 6.2) and with no other indication of AGN activity from the other methods, we obtain 57 sources (∼ 20% of the sample).As a final remark, we can compare the estimated fraction of RS-NIRdark galaxies hosting an AGN with a theoretical prediction based on analogous works present in the literature.Bonzini et al. (2013) and Novak et al. (2018), working on radio-selected catalogs of galaxies, reported in their work the fraction of AGN in the analysis of their sources as a function of the radio flux at 1.4 GHz.Considering these relations and integrating them on our 5 flux distribution at 1.4 GHz, we estimated an expected fraction of AGN in our sample lower than 28.3% (further divided in 27% of radio-quiet AGN and 1.3% of radio-loud AGN).We underline that this last estimation represents an upper limit on the expected number of AGN in our sample, since we expect that focusing on radio sources without an optical/NIR counterpart contributes to exclude a significant fraction of these objects.

Analysis of the physical properties
The results shown in Figure 4 and Table 2, allow us to infer some general properties of the RS-NIRdark galaxies: • The high value of the median stellar attenuation A V ∼ 3.5 confirms the initial hypothesis that our selection can provide a sample of highly-obscured galaxies.
5 Here, we employ the radio fluxes at 1.4GHz provided by Schinnerer et al. (2010).For the sources without a counterpart we conservatively assume a flux equal to the sensitivity of the survey • The high value of the infrared luminosity (L IR > 10 12 L ⊙ ) and -hence -the high median SFR (∼ 300 M ⊙ yr −1 ) suggests that our galaxies are actively forming stars (as expected from a radio selection).To determine the nature of the RS-NIRdark galaxies, we need to compare the SFRs obtained through the SED-fitting with those expected from main sequence galaxies in the same redshift range (e.g.Schreiber et al. 2015).As shown in Figure 8, most galaxies lie above the main sequence line at all the redshifts.Moreover, when computing the SF R/SF R M S (i.e. the ratio between the SFR and that expected from a main sequence galaxy of the same mass and in the same redshift bin), we obtain that more than 50% of the RS-NIRdark galaxies have a SF R/SF R M S > 3, being candidate starburst galaxies.
• An additional result, strictly connected to the previous one, concerns the high median stellar mass estimated through SED-fitting (M ⋆ ∼ 10 11 M ⊙ ).This quantity -albeit quite uncertain due to the weak constraints in the optical/NIR wavelengthssuggest that the RS-NIRdark galaxies are a population of massive star-forming galaxies, with a high-z tail suitable for playing a significant role in the evolution of the massive and passive galaxies at z ∼ 3 (see e.g.Straatman et al. 2014;Schreiber et al. 2018b;Valentino et al. 2020) All these findings, once combined, confirm the initial hypothesis that the RS-NIRdark galaxies represent a significant population of high-z DSFGs.

High-z tail
The redshift distribution of the RS-NIRdark galaxies is worth a deeper discussion.The median redshift around z ∼ 3 tells us that we are looking at a population whose bulk is located at the so-called cosmic noon.However, the presence of a significant tail of high-z sources (namely, 99 galaxies at z > 3 and 17 galaxies at z > 4.5, once excluded the possible AGNs selected in Section 6.) can provide some insights on the possible evolutionary path of these sources.
The main result concerns the number density of these sources.We compute this quantity through the V max method (Schmidt 1968), considering the galaxies located at z > 3.5 and that could -therefore -play a role in the evolution of the massive galaxies at z ∼ 3.5.We account for the uncertainties in the photo-z s and in the radio fluxes through a MonteCarlo integration.Specifically, we perform a large number of realizations (∼ 500) of the total number density, sampling each time .The gray shaded area represents the 1σ = 0.3 dex scatter, while the red dotted line reports our threshold for selecting star-burst galaxies (i.e. three times the SFR expected from MS galaxies).It is possible to notice how the vast majority of the galaxies in the sample lie above the main-sequence line.The median uncertainty on the SFR and on the stellar mass is reported in the lower corner of each panel.Further details are given in Section 7.1.
and for each source a couple of values for the redshift and for the radio flux from their distribution (namely, the p(z) computed by Magphys and the values and uncertainties reported in Smolčić et al. 2017).At the end of this procedure, we consider as our number density the median value of the distribution and as its uncertainty the symmetrized interval between the 16th and the 84th percentiles.We obtain a number density of n = (3.3±0.9)×10−6 Mpc −3 for the galaxies at z > 3.5.This quantity is by a factor 6 lower than those computed by Straatman et al. (2014), Schreiber et al. (2018b), andValentino et al. (2020) for the passive and massive galaxies at z ∼ 3.5 (∼ 2 × 10 −5 Mpc −3 ).It is important to notice -however -that the our estimation must be considered as a lower limit on the actual number density of the RS-NIRdark galaxies, since it is not corrected for the expected duty cycle of the galaxies and for the incompleteness of the selection.This issue will be discussed in detail in the forthcoming paper of the series (Gentile et al., in prep).Moreover, when looking at the stellar masses of our sources, we obtain that the RS-NIRdark galaxies located at z > 3.5 have a median stellar mass of log(M ⋆ ) = 11.0M⊙ with a 1σ dispersion of 0.45 dex.
Comparing this quantity with the expected properties of the progenitors of the massive galaxies at z ∼ 3.5 (see e.g. the forecasts by Valentino et al. 2020), we can notice how the low-mass end of those objects cannot be formed by the galaxies in our sample.This result also can be used to explain the difference between the number densities of the two populations.

Comparison with the literature
An interesting final point to discuss concerns the RS-NIRdark galaxies and their possible overlap with other important populations of DSFGs: • H-dropout: The first population is composed of the H-dropout galaxies (Wang et al. 2019).These galaxies are selected as H-dark sources in the CANDELS survey (the limiting magnitude at 5σ in the H-band is 27mag; Grogin et al. 2011;Koekemoer et al. 2011) with a counterpart in the second channel of the IRAC camera from the SEDS survey ([4.5]<24mag; 80% completeness limit; Ashby et al. 2013).When considering the full photometric catalog of 323 RS-NIRdark galaxies, we obtain that only 266 sources satisfy the selection criteria exposed in Wang et al. (2019), the others being too faint at 4.5 µm to be selected with the cut on the [4.5] flux performed by Wang et al. (2019).This result indicates that a significant fraction of the RS-NIRdark galaxies would not be selected as H-dropout.On the contrary, when examining the 18 H-dropout galaxies selected by Wang et al. (2019) in the COSMOS field, we find that only 2 sources have a significant (S/N > 5.5) radio counterpart at 3 GHz in the catalog by Smolčić et al. (2017).These two results, when combined, ensures that the two selections of RS-NIRdark galaxies and H-dropout are different, with just some sources belonging to both.
• SMGs: The second population is composed of the so-called sub-millimeter galaxies (SMGs).Since this definition can be applied to all the sources de-tected in a (sub)mm survey, it strongly depends on the considered instrument's sensitivity.To obtain results comparable with those reported by Wang et al. (2019), we consider the ALESS survey (Swinbank et al. 2014) targeting bright galaxies in the SCUBA survey (i.e.galaxies with a flux density at 870µm S 870 > 4.2 mJy).We can estimate the (sub)mm flux density of the RS-NIRdark galaxies at 870µm through the best-fitting SED provided by Magphys.In doing so, we obtain that only 76 galaxies have S 870 > 4.2 mJy, being consistent with the selection by Swinbank et al. (2014).On the other side, Thomson et al. (2014) pointed out that only ∼70% of SMGs are radio-bright galaxies, while Gruppioni et al. (2020) pointed out that only a tiny percentage of SMGs are NIR-dark galaxies.As before, these results suggest that the RS-NIRdark galaxies and SMGs are two different populations of galaxies with just some sources in common.
Finally, we can compare the physical properties estimated by Wang et al. (2019) andda Cunha et al. (2015) for the population of H-dropout and SMGs with those presented in this study.Figure 9 shows the comparison between the SFR and stellar mass of the RS-NIRdark galaxies in the redshift range 3.5 < z < 4.5 (i.e.around the median redshift z ∼ 4 reported by Wang et al. 2019 for the H-dropout) with those reported by da Cunha et al. ( 2015) and Wang et al. (2019) and with the main sequence at this redshift (Schreiber et al. 2018arescaled to a Chabrier 2003 IMF).For completeness, we also include in the plot the sample of "NIR-faint" (Ks > 23.9) SMGs selected in COSMOS by Smail et al. (2021).
We can notice that the RS-NIRdark galaxies areon average -more star-forming than the H-Dropout (∆ log(SFR/M ⊙ yr −1 )∼0.6), with a median SFR comparable with the SMGs (∆ log(SFR/M ⊙ yr −1 )∼0.02), but with a significant tail of less star-forming sources.This result can be explained by the radio-selection, cutting out most of the sources below the main sequence (belonging mainly to the H-dropout).Additionally, we can see how the RS-NIRdark galaxies cover most of the mass range covered by the other selections, being on average more massive than the H-dropout and less massive than the SMGs.This result is in agreement with the discussed overlap with both the other populations.

SUMMARY
In this paper, we presented the first panchromatic study of 263 Radio-Selected NIRdark galaxies discovered in the COSMOS field following the selection by Talia et al. (2021).The development of a new deblend-ing tool (PhoEBO: Photometry Extractor for Blended Objects) allowed us to extract accurate photometry in the optical-to-MIR regime, even for the sources with a close bright contaminant.This procedure, in particular, allowed us to analyze a wider sample of galaxies missed in previous studies on the RS-NIRdark galaxies in the COSMOS field (Talia et al. 2021;Behiri et al. 2023).The complete photometric catalog6 has been employed to estimate the photo-z s and physical properties of all the galaxies in the sample through an SEDfitting procedure performed with two complementary codes (Magphys and Cigale).The results obtained with these algorithms confirmed the initial hypothesis that the RS-NIRdark galaxies are a population of starburst DSFGs, lying above the main sequence in all the redshift bins and with a significant amount of dust absorbing their optical/NIR emission.Moreover, by studying in detail the redshift distribution of the galaxies in the sample and their number density, we obtain precious insights on the possible evolutionary path of these sources, collecting significant clues that the RS-NIRdark galaxies could play a key role in the evolution of the massive and passive galaxies discovered at z ∼ 3.In addition, the analysis of the multi-wavelength counterpart in all the wavelength regimes (from the X-rays to the radio) allowed us to estimate the possible AGN contribution in our sample.Finally, through a comparison with other populations of DSFGs, we confirmed that the radio-selection produces a population of galaxies with different physical properties with respect to the SMGs (e.g.da Cunha et al. 2015) and the H-Dropout (Wang et al. 2019).
The catalog presented in this work will be employed in the next papers of the series to derive the luminosity function of the RS-NIRdark galaxies, to estimate their contribution to the cosmic SFRD and to study their role in the evolution of the massive galaxies.In the near future, the results presented in this work will be reinforced by the plethora of new data coming from state-of-theart facilities.About half of the RS-NIRdark galaxies will be observed by the James Webb Space Telescope as a part of the COSMOS-Web survey (Casey et al. 2022).Moreover, a first spectral analysis of a pilot sample of these galaxies performed with ALMA will allow a robust determination of the spectroscopic redshifts (Gentile et al, subm.).
We warmly thank the anonymous referee for his/her comments on the first version of this paper, which really allowed us to increase the overall quality of our study.FG and MT thank Eleni As previously discussed in Section 1, this paper continues the previous studies on the RS-NIRdark galaxies in the COSMOS field started by Talia et al. (2021) and Behiri et al. (2023).Here, we discuss the main differences between this work and the previous ones.The main upgrade of this work concerns the selection of the sample of RS-NIRdark galaxies.Firstly, by crossmatching the radio catalog with the NIR-selected COS-MOS2020 (in place of the COSMOS2015 employed by Talia et al. 2021 andBehiri et al. 2023), we have been able to exclude from the sample 153 NIR-faint galaxies missed in the previous version of the COSMOS catalog.These galaxies were undetected in the COSMOS2015 because of the brighter limiting magnitude reached in Table 3.Comparison between the median properties obtained through Magphys in this paper and in the previous works based on sub-samples of the RS-NIRdark galaxies (Talia et al. 2021;Behiri et al. 2023 the NIR bands (∼ 1 mag in the Ks band) and in the shallower χ 2 -image employed as the detection image (see Laigle et al. 2016;Weaver et al. 2022).Secondly, both Talia et al. (2021) and Behiri et al. (2023) employed an additional selection step after the three discussed in this paper (Section 2).These step were aimed to only analyze the sources without contamination from a nearby source.More in detail, Talia et al. (2021) removed from the sample all the galaxies with a nearby source whose 3σ isophote in the χ 2 -image overlapped the 3σ isophote in the radio map at 3 GHz.Thanks to this selection, Talia et al. (2021) analyzed a sub-sample of the full catalog composed by 197 galaxies.Similarly, Behiri et al. (2023) slightly enlarged the sample by including 75 "slightly contaminated" sources (i.e.galaxy with a NIR-bright contaminant overlapping the 5σ radio isophote).Differently from these works, in this paper we developed an ad-hoc pipeline (PhoEBO; see Section 4) to extract the photometry of all the galaxies, regardless of the possible blending with nearby contaminants.Other important differences are: • Differently from Talia et al. (2021) and Behiri et al. (2023), in this paper we employ the IRAC maps from the DAWN survey (Section 5.1).These maps are significantly deeper than those (from the SPLASH survey 7 ) used in the previous studies, allowing us to obtain more detection in the IRAC bands for our galaxies and put more stringent constraints on the upper limits for the IRACundetected sources.Similarly, the flux densities in the (sub)mm regime are obtained from the most up-to-date version of the A3COSMOS catalog, re-7 https://splash.caltech.edu/index.htmlflecting in a more significant number of counterparts for our sample of galaxies.
• Finally, In Talia et al. (2021), the photometry is extracted ex-novo only for the IRAC bands, while the flux densities and the upper limits in the optical and NIR bands are obtained from preexisting catalogs or the average depth of the analyzed maps.In this work, the pipeline presented in Section 4 allows us to perform forced photometry on the exact position of the sources (thanks to the radio prior).This procedure (analogous to that followed by Behiri et al. 2023 for the NIR-dark galaxies) enable us to pose more stringent upper limits in the optical and NIR band, reducing the degrees of freedom in the SED-fitting.
These improvements allow us to build a more complete photometric catalog employed in the SED-fitting procedure.Analyzing the results obtained with Magphys (i.e. the same software used in the previous works; see Table 3), we can notice some interesting differences: • The median redshift of the whole sample of RS-NIRdark galaxies is slightly lower (see Table 3 than that reported in Talia et al. (2021) and Behiri et al. (2023).
• The median stellar attenuation is slightly lower than that reported in the previous studies (Table 3).This result can be explained by the more stringent upper limits obtained by the newdeblending algorithm in the optical and NIR bands and in the MIR regime thanks to the deeper IRAC maps employed in this work.
The main difference between the results of the various studies concerns the number density of the RS-NIRdark galaxies.Following the procedure discussed in Section 7.2, we can compute this quantity for the galaxies located at z > 4.5 (i.e. the highest redshift bin adopted in Talia et al. 2021).We obtain n = (1.1 ± 0.6) × 10 −6 Mpc −3 .This quantity is in moderate tension with that obtained by Talia et al. (2021): n = (5.2± 1.3) × 10 −5 Mpc −3 .This result can be explained by the lower number of galaxies analysed in this study with respect to the previous one and with the lower median redshift obtained here..

B. DATA RELEASE
The full photometric catalog of the RS-NIRdark galaxies in the COSMOS field, together with the other data supporting this study can be freely accessed on the website of the collaboration 8 .More in detail, the data release consists in: • One catalog containing all the photometry from optical to radio wavelengths: The catalog has follows this structure: -Column 1: A progressive ID.For consistency with the previous studies by Talia et al. (2021) and Behiri et al. (2023), the IDs refer to the total sample of 476 RS-NIRdark galaxies selected in the COSMOS field by Talia et al. (2021).
-Columns 2-27: Fluxes and photometric uncertainties obtained through the PhoEBO for all the bands discussed in Table 1.
-Columns 28-85: Fluxes and photometric uncertainties at FIR/(sub)mm/radio wavelengths obtained through cross-matching with the pre-existing catalogs listed in Section 3.2.
All the fluxes and uncertainties are reported in Jy.
The catalog can be directly employed to perform the SED-fitting with Magphys, therefore all the entries with a S/N < 1 are treated as upper limits (i.e. with a missing flux and with a positive uncertainty equal to the 1σ upper limit).
• Two catalogs containing all the properties estimated through SED-fitting with Magphys and Cigale: A complete description of all the columns can be found in in the FITS header of each file.
• A catalog with the AGN classifications performed in Section 6: The description of the columns is contained in the FITS header of the file.
• A series of pdf files showing the best-fitting SEDs of all the RS-NIRdark galaxies as computed by Magphys and Cigale and the cutouts in most of the bands analysed with PhoEBO

Figure 1 .
Figure 1.Scheme summarizing the workflow of the deblending algorithm PhoEBO employed to extract the optical/NIR/MIR photometry for the RS-NIRdark galaxies in the sample.Further details on the numbered steps and on the full procedure are given in Section 4.1.

Figure 2 .
Figure2.Accuracy of the PhoEBO pipeline as a function of the IRAC flux of the RS-NIRdark galaxies and of the angular distance between these sources and the blended contaminant (a proxy for the blending between the sources, once parameterized through the sum of the FWHMs of the two sources).These results are obtained by applying the PhoEBO pipeline to the set of simulated images described in Section 4.2.
4. We rescale the flux of the convolved radio source to obtain an integrated flux in the range [21,M] mag, where M is the limiting magnitude in each IRAC channel.Consequently, we rescale the flux of the convolved contaminant to obtain a flux ratio uniformly sampled in the range [0.1,1] (with the radio source brighter than the other galaxy).

Figure 3 .
Figure 3.Some examples of SEDs of high-z candidates fitted with Magphys (blue solid lines; da Cunha et al. 2008; Battisti et al. 2019) and Cigale (orange solid lines; Boquien et al. 2019).The photometry is represented by the black dots (for the detections at S/N > 3) and the black arrows (for the upper limits).The yellow boxes report the photo-z s computed with the two codes and the relative uncertainties.
empirical template (see da Cunha et al. 2008 for details).The code computes the dust temperature (T d ) as the luminosity-average of these three components.Finally, Magphys includes the radio emission from star formation as prescribed by da Cunha et al. (2015).All the ranges of the free parameters employed in the templates can be found in da Cunha et al. (2008) and da Cunhaet al. (2015).In this work, we use the photo-z version of Magphys(Battisti et al. 2019), able to estimate the photometric redshift together with the physical properties of the analyzed galaxies.
The choice is limited to the models by Draine & Li 2007 and their updated version by Draine et al. 2014 parameterized through the intensity of the radiation field and including the PAHs emission.A second possibility is the analytical model by Casey

Figure 4 .
Figure 4. Distribution of the photo-z s and of the main physical properties (in order: stellar mass, star formation rate, extinction, dust luminosity, dust temperature -both luminosity-and mass-weighted -and dust mass) of the RS-NIRdark as estimated by SED-fitting with Magphys+photo-z (da Cunha et al. 2015; Battisti et al. 2019) and Cigale (Boquien et al. 2019) in blue and orange, respectively.The dashed lines of the same colors mark the median values of the distributions.Further details are given in Section 5.(2012) parameterized through the dust temperature but not including the PAHs.Since our photometric catalog includes a point at 24µm from the SuperDeblended catalog, sampling the typical PAHs emission at z ∼ 3 and since this feature is generally crucial for determining a robust redshift, we decide to employ theDraine et al. (2014) model.We compute the mass-weighted dust temperature starting from the intensity of the radiation field reported in output by Cigale following the framework described inDraine & Li (2007) andDraine et al. (2014).

Figure 5 .
Figure 5.Some examples of the sources marked as likely AGN in the final catalog because of their morphology during the visual inspection of the radio maps at 3GHz.The postage have a 15 arcsec side.Further details are given in Section 6.

Figure 6 .
Figure 6.Behavior of the qTIR as a function of the redshift (left panel) and stellar mass (right panel).The blue solid lines show the relation by Delhaize et al. (2017) and Delvecchio et al. (2021), respectively.The objects distant more than 3σ from each relation are highlighted in red.

Figure 7 .
Figure 7. Diagram reporting the number of likely AGN unveiled by the different methods described in Section 6 and the relative overlap.Further details in Section 6.

Figure 8 .
Figure 8.Comparison between the SFRs and stellar masses of the RS-NIRdark galaxies (as computed by Magphys) and those expected from the main sequence of the star-forming galaxies (red solid line; Schreiber et al. 2015 rescaled to a Chabrier 2003 IMF).The gray shaded area represents the 1σ = 0.3 dex scatter, while the red dotted line reports our threshold for selecting star-burst galaxies (i.e. three times the SFR expected from MS galaxies).It is possible to notice how the vast majority of the galaxies in the sample lie above the main-sequence line.The median uncertainty on the SFR and on the stellar mass is reported in the lower corner of each panel.Further details are given in Section 7.1.

Figure 9 .
Figure 9.Comparison between the SFRs and stellar masses of several populations of DSFGs.The different symbols show the RS-NIRdark galaxies (blue dots), the H-dropout (Wang et al. 2019; red crosses), the SMGs found in ALESS (da Cunha et al. 2015; orange triangles) and the NIR-faint SMGs discovered by Smail et al. (2021) in the UKIDSS Ultra Deep Survey (purple triangles).The solid gray line represents the main sequence at z ∼ 4 (Schreiber et al. 2015 rescaled to a Chabrier 2003 IMF) and the gray shaded area its 1σ = 0.3 dex scatter.The dotted gray line reports our threshold to define starburst galaxies (i.e. three times the SFR of a main sequence gaalxy).The median uncertainty on the SFR and on the stellar mass is reported in the lower corner of each panel.For consistency with the cited studies, we only report galaxies with 3.5 < z < 4.5.Further details are given in Section 7.1.

Table 1 .
Main properties of the maps analyzed in Section 3. Part of the data are reproduced from Weaver et al. (2022) with author's permission.

Table 2 .
Comparison between the median properties estimated by Magphys and Cigale Vardoulaki for sharing her sample of AGN in COSMOS.FG and MT thank Ian Smail for sharing his sample of NIR-faint SMGs in COSMOS.FG thanks Alberto Traina for the valuable support on the SED-fitting.FG, MT, AL, MM and AC acknowledge the support from grant PRIN MIUR 2017-20173ML3WW 001.'Opening the ALMA window on the cosmic evolution of gas, stars, and supermassive black holes'.AL is supported by the EU H2020-MSCA-ITN-2019 project 860744 'BiD4BEST: Big Data applications for Black hole Evolution STudies'.MV acknowledges financial support from the Inter-University Institute for Data Intensive Astronomy (IDIA), a partnership of the University of Cape Town, the University of Pretoria, the University of the Western Cape and the South African Radio Astronomy Observatory, and from the South African Department of Science and Innovation's )