# A panchromatic view of N2CLS GOODS-N: The evolution of the dust cosmic density since $z \sim 7$

S. Berta, G. Lagache, A. Beelen, R. Adam, P. Ade, H. Ajeddig, S.

Amarantidis, P. André, H. Aussel, A. Benoît, et al.

► To cite this version:

S. Berta, G. Lagache, A. Beelen, R. Adam, P. Ade, et al. A panchromatic view of N2CLS GOODS-N: The evolution of the dust cosmic density since $z \sim 7$. Astronomy & Astrophysics - A&A, 2025, 696, pp.A193. 10.1051/0004-6361/202452894. hal-05040382

HAL Id: hal-05040382

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

Submitted on 19 Apr 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# A panchromatic view of N2CLS GOODS-N: The evolution of the dust cosmic density since $z \sim 7$

S. Berta¹,*, G. Lagache², A. Beelen², R. Adam³, P. Ade⁴, H. Ajeddig⁵, S. Amarantidis⁶, P. André⁵, H. Aussel⁵, A. Benoît⁷, M. Bethermin⁸, L.-J. Bing⁹, A. Bongiovanni⁶, J. Bounmy¹⁰, O. Bourrion¹⁰, M. Calvo⁷, A. Catalano¹⁰, D. Chérouvrier¹⁰, L. Ciesla², M. De Petris¹¹, F.-X. Désert¹², S. Doyle⁴, E. F. C. Driessen¹, G. Ejlali¹³, D. Elbaz⁵, A. Ferragamo¹¹,¹⁴, A. Gomez¹⁵, J. Goupy⁷, C. Hanser¹⁰, S. Katsioli¹⁶,¹⁷, F. Kéruzoré¹⁸, C. Kramer¹, B. Ladjelate⁶, S. Leclercq¹, J.-F. Lestrade¹⁹, J. F. Macías-Pérez¹⁰, S. C. Madden⁵, A. Maury²⁰,²¹,⁵, F. Mayet¹⁰, H. Messias²²,²³, A. Monfardini⁷, A. Moyer-Anin¹⁰, M. Muñoz-Echeverría²⁴, I. Myserlis⁶, R. Neri¹, A. Paliwal¹¹, L. Perotto¹⁰, G. Pisano¹¹, N. Ponthieu¹², V. Revéret⁵, A. J. Rigby²⁵, A. Ritacco²⁶,²⁷, H. Roussel²⁸, F. Ruppin²⁹, M. Sánchez-Portal⁶, S. Savorgnano¹⁰, K. Schuster¹, A. Sievers⁶, C. Tucker⁴, M.-Y. Xiao³⁰, and R. Zylka¹

*(Affiliations can be found after the references)*

Received 5 November 2024 / Accepted 27 February 2025

## ABSTRACT

To understand early star formation, it is essential to determine the dust mass budget of high-redshift galaxies. Sub-millimeter rest-frame emission, dominated by cold dust, is an unbiased tracer of dust mass. The New IRAM KID Arrays 2 (NIKA2) conducted a deep blank field survey at 1.2 and 2.0 mm in the GOODS-N field as part of the NIKA2 Cosmological Legacy Survey (N2CLS), detecting 65 sources with $S/N \ge 4.2$. Thanks to a dedicated interferometric program with NOEMA and other high-angular resolution data, we identified the multi-wavelength counterparts of these sources and resolved them into 71 individual galaxies. We built detailed spectral energy distributions (SEDs) and assigned a redshift to 68 of them over the range $0.6 < z < 7.2$. We fit these SEDs using modified blackbody and Draine & Li (2007, ApJ, 657, 810) models and the panchromatic approaches MAGPHYS, CIGALE, and SED3FIT, thus deriving their dust mass ($M_{dust}$), infrared luminosity ($L_{IR}$), and stellar mass ($M_*$). Eight galaxies require an active galactic nucleus torus component, and another six require an unextinguished young stellar population. A significant fraction of our galaxies are classified as starbursts based on their position on the $M_*$ versus star formation rate plane or their depletion timescales. We computed the dust mass function in three redshift bins ($1.6 < z \le 2.4, 2.4 < z \le 4.2$ and $4.2 < z \le 7.2$) and determined the Schechter function that best describes it. The dust cosmic density, $\rho_{dust}$, increases by at least an order of magnitude from $z \sim 7$ to $z \sim 1.5$, as predicted by theoretical works. At lower redshifts, the evolution flattens. Nonetheless, significant differences exist between results obtained with different selections and methods. The superb GOODS-N data set enabled a systematic investigation into the dust properties of distant galaxies. N2CLS holds promise for combining these deep field findings with the wide COSMOS field into a self-consistent analysis of dust in galaxies both near and far.

**Key words.** evolution – galaxies: evolution – galaxies: high-redshift – galaxies: luminosity function, mass function – galaxies: statistics – submillimeter: galaxies

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

<footnote>* Corresponding author; berta@iram.fr</footnote>
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 \sim 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 \sim 2$. The former computed the DMF of more than $5 \times 10^5$ 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 \le 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 \sim 7$ to $z \sim 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_0 = 67.4 \text{ km s}^{-1} \text{ Mpc}^{-1}$, $\Omega_m = 0.315$, and $\Omega_\Lambda = 0.685$ (Planck Collaboration VI 2020), a Chabrier (2003) stellar initial mass function (IMF), and the Draine (2003) frequency-dependent dust absorption coefficient, $\kappa_\nu$, renormalized as indicated by Draine et al. (2014, i.e., $\kappa_{850\mu m} = 0.047 \text{ m}^2 \text{ kg}^{-1}$ as reference).

All the products from this paper have been released online on the survey's home page<fnref n="1" />. 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

<footnote n="1">¹ Data release: https://data.lam.fr/n2cls/home</footnote>
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<fnref n="2" />) 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<fnref n="2" /> (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.

The NOEMA data were reduced and calibrated in the standard way with the GILDAS<fnref n="3" /> 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.

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 \sim 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

<footnote n="2">² Pointing and Imaging In Continuum, found at https://publicwiki.iram.es/PIIC and https://www.iram.fr/IRAMFR/GILDAS/</footnote>

<footnote n="3">³ Grenoble Image and Line Data Analysis Software, https://www.iram.fr/IRAMFR/GILDAS/</footnote>
<table><caption>**Table 1.** Statistics of the multi-wavelength counterparts adopted in the N2GN catalog.</caption><thead><tr><th>Survey</th><th>Instrument and band</th><th>Nr. sources</th><th>Ref.</th></tr></thead><tbody><tr><td rowspan="7">N2CLS GOODS-N</td><td>Total</td><td>71</td><td>This work</td></tr><tr><td>With z</td><td>68</td><td>This work</td></tr><tr><td>NIKA2</td><td>65</td><td>Bing et al. (2023)</td></tr><tr><td>1.2 mm</td><td>63</td><td>Bing et al. (2023)</td></tr><tr><td>2.0 mm</td><td>26</td><td>Bing et al. (2023)</td></tr><tr><td>NOEMA 1.2 or 2.0 mm</td><td>31 (*)</td><td>W21CV, this work</td></tr><tr><td>NOEMA 1.2 or 2.0 mm</td><td>2</td><td>W16EE, E16AI, Bing (2022)</td></tr><tr><td rowspan="5">Others sub-mm to radio</td><td>SMA 860 µm</td><td>22</td><td>Cowie et al. (2017)</td></tr><tr><td>SCUBA2 850 µm</td><td>46</td><td>Cowie et al. (2017)</td></tr><tr><td>SCUBA2 450 µm</td><td>31</td><td>Barger et al. (2022)</td></tr><tr><td>NOEMA</td><td>5</td><td>Aa.Vv. (Sect. 2.4)</td></tr><tr><td>VLA 1–2 GHz</td><td>53</td><td>Owen (2018)</td></tr><tr><td>CANDELS/GOODS</td><td>Optical U to K</td><td>59</td><td>Barro et al. (2019)</td></tr><tr><td>FRESCO, JADES, PANORAMIC, Congress</td><td>JWST F090W to F770W</td><td>24</td><td>New extraction</td></tr><tr><td>CANDELS/GOODS</td><td>IRAC 3.6 to 8.0 µm</td><td>59</td><td>Barro et al. (2019), Liu et al. (2018)</td></tr><tr><td>PeakUp IRS</td><td>IRS 16 µm</td><td>10</td><td>Liu et al. (2018)</td></tr><tr><td>GOODS</td><td>MIPS 24 µm</td><td>43</td><td>Barro et al. (2019), Liu et al. (2018)</td></tr><tr><td>PEP + GOODS-H</td><td>PACS 100, 160 µm</td><td>28</td><td>Magnelli et al. (2013), Liu et al. (2018)</td></tr><tr><td>GOODS-H/HerMES</td><td>SPIRE 250–500 µm</td><td>22</td><td>Elbaz et al. (2011), Oliver et al. (2012), Liu et al. (2018)</td></tr><tr><td>CDFN</td><td>Chandra 0.5 to 7.0 keV</td><td>5</td><td>Evans et al. (2024)</td></tr></tbody></table>

Notes. (*)26/27 NIKA2 detected/targeted; five are split in two components by NOEMA.

<figure><img src="image_2.png" /><figcaption>**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.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Comparison between NIKA2 and NOEMA flux densities for the sources observed in program W21CV. The green solid diagonal lines mark the 1:1 correspondence between the two instruments.</figcaption></figure>

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 \ge 4$, over an area of 450 arcmin$^2$. 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\sigma$ 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\sigma$. 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

**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 Multitiered 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).

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

**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<fnref n="4" />), 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 SEExtractor (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.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.35 is appropriate for the typical population of massive dusty sources at high redshift ($z > 3$). We observe that some of our $z \sim 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).

<footnote n="4">⁴ https://dawn-cph.github.io/dja/index.html</footnote>
<figure><img src="image_4.png" /><figcaption>**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.</figcaption></figure>

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<fnref n="5" />) 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<fnref n="6" />.

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

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

<footnote n="5">⁵ http://ned.ipac.caltech.edu/. NED is funded by the National Aeronautics and Space Administration and operated by the California Institute of Technology.</footnote>

<footnote n="6">⁶ These are N2GN_1_33 and N2GN_1_55, as reported in Appendix E.</footnote>
<figure><img src="image_5.png" /><figcaption>Fig. 4. Redshift distribution of the N2GN sources. The solid histogram includes spectroscopic redshifts, while the open histogram contains all sources in the sample.</figcaption></figure>

identified source. None of the SMA identifications is multi-ple. 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

$$d = \sqrt{\text{HWHM}_{\text{proxy}}^2 + \text{HWHM}_{\lambda}^2}, \quad (1)$$

where HWHMproxy 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 Extra-galactic Sky (SIDES) simulations (Béthermin et al. 2017, and in prep.), we clearly see an excess of galaxies at $z \sim 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.).

# 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_{\star}$ 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

$$L = \Omega \epsilon_{\nu} B_{\nu}(T_{\text{dust}}), \quad (2)$$

where $\Omega$ is the solid angle of emission, $\epsilon_{\nu}$ is the dust emissivity coefficient and $B_{\nu}(T_{\text{dust}})$ is the Planck function of temperature
$T_{\text{dust}}$. For a uniform medium of optical depth $\tau_{\nu}$, radiative transfer theory leads to $\epsilon_{\nu} = 1 - \exp(-\tau_{\nu})$. Relating $\tau_{\nu}$ to the mass absorption coefficient $\kappa_{\nu}$ 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):

$$L_{\nu} = 4\pi M_{\text{dust}}\kappa_{\nu}B_{\nu}(T_{\text{dust}}), \quad (3)$$

which is used to model the observed SEDs. The absorption coefficient depends on frequency as a power law: $\kappa_{\nu} = \kappa_{0}(\nu/\nu_{0})^{\beta}$. We adopt the values of $\kappa_{\nu}$ tabulated by Draine (2003), with the correction indicated by Draine et al. (2014). In this assumption, the reference value is $\kappa_{850\,\mu\text{m}} = 0.047\,\text{m}^2\,\text{kg}^{-1}$.

The MBB-thin fit includes three free parameters, namely normalization, $T_{\text{dust}}$ and $\beta$. The dust mass of the given galaxy, $M_{\text{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_{\text{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 $\beta$ is not necessarily the same as the one in the tabulated $\kappa_{\nu}$. 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(-\tau_{\nu})$, which can be expressed in terms of the wavelength $\lambda_{\text{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 $\lambda_{\text{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 $\lambda_{\text{thick}}$ as free parameter; the other based on the size A of the dust emitting region. The two are linked by the expression $\lambda_{\text{thick}} = \lambda_{0}(\kappa_{0}M_{\text{dust}}/A)^{1/\beta}$. With the use of mock simulated SEDs, Ismail et al. (2023) demonstrated that, adopting the first formulation, $\lambda_{\text{thick}}$ is hardly constrained. Furthermore, these authors showed that a reliable estimate of $T_{\text{dust}}$, $\beta$ and $M_{\text{dust}}$ with the general MBB requires an independent knowledge of A. The degeneracies that affect $\lambda_{\text{thick}}$ can induce a significant underestimation of $M_{\text{dust}}$ of 20% to 50%, depending on dust mass itself.

The MBB fit was limited to data at rest-frame wavelengths $\ge 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(\nu) \propto \nu^{\alpha}$ (spectral index $\alpha = -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_{\text{min}}$. The latter, representing a fraction $\gamma$ of the total amount of dust, is exposed to starlight with intensities within the range $U_{\text{min}}$ to $U_{\text{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_{\text{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 $-\alpha$.

The DL07 model thus has six free parameters: $U_{\text{min}}$, $U_{\text{max}}$, $\gamma$, $q_{\text{PAH}}$, $\alpha$, and $M_{\text{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 $\kappa_{\nu}$ 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 t^{-\alpha} \exp(-t/\tau)$, where $t$ is the model age and $\tau$ 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.
<figure><img src="image_6.png" /><figcaption>**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.</figcaption></figure>

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 \propto v^{-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\sigma$ 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_{\text{dust}}$, infrared luminosity, $L_{\text{IR}}$, and stellar mass, $M_{\star}$.

## 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^{11}$ 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.
<figure><img src="image_7.png" /></figure>

**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 $L_{850}$ versus z line adopting the minimum and maximum $M_{\text{dust}}/L_{850}$ ratio in the sample (Sect. 6.1).

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_{\text{dust}}$, the infrared luminosity, $L_{\text{IR}}$, and the stellar mass, $M_{\star}$, 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_{\text{IR}}$ (left) and by their $M_{\text{dust}}$ (right). The normalization by $M_{\text{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_{\text{dust}}/L_{\text{sub-mm}}$. Consequently, because of this relation between sub-millimeter emission and $M_{\text{dust}}$, dust masses derived by SED fitting without the knowledge of the actual rest-frame sub-millimeter emitted luminosity are subject to large uncertain-

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

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

ties (e.g., Berta et al. 2016) and/or need important assumptions in the SED modeling (e.g., fixing $\beta$ in a MBB analysis).

The two middle panels of Figure 6 present the integrated infrared luminosity ($L_{\text{IR,SF}}$, 8–1000 µm, excluding the possible AGN-torus contribution), and the monochromatic 850 µm luminosity of the N2GN galaxies ($L_{850}$), 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_{850}$: 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_{\text{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_{850}$, related to the variety of observed SED colors (Fig. 5).

The bottom panel of Fig. 6 reports the distribution of $M_{\text{dust}}$ as a function of redshift. Thanks to the negative sub-millimeter k-correction, there is no evident dependence of $M_{\text{dust}}$ on redshift. The dust masses of the N2GN galaxies span over the range between ~$10^8$ and ~$6 \times 10^9 M_{\odot}$.

## 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_{\text{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_{\text{dust}}$ on millimeter flux of all other N2GN galaxies.

Figure 7 shows the derived $M_{\text{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 $\pm 1.96 \times$ 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
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Study of the $M_{\text{dust}}/L_{850}$ ratio. Top-left panel: least squares fit to the correlation between $M_{\text{dust}}$ and $L_{850}$. Right-hand panels: $M_{\text{dust}}$ and $M_{\text{dust}}/L_{850}$ as a function of IR luminosity.</figcaption></figure>

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.

## 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_{\text{dust}}/L_{850\,\mu\text{m}}$ ratio, based on the $\kappa_{\nu}$ normalization by Draine et al. (2014) as usual.

We quantified the correlation between $M_{\text{dust}}$ and $L_{850}$ with a simple least squares fit to the data, leading to the relation $\log(M_{\text{dust}}) = (1.24 \pm 0.08)\log(L_{850}) + (12.11 \pm 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_{\text{dust}}$ and $L_{\text{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_{\text{dust}}/L_{850}$ and $L_{\text{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_{850}$ and $L_{\text{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_{\text{mol}}$) and molecular gas depletion timescale ($\tau_{\text{dep}}$), we computed $\tau_{\text{dep}}$ of the

<figure><img src="image_10.png" /><figcaption>**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 $\Delta_{\text{MS}} = \pm 0.6$ dex, black solid line); starburst galaxies ($\Delta_{\text{MS}} > 0.6$ dex, blue dashed line); extreme starbursts ($\Delta_{\text{MS}} > 1.2$ dex, purple dotted line); below-MS galaxies ($\Delta_{\text{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, 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).</figcaption></figure>

N2GN galaxies. Figure 9 compares the result to a collection of galaxies from the literature with $M_{\text{mol}}$ derived from CO data (Tacconi et al. 2020; Berta et al. 2023, and references therein).

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

Noteworthy, the stellar mass of several N2GN high-redshift galaxies is of the order of a few $10^{11} M_{\odot}$. 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).
<figure><img src="image_11.png" /><figcaption>Fig. 10. Position of the N2CLS GOODS-N galaxies in the M<sub>⋆</sub> 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))<sub>MS</sub> 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.</figcaption></figure>

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<sub>⋆</sub> 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 τ<sub>dep</sub> of adopting the M<sub>⋆</sub> 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<sub>⋆</sub> differences is a large difference of A<sub>V</sub> 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<sub>⋆</sub> > 10<sup>11</sup> M<sub>☉</sub> 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.

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

<figure><img src="image_12.png" /><figcaption>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<sub>⋆</sub> range is as large as 0.25 dex (2.5th to 97.5th percentiles).</figcaption></figure>

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<sub>dust</sub>, also known as DMF, is given by

$$ \Phi(M_{\text{dust}})\Delta M_{\text{dust}} = \sum_{i} \frac{1}{V_a^i} \Delta M_{\text{dust}}, \qquad (4) $$

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

$$ V_a = \frac{\Omega}{4\pi} \int_{z_{\text{min}}}^{\min(z_{\text{max}}, z_{\text{max}}^{\text{bin}})} \frac{dV}{dz}, \qquad (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
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Exemplification of how the dust mass completeness is computed in the redshift bin $1.6 < z \le 2.4$ (Sect. 6.1). The diagonal shaded areas represent tracks of the N2GN SED models, given the minimum and maximum $M_{dust}/L_{850}$ 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.</figcaption></figure>

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 N1KA2 map is not uniform, but it varies by up to a factor of three across the N2GN field. As a consequence, $\Omega$ 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_{850}$ versus $z$ line adopting the minimum and maximum $M_{dust}/L_{850}$ 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_{850}$ 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_{850}$ ratio in the N2GN sample, given the scat-

ter 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 \le 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_{850}$ 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.

In a given mass bin, $M_{dust}^{lo} < M_{dust} < M_{dust}^{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

$$f_{\text{compl}} = \frac{\int_{S_{\text{lim}}}^{S_{\text{up}}} N(M, S')dS'}{\int_{S_{\text{lo}}}^{S_{\text{up}}} N(M, S')dS'} \quad (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 $\ge 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 \le 2.4$, $2.4 < z \le 4.2$, and $4.2 < z \le 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.

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 $\kappa_y(250 \mu m) = 4.0 \text{ cm}^2 \text{ g}^{-1}$ (Bianchi 2013), which corresponds to the Draine (2003) dust models normalization. Therefore, we rescaled their DMF to our
<table><caption>**Table 2.** The N2GN DMF.</caption><thead><tr><th>log (M<sub>dust</sub>[M<sub>☉</sub>])</th><th>Center</th><th>Range</th><th>dN/dM<sub>dust</sub><br>[10<sup>-5</sup> h<sup>3</sup> Mpc<sup>-3</sup> dex<sup>-1</sup>]</th><th>N<br>srcs.</th><th>M<sub>dust</sub><br>Compl.</th></tr></thead><tbody><tr><td colspan="6" style="text-align: center;">1.6 &lt; z &le; 2.4</td></tr><tr><td>8.50</td><td>8.38</td><td>8.62</td><td>13.18 &plusmn; 8.98</td><td>5</td><td>0.55</td></tr><tr><td>8.75</td><td>8.63</td><td>8.88</td><td>12.41 &plusmn; 7.21</td><td>5</td><td>0.78</td></tr><tr><td>9.00</td><td>8.88</td><td>9.12</td><td>5.85 &plusmn; 3.12</td><td>6</td><td>1.00</td></tr><tr><td>9.25</td><td>9.12</td><td>9.38</td><td>2.93 &plusmn; 2.17</td><td>3</td><td>1.00</td></tr><tr><td colspan="6" style="text-align: center;">2.4 &lt; z &le; 4.2</td></tr><tr><td>8.30</td><td>8.15</td><td>8.45</td><td>6.91 &plusmn; 3.71</td><td>5</td><td>0.54</td></tr><tr><td>8.60</td><td>8.45</td><td>8.75</td><td>6.27 &plusmn; 2.73</td><td>10</td><td>0.76</td></tr><tr><td>8.90</td><td>8.75</td><td>9.05</td><td>3.16 &plusmn; 2.47</td><td>8</td><td>0.98</td></tr><tr><td>9.20</td><td>9.05</td><td>9.35</td><td>1.48 &plusmn; 1.10</td><td>4</td><td>1.00</td></tr><tr><td>9.50</td><td>9.35</td><td>9.65</td><td>0.69 &plusmn; 0.61</td><td>2</td><td>1.00</td></tr><tr><td colspan="6" style="text-align: center;">4.2 &lt; z &le; 7.2</td></tr><tr><td>8.18</td><td>8.04</td><td>8.32</td><td>2.64 &plusmn; 2.34</td><td>2</td><td>0.52</td></tr><tr><td>8.47</td><td>8.32</td><td>8.61</td><td>3.18 &plusmn; 1.73</td><td>6</td><td>0.86</td></tr><tr><td>8.76</td><td>8.61</td><td>9.91</td><td>1.28 &plusmn; 0.95</td><td>3</td><td>1.00</td></tr><tr><td>9.05</td><td>8.81</td><td>9.20</td><td>0.82 &plusmn; 0.55</td><td>3</td><td>1.00</td></tr></tbody></table>

**Notes.** The columns include: dust mass logarithmic bins (center, lower and upper boundaries); comoving number density; number of sources; and $M_{dust}$ fractional completeness (Sect. 6.1). The DMF uncertainties take into account the contribution of the $M_{dust}$ probability distribution function produced by the SED fitting, as well as redshift uncertainties and Poisson errors.

$\kappa_v$ of reference (Draine et al. 2014). The result is compared to the N2GN DMF at $1.8 < z \le 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^3$ 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^3$ 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

$$ \Phi(M_{\text{dust}}) \, \text{d} \log M_{\text{dust}} = \frac{\Phi_{\text{dust}}^*}{M_{\text{dust}}^*} \left(\frac{M_{\text{dust}}}{M_{\text{dust}}^*}\right)^{\alpha+1} e^{-\frac{M_{\text{dust}}}{M_{\text{dust}}^*}} \ln(10) \, \text{d} \log M_{\text{dust}} \quad (7) $$

Pozzi et al. (2020) fitted the observed DMF with a nonlinear 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 $\alpha$, $M_{\text{dust}}^*$ and $\Phi_{\text{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 $\alpha$ to the local value of 1.48, in order to derive the evolution of $M_{\text{dust}}^*$ and $\Phi_{\text{dust}}^*$.

We studied the variation of $M_{\text{dust}}^*$ as a function of cosmic age, based on Pozzi et al. (2020) results. By performing a simple least square fit to their $M_{\text{dust}}^*$ values as a function of cosmic time, we found the relation $\log M_{\text{dust}}^* = (-2.2 \pm 0.1) \times \log t_{\text{Univ}} + (30.1 \pm 1.0)$. The units are years for the age of the Universe and $M_{\odot}$ for $M_{\text{dust}}^*$.

While deriving this equation, we rescaled their values of $M_{\text{dust}}^*$ to the $\kappa_v$ normalization by Draine et al. (2014).

Using the derived equation, we evolved the Pozzi et al. (2020) $M_{\text{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 $\Phi_{\text{dust}}^*$ to match the observed N2GN DMF.

In both $1.6 < z \le 2.5$ and $2.5 < z \le 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 $\alpha = 1.48$ derived locally by Pozzi et al. (2020), suggesting that the $M_{\text{dust}}^*$ derived from their data is a good approximation at least up to $z \sim 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^3$ 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_{\text{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 ($\le 2-3\%$) on the mass function and on the value of $\rho$.

Following Berta et al. (2007, 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_{\text{dust}}^*$, $\Phi_{\text{dust}}^*$ and $\alpha$, 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.
<figure><img src="image_14.png" /><figcaption>**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 $\Delta M_{\text{dust}}$, including SED fitting parameters sampling and $\Delta 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 \le 2.5$, $3.5 < z \le 4.5$ and $4.5 < z \le 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\sigma$ uncertainty (shaded areas). The dashed green lines represent the $1.6 < z \le 2.4$ result, repeated in the other two redshift bins for comparison.</figcaption></figure>

<table><caption>**Table 3.** Results of the STY analysis of the N2GN DMF using a Schechter function.</caption><thead><tr><th>z</th><th>log M<sup>*</sup><sub>dust</sub></th><th>1σ range<br>log M<sub>☉</sub></th><th>α</th><th>1σ range</th><th>Φ<sup>*</sup><sub>dust</sub><br>10<sup>-5</sup> h Mpc<sup>-3</sup> dex<sup>-1</sup></th><th>1σ range<br>10<sup>5</sup> h M<sub>☉</sub> Mpc<sup>-3</sup></th><th>ρ<sub>dust</sub></th><th>1σ range<br>10<sup>5</sup> h M<sub>☉</sub> Mpc<sup>-3</sup></th></tr></thead><tbody><tr><td>1.6–2.4</td><td>9.35</td><td>9.22 to 9.64</td><td>-1.40</td><td>-1.74 to -1.26</td><td>4.36</td><td>1.35 to 7.47</td><td>1.43</td><td>0.98 to 2.87</td></tr><tr><td>2.4–4.2</td><td>9.61</td><td>9.46 to 9.73</td><td>-1.61</td><td>-1.83 to -1.41</td><td>0.69</td><td>0.37 to 1.39</td><td>0.64</td><td>0.40 to 1.11</td></tr><tr><td>4.2–7.2</td><td>9.56</td><td>9.31 to 9.79</td><td>-1.39</td><td>-1.91 to -1.26</td><td>0.34</td><td>0.07 to 1.15</td><td>0.18</td><td>0.13 to 0.37</td></tr></tbody></table>

## 7. Discussion

The variation of $M^*$ and $\Phi^*$, 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 \sim 0.5–1.0$ (Pozzi et al. 2020).

The characteristic dust mass, $M_{\text{dust}}^*$, does not evolve at least until redshift $z \sim 4$ and then is halved by $z \sim 2$. The evolution continues down to the local Universe, decreasing $M_{\text{dust}}^*$ by more than an order of magnitude until $z \sim 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{aligned} \rho_{\text{dust}} (M_{\text{dust}} > M_{\text{inf}}) &= \int_{M_{\text{inf}}}^{\infty} \Phi (M_{\text{dust}}) \, d \log M_{\text{dust}} \\ &= M_{\text{dust}}^* \Phi_{\text{dust}}^* \Gamma \left( \alpha + 2, \frac{M_{\text{dust}}}{M_{\text{dust}}^*} \right). \end{aligned} \qquad (8) $$

In its calculation, we adopted $M_{\text{inf}} = 10^4 M_{\odot}$. The resulting dust mass density of the N2GN survey is shown in Fig. 14 and compared to data found in the literature.

All literature data have been rescaled to the Draine et al. (2014) value of $\kappa_{850} = 0.047 \, \text{m}^2 \, \text{kg}^{-1}$. Appendix C lists the original assumptions of the different authors and shows the importance of rescaling all results to the same $\kappa_v$ in order to properly compare them to each other.

The local estimates of $\rho_{\text{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$^2$ 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 men-
tioned 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 unbi-
ased a M<sub>dust</sub> estimate as possible. Moreover these galaxies ben-
efit from a very rich multi-wavelength coverage, that allowed us
to avoid any strong assumption in the SED fitting that lead to
the M<sub>dust</sub> measurements (e.g., fixing the dust temperature). The
two samples lead to different results: in this work we detect a
strong evolution of ρ<sub>dust</sub> 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 ρ<sub>dust</sub> 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 ρ<sub>dust</sub> in Fig. 14,
and the apparently discordant trends of ρ<sub>dust</sub> 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 sam-

ples and in the methods adopted to derive M<sub>dust</sub> among dif-
ferent 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<sub>max</sub> 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 ρ<sub>dust</sub>
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).

Popping et al. (2017) included dust production and destruc-
tion 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 rela-
tion between stellar mass and dust mass in the local and high-
redshift Universe, as well as the dust-to-gas ratio of local galax-
ies 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 destruc-
tion processes. They explored different cosmological scenarios
of galaxy formation: pure luminosity evolution and number den-
sity 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 accre-
tion 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 hydro-
dynamic 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 cosmolog-
ical simulation code RAMSES-CUDATON (Ocvirk et al. 2016),
with the aim to study how dust affects the escape of ionizing pho-
tons 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 pre-
scriptions 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
<figure><img src="image_15.png" /><figcaption>**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 $\rho_{\text{dust}}$ are based on $\Omega_{\text{gas}}$ (Péroux & Howk 2020) and MgII absorbers (Ménard & Fukugita 2012). All data have been rescaled to our chosen cosmology and $\kappa_v$ (Draine et al. 2014). Figure C.1 presents a full-color version and Fig. D.1 shows $\rho_{\text{dust}}$ as a function of redshift.</figcaption></figure>

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 \rho_{\text{dust}}$ as determined by N2GN. The model by Gioannini et al. (2017) is consistent with higher measurements of $\rho_{\text{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 \sim 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 $\rho_{\text{dust}}(z)$, steadily decreasing from $z \sim 1$ to $\sim 0.2$ (Pozzi et al. 2020; Eales & Ward 2024). The model by Gioannini et al. (2017) fairly reproduces this decrease.

<figure><img src="image_16.png" /><figcaption>**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.</figcaption></figure>

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 $\rho_{\text{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 $\rho_{\text{dust}}$.

Given the large scatter in the data, after re-scaling all estimates to the same $\kappa_v$ assumptions, it is difficult to be conclusive on which model is most successful in reproducing the observed evolution of $\rho_{\text{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 $\kappa_v$), the eight methods produce consistent estimates of $M_{\text{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_{\text{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_{\text{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^{11} M_{\odot}$. 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 NIR 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_{\text{dust}}) = (1.24\pm0.08)\log(L_{850})+(12.11\pm0.20)$. This relation can be useful for galaxies lacking an SED coverage detailed enough to derive $M_{\text{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 $\Phi^*$, that doubles from $z \sim 7$ to $z \sim 2.5$, increases by a further factor of five at later epochs and ceases at $z \sim 0.5-1.0$; (b) a constant characteristic dust mass, $M_{\text{dust}}^*$, from $z \sim 7$ to $z \sim 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 $\rho_{\text{dust}}$ determinations to the same $\kappa_v$ assumption. The comoving $\rho_{\text{dust}}$ increased by at least one order of magnitude from $z \sim 7$ to $z \sim 1.5$. At lower redshift, a different data set leads to different trends, with $\rho_{\text{dust}}$ possibly remaining constant down to $z \sim 0$ or decreasing by a factor of approximately three. The N2GN sample leads to results systematically different from the recent ALMA $A^3$ COSMOS determination of $\rho_{\text{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 $\rho_{\text{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_{\text{dust}}$, and possibly in their underlying assumptions. A systematic and self-consistent study of $M_{\text{dust}}$ and $\rho_{\text{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 redshifts.
shift 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 \gg 7$ if dust already existed in galaxies in amounts similar to those in the N2GN sample. The natural next step in the exploration of $\rho_{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

**Acknowledgements.** 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 30m 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 et 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 - PEB-BLES). 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 Millimétrique (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 cos-

mic 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

Adam, R., Adane, A., Ade, P., et al. 2014, Proceedings of the 49th Rencontres de Morion on Cosmology: La Thuile, Italy, March 15-22, 2014, 73

Adam, R., Adane, A., Ade, P. A. R., et al. 2018, A&A, 609, A115

Adscheid, S., Magnelli, B., Liu, D., et al. 2024, A&A, 685, A1

Alaghband-Zadeh, S., Chapman, S. C., Swinbank, A. M., et al. 2013, MNRAS, 435, 1493

Alexander, D. M., Bauer, F. E., Brandt, W. N., et al. 2003, AJ, 126, 539

Andrews, S. K., Driver, S. P., Davies, L. J. M., et al. 2017, MNRAS, 464, 1569

Aoyama, S., Hou, K.-C., Hirashita, H., Nagamine, K., & Shimizu, I. 2018, MNRAS, 478, 4905

Aravena, M., Murphy, E. J., Aguirre, J. E., et al. 2013, MNRAS, 433, 498

Aravena, M., Hodge, J. A., Wagg, J., et al. 2014, MNRAS, 442, 558

Aravena, M., Spilker, J. S., Bethermin, M., et al. 2016, MNRAS, 457, 4406

Arrabal Haro, P., Rodríguez Espinosa, J. M., Muñoz-Tuñón, C., et al. 2018, MNRAS, 478, 3740

Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33

Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123

Astropy Collaboration (Price-Whelan, A. M., et al.) 2022, ApJ, 935, 167

Bakx, T. J. L. C., Dannerbauer, H., Frayer, D., et al. 2020, MNRAS, 496, 2372

Barger, A. J., Cowie, L. L., & Wang, W. H. 2008, ApJ, 689, 687

Barger, A. J., Wang, W. H., Cowie, L. L., et al. 2012, ApJ, 761, 89

Barger, A. J., Cowie, L. L., Blair, A. H., & Jones, L. H. 2022, ApJ, 934, 56

Barro, G., Pérez-González, P. G., Cava, A., et al. 2019, ApJS, 243, 22

Battisti, A. J., Cunha, E. D., Shivaei, I., & Calzetti, D. 2020, ApJ, 888, 108

Beeston, R. A., Wright, A. H., Maddox, S., et al. 2018, MNRAS, 479, 1077

Beeston, R. A., Gomez, H. L., Dunne, L., et al. 2024, MNRAS, 535, 3162

Berta, S., & Zylka, R. 2019-2024, Welcome to the PIIC, https://www.iram.fr/~gildas/dist/piic.pdf

Berta, S., Lonsdale, C. J., Polletta, M., et al. 2007, A&A, 476, 151

Berta, S., Lutz, D., Santini, P., et al. 2013a, A&A, 551, A100

Berta, S., Lutz, D., Nordon, R., et al. 2013b, A&A, 555, L8

Berta, S., Lutz, D., Genzel, R., Förster-Schreiber, N. M., & Tacconi, L. J. 2016, A&A, 587, A73

Berta, S., Stanley, F., Ismail, D., et al. 2023, A&A, 678, A28

Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393

Béthermin, M., De Breuck, C., Sargent, M., & Daddi, E. 2015, A&A, 576, L9

Béthermin, M., Wu, H.-Y., Lagache, G., et al. 2017, A&A, 607, A89

Bianchi, S. 2013, A&A, 552, A89

Bianchi, S., & Ferrara, A. 2005, MNRAS, 358, 379

Bing, L. J. 2022, Ph.D. Thesis, Mapping the Dusty Star Formation at High Redshift with the NIKA2 Cosmological Surveys, Aix-Marseille, France

Bing, L., Béthermin, M., Lagache, G., et al. 2023, A&A, 677, A66

Blain, A. W., Smail, I., Ivison, R. J., Kneib, J. P., & Frayer, D. T. 2002, Phys. Rep., 369, 111

Boquien, M., Burgarella, D., Roehlly, Y., et al. 2019, A&A, 622, A103

Borys, C., Scott, D., Chapman, S., et al. 2004, MNRAS, 355, 485

Bothwell, M. S., Chapman, S. C., Tacconi, L., et al. 2010, MNRAS, 405, 219

Bothwell, M. S., Smail, I., Chapman, S. C., et al. 2013, MNRAS, 429, 3047

Bothwell, M. S., Aguirre, J. E., Aravena, M., et al. 2017, MNRAS, 466, 2825

Bouché, N., Lehnert, M. D., Aguirre, A., Péroux, C., & Bergeron, J. 2007, MNRAS, 378, 525

Bourrion, O., Benoit, A., Bouly, J. L., et al. 2016, J. Instrum., 11, P11001

Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2015, ApJ, 803, 34

Brandt, W. N., Alexander, D. M., Hornschemeier, A. E., et al. 2001, AJ, 122, 2810

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

Burgarella, D., Buat, V., & Iglesias-Páramo, J. 2005, MNRAS, 360, 1413

Calvo, M., Benoît, A., Catalano, A., et al. 2016, J. Low Temp. Phys., 184, 816

Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682

Carilli, C. L., Daddi, E., Riechers, D., et al. 2010, ApJ, 714, 1407
Catalano, A., Calvo, M., Ponthieu, N., et al. 2014, A&A, 569, A9

Chabrier, G. 2003, ApJ, 586, L133

Chapin, E. L., Pope, A., Scott, D., et al. 2009, MNRAS, 398, 1793

Chapman, S. C., Blain, A. W., Smail, I., & Ivison, R. J. 2005, ApJ, 622, 772

Charlot, S., & Fall, S. M. 2000, ApJ, 539, 718

Chung, E. J., Rhee, M.-H., Kim, H., et al. 2009, ApJS, 184, 199

Ciesla, L., Charmandaris, V., Georgakakis, A., et al. 2015, A&A, 576, A10

Colina, L., Crespo Gómez, A., Álvarez-Márquez, J., et al. 2023, A&A, 673, L6

Combes, F., García-Burillo, S., Braine, J., et al. 2011, A&A, 528, A124

Combes, F., García-Burillo, S., Braine, J., et al. 2013, A&A, 550, A41

Cortzen, I., Magdis, G. E., Valentino, F., et al. 2020, A&A, 634, L14

Cowie, L. L., Barger, A. J., Hu, E. M., Capak, P., & Songaila, A. 2004, AJ, 127, 3137

Cowie, L. L., Barger, A. J., Hsu, L. Y., et al. 2017, ApJ, 837, 139

Cox, P., Neri, R., Berta, S., et al. 2023, A&A, 678, A26

Crespo Gómez, A., Colina, L., Álvarez-Márquez, J., et al. 2024, A&A, 691, A325

Cutri, R. M., Wright, E. L., Conrow, T., et al. 2013, Explanatory Supplement to the AllWISE Data Release Products, 1

da Cunha, E., Charlot, S., & Elbaz, D. 2008, MNRAS, 388, 1595

da Cunha, E., Walter, F., Decarli, R., et al. 2013, ApJ, 765, 9

da Cunha, E., Walter, F., Smail, I. R., et al. 2015, ApJ, 806, 110

Daddi, E., Dannerbauer, H., Stern, D., et al. 2009, ApJ, 694, 1517

Dale, D. A., Helou, G., Magdis, G. E., et al. 2014, ApJ, 784, 83

Dannerbauer, H., Harrington, K., Díaz-Sánchez, A., et al. 2019, AJ, 158, 34

Davé, R., Anglés-Alcázar, D., Narayanan, D., et al. 2019, MNRAS, 486, 2827

Davidzon, I., Ilbert, O., Laigle, C., et al. 2017, A&A, 605, A70

Davies, L. J. M., Driver, S. P., Robotham, A. S. G., et al. 2015, MNRAS, 447, 1014

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

Decarli, R., Walter, F., Gónzalez-López, J., et al. 2019, ApJ, 882, 138

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

Dickinson, M., & GOODS Legacy Team 2001, American Astronomical Society Meeting Abstracts, 198, 25.01

Dickinson, M., Bergeron, J., Casertano, S., et al. 2003a, Great Observatories Origins Deep Survey (GOODS), Validation Observations, Spitzer Proposal

Dickinson, M., Papovich, C., Ferguson, H. C., & Budavári, T. 2003b, ApJ, 587, 25

Dole, H., Lagache, G., Puget, J. L., et al. 2006, A&A, 451, 417

Draine, B. T. 2003, ARA&A, 41, 241

Draine, B. T., & Li, A. 2001, ApJ, 551, 807

Draine, B. T., & Li, A. 2007, ApJ, 657, 810

Draine, B. T., & Salpeter, E. E. 1979, ApJ, 231, 438

Draine, B. T., Dale, D. A., Bendo, G., et al. 2007, ApJ, 663, 866

Draine, B. T., Aniano, G., Krause, O., et al. 2014, ApJ, 780, 172

Driver, S. P., Popescu, C. C., Tuffs, R. J., et al. 2008, ApJ, 678, L101

Driver, S. P., Norberg, P., Baldry, I. K., et al. 2009, Astron. Geophys., 50, 5.12

Driver, S. P., Hill, D. T., Kelvin, L. S., et al. 2011, MNRAS, 413, 971

Driver, S. P., Andrews, S. K., da Cunha, E., et al. 2018, MNRAS, 475, 2891

Dudzevičiūtė, U., Smail, I., Swinbank, A. M., et al. 2020, MNRAS, 494, 3828

Dudzevičiūtė, U., Smail, I., Swinbank, A. M., et al. 2021, MNRAS, 500, 942

Dunne, L., & Eales, S. A. 2001, MNRAS, 327, 697

Dunne, L., Eales, S., Edmunds, M., et al. 2000, MNRAS, 315, 115

Dunne, L., Maddox, S. J., Ivison, R. J., et al. 2009, MNRAS, 394, 1307

Dunne, L., Gomez, H. L., da Cunha, E., et al. 2011, MNRAS, 417, 1510

Dunne, L., Bonavera, L., Gonzalez-Nuevo, J., Maddox, S. J., & Vlahakis, C. 2020, MNRAS, 498, 4635

Dunne, L., Maddox, S. J., Vlahakis, C., & Gomez, H. L. 2021, MNRAS, 501, 2573

Eales, S., & Ward, B. 2024, MNRAS, 529, 1130

Eales, S., Dunne, L., Clements, D., et al. 2010, PASP, 122, 499

Eisenstein, D. J., Johnson, B. D., Robertson, B., et al. 2023, ApJS, submitted [arXiv:2310.12340]

Elbaz, D., Dickinson, M., Hwang, H. S., et al. 2011, A&A, 533, A119

Enia, A., Talia, M., Pozzi, F., et al. 2022, ApJ, 927, 204

Evans, I. N., Evans, J. D., Martínez-Galarza, J. R., et al. 2024, ApJS, 274, 22

Feltre, A., Hatziminaoglou, E., Fritz, J., & Franceschini, A. 2012, MNRAS, 426, 120

Ferrarotti, A. S., & Gail, H. P. 2006, A&A, 447, 553

Finkelstein, S. L., Ryan, R. E., Papovich, C., et al. 2015, ApJ, 810, 71

Fixsen, D. J., Dwek, E., Mather, J. C., Bennett, C. L., & Shafer, R. A. 1998, ApJ, 508, 123

Fontana, A., Pozzetti, L., Donnarumma, I., et al. 2004, A&A, 424, 23

Fontana, A., Salimbeni, S., Grazian, A., et al. 2006, A&A, 459, 745

Frayer, D. T., Huynh, M. T., Chary, R., et al. 2006, ApJ, 647, L9

Freundlich, J., Combes, F., Tacconi, L. J., et al. 2019, A&A, 622, A105

Fritz, J., Franceschini, A., & Hatziminaoglou, E. 2006, MNRAS, 366, 767

Fujimoto, S., Ouchi, M., Shibuya, T., & Nagai, H. 2017, ApJ, 850, 83

Fujimoto, S., Brammer, G. B., Watson, D., et al. 2022, Nature, 604, 261

Fukugita, M. 2011, ArXiv e-prints [arXiv:1103.4191]

Geach, J. E., Smail, I., Moran, S. M., et al. 2011, ApJ, 730, L19

Gehrz, R. 1989, IAU Symp., 135, 445

Genzel, R., Baker, A. J., Tacconi, L. J., et al. 2003, ApJ, 584, 633

Genzel, R., Tacconi, L. J., Lutz, D., et al. 2015, ApJ, 800, 20

George, R. D., Ivison, R. J., Hopwood, R., et al. 2013, MNRAS, 436, L99

Gioannini, L., Matteucci, F., & Calura, F. 2017, MNRAS, 471, 4615

Granato, G. L., Ragone-Figueroa, C., Taverna, A., et al. 2021, MNRAS, 503, 511

Greve, T. R., Pope, A., Scott, D., et al. 2008, MNRAS, 389, 1489

Griffin, M. J., Abergel, A., Abreu, A., et al. 2010, A&A, 518, L3

Guidetti, D., Bondi, M., Prandoni, I., et al. 2017, MNRAS, 471, 210

Hagimoto, M., Bakx, T. J. L. C., Serjeant, S., et al. 2023, MNRAS, 521, 5508

Harris, A. I., Baker, A. J., Zonak, S. G., et al. 2010, ApJ, 723, 1139

Harris, A. I., Baker, A. J., Frayer, D. T., et al. 2012, ApJ, 752, 152

Henriques, B. M. B., Yates, R. M., Fu, J., et al. 2020, MNRAS, 491, 5795

Hezaveh, Y. D., Marrone, D. P., Fassnacht, C. D., et al. 2013, ApJ, 767, 132

Hollenbach, D., & Salpeter, E. E. 1971, ApJ, 163, 155

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

Inoue, A. K., Shimizu, I., Iwata, I., & Tanaka, M. 2014, MNRAS, 442, 1805

Ismail, D., Beelen, A., Buat, V., et al. 2023, A&A, 678, A27

Ivison, R. J., Greve, T. R., Smail, I., et al. 2005, MNRAS, 337, 1

Ivison, R. J., Greve, T. R., Dunlop, J. S., et al. 2007, MNRAS, 380, 199

Ivison, R. J., Swinbank, A. M., Swinyard, B., et al. 2010, A&A, 518, L35

Ivison, R. J., Papadopoulos, P. P., Smail, I., et al. 2011, MNRAS, 412, 1913

Ivison, R. J., Swinbank, A. M., Smail, I., et al. 2013, ApJ, 772, 137

Jin, S., Daddi, E., Magdis, G. E., et al. 2022, A&A, 665, A3

Jones, A. P., Tielens, A. G. G. M., Hollenbach, D. J., & McKee, C. F. 1994, ApJ, 433, 797

Kocevski, D. D., Finkelstein, S. L., Barro, G., et al. 2024, ApJ, accepted [arXiv:2404.03576]

Kodra, D., Andrews, B. H., Newman, J. A., et al. 2023, ApJ, 942, 36

Labbe, I., Greene, J. E., Bezanson, R., et al. 2025, ApJ, 978, 92

Lagache, G., Puget, J.-L., & Dole, H. 2005, ARA&A, 43, 727

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

Lewis, J. S. W., Ocvirk, P., Dubois, Y., et al. 2023, MNRAS, 519, 5987

Li, A., & Draine, B. T. 2001, ApJ, 554, 778

Li, Q., Narayanan, D., & Davé, R. 2019, MNRAS, 490, 1425

Lim, C.-F., Wang, W.-H., Smail, I., et al. 2020, ApJ, 889, 80

Liske, J., Baldry, I. K., Driver, S. P., et al. 2015, MNRAS, 452, 2087

Liu, D., Daddi, E., Dickinson, M., et al. 2018, ApJ, 853, 172

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

Lutz, D. 2014, ARA&A, 52, 373

Lutz, D., Poglitsch, A., Altieri, B., et al. 2011, A&A, 532, A90

Maudu, P. 1995, ApJ, 441, 18

Magdis, G. E., Daddi, E., Béthermin, M., et al. 2012, ApJ, 760, 6

Magnelli, B., Lutz, D., Santini, P., et al. 2012, A&A, 539, A155

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

Magnelli, B., Ivison, R. J., Lutz, D., et al. 2015, A&A, 573, A45

Magnelli, B., Boogaard, L., Decarli, R., et al. 2020, ApJ, 892, 66

Mancini, C., Matute, I., Cimatti, A., et al. 2009, A&A, 500, 705

Maraston, C. 2005, MNRAS, 362, 799

Maseda, M. V., de Graaff, A., Franx, M., et al. 2024, A&A, 689, A73

Mathis, J. S. 1990, ARA&A, 28, 37

Matthee, J., Naidu, R. P., Brammer, G., et al. 2024, ApJ, 963, 129

Mckee, C. F., Hollenbach, D. J., Seab, G. C., & Tielens, A. G. G. M. 1987, ApJ, 318, 674

Ménard, B., & Fukugita, M. 2012, ApJ, 754, 116

Millard, J. S., Eales, S. A., Smith, M. W. L., et al. 2020, MNRAS, 494, 293

Momcheva, I. G., Brammer, G. B., van Dokkum, P. G., et al. 2016, ApJS, 225, 27

Monfardini, A., Adam, R., Adane, A., et al. 2014, J. Low Temp. Phys., 176, 787

Nanni, A., Bressan, A., Marigo, P., & Girardi, L. 2013, MNRAS, 434, 2390

Neri, R., Downes, D., Cox, P., & Walter, F. 2014, A&A, 562, A35

Noll, S., Burgarella, D., Giovannoli, E., et al. 2009, A&A, 507, 1793

Nozawa, T., Kozasa, T., Habe, A., et al. 2007, ApJ, 666, 955

Ocvirk, P., Gillet, N., Shapiro, P. R., et al. 2016, MNRAS, 463, 1462

Oesch, P. A., Brammer, G., Naidu, R. P., et al. 2023, MNRAS, 525, 2864

Oliver, S. J., Bock, J., Altieri, B., et al. 2012, MNRAS, 424, 1614

Ostriker, J., & Silk, J. 1973, ApJ, 184, L113

Owen, F. N. 2018, ApJS, 235, 34

Parente, M., Ragone-Figueroa, C., Granato, G. L., & Lapi, A. 2023, MNRAS, 521, 6105
Pearson, C. P., Serjeant, S., Negrello, M., et al. 2010, A&A, 514, A9

Penner, K., Pope, A., Chapin, E. L., et al. 2011, MNRAS, 410, 2749

Penney, J. I., Blain, A. W., Assef, R. J., et al. 2020, MNRAS, 496, 1565

Perera, T. A., Chapin, E. L., Austermann, J. E., et al. 2008, MNRAS, 391, 1227

Pérez-González, P. G., Rieke, G. H., Egami, E., et al. 2005, ApJ, 630, 82

Pérez-González, P. G., Rieke, G. H., Villar, V., et al. 2008, ApJ, 675, 234

Pérez-González, P. G., Barro, G., Rieke, G. H., et al. 2024, ApJ, 968, 4

Perotto, L., Ponthieu, N., Macías-Pérez, J. F., et al. 2020, A&A, 637, A71

Péroux, C., & Howk, J. C. 2020, ARA&A, 58, 363

Pilbratt, G. L., Riedinger, J. R., Passvogel, T., et al. 2010, A&A, 518, L1

Planck Collaboration VI. 2020, A&A, 641, A6

Poglitsch, A., Waelkens, C., Geis, N., et al. 2010, A&A, 518, L2

Pope, A., Borys, C., Scott, D., et al. 2005, MNRAS, 358, 149

Pope, A., Scott, D., Dickinson, M., et al. 2006, MNRAS, 370, 1185

Pope, A., Chary, R.-R., Alexander, D. M., et al. 2008, ApJ, 675, 1171

Popesso, P., Concas, A., Cresci, G., et al. 2023, MNRAS, 519, 1526

Popping, G., Somerville, R. S., & Galametz, M. 2017, MNRAS, 471, 3152

Pozzi, F., Calura, F., Zamorani, G., et al. 2020, MNRAS, 491, 5073

Pozzi, F., Calura, F., Fudamoto, Y., et al. 2021, A&A, 653, A84

Puget, J. L., Abergel, A., Bernard, J. P., et al. 1996, A&A, 308, L5

Reuter, C., Vieira, J. D., Spilker, J. S., et al. 2020, ApJ, 902, 78

Rho, J., Kozasa, T., Reach, W. T., et al. 2008, ApJ, 673, 271

Richards, E. A., Kellermann, K. I., Fomalont, E. B., Windhorst, R. A., & Partridge, R. B. 1998, AJ, 116, 1039

Riechers, D. A., Hodge, J., Walter, F., Carilli, C. L., & Bertoldi, F. 2011, ApJ, 739, L31

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

Rovilos, E., Georgantopoulos, I., Akylas, A., & Fotopoulou, S. 2010, A&A, 522, A11

Rudnick, G., Hodge, J., Walter, F., et al. 2017, ApJ, 849, 27

Saintonge, A., Lutz, D., Genzel, R., et al. 2013, ApJ, 778, 2

Sandage, A., Tammann, G. A., & Yahil, A. 1979, ApJ, 232, 352

Santini, P., Maiolino, R., Magnelli, B., et al. 2014, A&A, 562, A30

Sargent, B. A., Srinivasan, S., Meixner, M., et al. 2010, ApJ, 716, 878

Schechter, P. 1976, ApJ, 203, 297

Schmidt, M. 1968, ApJ, 151, 393

Schneider, R., Valiante, R., Ventura, P., et al. 2014, MNRAS, 442, 1440

Scoville, N., Aussel, H., Brusa, M., et al. 2007, ApJS, 172, 1

Sharon, C. E., Riechers, D. A., Hodge, J., et al. 2016, ApJ, 827, 18

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

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

Skelton, R. E., Whitaker, K. E., Momcheva, I. G., et al. 2014, ApJS, 214, 24

Smail, I., Ivison, R. J., Owen, F. N., Blain, A. W., & Kneib, J. P. 2000, ApJ, 528, 612

Solomon, P. M., Downes, D., Radford, S. J. E., & Barrett, J. W. 1997, ApJ, 478, 144

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

Stach, S. M., Dudzevičiūtė, U., Smail, I., et al. 2019, MNRAS, 487, 4648

Staguhn, J. G., Kovács, A., Arendt, R. G., et al. 2014, ApJ, 790, 77

Stalevski, M., Fritz, J., Baes, M., Nakos, T., & Popović, L. Č. 2012, MNRAS, 420, 2756

Stalevski, M., Ricci, C., Ueda, Y., et al. 2016, MNRAS, 458, 2288

Steidel, C. C., Adelberger, K. L., Shapley, A. E., et al. 2003, ApJ, 592, 728

Sun, F., Helton, J. M., Egami, E., et al. 2024, ApJ, 961, 29

Swinbank, A. M., Smail, I., Chapman, S. C., et al. 2004, ApJ, 617, 64

Tacconi, L. J., Neri, R., Genzel, R., et al. 2013, ApJ, 768, 74

Tacconi, L. J., Genzel, R., Saintonge, A., et al. 2018, ApJ, 853, 179

Tacconi, L. J., Genzel, R., & Sternberg, A. 2020, ARA&A, 58, 157

Teplitz, H. I., Chary, R., Elbaz, D., et al. 2011, AJ, 141, 1

Thomson, A. P., Ivison, R. J., Smail, I., et al. 2012, MNRAS, 425, 2203

Tielens, A. G. G. M. 1998, ApJ, 499, 267

Todini, P., & Ferrara, A. 2001, MNRAS, 325, 726

Traina, A., Magnelli, B., Gruppioni, C., et al. 2024a, A&A, 690, A84

Traina, A., Gruppioni, C., Delvecchio, I., et al. 2024b, A&A, 681, A118

Übler, H., D'Eugenio, F., Perna, M., et al. 2024, MNRAS, 533, 4287

Valentino, F., Magdis, G. E., Daddi, E., et al. 2018, ApJ, 869, 27

Villanueva, V., Ibar, E., Hughes, T. M., et al. 2017, MNRAS, 470, 3775

Vlahakis, C., Dunne, L., & Eales, S. 2005, MNRAS, 364, 1253

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

Wang, W. H., Cowie, L. L., & Barger, A. J. 2004, ApJ, 613, 655

Wang, S., Liu, J., Qiu, Y., et al. 2016a, ApJS, 224, 40

Wang, T., Elbaz, D., Schreiber, C., et al. 2016b, ApJ, 816, 84

Wang, W.-H., Lin, W.-C., Lim, C.-F., et al. 2017, ApJ, 850, 37

Wang, T., Elbaz, D., Daddi, E., et al. 2018, ApJ, 867, L29

Weibel, A., Oesch, P. A., Barrufet, L., et al. 2024, MNRAS, 533, 1808

Weingartner, J. C., & Draine, B. T. 2001, ApJ, 548, 296

Williams, R. E., Blacker, B., Dickinson, M., et al. 1996, AJ, 112, 1335

Wolfire, M. G., Hollenbach, D., McKee, C. F., Tielens, A. G. G. M., & Bakes, E. L. O. 1995, ApJ, 443, 152

Xiao, M., Oesch, P., Elbaz, D., et al. 2023, ArXiv e-prints [arXiv:2309.02492]

Yang, C., Omont, A., Beelen, A., et al. 2017, A&A, 608, A144

Yates, R. M., Hendriks, D., Vijayan, A. P., et al. 2024, MNRAS, 527, 6292

Zhukovska, S., Gail, H. P., & Trieloff, M. 2008, A&A, 479, 453

Zylka, R. 2013, Astrophysics Source Code Library [record ascl:1303.011]

¹ Institut de Radioastronomie Millimétrique (IRAM), 300 Rue de la Piscine, 38400 Saint-Martin-d'Hères, France

² Aix Marseille Univ., CNRS, CNES, LAM (Laboratoire d'Astrophysique de Marseille), Marseille, France

³ Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, France

⁴ School of Physics and Astronomy, Cardiff University, Queen's Buildings, The Parade, Cardiff CF24 3AA, UK

⁵ Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, 91191 Gif-sur-Yvette, France

⁶ Institut de Radioastronomie Millimétrique (IRAM), Avenida Divina Pastora 7, Local 20, E-18012 Granada, Spain

⁷ Institut Néel, CNRS, Université Grenoble Alpes, 25 Av. des Martyrs, 38000 Grenoble, France

⁸ Université de Strasbourg, CNRS, Observatoire Astronomique de Strasbourg, UMR 7550, 67000 Strasbourg, France

⁹ Astronomy Centre, Department of Physics and Astronomy, University of Sussex, Brighton BN1 9QH, UK

¹⁰ Univ. Grenoble Alpes, CNRS, Grenoble INP, LPSC-IN2P3, 53, Avenue des Martyrs, 38000 Grenoble, France

¹¹ Dipartimento di Fisica, Sapienza Università di Roma, Piazzale Aldo Moro 5, I-00185 Roma, Italy

¹² Université Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France

¹³ Institute for Research in Fundamental Sciences (IPM), School of Astronomy, Tehran, Iran

¹⁴ Physics Department "Ettore Pancini", Università degli Studi di Napoli "Federico II", Via Cintia 21, I-80126 Napoli, Italy

¹⁵ Centro de Astrobiología (CSIC-INTA), Torrejón de Ardoz 28850, Madrid, Spain

¹⁶ National Observatory of Athens, Institute for Astronomy, Astrophysics, Space Applications and Remote Sensing, Ioannou Metaxa and Vasileos Pavlou, GR-15236 Athens, Greece

¹⁷ Department of Astrophysics, Astronomy & Mechanics, Faculty of Physics, University of Athens, Panepistimiopolis, GR-15784 Zografos, Athens, Greece

¹⁸ High Energy Physics Division, Argonne National Laboratory, 9700 South Cass Avenue, Lemont, IL 60439, USA

¹⁹ LERMA, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Université, UPMC, 75014 Paris, France

²⁰ Institute of Space Sciences (ICE), CSIC, Campus UAB, Carrer de Can Magrans s/n, E-08193 Barcelona, Spain

²¹ ICREA, Pg. Lluís Companys 23, Barcelona, Spain

²² Joint ALMA Observatory, Alonso de Córdova 3107, Vitacura 763-0355, Santiago, Chile

²³ European Southern Observatory, Alonso de Córdova 3107, Vitacura, Casilla 19001, Santiago de Chile, Chile

²⁴ IRAP, CNRS, Université de Toulouse, CNES, UT3-UPS, Toulouse, France

²⁵ School of Physics and Astronomy, University of Leeds, Leeds LS2 9JT, UK

²⁶ Dipartimento di Fisica, Università di Roma "Tor Vergata", Via della Ricerca Scientifica 1, I-00133 Roma, Italy

²⁷ Laboratoire de Physique de l'École Normale Supérieure, ENS, PSL Research University, CNRS, Sorbonne Université, Université de Paris, 75005 Paris, France

²⁸ Institut d'Astrophysique de Paris, CNRS (UMR7095), 98 Bis boulevard Arago, 75014 Paris, France

²⁹ University of Lyon, UCB Lyon 1, CNRS/IN2P3, IP21, 69622 Villeurbanne, France

³⁰ Department of Astronomy, University of Geneva, Chemin Pegasi 51, 1290 Versoix, Switzerland
# 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_{\text{dust}}$, $M_{\star}$ and $L_{\text{IR}}(8-1000\,\mu\text{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_{\text{dust}}$ or $L_{\text{IR}}(8-1000\,\mu\text{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 $\kappa_v$ normalization by Draine et al. (2014). We remind that the first DL07 modeling came with the parametrization of the dust absorption coefficient $\kappa_v$ by Li & Draine (2001), later revised by Draine (2003) ($\kappa_{850\,\mu\text{m}} = 0.038\,\text{m}^2\,\text{kg}^{-1}$). Draine et al. (2014) re-normalized their assumption on dust properties such that $\kappa_{850\,\mu\text{m}} = 0.047\,\text{m}^2\,\text{kg}^{-1}$. On the other hand, da Cunha et al. (2008, 2015) adopted the normalization by Dunne et al. (2000), that is $\kappa_{850\,\mu\text{m}} = 0.077\,\text{m}^2\,\text{kg}^{-1}$.

Figure A.1 compares the $M_{\text{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_{\text{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\times$ 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_{\text{dust}}$, to be at least partially ascribed to the missing warmer dust component that dominates the $\lambda_{\text{rest}} < 50\,\mu\text{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_{\star}$-SFR plane. Appendix B investigates the effect on depletion timescales. The most striking cases are very red galaxies at high redshift, with $M_{\star}$ 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(\text{MAGP}) - A_V(\text{CIG}) > 1.3$. These differences testify how critical and uncertain can the $M_{\star}$ 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 $\lambda_{\text{rest}} > 50\,\mu\text{m}$. Consequently, the integrated luminosity of these models is $L_{\text{FIR}}(50-1000\,\mu\text{m})$, rather than $L_{\text{IR}}$. By integrating the templates family of Herschel star forming galaxies defined by Berta et al. (2013a), we recovered a typical ratio $L_{\text{FIR}}/L_{\text{IR}} \approx 0.7$ for this kind of galaxies (Berta et al. 2023). We therefore rescaled the MBB $L_{\text{FIR}}$ estimate to the galaxies’ $L_{\text{IR}}$ using this ratio.

Figure A.3 compares the $L_{\text{IR}}$ results. In the cases when an AGN-torus was necessary in the fit, its contribution to the $L_{\text{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.
<figure><img src="image_17.png" /><figcaption>**Fig. A.1.** Comparison of $M_{\text{dust}}$ 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).</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. A.2.** Comparison of $M_*$ obtained with the five different UV-to-radio SED fitting methods considered here.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Fig. A.3.** Comparison of $L_{IR}$ obtained with the eight different SED fitting methods considered here.</figcaption></figure>

## Appendix B: Effect of $M_\star$ on the $\tau_{\text{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_\star$, redshift and distance from the MS to $\tau_{\text{dep}}$. This approach takes as reference the MS parametrization by Speagle et al. (2014), shown in Fig. 10. The $\tau_{\text{dep}}$ thus obtained depends on the value of $M_\star$. 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_\star$ 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_\star$ obtained with MAGPHYS/SED3FIT or with CIGALE. Few sources migrate to lower values of $\tau_{\text{dep}}$, and from the locus of the MS to the region occupied by starbursts. This happens specifically for the group of galaxies at $z \sim 5$, that is going to be the subject of Lagache et al. (in prep.).

<figure><img src="image_20.png" /><figcaption>**Fig. B.1.** Consequence of changing $M_\star$ when computing $\tau_{\text{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_\star$. The references of the literature data and lines are listed in the caption of Fig. 9</figcaption></figure>
## Appendix C: The importance of $\kappa_v$ 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 compared to literature data, pointing out the need to reconcile all data to the same dust model normalization, namely to the same $\kappa_v$. Our choice of reference is $\kappa_{850\mu m} = 0.047 \text{ m}^2 \text{ kg}^{-1}$ (Draine et al. 2014).

If a different dust model was adopted, a re-normalization of dust masses is needed, because $M_{\text{dust}}$ depends inversely on $\kappa_v$. 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 $\kappa_v$ 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) $\kappa_{850} = 0.038 \text{ m}^2 \text{ kg}^{-1}$ value, in order to compare the literature dust masses to N2GN's, they would need to be rescaled by a factor $\kappa_{850}(D03)/\kappa_{850}(D14)$, because dust mass is inversely proportional to $\kappa_v$. In other words, if the adopted $\kappa_{850}$ 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 $\kappa_{850}$ 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 $\kappa_v$ 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 $\kappa_{850} = 0.077 \text{ m}^2 \text{ kg}^{-1}$. 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 $\kappa_{250} = 0.89 \text{ m}^2 \text{ kg}^{-1}$, which they state corresponds to $\kappa_{850} = 0.077 \text{ m}^2 \text{ kg}^{-1}$.

Pozzi et al. (2020, 2021) used the value $\kappa_{1.2\text{THz}} = 4.0 \text{ cm}^2 \text{ g}^{-1}$, corresponding to Draine (2003, $\kappa_{850} = 0.038 \text{ m}^2 \text{ kg}^{-1}$). Magnelli et al. (2020) adopted $\kappa_{850} = 0.043 \text{ m}^2 \text{ kg}^{-1}$. 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_{\text{dust}}$ assumptions of 25 and 35 K. In all three cases, they adopted the Draine et al. (2014) $\kappa_v$ 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 $\kappa_v$ different from the one originally adopted is a shift of the DMF in logarithmic space to larger or smaller masses. The net result on $\rho_{\text{dust}}$ (the integral of the DMF) is therefore a simple scaling by the same amount as for $M_{\text{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 $\kappa_v$ assumptions. In the left panel the data are shown “as they are” computed by each work, each with its specific value of $\kappa_v$. 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.

<table><caption>**Table C.1.** Summary of the $\kappa_v$ values adopted in the literature.</caption><thead><tr><th>Adopted κ<sub>v</sub><br>[m<sup>2</sup> kg<sup>-1</sup>]</th><th>κ<sub>v</sub> reference</th><th>Work</th></tr></thead><tbody><tr><td>MW κ<sub>850</sub> = 0.047</td><td>Draine et al. (2014)</td><td>This work<br>Péroux & Howk (2020)<br>Traina et al. (2024a)</td></tr><tr><td>MW κ<sub>850</sub> = 0.038</td><td>Draine (2003)</td><td>Pozzi et al. (2020)<br>Pozzi et al. (2021)</td></tr><tr><td>MW κ<sub>850</sub> = 0.077</td><td>Dunne et al. (2000)</td><td>MAGPHYS<br>Beeston et al. (2018)<br>Beeston et al. (2024)<br>Driver et al. (2018)<br>Dudzevičiūtė et al. (2021)<br>Dunne et al. (2011)<br>Eales & Ward (2024)<br>Vlahakis et al. (2005)</td></tr><tr><td>MW κ<sub>850</sub> = 0.043</td><td>Li & Draine (2001)</td><td>Magnelli et al. (2020)</td></tr><tr><td>SMC κ<sub>V</sub> = 1.54 × 10<sup>3</sup></td><td>Weingartner & Draine (2001)</td><td>Ménard & Fukugita (2012)</td></tr></tbody></table>
<figure><img src="image_21.png" /><figcaption>**Fig. C.1.** Effect of rescaling $\rho_{\text{dust}}$ to the same $\kappa_{\nu}$. Left: data without rescaling; each data set is shown as obtained with its $\kappa_{\nu}$ underlying assumption. Right: the same data after rescaling to $\kappa_{850} = 0.047 \, \text{m}^2 \, \text{kg}^{-1}$ (Draine et al. 2014).</figcaption></figure>

## Appendix D: $\rho_{\text{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.

<figure><img src="image_22.png" /><figcaption>**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.</figcaption></figure>
## 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 $^{12}$CO spectral line energy distribution (SLED) of this source and its dynamical properties, by combining NOEMA and VLA observations of $^{12}$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 \times 10^{10} M_{\odot}$. 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 \times 10^{-5}$ µ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 \pm 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\alpha$ 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 \pm 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\alpha$ 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 $^{12}$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 \sim 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 $\mu \sim 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 \pm 0.16$.

N2GN_1_09 is the ID12646 in the Spitzer-driven superdeblended 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 \times 10^{-5}$ µ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*<sub>spec</sub> = 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*<sub>spec</sub> = 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*<sub>phot</sub> = 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*<sub>spec</sub> = 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*<sub>phot</sub> = 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*<sub>3.6</sub>*S*<sub>Ks</sub> ~ 13, interpreted as a deep Balmer break combined to a large extinction, that leads to a *z*<sub>phot</sub> = 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*<sub>spec</sub> = 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*<sub>spec</sub> = 2.098.

**N2GN_1_23** is an optically dark galaxy detected by JWST/FRESCO at *z*<sub>spec</sub> = 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*<sub>phot</sub> = 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*<sub>phot</sub> = 2.697 and 2.817 (Kodra et al. 2023).

**N2GN_1_25** was observed with NOEMA and identified as an optical galaxy at *z*<sub>spec</sub> = 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*<sub>phot</sub> = 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*<sub>spec</sub> = 1.988, 1.989; Kodra et al. 2023). The system lies at the posi-
tion 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_{\text{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_{\text{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_{\text{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 \times 10^{-4}$ µ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_{\text{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_{\text{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_{\text{phot}} = 4.12 \pm 0.53$ taking into account the available photometry (Sect. 3.3).

**N2GN_1_35** is identified with VLA with an edge-on $z_{\text{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_{\text{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_{\text{phot}} = 3.93 \pm 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_{\text{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_{\text{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_{\text{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_{\text{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_{\text{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_{\text{phot}} = 2.228$ optically bright galaxy showing a distinct break at $\sim 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 $\sim 10\sigma$ in the F277W band, $\sim 3\sigma$ 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_{\text{phot}} = 5.87 \pm 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\alpha$ and CIII]λ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} \sim 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 $\sim 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 \times 10^{-5}$ µ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 ($\sim 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\alpha$ 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 an 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 $\lambda \sim 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 defined 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 superdeblended 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).
<table><caption>**Table E.1.** Details of the N2GN sources: identification through high spatial resolution millimeter or radio data; redshift (spectroscopic or photometric); results of SED fitting.</caption><thead><tr><th rowspan="2">N2CLS ID</th><th rowspan="2">Name IAU (NIKA2)</th><th rowspan="2">Alt. †</th><th rowspan="2">Method</th><th colspan="3">Identification</th><th rowspan="2">z</th><th rowspan="2">Redshift Nature</th><th rowspan="2">Ref.</th><th colspan="3">SED fitting results</th><th rowspan="2">log(L<sub>☉</sub>)</th><th rowspan="2">log(L<sub>☉</sub>)</th></tr><tr><th>Ref.</th><th>R.A.</th><th>[hh:mm:ss]</th><th>log(M<sub>☉</sub>)</th><th>2.5%</th><th>97.5%</th><th>L<sub>IR</sub> [10<sup>2</sup>L<sub>☉</sub>]</th><th>2.5%</th><th>97.5%</th></tr></thead><tbody><tr><td>N2GN_1_01</td><td>N2GN J123633+621408</td><td>GN10</td><td>SMA</td><td>1</td><td>12:36:33.44</td><td>62:14:08.7</td><td>5.303</td><td>spec</td><td>a</td><td>20.92</td><td>9.18</td><td>9.36</td><td>23.90</td><td>13.22</td><td>13.46</td></tr><tr><td>N2GN_1_02</td><td>N2GN J123730+621259</td><td>850-6</td><td>SMA</td><td>1</td><td>12:37:30.79</td><td>62:12:58.9</td><td>1.963</td><td>phot</td><td>b</td><td>60.46</td><td>9.78</td><td>9.78</td><td>4.73</td><td>12.68</td><td>12.68</td></tr><tr><td>N2GN_1_03</td><td>N2GN J123707+621408</td><td>GN12</td><td>SMA</td><td>1</td><td>12:37:07.21</td><td>62:14:08.3</td><td>2.484</td><td>spec</td><td>c</td><td>38.26</td><td>9.58</td><td>9.59</td><td>4.71</td><td>12.58</td><td>12.68</td></tr><tr><td>N2GN_1_04</td><td>N2GN J123711+622211</td><td>GN20</td><td>PdBI</td><td>2</td><td>12:37:11.90</td><td>62:22:12.1</td><td>4.055</td><td>spec</td><td>d</td><td>66.28</td><td>9.71</td><td>9.91</td><td>19.85</td><td>13.29</td><td>13.30</td></tr><tr><td>N2GN_1_05</td><td>N2GN J123711+621330</td><td>GN1200.3</td><td>NOEMA</td><td>3</td><td>12:37:11.33</td><td>62:13:30.9</td><td>1.996</td><td>spec</td><td>e</td><td>21.48</td><td>9.31</td><td>9.59</td><td>6.20</td><td>12.76</td><td>12.82</td></tr><tr><td>N2GN_1_06<sup>‡</sup></td><td>N2GN J123652+621226</td><td>HDF850.1</td><td>PdBI</td><td>4</td><td>12:36:51.99</td><td>62:12:25.7</td><td>5.185</td><td>spec</td><td>f</td><td>11.52</td><td>9.03</td><td>9.06</td><td>8.27</td><td>12.92</td><td>12.96</td></tr><tr><td>N2GN_1_07</td><td>N2GN J123645+621448</td><td>GN12</td><td>SMA</td><td>1</td><td>12:36:46.09</td><td>62:14:48.5</td><td>2.960</td><td>spec</td><td>g</td><td>21.07</td><td>8.99</td><td>9.41</td><td>7.23</td><td>12.85</td><td>13.00</td></tr><tr><td>N2GN_1_08</td><td>N2GN J123631+621714</td><td>AzGN2</td><td>SMA</td><td>1</td><td>12:36:31.92</td><td>62:17:14.7</td><td>4.340</td><td>spec</td><td>h</td><td>24.05</td><td>9.10</td><td>9.38</td><td>12.45</td><td>12.75</td><td>13.24</td></tr><tr><td>N2GN_1_09</td><td>N2GN J123627+621217</td><td>AzGN10</td><td>PdBI</td><td>5.6</td><td>12:36:27.50</td><td>62:12:17.8</td><td>4.147</td><td>spec</td><td>i</td><td>12.49</td><td>8.93</td><td>9.13</td><td>3.86</td><td>12.45</td><td>12.76</td></tr><tr><td>N2GN_1_10</td><td>N2GN J123635+621423</td><td></td><td>VLA</td><td>7</td><td>12:36:35.59</td><td>62:14:24.0</td><td>2.005</td><td>spec</td><td>j</td><td>11.21</td><td>9.01</td><td>9.11</td><td>6.25</td><td>12.73</td><td>12.80</td></tr><tr><td>N2GN_1_11</td><td>N2GN J123713+621826</td><td>GN40</td><td>VLA</td><td>7</td><td>12:37:13.89</td><td>62:18:26.3</td><td>3.990</td><td>phot</td><td>k</td><td>8.49</td><td>8.90</td><td>9.11</td><td>12.90</td><td>13.03</td><td>13.16</td></tr><tr><td>N2GN_1_12<sup>a</sup></td><td>N2GN J123636+621155</td><td>GN11</td><td>PdBI</td><td>6</td><td>12:36:37.50</td><td>62:11:56.5</td><td>2.590</td><td>phot</td><td>h</td><td>33.69</td><td>8.39</td><td>9.11</td><td>5.05</td><td>12.27</td><td>13.89</td></tr><tr><td>N2GN_1_13</td><td>N2GN J123636+621155</td><td>GN11</td><td>PdBI</td><td>6</td><td>12:36:36.10</td><td>62:11:54.4</td><td>4.163</td><td>phot</td><td>b</td><td>1.98</td><td>8.15</td><td>9.99</td><td>12.36</td><td>13.04</td><td>13.14</td></tr><tr><td>N2GN_1_14</td><td>N2GN J123701+621146</td><td>GN17</td><td>NOEMA</td><td>3</td><td>12:36:58.61</td><td>62:14:51.1</td><td>5.181</td><td>spec</td><td>l</td><td>6.34</td><td>8.80</td><td>8.88</td><td>18.13</td><td>12.71</td><td>13.26</td></tr><tr><td>N2GN_1_15</td><td>N2GN J123618+621550</td><td>GN6</td><td>SMA</td><td>1</td><td>12:37:01.60</td><td>62:11:45.9</td><td>1.760</td><td>spec</td><td>b</td><td>19.16</td><td>9.26</td><td>9.44</td><td>3.69</td><td>12.53</td><td>12.57</td></tr><tr><td>N2GN_1_16</td><td>N2GN J123622+621615</td><td>GN9</td><td>SMA</td><td>1</td><td>12:36:18.32</td><td>62:15:50.7</td><td>2.000</td><td>spec</td><td>b</td><td>17.89</td><td>9.25</td><td>9.25</td><td>5.77</td><td>12.77</td><td>12.77</td></tr><tr><td>N2GN_1_17<sup>a</sup></td><td>N2GN J123738+621734</td><td>GN37</td><td>SMA</td><td>1</td><td>12:36:22.08</td><td>62:16:16.2</td><td>3.470</td><td>phot</td><td>h</td><td>12.11</td><td>8.95</td><td>9.31</td><td>17.65</td><td>12.67</td><td>13.42</td></tr><tr><td>N2GN_1_17<sup>b</sup></td><td>N2GN J123738+621734</td><td>GN37</td><td>NOEMA</td><td>3</td><td>12:37:38.11</td><td>62:17:37.1</td><td>3.190</td><td>spec</td><td>m</td><td>13.49</td><td>8.90</td><td>9.26</td><td>2.39</td><td>12.36</td><td>12.49</td></tr><tr><td>N2GN_1_18</td><td>N2GN J123702+621425</td><td></td><td>NOEMA</td><td>3</td><td>12:37:02.59</td><td>62:14:26.7</td><td>4.390</td><td>phot</td><td>h</td><td>—</td><td>8.20</td><td>9.26</td><td>—</td><td>—</td><td>—</td></tr><tr><td>N2GN_1_19</td><td>N2GN J123616+621515</td><td>GN4</td><td>SMA</td><td>1</td><td>12:36:16.10</td><td>62:15:13.7</td><td>2.578</td><td>spec</td><td>b</td><td>6.07</td><td>8.62</td><td>9.03</td><td>7.78</td><td>12.53</td><td>13.14</td></tr><tr><td>N2GN_1_20</td><td>N2GN J123550+621041</td><td>AzGN4</td><td>SMA</td><td>1</td><td>12:35:50.22</td><td>62:10:42.4</td><td>5.333</td><td>phot</td><td>n</td><td>10.66</td><td>8.92</td><td>9.26</td><td>30.19</td><td>12.89</td><td>13.48</td></tr><tr><td>N2GN_1_21</td><td>N2GN J123634+621241</td><td>GN26</td><td>VLA</td><td>7</td><td>12:36:34.52</td><td>62:12:41.0</td><td>1.219</td><td>spec</td><td>b</td><td>8.84</td><td>8.94</td><td>8.94</td><td>3.61</td><td>12.56</td><td>12.56</td></tr><tr><td>N2GN_1_22</td><td>N2GN J123553+621340</td><td>AzGN13</td><td>SMA</td><td>1</td><td>12:35:53.23</td><td>62:13:37.9</td><td>2.098</td><td>spec</td><td>b</td><td>33.07</td><td>9.30</td><td>9.56</td><td>2.59</td><td>12.31</td><td>12.44</td></tr><tr><td>N2GN_1_23</td><td>N2GN J123656+621207</td><td>GN15</td><td>NOEMA</td><td>3</td><td>12:36:56.54</td><td>62:12:07.5</td><td>5.179</td><td>spec</td><td>o</td><td>3.84</td><td>8.52</td><td>8.76</td><td>10.42</td><td>12.87</td><td>13.31</td></tr><tr><td>N2GN_1_24<sup>a</sup></td><td>N2GN J123719+621219</td><td>AzGN27</td><td>NOEMA</td><td>3</td><td>12:37:18.94</td><td>62:12:17.5</td><td>2.697</td><td>phot</td><td>b</td><td>2.38</td><td>7.67</td><td>9.30</td><td>1.54</td><td>11.95</td><td>12.31</td></tr><tr><td>N2GN_1_24<sup>b</sup></td><td>N2GN J123719+621219</td><td>AzGN27</td><td>NOEMA</td><td>3</td><td>12:37:19.61</td><td>62:12:21.0</td><td>2.817</td><td>phot</td><td>b</td><td>8.18</td><td>7.73</td><td>9.68</td><td>1.23</td><td>11.68</td><td>12.66</td></tr><tr><td>N2GN_1_25</td><td>N2GN J123712+621212</td><td>GN1200.29</td><td>NOEMA</td><td>3</td><td>12:37:12.05</td><td>62:12:11.8</td><td>2.914</td><td>spec</td><td>b</td><td>22.95</td><td>8.96</td><td>9.36</td><td>2.44</td><td>12.37</td><td>12.55</td></tr><tr><td>N2GN_1_26</td><td>N2GN J123740+621221</td><td>GN38</td><td>SMA</td><td>1</td><td>12:37:41.14</td><td>62:12:20.4</td><td>5.310</td><td>phot</td><td>b</td><td>13.20</td><td>8.68</td><td>9.33</td><td>17.26</td><td>13.23</td><td>13.25</td></tr><tr><td>N2GN_1_27<sup>a</sup></td><td>N2GN J123621+621709</td><td>GN7</td><td>NOEMA</td><td>3</td><td>12:36:21.26</td><td>62:17:08.5</td><td>1.988</td><td>spec</td><td>b</td><td>6.85</td><td>8.76</td><td>9.03</td><td>3.73</td><td>12.49</td><td>15.58</td></tr><tr><td>N2GN_1_27<sup>b</sup></td><td>N2GN J123621+621709</td><td>GN7</td><td>NOEMA</td><td>3</td><td>12:36:20.98</td><td>62:17:09.7</td><td>1.989</td><td>spec</td><td>b</td><td>16.70</td><td>8.77</td><td>9.25</td><td>6.30</td><td>12.79</td><td>12.80</td></tr><tr><td>N2GN_1_28</td><td>N2GN J123728+621920</td><td></td><td>VLA</td><td>7</td><td>12:37:28.12</td><td>62:19:20.2</td><td>3.222</td><td>spec</td><td>g</td><td>9.80</td><td>8.86</td><td>9.13</td><td>5.30</td><td>12.67</td><td>12.79</td></tr><tr><td>N2GN_1_29</td><td>N2GN J123702+621302</td><td>GN18</td><td>VLA</td><td>7</td><td>12:37:02.57</td><td>62:13:02.4</td><td>2.819</td><td>phot</td><td>b</td><td>6.97</td><td>8.60</td><td>8.86</td><td>3.77</td><td>12.51</td><td>12.62</td></tr><tr><td>N2GN_1_30</td><td>N2GN J123723+621714</td><td></td><td>NOEMA</td><td>3</td><td>12:37:23.52</td><td>62:17:15.0</td><td>2.613</td><td>phot</td><td>b</td><td>3.80</td><td>8.57</td><td>8.94</td><td>4.89</td><td>12.53</td><td>12.71</td></tr><tr><td>N2GN_1_32</td><td>N2GN J123712+621035</td><td>AzGN20</td><td>VLA</td><td>7</td><td>12:36:29.45</td><td>62:15:13.2</td><td>3.652</td><td>spec</td><td>b</td><td>25.38</td><td>9.09</td><td>9.45</td><td>3.89</td><td>12.33</td><td>12.71</td></tr><tr><td>N2GN_1_33</td><td>N2GN J123628+621044</td><td>GN25</td><td>JWST</td><td>8</td><td>12:36:29.03</td><td>62:10:45.5</td><td>5.388</td><td>phot</td><td>b</td><td>4.46</td><td>8.60</td><td>8.77</td><td>19.66</td><td>13.26</td><td>13.47</td></tr><tr><td>N2GN_1_34<sup>a</sup></td><td>N2GN J123644+621938</td><td>AzGN28</td><td>NOEMA</td><td>3</td><td>12:36:45.67</td><td>62:19:39.9</td><td>—</td><td>—</td><td>—</td><td>—</td><td>8.57</td><td>9.47</td><td>—</td><td>—</td><td>—</td></tr><tr><td>N2GN_1_34<sup>b</sup></td><td>N2GN J123644+621938</td><td>AzGN28</td><td>NOEMA</td><td>3</td><td>12:36:44.04</td><td>62:19:38.5</td><td>4.120</td><td>phot</td><td>h</td><td>15.08</td><td>8.79</td><td>9.57</td><td>19.67</td><td>12.91</td><td>13.37</td></tr><tr><td>N2GN_1_35</td><td>N2GN J123637+620851</td><td></td><td>VLA</td><td>7</td><td>12:36:37.04</td><td>62:08:52.4</td><td>1.950</td><td>phot</td><td>b</td><td>15.62</td><td>9.19</td><td>9.19</td><td>2.13</td><td>12.32</td><td>12.32</td></tr><tr><td>N2GN_1_36</td><td>N2GN J123658+620930</td><td></td><td>SMA</td><td>1</td><td>12:36:58.53</td><td>62:09:31.6</td><td>3.930</td><td>phot</td><td>h</td><td>5.85</td><td>8.66</td><td>9.03</td><td>8.25</td><td>12.90</td><td>12.94</td></tr></tbody></table>