© The Authors 2025

1. Introduction

Dust plays a central role in galaxy evolution. It functions as a catalyst in transforming atomic hydrogen into the molecular hydrogen from which stars are formed, while chemical reactions on the surface of dust grains define the structure of the interstellar medium (ISM; e.g., Hollenbach & Salpeter 1971; Mathis 1990; Wolfire et al. 1995; Weingartner & Draine 2001; Draine 2003). Dust is mainly produced in the envelopes of asymptotic giant branch stars (AGB; Gehrz 1989; Ferrarotti & Gail 2006; Sargent et al. 2010; Nanni et al. 2013; Schneider et al. 2014) and at the end of the life of massive stars during supernovae (SNe) explosions (Todini & Ferrara 2001; Rho et al. 2008; Dunne et al. 2009). On the other hand, SNe shocks destroy dust grains (Draine & Salpeter 1979; McKee et al. 1987; Jones et al. 1994; Bianchi & Ferrara 2005; Nozawa et al. 2007), which can form again in the ISM through an accretion process (Tielens 1998; Zhukovska et al. 2008). Dust can also be ejected from galaxies by powerful winds, providing an additional cooling channel and contributing to the metal enrichment of the intergalactic medium (IGM; e.g., Ostriker & Silk 1973; Bouché et al. 2007).

Dust absorbs the ultraviolet (UV) emission of young OB-type stars, further acting as a coolant to condense gas and form new stars. This radiation is reemitted at long wavelengths in the mid-infrared (MIR) and far-infrared (FIR) spectral domains up to the millimeter-wavelength regime (e.g., Puget et al. 1996; Fixsen et al. 1998; Dole et al. 2006; Driver et al. 2008). Thus, it has become clear that studying the cosmic dust budget and how it has evolved from the primordial Universe to the current epoch plays a key role in understanding the evolution of galaxies through cosmic history. So far, not many studies have been dedicated to the statistics of the dust content of galaxies. Direct measurements of the dust mass function (DMF) of galaxies are hindered by the difficulty of performing an unbiased selection on the basis of their dust mass.

The first attempts to measure the DMF in the local Universe were those by Dunne et al. (2000) and Dunne & Eales (2001), who used a sample of IRAS bright galaxies observed with the Sub-millimeter Common User Bolometer Array (SCUBA) at 450 and 850 μm, and by Vlahakis et al. (2005), who added a sample of optically selected galaxies. These early studies were extended up to redshift z ∼ 0.5 by Dunne et al. (2011), exploiting a sample of 250 μm-selected Herschel-sources with a reliable counterpart in the Sloan Digital Sky Survey catalog. It was not until the works by Driver et al. (2018) and Pozzi et al. (2020) that the study of the DMF reached redshift z ∼ 2. The former computed the DMF of more than 5 × 10⁵ optical and near-infrared (NIR) selected galaxies at z < 1.75 in Herschel extragalactic surveys. The latter studied the DMF of 160 μm selected Herschel galaxies in the redshift range from z = 0.2 to 2.5.

All the works mentioned so far based the determination of the dust mass in galaxies, the DMF, and the dust cosmic density on fitting the FIR-millimeter spectral energy distributions (SEDs) of galaxies with a model describing their dust emission. Several alternative approaches have also been employed. Fukugita (2011) computed the local dust cosmic density by integrating the star formation rate density (SFRD). Ménard & Fukugita (2012) measured the amount of dust residing in MgII absorbers along the lines of sight to distant quasars. Péroux & Howk (2020) derived the comoving dust density of galaxies from the gas cosmic density, adopting a dust-to-gas mass ratio.

An unbiased and direct dust-driven selection implies using an observable tracer univocally linked to the amount of dust in galaxies. Such an indicator must trace the bulk of the dust mass, minimizing losses or biases. A natural choice is the rest-frame sub-millimeter or millimeter emission of galaxies tracing the cold dust component that dominates their dust mass budget. So far, few studies have attempted to determine the DMF and comoving dust density in this way, especially at high redshift. Dudzevičiūtė et al. (2021) based their study on an SED fitting of 450 and 850 μm selected galaxies detected by SCUBA2 and the Atacama Large Millimeter/submillimeter Array (ALMA). Magnelli et al. (2020) studied the contribution of NIR H band selected galaxies to the cosmic dust density up to z = 3.2 by performing stacking on deep 1.2 mm ALMA maps. Pozzi et al. (2021) followed a similar approach for UV-selected galaxies at z = 4.4 − 5.9 in the ALMA/ALPINE survey (Le Fèvre et al. 2020). Recently, Eales & Ward (2024) exploited a stacking of optically selected galaxies on 850 μm SCUBA-2 maps (Millard et al. 2020) to study the evolution of dust up to z = 5.5.

In order to access the rest-frame sub-millimeter emission of dust in high-z galaxies, it is necessary to carry out observations in the millimeter wavelength domain. Traina et al. (2024a) derived the DMF of 189 galaxies from the ALMA A³ COSMOS collection (Liu et al. 2019; Adscheid et al. 2024), finding a smooth evolution of the dust mass density over the whole redshift range 0.5 < z ≤ 6.0, which is at odds with the sudden drop at z > 3 predicted by models.

The New IRAM KID Arrays 2 (NIKA2; Perotto et al. 2020; Adam et al. 2018; Calvo et al. 2016; Bourrion et al. 2016; Monfardini et al. 2014) is a dual-band millimeter continuum camera operating at 1.2 and 2.0 mm (260 and 150 GHz) installed at the IRAM 30 m telescope in Spain. It consists of three arrays of kinetic inductance detectors, two at 1.2 mm (with 1140 detectors each) and one at 2.0 mm (616 detectors), covering an effective field of view of ∼6.5 arcmin in diameter.

As part of the guaranteed time program, the NIKA2 Cosmological Legacy Survey (N2CLS) dedicated 300 hours to the observation of the GOODS-N and COSMOS fields (Dickinson & GOODS Legacy Team 2001; Scoville et al. 2007) to obtain a systematic census of dusty star forming galaxies (DSFGs) at 1.2 and 2.0 mm (Bing et al. 2023). Because of the long wavelength sampled by this blind survey, the detection of high redshift galaxies is favored (Blain et al. 2002; Lagache et al. 2005; Lutz 2014; Béthermin et al. 2015). The GOODS-N field benefits from a very rich multi-wavelength database spanning from the X-rays to radio frequencies and including observations with space telescopes such as Hubble, James Webb, Spitzer, Herschel, Chandra, and ground-based facilities such as the Gran Telescopio Canarias (GTC), Keck, the Large Binocular Telescope (LBT), Subaru, the Canada-France Hawaii Telescope (CFHT), the Kitt Peak National Observatory (KPNO), the Very Large Array (VLA), and others. We matched the NIKA2 catalog to the available multi-wavelength data, starting from high-resolution sub-millimeter and radio observations of the Sub-millimeter Array (SMA) and VLA down to NIR and optical through FIR and MIR. In case of missing high-resolution long-wavelength data or redshift estimates, in order to pinpoint the positions of the NIKA2 sources and measure their redshifts, we performed further millimeter interferometric observations with the IRAM NOrthern Extended Millimeter Array (NOEMA).

We exploited this exceptional data set, driven by the NIKA2 1.2 and 2.0 mm observations, and performed a panchromatic SED fitting, with the ultimate goal being to derive the dust mass of the N2CLS galaxies. In this way, we studied the evolution of the DMF and of the dust cosmic density from z ∼ 7 to z ∼ 1.5 in an unbiased way, as close as possible to a direct dust mass selection. This paper is organized as follows. Section 2 describes the NIKA2 and NOEMA data as well as all the ancillary data sets available in the GOODS-N region. The details about how the NIKA2 catalogs have been matched to all other multi-wavelength photometry are given in Sect. 3. In Sect. 4, the eight methods adopted to fit the SEDs of our sources are presented. A comparison between the results of the different methods is in Appendix A. The overall properties of the N2CLS GOODS-N galaxies, as found with the SED fitting, are summarized in Sect. 5. The comoving number density of these galaxies is computed in Sect. 6, and the evolution of the DMF and of the dust cosmic density are discussed in Sect. 7. Finally, Sect. 8 summarizes our findings. Throughout this manuscript, we adopt a ΛCDM cosmology with H₀ = 67.4 km s⁻¹ Mpc⁻¹, Ω_m = 0.315, and Ω_Λ = 0.685 (Planck Collaboration VI 2020), a Chabrier (2003) stellar initial mass function (IMF), and the Draine (2003) frequency-dependent dust absorption coefficient, κ_ν, renormalized as indicated by Draine et al. (2014, i.e., κ_850 μm = 0.047 m² kg⁻¹ as reference).

All the products from this paper have been released online on the survey’s home page¹. These comprise the N2CLS final maps and catalogs, the NOEMA follow-up data, and the matched catalog, which includes the identification, multi-wavelength counterparts, and redshifts.

2. N2CLS and multi-wavelength data

The Great Observatories Origins Deep Survey Northern field (GOODS-N, RA 12:36:55, Dec +62:14:19; Dickinson & GOODS Legacy Team 2001) is centered on the now-historic Hubble Deep Field North (HDF; Williams et al. 1996) and benefits from a very rich multi-wavelength coverage, ranging from the X-rays to radio frequencies, including observations carried out with all major facilities in the northern hemisphere and in space. This section presents an overview of all data sets that have been considered when building the SEDs of the N2CLS sources in GOODS-N, starting from the driving sample detected by NIKA2 at the IRAM 30 m telescope and by NOEMA.

2.1. N2CLS observations and source extraction

The N2CLS observations of the GOODS-N field (∼160 arcmin²) were carried out with NIKA2 between October 2017 and January 2023, under project ID 192-16, for a total of 86.15 hours of telescope time (Bing et al. 2023). NIKA2 observes simultaneously in two photometric bands at wavelengths of 1.2 and 2.0 mm, with half power beam widths (HPBW) of ∼12 and 18 arcsec, respectively. The data have been reduced with the IRAM PIIC software² (Zylka 2013; Berta & Zylka 2019-2024) following the standard iterative procedure for deep fields. We defer to Bing et al. (2023) for a detailed description of the data set and its reduction.

The final NIKA2 maps of GOODS-N reach noise levels of 0.11 mJy/beam, and 0.031 mJy/beam in the deepest, central area, at 1.2 and 2.0 mm, respectively. At these depths, NIKA2 hits the photometric confusion limit at 2.0 mm and reaches within a factor of 2 from the confusion limit at 1.2 mm (Ponthieu et al., in prep.). The noise in the map is not uniform and increases toward the outer regions of the field.

Source extraction was performed with a dedicated software on the matched-filter PIIC maps, using circular 2D Gaussian kernels matching the HPBW. Sources were identified as peaks on the signal to noise ratio (S/N) maps and measured using PSF fitting on the signal maps (Bing et al. 2023). Possible systematic effects introduced by the data reduction were evaluated with simulations, by injecting an artificial sky model into the NIKA2 data timelines and performing the data reduction in the same way as for the original data. In this way the purity and photometric completeness of the extracted catalog were determined. The GOODS-N catalog reaches a purity of 80% at S/N = 3.0 and 2.9, at 1.2 and 2.0 mm, respectively, and > 95% at S/N > 4.2 and > 4.1. This procedure allowed us also to correct the effect of flux boosting/filtering, that is the impact of noise and data reduction on the measured source fluxes, as well as to determine the “effective area” of the survey “seen” by each detected object, given the non-uniform noise level across the maps.

2.2. N2CLS sample selection

The sample analyzed in this work consists of the N2CLS sources extracted from the NIKA2 1.2 and 2.0 mm GOODS-N maps, with S/N ≥ 4.2 in at least one NIKA2 band. This translates in a flux cut of ∼0.7 mJy at 1.2 mm, with a few sources that are slightly fainter because lying in the deepest central region of the maps. In the statistical analysis that follows, the non-uniform coverage and noise level are taken into account in a natural way by using the effective area mentioned above.

The NIKA2 sample analyzed here includes a total of 65 sources selected with S/N ≥ 4.2 at either 1.2 or 2.0 mm (corresponding to > 95% purity; Bing et al. 2023). Out of these, 63 are detected at 1.2 mm and 26 at 2.0 mm, with S/N ≥ 4.2. Two sources are detected only at 2.0 mm. We complement the flux at either 1.2 or 2 mm with our blind catalog of S/N > 3 detected sources. Table 1 summarizes the number of sources and the number of matches to all multi-wavelength catalogs considered here. In what follows, we name this sample and catalog “N2GN”.

2.3. NOEMA follow-up observations

With the current sensitivity, N2CLS is constructing the most complete sample of high-z IR-luminous massive galaxies in GOODS-N. However, the angular resolution of NIKA2 at 1.2 and 2 mm makes it difficult to unambiguously identify the multi-wavelength counterparts of the N2CLS sources. As already demonstrated by the follow-ups of SCUBA-2 sources with the SMA (e.g., Cowie et al. 2017) or ALMA (e.g., Stach et al. 2019), the combination of (sub-)mm single-dish and interferometer surveys is by far the most efficient way of constraining the dusty star formation at z > 2. The high resolution and sensitivity of NOEMA were thus used to provide accurate position measurements on N2CLS sources that had ambiguous proxy for multi-wavelength identification or no counterpart (Sect. 3).

We observed 27 N2GN sources (program W21CV, P.I. L. Bing, for 39.5 hours in C-configuration) either at 255 GHz (23 sources) or 150 GHz (4 sources) in the continuum, matching the frequencies of the two NIKA2 bands. We reached beam sizes of 0.84″ × 0.61″ and 1.61″ × 0.93″ at 255 and 150 GHz, respectively. The observations reveal at least one counterpart at S/N > 4 for all but one source (N2GN_1_59 which has a NOEMA source at the very edge of the beam, at 9.7″ from the NIKA2 position). Five N2CLS sources break into two millimeter sources (N2GN_1_17, N2GN_1_24, N2GN_1_27, N2GN_1_34 and N2GN_1_56). Our NOEMA data for N2GN_1_05 also reveal two broad CO lines, giving a redshift similar to that obtained by 3D-HST (z = 1.996). Figure 1 shows the resulting continuum maps of two N2GN example sources.

Fig. 1. NOEMA 150 GHz maps for N2GN_1_13 (left) and N2GN_1_17 (right) that reveal one and two counterparts, respectively. The beam shown at the bottom left has a size of 1.61″ × 0.93″. Multi-wavelength postage stamps, including NIKA2, are shown in Fig. E.1.

The NOEMA data were reduced and calibrated in the standard way with the GILDAS³ software, producing uv-tables of the lower and upper sidebands (LSB, USB). The main calibrators adopted were MWC349 and LkHα101. The absolute flux uncertainty is 10% and the positional error is 0.2 arcsec.

All channels of LSB and USB were combined together, taking care to exclude possible spectral lines, into double side band (DSB) continuum maps. These maps were cleaned using natural weights and support masks including all detected sources. Aperture fluxes were measured using ad hoc extraction polygons, taking into account primary beam losses. Flux statistical uncertainties were measured by rescaling the noise to the same apertures and finally source positions were computed as barycenters of the signal within the polygonal apertures. Figure 2 compares the 1.2 and 2.0 mm flux densities measured on the NIKA2 and NOEMA maps. No bandpass correction has been applied. In case of multiple sources, we have summed the contribution of the different NOEMA components.

Fig. 2. Comparison between NIKA2 and NOEMA flux densities for the sources observed in program W21CV. The green solid diagonal lines marks the 1:1 correspondence between the two instruments.

Furthermore, as a N2CLS pilot project, we previously observed with the NIKA pathfinder (Adam et al. 2014; Catalano et al. 2014) a small 3.5′ × 2.5′ field (proposal 230-14) centered at (RA, Dec) = (12:36:27, 62:12:18) on AzGN10 (Penner et al. 2011), a z ∼ 6.5 candidate galaxy (Liu et al. 2018). We followed up with the Plateau de Bure Interferometer (PdBI) two sources that are now included in N2CLS: N2GN_1_09 and N2GN_1_12. These two sources were observed using the WideX correlator in band3 (255 GHz) with the D configuration (project W16EE). We further followed-up N2GN_1_09 at 3 mm with the D configuration (project E16AI) which provided a continuum flux at 109 GHz (Bing 2022). N2GN_1_12 has two interferometric counterparts. Summarizing, our N2CLS selected sample comprises 71 individual millimeter galaxies.

2.4. Other (sub-)millimeter data

Further NOEMA observations by other teams allow us to get an accurate position and/or complementary data for N2GN_1_01 (GN10; Riechers et al. 2020), N2GN_1_04 (GN20; Daddi et al. 2009), N2GN_1_06 (HDF850.1; Neri et al. 2014), N2GN_1_09 (ID12646; Jin et al. 2022), and N2GN_1_44 (z = 7.2; Fujimoto et al. 2022). We also use the SCUBA-2 850 μm source catalog by Cowie et al. (2017). This catalog comprises 186 sources with S/N ≥ 4, over an area of 450 arcmin². In the central region, the 850 μm observations cover the field to near the confusion limit of ∼1.65 mJy, while over the wider region, they have a median 4σ limit of 3.5 mJy. The James Clerk Maxwell Telescope (JCMT) at 850 μm has an angular resolution of ∼14″ FWHM.

To find accurate positions and identify the true optical, NIR, and MIR counterparts of the SCUBA-2 sources, Cowie et al. (2017) conducted interferometric follow-up observations with the SMA at 860 μm. They observed nearly all of the bright SCUBA-2 sources with the SMA. Including archival data from other SMA programs, they have identifications for 33 SCUBA-2 sources with the SMA at 4σ. All sources but GN20 were observed in the compact mode of the SMA with a spatial resolution of 2″ at 860 μm.

Barger et al. (2022) presented the deepest SCUBA-2 450 μm data, achieving a central rms of 1.14 mJy. They detect 79 sources with an S/N above four, on 175 arcmin². The JCMT at 450 μm has an angular resolution of 7.5″ FWHM.

2.5. Radio data

Owen (2018) obtained new wide-band continuum observations in the 1−2 GHz band using the Karl G. Jansky Very Large Array. The best image with an effective frequency of 1525 MHz reaches an rms noise in the field center of 2.2 μJy, with 1.6″ angular resolution. We used their catalog, that contains 795 sources with an S/N = 5 or greater, covering a radius of 9 arcmin centered near the nominal center for the GOODS-N field. This area covers all our N2CLS sources.

2.6. Mid- and far-infrared data

2.6.1. Herschel:

The Herschel space telescope (Pilbratt et al. 2010) observed the GOODS-N field with the Photodetector Array Camera and Spectrometer (PACS; Poglitsch et al. 2010) and the Spectral and Photometric Imaging REceiver (SPIRE; Griffin et al. 2010). The field was included in the Guaranteed Time surveys PACS Evolutionary Probe (PEP, in the 100 and 160 μm bands; Lutz et al. 2011) and Herschel Multi-tiered Extragalactic Survey (HerMES, at 250, 350 and 500 μm; Oliver et al. 2012). The open time GOODS-Herschel survey (Elbaz et al. 2011) reached the confusion limit in GOODS-N at 160 μm. In this work, we make use of the combined PEP plus GOODS-Herschel PACS maps and catalogs by Magnelli et al. (2013) and the GOODS-Herschel SPIRE data from Elbaz et al. (2011).

2.6.2. Spitzer:

The GOODS-N field has been observed at MIR wavelengths with Spitzer’s Multiband Imaging Photometer (MIPS) at 24 and 70 μm as part of the Guaranteed Time Observers (GTO) and GOODS surveys (Dickinson et al. 2003a; Frayer et al. 2006). We retrieved the MIPS data from the Barro et al. (2019) multi-wavelength catalog. They use the photometric catalogs in both bands described in Pérez-González et al. (2005, 2008), which are based on the reduced and mosaicked data.

2.6.3. Deblended catalog:

Source confusion and clustering can introduce substantial biases in photometric works (e.g., Béthermin et al. 2017). Herschel SPIRE images have point spread functions (PSFs) that are several times larger (17.6−35.2″) than those of PACS images (7−12″), which are also several times higher than optical and NIR images. The fluxes of individual galaxies are often difficult to measure because of blending from close neighbors. When appropriate (i.e., when the N2CLS counterpart is clearly blended and we visually saw from the cutouts that an accurate flux could be measured using deblending techniques), we used the flux from the super-deblended catalog by Liu et al. (2018). Finally, the super-deblended catalog provides the additional 16 μm flux of ten sources in our N2CLS selected sample. The photometry at 16 μm was measured from the Spitzer IRS peak-up image by Teplitz et al. (2011).

2.7. Optical and near-infrared data

We used the multi-wavelength catalog from Barro et al. (2019), which is selected in the WFC3 F160W (H-band) in the GOODS-N CANDELS (Cosmic Assembly Near-IR Deep Extragalactic Legacy Survey) field. The multi-wavelength photometry includes broad-band data from UV (U-band from KPNO and LBT), optical (Hubble Space Telescope, Advanced Camera for Surveys, HST/ACS F435W, F606W, F775W, F814W, and F850LP), NIR-to-MIR (HST/WFC3 F105W, F125W, F140W and F160W, Subaru/MOIRCS Ks, CFHT/Megacam K, and Spitzer InfraRed Array Camera, IRAC 3.6, 4.5, 5.8, 8.0 μm). Similarly to the case of MIR to FIR data, for the K and IRAC bands we used instead fluxes from the super-deblended catalog by Liu et al. (2018), when appropriate.

2.8. JWST photometric data

We used the available public imaging of the James Webb Space Telescope (JWST) in the GOODS-N field (i.e., the mosaics v7.3 obtained from the DAWN JWST Archive, DJA⁴), including the surveys of FRESCO (First Reionization Epoch Spectroscopically Complete Observations; Oesch et al. 2023), JADES (JWST Advanced Deep Extragalactic Survey; Eisenstein et al. 2023), PANORAMIC (Williams et al., in prep.), and Congress (Egami et al., in prep.). From shortest to longest wavelengths, we used the NIRCam (Near Infrared Camera) and MIRI (Mid-InfraRed Instrument) filters F090W, F115W, F150W, F182M, F200W, F210M, F277W, F335M, F356W, F410M, F444W, and F770W. The photometry measurements are similar to those described by Weibel et al. (2024) and Xiao et al. (2023). Specifically, we used SExtractor (Bertin & Arnouts 1996) in dual-image mode and the longest wavelength NIRCam filter, F444W, as the detection image. The fluxes were measured in a circular aperture of radius $0\stackrel{″}{.}35$. Prior to flux measurement, all images were co-aligned and drizzled to a common grid of 40 mas/pixel. All NIRCam filters were PSF-matched to the F444W band. For MIRI/F770W, whose PSF is broader than the NIRCam/F444W PSF, we further computed the matching kernel from F444W to F770W and produced PSF-matched flux and rms images. We measured the flux of F770W based on the original F444W map and in accordance with the PSF-matched color between F770W and PSF-matched F444W. Finally, we performed an aperture correction based on the flux measured through the Kron aperture in F444W, and scaled the fluxes to total by computing the encircled energy of the Kron ellipse on the F444W PSF.

The circular aperture of radius $0\stackrel{″}{.}35$ is appropriate for the typical population of massive dusty sources at high redshift (z > 3). We observe that some of our z ∼ 2 dusty galaxies are larger, with very complex and highly wavelength-dependent morphologies. For these kinds of sources, we do not consider our JWST photometry and rather rely on previous HST+IRAC published catalogs. For HDF850.1 we adopted the fluxes published in Sun et al. (2024). In summary, we complemented our multi-wavelength SEDs with JWST fluxes for 24 sources, including 10 optically dark galaxies (notice that for N2GN_1_61 we only have one photometric data point, F444W, from Sun et al. 2024).

2.9. X-ray data

GOODS-N is included in the 2 Ms Chandra observations of the Chandra Deep Field North (CDFN), in the soft and hard X-ray bands (0.5 − 2.0 and 2.0 − 7.0 keV). Alexander et al. (2003) published the catalog of CDFN, including a total of 503 X-ray sources, 348 of which lie in the IRAC GOODS-N area (Rovilos et al. 2010).

The GOODS-N matched catalog by Barro et al. (2019) includes the X-ray ID in the Alexander et al. (2003) catalog. In addition, we also used the Chandra Source Catalog (CSC 2.0; Evans et al. 2024) and the NASA Extragalactic Database (NED⁵) to find additional X-ray counterparts of the N2GN sources and identify potential active galactic nuclei (AGN). As a result of this search, 17 N2GN sources have a X-ray counterpart (Evans et al. 2024; Barro et al. 2019; Wang et al. 2016a; Alexander et al. 2003; Brandt et al. 2001). Out of these, only six require an AGN-torus component to reproduce their observed optical-to-radio SEDs (Sect. 4), namely N2GN_1_03, 10, 11, 19, 32, and 46. Appendix E gives more details about the sources.

3. Identification of N2CLS sources

Due to the relatively low angular resolution of NIKA2, it is challenging to confidently match our N2CLS galaxies with any counterparts in the rich ancillary GOODS-N data set. This is a notoriously difficult problem for (sub-)millimeter galaxies found with blind single-dish surveys. In this Section, we describe the process of identification and the procedure followed to build the UV-to-radio SEDs of the N2GN sources.

3.1. Identification process

Based on the empirical relation observed between the FIR and radio emission in star-forming galaxies, the radio emission has often been used as a high-angular resolution proxy to get the position of the rest-frame FIR emission observed in the (sub-)mm (Smail et al. 2000; Ivison et al. 2002; Borys et al. 2004). The radio proxy has also been combined with MIR imaging from Spitzer (e.g., Ivison et al. 2007). However, when observed at these radio or 24 μm wavelengths, galaxies are subject to a positive k-correction and therefore become fainter with increasing redshift. This technique thus biases the counterpart identification toward lower redshift (e.g., Chapman et al. 2005).

As demonstrated by the follow-up observations of SCUBA-2 (e.g., Cowie et al. 2017; Simpson et al. 2017, 2020; Stach et al. 2019) or South Pole Telescope (SPT; e.g., Hezaveh et al. 2013; Reuter et al. 2020) sources with the SMA or ALMA, the synergy between single-dish and interferometer surveys at similar wavelengths is the most effective approach to obtain large unbiased samples of dusty star-forming galaxies, pinpoint their position and obtain their redshift. Following these previous studies, we located the precise positions of our galaxies using, by order of priority:

• Known millimeter counterparts from previous PdBI (e.g., HDF850.1) or NOEMA observations: six sources.

• Our own NOEMA follow-up observations (see Sect. 2.3): 31 sources.

• The SMA follow-up observations of SCUBA-2 sources from Cowie et al. (2017): 15 sources.

• Two forced identifications: one with a high-redshift optical galaxy revealed by JWST and one with the closest VLA counterpart to the N2CLS position⁶.

• The very deep VLA radio catalog by Owen (2018): 17 sources.

The list of sources with their proxy for identification is given in Table E.1. In the left panel of Fig. 3, we show the distribution of the distances between the N2CLS coordinates and the precise position of the matched proxy, as obtained through our identification process. Given the small angles involved, the total distance is computed with the Pythagorean theorem. The median total distance is 1.56″. Of the ten sources with a distance higher than 4″, six are double, three are located with VLA, and one is located with NOEMA.

Fig. 3. Differences between the coordinates of the N2GN sources and their counterparts. Left panel: Distance between the NIKA2 coordinates (all from the 1.2 mm position but two) and the matched proxy (Sect. 3.1). Right panel: Distance between the matched proxy and the N2GN optical counterparts. The green dashed lines and boxes mark the rms of the distribution.

When a NIKA2 source has multiple associations in high-resolution NOEMA data, all such associated objects are used. Their low spatial resolution, blended photometry (e.g. NIKA2 and Herschel) is excluded and different SEDs are built for each identified source. None of the SMA identifications is multiple. Finally, in the case of radio identifications, 16 out of the 17 NIKA2 sources identified through VLA data have only one radio counterpart. The remaining one, N2GN_1_55, turns out to be an optically dark galaxy with no redshift information yet (Appendix E).

3.2. Multi-wavelength spectral energy distribution

We used the multi-wavelength data presented in Sect. 2 to build the SED of each galaxy. We automatically searched for multi-wavelength counterparts at a given distance d given by

$$\begin{array}{c}\hfill d=\sqrt{{{\mathrm{HWHM}}_{\mathrm{proxy}}}^2+{{\mathrm{HWHM}}_{\lambda }}^2},\end{array}$$(1)

where HWHM_proxy is the half width at half maximum (HWHM) of the beam of the proxy used to get the precise position (e.g., NOEMA) and HWHM_λ is the HWHM of the beam of the multi-wavelength observations. The right-hand panel of Fig. 3 presents the distances between the matched proxy (Sect. 3.1) and the optical identification of the N2GN sources. The median distance is 0.25″. Five sources out of 71 have a double entry in the optical catalogs when applying Eq. (1): N2GN_1_05, 15, 24a, 26, and 32. However, the two optical counterparts correspond to: (i) the same galaxy in the first two cases; (ii) a very-close merger in the second two cases. For the last source (N2GN_1_32), we chose the closest optical counterpart (0.22 arcsec against 0.93 arcsec), which is also the brightest one.

Using the multi-wavelength cutouts, we checked one by one the multi-wavelength counterparts and removed the photometric point for blended sources (e.g., SPIRE “blobs”) or replace their fluxes with the deblended photometry by Liu et al. (2018), when appropriate. The SEDs of all individual galaxies are shown in Appendix E. The numbers of the multi-wavelength counterparts used for the SED fitting are given in Table 1.

3.3. Redshifts

We gathered all the redshift information found in the dedicated surveys and catalogs from the literature. These are from Steidel et al. (2003), Skelton et al. (2014), Bouwens et al. (2015), Cowie et al. (2017), Arrabal Haro et al. (2018), Owen (2018), Barro et al. (2019), Kodra et al. (2023). We complemented this survey approach with a dedicated search for individual sources, to obtain the redshifts measured from e.g., PdBI or NOEMA, VLA, or JWST observations. We also used the NED to check some individual objects.

For four sources (N2GN_1_08, N2GN_1_12_a, N2GN_1_16, N2GN_1_36), we found some discrepancies between different photometric redshifts from the literature. We also have three sources (N2GN_1_18, N2GN_1_34_b, and N2GN_1_43) that previously lacked known redshift values but now have sufficient photometric data points to determine a redshift. For these seven sources, we used The Code Investigating GALaxy Emission (CIGALE) software (Burgarella et al. 2005; Noll et al. 2009; Boquien et al. 2019) to determine a photometric redshift.

The redshift of each source is given in Table E.1. We have spectroscopic redshifts for 29 sources (∼43% of the sample) and photometric redshift for 39 sources. Three galaxies (N2GN_1_17_b, N2GN_1_34_a and N2GN_1_55), identified with NOEMA (the first two) and VLA (the third), have no redshift.

Figure 4 presents the redshift distribution of the N2GN sources, distinguishing between spectroscopic (solid histogram) and photometric redshifts. The median redshift of the N2GN sample is z = 2.819, with a median absolute deviation (M.A.D.) of 0.831. Compared to the redshift distribution expected for the N2CLS galaxies from the Simulated Infrared Dusty Extragalactic Sky (SIDES) simulations (Béthermin et al. 2017, and in prep.), we clearly see an excess of galaxies at z ∼ 5.1 − 5.3, linked to the complex overdense environment hosting the famous HDF850.1 (N2GN_1_06) and GN10 (N2GN_1_01) sub-millimeter galaxies (Sun et al. 2024; Xiao et al. 2023). A dedicated analysis of the N2GN galaxies in this overdensity will be presented in Lagache et al. (in prep.).

Fig. 4. Redshift distribution of the N2GN sources. The solid histogram includes spectroscopic redshifts, while the open histogram contains all sources in the sample.

4. Fitting the spectral energy distributions

The SEDs of the N2GN sources consist of up to 34 broad-band photometric measurements, spanning over the wavelength range from 3600 Å to 21 cm. We reproduced the SEDs with different modeling approaches, with the main goal to derive as robust dust properties as possible.

The observed FIR to millimeter SEDs were studied with a modified blackbody (MBB; Sect. 4.1) in the optically thin approximation and in its general form, as well as with the Draine & Li (2007, DL07, Sect. 4.2) model. The panchromatic optical to radio SEDs were modeled with the MAGPHYS and SED3FIT codes, in their original and high-z versions (da Cunha et al. 2008, 2015; Battisti et al. 2020; Berta et al. 2013a, Sects. 4.3 and 4.5), as well as with CIGALE (Burgarella et al. 2005; Noll et al. 2009; Boquien et al. 2019, Sect. 4.4).

Appendix A compares the results of the eight SED fitting methods presented here, in terms of the derived M_dust, L_IR and M_⋆ of the N2GN sources. The different approaches lead to consistent results, within their limitations. In the main analysis of this work we adopt the results obtained with MAGPHYS high-z for purely star-forming galaxies and with SED3FIT-hz in case an additional component (AGN or a young stellar population) is needed to reproduce the observed SED.

4.1. Modified blackbody

The emitted FIR-to-millimeter dust SED of a galaxy can be approximated with a single-temperature MBB, described as

$$\begin{array}{c}\hfill L=\mathrm{\Omega }{ϵ}_{\nu }{B}_{\nu }(T_{\mathrm{dust}}),\end{array}$$(2)

where Ω is the solid angle of emission, ϵ_ν is the dust emissivity coefficient and B_ν(T_dust) is the Planck function of temperature T_dust. For a uniform medium of optical depth τ_ν, radiative transfer theory leads to ϵ_ν = 1 − exp(−τ_ν). Relating τ_ν to the mass absorption coefficient κ_ν and assuming that the dusty medium is optically thin, the broadly used expression of the MBB is obtained (see Berta et al. 2016, for a thorough derivation):

$$\begin{array}{c}\hfill L_{\nu }=4\pi M_{\mathrm{dust}}{\kappa }_{\nu }{B}_{\nu }(T_{\mathrm{dust}}),\end{array}$$(3)

which is used to model the observed SEDs. The absorption coefficient depends on frequency as a power law: κ_ν = κ₀(ν/ν₀)^β. We adopt the values of κ_ν tabulated by Draine (2003), with the correction indicated by Draine et al. (2014). In this assumption, the reference value is κ_850 μm = 0.047 m² kg⁻¹.

The MBB-thin fit includes three free parameters, namely normalization, T_dust and β. The dust mass of the given galaxy, M_dust, was computed by inverting Eq. (3) and evaluating the model at the longest wavelength covered by the observed SED, such to minimize the effect of possible optically thick dust at shorter wavelengths and to avoid extrapolations beyond the available data. In computing M_dust, we took care to account for the correction described by Bianchi (2013) and Berta et al. (2016), related to the fact that the best fit value of β is not necessarily the same as the one in the tabulated κ_ν. We defer to these works for more details.

The general form of the MBB differs from the optically thin approximation by the factor 1 − exp(−τ_ν), which can be expressed in terms of the wavelength λ_thick at which the medium becomes optically thick or in terms of the size of the absorbing medium (Ismail et al. 2023). Here we adopted the former approach: a fourth additional free parameter is hence λ_thick.

Ismail et al. (2023) modeled the FIR-to-millimeter SEDs of a sample of 125 high-z dusty galaxies with an exquisite wavelength coverage with up to 12 bands between 250 μm and 3.5 mm in the observed frame (Cox et al. 2023). They studied the general MBB in two different formulations: one adopting λ_thick as free parameter; the other based on the size A of the dust emitting region. The two are linked by the expression λ_thick = λ₀(κ₀M_dust/A)^1/β. With the use of mock simulated SEDs, Ismail et al. (2023) demonstrated that, adopting the first formulation, λ_thick is hardly constrained. Furthermore, these authors showed that a reliable estimate of T_dust, β and M_dust with the general MBB requires an independent knowledge of A. The degeneracies that affect λ_thick can induce a significant underestimation of M_dust of 20% to 50%, depending on dust mass itself.

The MBB fit was limited to data at rest-frame wavelengths ≥50 μm, because below this limit the contribution of warm dust to the SED becomes non-negligible, and is not included in the MBB models. The effect of the cosmic microwave background (CMB) was taken into account as described by da Cunha et al. (2013). Radio data were not included in the fit, but they were compared to a synchrotron model obtained with a simple power law S(ν)∝ν^α (spectral index α = −0.8) normalized such to obey the radio-FIR correlation for star formation (Magnelli et al. 2015; Delhaize et al. 2017). Sources that show a significant radio excess with respect to the radio-FIR correlation are very likely to host a radio-AGN.

4.2. The Draine & Li (2007) dust model

The Draine & Li (2007) models are an upgrade of those originally developed by Draine & Li (2001), Li & Draine (2001) and Weingartner & Draine (2001). Interstellar dust is described as a mixture of carbonaceous and amorphous silicate grains, whose size distributions are chosen to mimic the observed extinction law in the Milky Way (MW), Large Magellanic Cloud (LMC) or Small Magellanic Cloud (SMC) bar region.

The dust distribution is divided in two components: the diffuse ISM, responsible for the bulk of the dust mass budget; and dust enclosed in photo-dissociation regions (PDRs). The former is heated by a radiation field of constant intensity U_min. The latter, representing a fraction γ of the total amount of dust, is exposed to starlight with intensities within the range U_min to U_max. Although PDRs usually provide a small fraction of the total dust mass, they can contribute to the majority of the MIR dust emission. The properties of grains are parameterized by the index q_PAH, defined as the fraction of dust mass in the form of PAH (polycyclic aromatic hydrocarbons) molecules. Finally the fraction of dust in PDRs heated by starlight with an intensity U is a power law of U with index −α.

The DL07 model thus has six free parameters: U_min, U_max, γ. q_PAH, α, and M_dust. Studying local Spitzer galaxies, Draine et al. (2007) demonstrated that some of the parameters can be limited to a restricted range of values. We adopted these as commonly done in the literature (e.g., Magdis et al. 2012; Magnelli et al. 2012; Santini et al. 2014; Berta et al. 2016). Also in the case of DL07 modeling, we adopted the renormalization of κ_ν prescribed by Draine et al. (2014).

4.3. MAGPHYS

The Multi-wavelength Analysis of Galaxy Physical Properties (MAGPHYS) software reproduces the observed SEDs of galaxies linking together their optical stellar emission to their dust component through energy balance (no radiation transfer involved). The energy absorbed by dust in stellar birth clouds and in the ISM is re-distributed to the dust emission at infrared wavelengths. Here we used two different versions: the original by da Cunha et al. (2008); and one dedicated to z > 1 galaxies (named “high-z v2”; da Cunha et al. 2015; Battisti et al. 2020).

The code combines the Bruzual & Charlot (2003, BC03) optical/NIR stellar library, with the MIR/FIR dust emission computed by da Cunha et al. (2008). The adopted star formation history (SFH) is a continuous delayed exponential function of the form $\psi (t)\propto \frac{t}{{\tau }^2}exp(-t/\tau )$, where t is the model age and τ the star formation timescale in units of Gyr. Superimposed to this continuous SFH are bursts of star formation of random duration and age. Dust attenuation is described with the recipe by Charlot & Fall (2000). The main assumptions are that the energy re-radiated by dust is equal to that absorbed, and that starlight is the only significant source of dust heating.

The power re-radiated by dust in stellar birth clouds is computed as the sum of three components: PAHs; a MIR continuum describing the emission of hot grains with temperatures T = 130 to 250 K; and grains in thermal equilibrium with T = 30 to 60 K. The “ambient” ISM is modeled by fixing the relative proportions of these three components to reproduce the cirrus emission of the Milky Way, and adding a component of cold grains in thermal equilibrium, with adjustable temperature in the range T = 15 to 25 K.

Different combinations of star formation histories, metallicities and dust content can lead to similar amounts of energy absorbed by dust in the stellar birth clouds, and these energies can be distributed in wavelength using different combinations of dust parameters. Consequently, in the process of fitting, a wide range of optical models is associated with a wide range of infrared spectra. We defer to da Cunha et al. (2008) for a formal description of how the SED models are built.

In the high-z version of the code, da Cunha et al. (2015) extended the parameter space of the models, such that they include higher dust optical depths, higher SFRs, and younger ages. Also the dust properties have been modified to allow for higher values of dust attenuation both in the stellar birth clouds and in the diffuse ISM. Moreover, the high-z MAGPHYS takes into account the UV absorption by the IGM including Lyman series line blanketing and Lyman-continuum absorption.

At radio frequencies, a thermal and a non-thermal components have been added to the model. The former (Bremsstrahlung, or free-free) is a power-law L ∝ ν^−0.1, normalized such to have a fixed contribution of 10% to the rest-frame 20 cm emission. The latter (synchrotron) is a power-law with spectral index −0.8, linked to the FIR emission by the local radio-FIR correlation (q_FIR = 2.34 with a 1σ dispersion of 0.25).

Finally, Battisti et al. (2020) introduced a flexible 2175 Å absorption feature to the diffuse ISM attenuation curve, with strength depending on the actual A_V/E(B − V) of the model. These authors also upgraded the IGM absorption from the Madau (1995) to the Inoue et al. (2014) prescription.

The code produces the marginalized probability distribution of the derived quantities. To our aims, among the many products available, we focus the analysis on dust mass, M_dust, infrared luminosity, L_IR, and stellar mass, M_⋆.

4.4. CIGALE

The CIGALE software (Burgarella et al. 2005; Noll et al. 2009; Boquien et al. 2019) is a public and versatile code that models galaxy SEDs from X-rays to radio taking into account the balance between the energy absorbed by dust in the UV-optical and reemitted in the IR. The code can be used to fit photometry and spectroscopy, or to create mock SEDs thanks to its large library of models. CIGALE allows for different kinds of star formation histories (SFH, parametric and non-parametric), stellar population models (e.g., Bruzual & Charlot 2003; Maraston 2005), dust attenuation laws (e.g., Calzetti et al. 2000; Charlot & Fall 2000), IR emission models (e.g., Draine & Li 2007; Dale et al. 2014), AGN components (e.g., Fritz et al. 2006; Ciesla et al. 2015), X-rays and radio emission.

SED fitting was performed on the N2GN sources adopting the Bruzual & Charlot (2003) stellar models, with a Chabrier (2003) IMF and solar metallicity, the Charlot & Fall (2000) attenuation law, the Draine & Li (2007) dust emission models, renormalized as in Draine et al. (2014), and Stalevski et al. (2016, 2012) AGN models when hinted by the X-ray, UV, MIR or radio emission. A delayed and truncated SFH was used, with e-folding time of the main stellar population of 0.1, 0.5, 1, 5, 10 Gyr, and truncation ages of 50 or 100 Myr.

4.5. SED3FIT

Originally inspired by MAGPHYS, SED3FIT (Berta et al. 2013a) combines three emission components simultaneously to reproduce the observed SEDs of galaxies. It adds a third component to the MAGPHYS stellar and dust models (Bruzual & Charlot 2003; da Cunha et al. 2008; Charlot & Fall 2000). The third component can be an AGN-torus model or any other kind of emission that might need to be considered in addition to the stellar+dust model. Also SED3FIT exists in the original version by Berta et al. (2013a) and in a new high-z implementation (named SED3FIT-hz).

Due to the huge number of possible model combinations that arise from adding a third library to the already big stellar and dust collections, as well as from the free normalization of the stars+dust component, the code samples the parameters space and the libraries of models randomly instead of systematically test all possibilities. A total of 10¹¹ models are compared to the observed SED of each N2GN galaxy. SED3FIT is here used in two cases:

• (a)

  Sources with evidence of AGN activity, such as radio excess with respect to the radio-FIR correlation, X-ray emission, or MIR excess or power-law SED. In this case, we adopted the AGN-torus library by Fritz et al. (2006), updated by Feltre et al. (2012), that includes both the emission of the dusty torus, heated by the central AGN engine, and the emission of the accretion disc.

• (b)

  Sources with the short-wavelength photometry not reproduced by the standard MAGPHYS or CIGALE fit, thus having a rest-frame UV excess. In this case, we used SED3FIT with a library of non-extinguished young simple stellar populations (SSPs), of solar metallicity, Chabrier (2003) IMF, and with ages between 10 and 100 Myr, drawn from the stellar library by Bruzual & Charlot (2003). This choice is aimed at schematically simulating a young stellar component producing the observed excess and does not affect the estimate of dust masses.

Out of the whole N2GN sample, eight sources require an AGN-torus component to fit their MIR and UV emission. Six of them are also detected in the X-rays (Evans et al. 2024; Barro et al. 2019; Alexander et al. 2003). More objects are matched to X-rays sources, for a total number of 17, and in six cases the observed radio emission is in excess to the radio-FIR correlation of star forming galaxies. A detailed description of each individual source is given in Appendix E.

5. Derived properties of the N2CLS sources

The quantities produced by the SED fitting, of interest for this work, are the dust mass, M_dust, the infrared luminosity, L_IR, and the stellar mass, M_⋆, of the N2GN galaxies. The best fit of each individual galaxy is shown in Fig. E.1.

The best models of all N2GN sources (obtained with the UV-to-radio SED fitting) are also collected in Fig. 5, normalized by their L_IR (left) and by their M_dust (right). The normalization by M_dust is very similar to a sub-millimeter normalization along the Raileigh-Jeans (RJ) tail of the dust emission (e.g., at 850 μm in the rest frame). This is easily explained by the bulk of the dust mass budget being locked into the cold dust that dominates the RJ tail. The scatter is driven by the fact that the mass-to-light ratio is not unique, since each galaxy comes with its own value of M_dust/L_{sub − mm}. Consequently, because of this relation between sub-millimeter emission and M_dust, dust masses derived by SED fitting without the knowledge of the actual rest-frame sub-millimeter emitted luminosity are subject to large uncertainties (e.g., Berta et al. 2016) and/or need important assumptions in the SED modeling (e.g., fixing β in a MBB analysis).

The two middle panels of Figure 6 present the integrated infrared luminosity (L_IR, SF, 8−1000 μm, excluding the possible AGN-torus contribution), and the monochromatic 850 μm luminosity of the N2GN galaxies (L₈₅₀), as a function of redshift. The effects of the strong negative k-correction induced by the steep RJ tail of the dust SED are evident in both panels, and are more accentuated for L₈₅₀: on average the sub-millimeter luminosity of the sources keeps almost constant (and even slightly decreases) as redshift increases. The blue line in the L_IR, SF panel shows the expected luminosity limit, given the median SED template based on the N2GN sources and the 0.7 mJy flux cut at 1.2 mm. The infrared luminosity is affected by a much larger scatter than L₈₅₀, related to the variety of observed SED colors (Fig. 5).

The bottom panel of Fig. 6 reports the distribution of M_dust as a function of redshift. Thanks to the negative sub-millimeter k-correction, there is no evident dependence of M_dust on redshift. The dust masses of the N2GN galaxies span over the range between ∼10⁸ and ∼6 × 10⁹ M_⊙.

Fig. 5. Spectral energy distributions of all N2GN sources obtained with the UV-to-radio SED fitting (gray lines). The red and blue lines represent the median and average SEDs, respectively, as obtained by combining all models.

Fig. 6. Distribution of the N2GN galaxies, their luminosity, and dust mass as a function of redshift. Top panel: Redshift distribution. The filled histogram includes the sources with spectroscopic redshift, while the open histogram includes all sources. Second panel: Infrared luminosity of the star forming component, integrated between 8 and 1000 μm. The blue solid line represents the expected luminosity for a source of 0.7 mJy in the NIKA2 1.2 mm band, assuming its emission is described by the median SED shown in the left panel of Fig. 5. The shaded area is given by the best fit models of all sources. Third: 850 μm monochromatic luminosity. Bottom panel: Dust mass distribution as a function of redshift. The dashed and dotted lines represent the minimum dust mass and maximal-completeness dust mass obtained rescaling the L850 versus z line adopting the minimum and maximum Mdust/L850 ratio in the sample (Sect. 6.1).

5.1. Sources without redshift

Three N2GN sources have no match to multi-wavelength catalogs, except to radio data, namely N2GN_1_17b, 34a and 55. Therefore, these sources have no M_dust, nor a redshift determination. N2GN_1_17b and 34a are sub-components of NIKA2 sources and are detected by NOEMA at 2.0 mm and 1.2 mm, respectively. N2GN_1_55 is a single source, detected by NIKA2 at 1.2 mm. We derived the dust mass of these objects by analyzing the dependence of M_dust on millimeter flux of all other N2GN galaxies.

Figure 7 shows the derived M_dust versus the observed 1.2 mm flux of the galaxies. The red dashed line represents a simple least square fit to the sample and the dotted lines are the same fit, rescaled ±1.96× the residuals rms (corresponding to the 2.5−97.5 percentile range) around the fit. The red, solid vertical lines mark the 1.2 mm flux of the two sources without a redshift determination detected at 1.2 mm. A similar analysis was carried out at 2.0 mm for N2GN_1_17b. As a result, the range of dust masses associated to these three sources is roughly as broad as one order of magnitude. Table E.1 includes these values.

Fig. 7. Distribution of the N2GN galaxies as a function of dust mass and 1.2 mm flux density, color coded on the basis of their redshift. Median error bars are shown in the bottom-right corner. The red dashed line is a least squares fit to the data. The red dotted lines represent the same line scaled to the 2.5th and 97.5th percentiles of the residuals distribution. The red, solid vertical lines mark the flux of the sources with no redshift measurement and multi-wavelength counterparts available.

5.2. Dust mass to light ratio

For many dusty infrared galaxies, the absence of a fully sampled SED hinders the possibility to perform a detailed SED fitting and to derive a reliable estimate of – among other quantities – its dust mass. Thanks to the exquisite SED sampling of the N2GN sources, here we study the dust mass to rest-frame sub-millimeter luminosity ratio, with the goal to define a relation between these two quantities, to be used in cases that do not benefit from such a good wavelength coverage as the GOODS-N field. To this aim, Fig. 8 shows the M_dust/L_850 μm ratio, based on the κ_ν normalization by Draine et al. (2014) as usual.

Fig. 8. Study of the Mdust/L850 ratio. Top-left panel: least squares fit to the correlation between Mdust and L850. Right-hand panels: Mdust and Mdust/L850 as a function of IR luminosity.

We quantified the correlation between M_dust and L₈₅₀ with a simple least squares fit to the data, leading to the relation log(M_dust) = (1.24 ± 0.08)log(L₈₅₀)+(12.11 ± 0.20). The Spearman rank correlation coefficient is r_s = 0.83, with 66 degrees of freedom and a Student’s t distribution of t = 12.1, translating in a probability of correlation larger than 99.9%.

No correlation is seen between M_dust and L_IR (top-right panel), thus testifying that no artificial dependence between the two quantities has been introduced by the SED fitting underlying assumptions. On the other hand, there seem to be an apparent anti-correlation between M_dust/L₈₅₀ and L_IR (bottom-right panel), but it is driven by the avoidance zone dictated by the N2GN selection and by the more obvious correlation between L₈₅₀ and L_IR.

5.3. Starbursty nature of the N2GN galaxies

Adopting the scaling relation by Tacconi et al. (2020), that links M_⋆, redshift and distance from the main sequence (MS) of star forming galaxies to molecular gas mass (M_mol) and molecular gas depletion timescale (τ_dep), we computed τ_dep of the N2GN galaxies. Figure 9 compares the result to a collection of galaxies from the literature with M_mol derived from CO data (Tacconi et al. 2020; Berta et al. 2023, and references therein).

Fig. 9. Depletion timescale of the N2GN galaxies obtained applying the Tacconi et al. (2020) scaling relation. The gray shaded area is the trend found by Saintonge et al. (2013). The different lines represent the trends found by Tacconi et al. (2020) for: main sequence galaxies (selected within ΔMS = ±0.6 dex, black solid line); starburst galaxies (ΔMS > 0.6 dex, blue dashed line); extreme starbursts (ΔMS > 1.2 dex, purple dotted line); below-MS galaxies (ΔMS < 0.4 dex, red long-dashed line). Literature data are from Alaghband-Zadeh et al. (2013), Aravena et al. (2016, 2014, 2013), Bakx et al. (2020), Berta et al. (2023), Bothwell et al. (2017, 2013), Carilli et al. (2010), Chung et al. (2009), Combes et al. (2011, Combes et al. (2013, Dannerbauer et al. (2019), Decarli et al. (2016, 2019), Dunne et al. (2021, 2020), Fujimoto et al. (2017), Freundlich et al. (2019), Geach et al. (2011), Genzel et al. (2015, 2003), George et al. (2013), Hagimoto et al. (2023), Harris et al. (2012, 2010), Ivison et al. (2013, 2011, 2010), Penney et al. (2020), Riechers et al. (2020, 2011), Rudnick et al. (2017), Sharon et al. (2016), Solomon et al. (1997), Tacconi et al. (2018, 2013), Thomson et al. (2012), Valentino et al. (2018), Villanueva et al. (2017), Wang et al. (2018), and Yang et al. (2017).

The N2GN galaxies occupy mainly the locus of starbursts, with depletion timescales of the order of 0.1 to 1.0 Gyr. It is worth recalling that the Tacconi et al. (2020) relation uses the MS parametrization by Speagle et al. (2014) as reference, that does not include the bending of the MS at high stellar masses.

Figure 10 shows the position of the N2GN sources in the M_⋆ versus SFR space, as obtained with our fiducial UV-to-radio SED fitting (MAGPHYS or SED3FIT high-z) and for comparison with CIGALE. Two different parametrizations of the MS are overlaid to the data: those by Speagle et al. (2014, red lines) and Popesso et al. (2023, blue lines). The latter includes the well-known flattening at large stellar masses, while the former does not. Depending on the adopted reference MS, most of the N2GN sources are outliers of the MS (i.e., starbursts, case of Popesso et al. 2023) or include also a non-negligible number of MS galaxies (case of Speagle et al. 2014).

Fig. 10. Position of the N2CLS GOODS-N galaxies in the M⋆ versus SFR plane, as obtained with the SED fitting results, in different redshift bins. Each galaxy is positioned on the basis of its Δ(log(SFR))MS computed at its actual redshift with respect to the Popesso et al. (2023) MS, translated to the average redshift of each bin. The solid blue lines represent the parametrization of the MS of star forming galaxies by Popesso et al. (2023). For comparison the solid red lines refer to the one by Speagle et al. (2014). The dotted lines are placed at ±0.5 dex from the solid ones. Top row: results obtained with the fiducial fit (MAGPHYS or SED3FIT high-z). Bottom row: results obtained with CIGALE for comparison.

Noteworthy, the stellar mass of several N2GN high-redshift galaxies is of the order of a few 10¹¹ M_⊙. Similar very large values of M_⋆ at high redshift have been found in the recent literature, when exploiting deep JWST data (e.g., Xiao et al. 2023). They represent a challenge to galaxy formation models, as such high stellar masses would require very efficient star formation in the early stages of galaxy evolution, and would account for about 17% of the cosmic SFRD at z = 5 − 6 (Xiao et al. 2023).

For some of them, significant differences between MAGPHYS/SED3FIT and CIGALE are evident in Fig. 10. Appendix A discusses the direct comparison of the derived stellar masses, revealing a good overall agreement but also the presence of a few important outliers, with differences in the best fit M_⋆ estimate as large as one order of magnitude. The most critical cases are very red galaxies at z > 4. Appendix B shows the consequence on τ_dep of adopting the M_⋆ value derived with CIGALE instead of the one by MAGPHYS/SED3FIT. The main result of short depletion timescales is not affected. The most likely cause of these M_⋆ differences is a large difference of A_V for very extinguished galaxies.

The N2GN data set can arguably be considered the best photometric collection of an extragalactic blank field to date, including observations obtained with all major space-borne facilities available and covering the whole electromagnetic spectrum. Despite this opulence, it turns out that the majority of these M_⋆ > 10¹¹ M_⊙ high-z galaxies are affected by a large, uncovered gap between their rest-frame optical and FIR data. One example is shown in Fig. 11, where the 2.5%, 16% 84% and 97.5% percentiles of the modeled photometry are shown. Very deep MIR observations (in the observed frame, for example with JWST/MIRI) are necessary to constrain the NIR emission of these sources and their stellar mass.

Fig. 11. Example of fit uncertainty when no NIR and MIR data are available and a big gap in wavelength exists between the rest-frame optical and FIR spectral domains. The two different tones of red shaded areas represent the 2.5%, 16% 84% and 97.5% percentiles of the models, as computed in the photometric bands of the input catalog; the corresponding M⋆ range is as large as 0.25 dex (2.5th to 97.5th percentiles).

5.4. Blue excess

A few N2GN sources are characterized by a blue excess detected in the short-wavelength JWST and HST bands, with respect to the extinguished stellar models fitting the optical-NIR SED and providing the power to heat the bright FIR dust emission. At the redshift of the N2GN galaxies, this excess lies in the rest-frame UV or blue optical domain. For these dusty powerful star forming galaxies – hence excluding the case of passive galaxies with UV emission from evolved stellar populations – two possible components might contribute to this excess: a type-1 AGN component, or a young stellar population with no (or low) extinction.

As mentioned before, in absence of any other evidence of AGN activity, such as a MIR excess, optical/blue point-like morphology, (hard) X-ray emission or a radio excess with respect to the radio-FIR correlation, we opted to reproduce the observed blue data with a young SSP. Six N2GN sources were treated in this way, namely N2GN_1_06 and 13 (at z > 5) and 15, 49, 53, and 56b (at “cosmic noon”; Fig. E.1).

A similar excess of emission at short wavelengths was found for “little red dots” observed with JWST: compact red galaxies at z > 5, with F277W–F444W > 1, and F150W–F200W < 0.5 mag and F444W ≤ 28 mag (Labbe et al. 2025; Pérez-González et al. 2024; Matthee et al. 2024; Kocevski et al. 2024). Based on the rest-frame optical-NIR JWST photometry, the very red color testifies the presence of a large amount of dust, while the blue excess can be explained by either a QSO, a clumpy ISM surrounding the star forming regions, allowing us to see unobscured very young star formation, or a gray attenuation law, typically linked to significant scattering (Pérez-González et al. 2024). The two z > 5 such cases present in N2CLS GOODS-N (N2GN_1_06 and 13) benefit from a very complete multi-wavelength coverage and do not show any hints of AGN activity over the whole X-rays to radio broad-band SED.

6. Comoving number density

The comoving number density of the sources in the N2GN survey, in intervals of M_dust, also known as DMF, is given by

$$\begin{array}{c}\hfill \mathrm{\Phi }(M_{\mathrm{dust}})\mathrm{\Delta }{M}_{\mathrm{dust}}=\sum _i\frac{1}{V_a^i}\mathrm{\Delta }{M}_{\mathrm{dust}},\end{array}$$(4)

where the sum is computed over all galaxies in the given mass bin of width ΔM_dust. The volume within which each source is accessible to the N2CLS GOODS-N survey is a spherical shell:

$$\begin{array}{c}\hfill V_a=\frac{\mathrm{\Omega }}{4\pi }{\displaystyle {\int }_{z_{\min }^{\mathrm{bin}}}^{\min (z_{\max },z_{\max }^{\mathrm{bin}})}}\frac{\mathrm{d}V}{\mathrm{d}z},\end{array}$$(5)

where the minimum redshift is given by the lower boundary of the given redshift bin considered and the maximum redshift is the minimum value between the upper boundary of the bin and the maximum accessible redshift. The latter is the highest redshift at which a galaxy would be observable in the survey (Schmidt 1968), given the N2GN 1.2 mm flux limit of 0.7 mJy.

The total area of the N2GN field is 159 arcmin² (Bing et al. 2023). Nevertheless, the depth of the NIKA2 map is not uniform, but it varies by up to a factor of three across the N2GN field. As a consequence, Ω is the effective area associated to each galaxy, as computed by Bing et al. (2023).

Because of the strong negative k-correction along the steep RJ tail of the dust SED, z_max turns out to be very large (typically z_max > 20). Hence, the accessible volume V_a is limited by the upper boundary of the given redshift bin.

6.1. Dust mass completeness

The N2GN sample selection is based on a 1.2 mm flux cut in the observed frame. Flux completeness has already been taken into account by using the effective area of each source in the computation of the accessible volume.

The dust mass budget is dominated by the rest-frame sub-millimeter emission. Therefore, there is an almost-direct link between the flux selection and the derived quantity M_dust. Nevertheless, the relation between S(1.2 mm) and M_dust is not univocal, because the latter is derived from the former by means of a proper SED fitting of each galaxy, and therefore the dust mass to sub-millimeter light ratio varies in the sample (Fig. 8 and associated text). Consequently, it is not possible to define a sharp mass limit encompassing the whole sample.

The bottom panel of Fig. 6 shows the distribution of M_dust as a function of redshift. The dashed and dotted lines represent the minimum dust mass and maximal-completeness dust mass obtained rescaling the L₈₅₀ versus z line adopting the minimum and maximum M_dust/L₈₅₀ ratio in the sample. The shaded areas are based on the best fit SED models of all N2GN sources.

The minimal area represents the absolute minimum M_dust that objects in the sample could have, at the given redshift, if they had the minimum “observed” M_dust/L₈₅₀ and a flux of 0.7 mJy at 1.2 mm, accounting for the scatter due to the different SED shapes of the sample. The maximal area, instead, represents the M_dust that an object would have, if its flux were 0.7 mJy at 1.2 mm, and if it were characterized by the maximum “observed” M_dust/L₈₅₀ ratio in the N2GN sample, given the scatter of SED shapes. These maximal M/L tracks do not indicate a maximum M_dust limit. Schematically, above this area/line the sample is maximally complete in terms of dust mass. Below it, dust mass incompleteness must be corrected.

This effect has been first examined in the past for the stellar mass function (e.g., Dickinson et al. 2003b; Fontana et al. 2004, 2006; Berta et al. 2007) and later on for the gas mass function and the DMF (e.g., Berta et al. 2013b; Pozzi et al. 2020). Here we applied the recipe by Fontana et al. (2004) to the case of dust masses and the 1.2 mm flux cut.

Figure 12 shows the distribution of dust masses as a function of the observed 1.2 mm flux for the N2GN galaxies in the 1.6 < z ≤ 2.4 redshift bin, as an example. A similar procedure was applied also to the other redshift bins. The vertical line represents the adopted 1.2 mm flux limit. The diagonal shaded areas are the tracks described by the minimum and maximum M_dust/L₈₅₀ ratios in the sample, at the central redshift, taking into account the scatter due to the different SED shapes of the N2GN sources. The horizontal dashed line represents the dust mass level above which the sample is definitely complete, for a given SED.

Fig. 12. Exemplification of how the dust mass completeness is computed in the redshift bin 1.6 < z ≤ 2.4 (Sect. 6.1). The diagonal shaded areas represent tracks of the N2GN SED models, given the minimum and maximum Mdust/L850 ratios in the sample. The vertical dashed line is the 1.2 mm flux limit of the survey, and the horizontal dashed line is the corresponding mass maximal completeness threshold. The dotted horizontal lines mark the boundaries of an example mass bin, in which the hatched area would contain the sources missed by the flux cut.

In a given mass bin, $M_{\mathrm{dust}}^{\mathrm{lo}}<M_{\mathrm{dust}}<M_{\mathrm{dust}}^{\mathrm{up}}$, below the maximal completeness mass, the observed 1.2 mm flux corresponding to a given mass is encompassed between the maximal and minimum diagonal areas (S^lo, S^up), but the flux limit of the survey cuts the distribution. The sources in the shaded area are missed by the survey. The fraction of source actually observed, with respect to the total in that bin (i.e., the mass completeness) is

$$\begin{array}{c}\hfill f_{\mathrm{compl}}=\frac{{\displaystyle {\int }_{S_{\lim }}^{S_{\mathrm{up}}}}N(M,S^{\prime })\mathrm{d}{S}^{\prime }}{{\displaystyle {\int }_{S_{\mathrm{lo}}}^{S_{\mathrm{up}}}}N(M,S^{\prime })\mathrm{d}{S}^{\prime }}·\end{array}$$(6)

Since the actual distribution of the galaxies below the flux limit is not known, we assumed that it does not depend of M_dust, implying that the distribution observed above the maximal mass threshold still holds below.

We computed the observed DMF only in those bins with dust mass completeness ≥0.5. Table 2 lists the completeness values for each mass bin considered. We defer the reader to Fontana et al. (2004). Berta et al. (2007). Pozzi et al. (2020) for further details.

6.2. Dust mass function

The comoving number density of galaxies as a function of dust mass (Eq. 4) was corrected for mass incompleteness as described above. The incompleteness due to the non homogeneous map coverage and the limited S/N of the survey was already taken into account in the computation of V_a, when accounting for the effective area associated to each source.

Figure 13 and Tab. 2 present the resulting DMF, in three redshift bins: 1.6 < z ≤ 2.4, 2.4 < z ≤ 4.2, and 4.2 < z ≤ 7.2. The choice of these redshift bins comes from an optimization of the DMF signal and maximizes the number of sources in each bin. Only a few mass bins are populated because of the overall relatively small number of sources. Error bars account for M_dust uncertainties, based on the probability distribution function (PDF) produced by SED fitting, redshift uncertainties, and Poisson errors.

Fig. 13. Dust mass function of the N2GN sources (red symbols). Only the mass bins containing at least two galaxies and with a completeness ≥50% are plotted. Uncertainties take into account ΔMdust, including SED fitting parameters sampling and Δz contributions, and Poisson errors. Left panel: comparison to literature results by Pozzi et al. (2020, blue filled squares and long-dashed line) and Traina et al. (2024a, light blue triangles and short-dashed lines, fiducial results in the 2.0 < z ≤ 2.5, 3.5 < z ≤ 4.5 and 4.5 < z ≤ 6.0 redshift bins). The blue dotted lines are the same Pozzi et al. (2020) DMF, evolved to z = 2.0 and 3.3 (see Sect. 6.2). Right panel: results of the STY analysis. The green solid lines and shaded areas are the result of the STY parametric analysis: most probable model (solid lines) and 1σ uncertainty (shaded areas). The dashed green lines represent the 1.6 < z ≤ 2.4 result, repeated in the other two redshift bins for comparison.

Pozzi et al. (2020) derived the DMF of 5546 galaxies detected by Herschel in the COSMOS field, in the redshift range 0.1 < z < 2.5. They adopted a dust κ_ν(250 μm) = 4.0 cm² g⁻¹ (Bianchi 2013), which corresponds to the Draine (2003) dust models normalization. Therefore, we rescaled their DMF to our κ_ν of reference (Draine et al. 2014). The result is compared to the N2GN DMF at 1.8 < z ≤ 2.5 in the left-hand panel of Fig. 13. The two estimates are consistent within the errors, the Herschel-based DMF being systematically lower than the NIKA2 one.

Traina et al. (2024a) computed the DMF of 189 galaxies detected at millimeter wavelengths by ALMA in the COSMOS field as part of the A³ COSMOS data collection (Liu et al. 2019; Adscheid et al. 2024; Traina et al. 2024b). Their analysis covers the range from z = 0.5 to z = 6.0, split in eight redshift bins. In the left panel of Fig. 13, the data by Traina et al. (2024a) are shown for only three z bins, intersecting those of N2GN. The A³ COSMOS DMF of the other intersecting bins is very similar and we omit it for sake of figure readability. The DMF determined by Traina et al. (2024a) is very similar to our results, and extends to larger dust masses.

The DMF can be reproduced with a Schechter (1976) function. For logarithmic mass bins, this is parametrized as

$$\begin{array}{c}\hfill \mathrm{\Phi }(M_{\mathrm{dust}})\mathrm{d}log{M}_{\mathrm{dust}}=\frac{{\mathrm{\Phi }}_{\mathrm{dust}}^{\ast }}{M_{\mathrm{dust}}^{\ast }}{\left(\frac{M_{\mathrm{dust}}}{M_{\mathrm{dust}}^{\ast }}\right)}^{\alpha +1}{e}^{-\frac{M_{\mathrm{dust}}}{M_{\mathrm{dust}}^{\ast }}}ln(10)\mathrm{d}log{M}_{\mathrm{dust}}.\end{array}$$(7)

Pozzi et al. (2020) fitted the observed DMF with a non-linear least squares approach to determine the Schechter parameters over six redshift bins. At z < 0.25 they sampled the DMF over a mass range sufficiently large to constrain the three parameters α, $M_{\mathrm{dust}}^{*}$ and ${\mathrm{\Phi }}_{\mathrm{dust}}^{*}$. At higher redshift, given the bright flux cut of their survey, they sampled only the massive tail of the DMF, and they needed to fix α to the local value of 1.48, in order to derive the evolution of $M_{\mathrm{dust}}^{*}$ and ${\mathrm{\Phi }}_{\mathrm{dust}}^{*}$.

We studied the variation of $M_{\mathrm{dust}}^{*}$ as a function of cosmic age, based on Pozzi et al. (2020) results. By performing a simple least square fit to their $M_{\mathrm{dust}}^{*}$ values as a function of cosmic time, we found the relation $log{M}_{\mathrm{dust}}^{*}=(-2.2\pm 0.1)\times log{t}_{\mathrm{Univ}}+(30.1\pm 1.0)$. The units are years for the age of the Universe and M_⊙ for $M_{\mathrm{dust}}^{*}$. While deriving this equation, we rescaled their values of $M_{\mathrm{dust}}^{*}$ to the κ_ν normalization by Draine et al. (2014).

Using the derived equation, we evolved the Pozzi et al. (2020) $M_{\mathrm{dust}}^{*}$ to the central redshift of our two lower redshift bins (z = 2.0, 3.3), but we did not extrapolate it to the third (z = 5.7), because it is too far from the range of redshift actually covered by Pozzi et al. (2020). The light-blue long-dashed line in the top panel of Fig. 13 is the original 1.8 < z < 2.5 DMF by Pozzi et al. (2020). The blue dotted lines in each panel represent the Pozzi et al. (2020) DMF evolved as described above, rescaling their value of ${\mathrm{\Phi }}_{\mathrm{dust}}^{*}$ to match the observed N2GN DMF.

In both 1.6 < z ≤ 2.5 and 2.5 < z ≤ 4.2 redshift bins, the result (dotted lines in Fig. 13) is consistent with the N2GN observed DMF. Notably so is also the low-mass end slope α = 1.48 derived locally by Pozzi et al. (2020), suggesting that the $M_{\mathrm{dust}}^{*}$ derived from their data is a good approximation at least up to z ∼ 4.2.

Traina et al. (2024a) performed a Markov Chain Monte Carlo (MCMC) fitting of the observed DMF. Their best fiducial model was obtained by fitting a combination of the A³ COSMOS data and the Herschel results by Pozzi et al. (2020), simultaneously over all redshift bins. While doing so, they fixed the low-mass slope to the value determined locally by Pozzi et al. (2020). Their result is overplotted to the N2GN observed DMF in Fig. 13. Their observed DMF is well consistent with the NIKA2 DMF at all redshifts, over the mass ranges in common, but at z > 2.4 their Schechter parametrization differs significantly from ours (Sect. 6.3).

6.3. STY analysis

We applied the STY (Sandage, Tammann, & Yahil) method (Sandage et al. 1979) to the N2GN galaxies, with the aim to derive new Schechter parameters describing the evolution of the DMF, based on the N2GN data themselves. It is important to remind that this is not a fit to the 1/V_a mass density using a Schechter function, but a statistical method taking into account the dust properties of each individual galaxy that contributes to the mass function. In this way, the method naturally includes all sources and accounts for the Poisson noise and the M_dust uncertainty (including the errors on z).

The objects without a known redshift are included in this piece of analysis such that they contribute to the uncertainty of the Schechter parameters and of the integrated dust mass density (Sect. 7.1) in each redshift bin considered. Nevertheless, they have a very small impact (≤2 − 3%) on the mass function and on the value of ρ.

Following Berta et al. (2007); Berta et al. (2013b), we used a MCMC sampling of the parameter space to explore the posterior probability function of the Schechter model, with parameters comprising both the three Schechter $M_{\mathrm{dust}}^{*}$, ${\mathrm{\Phi }}_{\mathrm{dust}}^{*}$ and α, and additional hyper-parameters to represent the dust mass of each individual galaxy dust masses. The MCMC engine was used with a Metropolis-Hastings sampling algorithm to explore the posterior probability function of the model. The hyper-parameters are constrained solely by the prior knowledge of the dust mass probability distribution function of each galaxy (produced by the SED fitting). Therefore, their marginalization automatically accounts for the mass uncertainties. We defer to Berta et al. (2007) for a thorough description of the method and of the adopted algorithm.

Figure 13 shows the results of the STY analysis and Tab. 3 summarizes the values found for the Schechter parameters. Despite the range of dust masses covered by the data is limited, we are able to set a constraint on all three Schechter parameters and their uncertainty.

7. Discussion

The variation of M^* and Φ^*, obtained with the STY study of the DMF (Sect. 6.2, Tab. 3), reveals a rapid evolution of the comoving density of DSFGs, roughly doubling in the first Gyr from z = 4.2 − 7.2 to z = 2.4 − 4.2, and increasing by more than a factor of 5 during the next two Gyr. Based on previous works, at more recent epochs, the comoving density evolution slows down and practically ceases by z ∼ 0.5 − 1.0 (Pozzi et al. 2020).

The characteristic dust mass, $M_{\mathrm{dust}}^{*}$, does not evolve at least until redshift z ∼ 4 and then is halved by z ∼ 2. The evolution continues down to the local Universe, decreasing $M_{\mathrm{dust}}^{*}$ by more than an order of magnitude until z ∼ 0.

These trends indicate a scenario in which DSFGs at “cosmic noon” are more numerous than at earlier epochs and contain more dust than their low-z cousins. This two-fold evolution can be explained with dust being rapidly produced in the distant Universe, heating up and emitting at FIR-millimeter wavelengths, thus peaking at z = 2 − 3, and finally being slowly and gradually consumed by star formation or expelled by galactic winds in more recent epochs, when dust production by AGB stars and Supernovae does not balance the “losses” anymore.

7.1. Evolution of the dust mass cosmic density

The integral of the Schechter function is given by the upper incomplete Gamma function:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\rho }_{\mathrm{dust}}(M_{\mathrm{dust}}>M_{\inf })& ={\displaystyle {\int }_{M_{\inf }}^{\infty }}\mathrm{\Phi }\left(M_{\mathrm{dust}}\right)\mathrm{d}log{M}_{\mathrm{dust}}\hfill \\ \hfill & =M_{\mathrm{dust}}^{\ast }{\mathrm{\Phi }}_{\mathrm{dust}}^{\ast }\mathrm{\Gamma }(\alpha +2,\frac{M_{\mathrm{dust}}}{M_{\mathrm{dust}}^{\ast }})·\hfill \end{array}\end{array}$$(8)

In its calculation, we adopted M_inf = 10⁴ M_⊙. The resulting dust mass density of the N2GN survey is shown in Fig. 14 and compared to data found in the literature.

Fig. 14. Evolution of the dust mass density as a function of look-back time. The N2GN results are depicted as red filled circles. All other symbols belong to literature data, as indicated in the right-hand panel. Local estimates are by Vlahakis et al. (2005) and Beeston et al. (2018); combinations of optical-NIR and FIR-sub-millimeter data are by Eales & Ward (2024) and Driver et al. (2018); Herschel-based studies are by Beeston et al. (2024), Pozzi et al. (2020) and Dunne et al. (2011); sub-millimeter and millimeter driven results are by Pozzi et al. (2021), Dudzevičiūtė et al. (2021), and Magnelli et al. (2020); and finally derivations of ρdust are based on Ωgas (Péroux & Howk 2020) and MgII absorbers (Ménard & Fukugita 2012). All data have been rescaled to our chosen cosmology and κν (Draine et al. 2014). Figure C.1 presents a full-color version and Fig. D.1 shows ρdust as a function of redshift.

All literature data have been rescaled to the Draine et al. (2014) value of κ₈₅₀ = 0.047 m² kg⁻¹. Appendix C lists the original assumptions of the different authors and shows the importance of rescaling all results to the same κ_ν in order to properly compare them to each other.

The local estimates of ρ_dust shown in Fig. 14 were established by Vlahakis et al. (2005) and Beeston et al. (2018). The former observed a sample of 81 optically selected galaxies with SCUBA at 450 and 850 μm. The latter studied the local DMF of 15 750 galaxies drawn from the 145 deg² intersection of the Galaxy And Mass Assembly (GAMA; Driver et al. 2009; Liske et al. 2015) and the Herschel Astrophysical Terahertz Large Area Survey (H-ATLAS; Eales et al. 2010) surveys, detected in the optical r band and in the Herschel 250 μm SPIRE band.

In this collection, FIR Herschel-based studies are included: the z < 0.5 Dunne et al. (2011) result, comprising 1867 H-ATLAS 250 μm sources with a reliable SDSS counterpart; the z < 0.5 Beeston et al. (2024) data, based on a sample of 29 241 galaxies in H-ATLAS; and Pozzi et al. (2020), already mentioned earlier, including 5546 PEP 160 μm selected galaxies with multi-wavelength counterparts, in the redshift range z = 0.2 − 2.5.

Combining optically selected galaxies and long wavelength data, Driver et al. (2018) brought together three complementary data sets, GAMA (Driver et al. 2011; Liske et al. 2015), G10-COSMOS (Davies et al. 2015; Andrews et al. 2017), and 3D-HST (Momcheva et al. 2016) for a total of more than 500 000 optical and NIR selected galaxies up to z ≤ 1.75. Eales & Ward (2024) combined the stellar mass function by Davidzon et al. (2017) in the COSMOS field (Scoville et al. 2007) with the SCUBA-2 850 μm stacking analysis by Millard et al. (2020) in the same field to study the evolution of the cosmic density of dust from z = 0.1 to z = 5.5.

Few sub-millimeter-driven studies are included, based on SCUBA-2 and ALMA data. Dudzevičiūtė et al. (2021) derived ρ_dust by fitting the UV-to-radio SEDs of 450 μm galaxies selected from the SCUBA-2 STUDIES survey (Wang et al. 2017; Lim et al. 2020). They compared and combined the results to a similar modeling carried out on 850 μm ALMA/SCUBA-2 AS2UDS (Stach et al. 2019; Dudzevičiūtė et al. 2020) galaxies. Magnelli et al. (2020) used the ALMA spectroscopic survey (ASPECS) in the Hubble Ultra Deep Field (HUDF) to study the contribution of H band-selected galaxies to the dust cosmic density. These authors performed stacking on the 1.2 mm ALMA maps of the selected galaxies in bins of stellar mass, over the redshift range z = 0.3 − 3.2. The resulting ρ_dust is limited by the lower limit of the stellar mass selection and can be considered a lower limit. Also applying a stacking technique, Pozzi et al. (2021) derived the contribution to ρ_dust of 113 UV-selected galaxies at redshift z = 4.4 − 5.9 in the ALMA/ALPINE survey (Le Fèvre et al. 2020). They also estimated the dust density given by additional 23 serendipitously detected sub-millimeter galaxies. Traina et al. (2024a) measured the evolutionary ρ_dust from z = 6.0 to z = 0.5 analyzing the ALMA-selected A³ COSMOS sample mentioned in Sect. 6.2. In Fig. 14, we show their fiducial results, based on CIGALE SED fitting and the Schechter fit to ALMA plus Herschel data. Finally, Ménard & Fukugita (2012) computed the dust density studying MgII absorbing clouds in the lines of sight of quasars, and deriving the amount of intervening dust by fitting the corresponding extinction curve, and Péroux & Howk (2020) derived the dust density Ω_dust = ρ_dust/ρ_crit, 0 from the gas density Ω_gas, making an assumption on the dust to gas ratio.

The N2GN and A³ COSMOS samples are practically the only ones consisting of galaxies directly selected at observed millimeter wavelengths, that translate into the rest-frame sub-millimeter domain. As pointed out in Sect. 5, selecting along the RJ tail of the dust emission is paramount to have as unbiased a M_dust estimate as possible. Moreover these galaxies benefit from a very rich multi-wavelength coverage, that allowed us to avoid any strong assumption in the SED fitting that lead to the M_dust measurements (e.g., fixing the dust temperature). The two samples lead to different results: in this work we detect a strong evolution of ρ_dust from z ∼ 7 to z > 2, while Traina et al. (2024a) find a smoother increase of the dust density from z = 6 to z = 0.5. The A³ COSMOS values of ρ_dust are larger than those derived here for N2GN by a factor ∼5 in the largest redshift bin (i.e., ∼3.3σ higher) and by a factor 1.5 to 2.5 at 2.4 < z ≤ 4.2 (i.e., 2−3σ).

The large spread in the observed values of ρ_dust in Fig. 14, and the apparently discordant trends of ρ_dust at z < 1.5 in the different literature works shown here, are likely to be (at least in part) ascribed to differences in the selection of the samples and in the methods adopted to derive M_dust among different authors. More systematic and extensive measurements, possibly carried out self-consistently over the whole z = 0 − 7 redshift range, and over larger sky areas (e.g., COSMOS), are needed to be more conclusive about the actual cosmic evolution of dust.

A larger survey area will also give the chance to detect rare objects at even higher redshift than in N2GN. The analysis of z_max in Sect. 6.2 has in fact shown that – if dust already existed in such amounts – at the N2CLS depth a dusty star forming galaxy would be detected at z ≫ 7.

7.2. Comparison to models predictions

In Fig. 15, the N2GN data and the collection of literature ρ_dust results are compared to model predictions by Popping et al. (2017), Gioannini et al. (2017), Aoyama et al. (2018), Li et al. (2019), Lewis et al. (2023), Parente et al. (2023), and Yates et al. (2024).

Fig. 15. Comparison of the dust mass density to model predictions. The N2GN results are depicted as red filled circles. The different lines represent the semi-analytical models by Popping et al. (2017), Parente et al. (2023) and Yates et al. (2024), the chemical evolution model by Gioannini et al. (2017), and the cosmological hydrodynamical simulations by Lewis et al. (2023), Li et al. (2019) and Aoyama et al. (2018). Figure D.1 presents a version of this figure as a function of redshift.

Popping et al. (2017) included dust production and destruction in a semi-analytic model of galaxy formation, taking into account condensation of dust in stellar ejecta, growth of dust in the ISM, the destruction of dust by supernovae and in hot halos, and dusty winds and inflows. Their model reproduces the relation between stellar mass and dust mass in the local and high-redshift Universe, as well as the dust-to-gas ratio of local galaxies as a function of their stellar mass.

Gioannini et al. (2017) studied the evolution of interstellar dust with chemical evolution models of galaxies of different morphological types, including dust production from supernovae and asymptotic giant branch stars, dust accretion and destruction processes. They explored different cosmological scenarios of galaxy formation: pure luminosity evolution and number density evolution. The latter is plotted in Fig. 15.

Aoyama et al. (2018) investigated the evolution of dust by means of hydrodynamic simulations, in which the enrichment of metals and dust is treated self-consistently with star formation and stellar feedback. Dust evolution is driven by production in stellar ejecta, destruction by sputtering, grain growth by accretion and coagulation, and grain disruption by shattering. The model is successful at reproducing the local relation between dust-to-gas ratio and metallicity of galaxies. In Fig. 15 only the dust component locked in galaxies is plotted.

Li et al. (2019) implemented the production, growth, and destruction of dust grains into the SIMBA cosmological hydrodynamic galaxy formation simulation (Davé et al. 2019). Dust forms in stellar ejecta, grows by the accretion of metals, and is destroyed by thermal sputtering and supernovae. This simulation reproduces the observed local DMF but under-predicts the DMF by a factor ∼3 in the redshift range z = 1 to 2. In Fig. 15 we show the prediction obtained including and excluding the dust ejected out of the galaxies via galactic winds.

The Lewis et al. (2023) model, called DUSTiER (DUST in the Epoch of Reionization), couples a physical model for dust production to the fully radiation-hydrodynamics cosmological simulation code RAMSES-CUDATON (Ocvirk et al. 2016), with the aim to study how dust affects the escape of ionizing photons into the IGM during the Epoch of Reionization. In Fig. 15 the fiducial model is shown (at z ≥ 5); the uncertainty range is ∼0.3 dex (not shown).

Parente et al. (2023) implemented some standard prescriptions for dust production and evolution on top of the L-GALAXIES2020 semi-analytic model (Henriques et al. 2020). Their approach adopts a new disc instability criterion to trigger bulge and central black hole growth. As for the dust properties, the model includes two populations of large and small grains (e.g., Granato et al. 2021), production by stars, shattering and coagulation, grain growth in molecular clouds, destruction in SNe shocks, sputtering in the hot phase and ejection from the host galaxy. The model is successful at reproducing the local stellar mass function and DMF, but it is known to lack dusty highly star forming galaxies at z > 1 (Traina et al. 2024a). The Yates et al. (2024) simulation is a further implementation of the L-GALAXIES model, taking into account binary stars evolution and dust production and destruction. This implementation reproduces observed dust-to-metal and dust-to-gas ratios at z = 0 to 4 but systematically underpredicts the dust masses of galaxies at z > 4. In Fig. 15 the fiducial model of Parente et al. (2023) and the total ISM dust component of Yates et al. (2024) are shown.

All models predict a steep fall of the dust mass cosmic density above z = 2 − 3, despite significant differences. The models by Popping et al. (2017) and Li et al. (2019) are the most successful at predicting the z = 2.4 − 7.2 ρ_dust as determined by N2GN. The model by Gioannini et al. (2017) is consistent with higher measurements of ρ_dust as those by Eales & Ward (2024).

At lower redshift, the differences between the different simulations become critical. In the Popping et al. (2017) model, the dust density is roughly constant from z ∼ 2 to the local Universe, suggesting a quasi-equilibrium between the production and destruction of dust grains (Magnelli et al. 2020). This trend is consistent with the observed density by Driver et al. (2018), based on optically selected galaxies. On the other hand, a FIR selection leads to an evolving ρ_dust(z), steadily decreasing from z ∼ 1 to ∼0.2 (Pozzi et al. 2020; Eales & Ward 2024). The model by Gioannini et al. (2017) fairly reproduces this decrease.

The hydrodynamical simulations by Aoyama et al. (2018) and Li et al. (2019) significantly over-predict the low-redshift dust density. Notably, when excluding the fraction of dust ejected from the galaxies into the IGM, the Li et al. (2019) model predicts a decrease of ρ_dust similar to the one observed by Pozzi et al. (2020). Arguably, the expelled dust residing in the IGM is too faint to contribute to the observed (sub-)mm flux of the N2GN galaxies, and therefore it is not included in the measured DMF and ρ_dust.

Given the large scatter in the data, after re-scaling all estimates to the same κ_ν assumptions, it is difficult to be conclusive on which model is most successful in reproducing the observed evolution of ρ_dust. Once again, selection effects and differences in the assumptions governing the dust mass derivation produce too large of a scatter to allow either of the models to be favored.

8. Summary and conclusions

We have conducted a panchromatic study of the sources detected by NIKA2 in the GOODS-N field, combining multi-wavelength data obtained with NIKA2 at the IRAM 30 m telescope, NOEMA, SMA, SCUBA2, VLA, Herschel, Spitzer, the Hubble and James Webb space telescopes, and ground-based optical facilities. The sample is named N2GN.

Out of the 71 individual galaxies associated with the N2GN detections, 68 have a detailed SED coverage and a redshift measurement. Redshifts span the range 0.6 < z < 7.2, with a median of 2.819. We performed SED fitting of these sources with eight different methods, including a MBB in the optically thin approximation and in its general form, the Draine & Li (2007) models, and UV-to-radio fits with MAGPHYS (da Cunha et al. 2008, 2015; Battisti et al. 2020), CIGALE (Burgarella et al. 2005; Noll et al. 2009; Boquien et al. 2019), and SED3FIT (Berta et al. 2013a). The end product of this analysis is a robust estimate of their dust masses and of the evolution of the DMF and the dust cosmic density as a function of redshift. The main results of this analysis are the following.

• Once the differences in their underlying assumptions are taken into account (e.g., in the reference value of the dust absorption coefficient κ_ν), the eight methods produce consistent estimates of M_dust. The only exception is the MBB in general form, which tends to underestimate the dust mass. This model is known to fail when no information about the size of the dust emitting region, or the frequency at which the medium becomes optically thick, are known (Ismail et al. 2023). We also found a similar good consistency between the different methods for L_IR. As for M_⋆, the different codes produce results in general agreement, but significant differences exist for individual sources, especially those with very red optical SEDs.

• Combining the SFR derived from L_IR and the stellar mass, we compared the N2GN galaxies to the “main sequence” locus between redshift 1.6 and 7.2. When using the parametrization of the MS given by Popesso et al. (2023), the majority of the sources turn out to be outliers that are likely powered by an intense starburst.

• The stellar mass of several z > 4 N2GN galaxies is of the order of a few 10¹¹ M_⊙. If confirmed, such a high M_⋆ would challenge galaxy formation models and imply very high star formation efficiencies (Xiao et al. 2023). Nevertheless, most of these sources have a very broad gap in their observed SEDs, as they lack data between rest-frame optical and FIR wavelengths. Moreover significant differences exist between the M_⋆ estimates obtained by different codes, which is possibly driven by a different SFH or attenuation law details. Only future very deep MIR observations (e.g., with JWST/MIRI) will be able to constrain the rest-frame NIR emission and the stellar mass of these galaxies.

• A few sources were detected in the X-rays or showed other evident signs of AGN activity (power-law-like MIR SED, radio excess, optical point-like morphology). These sources required an AGN-torus component in the SED fitting (Berta et al. 2013a). In some other cases, the rest-frame UV-blue optical SED had a strong excess with respect to the extinguished stellar component. This excess can be explained with the contribution of a non-extinguished young stellar population to the SED. Two of these blue-excess galaxies lie at z > 5, and their SEDs are similar to those of the “little red dots” recently identified by JWST (Pérez-González et al. 2024).

• The millimeter selection, possible thanks to NIKA2, closely resembles a dust mass selection modulo the spread in redshift and in the dust mass to light ratio of the sample. The strong correlation between the dust mass and the rest-frame 850 μm luminosity of the N2GN galaxies is described by the equation: log(M_dust) = (1.24 ± 0.08)log(L₈₅₀)+(12.11 ± 0.20). This relation can be useful for galaxies lacking an SED coverage detailed enough to derive M_dust.

• We derived the DMF of the N2GN galaxies over the redshift range between z = 1.6 and z = 7.2, and we reproduced it with a Schechter function using the STY method. Combined with literature results (e.g., Pozzi et al. 2020; Eales & Ward 2024), our data reveal a two-fold evolution of the DMF of millimeter galaxies: (a) a rapid evolution of the characteristic comoving density Φ^*, that doubles from z ∼ 7 to z ∼ 2.5, increases by a further factor of five at later epochs and ceases at z ∼ 0.5 − 1.0; (b) a constant characteristic dust mass, $M_{\mathrm{dust}}^{*}$, from z ∼ 7 to z ∼ 4 and then a smooth decrease by more than an order of magnitude down to the local Universe.

• By integrating the DMF, we derived the comoving dust mass density of millimeter galaxies and compared it with previous studies. In this comparison, we showed that it is of uttermost importance to rescale all ρ_dust determinations to the same κ_ν assumption. The comoving ρ_dust increased by at least one order of magnitude from z ∼ 7 to z ∼ 1.5. At lower redshift, a different data set leads to different trends, with ρ_dust possibly remaining constant down to z ∼ 0 or decreasing by a factor of approximately three. The N2GN sample leads to results systematically different from the recent ALMA A³ COSMOS determination of ρ_dust (Traina et al. 2024a), with the latter showing a smoother evolution and a larger dust mass density at the highest redshifts. A comparison with simulations of dust and galaxy evolution favors semi-analytical and chemical evolution models that take into account dust growth and destruction (Popping et al. 2017; Gioannini et al. 2017). Hydrodynamical simulations are also successful when excluding the fraction expelled into the IGM from the dust budget (Li et al. 2019).

The deep NIKA2 maps of GOODS-N observed by the N2CLS survey have enabled for the first time the determination of the dust mass comoving density to be extended up to ∼7, that is, an epoch when the Universe was less than 1 Gyr old. A very large spread in the determinations of ρ_dust at all redshifts (z = 0 − 7) still exists. The discrepancies between different works likely arise from the differences in sample selections, in the methods employed to derive M_dust, and possibly in their underlying assumptions. A systematic and self-consistent study of M_dust and ρ_dust over as broad a redshift range as possible, from the local Universe to remote epochs, is envisaged and recommended in order to reach a consensus about the actual cosmic evolution of dust.

Future studies will need to be carried out on larger areas in order to cover larger volumes and be sensitive to both lower redshift regimes and rare very distant dusty galaxies. Thanks to the negative k-correction at millimeter wavelengths, at the unprecedented depth of N2CLS, a dusty star forming galaxy would be easily detected at z ≫ 7 if dust already existed in galaxies in amounts similar to those in the N2GN sample. The natural next step in the exploration of ρ_dust will exploit the N2CLS COSMOS field (N2CLS team, in prep.).

Data availability

The N2CLS final maps and catalogs, the NOEMA follow-up data, and the multi-wavelength catalog are available on line on the survey’s home page: https://data.lam.fr/n2cls/home. The JWST data are available from the Mikulski Archive for Space Telescopes (MAST: https://mast.stsci.edu). Figure E.1 is available on Zenodo, at the address https://doi.org/10.5281/zenodo.14965423

Acknowledgments

We thank the anonymous referee for their insightful comments that improved the quality and the presentation of this work. We recognize the contribution of Mael Voyer, Jean Anquetil and Aymeric Garnier to the early identification process of the N2CLS sources. We are grateful to J. Lewis for providing the DUSTIER model. This work is based on observations carried out under project numbers 192-16 with the IRAM 30 m telescope, and projects W16EE, E16AI, W21CV and W23CX with NOEMA. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). We would like to thank the IRAM staff for their support during the observation campaigns. This work is also based on observations made with the NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. These observations are associated with program nos. 1181, 1895, 2514, and 3577. The NIKA2 dilution cryostat has been designed and built at the Institut Néel. In particular, we acknowledge the crucial contribution of the Cryogenics Group, and in particular Gregory Garde, Henri Rodenas, Jean-Paul Leggeri, Philippe Camus. We acknowledge financial support from the “Programme National de Cosmologie and Galaxies” (PNCG) funded by CNRS/INSU-IN2P3-INP, CEA and CNES, France, from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (project CONCERTO, grant agreement No 788212) and from the Excellence Initiative of Aix-Marseille University-A*Midex, a French “Investissements d’Avenir” programme. This work has been partially funded by the Foundation Nanoscience Grenoble and the LabEx FOCUS ANR-11-LABX-0013. This work is supported by the French National Research Agency under the contracts “MKIDS”, “NIKA” and ANR-15-CE31-0017 and in the framework of the “Investissements d’avenir” program (ANR-15-IDEX-02). This work has been supported by the GIS KIDs. This work has benefited from the support of the European Research Council Advanced Grant ORISTARS under the European Union’s Seventh Framework Programme (Grant Agreement no. 291294). A. M. acknowledges support the funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 101098309 – PEBBLES) A. R. acknowledges financial support from the Italian Ministry of University and Research – Project Proposal CIR01_00010. E. A. acknowledges funding from the French Programme “Investissements d’avenir” through the Enigmass Labex. M. M. E. acknowledges the support of the French Agence Nationale de la Recherche (ANR), under grant ANR-22-CE31-0010. R. A. acknowledges support from the Programme National Cosmology et Galaxies (PNCG) of CNRS/INSU with INP and IN2P3, co-funded by CEA and CNES. R. A. was supported by the French government through the France 2030 investment plan managed by the National Research Agency (ANR), as part of the Initiative of Excellence of Université Côte d’Azur under reference number ANR-15-IDEX-01. The NIKA2 data were processed using the Pointing and Imaging In Continuum software (PIIC; Zylka 2013; Berta & Zylka 2019-2024), developed by Robert Zylka at the Institut de Radioastronomie Millimetrique (IRAM) and distributed by IRAM via the GILDAS pages. PIIC is the extension of the MOPSIC data reduction software to the case of NIKA2 data. This work made use of Astropy: (http://www.astropy.org) a community-developed core Python package and an ecosystem of tools and resources for astronomy (Astropy Collaboration 2013, 2018, 2022).

Note added in proof. In this note, we add an important clarification to what is described in Appendix C about the dust cosmic density computed by Péroux & Howk (2020). On the basis of relative metal abundances of damped Lyα systems (DLAs), these authors estimated the number of metal atoms missing from the gas, which were assumed to be incorporated into dust. By integrating the mass of these “depleted” elements relative to hydrogen, they determined the dust-to-gas ratio of each DLA absorber. The dust cosmic density was then estimated by applying to the gas density an average of these ratios, weighted by the HI column density. This approach allows for the possible variety of the dust-to-gas ratio of the absorbers. No correction to the Péroux & Howk (2020) values has been applied in Figs. 14 and C.1.

References

Appendix A: Comparison between the different fitting methods

The eight SED fitting codes used in this work are here compared in terms of their M_dust, M_⋆ and L_IR(8 − 1000μm) outputs. For each of these quantities, we directly compare the results of each method to those of each other approach.

For each object, we compute the dispersion of all eight determinations around their median value (i.e., computing the median absolute deviation, M.A.D., per object). To quantify the scatter in each combination of two methods, we compute the distribution of the median absolute deviation of the quantity (y − x)/x, computed for each pair of codes using all sources. The quantities x and y simply represent the M_dust or L_IR(8 − 1000μm) values computed with methods in the abscissa or ordinate axis in the comparison diagrams.

A.1. Dust mass

In order to produce a meaningful comparison, all dust masses have been re-scaled to the κ_ν normalization by Draine et al. (2014). We remind that the first DL07 modeling came with the parametrization of the dust absorption coefficient κ_ν by Li & Draine (2001), later revised by Draine (2003) (κ_850μm = 0.038 m² kg⁻¹). Draine et al. (2014) re-normalized their assumption on dust properties such that κ_850μm = 0.047 m² kg⁻¹. On the other hand, da Cunha et al. (2008, 2015) adopted the normalization by Dunne et al. (2000), that is κ_850μm = 0.077 m² kg⁻¹.

Figure A.1 compares the M_dust results of each method to those of each other approach. The comparison shows a general agreement between the different methods, modulo their natural limitations (especially those of the MBB models).

The most striking exception is the MBB-gen model, which systematically underestimates M_dust with respect to the other methods. The effect varies between 0.2 and 0.5 dex (median offset over all sources, with respect to different codes), corresponding to up to 4× the M.A.D. of the quantity (y − x)/x. This confirms the results by Ismail et al. (2023), who showed that a more thorough knowledge of the source properties (specifically the size of the dusty IR-emitting region) is needed to properly apply this method (Sect. 4.1).

The MBB-thin model shows a milder systematic under-estimation of M_dust, to be at least partially ascribed to the missing warmer dust component that dominates the λ_rest < 50 μm emission. A systematic offset between the original and the high-z versions of MAGPHYS and SED3FIT exists, most likely caused by the broader range of dust properties allowed in the latter, which leads to a better fit to the observed SEDs than in the original version. Moreover, MAGPHYS and SED3FIT (which are tightly related to each other) show the largest scatter when compared to the other methods. This drives the preference to use the result of the high-z version of these codes for the N2CLS sample.

A.2. Stellar mass

For the stellar component, the MAGPHYS-based codes adopt Bruzual & Charlot (2003, BC03) stellar populations with a Chabrier (2003) IMF and apply the Charlot & Fall (2000, CF00) two-components attenuation law. The CIGALE fits have been performed adopting a similar setup: the BC03 stellar populations with the Chabrier (2003) IMF and the same CF00 attenuation law.

Figure A.2 presents the comparison of stellar masses and shows a very good agreement, although a few important outliers exist, especially in the high-z versions of MAGPHYS and SED3FIT with respect to CIGALE. Section 5.3 show the consequence of these differences on the position of the N2GN galaxies in the M_⋆-SFR plane. Appendix B investigates the effect on depletion timescales. The most striking cases are very red galaxies at high redshift, with M_⋆ differences up to 1 dex between MAGPHYS and CIGALE. In these cases, despite adopting the same stellar models and the same dust attenuation law, MAGPHYS converges to best fit models with NIR luminosities much brighter than those found by CIGALE. The SEDs of the majority of these cases miss data points to constrain the rest-frame NIR emission. The most likely culprits are differences in the adopted SFH and in the retrieved A_V. For these very red outliers, A_V > 3 and A_V(MAGP)−A_V(CIG) > 1.3. These differences testify how critical and uncertain can the M_⋆ derivation be in the case of extreme, high-z SMGs.

A.3. Infrared luminosity

As for the infrared luminosity, the single-temperature MBB models do not include any warm dust component and are therefore limited to data points at λ_rest > 50 μm. Consequently, the integrated luminosity of these models is L_FIR(50 − 1000μm), rather than L_IR. By integrating the templates family of Herschel star forming galaxies defined by Berta et al. (2013a), we recovered a typical ratio L_FIR/L_IR ≃ 0.7 for this kind of galaxies (Berta et al. 2023). We therefore rescaled the MBB L_FIR estimate to the galaxies’ L_IR using this ratio.

Figure A.3 compares the L_IR results. In the cases when an AGN-torus was necessary in the fit, its contribution to the L_IR has been excluded. A good overall agreement between the different approaches is found. CIGALE behaves naturally similarly to DL07 and is in good agreement with SED3FIT-hz. The extreme outliers in some of the panels are those sources for which the codes diverge because not enough data are available to constrain all modeled components.

Fig. A.1. Comparison of Mdust obtained with the eight different SED fitting methods considered here. The top histogram shows the distribution of the median absolute deviation (M.A.D.) of the quantity (y − x)/x, computed for each pair of codes using all sources. The bottom histogram depicts the dispersion of all eight determinations around their median value (i.e., computing the M.A.D. per object).

Fig. A.2. Comparison of M⋆ obtained with the five different UV-to-radio SED fitting methods considered here.

Fig. A.3. Comparison of LIR obtained with the eight different SED fitting methods considered here.

Appendix B: Effect of M_⋆ on the τ_dep estimate

In Sect. 5.3, the results of SED fitting have been used to derive the molecular gas depletion timescale of the N2GN galaxies, by applying the Tacconi et al. (2020) relation that links M_⋆, redshift and distance from the MS to τ_dep. This approach takes as reference the MS parametrization by Speagle et al. (2014), shown in Fig. 10. The τ_dep thus obtained depends on the value of M_⋆. In the same Sect. 5.3, it has also been shown that for a few extreme galaxies at z > 4, with very red rest-frame optical colors, MAGPHYS/SED3FIT gives a significantly different estimate of M_⋆ with respect to CIGALE. These sources are outliers in Fig. A.2.

Here we show in Fig. B.1 the consequence of using the value of M_⋆ obtained with MAGPHYS/SED3FIT or with CIGALE. Few sources migrate to lower values of τ_dep, and from the locus of the MS to the region occupied by starbursts. This happens specifically for the group of galaxies at z ∼ 5, that is going to be the subject of Lagache et al. (in prep.).

Fig. B.1. Consequence of changing M⋆ when computing τdep with the Tacconi et al. (2020) scaling relations. Left panel: computation performed adopting the MAGPHYS/SED3FIT fiducial results. Right panel: same result obtained adopting the CIGALE estimate of M⋆. The references of the literature data and lines are listed in the caption of Fig. 9

Appendix C: The importance of κ_ν when comparing different works

Section 7.1 presented the evolution of the dust cosmic density computed by integrating the STY results of the N2GN DMF in the redshift range z = 1.6 to 7.2. The N2GN results have been there compared to literature data, pointing out the need to reconcile all data to the same dust model normalization, namely to the same κ_ν. Our choice of reference is κ_850μm = 0.047 m² kg⁻¹ (Draine et al. 2014).

If a different dust model was adopted, a re-normalization of dust masses is needed, because M_dust depends inversely on κ_ν. This is easy to visualize in the simple case of the MBB (Eq. 3). The re-normalization factor is simply the ratio between the two different κ_ν values of reference (at the same reference frequency). Consequently, the observed DMF is shifted (in log space) by the same amount along the mass axis and so does its integral.

For example, if a work found in the literature adopted the Draine (2003)κ₈₅₀ = 0.038 m² kg⁻¹ value, in order to compare the literature dust masses to N2GN’s, they would need to be rescaled by a factor κ₈₅₀(D03)/κ₈₅₀(D14), because dust mass is inversely proportional to κ_ν. In other words, if the adopted κ₈₅₀ is smaller than the Draine et al. (2014) value, the corresponding dust masses need to be scaled down, in order to be compared to those obtained with the Draine et al. (2014) normalization. Vice versa if the adopted κ₈₅₀ is larger than the Draine et al. (2014) value, they will need to be scaled up. This holds if all adopted dust models differ only by their normalization but have the same functional dependence on wavelength (e.g. the R_V = 3.1 MW models by Draine and collaborators). On the other hand, if more complex differences exist in the dust model assumptions, further transformations might be needed. The specific κ_ν assumptions of the different works shown in Fig. 14 are as follows.

The SED fitting code MAGPHYS (da Cunha et al. 2008) was used by Beeston et al. (2018, local galaxies), Beeston et al. (2024, 0.0 < z < 0.5), Driver et al. (2018, optically selected galaxies at z < 1.75), and Dudzevičiūtė et al. (2021, sub-millimeter galaxies at z = 1 − 2 and z = 3 − 4). MAGPHYS adopts κ₈₅₀ = 0.077 m² kg⁻¹. Other works in the literature collection shown in Fig. 14 also adopted the same value, namely Vlahakis et al. (2005) and Eales & Ward (2024). Dunne et al. (2011) adopted κ₂₅₀ = 0.89 m² kg⁻¹, which they state corresponds to κ₈₅₀ = 0.077 m² kg⁻¹.

Pozzi et al. (2020, 2021) used the value κ_1.2THz = 4.0 cm² g⁻¹, corresponding to Draine (2003, κ₈₅₀ = 0.038 m² kg⁻¹). Magnelli et al. (2020) adopted κ₈₅₀ = 0.043 m² kg⁻¹. Traina et al. (2024a) estimated the dust mass of their 189 galaxies in three ways: a multi-wavelength SED fit with CIGALE, and a RJ flux scaling based on two different T_dust assumptions of 25 and 35 K. In all three cases, they adopted the Draine et al. (2014)κ_ν re-normalization.

Ménard & Fukugita (2012) computed the dust density studying MgII absorbing clouds in the lines of sight of quasars, applying the Weingartner & Draine (2001) SMC (Small Magellanic Cloud) dust model. Transforming to the MW dust model implies multiplying by a factor 1.8 (see their Sect. 4). A further rescaling by 0.93 brings the value to the Draine (2003) normalization and finally this needs to be transformed to the Draine et al. (2014) assumption. Finally, Péroux & Howk (2020) derived the dust cosmic density by applying the dust to gas ratio by Draine (2003) – with the correction prescribed by Draine et al. (2014) – to the gas cosmic density. Therefore, their data do not need any further correction.

The consequence of rescaling some results to a value of κ_ν different from the one originally adopted is a shift of the DMF in logarithmic space to larger or smaller masses. The net result on ρ_dust (the integral of the DMF) is therefore a simple scaling by the same amount as for M_dust.

Table C.1 summarizes the values adopted in literature, as described above. Figure C.1 shows the importance of reconciling all the data to the same κ_ν assumptions. In the left panel the data are shown “as they are” computed by each work, each with its specific value of κ_ν. The right panel is a repetition of Figs. 14 and 15, in full color and excluding the models, after rescaling all results to the Draine et al. (2014) choice of reference.

Fig. C.1. Effect of rescaling ρdust to the same κν. Left: data without rescaling; each data set is shown as obtained with its κν underlying assumption. Right: the same data after rescaling to κ850 = 0.047 m2 kg−1 (Draine et al. 2014).

Appendix D: ρ_dust versus z

In this brief Appendix, Figs. 14 and 15 are presented with redshift on the main x-axis instead of look-back time, as is often found in the literature.

Fig. D.1. Cosmic density of dust in star forming galaxies as a function of redshift. Figures 14 and 15 present the same data as a function of look-back time.

Appendix E: Notes on individual sources

Figure E.1 (available on Zenodo, as mentioned in the “Data availability” Section) presents the multi-wavelength postage stamps, the FIR SED fitting and the panchromatic optical-to-radio modeling of all N2GN sources. Table E.1 lists their main properties: names, identification, redshift, and SED fitting results. In this Appendix, some notes about the N2GN individual sources are given.

N2GN_1_01 is an optically dark galaxy, identified with the GN10 SCUBA source (Pope et al. 2006) at z = 5.303. It was detected at 2 mm by GISMO on the IRAM 30 m telescope (GDF2000.6; Staguhn et al. 2014). It is detected by the SMA (Cowie et al. 2017) and NOEMA (Riechers et al. 2020). The latter authors studied the ¹²CO spectral line energy distribution (SLED) of this source and its dynamical properties, by combining NOEMA and VLA observations of ¹²CO and [CII] 158 μm. GN10 is a compact starburst (∼1.6 kpc in diameter) hosted in a ∼6.4 kpc rotating cold gas disk, with a total dynamical mass of 8.6 × 10¹⁰ M_⊙. The observed SED includes data from JWST, HST, Spitzer, Herschel, NIKA2, NOEMA and VLA.

N2GN_1_02’s SMA position is ∼2 arcsec away from the SCUBA source GOODS 850.6 (Wang et al. 2004; Barger et al. 2012) and 0.25 arcsec from AzGN5 (Chapin et al. 2009). It was previously assigned a low photometric redshift (Cowie et al. 2017), likely related to the nearby galaxy lying at only ∼6 arcsec away. We adopt the z_phot = 1.963 reported by Kodra et al. (2023) and based on optical-NIR SED fitting.

N2GN_1_03 was detected by SCUBA2 at 450 μm and corresponds to GOODS 850-9 and GN19 (Pope et al. 2006), with a spectroscopic redshift z = 2.484. It lies close to a pair of radio sources. The SMA position identifies it with the southwestern one. It was also detected by IRS at 16 μm with a flux of 32.6 μJy Teplitz et al. (2011). At its position, the 2 Ms Chandra map of the CDFN reveals a X-ray source with a hard spectrum and fluxes of 5.94 and 1.56 10⁻⁵μJy in the hard (2-7 keV) and soft (1.2-2 keV) band, respectively (Alexander et al. 2003; Evans et al. 2024). We fit the SED including an AGN-torus component, which contributes to the UV (and X-rays) emission of the galaxy but is negligible in the FIR-millimeter domain. The VLA radio measurements are in agreement with the radio-FIR correlation, given the Herschel, SCUBA2 and NIKA2 flux densities.

N2GN_1_04 is identified with GN20 (Pope et al. 2005; Daddi et al. 2009), a well studied z = 4.055 starburst system of two interacting galaxies. The NIKA2 source coincides with the brightest of the two (GN20 proper), while GN20.2 is too far for contributing to the NIKA2 flux density. Cortzen et al. (2020) studied the [CI](2-1) to (1-0) ratio, as measured with NOEMA, finding that its excitation temperature points to an optically thick dusty medium with optical depth reaching unity at 170 ± 23 μm. Recent JWST MIRI and NIRSpec observations reveal a conspicuous central source hosted by an extended envelope, faint stellar clumps associated to CO similar structures, a turbulent rotating disk (500 km/s) with additional non circular motions, and a very broad Hα line in the central regions, suggestive of strong, possibly AGN driven, winds (Colina et al. 2023; Übler et al. 2024; Crespo Gómez et al. 2024).

N2GN_1_05 corresponds to the source HDF254 observed at 1.4 GHz at high resolution with MERLIN by Bothwell et al. (2010). The radio observations show five peaks within a 1.5 arcsec maximum distance from each other, with a large velocity dispersion (249 ± 18 km/s) indicative of a possible merger. A strong MERLIN point source 0.5 arcsec to the west might be an AGN component, although the optical and MIR spectroscopy show a normal Hα to [NII] line ratio and strong PAH features (Swinbank et al. 2004; Pope et al. 2008). Our W21CV NOEMA program detected an emission line at 154.2 GHz, consistent with being ¹²CO(4-3) at z = 1.996, the same redshift measured by optical spectroscopy (3D-HST; Skelton et al. 2014, and those previously mentioned). The proxy position lies between the two GISMO 2 mm sources GDF2000.3 and GDF2000.8 (Staguhn et al. 2014). Although the source is detected in the 2 Ms CDFN X-ray map (Barro et al. 2019; Alexander et al. 2003; Bothwell et al. 2010), no AGN-torus component is needed to fit its optical to millimeter SED, nor the VLA continuum shows any excess with respect to the radio-FIR correlation.

N2GN_1_06 is HDF850.1, the first and brightest sub-millimeter galaxy discovered by SCUBA in the Hubble Deep Field (Hughes et al. 1998). It was also detected at 2 mm by GISMO (GDF2000.1; Staguhn et al. 2014). It is an optically dark galaxy and it took nearly 14 years to determine its redshift. Walter et al. (2012) and Neri et al. (2014) placed it at z = 5.184, detecting CO and [CII] 158 μm at millimeter wavelengths with the PdBI, precursor of NOEMA. HDF850.1 is a starburst galaxy belonging to the known z ∼ 5.2 overdensity in the GOODS-N field. JWST detected its rest-frame UV-optical emission, solving the long-lasting puzzle of its identification (Sun et al. 2024). The galaxy is mildly lensed by a z = 1.224 elliptical galaxy, with a magnification factor μ ∼ 2.5 The NIRCAM photometry reveals a rest-frame UV excess with respect to the light emitted from the extinguished stars that power its bright IR emission. We fit the SED of this galaxy adding a unextinguished young SSP to reproduce this excess (Sects. 4.5 and 5.4).

N2GN_1_07 corresponds to GN12 (Pope et al. 2005; Daddi et al. 2009; Mancini et al. 2009). It was detected at 2 mm also by GISMO (GDF2000.11; Staguhn et al. 2014). NIRSpec Wide GTO data give z_spec = 2.9601. The source lies in a very crowded field, with the consequence that the Herschel/SPIRE photometry possibly includes the contribution of a few strongly blended sources and was therefore not used in our SED modeling. We kept instead the PACS photometry, that has a PSF size matching well the NIKA2 beam, as well as the 450 μm SCUBA2 measurement (Barger et al. 2022). The SED is well reproduced with a starburst model; the available VLA data are in slight excess with respect to the radio-FIR correlation.

N2GN_1_08 is detected by SMA at a position 0.15 arcsec away from AzGN2 (Chapin et al. 2009) and ∼4 arcsec away from JWST JADES-GN-189.13100+62.28684 (z = 5.267; Sun et al. 2024). A new JWST extraction reveals a large obscured galaxy at a projected distance of 0.04 arcsec, composed of four “blobs” (M. Xiao, priv. comm.). We used CIGALE to derive the photometric redshift of the latter, z_phot = 4.34 ± 0.16.

N2GN_1_09 is the ID12646 in the Spitzer-driven super-deblended catalog by Liu et al. (2018), and the AzGN10 source detected by Atzec at 1.1 mm (Perera et al. 2008; Chapin et al. 2009), pinpointed by VLA. Its position is 0.465 arcsec away from the IRAC J123627.54+621218.4 object reported in Enia et al. (2022). It is an optically dark galaxy and was observed with NOEMA by Jin et al. (2022) and Bing (2022) at z_spec = 4.147. Its optical counterpart is finally identified on JWST maps with a very red galaxy at ∼0.3 arcsec distance.

N2GN_1_10 is a known bright AGN at z = 2.005, detected in the Chandra 2 Ms map with fluxes of 2.92, 3.66 and 18.10 × 10⁻⁵ μJy in the 0.5-1.2, 1.2-2.0 and 2.0-7.0 keV bands, respectively (Evans et al. 2024). It is also detected by WISE and Akari in the MIR and at 8.5, 5.5 and 1.4 GHz in the radio (Cutri et al. 2013; Pearson et al. 2010; Richards et al. 1998; Guidetti et al. 2017; Cowie et al. 2017). The MIR SED is a well defined power-law and is reproduced with an important AGN-torus component in the SED fitting. The radio emission, nevertheless, is consistent with the radio-FIR correlation, indicating that the AGN does not contribute significantly in this spectral domain.

N2GN_1_11 is an optically dark source not present in the Barro et al. (2019) catalogs. It lies at 0.223 arcsec from the IRAC J123713.86+621826.6 object reported by Enia et al. (2022). It has a bright radio counterpart, which is here used to pinpoint its precise position and corresponds to GN40 or AzGN26. To date, no optical photometry is available. The photometric redshift z = 3.99 was estimated by Wang et al. (2016a) based on F160W, Ks and Spitzer photometry. The object is a bright X-ray source with fluxes of 5.1, 7.9 and 10.7 10⁻⁵ μJy in the 0.2-0.5, 0.5-1.2 and 2.0-7.0 Chandra bands. The SED fitting includes a bright AGN-torus component that dominates the IRAC 8 μm, IRS 16 μm and MIPS 24 μm observed fluxes. The 21 cm emission is in excess with respect to the radio-FIR correlation by at least a factor ∼5.

N2GN_1_12 is a blend of two PdBI 255 GHz sources (Bing 2022). The two components, named 12a and 12b, correspond to ID13161 and ID13207 in Liu et al. (2018) and are detected by both HST and JWST. The coordinates of 12a are 0.7 arcsec away from those of GN11 (Pope et al. 2005, 2006; Mancini et al. 2009). The high-resolution optical data show two peaks each, belonging to the same galaxies. The photometric redshift of 12a is determined here with CIGALE (z = 2.59 ± 0.08) and that of 12b was estimated by Kodra et al. (2023, z = 4.163). Because of blending, the FIR-millimeter emission is defined by the NOEMA fluxes only. Source 12b is also detected at 21 cm and is consistent with the radio-FIR correlation.

N2GN_1_13 is an optically dark galaxy. Its position and spectroscopic redshift has been determined by recent NOEMA observations with ¹²CO(4-3) and (5-4) transitions at z = 5.181 (this work and Lagache et al. in prep.). It is now detected by JWST/FRESCO (M. Xiao, priv. comm.). No MIR data are available to verify the possible presence of an AGN-torus component contributing to its MIR SED. The optical photometry reveals a rest-frame UV excess, that we reproduce by adding a unextinguished young SSP to the SED modeling (Sect. 5.4).

N2GN_1_14 is identified with a faint optical galaxy at z_spec = 1.76 thanks to the SMA observations (Cowie et al. 2017). Its coordinates match those of GN17. The SED is very complete and finely sampled with 27 data points, a well-defined 1.6 μm (rest-frame) stellar “bump”, and data on both sides of the FIR emission peak. It is detected in the X-rays (Barro et al. 2019; Alexander et al. 2003), but no AGN-torus component is needed to reproduce its optical-to-radio emission.

N2GN_1_15, similarly to the previous one, is identified by SMA and has a very well defined SED at z_spec = 2.0. It corresponds to GN6. The deep optical data show an excess that we model with an unextinguished young SSP.

N2GN_1_16 is a high-redshift galaxy, identified by SMA and detected in the optical by HST (Kodra et al. 2023; Liu et al. 2018). It lies at the position of GN9. Including JWST data (one source at a distance of 0.23 arcsec, M. Xiao priv. comm.) we obtain a photometric redshift of z_phot = 3.47 ± 0.45.

N2GN_1_17 is split in two components by NOEMA (Sect. 2.3), named 17a and 17b. The component 17b lies very close to a bright star, and therefore no optical-IR photometry is available, nor a redshift estimate. Because of blending, the SPIRE, SCUBA2 and NIKA2 data are not considered and the FIR emission of 17a is defined by NOEMA alone. Liu et al. (2018) deblended the MIPS and IRAC emission. Its position matches that of GN37. The spectroscopic redshift, z_spec = 3.190, was measured by Cowie et al. (2004).

N2GN_1_18 is an optically dark galaxy. Cowie et al. (2017) associated it to a Keck/LRIS and DEIMOS redshift but it has no K band counterpart on the CFHT maps; therefore we ignored it. In the JWST grism spectra, it has a bright neighboring source in the same dispersion direction. Therefore, these spectra were of no use either. Thanks to the NOEMA data, we identify it on the JWST maps and derived a new z_phot = 4.39 ± 0.15 with CIGALE.

N2GN_1_19 is matched through SMA to a X-ray source at z = 2.578. Its integrated 0.3-8.0 keV flux is 1.68 × 10⁻⁴ μJy (Wang et al. 2016a). It lies at the same position of GN4 and AzGN16 and shows an 8.0 μm excess with respect to its stellar emission but no radio excess. We model the SED including an AGN-torus component.

N2GN_1_20 is another optically dark galaxy, detected by SMA. It lies at a distance of 1.3 arcsec from AzGN4 (Chapin et al. 2009). It corresponds to ID14914 by Liu et al. (2018) and has a Ks band flux of 0.75 μJy, resulting in a very red Ks to 3.6 μm color, S_3.6/S_Ks ∼ 13, interpreted as a deep Balmer break combined to a large extinction, that leads to a z_phot = 5.333 (Liu et al. 2018). In our SED fitting, the bright 24 μm flux is contributed by a strong 3.3 μm PAH feature.

N2GN_1_21 is identified with a z_spec = 1.219 radio source at a projected distance of 0.25 arcsec from the N2CLS position, with a rich and very well sampled SED. Barro et al. (2019) matched it to an X-ray source (Alexander et al. 2003), but no AGN-torus component is needed in the SED fitting. A JADES z = 5.446 source (Sun et al. 2024) lies at a projected distance of 1.05 arcsec from the N2CLS position.

N2GN_1_22 lies in a crowded area. The SMA counterpart seems blended with a different bright 24 μm object. Hence, we adopted the MIPS and IRS fluxes deblended by Liu et al. (2018) and we do not consider the Herschel photometry. Kodra et al. (2023) give a spectroscopic redshift z_spec = 2.098.

N2GN_1_23 is an optically dark galaxy detected by JWST/FRESCO at z_spec = 5.179 (Xiao et al. 2023). The NOEMA identification is at a distance of only 0.167 arcsec. The MIR and FIR corresponding sources are blends of several objects, and therefore we ignored them for the SED. The NIKA2 position is 0.336 arcsec away from the GOODS J123656.5+621207 objects mentioned by Enia et al. (2022, with z_phot = 3.03).

N2GN_1_24 is a blend of two millimeter sources detected by NOEMA at ∼10 arcsec distance from each other and named 24a and 24b. In their SEDs, the FIR emission is defined by NOEMA alone because of blending. The NOEMA positions allow the sources to be matched with two optical galaxies at z_phot = 2.697 and 2.817 (Kodra et al. 2023).

N2GN_1_25 was observed with NOEMA and identified as an optical galaxy at z_spec = 2.914. The exquisite SED shows clear Lyman and Balmer breaks. In the Barro et al. (2019) catalog it is matched to an X-ray source, but neither the MIR or the radio emission is in excess with respect to the SED of a star forming galaxy.

N2GN_1_26 is identified with a high-z HST galaxy thanks to SMA data. It corresponds to AzGN18 or GN38. Kodra et al. (2023) give a photometric redshift z_phot = 5.31. The FIR SED is well defined by the PACS, SCUBA2, SMA, and NIKA2 data. The 24 μm flux is well matched by an intense 3.3 μm PAH feature.

N2GN_1_27 is the merger of two galaxies (27a and 27b), identified with NOEMA and at nearly the same redshift (z_spec = 1.988, 1.989; Kodra et al. 2023). The system lies at the position of GN7. It has been detected by IRS at 16 μm (Teplitz et al. 2011), for which we adopt the deblended fluxes by Liu et al. (2018).

N2GN_1_28 is identified with VLA 21 cm source (Owen 2018) and has a spectroscopic redshift z_spec = 3.2223 from the NIRSpec Wide GTO. The very red optical SED is well fitted by an intense starburst model and the VLA flux is consistent with the radio-FIR correlation.

N2GN_1_29 is also identified through VLA data and is at the position of GN18. A JADES z = 5.432 source lies at a projected distance of 3.9 arcsec (Sun et al. 2024), but the z = 2.819 radio source has a higher match probability. The SED is well sampled and is well fitted with a starburst model.

N2GN_1_30 is very close to N2GN_1_51; therefore we excluded the SPIRE photometry that is partially blended. An AGN-torus component is required to fit the optical blue slope, as a pure stellar model does not reproduce successfully both the 8 μm and the observed U, F435W and F606W photometry at the same time; in the best fit solution, the AGN emission supplies the optical and the stellar component dominates in the NIR.

N2GN_1_31 has a bright VLA 1.4 GHz counterpart, also detected by HST and JWST. It corresponds to AzGN20 and its redshift is z_phot = 2.638 (Kodra et al. 2023). In the SPIRE 500 μm channel, the source is heavily blended with another, brighter, object; therefore we ignore this band.

N2GN_1_32 is a radio source at z_spec = 3.652 (Owen 2018; Barger et al. 2008; Kodra et al. 2023). It is detected in the Chandra 2 Ms map of the CDFN with fluxes of 2.42, 2.23, 0.95, and 1.11 ×  10⁻⁴ μJy in the 0.2-0.5, 0.5-1.2, 1.2-2.0, and 2.0-7.0 keV bands, respectively (Evans et al. 2024). The observed rest-frame NIR-MIR SED is a power law and is dominated by an AGN-torus component in the IRAC and MIPS bands. A stellar component contributes the bulk of the optical emission. The radio emission does not show any excess and is consistent with the radio-FIR correlation.

N2GN_1_33 corresponds to the SMA GOODS-N source nr. 13 and SCUBA2 nr. 63 (Cowie et al. 2017), as well as to GN25. Three sources lie inside the SMA beam: a bright z_spec = 1.013 optical galaxy, a very faint optical galaxy revealed by JWST, and a radio source which lies in the middle of the two optical sources. The bright and faint optical galaxies are respectively at 1.21 and 0.56 arcsec from the SMA position. We chose the closest one as the identification. Kodra et al. (2023) reports a redshift z_phot = 5.388. Because of blending, we ignore the Spitzer and Herschel data in the SED fitting. Therefore, a large gap in the wavelength coverage exists between the rest-frame optical and sub-millimeter spectral domains.

N2GN_1_34 corresponds to AzGN28 (Chapin et al. 2009) but breaks in two components (34a and 34b), when observed with NOEMA. Both are optically dark sources. The 34a component lies very close to a bright nearby galaxy, it might be lensed, but no counterpart is found in any catalog or map. On the other hand, 34b is an isolated source detected in the IRAC, IRS 16 μm, MIPS 24 μm channels, and at 1.4 GHz by VLA. Enia et al. (2022) report an optically dark source at 0.233 arcsec from its position (GOODS J123644.0+621938). CIGALE gives a redshift z_phot = 4.12 ± 0.53 taking into account the available photometry (Sect. 3.3).

N2GN_1_35 is identified with VLA with an edge-on z_phot = 1.95 galaxy. The SED is superb, with 26 bands covering the whole electromagnetic spectrum from 3500 Å to 21 cm and allows a textbook fit.

N2GN_1_36 is an optically dark source that lies very close (< 5 arcsec) to a bright star. The SMA data pinpoint its position. It lies at 0.233 arcsec from the GOODS J123658.5+620931 source discussed in Enia et al. (2022). Liu et al. (2018) derive a redshift z_phot = 3.55, but their photometry is affected by the nearby star. Here we adopt a new HST extraction by the FRESCO team (M. Xiao, priv. comm.) and derive a new estimate of z_phot = 3.93 ± 0.17 with CIGALE. Similarly to the case of N2GN_1_30, an AGN-torus component is needed: the bright NIR emission is dominated by stars and the blue excess by the AGN.

N2GN_1_37 is an isolated galaxy at z_spec = 2.302 (Kodra et al. 2023), whose accurate position is given by the VLA data. Its observed SED consists of 24 data points. Barro et al. (2019) match it to an X-ray source in the 2 Ms Chandra map Alexander et al. (2003), but its optical-to-radio SED is well reproduced by a simple starburst model.

N2GN_1_38 lies in a crowded region of GOODS-N. The NOEMA position matches that of a z_phot = 3.126 optical galaxy detected by HST and Spitzer/IRAC (Kodra et al. 2023; Liu et al. 2018). The FIR-millimeter emission is defined by the 850 μm to 2.0 mm data only, because of the heavy blending affecting the Herschel data.

N2GN_1_39 is a z_spec = 2.301 galaxy, identified thanks to the NOEMA data. Its SED shows a clear 1.6 μm stellar peak and a well defined FIR dust emission sampled by PACS, SCUBA2, NIKA2 and NOEMA. It is detected in the X-rays (Barro et al. 2019; Alexander et al. 2003) but no AGN torus is needed in the SED fitting. The radio measurement is consistent with being powered by star formation only.

N2GN_1_40 is a galaxy at z_phot = 2.764 (Kodra et al. 2023). In the 1.2 mm NOEMA map, two millimeter sources are detected, but the fainter one lies 7.87 arcsec apart and does not contribute significantly to the NIKA2 flux. The Herschel and MIPS photometry is affected by a significant blending and is therefore ignored. It is clearly detected in the HST and JWST maps at the same NOEMA position. The observed SED has a large wavelength gap between the observed NIR (IRAC) and sub-millimeter data, implying that the rest-frame MIR and the peak of the FIR emission are not well sampled.

N2GN_1_41 has a very red, but complex optical SED, with a slight flux offset between the HST and JWST measurements. The redshift z_phot = 4.101 by Kodra et al. (2023) is based on the HST, Ks and IRAC photometry. Its position is pinpointed by our NOEMA data. It was not detected by Herschel, nor by MIPS.

N2GN_1_42 is in a crowded area and the actual multi-wavelength counterparts are identified thanks to the NOEMA data. The IR deblended photometry by Liu et al. (2018) is used. It is a z_phot = 2.228 optically bright galaxy showing a distinct break at ∼1.5 μm in the observed frame and a strong 1.6 μm stellar peak in the IRAC 5.8 μm band. The 1.4 GHz measurement is in excess to the radio synchrotron powered by star formation by more than an order of magnitude, indicating the presence of an AGN, but no signs are detected at MIR or optical wavelengths.

N2GN_1_43 is an optically and K-band dark galaxy, detected by JWST at ∼10σ in the F277W band, ∼3σ in F200W, and tentatively by HST at shorter wavelengths (M. Xiao, priv. comm.). NOEMA matches it with a clear detection in the Spitzer IRAC bands, but the object is not visible in the MIPS and Herschel maps. CIGALE places it at z_phot = 5.87 ± 0.38 with the Balmer break between the two JWST measurements mentioned here.

N2GN_1_44 is identified as the GNz7q source studied by Fujimoto et al. (2022). These authors report on the results by NOEMA and HST, highlighting its compact morphology at rest-frame UV wavelength, a bright FIR-millimeter emission, and a MIR excess in the observed 24 μm band. The object is not detected in the Chandra 2 Ms X-ray map. Its rest-frame UV spectrum does not show any bright, nor broad, line between Lyα and C III]λ1909. NOEMA detected the [CII] 158 μm line at redshift z = 7.1899, consistent with the observed Lyman break. Two possible interpretations of the emission of GNz7q are proposed by Fujimoto et al. (2022): a UV-compact star-forming object or (their favorite scenario) a red quasar explaining the MIR excess at odds with the lack of broad UV MgII and CIV lines. Our SED fitting reproduces the source with a blue star forming galaxy model, characterized by a prominent 3.3 μm PAH feature that accounts for the observed MIR excess. On the other hand, SED3FIT and CIGALE failed at finding a viable solution including a type-1 AGN torus, as suggested by Fujimoto et al. (2022, their Fig. 2).

N2GN_1_45 is pinpointed by the VLA data (Owen 2018) and matches a z_phot ∼ 3 galaxy detected by HST, JST, Spitzer, Herschel and SCUBA2. It corresponds to AzGN24 and GN23 (Chapin et al. 2009; Mancini et al. 2009; Pope et al. 2006). Barro et al. (2019) identify it with an X-ray source. No AGN torus is needed to describe the observed SED, because the 1.6 μm stellar peak falls in the Spitzer/IRAC 8.0 μm channel and the stellar component fits well the IRAC monotonically increasing SED. The VLA radio flux is in significant excess to the radio-FIR correlation, as applied by MAGPHYS, but consistent with the Delhaize et al. (2017) correlation scaled from the MBB fit.

N2GN_1_46 is matched by NOEMA to one of three tightly grouped optical sources. It was previously erroneously associated with a z_spec = 0.55 galaxy (Kodra et al. 2023; Barger et al. 2012; Owen 2018), but the high-resolution HST and JWST optical photometry reveals that this is a nearby galaxy at ∼1.5 arcsec distance and not the NOEMA source. Moreover this low redshift is not consistent with the observed 1.6 μm stellar peak detected between the IRAC 4.5 and 5.8 μm bands. Therefore we adopted the z_phot = 1.894 redshift estimate by Kodra et al. (2023). The source is detected in the 2 Ms Chandra map with fluxes of 0.99 and 7.6 × 10⁻⁵ μJy in the 0.5-2.0 and 2.0-8.0 keV bands. In the SED fitting we include an AGN-torus component, that turns out to contribute significantly only at the bluest optical wavelengths.

N2GN_1_47 is a high-z galaxy detected by HST and identified via the NOEMA observations. The observed photometry shows a break between the Ks and IRAC 3.6 μm band, placing it at z_phot = 4.87 (Kodra et al. 2023). Barro et al. (2019) associate it with an X-ray source (Alexander et al. 2003), but no AGN component is required to fit its SED.

N2GN_1_48 is detected by HST and is identified by the SMA, It corresponds to GN21. JWST spectroscopy sets the redshift to z_spec = 2.6 (JADES, DJA). The SED is reproduced by a very extinguished starburst model.

N2GN_1_49 lies close to a very bright star. Identified on the basis of its VLA radio counterpart (Owen 2018), is blended to another galaxy in the IRAC bands (Liu et al. 2018). Kodra et al. (2023) derive a z_phot = 2.592. The optical SED is characterized by a blue excess, that we fit with an unextinguished young SSP (Sect. 5.4).

N2GN_1_50 corresponds to AzGN22, so far associated to a radio galaxy at z = 1.0662. The VLA position is at 0.19 arcsec from a JWST source at z_spec = 3.1312, that nicely matches the observed 1.6 μm stellar peak detected in the IRAC 5.8 μm channel and the Balmer break in the F105W JWST band. The adopted IRAC photometry is from the super-deblended catalog by Liu et al. (2018), that removes the contribution of a close by (∼3 arcsec away) face-on spiral galaxy (see HST postage stamps, arguably the z = 1.0662 galaxy). Closer to the N2CLS position, there is a JWST galaxy at z_spec = 5.186 (Sun et al. 2024),

N2GN_1_51 lies roughly two NIKA2 HPBW away from another, much brighter 1.2 mm source, in a crowded region with several MIR objects (on the IRAC and MIPS maps). NOEMA confirms the NIKA2 detection and identifies it with a HST object at z_phot = 2.013 that becomes progressively more point-like in the bluest HST bands. Barro et al. (2019) matches it to an X-ray source in the Alexander et al. (2003) CDFN catalog. Despite the compact blue-optical morphology and the X-rays detection, no AGN component is needed to fit the optical-to-millimeter SED of this source.

N2GN_1_52 is an optically dark source and lies outside of the current JWST covered area. Only a very faint Spitzer photometry is available, identified with the NOEMA positioning. Liu et al. (2018) give a redshift z_phot = 4.06. The IRAC 3.6 μm flux is brighter than the one at 4.6 μm, which could be explained with a bright Hα line at this redshift, contributing to the broad-band photometry.

N2GN_1_53 is an isolated IRAC galaxy at z_phot = 1.994, located by NOEMA. The rich optical SED has a remarkable blue excess, that we reproduce adding an unextinguished young SSP to the dusty starburst model.

N2GN_1_54 is confirmed by NOEMA, that defined its precise position. It is detected on the Spitzer IRAC and MIPS maps, as well as by HST and JWST. Its SED shows a remarkable Lyman break at z_phot = 2.885. The Barro et al. (2019) catalog includes a X-ray identification (Alexander et al. 2003). A big gap between the NIR and sub-millimeter observed photometry still exists and does not allow us to verify if an AGN torus is needed to fit the MIR SED of this source.

N2GN_1_55 lies on top of three radio sources (Owen 2018): two very bright nearby galaxies (at 3.8 and 4.6 arcsec from the NIKA2 position) and one optically and IRAC-dark object (at 0.08 arcsec from the NIKA2 position). The SCUBA2 data were matched to one of the nearby galaxies because the SCUBA2 position is shifted to the north with respect to the NIKA2 position. We adopt the dark source as the counterpart. The Spitzer map being highly blended in this region, currently no photometry other than NIKA2 and VLA is available, nor any redshift estimate is at hand.

N2GN_1_56 splits in two NOEMA sources at a distance of 13.5 arcsec from each other, named 56a and 56b. The former corresponds to AzGN21 (Chapin et al. 2009). Both sources have IRAC, MIPS and HST counterparts. Their redshifts are z_phot = 1.866 and 2.598, respectively (Kodra et al. 2023). Because of blending, the Herschel photometry is ignored. Barro et al. (2019) matched the 56a counterpart to an X-ray source. The 56b galaxy shows a significant optical-blue excess that is modeled with a unextinguished, young SSP as usual. The radio emission is consistent with the radio-FIR correlation in both cases.

N2GN_1_57 is detected by SCUBA at 450 and 850 μm (Barger et al. 2012). The corresponding SCUBA2 source was matched to a z = 0.9 radio galaxy. In the NIKA2 beam lies the high-z candidate HRG14 J123731.66+621616.7 (Finkelstein et al. 2015), also detected by the SHARDS middle-band survey at z_phot = 5.3 (Kodra et al. 2023). This is our counterpart choice. It is detected in the four IRAC channels. At longer wavelengths it is strongly blended to a couple of nearby objects, and therefore we ignored the MIPS and Herschel photometry. Summarizing: the observed SED includes HST and IRAC photometry, with a pronounced Balmer break at λ ∼ 2 μm, as well as SCUBA2 and NIKA2 data.

N2GN_1_58 is identified by VLA with a z_phot = 2.723 galaxy with an exceptional SED coverage. A total of 20 data points cover the wavelength range from 3500 Å to 21 cm, excluding the SPIRE data that are a blend of different objects.

N2GN_1_59 is also pinpointed by VLA, although the radio source is at 5.74 arcsec from the N2CLS position. It was observed by NOEMA as part of our follow-up program, but only a 5σ source at 9.7 arcsec distance from the NIKA2 peak was detected. The high-z candidate HRG14 J123702.11+621733.4 (z_phot = 3.434; Arrabal Haro et al. 2018) lies at 0.79 arcsec from the radio position, and CANDELS-GDN J123702.1+621734.2 is 0.3 arcsec away. The final, observed SED combines the HST, Spitzer, NIKA2 and VLA data and is best fitted with a very young starburst (1.3 × 10⁷ yr).

N2GN_1_60 is blended with a few large nearby galaxies and is identified by SMA at z_phot = 1.811 (Kodra et al. 2023). It is detected both in the Ks band and at 1.4 GHz (Owen 2018). The IRAC photometry is finely deblended by Liu et al. (2018), but the MIPS, PACS and SPIRE data are not recoverable. HST constrains the optical emission of this red, dusty galaxy.

N2GN_1_61 is an optical- and IRAC-dark source, confirmed by NOEMA. JWST grism observations detected one emission line, consistent with Hα at z_spec = 5.201 (M. Xiao, priv. comm.). A faint source is visible on the IRAC 3.6 μm and 4.5 μm maps at ∼3 arcsec distance. The observed SED consists of only three data points (NIKA2 1.2 mm, SCUBA2 850 μm and JWST F444W), but the availability of the redshift allowed us to perform a fit with the sole aim of deriving M_dust.

N2GN_1_62 is a nearby spiral galaxy at z_spec = 0.641, confirmed by NOEMA. Its SED is well sampled by all available instruments.

N2GN_1_63 is identified by VLA with a galaxy at z_spec = 2.21 (Kodra et al. 2023) and corresponds to GN5. The SED is defined by 23 bands from the optical to the radio domain and shows a clear 1.6 μm stellar peak, a very red NIR-to-MIR color (L_7.5μm/L_2.5μm ∼ 12 in the rest frame) and a well define FIR-millimeter emission with data on both the Wien and RJ tails.

N2GN_2_13 is a NIKA2 2 mm source with no 1.2 mm detection (in the blind catalog with S/N> 3σ). It was observed by NOEMA and located in a crowded region with several optical sources. The NOEMA position corresponds to a z_phot = 2.071 galaxy (Kodra et al. 2023) detected by HST. We use the super-deblended IRAC catalog by Liu et al. (2018), but MIPS and SPIRE fluxes are not recoverable.

N2GN_2_20, finally, is an optically dark source, identified by VLA with a z_phot = 6.67 isolated galaxy. The photometric redshift is constrained mainly by the Balmer break observed between the Ks band and the IRAC 3.6 μm band. Its position corresponds to GN3 (Pope et al. 2006).

All Tables

All Figures

	Fig. 1. NOEMA 150 GHz maps for N2GN_1_13 (left) and N2GN_1_17 (right) that reveal one and two counterparts, respectively. The beam shown at the bottom left has a size of 1.61″ × 0.93″. Multi-wavelength postage stamps, including NIKA2, are shown in Fig. E.1.
In the text

	Fig. 2. Comparison between NIKA2 and NOEMA flux densities for the sources observed in program W21CV. The green solid diagonal lines marks the 1:1 correspondence between the two instruments.
In the text

	Fig. 3. Differences between the coordinates of the N2GN sources and their counterparts. Left panel: Distance between the NIKA2 coordinates (all from the 1.2 mm position but two) and the matched proxy (Sect. 3.1). Right panel: Distance between the matched proxy and the N2GN optical counterparts. The green dashed lines and boxes mark the rms of the distribution.
In the text

	Fig. 4. Redshift distribution of the N2GN sources. The solid histogram includes spectroscopic redshifts, while the open histogram contains all sources in the sample.
In the text

	Fig. 5. Spectral energy distributions of all N2GN sources obtained with the UV-to-radio SED fitting (gray lines). The red and blue lines represent the median and average SEDs, respectively, as obtained by combining all models.
In the text

Fig. 6. Distribution of the N2GN galaxies, their luminosity, and dust mass as a function of redshift. Top panel: Redshift distribution. The filled histogram includes the sources with spectroscopic redshift, while the open histogram includes all sources. Second panel: Infrared luminosity of the star forming component, integrated between 8 and 1000 μm. The blue solid line represents the expected luminosity for a source of 0.7 mJy in the NIKA2 1.2 mm band, assuming its emission is described by the median SED shown in the left panel of Fig. 5. The shaded area is given by the best fit models of all sources. Third: 850 μm monochromatic luminosity. Bottom panel: Dust mass distribution as a function of redshift. The dashed and dotted lines represent the minimum dust mass and maximal-completeness dust mass obtained rescaling the L850 versus z line adopting the minimum and maximum Mdust/L850 ratio in the sample (Sect. 6.1).

In the text

Fig. 7. Distribution of the N2GN galaxies as a function of dust mass and 1.2 mm flux density, color coded on the basis of their redshift. Median error bars are shown in the bottom-right corner. The red dashed line is a least squares fit to the data. The red dotted lines represent the same line scaled to the 2.5th and 97.5th percentiles of the residuals distribution. The red, solid vertical lines mark the flux of the sources with no redshift measurement and multi-wavelength counterparts available.

In the text

	Fig. 8. Study of the Mdust/L850 ratio. Top-left panel: least squares fit to the correlation between Mdust and L850. Right-hand panels: Mdust and Mdust/L850 as a function of IR luminosity.
In the text

Fig. 9. Depletion timescale of the N2GN galaxies obtained applying the Tacconi et al. (2020) scaling relation. The gray shaded area is the trend found by Saintonge et al. (2013). The different lines represent the trends found by Tacconi et al. (2020) for: main sequence galaxies (selected within ΔMS = ±0.6 dex, black solid line); starburst galaxies (ΔMS > 0.6 dex, blue dashed line); extreme starbursts (ΔMS > 1.2 dex, purple dotted line); below-MS galaxies (ΔMS < 0.4 dex, red long-dashed line). Literature data are from Alaghband-Zadeh et al. (2013), Aravena et al. (2016, 2014, 2013), Bakx et al. (2020), Berta et al. (2023), Bothwell et al. (2017, 2013), Carilli et al. (2010), Chung et al. (2009), Combes et al. (2011, Combes et al. (2013, Dannerbauer et al. (2019), Decarli et al. (2016, 2019), Dunne et al. (2021, 2020), Fujimoto et al. (2017), Freundlich et al. (2019), Geach et al. (2011), Genzel et al. (2015, 2003), George et al. (2013), Hagimoto et al. (2023), Harris et al. (2012, 2010), Ivison et al. (2013, 2011, 2010), Penney et al. (2020), Riechers et al. (2020, 2011), Rudnick et al. (2017), Sharon et al. (2016), Solomon et al. (1997), Tacconi et al. (2018, 2013), Thomson et al. (2012), Valentino et al. (2018), Villanueva et al. (2017), Wang et al. (2018), and Yang et al. (2017).

In the text

Fig. 10. Position of the N2CLS GOODS-N galaxies in the M⋆ versus SFR plane, as obtained with the SED fitting results, in different redshift bins. Each galaxy is positioned on the basis of its Δ(log(SFR))MS computed at its actual redshift with respect to the Popesso et al. (2023) MS, translated to the average redshift of each bin. The solid blue lines represent the parametrization of the MS of star forming galaxies by Popesso et al. (2023). For comparison the solid red lines refer to the one by Speagle et al. (2014). The dotted lines are placed at ±0.5 dex from the solid ones. Top row: results obtained with the fiducial fit (MAGPHYS or SED3FIT high-z). Bottom row: results obtained with CIGALE for comparison.

In the text

Fig. 11. Example of fit uncertainty when no NIR and MIR data are available and a big gap in wavelength exists between the rest-frame optical and FIR spectral domains. The two different tones of red shaded areas represent the 2.5%, 16% 84% and 97.5% percentiles of the models, as computed in the photometric bands of the input catalog; the corresponding M⋆ range is as large as 0.25 dex (2.5th to 97.5th percentiles).

In the text

Fig. 12. Exemplification of how the dust mass completeness is computed in the redshift bin 1.6 < z ≤ 2.4 (Sect. 6.1). The diagonal shaded areas represent tracks of the N2GN SED models, given the minimum and maximum Mdust/L850 ratios in the sample. The vertical dashed line is the 1.2 mm flux limit of the survey, and the horizontal dashed line is the corresponding mass maximal completeness threshold. The dotted horizontal lines mark the boundaries of an example mass bin, in which the hatched area would contain the sources missed by the flux cut.

In the text

Fig. 13. Dust mass function of the N2GN sources (red symbols). Only the mass bins containing at least two galaxies and with a completeness ≥50% are plotted. Uncertainties take into account ΔMdust, including SED fitting parameters sampling and Δz contributions, and Poisson errors. Left panel: comparison to literature results by Pozzi et al. (2020, blue filled squares and long-dashed line) and Traina et al. (2024a, light blue triangles and short-dashed lines, fiducial results in the 2.0 < z ≤ 2.5, 3.5 < z ≤ 4.5 and 4.5 < z ≤ 6.0 redshift bins). The blue dotted lines are the same Pozzi et al. (2020) DMF, evolved to z = 2.0 and 3.3 (see Sect. 6.2). Right panel: results of the STY analysis. The green solid lines and shaded areas are the result of the STY parametric analysis: most probable model (solid lines) and 1σ uncertainty (shaded areas). The dashed green lines represent the 1.6 < z ≤ 2.4 result, repeated in the other two redshift bins for comparison.

In the text

Fig. 14. Evolution of the dust mass density as a function of look-back time. The N2GN results are depicted as red filled circles. All other symbols belong to literature data, as indicated in the right-hand panel. Local estimates are by Vlahakis et al. (2005) and Beeston et al. (2018); combinations of optical-NIR and FIR-sub-millimeter data are by Eales & Ward (2024) and Driver et al. (2018); Herschel-based studies are by Beeston et al. (2024), Pozzi et al. (2020) and Dunne et al. (2011); sub-millimeter and millimeter driven results are by Pozzi et al. (2021), Dudzevičiūtė et al. (2021), and Magnelli et al. (2020); and finally derivations of ρdust are based on Ωgas (Péroux & Howk 2020) and MgII absorbers (Ménard & Fukugita 2012). All data have been rescaled to our chosen cosmology and κν (Draine et al. 2014). Figure C.1 presents a full-color version and Fig. D.1 shows ρdust as a function of redshift.

In the text

Fig. 15. Comparison of the dust mass density to model predictions. The N2GN results are depicted as red filled circles. The different lines represent the semi-analytical models by Popping et al. (2017), Parente et al. (2023) and Yates et al. (2024), the chemical evolution model by Gioannini et al. (2017), and the cosmological hydrodynamical simulations by Lewis et al. (2023), Li et al. (2019) and Aoyama et al. (2018). Figure D.1 presents a version of this figure as a function of redshift.

In the text

Fig. A.1. Comparison of Mdust obtained with the eight different SED fitting methods considered here. The top histogram shows the distribution of the median absolute deviation (M.A.D.) of the quantity (y − x)/x, computed for each pair of codes using all sources. The bottom histogram depicts the dispersion of all eight determinations around their median value (i.e., computing the M.A.D. per object).

In the text

	Fig. A.2. Comparison of M⋆ obtained with the five different UV-to-radio SED fitting methods considered here.
In the text

	Fig. A.3. Comparison of LIR obtained with the eight different SED fitting methods considered here.
In the text

	Fig. B.1. Consequence of changing M⋆ when computing τdep with the Tacconi et al. (2020) scaling relations. Left panel: computation performed adopting the MAGPHYS/SED3FIT fiducial results. Right panel: same result obtained adopting the CIGALE estimate of M⋆. The references of the literature data and lines are listed in the caption of Fig. 9
In the text

	Fig. C.1. Effect of rescaling ρdust to the same κν. Left: data without rescaling; each data set is shown as obtained with its κν underlying assumption. Right: the same data after rescaling to κ850 = 0.047 m2 kg−1 (Draine et al. 2014).
In the text

	Fig. D.1. Cosmic density of dust in star forming galaxies as a function of redshift. Figures 14 and 15 present the same data as a function of look-back time.
In the text