© The Authors 2026

1. Introduction

At high redshifts (z ≳ 4), and despite recent progress, the contribution of dust-obscured galaxies to the young Universe remains quite uncertain (e.g. Traina et al. 2024). Recent Atacama Large Millimetre/sub-millimetre Array (ALMA) surveys (such as ALPINE and REBELS) show that dust-obscured (far-IR and millimetre) and unobscured (UV rest-frame) star formation could have made a similar contribution to the cosmic budget even at z ∼ 6 (e.g. Gruppioni et al. 2020; Fudamoto et al. 2020; Khusanova et al. 2021; Barrufet et al. 2023; Fujimoto et al. 2024). In addition, recent JWST observations have reported candidate heavily obscured sources at z ≳ 7 (Akins et al. 2023; Rodighiero et al. 2023; Gandolfi et al. 2025; Martis et al. 2025), although spectroscopic confirmation and direct detection of the thermal dust continuum are still lacking. These studies are not based on volume-complete surveys and/or rely on complex corrections and selections. The population of luminous and dusty galaxies at z > 4 is largely unconstrained, and accurately quantifying the total star formation rate density (SFRD) at z ≳ 4 from a well-controlled volume-limited survey remains an ambitious goal. One of the main contributors to the obscured SFRD at such high redshifts are millimetre-bright, ultra-massive, and optically dark galaxies¹. A recent JWST study based on three millimetre-bright sources (Xiao et al. 2024) at z ∼ 5 − 6 shows that they account for ∼20% of the SFRD derived from rest-frame far-UV to optical observations and corrected for dust attenuation at z ∼ 5.8. These three galaxies are also ultra-massive (log(M_★/M_⊙) > 11.0) and would require, on average, about 50% of the baryons in their halos to be converted into stars. This is two to three times higher than even the most efficient galaxies at later times (Xiao et al. 2024, and references therein).

The two ultra-massive galaxies in the great observatories origins deep surveys north (GOODS-N) field (S3 and GN10) are located in a large-scale structure in the process of formation (Herard-Demanche et al. 2025), which also contains the well known HDF850.1 galaxy. HDF850.1 is forming in one of the richest environments identified to date at z > 5, with ∼100 z = 5.17 − 5.20 galaxies (Herard-Demanche et al. 2025; Sun et al. 2024; see also e.g. Calvi et al. 2021). Even before the extent of this structure was known, the redshifts of HDF850.1 and GN10 alone already suggested a significant overdensity of dusty star-forming galaxies at z > 5 in this field compared to expectations, though it remained unclear whether the individual substructures traced by these sources were physically connected (e.g. Riechers et al. 2020).

Observations suggest that such (massive) overdensities are not rare at z ∼ 5 (e.g. Guaita et al. 2022; Helton et al. 2024a). They are likely progenitors of today’s massive galaxy clusters. Understanding their distribution and properties shortly after the epoch of reionisation helps bridge observations of the early Universe with modern-day structures. Given the existence of a similar galaxy overdensity in both the GOODS-N and GOODS-S fields, Sun et al. (2024) suggest that 50% ± 20% of the cosmic star formation at z = 5.1 − 5.5 occurs in protocluster environments.

Galaxies in high-redshift overdensities often have different properties compared to those in less dense regions. For example, Lemaux et al. (2022) find a highly significant positive relation between the star formation rate (SFR) and galaxy densities at 3 ≲ z ≲ 5, with galaxies in the densest protocluster environments exhibiting average SFRs > 0.5 dex higher than those in the coeval field. At higher redshifts (z > 4.5), Li et al. (2025) observe that galaxies in dense environments exhibit a slightly elevated SFR at a given mass compared to those in less dense regions. Additionally, they find that galaxies in high-density regions tend to have redder UV slopes, which may indicate increased dust extinction (see also Helton et al. 2024b). Herard-Demanche et al. (2025) compared the properties of galaxies forming within the GOODS-N z ∼ 5.2 overdense region to those outside and found that the masses, SFRs, and UV luminosities of galaxies inside the overdensity are significantly higher than those outside. Morishita et al. (2025) identify two overdensities of galaxies at z ∼ 5.7 along the sightline of the galaxy cluster Abell 2744, where some members exhibit properties suggestive of accelerated maturation. These studies show that the environment influences galaxy growth early in cosmic history even if high-z galaxies are also much more gas-rich than local ones, for the same total mass, which could also give rise to very high SFRs.

As part of the NIKA2 Cosmological Legacy Survey (N2CLS), we conducted very deep observations of GOODS-N at both 1.2 and 2 mm (Bing et al. 2023; Berta et al. 2025; Béthermin et al. 2025). We covered 160 arcmin² and to a depth close to the confusion noise of the Institut de Radioastronomie Millimétrique (IRAM) 30 m telescope (Ponthieu et al. 2025). The large-scale overdense region at z = 5.1 − 5.3, as observed by the FRESCO and JADES JWST programmes (Herard-Demanche et al. 2025; Sun et al. 2024), is completely covered by N2CLS. The survey provides a unique opportunity to study in detail the obscured star formation in the over-dense environment and analyse the properties of its millimetre-bright dusty star-forming galaxy (DSFG). This paper presents the discovery of N2GN_1_13 at z_spec = 5.182, a massive DSFG belonging to the overdensity that we characterised using Northern Extended Millimetre Array (NOEMA) observations. Furthermore, by combining our analysis with JWST/FRESCO data, we revealed another N2CLS DSFG at z_spec = 5.201. We additionally identify three N2CLS galaxies with a photometric redshift of z ∼ 5.3, one of which is optically dark.

The paper is organised as follows. Section 2 gives a summary of the N2CLS survey of GOODS-N. In Sect. 3 we present the NOEMA observations of N2GN_1_13 and give its properties, which we then compare with those of HDF850.1 and GN10 in Sect. 4. In Sect. 5 (complemented by Appendices B and C) we present the remaining five N2CLS DSFGs identified at z = 5.1 − 5.3 (two with z_spec, three with z_phot) and discuss their physical properties. In Sect. 6 we discuss the efficiency of the DSFGs in converting baryonic gas into stars. In Sect. 7 we report their contribution to the cosmic SFRD and present their localisation in the overdense structure. In Sect. 8 we compare our results with those expected from empirical, semi-analytical, and hydrodynamical models. We finally provide a summary in Sect. 9. We also examine the dust production in these high-redshift galaxies in Appendix A, and present another N2CLS galaxy candidate at z ∼ 5.2 (N2GN_1_20) in Appendix D.

Throughout the paper, we adopt a Λ cold dark matter cosmology with H₀ = 67.4 km s⁻¹ Mpc⁻¹, Ω_m = 0.315, and Ω_Λ = 0.685 (Planck Collaboration VI 2020), a Chabrier (2003) stellar initial mass function (IMF), and the Draine (2003) frequency-dependent dust absorption coefficient, κ_ν, re-normalised as indicated by Draine et al. (2014, i.e. κ_850 μm = 0.047 m² kg⁻¹, as reference).

2. DSFGs from the NIKA2 cosmological legacy survey in GOODS-N

N2CLS is the deepest and largest survey ever made at 1.2 and 2 mm, covering COSMOS and GOODS-N (Table 1 in Bing et al. 2023). In GOODS-N, it reaches an average 1σ noise level of 0.17 and 0.048 mJy over 160 sq. arcmin. Our complete sample with a purity threshold > 95% comprises 65 galaxies detected with S/N > 4.2 in a least one of the NIKA2 bands. Thanks to a very successful NOEMA continuum observation programme, combined with sub-millimetre array (SMA) 860 μm (Cowie et al. 2017) and deep radio observations (Owen 2018), all GOODS-N sources have accurate positions that allow us to search for their multi-wavelength counterparts (Berta et al. 2025, hereafter B25). The final sample (including multiplicities) comprises 71 individual millimetre sources. Three galaxies (two identified with NOEMA, one identified with the VLA) have no multi-wavelength counterpart and thus no redshift. We refer the reader to B25 for a complete description of the identification process and construction of multi-wavelength spectral energy distributions (SEDs). The paper also provides a full description of the data employed for the SED fitting, which have been released and are publicly accessible². We have spectroscopic redshifts for 29 sources³ (i.e. 41% of the sample) and photometric redshifts for 39 sources.

Among the 71 individual galaxies, 16 are classified as Hubble Space Telescope (HST)-dark (some of which are also IRAC-dark). Two HST-dark galaxies remain unidentified despite their NOEMA (N2GN_1_17_b) and Very Large Array (VLA; N2GN_1_55) counterpart. For the other 14 galaxies (6 with spectroscopic redshifts and 8 with photometric redshifts), the average redshift is z = 4.9 and the average stellar mass is M_★ = 3.9 × 10¹¹ M_⊙ (stellar masses are derived using CIGALE; see Sect. 3.2). This is in agreement with the general consensus in the literature that the majority of these optically dark galaxies are heavily dust-obscured, star-forming sources at z ∼ 3 − 6 that dominate the high-mass end of the galaxy population at these cosmic times (e.g. Alcalde Pampliega et al. 2019; Wang et al. 2019; Barrufet et al. 2023; Xiao et al. 2023; Gottumukkala et al. 2024).

3. N2GN_1_13: Identifying a new DSFG in the overdense environment of HDF850.1

We selected one of these optically dark N2CLS galaxies, N2GN_1_13, observed with JWST/NIRCam as part of the FRESCO programme (Oesch et al. 2023) but without a confirmed redshift, as a secure target for NOEMA follow-up observations. N2GN_1_13 is a bright DSFG, with NIKA2 fluxes of 1.95 ± 0.25 and 0.49 ± 0.08 mJy at 1.2 and 2 mm, respectively. The accurate position of N2GN_1_13 has been secured thanks to our previous NOEMA programme (W21CV, P.I. L. Bing), which detected at 2 mm a single counterpart with a flux of 393 ± 47 μJy and a beam size of 1.61″ × 0.93″. As illustrated in Fig. 1 (see also Fig. C.1), the galaxy is absent in the HST images and emerges in JWST images at λ ≳ 3.3 μm. The galaxy has a 6.2σ single line detection in the NIRCam grism spectrum (Xiao, priv. comm.) but no conclusive redshift measurement. Assuming the NIRCam line is the Hα line, then z ∼ 5.196 and N2GN_1_13 would have the same redshift as the gigantic overdensity (Herard-Demanche et al. 2025; Sun et al. 2024).

Fig. 1. N2GN_1_13. Left: False colour image produced using HST data, with F814W, F125W, and F160W images in blue, green, and red, respectively. Overlaid in red are the NOEMA continuum contours at 2 mm, ranging from 4σ to 14σ in steps of 2. The white contours are JWST F444W brightness. The scale bar in the bottom left has a length of 1″. Middle: NOEMA CO(4−3) and CO(5−4) line spectra. The red line shows the Gaussian fit to the line profile (including the continuum). Notice that the CO(5−4) line lies precisely at the interface between two basebands, resulting in the partial fragmentation of the spectral profile. Right: NOEMA 3.5 mm continuum map. Contours are displayed with a step of 12 μJy/beam (corresponding to 2σ).

3.1. NOEMA observations

We observed N2GN_1_13 with NOEMA in band-1 using the C-configuration as part of programme W23CX (PI: G. Lagache). The total on-source integration time was 6.6 hours. The receivers were tuned to cover a frequency range suitable for detecting the CO(4−3) and CO(5−4) transitions at z ∼ 5.196.

The NOEMA data were processed using the standard calibration and reduction procedures implemented in the GILDAS⁴ software package. This resulted in calibrated uv tables for the lower and upper sidebands, sampled with a spectral resolution of 20 MHz. The source 1125+596 was used as the phase and amplitude calibrator. The absolute flux uncertainty is estimated to be within 10%, while the positional accuracy is better than 0.16″. The final dataset achieves a synthesised beam of 2.3″ × 2.1″, with a line sensitivity of 0.7 mJy/MHz and a continuum sensitivity of 5.7 μJy/beam.

We detect two spectral lines with S/N > 6, which we unambiguously identify as CO(4−3) and CO(5−4) at z = 5.182. As expected from the CO/[CI] flux line ratio in distant galaxies and our sensitivity, the [CI](1−0) at 79.6 GHz remains undetected with S​_[CI] < 0.89 Jy km s⁻¹ (assuming the same full width at half maximum as CO(4−3))⁵. Figure 1 presents the two CO detections, fitted along with the continuum. The best-fit Gaussian parameters derived from the line profiles are listed in Table 1.

To characterise the continuum emission, we combined all channels from the lower and upper sidebands, excluding the spectral windows containing line emission, and produced a double sideband continuum map (Fig. 1). The map was cleaned using natural weighting with a limited number of clean components. The continuum flux was measured using two methods: integration within the 2σ contour and direct measurement at the position of the emission peak. Both approaches yield a flux of 57 ± 6 μJy at 84.6 GHz, suggesting that the continuum source is unresolved.

3.2. Properties of N2GN_1_13

The total IR luminosity (L_IR; integrated between 8 and 1000 μm), SFR, stellar mass (M_★), dust mass (M_dust), and dust spectral index (β) have been derived in B25 using multi-wavelength SED fitting. The dust mass has been obtained using Draine (2003) frequency-dependent dust absorption coefficient (κ_ν), re-normalised as indicated by Draine et al. (2014), i.e. κ_850 μm = 0.047 m² kg⁻¹ as reference. The panchromatic optical to radio SEDs are modelled 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. 2013), as well as with CIGALE (Burgarella et al. 2005; Noll et al. 2009; Boquien et al. 2019). The spectral index is obtained using a modified black-body fit to the far-IR to millimetre data. The effect of the cosmic microwave background is taken into account as described by da Cunha et al. (2013).

The stellar mass has been obtained using the same assumptions for CIGALE and the high-z versions of MAGPHYS and SED3FIT, i.e. Bruzual & Charlot (2003, BC03) stellar populations with a Chabrier (2003) IMF and the Charlot & Fall (2000) two-components attenuation law. B25 presents the comparison of stellar masses obtained using the different codes, and shows a very good agreement, but for few outliers lying at high-z. Unfortunately, the sources considered in this paper are among the outliers. For example, considering the 6 optically dark sources at z ∼ 5.2 (see Table 3), the ratio between the high-z versions of MAGPHYS or SED3FIT and CIGALE stellar mass determination ranges from 1.4 to 9.2, with an average ratio of 5.8. Masses up to 10¹² M_⊙ are derived using MAGPHYS and SED3FIT (and CIGALE for N2GN_1_20), indicating the presence of extremely massive objects at this redshift. Excluding any active galactic nucleus (AGN) contamination, the high flux observed in the long-wavelength NIRCam/JWST F444W filter (∼0.7−1.3 μJy) supports by itself the very large stellar masses, i.e. M_★ ∼ 10¹¹ − 10¹² M_⊙. Such very high stellar masses go with high extinction. For N2GN_1_13, the stellar masses from SED3FIT and CIGALE are log(M_★) = 11.68 and log(M_★) = 11.03, respectively. In the following, and to be conservative, we adopt the lower stellar mass values obtained using CIGALE. To be consistent, all the other parameters (e.g. SFR, dust mass) are taken from the CIGALE SED fitting.

The total molecular gas mass is obtained using $M_{\mathrm{Mol}}={\alpha }_{\mathrm{CO}}\times L{\prime }_{\mathrm{CO}(1-0)}$. Following Berta et al. (2023, see also Dunne et al. 2022), we adopted a conversion factor α_CO = 4.0 M_⊙ (K km s⁻¹ pc²)⁻¹, including the contribution of helium. We transformed the luminosity of the CO lines detected by NOEMA into $L{\prime }_{\mathrm{CO}(1-0)}$ by adopting the average ratios given by Carilli & Walter (2013) for sub-millimetre galaxies (SMGs), r_43/10 = 0.46, and r_54/10 = 0.39, which give r_54/43 (and also r_54/21; see Sect. 4) consistent with our measurement. We obtain M_Mol = (2.3 ± 0.6) × 10¹⁰ M_⊙ and (3.1 ± 0.9) × 10¹⁰ M_⊙ from the CO(5−4) and CO(4−3) line luminosity, respectively. We adopted the molecular gas mass from the CO(4−3), which is detected with a higher S/N and with a measurable full width at half maximum (Table 1).

We could then derive the depletion time of the molecular gas using the molecular gas mass and the SFR:

$$\begin{array}{c}\hfill {\tau }_{\mathrm{dep}}=\frac{M_{\mathrm{mol}}}{\mathrm{SFR}}·\end{array}$$

This is the time it would take for a galaxy to consume the entirety of its molecular gas reservoir, if it formed stars at the rate derived from far-IR and sub-millimetre data in a single event (not considering any mass return to the interstellar medium (ISM) from stellar winds and supernovae). We also finally computed the gas-to-dust mass ratio (GDR) as M_Mol/M_dust.

We summarise the physical properties of N2GN_1_13 in Table 2 and show JWST/NIRCam images in Fig. C.1. Based on these findings, we can conclude that:

• N2GN_1_13 is a very massive galaxy (1.07 × 10¹¹ M_⊙). Using the same survey volume as FRESCO (Xiao et al. 2024), such a high mass requires a high fraction of the baryons (∼30%) to be converted into stars (ϵ_★; as discussed in Sect. 6). This is significantly higher than what is commonly observed in lower-z galaxies (Xiao et al. 2024; Gentile et al. 2024).

• N2GN_1_13 has a very red continuum with F182M − F444W = 3.8 ABmag. As shown in B25, the optical photometry reveals a rest-frame UV excess, which is reproduced by adding an unextinguished young (ages between 10 and 100 Myr) simple stellar population to the SED modelling. As can be seen in Fig. C.1, it has an extended morphology. The source is very red simply because it is dusty with (i) a very high A_V = 3.96 ± 0.13 mag, indicating a large dust content and (ii) rest-frame dust emission detected by NOEMA at both 324 and 574 μm. Furthermore, this source cannot be classified as a broad-line AGN, as the Hα line does not show any additional broadening beyond what is expected from morphological and kinematic effects. It has no X-ray counterpart and no rest-frame mid-IR data are available to verify the possible presence of an AGN-torus component contributing to its mid-IR SED.

• The very high SFR relative to the available gas reservoir results in a short depletion time of 43 Myr, indicative of a brief (≤100 Myr) starburst phase⁶. Such a low τ_dep points to an extreme case of star formation. Assuming no gas accretion or replenishment, and modelling the gas consumption as an exponential decay M_gas(t) = M_gas, 0 × e^−t/τ_dep, the galaxy’s SFR would drop to 10% of its initial value within approximately 100 Myr, effectively marking the quenching of star formation.

• With a dust mass of (7.7 ± 0.8) × 10⁸ M_⊙, and given its stellar mass, N2GN_1_13 lies near the maximal theoretical dust yield expected from the combined contributions of supernovae (assuming no reverse-shock destruction) and asymptotic giant branch (AGB) stars (see Fig. A.1 and Witstok et al. 2023). However, the short gas depletion time suggests that dust production from AGB stars, whose contribution operates on longer timescales, is likely minimal. This implies that additional dust formation mechanisms, such as rapid grain growth in the ISM, may be required to account for the observed dust mass.

• Our observed gas-to-dust ratio leads to a galaxy metallicity consistent with solar, Z ≈ Z_⊙ (e.g. Li et al. 2019). However, this estimate should be interpreted with caution, as our adopted values for α_CO and κ_ν are themselves based on the assumption of solar-like metallicity.

• With a CO(5−4)/CO(4−3) flux ratio of 1.37, N2GN_1_13 closely matches the ratio observed in the galaxy CGG-z4.b at z = 4.3 (Brinch et al. 2025), and falls within the upper range of values reported in the literature (see Fig. 5 in Brinch et al. 2025). This elevated ratio suggests the presence of high gas temperatures, densities, and/or pressures.

• N2GN_1_13 has a very high SFR and is thus in the process of very efficient stellar mass build-up. Compared to the star-forming main sequence (SFMS), N2GN_1_13 lies significantly above the relation if we adopt the parametrisation of Popesso et al. (2023), with ΔSFR = log(SFR/SFR_MS) = 0.55 ± 0.02 (for a scatter of 0.09 dex). Using the parametrisation of Speagle et al. (2014), however, the offset is smaller, ΔSFR = 0.29 ± 0.04 (for a scatter of 0.2 dex). These SFMS relations were calibrated with pre-JWST data, and it is now well established that massive galaxies at z ∼ 5 can remain undetected in HST observations. As a result, the SFMS at the high-mass end (M_★ ≳ 10¹¹ M_⊙) is not necessarily well constrained. Considering the specific star-formation rate (sSFR), N2GN_1_13 has log(sSFR) = −8.18 ± 0.13, which is close to the empirical limit of log(sSFR) ≲ −8.05 for main-sequence (MS) galaxies reported by Caputi et al. (2017, see also Rinaldi et al. 2025), based on star-forming galaxies selected from optical and near-IR surveys. Taken together, these results suggest that N2GN_1_13 is broadly consistent with, though slightly above, the SFMS, within the current uncertainties at the high-mass end.

N2GN_1_13 is thus a massive, extended, and highly obscured galaxy caught in a short-lived starburst phase at z ∼ 5.2. Its high stellar mass, red optical continuum, and gas reservoir point to rapid evolution and imminent quenching.

4. Gas masses, gas-to-dust ratios, and gas depletion times for the three brightest DSFGs in the overdense environment

In addition to N2GN_1_13, two other galaxies in the overdense environment have CO measurements: GN10 and HDF850.1 (see the compilation in Riechers et al. 2020, R20). We used the recipes described in Sect. 3 to compute gas masses, GDR and gas depletion times for these two galaxies and use the other source properties (e.g. SFR) as derived from the SED fitting with CIGALE. The results are given in Table 2.

We used the CO(2−1) line measurements for HDF850.1 and GN10 and the conversion factor from Carilli & Walter (2013), which gives M_Mol consistent with R20 if we consider α_CO = 1 rather than α_CO = 4, and assuming a mean gravitational lensing magnification factor ${\overline{\mu }}_L=1.6$ for HDF850.1 as in R20. If we further correct for κ_850 μm (which is equal to 0.057 m² kg⁻¹ in R20 versus 0.047 here), we would obtain GDR = 46 and 31 w.r.t. 65 and 130 in R20, and τ_dep = 39 and 34 Myr w.r.t. to 70 and 40 Myr in R20, for GN10 and HDF850.1, respectively. The large difference in the GDR for HDF850.1 arises from a factor of 3.7 difference in M_dust (assuming the same κ_850 μm), while the discrepancy in τ_dep for GN10 is mainly due to a factor of 1.6 difference in the SFR. Notice that we have the same spectral index β as in previous studies (R20, Walter et al. 2012).

We derive a gas mass fraction (f_gas = M_Mol/(M_Mol + M_★)) of ∼0.6 for HDF850.1 and GN10, which is typical of gas fraction of z ∼ 5 SMGs (Dudzevičiūtė et al. 2020). For N2GN_1_13, the gas mass fraction is ∼3 times lower but remains consistent with certain z ∼ 5 SMGs reported in the literature (e.g. Zavala et al. 2022). However, we warn that such measurements have to be taken with caution. For instance, Sun et al. (2024) reported a gas fraction of 0.18 for HDF850.1, which is significantly lower than the value obtained in this study. This discrepancy primarily arises from the considerably lower gas mass reported in Neri et al. (2014), attributed to the adoption of a much lower α_CO (this is also the case in Zavala et al. 2022 where α_CO = 1.0 M_⊙ (K km s⁻¹ pc²)⁻¹).

Tacconi et al. (2020) give updated scaling relations that links M_★, redshift and distance from the MS of star-forming galaxies (from Speagle et al. 2014) to molecular gas depletion timescale (τ_dep) and thus molecular gas mass (M_mol). We used these relations to compute τ_dep and M_mol (see Table 4) and obtain τ_dep = 206, 142, 161 Myr and M_mol = 1.4 × 10¹¹, 5.1 × 10¹⁰, 2.6 × 10¹¹ M_⊙ for N2GN_1_13, HDF850.1 (de-lensed) and GN10, respectively. While this compares quite well to the results given in Table 2 for HDF850.1 and GN10 (∼5% level), M_mol (and thus τ_dep) is overestimated by a factor of ∼4.5 for N2GN_1_13. If we used higher J transitions to get the molecular mass for HDF850.1 and GN10, we would have obtained the same discrepancy (i.e. molecular gas mass of 3.1 and 5.4 × 10¹⁰ M_⊙ using J_up = 5 rather than 12 and 25 × 10¹⁰ M_⊙). This shows how these determinations have to be taken with caution.

5. Bright dusty galaxies in the overdense environment

5.1. N2CLS galaxies at similar redshifts

In addition to N2GN_1_01 (GN10), N2GN_1_06 (HDF850.1), and N2GN_1_13, we have five more galaxies with redshifts and position compatible with the overdense structure⁷ (see Table 3). Two are optically dark N2CLS galaxies with spectroscopic redshift from JWST: N2GN_1_23 (also known as S3 in Xiao et al. 2024) and N2GN_1_61. For the three additional galaxies with previously reported photometric redshifts of z ∼ 5.2 (see B25 and references therein), independent photometric redshifts were also derived using the enhanced JWST photometric coverage (when available) and the SED-fitting code New-Hyperz. This is a modified version of the public code Hyperz originally developed by Bolzonella et al. (2000)⁸. The procedure and its validation is discussed in Appendix B. Results are shown in Table 3. Our new photometric redshifts are within one-sigma uncertainty of the literature based z_phot, for all but N2GN_1_26 (2σ). The new photometric redshifts therefore support their possible association with the z ∼ 5.1 − 5.4 overdensity. While photometric redshifts provide the best constraints currently available, spectroscopic confirmation will be required to secure their membership. For example, a photometric-redshift-selected galaxy, N2GN_1_33, initially estimated at z_phot ∼ 5.3 (B25) was later confirmed with JWST spectroscopy to lie at z = 1.78, illustrating the uncertainties inherent in photometric redshift estimates. In the following, we therefore clearly distinguish results for galaxies with spectroscopic and photometric redshifts. Among the three N2CLS galaxies with z_phot, one is optically dark (N1_1_43), while the remaining two are optically bright (N2GN_1_26 and N2GN_1_57). Details on each of the five sources are provided in Appendix C.

Of the eight sources (five with z_spec and three with z_phot), N2GN_1_01, 06, and 23 belong to the overdensity according to Sun et al. (2024), Herard-Demanche et al. (2025), and Xiao et al. (2024); N2GN_1_13 and N2GN_1_61 fall into their area of study but were not identified. N2GN_1_43 is in the JADES footprint but was not reported. N2GN_1_57 lies at the border of the JADES footprint and N2GN_1_26 is farther away.

5.2. Source properties

The multi-wavelength SEDs for these eight sources are shown in B25. Here, we used CIGALE to derive their properties. As in B25, we adopted the Bruzual & Charlot (2003) stellar models, with a Chabrier (2003) IMF and solar metallicity, the Charlot & Fall (2000) attenuation law, and the Draine & Li (2007) dust emission models, re-normalised as in Draine et al. (2014). A delayed and truncated star formation history 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.

To assess the impact of photometric-redshift uncertainties on the derived physical parameters, we repeated the fits at z_phot ± 1σ for the three sources with photometric redshifts. The additional uncertainties were estimated from the standard deviation of the resulting physical parameters and added in quadrature to the fitting errors returned by CIGALE at z_phot.

We summarise the results in Table 4. Our sample is characterised by high V-band attenuation (A_V), SFR and stellar mass, with median values of A_V = 3.26, SFR = 587 M_⊙ yr⁻¹ and log(M_★/M_⊙) = 10.93. These values are comparable to those reported for the three ultra-massive galaxies in Xiao et al. (2024). However, they are systematically much higher than those found in the larger DSFG sample from JWST (36 galaxies) analysed in the same study.

Combining the derived SFR and M_★, and the known redshift, we applied the Tacconi et al. (2020) scaling relations to derive the molecular gas depletion timescale, τ_dep, and mass, M_gas = τ_dep × SFR of these N2GN sources. Figure 2 locates the sources in the M_★ versus SFR plane and compares their position to the MS loci defined by Speagle et al. (2014) and Popesso et al. (2023). The Tacconi et al. (2020) scalings are based on the distance of the sources from the Speagle et al. (2014) MS, which does not include the bending at large M_★, appearing clearly instead in Popesso et al. (2023).

The N2GN z ∼ 5.2 galaxies lie significantly above the SFMS as defined by Popesso et al. (2023), indicating a bursty mode of star formation. Although they appear closer to the SFMS when adopting the Speagle et al. (2014) parametrisation, the eight galaxies still lie above the 1σ envelope of the relation. This result will need to be confirmed for such high-mass and high-SFR galaxies at such redshifts once the SFMS is better constrained with JWST data.

The median τ_dep obtained for these galaxies is of 0.18 Gyr, with individual values spanning over the range between ∼0.12 and ∼0.23 Gyr. These short timescales (or equivalently a very high star formation efficiency, ϵ = 1/τ_dep ∼ 5.6 × 10⁻⁹ yr⁻¹) confirm the intense ongoing starbursts activity that powers their far-IR and millimetre emission. Despite their short depletion timescales, the starbursts are still operating over periods significantly longer than their dynamical timescales, suggesting that the activity is sustained through multiple bursts or extended starburst phases. The fact that these galaxies reside in a confirmed overdense environment likely provides both the fuel and favourable conditions to maintain such prolonged star formation.

Given their short gas depletion timescales, these DSFGs may transition into a quiescent phase by z ∼ 4.5, assuming a constant star formation efficiency, complete molecular gas consumption, and no external gas replenishment. This evolutionary pathway would be consistent with the emergence of massive quiescent galaxies observed in other overdense environments at z > 4 (e.g. Carnall et al. 2024; Kakimoto et al. 2024; Tanaka et al. 2024).

We finally examined the rest-frame optical morphologies of our sources using JWST/NIRCam imaging (Fig. C.1). Most galaxies appear extended or clumpy, with no strong evidence of point-like or compact morphologies; only one galaxy in the sample appears marginally compact. Several galaxies show irregular or asymmetric light distributions, suggestive of ongoing or recent dynamical activity, such as minor mergers, interactions, or internal disk instabilities. These features are broadly consistent with recent JWST results for high-redshift, dust-obscured galaxies (e.g. McKinney et al. 2025), which reveal a wide range of morphologies, from clumpy and irregular to disk-like systems, among massive DSFGs.

5.3. Impact of missing far-IR data on SFR and stellar mass estimates for our sample

We also display in Fig. 2 the locations of the DSFGs identified by Xiao et al. (2024), selected based on their JWST colours F182M − F444W > 1.5 mag, which suggest significant dust attenuation. Notably, these JWST-selected galaxies lie systematically below the SFMS. This offset may stem from underestimated SFRs and/or overestimated M_★. In the absence of far-IR or (sub-)millimetre constraints (‘No IR’ case), UV-based SFRs, corrected for dust attenuation, frequently fail to capture the bulk of star formation in dusty systems. GOODS-ALMA studies (Elbaz et al. 2018; Xiao et al. 2023) report typical ratios of total to dust-corrected UV SFRs of the order of ∼8. Conversely, Pacifici et al. (2023) find that including far-IR data in SED fitting yield systematically lower SFRs compared to fits that omit such data. To investigate this further, we performed SED fitting on a subset of our N2CLS DSFGs using only optical to near-IR photometry (λ_obs < 10 μm). Using BAGPIPES (Carnall et al. 2018), we modelled the UV to near-IR SEDs of five sources with available HST+JWST photometry (N2GN_1_01, 06, 13, 23 and 43). For N2GN_1_01 (S2) and N2GN_1_23 (S3), we confirm the results of Xiao et al. (2024), with SFRs significantly underestimated. However, this is not the case for N2GN_1_43 where higher SFRs are recovered. To further probe this discrepancy, we repeated the analysis using the CIGALE code. The inferred SFRs from the No IR fits differ notably from those obtained with the full UV to (sub-)millimetre coverage, with differences reaching up to a factor of ∼20. SFRs are overestimated in the No IR scenario, except for two sources. Importantly, the stellar masses derived in the No IR case are systematically overestimated, with offsets ranging from a factor of 1.2 up to 7.5. These comparisons underscore that, even with the unprecedented depth of JWST, (sub-)millimetre data remain essential to accurately determining SFRs and M_★ (and A_V) in dusty galaxies. However, such data require very deep observations. Indeed, none of the 25 DSFGs in the GOODS-N field reported by Xiao et al. (2024, except for the brightest sources N2GN_1_01, 06, and 23) are detected at > 3σ in N2CLS, which already nearly reaches the confusion limit in GOODS-N (Ponthieu et al. 2025).

Fig. 2. Position of our eight N2CLS DSFGs in the SFR–M★ plane at z = 5.2, with sources at spectroscopic and photometric redshifts shown in red and blue, respectively. The SFMS from Popesso et al. (2023) and Speagle et al. (2014) is shown, with the 1σ scatter indicated by the light grey lines. For comparison, we include DSFGs from the GOODS-N field identified by Xiao et al. (2024) using JWST at 5.1 < z < 5.4 (excluding HDF850.1, GN10, and S3, which are part of our N2CLS sample). Galaxies from the 2 deg2 simulated infrared dusty extragalactic sky (SIDES) simulation (Béthermin et al. 2017), selected in a similar way as N2CLS with S​1.2 mm > 0.7 mJy at 5.1 < z < 5.4, are shown as downward orange triangles. The blue area indicates the location of the z = 5.2 DSFGs from Kumar et al. (2025), selected based on their 850 μm flux, using a flux threshold consistent with our N2CLS sample.

6. Stellar baryon budget and the efficiency of star formation

Recent observations have revealed a substantial population of massive galaxies in the early Universe, implying star-formation efficiencies that exceed the highest values inferred from lower-redshift studies. This may indicate a distinct and more intense mode of star formation during the earliest phases of galaxy evolution, which remains challenging for current models to reproduce (Shen et al. 2025, and references therein).

Studies leveraging data from the JWST have identified galaxies at redshifts z ∼ 5 to z ∼ 9 with stellar masses approaching 10¹¹ M_⊙. For instance, utilising JWST/Mid-infrared instrument (MIRI) data, Wang et al. (2025) found high densities of massive galaxies at z > 5, implying elevated stellar baryon efficiencies ϵ_★. Specifically, they report that the most massive galaxies at z ∼ 8 require ϵ_★ ∼ 0.3, compared to ϵ_★ ∼ 0.14, for typical low-redshift galaxies. Similarly, Xiao et al. (2024) examined 36 massive, dust-obscured galaxies with spectroscopic redshifts between z ∼ 5 and z ∼ 9 from the JWST FRESCO survey. They found that three ultra-massive galaxies (log(M_★/M_⊙)≳11.0, among them S2 ≡ N2GN_1_01 and S3 ≡ N2GN_1_23) require an exceptional ϵ_★ of approximately 50%, which is two to three times higher than the efficiencies observed in galaxies at later epochs. Further, Gentile et al. (2024) reported two massive, dusty starburst galaxies discovered in the COSMOS-Web survey, with star-formation efficiencies of ϵ_★ ∼ 0.25 at z_spec = 5.051 and ϵ_★ ∼ 0.5 at z_phot = 6.7 ± 0.3, respectively. In numerical simulations, ϵ_★ is found to increase significantly with both stellar mass and redshift, indicating that, at fixed mass, halos at (very) high redshifts convert gas into stars more efficiently than their lower-redshift counterparts (e.g. Ceverino et al. 2024). All these results imply that certain populations of high-redshift galaxies experience very efficient star formation at early times.

Building upon these studies, we computed the stellar baryon conversion efficiency (ϵ_★) for our sample of galaxies at z ∼ 5.2, following the methodology adopted in Boylan-Kolchin (2023) and Xiao et al. (2024). Using the cumulative halo mass function from Tinker et al. (2008), we computed the comoving number density of halos above a given halo mass threshold, at redshift z. The maximum stellar mass ($M_{★}^{\max }$) that is statistically allowed within an observed volume is given by $M_{★}^{\max }={ϵ}_{★}{f}_{\mathrm{b}}{M}_{\mathrm{halo}}^{\max }$, where ϵ_★ ≤ 1 represents the (integrated) efficiency with which baryons are converted into stars and f_b ≡ Ω_b/Ω_m = 0.156 is the cosmic baryon fraction from Planck Collaboration VI (2020). We show in Fig. 3 our 8 galaxies in the stellar mass versus redshift plane. The lines correspond to different number densities: 2.2 × 10⁻⁶ Mpc⁻³ (i.e. one source in the cosmic volume explored by N2CLS, assuming z = 5.2 ± 0.5; with ϵ_★ = 0.2 and ϵ_★ = 1 displayed as blue-dotted and grey-dashed lines, respectively) and 2.4 × 10⁻¹² Mpc⁻³ (i.e. a single galaxy in the whole sky; with ϵ_★ = 0.2 displayed as a dotted grey line).

Fig. 3. Stellar masses of our galaxies at z ∼ 5.2 compared to expectations from the simple relation that bounds the total stellar content of a dark matter halo to M★ ≤ ϵ★ fb Mhalo, where ϵ★ (which should be ≤1) is the efficiency with which baryons are converted into stars and fb is the cosmic baryonic fraction. We indicate the maximum stellar mass ($M_{★}^{\max }$) calculated from the maximum halo mass given our survey volume (160 arcmin2, z = 5.2 ± 0.5) for ϵ★ = 0.2 (dotted line) and ϵ★ = 1 (dashed grey line). We also show $M_{★}^{\max }$ for ϵ★ = 0.2 using the full sky area (dotted line).

Our analysis reveals notably high ϵ_★ values (ϵ_★ > 0.2 for all but N2GN_1_06 and 61), indicating that a significant fraction of baryonic mass in these galaxies has been converted into stars. The elevated values of ϵ_★ observed in our sample suggest that these galaxies are undergoing unusually efficient star formation at z = 5.2. Given their probable high metallicities and host halo masses exceeding 10¹² M_⊙, it is unlikely that cold gas inflows are the primary driver of this efficiency. Instead, the observed efficiency may reflect merger-driven inflows that channel gas into central regions and ignite intense starbursts, or compaction events that funnel enriched gas into central regions, triggering intense, short-lived starbursts.

7. Spatial distribution and abundance of N2CLS DSFGs in the GOODS-N overdensity

7.1. Contribution to the star formation rate density at z ∼ 5.2

Recent studies have shown that a significant fraction of the cosmic SFRD at high redshifts can be driven by rare, massive galaxies with high star formation efficiencies. For instance, Xiao et al. (2024) found that at z ∼ 5.8, ultra-massive galaxies with ϵ_★ > 0.2 can account for up to ${17}_{-3}^{+8}$% of the total SFRD (from Madau & Dickinson 2014). These findings highlight the critical role of massive, DSFGs in early cosmic star formation, motivating a closer look at their contribution in overdense environments such as GOODS-N.

Based on the five N2CLS galaxies with spectroscopic redshifts, we derive a dust-obscured cosmic SFRD⁹ of ${\psi }_{\mathrm{spec}}^{\mathrm{IR}}=(2.54\pm 0.10)\times {10}^{-2}$ M_⊙ yr⁻¹ Mpc⁻³. Including the additional contribution from sources with photometric redshifts raises the SFRD to ${\psi }_{\mathrm{spec}+\mathrm{phot}}^{\mathrm{IR}}=(4.77\pm 0.18)\times {10}^{-2}$ M_⊙ yr⁻¹ Mpc⁻³. The quoted uncertainties combine Poisson errors and uncertainties on individual SFR estimates in quadrature.

Figure 4 illustrates the contribution of the N2CLS galaxies in the z ∼ 5.2 overdensity to the cosmic SFRD. For reference, we compared our measurements to the widely used compilation of Madau & Dickinson (2014, MD14). The SFRD derived from spectroscopically confirmed galaxies, ${\psi }_{\mathrm{spec}}^{\mathrm{IR}}$, is 1.39 ± 0.06 times higher than the total SFRD reported by MD14 at similar redshift. When including galaxies with photometric redshifts, the excess increases to a factor of 2.62 ± 0.10. However, it is now well established that the MD14 compilation underestimates the SFRD at z ∼ 4 − 6. ALMA studies (e.g. Fujimoto et al. 2024; Traina et al. 2024; Sun et al. 2025; Liu et al. 2025) and JWST Hα luminosity function measurements (Covelo-Paz et al. 2025; Fu et al. 2025) both point to a systematically higher cosmic SFRD at these epochs. Using the updated estimate from Fujimoto et al. (2024), the contrast between our overdensity and the global average decreases to about unity for ${\psi }_{\mathrm{spec}}^{\mathrm{IR}}$ and to ∼2 for ${\psi }_{\mathrm{spec}+\mathrm{phot}}^{\mathrm{IR}}$. Restricting the comparison to the obscured component, our ${\psi }_{\mathrm{spec}}^{\mathrm{IR}}$ is already consistent with the levels reported by Traina et al. (2024) and Sun et al. (2025), and rises substantially above them once photometric redshift members are included. Relative to the dust-obscured SFRD from Zavala et al. (2021) at similar redshifts, ${\psi }_{\mathrm{spec}}^{\mathrm{IR}}$ alone is more than five times higher. If all galaxies with photometric redshifts are spectroscopically confirmed, the elevated SFRD in our sample will reflect the fact that these sources reside in a confirmed overdense environment, where star formation activity is strongly boosted compared to average field regions at the same epoch.

Fig. 4. Contribution of the dusty galaxies in the GOODS-N overdensity to the cosmic SFRD, ψ, at z ∼ 5.2. The red star represents the SFRD derived from spectroscopically confirmed galaxies, while the blue star also includes the contribution from galaxies with photometric redshifts. The purple curves show the total SFRD as a function of redshift: the dashed line corresponds to the Madau & Dickinson (2014) result, obtained by correcting the UV SFRD for dust attenuation, and the solid line to Fujimoto et al. (2024), obtained by summing the obscured (IR) and unobscured (UV) contributions. At z ∼ 5, the dust-obscured SFRD estimates from Traina et al. (2024) and Sun et al. (2025) are shown as an orange diamond and circle, respectively. The orange and blue shaded areas represent the obscured and unobscured SFRD, respectively, from Zavala et al. (2021).

Beyond our specific measurement, current research suggests that massive galaxy overdensities at high redshifts may be relatively common (e.g. Kashino et al. 2023; Wang et al. 2023; Helton et al. 2024a; Brinch et al. 2024) and may thus contribute significantly to the SFRD. In the GOODS-N field, Lin et al. (2025) estimate the fraction of the cosmic SFRD hosted within the protocluster environment at z ∼ 5.2, based on the number density of Hα emitters. They find that approximately 55% of the SFRD at this redshift is associated with the protocluster, somewhat higher than the ∼30% predicted by simulations (Chiang et al. 2017) but in good agreement with recent results based on the FLAMINGO simulation suite (Lim et al. 2024). Within the same volume considered in Lin et al. (2025, i.e. 62 arcmin² and 5.15 < z < 5.29), our four DSFGs yield a SFR density of log(ρ_SFR [M_⊙ yr⁻¹ Mpc⁻³]) = $-1.{39}_{-0.027}^{+0.026}$. This value is 1.6 times higher than the estimate based on the > 150 Hα emitters alone, suggesting that a substantial fraction of star formation activity is missed even in deep JWST surveys of highly overdense environments.

7.2. Environmental localization of N2CLS DSFGs in the GOODS-N overdensity

Recent JWST observations from the JADES and FRESCO programmes have revealed that the overdensity surrounding HDF850.1 represents one of the richest environments identified to date at z > 5, with approximately 100 galaxies at z = 5.17−5.20 distributed across 10 substructures within a ∼(15 cMpc)³ volume (Herard-Demanche et al. 2025; Sun et al. 2024). This region is particularly extreme, as ∼28% of all z = 5.0 − 6.0 Hα emitters detected across the GOODS-N field are concentrated within a narrow redshift slice of Δz ∼ 0.03. The spatial extent of this protocluster is consistent with the typical size of such structures predicted by Chiang et al. (2017) using a set of N-body simulations and semi-analytic models. In addition, Calvi et al. (2023) identified the galaxy GN10 (N2GN_1_01) as residing in an overdense region based on photometric redshifts, a result that has now been confirmed spectroscopically by Sun et al. (2024).

We show the locations of the N2CLS galaxies within the overdensity in Fig. 5. Cosmological simulations suggest that at z ≳ 5, approximately 90% of baryonic and dark matter resides in unvirialised, filamentary large-scale structures (e.g. Haider et al. 2016). Sun et al. (2024) interpret the 3D spatial distribution of group member galaxies as being consistent with three distinct filamentary structures. For simplicity, we compared the positions of our galaxies relative to the filaments identified in Sun et al. (2024), rather than to the 18 substructures detailed in Herard-Demanche et al. (2025). The z = 5.17 − 5.20 overdensity, traced by filament-1 and filament-2, spans the full extent of the GOODS-North FRESCO field (62 arcmin²). Within this region lie four N2CLS DSFGs with spectroscopic redshifts: N2GN_1_06 and N2GN_1_23 are located in filament-1, while the newly identified N2GN_1_13 and N2GN_1_61 fall within filament-2. In filament-3, N2GN_1_26 (photometric redshift) is found along with GN10, and N2GN_1_43 lies nearby. N2GN_1_57 may trace a distinct, possibly new, structure. We note that some previously confirmed spectroscopic members of the z ∼ 5.2 protocluster (Walter et al. 2012; Calvi et al. 2021) fall outside the JADES and FRESCO joint footprint, implying that the overdensity exhibits more complex substructures than currently captured.

Fig. 5. Large-scale structure of the overdense environment in the GOODS-N field at z = 5.17 − 5.40. Galaxy positions are shown relative to HDF850.1, in units of comoving Mpc. All galaxies are colour-coded by redshift. The eight N2CLS galaxies are indicated by crosses (for spectroscopic redshifts, zspec) and upward pointing triangles (for photometric redshifts, zphot); their NIKA2 identifiers are labelled in red, with source 20 discussed in Appendix D. Galaxy groups and filaments are taken from the JADES survey (Sun et al. 2024), with galaxies in less dense regions shown as transparent symbols. Additional galaxies, shown as dots, are drawn from the FRESCO survey (Herard-Demanche et al. 2025).

We also show the location of N2GN_1_20, another N2CLS galaxy at z_phot = 5.33 (see Appendix D), which lies farther from the core overdensity. Interestingly, Jiang et al. (2018) report a giant protocluster at z ∼ 5.7 occupying a volume of ∼35³ cMpc³, embedded within a larger overdense region spanning 53 × 41 cMpc² and extending 34 cMpc along the line of sight. Further investigation is needed to determine whether the overdensity could extend towards the region traced by N2GN_1_20. Finally, we note that the spatial distribution of N2CLS galaxies at z ∼ 5.2, when viewed in two dimensions (RA, Dec), does not exhibit a recognisable filamentary pattern, unlike the tentative identification of a filament with NIKA2 in the GJ526 field (Lestrade et al. 2022).

8. Accelerated galaxy evolution in the z ∼ 5.2 overdense structure

Chiang et al. (2017) outlined three evolutionary stages of cluster formation: an early inside-out growth phase from z = 10 to z ∼ 5, an extended period of vigorous star formation from z ∼ 5 to z ∼ 1.5, and a final stage of rapid mass assembly and quenching from z ∼ 1.5 to z = 0. Our results are consistent with the second phase of this sequence, during which intense star formation is expected to occur across 10−20 cMpc-scale structures, contributing significantly to the stellar mass of present-day clusters. Interestingly, Lim et al. (2024), using the FLAMINGO simulation suite, suggest that the suppression of star formation in protocluster environments already began by z ≈ 5, indicating that this epoch could represent a transitional phase in cluster evolution.

Recent results from the FOREVER22 project (Morokuma-Matsui et al. 2025) indicate that the elevated star formation activity observed in protocluster galaxies at z ≥ 6, compared to field galaxies, is primarily driven by a higher abundance of massive galaxies rather than enhanced sSFRs. This trend is likely a consequence of the denser environments in protoclusters, where increased galaxy interactions and sustained gas accretion promote the rapid assembly of massive systems. These massive galaxies tend to exhibit stronger dust attenuation, resulting in redder UV continuum slopes. This is consistent with the findings of Helton et al. (2024b), who reported that the UV-brightest and UV-reddest galaxies are preferentially located in denser environments with a higher number of neighbouring galaxies. Together, these observations provide compelling evidence that galaxy evolution is accelerated within high-redshift overdense structures and may help explain the record number of bright DSFGs we identified in this particular structure.

In the FoREVER22 project, protocluster galaxies show a combination of stellar masses and SFRs that are somewhat lower, but still comparable to those of the DSFGs identified in our study. This places our DSFGs at the upper end of the FoREVER22 population in terms of total star formation activity. In contrast, the FLAMINGO simulations from Lim et al. (2024) produce protocluster galaxies with high SFRs that are consistent with the SFRs measured for our sample. This is not the case for the general DSFG population modelled in FLAMINGO (Kumar et al. 2025), selected based on their 850 μm flux (S​_850 μm > 0.3 mJy, corresponding to the N2CLS flux cut S_1.2,mm > 0.7 mJy at z = 5.2), which predominantly lies on the SFMS (see Fig. 2). Furthermore, Chiang et al. (2017) estimate that the total SFR in protoclusters at z ∼ 5 is typically below 1000 M_⊙ yr⁻¹, which is about four times lower than the combined SFR of our five spectroscopically confirmed DSFGs. The unusually high level of star formation highlights the exceptional nature of this structure.

In the context of the SIDES simulation (Béthermin et al. 2017; Gkogkou et al. 2023), where no environmental effects are considered, we can use this dataset to address whether the enhanced dusty star formation observed here is driven by clustering or by the presence of large gas masses. To explore this, we selected SIDES galaxies in the 2-deg² simulation in the redshift range 5.1 < z < 5.4 and with a NIKA2 1.2 mm flux greater than 0.7 mJy, which mimics the selection criteria for our sample (B25). The SFR–M_★ relation for a NIKA2-like population is shown in Fig. 2. With the exception of N2GN_1_61, SIDES galaxies generally exhibit lower SFR and M_★. In SIDES, the SFR is capped at a hard limit of 1000 M_⊙/yr, and each galaxy’s SFR is redrawn until it falls below this threshold. This is explaining the lack of SFR > 1000 M_⊙/yr galaxies in SIDES. However, we can also see that SIDES does not reach high stellar masses. The most massive galaxy in SIDES has a stellar mass of log(M_★/M_⊙) = 10.92. When considering the larger 117-deg² SIDES-UCHUU simulation (Gkogkou et al. 2023), the most massive galaxy reaches log(M_★/M_⊙) = 11.29. SIDES employs abundance matching to obtain the stellar mass, which connects halo to galaxy stellar mass via mass functions. This approach may suppress the presence of high-stellar mass objects if they are not represented in the observed stellar mass function used to build the SIDES simulation. Indeed, using the stellar mass function from SIDES, we expect 0.4 galaxies with log(M_★/M_⊙) > 10.7 at 5.1 < z < 5.4 in the GOODS-N area (160 arcmin²), whereas we observe six (or seven if N2GN_1_20 is included).

Given the non-uniform distribution of sources, we investigated the impact of field-to-field variance by dividing the SIDES-UCHUU simulation into 2562 independent sub-fields, each covering 160 arcmin². Among these, 13 sub-fields contain 4 galaxies above the stellar mass threshold, a small fraction but consistent with our spectroscopic (z_spec) sample. To evaluate the effect of stellar mass uncertainties, we added a random scatter of 0.5 dex to the SIDES stellar masses and repeated the sub-field analysis. With this added uncertainty, the average stellar mass of the most massive galaxy in the redshift range 5.1 < z < 5.4 across all sub-fields is log(M_★/M_⊙) = 11.28, consistent with the most massive spectroscopically confirmed galaxy in our sample. Furthermore, we now identify 25 sub-fields (1% of the total) that contain 8 galaxies with log(M_★/M_⊙) > 10.7 within 5.1 < z < 5.4. In conclusion, while such overdensities are rare, their existence is plausible when considering both the strong clustering of massive galaxies at z > 5 and the significant uncertainties in stellar mass estimates.

Based on their stellar masses, the five z_spec galaxies are each expected to reside in massive dark matter halos, with a combined halo mass of M_halo ∼ 5 × 10¹³ M_⊙, as inferred from the empirical M_★ − M_halo relation at z = 5 (Behroozi et al. 2019). While it is uncertain whether all these galaxies will eventually merge into a single structure by z = 0, such an outcome would imply a descendant halo more massive than that of a typical Coma-like cluster progenitor at the same redshift, which simulations suggest is closer to ∼10¹³ M_⊙ (Chiang et al. 2013). This highlights the exceptional mass scale of the environment traced by these galaxies. However, it is important to note that total mass estimates for protoclusters are sensitive to the choice of aperture and redshift interval used for their definition, influencing such comparisons across studies (Lim et al. 2024).

9. Conclusion

Using 1.2 mm observations of GOODS-N taken with NIKA2, we identified eight N2CLS galaxies at z ∼ 5.2 (five with z_spec and three with z_phot), many of which are extremely faint or even undetected in deep HST and some JWST imaging. These galaxies have remarkably high SFRs, substantial dust content, and significant stellar masses, placing them among the most extreme star-forming systems known at this epoch. Six are securely associated with the exceptional overdensity reported in Sun et al. (2024) and Herard-Demanche et al. (2025).

We first identified a new DSFG in the overdense environment of HDF850.1 through NOEMA observations, which detect CO(4−3) and CO(5−4) lines at z = 5.182. In addition, another N2CLS galaxy shows a 7σ emission line in FRESCO data, consistent with Hα at z_spec = 5.201. Combined with previously known galaxies (HDF850.1, GN10, and S3), this brings the number of spectroscopically confirmed DSFGs within the overdense structure to five. We further identified three DSFGs with photometric redshifts consistent with the overdensity.

These galaxies exhibit exceptionally high SFRs (ranging from ∼200 to 2000 M_⊙ yr⁻¹) and high stellar masses (log(M_★/M_⊙) = 10−11), placing most of them significantly above the SFMS at z ∼ 5.2, especially under the Popesso et al. (2023) parametrisation. Their properties are comparable to those of the ultra-massive galaxies analysed by Xiao et al. (2024), suggesting an early and efficient buildup of both stellar and dust mass. They show short gas depletion timescales (with a median of 0.18 Gyr) and remarkably high stellar baryon conversion efficiencies (ϵ_★ > 20%, compared to ≲10% locally). The observed dust masses further support the scenario of efficient dust formation. Due to their short gas depletion timescales, these DSFGs could evolve into quiescent galaxies by z ∼ 4.5, assuming no new gas inflow. This scenario aligns with the early formation of massive quiescent galaxies seen in similar dense environments.

Their contribution to the cosmic SFRD is significant: even considering that it represents only the obscured component, the five spectroscopically confirmed galaxies exceed the total SFRD from Madau & Dickinson (2014) by a factor of 1.4, which rises to 2.6 when including all eight sources. Compared to more recent SFRD estimates at z ∼ 5 (e.g. Fujimoto et al. 2024), the excess is reduced to roughly unity for the spectroscopic sample and by about a factor of 2 when including photometric redshift members. These few dusty galaxies dominate the star formation budget within the overdensity, contributing more than the much larger population of Hα emitters.

Our DSFGs either lie at the upper extremes of stellar mass, dust mass, and SFR or are entirely absent from current simulations and models. This highlights the difficulty existing models face in reproducing such extreme systems. Their properties may be driven by efficient gas inflows along cosmic filaments in protocluster environments fuelling intense star formation despite short depletion timescales (e.g. Narayanan et al. 2015; Umehata et al. 2019).

We examined the spatial distribution of N2CLS galaxies within the overdensity, comparing their positions to the three filamentary structures identified by Sun et al. (2024). Four spectroscopically confirmed DSFGs lie along filament-1 and filament-2, two others align with filament-3, and two may trace new separate structures. We also identified an extremely massive (∼3 × 10¹² M_⊙) galaxy candidate at z = 5.33 that lies farther away (see Appendix D). The protocluster may therefore be more extended than previously mapped, as some galaxies lie beyond current survey footprints. Confirming this will require spectroscopic redshifts for sources currently identified via photometric estimates. If all these galaxies merged by z = 0, they would form a halo significantly more massive than a typical Coma-like cluster progenitor at that epoch, emphasising the exceptional mass scale of the environment traced by these galaxies.

The presence of such massive, dusty, and intensely star-forming galaxies at z ∼ 5.2 has significant implications for our understanding of galaxy formation and evolution. It suggests that the processes leading to rapid stellar and dust mass assembly were already in place within the first billion years of cosmic history in overdense environments. Future observations, particularly with facilities capable of probing the ISM properties and kinematics of these galaxies, will be crucial in unravelling the mechanisms driving their evolution.

Acknowledgments

We acknowledge financial support from the Programme National de Cosmologie and Galaxies (PNCG) funded by CNRS/INSU-IN2P3-INP, CEA and CNES, France, and from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (project CONCERTO, grant agreement No 788212). This work is based on observations carried out under project numbers 192-16 with the IRAM 30-m telescope, and projects W21CV, W23CX, and S24CF with NOEMA. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). 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. The observations are associated with programs #1181, 1895 and 3577. Some of the data products used here were retrieved from the Dawn JWST Archive (DJA). DJA is an initiative of the Cosmic Dawn Center (DAWN), which is funded by the Danish National Research Foundation under grant DNRF140. We thank Eric F. Jiménez-Andrade for sharing their data products from the 10 GHz Karl G. Jansky Very Large Array survey of GOODS-N. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. GL thanks Joris Witstok and Thomas Henning for enlightening discussions on dust production, and warmly thanks Fengwu Sun for sharing his Python code used to generate Figure 4.

References

Appendix A: Dust production

Observations have shown a relation between dust mass and stellar mass across an extensive redshift range. At all redshifts, the dust masses increase linearly with stellar masses. However, at fixed stellar mass, there is a large scatter of dust masses around the median, revealing the existence of a dust-poor galaxy population across all stellar masses (Triani et al. 2020). We compare in Fig. A.1 our stellar-mass estimates with our inferred dust masses for the eight z ∼ 5.2 DSFG. We also show the stellar mass–dust mass relations derived by Witstok et al. (2023) based on dust production from stellar sources (including supernovae and AGB stars; with and without reverse shock destruction, which can lead to ∼95% dust destruction). Additionally, we include the relation corresponding to the highest possible metal yield (i.e. assuming all metals are locked up in dust grains), representing the theoretical maximum dust mass. Note that the exact normalisation of these relations remains uncertain due to systematic effects, including the choice of IMF used to compute dust yields for a given stellar mass, and the range of plausible dust absorption cross-sections.

The dust masses observed in our galaxies point to efficient dust formation. Indeed our galaxies are consistent with the dust masses expected from a simple, IMF-averaged yield of stellar sources alone (i.e. supernovae and AGB stars), assuming no dust destruction in the ISM or ejection from the galaxy. In general, it is challenging to populate our observed M_★–M_dust relation using cosmological simulations, as our dusty galaxies occupy the upper envelope in both stellar and dust mass. They are, for instance, only sparsely represented in Di Cesare et al. (2023) and entirely absent in the sample of Triani et al. (2020).

We also show in Fig. A.1 the location of N2GN_1_20, the most massive galaxy in our z∼5.2 sample, which is discussed in Appendix D. It falls into the region where reverse-shock destruction of dust becomes significant, implying that additional dust growth mechanisms, such as grain growth in the ISM, may be necessary to account for their observed dust masses. However, we caution that this source currently lack spectroscopic redshift, and robust spectroscopic confirmation is essential before drawing firm conclusions about its dust production.

For comparison, we also present the complete sample of N2CLS galaxies at 2 < z < 5 from Berta et al. (2025), which occupy the same region i.e. near the supernova+AGB dust production line. Notably, three galaxies lie above the maximum dust mass allowed: N2GN_1_17a (z_spec = 3.19) and N2GN_1_32 (z_spec = 3.65; detected in the Chandra 2 Ms map and requiring an AGN-torus component in the IRAC and MIPS bands) remain consistent within 1.3σ (in M_dust), while N2GN_1_04 (z_spec = 4.056, also known as GN20) is significantly offset.

The distribution of our measurements in the M_dust–M_★ plane appears relatively flat. Sommovigo et al. (2022) show that the M_dust–M_★ relation differs between z ∼ 7 REBELS galaxies and z ∼ 5 ALPINE galaxies (even when using consistent [CII]-based methods to estimate M_dust), with REBELS galaxies exhibiting a flat trend, while ALPINE galaxies, as well as theoretical models (see Fig. A.1 and their Fig. 3), suggest that more massive galaxies are dustier. Both samples probe lower dust and stellar masses than our NIKA2 sources. The observed flatness of the M_dust–M_★ relation could partly arise from observational biases: the REBELS survey covers a limited stellar mass range, while for NIKA2, the 1.2 mm flux selection imposes a dust mass detection threshold.

Fig. A.1. Dust mass versus stellar mass for N2CLS galaxies. Filled red circles represent the five DSFGs with spectroscopic redshifts at z ∼ 5.2, while open blue circles indicate those with photometric redshifts (including the very massive N2GN_1_20 galaxy discussed in Appendix D). For comparison, we also show all N2CLS galaxies at 2 < z < 5 from Berta et al. (2025). The region between the orange and purple lines is the expected dust mass range assuming a maximal dust production from both AGB stars and supernovae (SN), without significant reverse-shock destruction. The area below the orange line represents the same scenario, but accounting for ∼95% reverse-shock destruction. The blue line shows the upper limit on dust mass based on the highest metal yields from Witstok et al. (2023), i.e. the extreme case of metal accretion onto dust grains resulting in maximal grain growth.

Appendix B: Photometric redshifts with New-Hyperz

For the four sources with literature-based photometric redshifts at z ∼ 5.2 (N2GN_1_20, N2GN_1_26, N2GN_1_43, N2GN_1_57), we re-derived their redshifts independently using the updated photometry and the New-Hyperz SED-fitting code. New-Hyperz is a modified version of the public code Hyperz originally developed by Bolzonella et al. (2000). The library of templates used for this need includes: (i) Eight evolutionary synthetic SEDs computed with the Bruzual & Charlot code (Bruzual & Charlot 2003), with Chabrier IMF (Chabrier 2003) and solar metallicity, matching the observed colours of local galaxies from E to Im types (namely a delta burst single stellar population, a constant star-forming system, and six τ-models with exponentially decaying SFR); (ii) Six empirical templates for galaxies in the local Universe: four SEDs compiled by Coleman et al. (1980), and two starburst galaxies from the Kinney et al. (1996) library; (iii) Ten Starbursts99 evolutionary models (Leitherer et al. 1999), including e-lines, for single bursts and constant SFR models, each spanning five metallicities (Z = 0.04, 0.02, 0.008, 0.004 and 0.001). The explored redshift domain goes from z = 0 to 16, without any prior in luminosity. The only restriction is that the age of the stellar population should be smaller than the age of the Universe for a given redshift. The internal dust extinction is a free parameter that is allowed to take values between A_V = 0 and 6 magnitudes following the Calzetti et al. (2000) extinction law. During the SED-fitting process, to account in a realistic way for calibration uncertainties related to the combination of data coming from different facilities, across a wide range in wavelengths, we imposed a minimum photometric error of 0.2 magnitudes, which is quadratically added to the original error bars. Indeed, in our experience this is needed to prevent artificially small error bars leading to low-quality fits, by correcting from effects such as biases between different instruments due to (absolute) zero-point uncertainties, photometry computed with an effectively different spatial sampling and spatial resolution, and contribution to the total flux coming from different spatial regions.

The quality of the photometric redshifts with New-Hyperz and the settings described above for the SED-fitting has been checked using the 29 N2CLS galaxies with spectroscopic redshifts (Sect. 2), excluding N2GN_1_61 because it is only detected in one filter. Spectroscopic redshifts are in the range between z = 0.6 and 7.2. The results are as follows, with Δz = z_phot − z_spec: σ(Δz/(1 + z_spec)) = 0.058, the median of (Δz/(1 + z_spec)) = -0.016, and the normalised median absolute deviation, which is less sensitive to outliers: σ_z, MAD = 1.48 × median (|Δz|/(1 + z_spec)) = 0.018. Regarding the fraction of outliers, defined in a conservative way as sources with Δz > 0.15(1 + z_spec), it is found to be 21% (6 of 28 objects), and only two galaxies (7%) exhibit |Δz|> 0.30(1 + z_spec). Among these six outliers, four exhibit a degenerate solution with a secondary solution consistent with z_spec. One of them is an X-ray AGN at z_spec = 3.65 with a near-IR to mid-IR SED dominated by a torus component, exhibiting a degenerate solution between z${}_{\mathrm{phot}}=0.{085}_{-0.085}^{+0.038}$ and z${}_{\mathrm{phot}}=2.{960}_{-0.04}^{+0.11}$. Also, a poor z_phot estimate is obtained for an obscured quasar at z_spec = 7.19 (Fujimoto et al. 2022).

Appendix C: Notes about individual sources and JWST cutouts

Below, we provide a detailed description of each individual source identified at z∼5.2 in the overdensity (but HDF850.1, GN10, and N2GN_1_13). The JWST photometry measurements follow a similar approach to that described by Weibel et al. (2024), with some modifications as outlined in Sect. 2.8 of B25. The multi-wavelength cutouts of all sources are shown in B25. We complement those cutouts with the JWST postage stamps, where available, in Fig. C.1.

• N2GN_1_23: is detected by JWST FRESCO at z_spec = 5.179 and is the so-called S3 source extensively discussed in Xiao et al. (2024). We have a firm identification of our N2CLS source through our NOEMA follow-up programme (see B25) with the NOEMA counterpart lying at a distance of 0.167 arcsec. The NIKA2 source has both a SCUBA2 450 and a 850 μm counterpart but the NOEMA position lies at 3.6 arcsec from the 450 μm source Number 78 of Barger et al. (2022), identified with a galaxy at z = 0.29 (see their Table 2). The 850 μm source is identified with SMM123656621205 (Number 80), which is at 2.8 arcsec from the NOEMA position, and with a redshift estimate of 3.75 (Cowie et al. 2017). The high-z source is clearly blended with a nearby object and flux measurements from Herschel and SCUBA2 450 μm cannot be used to construct the multi-wavelength SED. The NOEMA source has both a 1.4 GHz (Owen 2018) and 10 GHz (Jiménez-Andrade et al. 2024) radio counterpart.

• N2GN_1_26 is identified with a high-z HST galaxy thanks to SMA data. Kodra et al. (2023) give a photometric redshift z_phot = 5.31 ± 0.08. The far-IR 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 polycyclic aromatic hydrocarbon feature. There is a 10 GHz source at the HST position. This source is newly observed and detected by the PANORAMIC survey with F277W, F356W. Using New-Hyperz, it has a photometric redshift of 4.87${}_{-0.19}^{+0.21}$, consistent with previous findings at 1σ level.

• N2GN_1_43 is identified with NOEMA. It is a K-band dark galaxy. It is detected by JWST at ∼10σ in the F277W band, ∼3σ in F200W, and tentatively by HST at shorter wavelengths. Although not in the IRAC catalogue, it is clearly seen in the maps at 4.5, 5.8 and 8 μm. It has also a clear 10 GHz signal (Jiménez-Andrade et al. 2024) at the JWST source position. A fit with CIGALE in B25 places it at z_phot = 5.87 ± 0.38 with a ∼4000 Å break between the two JWST measurements mentioned here. The photometric redshift from New-Hyperz is 5.31${}_{-0.86}^{+1.33}$, consistent with previous estimates, with a degenerated solution in the parameter space between redshift and extinction leading to broad error bars.

• N2GN_1_57: The VLA source, coinciding with GN4 26261 (MIOP programme) at z = 0.9, is likely not the right counterpart as it lies at 4.9 arcsec from the N2CLS position. Closer to the N2CLS position (2.2 arcsec), the source HRG14 J123731.66+621616.7 (also detected by the SHARDS middle-band survey; Finkelstein et al. 2015) is at a photometric redshift of z = 5.2±0.06 in Arrabal Haro et al. (2018) and z = 5.30${}_{-0.19}^{+017}$ in Kodra et al. (2023). It is detected in the four IRAC channels and has a pronounced ∼4000 Å break at λ ∼ 2 μm. The photometric redshift obtained with New-Hyperz is z = 5.18${}_{-0.20}^{+0.19}$, consistent with previous findings.

• N2GN_1_61 is not only optically dark but also IRAC-dark ([4.5]> 24.6 mag at 5σ). The accurate position is given by NOEMA (see B25) revealing a single source with a flux of 698±84 μJy at 1.18 mm, in perfect agreement with the NIKA2 flux. JWST grism observations from FRESCO detected one ∼7 σ emission line, consistent with Hα at z_spec = 5.201 (M. Xiao, priv. comm.). It is observed by JWST with F090W, F115W, F182M, F210M, F356W, F444W filters, but only detected by F356W and F444W (see Fig. C.1). 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 allows us to perform a fit. Derived parameters have large uncertainties (but are consistent between MAGPHYS-highz and CIGALE).

Fig. C.1. JWST/NIRCam 10"×10" cutout images centred on N2GN_1_01, 06, 13, 23, 43, and 61 (from top to bottom). The filters used are indicated at the top of each corresponding multi-colour image (on the right). White contours represent the NOEMA signal-to-noise ratio, ranging from 3σ to 19σ in steps of 2. N2GN_1_20, 26, and 57 are not covered by current JWST data.

Appendix D: Another N2CLS candidate at z∼5.2

We identified an additional N2CLS galaxy at z∼5.2 that lies close to the overdense structure. N2GN_1_20 is detected by SMA. SMA position (RA = 12:35:50.22; Dec = 62:10:42.4) lies at 1.2 arcsec from S-CANDELS J123550.35+621041.8, which is an optically dark galaxy. It has a Ks band flux of 0.75 μJy, resulting in a very red 3.6 μm to Ks colour (S​_3.6/S​_Ks ∼ 13), interpreted as a deep break at ∼4000 Å combined to a large extinction, that leads to a z_phot = 5.333 ± 0.317 (Liu et al. 2018). Using New-Hyperz, we obtained z_phot = 5.16${}_{-1.10}^{+2.8}$, the broad error bars being the result of a degenerated solution in the parameter space between redshift, extinction and age of the stellar population. This source is not covered by JWST.

We show in Fig. D.1 the best fit from CIGALE, assuming z = 5.33. The parameters derived from the fit are: SFR = 2160±611 M_⊙ yr⁻¹; A_V = 2.7±0.3; M_★ = (303±66)×10¹⁰ M_⊙; and M_dust = (23.3±19.70)×10⁸ M_⊙. Combining the derived SFR and M_★, and applying the Tacconi et al. (2020) scaling relations, we obtain M_gas= 103×10¹⁰ M_⊙ and τ_dep = 476 Myr. Assuming its photometric redshift, the galaxy is extremely massive. As such, it is lying between the Popesso et al. (2023) and Speagle et al. (2014) parametrisations of the MS.

Fig. D.1. N2GN_1_20: best fit from CIGALE assuming z = 5.33 (top) and z = 3 (bottom). Observed and modelled fluxes are shown with purple and red circles, respectively.

It is also lying above the model expectations for the GOODS-N survey volume (Fig. 3), requiring ϵ_★ = 0.7 for the full sky area. This discrepancy could arise from assumptions related to the IMF. For example, a top-heavy IMF could yield lower stellar masses by ≳0.3 dex, resulting in slightly reduced values of ϵ_★. It could also arise from an overestimation of its photometric redshifts. A substantially lower true redshift compared to the photometric estimate may result in a lower best-fit stellar mass, thereby reducing the required value of ϵ_★. To test this hypothesis, we modelled the SED of N2GN_1_20 with CIGALE assuming z = 3 (see Fig. D.1). We obtain a stellar mass of (3.15±0.34)×10¹² M_⊙, essentially comparable to that inferred at z = 5.33. It is important to note that the fit at z = 3 is not merely a rescaling of the same template at a different redshift: CIGALE recomputes the fit and searches for an optimal solution at the imposed z. The reduced χ² is 30% higher at z = 3 than at z = 5.33. Such a high stellar mass at z = 3 remains above model expectations for the GOODS-N field area, and would still require ϵ_★ = 0.12 for the full sky. If confirmed at z = 5.33, the extremely high mass could reflect the presence of an unusually dense concentration of massive halos in this region.

We also compare the SED of N2GN_1_20 with those of the eight galaxies discussed in the main text, using different normalisations (Fig. D.2). We do not find any indication that N2GN_1_20 deviates significantly from the others.

As shown in Fig. 5, N2GN_1_20 lies outside the core of the overdensity. Its distance to GN10 is ΔRA = 11.8 cMpc, ΔDec=-8.1 cMpc, and ΔLOS_distance= 14.9 cMpc. Further investigation is required to determine whether this optically dark galaxy truly lies at high redshifts, to understand its unusually high stellar mass and efficiency, and to assess whether it may trace an extended region of the overdensity.

Fig. D.2. SED of N2GN_1_20 at z = 5.33 (red stars), compared to the eight N2CLS galaxies in the overdensity at z∼5.2 (circles), shown as a function of rest-frame wavelength. From top to bottom, the SEDs are normalised at 200, 3.6, and 0.9 μm (rest frame).

All Tables

All Figures

Fig. 1. N2GN_1_13. Left: False colour image produced using HST data, with F814W, F125W, and F160W images in blue, green, and red, respectively. Overlaid in red are the NOEMA continuum contours at 2 mm, ranging from 4σ to 14σ in steps of 2. The white contours are JWST F444W brightness. The scale bar in the bottom left has a length of 1″. Middle: NOEMA CO(4−3) and CO(5−4) line spectra. The red line shows the Gaussian fit to the line profile (including the continuum). Notice that the CO(5−4) line lies precisely at the interface between two basebands, resulting in the partial fragmentation of the spectral profile. Right: NOEMA 3.5 mm continuum map. Contours are displayed with a step of 12 μJy/beam (corresponding to 2σ).

In the text

Fig. 2. Position of our eight N2CLS DSFGs in the SFR–M★ plane at z = 5.2, with sources at spectroscopic and photometric redshifts shown in red and blue, respectively. The SFMS from Popesso et al. (2023) and Speagle et al. (2014) is shown, with the 1σ scatter indicated by the light grey lines. For comparison, we include DSFGs from the GOODS-N field identified by Xiao et al. (2024) using JWST at 5.1 < z < 5.4 (excluding HDF850.1, GN10, and S3, which are part of our N2CLS sample). Galaxies from the 2 deg2 simulated infrared dusty extragalactic sky (SIDES) simulation (Béthermin et al. 2017), selected in a similar way as N2CLS with S​1.2 mm > 0.7 mJy at 5.1 < z < 5.4, are shown as downward orange triangles. The blue area indicates the location of the z = 5.2 DSFGs from Kumar et al. (2025), selected based on their 850 μm flux, using a flux threshold consistent with our N2CLS sample.

In the text

Fig. 3. Stellar masses of our galaxies at z ∼ 5.2 compared to expectations from the simple relation that bounds the total stellar content of a dark matter halo to M★ ≤ ϵ★ fb Mhalo, where ϵ★ (which should be ≤1) is the efficiency with which baryons are converted into stars and fb is the cosmic baryonic fraction. We indicate the maximum stellar mass ($M_{★}^{\max }$) calculated from the maximum halo mass given our survey volume (160 arcmin2, z = 5.2 ± 0.5) for ϵ★ = 0.2 (dotted line) and ϵ★ = 1 (dashed grey line). We also show $M_{★}^{\max }$ for ϵ★ = 0.2 using the full sky area (dotted line).

In the text

Fig. 4. Contribution of the dusty galaxies in the GOODS-N overdensity to the cosmic SFRD, ψ, at z ∼ 5.2. The red star represents the SFRD derived from spectroscopically confirmed galaxies, while the blue star also includes the contribution from galaxies with photometric redshifts. The purple curves show the total SFRD as a function of redshift: the dashed line corresponds to the Madau & Dickinson (2014) result, obtained by correcting the UV SFRD for dust attenuation, and the solid line to Fujimoto et al. (2024), obtained by summing the obscured (IR) and unobscured (UV) contributions. At z ∼ 5, the dust-obscured SFRD estimates from Traina et al. (2024) and Sun et al. (2025) are shown as an orange diamond and circle, respectively. The orange and blue shaded areas represent the obscured and unobscured SFRD, respectively, from Zavala et al. (2021).

In the text

Fig. 5. Large-scale structure of the overdense environment in the GOODS-N field at z = 5.17 − 5.40. Galaxy positions are shown relative to HDF850.1, in units of comoving Mpc. All galaxies are colour-coded by redshift. The eight N2CLS galaxies are indicated by crosses (for spectroscopic redshifts, zspec) and upward pointing triangles (for photometric redshifts, zphot); their NIKA2 identifiers are labelled in red, with source 20 discussed in Appendix D. Galaxy groups and filaments are taken from the JADES survey (Sun et al. 2024), with galaxies in less dense regions shown as transparent symbols. Additional galaxies, shown as dots, are drawn from the FRESCO survey (Herard-Demanche et al. 2025).

In the text

Fig. A.1. Dust mass versus stellar mass for N2CLS galaxies. Filled red circles represent the five DSFGs with spectroscopic redshifts at z ∼ 5.2, while open blue circles indicate those with photometric redshifts (including the very massive N2GN_1_20 galaxy discussed in Appendix D). For comparison, we also show all N2CLS galaxies at 2 < z < 5 from Berta et al. (2025). The region between the orange and purple lines is the expected dust mass range assuming a maximal dust production from both AGB stars and supernovae (SN), without significant reverse-shock destruction. The area below the orange line represents the same scenario, but accounting for ∼95% reverse-shock destruction. The blue line shows the upper limit on dust mass based on the highest metal yields from Witstok et al. (2023), i.e. the extreme case of metal accretion onto dust grains resulting in maximal grain growth.

In the text

	Fig. C.1. JWST/NIRCam 10"×10" cutout images centred on N2GN_1_01, 06, 13, 23, 43, and 61 (from top to bottom). The filters used are indicated at the top of each corresponding multi-colour image (on the right). White contours represent the NOEMA signal-to-noise ratio, ranging from 3σ to 19σ in steps of 2. N2GN_1_20, 26, and 57 are not covered by current JWST data.
In the text

	Fig. D.1. N2GN_1_20: best fit from CIGALE assuming z = 5.33 (top) and z = 3 (bottom). Observed and modelled fluxes are shown with purple and red circles, respectively.
In the text

	Fig. D.2. SED of N2GN_1_20 at z = 5.33 (red stars), compared to the eight N2CLS galaxies in the overdensity at z∼5.2 (circles), shown as a function of rest-frame wavelength. From top to bottom, the SEDs are normalised at 200, 3.6, and 0.9 μm (rest frame).
In the text