© The Authors 2026

1. Introduction

Since the discovery of high-redshift bright submillimeter (submm) galaxies (e.g., Smail et al. 1997; Hughes et al. 1998; Blain et al. 2002; Chapman et al. 2003), the (sub)mm atmospheric window is the prime wavelength range to probe dusty star-forming systems in the high-z Universe (z > 2; see, e.g., Casey et al. 2014 and Hodge & da Cunha 2020 for reviews). Indeed, the peak of cold dust emission in the far-infrared (FIR, ∼80 μm) tracing obscured star formation is shifted to the submm windows by the redshift. The observed submm flux from galaxies at higher redshifts comes from rest-frame wavelengths closer to their peak of emission, compensating for the increasing luminosity distance and leading to an almost constant flux-versus-luminosity ratio at z≳ 2 (e.g., Blain et al. 2002; Lagache et al. 2004). This effect is commonly referred to as the negative K-correction.

Observations have unveiled a numerous population of these submm galaxies that pose a challenge to galaxy evolution models, which struggle to produce enough objects with such high star formation rates in the early Universe (e.g., Baugh et al. 2005; Lovell et al. 2021; Hayward et al. 2021; Kumar et al. 2025). The Atacama large millimeter array (ALMA) showed that massive dusty galaxies are already present up to z∼ 7 (e.g., Strandet et al. 2017; Fudamoto et al. 2021). Several surveys have found that the obscured star formation density is already on the same order of magnitude as the unobscured one at z∼ 5 (e.g., Gruppioni et al. 2020; Zavala et al. 2021; Khusanova et al. 2021; Barrufet et al. 2023). Because of their high dust content, which absorbs a large fraction of the light produced by their stars and re-emits it in the rest-frame FIR, dusty galaxies are very faint in the optical (e.g., Wang et al. 2019; Talia et al. 2021), but very bright in the observed (sub)mm.

While the first surveys were performed at 850 μm by the submm common-user bolometer array (SCUBA, Holland et al. 1999), other instruments on single-dish telescopes were developed to probe longer wavelengths¹ such as Bolocam (Laurent et al. 2005), AzTEC (e.g., Scott et al. 2012), or MAMBO (Lindner et al. 2011) at 1.1/1.2 mm, and GISMO (Staguhn et al. 2014; Magnelli et al. 2019) at 2 mm. In addition, a population of very bright sources at 1.4 mm magnified by gravitational lensing was found by the south pole telescope (SPT, Vieira et al. 2013, see also Lestrade et al. 2010 about an earlier single detection with MAMBO). Finally, ALMA offers the capacity to perform interferometric surveys through large mosaics of pointings and blind photometric surveys have been performed at 1 mm (e.g., Dunlop et al. 2017; Franco et al. 2018) and 2 mm (e.g., Casey et al. 2021; Long et al. 2024). The survey performed at longer wavelengths found sources at a higher average redshift (e.g., Simpson et al. 2014; Weiß et al. 2013; Staguhn et al. 2014; Casey et al. 2021; Cowie et al. 2023). This effect can easily be modeled considering the combination of the different (sub)mm flux versus redshift relation at different wavelengths (Zavala et al. 2014), the strong evolution of the dust luminosity function, and the impact of lensing at high flux (Béthermin et al. 2015a).

Performing deep surveys at longer wavelengths may appear an effective strategy for maximizing the number of high-redshift detections. However, things are more complex than they seem. At 2 mm, dusty galaxies even up to z ∼ 7 are observed in the Rayleigh-Jeans regime with a flux density decreasing steeply with increasing wavelength. If the survey sensitivity is not significantly better than at shorter wavelengths, this could lead to fewer detections even at high redshift. However, surveys at 2 mm could be potentially powerful tools for detecting hypothetical rare dusty objects at very high redshifts. They could also be combined with shorter wavelength surveys to identify follow-up targets at high redshift (e.g., Casey et al. 2019; Cooper et al. 2022; Cowie et al. 2023) or with a very cold apparent dust temperature (e.g., Jin et al. 2022; Bing et al. 2024).

The new IRAM KIDS array (NIKA2, Monfardini et al. 2014; Calvo et al. 2016; Bourrion et al. 2016) was installed on the IRAM 30-meter telescope in October 2015 (Adam et al. 2018; Perotto et al. 2020) and can observe its circular 6.5 arcmin-diameter field of view simultaneously at 1.2 and 2 mm. Thanks to its combined high sensitivity and large field of view, NIKA2 is particularly efficient to perform blind surveys at these two wavelengths. The NIKA2 cosmological legacy survey (N2CLS) observed the GOODS-N and COSMOS fields. It provides both an unprecedented statistical sample of 2 mm sources and an unique opportunity to combine them with their 1.2 mm counterparts. First results based on data obtained before May 2021 are shown in Bing et al. (2023), which extensively discusses the number counts. At that time, the GOODS-N survey was complete and COSMOS was still missing a factor of 2 in integration time. Since then, the depth of COSMOS has increased, but this flux regime had already been probed by GOODS-N. Therefore, the new data do not significantly improve the constraints on the number counts. For the GOODS-N field, the detailed analysis of individual sources is presented in Berta et al. (2025) and Lagache et al. (2026), while the confusion due to undetected sources is covered in Ponthieu et al. (2026). A similar study was also performed on the serendipitous point-source detections in the NIKA2 large program targeting the Sunyaev-Zeldovich effect in galaxy clusters (LPSZ, Désert et al. in prep.).

In this paper, we present the final N2CLS 2 mm data, interpret their color and redshift distributions, and discuss the consequences for 2 mm surveys. In Sect. 2, we introduce the observed (N2CLS) and simulated data used in the paper. We describe the 2 mm catalog together with the methods used to produce it. In Sect. 3, we discuss the properties of the sources, such as their color and their redshift. We then discuss the relevance of 2 mm single-dish and interferometric photometric surveys in identifying the high-redshift populations of dusty star-forming galaxies in Sect. 4. Finally, we present our conclusions in Sect. 5.

2. Data

2.1. The NIKA2 cosmological legacy survey

The N2CLS is a 300 h guaranteed-time large program of the NIKA2 camera at the 30 m IRAM telescope (192-16, PIs: Lagache, Beelen, and Ponthieu). Thanks to its unprecedented mapping speed for a single-dish mm instrument, NIKA2 allowed us to map two popular fields with unprecedented depth, compared to other mm surveys of the same size in GOODS-N, and a similar depth, but a twice larger area than the Ex-MORA survey in COSMOS (Long et al. 2024). For both fields, the data reduction and map making were performed using the pointing and imaging in continuum (PIIC) data reduction pipeline developed and supported by IRAM (Zylka 2013; Berta & Zylka 2019). The beam FWHM is ∼12 arcsec and ∼18 arcsec at 1.2 mm and 2 mm, respectively.

The GOODS-N field covers 159 arcmin² with a 1σ depth of 0.17 mJy/beam and 0.047 mJy/beam at 1.2 mm and 2 mm, respectively. All the data were taken between October 2017 and May 2021. After scan selection (i.e. removing scans with higher noise linked to weather instabilities), the total time on source reached 78.0 h and 72.8 h at 1.2 and 2 mm, respectively. The data reduction is described in Bing et al. (2023). These data are close to the confusion limit (Ponthieu et al. 2026).

In the COSMOS field, we achieve a 1σ sensitivity of 0.3 mJy/beam at 1.2 mm and 0.09 mJy/beam at 2 mm, covering an area of 1130 arcmin². The data were taken between October 2017 and January 2023. After the scan selection, the total observing time reached 195 h. Bing et al. (2023) presented the analysis of the data obtained before May 2021 (representing on-source time of 78.7 h at 1.2 mm and and 79 h at 2 mm). We used the same method but with an updated version of PIIC to reduce the full dataset (Carvajal-Bohorquez et al. in prep. for the details).

The absolute calibration uncertainty coming from the planet models used to calibrate the NIKA2 data is 5%. The point-source root-mean-square (RMS) calibration uncertainties are about 6% at 1.2 mm and about 3% at 2 mm (Perotto et al. 2020).

2.2. The GOODS-N catalog

We used the 1.2 mm and 2 mm catalogs produced from NIKA2 and described in Bing et al. (2023). The sources were extracted from match-filtered² signal-to-noise (S/N) maps by applying a detection threshold. The source fluxes were then measured via point spread function (PSF) fitting photometry. End-to-end simulations (see Sect. 2.5) provided estimates of the purity (fraction of real sources in the extracted sample) as a function of this threshold. Bing et al. (2023) used a S/N threshold corresponding to 80% purity (3.0 at 1.2 mm and 2.9 at 2 mm). The biases impacting the photometry, such as flux boosting and filtering effects, were corrected using the median ratio between the intrinsic flux in the simulation and the output flux measured using the simulated map. The simulations were also used to compute the flux uncertainties, which are significantly larger than what we could expect from pure uncorrelated noise.

In this paper, we focused on 2 mm sources, selecting only high-reliability 2 mm sources corresponding to a sample purity of 95% (S/N ≥ 4.2). We obtained a subsample of 25 sources³ that we then matched with the full 1.2 mm catalog (S/N ≥ 3.0 corresponding to 80% purity) using a matching radius of 6.5 arcsec. Only one object had no 1.2 mm counterpart (N2GN_2_13)⁴. This source is close to the noisy edge of the field. We estimated a 5σ upper limit on its flux based on the flux uncertainties of sources located in areas with a similar depth. The catalog is given in Table B.1.

The counterparts of all these sources were identified using a large set of ancillary data and dedicated NOEMA follow up observations. This counterpart identification and the redshift determination is described in a companion paper (Berta et al. 2025). This procedure also shows that 4 of our 25 sources are multiple (N2GN_1_12, N2GN_1_17, N2GN_1_24, N2GN_1_34).

2.3. The COSMOS catalog

The observing time in the COSMOS field has been more than doubled since the source extraction of Bing et al. (2023). Thus, we produced new 1.2 mm and 2 mm catalogs using the full-depth maps, together with updated end-to-end simulations to estimate the sample purity and correct for the flux boosting and filtering effects. The source extraction and flux measurement were performed following the same approach as in GOODS-N. To estimate the purity, we compared the number of sources between end-to-end simulations and pure-noise maps to derive the fraction of spurious sources at a given S/N (see Carvajal-Bohorquez et al. for a full description of the method). We corrected the raw flux densities and estimated uncertainties using the same approach as in GOODS-N, with only minor changes. We estimate that a purity of 95% is reached for a S/N threshold of 4.6, which we then use to produce our secure 2 mm catalog containing 90 objects. The details will be given in Carvajal-Bohorquez et al. (in prep.).

We used a similar matching procedure as the one described for the GOODS-N field to associate the 2 mm source with their 1.2 mm counterparts. The 1.2 mm catalog used for this procedure was cut off at S/N ≥ 3.9, corresponding to 80% purity. We find a 1.2 mm counterpart for 82 of our 90 sources. Two of the 2 mm sources (N2CO_1_8 and N2CO_1_41) split into two 1.2 mm counterparts each thanks to the sharper beam at shorter wavelength. Two other sources (N2CO_2_29 and N2CO_2_61) are complex blends with a clear 1.2 mm emission near the 2 mm source. However, it is unclear what fraction of the 1.2 mm emission comes from the main 2 mm source. For instance, N2CO_2_61 is detected at 1.2 mm, but it is an intriguing case of a complex blend between the faint 1.2 mm emission associated with the radio source COSMOSVLADP_J095951.93+020542.6 and a brighter 1.2 mm dusty neighbor 8.8 arcsec away (Carvajal-Bohorquez et al., in prep.). The final catalog is provided in Tables B.2 and B.3.

To identify the multi-wavelength counterpart(s) corresponding to each N2CLS source, we first searched for a counterpart in ancillary high-resolution (sub)mm interferometric catalogs, which serve as a proxy for accurate source position (following Berta et al. 2025). We used the following datasets, in decreasing order of priority:

• AS2COSMOS (Simpson et al. 2020), which was given highest priority because of the dedicated follow-ups with ALMA of the bright 850 μm SMGs;

• A³COSMOS (Liu et al. 2019), a (sub)mm ALMA archival survey comprising a prior catalog (S/N ≥4.35) and a complementary blind catalog (corresponding to sources detected blindly with S/N ≥ 5.4 and not present in the prior catalog). We used the latest available version of the catalogs⁵ (2025/03/12).

• Radio JVLA 3 GHz catalog from Smolčić et al. (2017) if no ALMA counterpart was found.

We visually inspected all our sources to verify the proxy positions and to check for any missing counterparts. During this process, we identified three MeerKAT sources (MIGHTEE – MeerKAT International Gigahertz Tiered Extragalactic Exploration survey; Heywood et al. 2022). After pinpointing the precise location of the N2CLS sources using these high-resolution data, we searched for their short-wavelength counterparts to determine a spectroscopic redshift (when available) or a photometric redshift otherwise. We obtained eight spectroscopic redshifts from Chen et al. (2022), complemented by Mitsuhashi et al. (2021) from the brightest AS2COSMOS sources. We also collected 33 spectroscopic redshifts from the COSMOS spectroscopic redshift compilation by Khostovan et al. (2026), three from Jin et al. (2019), three from Brisbin et al. (2017), two from Daddi et al. (2021), one from Hasinger et al. (2018), one from Casey et al. (2015), one from Ikeda et al. (2026), one from Sillassen et al. (2025), one from Jin et al. (2022), one from Jin et al. (2024), and one from Jin (private communication). We also obtained 26 photometric redshifts from the COSMOS2020 catalog (Weaver et al. 2022), 17 from the COSMOS2025 (COSMOS-Web, Shuntov et al. 2025), three from Liu et al. (2019), one from Miettinen et al. (2017), and one from van der Vlugt et al. (2023). The full procedure will be described in Carvajal-Bohorquez et al. (in prep.).

We identified counterparts for 80 2-mm sources (with ten sources having no proxy detections), of which 18 have multiple counterparts. These 18 multiple sources break into 43 galaxies, including 19 with photometric redshifts and 24 with spectroscopic redshifts. Among the remaining 62 sources with a single counterpart, we have 29, 32, and one galaxies with photometric, spectroscopic, and no redshift, respectively. Finally, 11 galaxies were left without a redshift (one galaxy identified through a proxy but without redshift and ten 2-mm sources without proxy). Compared to GOODS-N, the redshift identification completeness is lower, which limits the interpretation of redshift-dependent statistics.

2.4. The SIDES simulations

The simulated IR dusty extragalactic sky (SIDES, Béthermin et al. 2017, 2022) is a semi-empirical realistic simulation of the sky from the mid-IR (MIR) to the mm. It is based on a dark-matter halo catalog from N-body numerical simulations projected into a light cone. The stellar mass of the galaxies hosted by the (sub)halos is derived using the abundance matching technique (see, e.g., Behroozi et al. 2010 or Moster et al. 2013 for a description of the method), which reproduces the observed stellar mass function by construction. The star-forming properties were then drawn following empirical scaling relations that included their scatter, such as the relation between the star formation rate (SFR) and the stellar mass (Schreiber et al. 2015). The galaxy population on this relation was tagged as main sequence and the positive outliers as starbursts.

A spectral energy distribution (SED) template was then attributed to each galaxy, depending on whether it is a main sequence or starbursting object. These templates are based on observations, following the observed evolution of the mean radiation field ⟨U⟩ with redshift (Magdis et al. 2012; Béthermin et al. 2015a). The ⟨U⟩ parameter is related to the dust temperature (approximately $T_{\mathrm{dust}}\propto {⟨U⟩}^{\frac{1}{4+\beta }}$, where β is the dust spectral index) and it is drawn by assuming an observationally motivated 0.2 dex scatter around the mean evolution. Finally, the SEDs were integrated into broad-band filters. To properly evaluate the impact of cosmic variance on our observations, we used the new 117 deg² version of SIDES introduced in Gkogkou et al. (2023, hereafter SIDES-Uchuu).

SIDES includes the free-free and the synchrotron emission of star-forming galaxies in their SED template, assuming the IR-radio correlation of Sargent et al. (2010) and an equality between synchrotron and free-free emission at 30 GHz. The radio galaxies were not included in SIDES, however, they are discussed in Sect. 3.4.

In Bing et al. (2023), we showed that SIDES reproduces well the number counts at both 1.2 and 2 mm from single-dish and interferometric surveys, if we take resolution effects into account. In this paper, we compare SIDES with measured colors and redshift distributions (see Sect. 3) to further validate it in the mm regime, before using it to discuss the relevance of 2 mm surveys (see Sect. 4).

2.5. End-to-end N2CLS simulation based on SIDES

The source fluxes measured by single-dish mm surveys are affected by a large set of effects such as blending, data filtering, correlated noise, or flux boosting. To simulate only the blending effect, we generated noiseless maps using a Gaussian beam for the sources and we extracted the sources to produce the so-called blob catalog.

To include the other effects, we produced end-to-end (hereafter, E2E) simulations based on SIDES and real N2CLS data. The GOODS-N simulations are described in Bing et al. (2023). For COSMOS, we generated a new simulation for the full dataset (see Carvajal-Bohorquez et al. in prep.) using a similar method. For each detector in the arrays, the incident flux is computed as a function of time using the SIDES simulated map and the real pointing coordinates, which vary with time. The resulting detector timelines are then combined with the real detector timelines after reversing the sign of half of them. These new timelines are used to produce simulated maps using PIIC, in which the true sky signal cancels out while the SIDES signal and the true noise remain. We then run the same source extraction procedure, as described in Sects. 2.2 and 2.3. A similar procedure was applied to each of the 117 SIDES-Uchuu 1-deg² tiles. For the purposes of this paper, we matched the simulated 2 mm catalog with the 1.2 mm one, using the same procedure as the one applied to the observational catalogs. The result is hereafter called the E2E catalog.

We checked that the number of robust 2 mm sources per simulated COSMOS field (71 ± 16) is compatible at the 1.2-σ level with the N2CLS (90). Bing et al. (2023) also found marginally higher number counts in their shallower COSMOS data than in GOODS-N and SIDES E2E. A similar result was found by Adscheid et al. (2024) using the A3COSMOS archival survey.

3. Colors and redshifts of the 2 mm sources

3.1. 2 mm versus 1.2 mm colors

Thanks to its simultaneous observations at 1.2 mm and 2 mm, the N2CLS is ideal to study the mm colors of extragalactic sources. We focus our analysis on the COSMOS field for which we have a significantly larger number of detections. Since the color is a ratio, it is affected by the uncertainties on both the numerator and the denominator. Thus, we limited our analysis to high-reliability (S/N ≥4.6) sources at both 1.2 mm and 2 mm. However, in COSMOS, this choice has no impact compared to a S/N ≥3.9 cut at 1.2 mm (80% purity), since there are no robust 2 mm sources with an S/N between 3.9 and 4.6 at 1.2 mm. The measured 2 mm versus 1.2 mm color distribution is shown in Fig. 1. This distribution is asymmetric with a long tail of high color values, corresponding to red objects.

Fig. 1. Distribution of the 2–1.2 mm color for the true N2CLS COSMOS catalog (red filled histogram, only S/N ≥4.6 sources at both wavelengths and various simulated catalogs based on SIDES. The solid black histogram corresponds to the simulated SIDES galaxy catalog after applying a flux cut similar to COSMOS data (see Sect. 3.1). The yellow histogram is based on a similar selection, but applied to the blob catalog extracted from the noiseless simulated map. The grey histogram is the distribution obtained using the E2E simulation described in Sect. 2.5. The error bars on the N2CLS histogram are computed assuming a Poisson law, and those from the E2E are the standard deviation between the various simulated fields. The histograms are normalized to have a unity area (∫n(c) dc = 1, where c is the color). The lower limits on the color of 1.2 mm non-detections are shown as right-pointing triangles. Their y-axis position is arbitrary and was chosen to reduce overcrowding in the figure.

To interpret these data, we computed several distributions based on the SIDES simulation. In this work, we used the 117 deg² version of the simulation (Gkogkou et al. 2023) up to z∼ 7 to estimate the field-to-field variance. We estimated the impact of this redshift limit using the older 2 deg² version with redshifts up to z = 10 (Béthermin et al. 2017), and found that it is negligible. For instance, the mean color varies by only 0.3%. For the individual galaxy color distribution, we just cut the SIDES simulated galaxy catalog at the typical flux cut of the N2CLS catalog. For this purpose, we use the S/N limit of 4.6 and multiply it by the mean 1-σ depth of our maps, and obtained the following selection: S_2 mm ≥ 0.414 mJy and S_1.2 mm ≥ 1.38 mJy. To evaluate the impact of the NIKA2 angular resolution, we applied similar cuts to the blob catalog presented in Sect. 2.5. Finally, we use the catalogs resulting from our E2E simulations (including real noise in the simulated timelines, the instrumental beam, and map-making effects), which were also introduced in Sect. 2.5 to compare accurately SIDES to N2CLS. Since these simulations produced 117 independent fields with the same field geometry as N2CLS COSMOS, we could use them to estimate the field-to-field variance. These uncertainties are larger than the ones obtained from a bootstrap method, since galaxies in a given region of the sky tend to accumulate at similar redshift and thus have similar colors. These correlations between sources increase the variations between fields of the various summary statistics. The three distributions obtained from SIDES are presented in Fig. 1.

The means and medians of all these distributions are listed in Table 1, and are roughly compatible. The SIDES galaxy catalog have a slightly lower mean but the same median as N2CLS, while the blob and E2E catalogs have higher means and medians. This shows that the resolution effects have a small impact (∼10%) on the average colors. To evaluate the quality of SIDES predictions, the most relevant comparison is between the N2CLS and the E2E, since they should be affected by exactly the same instrumental and source extraction effects. The difference between the two means and medians corresponds to 1.2 and 1.7 times the field-to-field standard deviation. We also have to take into account the calibration uncertainties. Because the 1.2 and 2 mm fluxes are affected in the same way (i.e. systematic uncertainty), the absolute calibration uncertainty does not have to be taken into account for the color uncertainty. On the other hand, the point-source rms calibration uncertainties (6% at 1.2 mm and about 3% at 2 mm, Perotto et al. 2020), lead to a 6.7% statistical uncertainty on the color. By quadratically combining these uncertainties with the field-to-field variance, we obtained uncertainties of 0.020 and 0.018 on the mean and the median color, corresponding for both to a 0.8 and 1.1σ difference. Thus, there is no significant tension observed with respect to the average colors of the N2CLS sources and the SIDES predictions.

However, the distributions showed in Fig. 1 exhibit some difference of shapes. The most striking one is how narrow is the SIDES galaxy color distribution compared to the others with a standard deviation which is ∼3 times smaller. This is not the case for the SIDES blob catalog, suggesting that the resolution effects drive the dispersion in color rather than the noise. In contrast, the N2CLS and the SIDES E2E catalogs agree on the dispersion at the 0.5σ level, demonstrating how efficiently our E2E simulation can predict the observations.

Despite the agreement between N2CLS and SIDES E2E on these basic metrics, the distributions are not exactly similar. For instance, the probabilities returned by the Kolmogorov-Smirnoff (KS) null-hypothesis test comparing N2CLS and SIDES are always lower than 0.01. Above a p-value of 0.05, we can consider that the two distributions are drawn from the same probability density function. In contrast, a low probability indicates that the two distributions are not similar. This test does not take into account the field-to-field variance and the calibration uncertainties. It may thus be improper and must be interpreted with caution. However, there is a 2.5-σ excess in the 0.16–0.18 bin after taking into account the field-to-field variance. If we exclude this color bin from the KS test, we end up obtaining a p-value of 0.23, confirming that the disagreement is caused by this single bin. The origin of this peak in the color distribution remains unclear.

Finally, we also checked the mean color in the GOODS-N field. This field is limited by its small number of 2 mm detection (25), but it is deeper. We find a mean color of 0.196 ± 0.014, which is compatible with the 0.208 value found in the SIDES galaxy catalog after applying flux cuts similar to GOODS-N (S_2 mm ≥ 0.2 mJy and S_1.2 mm ≥ 0.71 mJy. Contrary to COSMOS, we do not observe any excess in the 0.16–0.18 color bin, suggesting this feature is a field variance effect.

3.2. Redshift distribution

The redshift distribution is an important constraint for models. For a flux-limited sample selected at a single wavelength, this distribution varies with both the flux cut and the wavelength (e.g., Smolčić et al. 2012; Zavala et al. 2014; Béthermin et al. 2015b). In the (sub)mm, the mean redshift is expected to increase with the wavelength. With its combination of size and depth, the N2CLS 2 mm survey is ideal to test model predictions in the long-wavelength and faint-flux regime. Since it has a complete redshift association, we first focused our analysis on the GOODS-N field where all NIKA2 sources have at least one associated galaxy with a known redshift (Berta et al. 2025).

In Fig. 2, we present the redshift distribution of 2 mm high-reliability sources with Poisson uncertainties. In the case of a multi-component source, we use the redshift of the main component (listed as “a” in Berta et al. 2025)⁶. We compared these results with the predictions from SIDES. The source identification is a complex process involving multi-wavelength data.

Fig. 2. Redshift distribution of the GOODS-N high-reliability 2 mm sources. The red histogram is the distribution measured in the N2CLS. When sources are multiple, the redshift of the main component is used (see Sect. 3.4). The uncertainties correspond to a Poisson law. The grey histogram is derived from the SIDES galaxy catalog after applying a flux cut similar to the GOODS-N field (see Sect. 3.2). The uncertainties correspond to the field-to-field variance of a 15.9 × 10 arcmin2 field. The yellow histogram is based on a similar selection, but applied to the SIDES blob catalog and using the redshift of the brightest galaxy in the beam. The histograms are normalized to have a unity area (∫n(z) dz = 1).

Since SIDES does not produce simulated optical and radio images, we had to rely on two simplified approaches. First, we selected a sample in SIDES by cutting the galaxy catalog above the survey flux limit (S_2 mm ≥ S/N_lim × σ = 0.2 mJy). This sharp flux limit is a simplification, since the survey depth varies across the field and the transfer function can also alter the true flux limit. However, the redshift distribution is only weakly dependent on the exact flux cut (Béthermin et al. 2015b), so an approximate flux limit is sufficient for our purpose. Our second approach takes resolution effects into accounts and is based on the SIDES blob catalog. There are always several galaxies with different redshifts present in the same NIKA2 beam, and we attribute the redshift of the brightest galaxy at 2 mm to the blob (see Béthermin et al. 2017, 2024 for a detailed presentation of the method). For both approaches, we compute the field-to-field variance by dividing the SIDES catalog into 15.9 × 10 arcmin² subfields (2541 in total). The galaxy (gray histogram) and blob (yellow histogram) are very similar (see Figs. 2 and 3). This shows that, contrary to the number counts, the redshift distribution is not sensitive to the resolution effects (see, e.g., Bing et al. 2023).

Fig. 3. Same as Fig. 2, but for the COSMOS field. The blue dashed histogram represents the Ex-MORA survey (Long et al. 2024).

There is an overlap between the 1σ error bars of the N2CLS/GOODS-N and both galaxy and blob SIDES redshift distributions in all the redshift bins except the 5.0 < z < 5.5 one, where N2CLS is ∼1.5σ above the SIDES predictions. This is due to a known galaxy overdensity at z∼ 5.2 around the HDF850.1 galaxy. A companion N2CLS paper shows that this overdensity contain numerous dusty star-forming galaxies (Lagache et al. 2026). There is also an excess in the 4.0 < z < 4.5 bin, but the 1-σ error bars overlap and, thus, it is not significant. We measure a mean redshift of 3.6 ± 0.3, while the SIDES simulation predicts 2.9 ± 0.2 for both the galaxy and blob selection. The redshift limit at z∼ 7 has a minor impact on this quantity and we found a mean redshift of 3.0 using the 2 deg² version of SIDES which is cut at z = 10 (Béthermin et al. 2017)⁷. This ∼2σ difference is likely due to the overdensity at z∼ 5.2.

To test this, we performed a KS test on both redshift distributions after excluding the 5.0 < z < 5.5 range. The test returns a p-value of 0.4 (same values for the SIDES galaxy and blob catalogs), and N2CLS in GOODS-N and SIDES are thus fully compatible after excluding this known overdensity, while it is only marginally compatible (p-value of 0.06 for both the galaxy and blob catalogs) if we include this redshift bin.

The COSMOS 2 mm catalog is slightly incomplete in redshift (only 88% have a main component with an associated redshift), but has a significantly larger number of detections. We thus performed a similar analysis as in GOODS-N, but using a S_2 mm flux cut of 0.414 mJy for the SIDES selection and a 34 × 34 arcmin² for the field-to-field variance. This field size matches the 1130 deg² area where sources are extracted in the real N2CLS data (Carvajal-Bohorquez et al. in prep.). However, the results should be taken with caution because of the mild redshift incompleteness. The result is shown in Fig. 3. We obtained a mean N2CLS redshift of 3.2 ± 0.2, which agrees with SIDES galaxy prediction of 3.2 ± 0.2 and the blob prediction of 3.1 ± 0.1. The 1σ error bars of the N2CLS COSMOS and the SIDES distribution overlap, and the KS test returns a p-value of 0.58 for the galaxy catalog and 0.47 for the blob catalog, compatible with the two samples being drawn from the same distribution.

We also compared our N2CLS results with those of the Ex-MORA survey (Long et al. 2024), which has a similar depth and a partial overlap. The redshift distributions of the two surveys agree all the 1-σ error bars overlapping except in the 1< z< 2 bin, where this is almost the case. The KS test p-value is 0.2, confirming the agreement between the two surveys.

Contrary to GOODS-N, there is no tension between the observed and predicted mean redshifts. This suggests that GOODS-N is a peculiar field with a strong overdensity at z∼ 5.2. However, COSMOS is slightly incomplete, and could be more incomplete at higher redshift, biasing the mean redshift toward lower values.

3.3. Colors versus redshift

Since the dust emission of galaxies peaks at ∼80 μm rest-frame, red colors in the (sub)mm can be used to find high-redshift candidates. For instance, Herschel surveys identified sources with red Herschel/SPIRE colors (usually 250 versus 500 μm) as z≳ 4 candidates (Dowell et al. 2014; Asboth et al. 2016; Donevski et al. 2018). However, Béthermin et al. (2017) and Donevski et al. (2018) showed using simulations that a large fraction of these sources have redder observed colors than their intrinsic ones because of fluctuations of the noise in the bands used to compute the color (positive in the long-wavelength band and negative in the short-wavelength one). Similarly, mm colors were proposed to identify candidates at even higher redshifts (e.g. Casey 2020; Cooper et al. 2022; Cowie et al. 2023). In this section, we study the (sub)mm color of N2CLS galaxies versus redshift to test this approach.

3.3.1. Colors versus redshift in the N2CLS

In the color-versus-redshift relation, we can only place sources which have all their flux coming from the same redshift. Thus, we selected only galaxies with a single counterpart to have no ambiguity. We derived both the N2CLS colors (S_2 mm/S_1.2 mm) and the 2 mm versus 850 μm colors after matching our catalogs with the S2COSMOS (Simpson et al. 2017) and SUPER GOODS (Cowie et al. 2017) SCUBA2 catalogs in COSMOS and GOODS-N, respectively. Contrary to Sect. 3.1, we use the full 1.2 mm catalog to compute the colors. The results are presented in Fig. 4.

Fig. 4. Color as a function of the redshift in the COSMOS (top) and GOODS-N (bottom) fields. The left plots show the internal NIKA2 color (S2 mm/S1.2 mm) and the right ones show the ratio between the N2CLS 2 mm and the SCUBA2 850 μm fluxes from S2COSMOS (Simpson et al. 2017) and SUPER GOODS (Cowie et al. 2017). The filled blue circle are the N2CLS sources with spectroscopic redshifts, while the open red circles corresponds to photometric redshifts. The grey dots are from the SIDES 2 deg2 simulated catalog after applying a 2 mm flux selection similar to SIDES (see Sect. 3.2). The yellow solid line corresponds to the median of the SIDES blob catalog, using the redshift of the brightest galaxy in the beam, and the colored area correspond to the 1σ confidence region. The dashed line is the mean color as a function of redshift from SIDES 2 deg2 without applying any selection. We also show the sensitivity ratio (σ2 mm/σ1.2 mm where σ2 mm and σ1.2 mm are the RMS of the noise at 2 mm and 1.2 mm, respectively) for ALMA (dotted line) and N2CLS (dot-dashed line). See the discussion in Sect. 4.1 for more details. The right-side panel of each plot show the color distribution of sources with a photometric redshift (red), a spectroscopic redshift (blue), and no redshift (green).

There is no clear trend of the observed N2CLS S_2 mm/S_1.2 mm color versus redshift (Fig. 4, left) and the data points exhibit a large scatter. For a linear fit of the color versus redshift, we find a slope of 0.009 ± 0.005 when combining GOODS-N and COSMOS. The color measurements are affected by the uncertainties in both the numerator and the denominator. Since we accounted for both the limited S/N and the uncertainties in the deboosting, this can lead to large relative uncertainties (> 1/3); thus, the N2CLS sample is too small to detect any significant trend.

The observed evolution of the 2 mm versus 850 μm color with redshift exhibits a tentative mildly increasing trend (Fig. 4, left) with a slope of 0.014 ± 0.006 in COSMOS, corresponding to a 2.3σ significance. This result agrees with Cooper et al. (2022), Cowie et al. (2023), and Long et al. (2024), who also found a weak trend without indicating its significance. We attempted to fit separately the trend below and above z = 2.5 (motivated by the trend predicted by SIDES in Sect. 3.3.2). The results are not statistically significant (0.002 ± 0.026 at z< 2.5 and 0.017 ± 0.011 at z> 2.5), but they are compatible with the flat trend below z = 2.5 and the rise above predicted by SIDES.

In the COSMOS field, a small fraction of the sources has no redshift (12%). We compared the color distributions of the sources depending of the type of information on their redshift (see the right side panels of each subfigure of Fig. 4). For the N2CLS S_2 mm/S_1.2 mm colors in COSMOS, the three distributions are similar, except a mild tension between the sources with a spectroscopic redshift and without redshift in COSMOS (p-value of 0.04 to have similar parent distributions according to the KS test). The redshift incompleteness of our sample should thus not affect significantly our results. Concerning the S_2 mm/S_850 μm colors in COSMOS, there is a tension between both the photometric and spectroscopic redshift samples and the objects without redshift (p-value of 0.005 and 0.003, respectively). There is no tension between the other distributions. The sources without redshift tend to be redder. It is thus possible that the sample of 11 COSMOS 2 mm galaxies without redshift (no proxy or no redshift available for the counterpart) is biased toward cold or high-z sources.

3.3.2. Colors versus redshift in SIDES

We compared our results with the intrinsic colors of the galaxies in SIDES after applying the same 2 mm flux cut as in Sect. 3.1. We used the 2 deg² version of the simulation (Béthermin et al. 2017) for this analysis, since it goes up to z∼ 10 and this makes it wide enough for this qualitative comparison. The results are shown in Fig. 4 as grey filled circles. We remark that the intrinsic scatter predicted by SIDES is much smaller than the observational scatter, which is driven by measurement uncertainties. We also computed the mean color of the full SIDES galaxy catalog without applying any flux selection (dashed line). Finally, we derived the evolution of the color in the SIDES blob catalog as a function of the redshift of the brightest component (yellow line and shaded area).

The SIDES simulation includes only a scatter on the mean intensity of the radiation field ⟨U⟩, which is directly connected to the dust temperature. In contrast, the slopes of the SEDs in the Rayleigh-Jeans regime have only small variations. Since most of the sources are observed close from this regime, the S_2 mm/S_1.2 mm color at fixed redshift has thus only a very small scatter. We could thus think that SIDES may severely under-predict the color dispersion. However, as shown in Sect. 3.1 and Fig. 1, the observed color scatter is 1-σ compatible with the SIDES scatter after applying all the observational effects. Thus, there is no evidence that the SIDES color scatter is underestimated.

The mean NIKA2 color in both SIDES blob and galaxy catalogs is nearly constant (∼0.2) between z = 0.5 and 4. It increases slightly at higher redshifts (∼0.3 at z = 8). The colors obtained for the 2 mm-selected SIDES galaxy (gray symbols) are redder (higher values) than the mean value of all individual galaxies in the SIDES simulation (dashed line). This selection bias will be discussed in Sect. 4.2. The color slightly increases at z< 0.5 (as shown by the dashed line for the SIDES mean and the gray dots for individual 2 mm-selected SIDES galaxies) due to the free-free and synchrotron contributions included in the SIDES templates (Béthermin et al. 2012). Additionally, we note that the blob selection produces slightly redder colors at low redshifts, compared to the 2 mm-selected galaxy catalog.

The results are similar for the S_2 mm/S_850 μm color (Fig. 4, right). However, the reddening of the colors with increasing redshift starts at lower redshift (z∼ 2.5). This is expected since the Rayleigh-Jeans power-law regime starts above 200 μm rest-frame. Below this rest-frame wavelength, the SED is flatter and this regime is thus reached at lower redshift for the observed 850 μm than the 1.2 mm. The blob catalog produces marginally redder colors at low z and bluer ones at z> 6 than the 2 mm-selected galaxy catalog. Cowie et al. (2023) proposed a power-law relation between the logarithm of the color and the redshift (log₁₀(S_2 mm/S_850 μm) = 0.055z − 1.287), which disagrees with SIDES forecast at z < 2 and z > 5. Unfortunately, we lack data in these redshift ranges to test whether SIDES is correct or if the power law is a good approximation.

3.4. 2 mm sources without 1.2 mm counterparts

Because of the expected mild evolution of the NIKA2 color with redshift at z≳ 4, the sources detected only at 2 mm could potentially turn out to be at very high redshift. However, our sample has a 95% purity and these sources could thus be spurious.

In the GOODS-N field, only one source (N2GN_2_13) is detected only at 2 mm. As explained in Berta et al. (2025), this source was detected by NOEMA and has a short-wavelength counterpart with a photometric redshift of 2.071. This suggests that their is no spurious source in the GOODS-N sample. Based on the 25 detections and the 5% spurious rate, we expect in average 1.25 spurious sources per GOODS-N-sized field, and the Poisson probability of having 0 spurious objects in such a field is 0.29. This absence of spurious source is thus compatible with our estimated sample purity.

In the COSMOS field, out of the 90 high-reliability 2 mm sources, we found 8 sources without 1.2 mm counterpart (S/N ≥ 3.9 cut). The two brightest 2 mm-only sources (N2CO_2_37 and N2CO_2_44) are associated (matching distance of 1.40 and 1.42 arcsec, respectively) with two radio sources (COSMOSVLADP_J100114.85+020208.8 and COSMOSVLADP_J095945.19+023439.3, respectively) with spectroscopic redshifts of 0.97 and 0.12, respectively. N2CO_2_37 has a 4 GHz flux density of 7.5 mJy (Schinnerer et al. 2010). Extrapolating this emission at 150 GHz (2 mm) by assuming a standard synchrotron spectrum in ν^−0.8, we obtained 0.41 mJy (and 0.59 mJy for ν^−0.7). The majority of its 2 mm N2CLS flux (0.71 mJy) is thus expected to come from synchrotron emission. For the second source, which is more than one order of magnitude fainter in radio, the synchrotron is less likely to contribute, but we cannot exclude a contamination of dust emission by some free-free emission.

To understand the nature of the six remaining sources, we used the end-to-end simulation introduced in Sect. 2.5 to estimate the expected number of S/N ≥ 4.6 2 mm sources without S/N ≥3.9 1.2 mm counterpart, and found 5.8 ± 3.0 objects. The observations are thus fully compatible with our simulation. We also estimated the probability using another approach by computing the probability to find 6 or more spurious sources assuming a mean expected value of 4.5 (5% spurious rate times 90 sources found in the N2CLS), and found a probability of 0.29. It is thus statistically possible that these 6 sources are all spurious. Searching for these 2 mm-only sources is thus not a reliable way to select very high-z candidates.

3.5. Comparison with the Ex-MORA survey

In the COSMOS field, the extended mapping obscuration to reionization with ALMA (Ex-MORA) survey (Long et al. 2024, see also Casey et al. 2021 and Zavala et al. 2021 for the first MORA survey) has produced an interferometric mosaic of 577 arcmin² at 2 mm, overlapping with the N2CLS footprint⁸ on 533 arcmin². In this section, we compare the Ex-MORA and N2CLS catalogs. The Ex-MORA survey is slightly deeper than N2CLS with flux limit of ∼0.3 mJy versus ∼0.4 mJy. Out of the 37 Ex-MORA sources, four of them are on the noisy edge of the N2CLS field and thus outside of the region where sources were extracted (eMORA.18, eMORA.30, eMORA.31, eMORA.36). Of the 33 remaining Ex-MORA sources, 22 have a corresponding N2CLS 2 mm detection in the full catalog with S/N ≥ 3.9⁹ For most of the sources without counterpart, we can observe a faint signal (1.5 ≤ S/N ≤ 3.9) in the N2CLS map at the Ex-MORA position: S/N = 3.4 for eMORA.5, S/N = 3.3 for eMORA.6, S/N = 1.5 for eMORA.9, S/N = 1.9 for eMORA.17, S/N = 2.8 for eMORA.19, S/N = 2.1 for eMORA.23, S/N = 3.6 for eMORA.28, S/N = 2.9 for eMORA.29, S/N = 2.9 for eMORA.32, S/N = 2.5 for eMORA.33. Only the source position of eMORA.34 has no positive fluctuation of the signal (S/N = − 0.7), and is among the lowest S/N of the Ex-MORA survey (5.1).

We also investigated the Ex-MORA counterparts of robust (S/N ≥ 4.6) N2CLS 2 mm sources. Of the 54 N2CLS sources within the Ex-MORA footprint, only 18 have corresponding Ex-MORA detections within the N2CLS beam (< 9 arcsec from the N2CLS centroid). In addition, we inspected the Ex-MORA maps and found at least one candidate source with an S/N > 4 signal in the beam of 13 other N2CLS sources. Four other sources are located in a noisy region at the edge of the Ex-MORA survey, where no detection is possible. Thus, we are left with 19 sources without any significant counterparts in the NIKA2 beam. Only four of them are 2 mm-only detections (see Sect. 3.4) and are potentially spurious. The 15 others are also detected at 1.2 mm and are likely true objects missed by Ex-MORA.

In Fig. 5, we show the comparison between the N2CLS and the Ex-MORA 2 mm flux of objects detected by both surveys. Most of the low-flux source (0.3–0.6 mJy) have compatible fluxes at 2σ, but we can observe a small average excess of N2CLS flux compared to Ex-MORA. This excess is striking at higher flux, where the relative uncertainties are much smaller. A similar excess has been found in GOODS-N for four 2 mm sources followed up by NOEMA (Berta et al. 2025). For sources with a single known component, the N2CLS fluxes have a mean excess ratio of 1.3, and it reaches 2.2 for multiple sources. In the case of multiple sources, the flux difference is straightforward to explain. Indeed, for our multiple sources, Ex-MORA is detecting only one of the counterparts, leading to a significantly smaller Ex-MORA flux than the total flux from the galaxies in the NIKA2 beam.

Fig. 5. Comparison between the exMORA (x axis, Long et al. 2024) and N2CLS (y axis) 2 mm flux. The N2CLS sources with a single counterpart are in blue and the multiple sources are in red. The black dashed line indicates the one-to-one relation and the dotted line includes the various corrections (bandpass, source blending, and source sizes) discussed in Sect. 3.5.

The origin of the non-unity flux ratio between Ex-MORA and single-component N2CLS sources is less clear. It could be due to the combination of several effects. A first systematic factor comes from the large bandpass of NIKA2 and its different central observed wavelengths compared to Ex-MORA. We compute the flux density ratio between the two instruments assuming a power-law spectrum in λ^−2.95 (or ν^2.95) based on the measured mean color between 1.2 mm and 2 mm (Sect. 3.1 and Table 1):

$$\begin{array}{c}\hfill \frac{S_{2\mathrm{mm}}^{\mathrm{NIKA}2}}{S_{2\mathrm{mm}}^{\mathrm{Ex}-\mathrm{MORA}}}=\frac{\int {\left(\frac{\nu }{{\nu }_{\mathrm{Ex}-\mathrm{MORA}}}\right)}^{2.95}{t}_{\nu }d\nu }{\int t_{\nu }d\nu }=1.085,\end{array}$$(1)

where the central frequency of Ex-MORA is ν_{Ex − MORA} = 147 GHz and t_ν is the transmission of the NIKA2 2 mm bandpass as a function of frequency. This ratio increase to 1.169 when assuming a modified black body in the Rayleigh-Jeans regime and with an emissivity index β of 2.4, as found by Cooper et al. (2022).

In addition, the flux measured in single-dish blobs can be significantly larger than the flux of the brightest galaxy in the beam (Hayward et al. 2013; Karim et al. 2013; Cowley et al. 2015; Scudder et al. 2016; Béthermin et al. 2017; Bing et al. 2023) due to the presence of other galaxies in the beam. To evaluate the impact of this effect, we compared the flux in the SIDES blob catalog (Sect. 2.5) with the flux of the brightest galaxy in the beam (see Béthermin et al. 2024 for a detailed description of the method), and found a median ratio of 1.14.

Finally, the Ex-MORA flux is measured at the peak of emission, but this method can slightly underestimate the flux if the sources are marginally resolved (e.g., Sect. 3.5 of Béthermin et al. 2020). Assuming Gaussian profiles for both the beam and the galaxy, the peak flux underestimates the total flux by the ratio between the synthesized beam area (Ω_beam ∝ a_beam b_beam, where a_beam and b_beam are the beam FWHM of the major and minor axes) and the area of its convolution with the galaxy profile (${\mathrm{\Omega }}_{\mathrm{beam}}\propto \sqrt{(a_{\mathrm{beam}}^2+s_{\mathrm{intr}}^2)(b_{\mathrm{beam}}^2+s_{\mathrm{intr}}^2)}$, where s_intr is the intrinsic FWHM of the galaxy). For a typical Ex-MORA beam size of 1.68 arcsec × 1.44 arcsec, this leads to an underestimate by a factor of 0.97, 0.91, and 0.71 for a source intrinsic FWHM of 0.25, 0.5, and 1 arcsec, respectively. Thus, this effect is non-negligible even if the galaxies are three times more compact than the beam. The size of the sources can vary from a radius of ∼0.5 kpc to ∼3 kpc (Gómez-Guijarro et al. 2022; Pozzi et al. 2024; Hodge et al. 2025), corresponding to a FWHM from ∼0.1 arcsec to ∼1 arcsec for galaxies between z = 1 to z = 5. Depending on the source, the effect can range from negligible to strong.

The dotted line in Fig. 5 shows the expected relation between the N2CLS and Ex-MORA flux after applying all these corrections for a source FWHM of 0.5 arcsec and a ν^−2.95 reference spectrum. When all these effects are included, there is no longer a discrepancy between the two surveys. The comparison between single-dish and interferometric surveys is thus not straightforward. At similar depth, interferometric surveys will tend to miss extended and multi-component objects, while single-dish surveys will be biased toward groups of galaxies or random alignments. The significant number of N2CLS sources without an Ex-MORA counterpart suggests that these extended and composite sources are not rare. In addition, depending on the exact noise realization in each survey, some sources close to the detection limit may be detected by only one survey. This emphasizes the need for future combined analyses of surveys of the same area at different resolutions.

4. Relevance of 2 mm continuum surveys in searching for high-z DSFGs

The colors and redshifts of the N2CLS 2 mm sources agree with the predictions from the SIDES model (Sects. 3.1, 3.2, and 3.3). In addition, the 2 mm-only sources (only 8 out of 115 N2CLS 2 mm sources for both GOODS-N and COSMOS) are either z≲ 1 radio sources or are compatible with the expected spurious detection numbers. This tends to rule out the possibility of an exotic population of numerous galaxies detected only at 2 mm at high z or with very cold dust. A similar conclusion is reached by Ponthieu et al. (2026) based the measured N2CLS confusion noise, which is similar to the SIDES predictions. This is also consistent the results of the MORA and Ex-MORA surveys (Casey et al. 2021; Zavala et al. 2021; Long et al. 2024), which also found no excess of exotic very dusty sources. These new data, combined with our model, can also provide insights on the scientific performance (e.g., number of detections, coverage of the SFR-z plane) expected from 2 mm surveys and help us to better understand the impact of the flux selection at this wavelength.

4.1. Based on N2CLS observations

In this section, we discuss the merits of surveys at different wavelengths for detecting galaxies at various redshifts. Here, we consider a wavelength to be more efficient than another if it can detect a source with a higher S/N, or equivalently with less telescope time. Surveys can have different integration times and areas. It is thus important to find a common metric for all of them. In this discussion, we first focus on the case of comparing surveys in different bands (whether simultaneous or not) with the same area Ω_survey and total time t_survey. Since the sensitivity σ in both bands scales as $t_{\mathrm{survey}}^{-1/2}$ at fixed area and as Ω_survey^−1/2 at fixed total time, the sensitivity ratio between the two bands will be invariant as long as the two bands have the same area Ω_survey and total time t_survey. Sources with a color S_A/S_B equal to the sensitivity ratio σ_A/σ_B will have the same expected S/N in both bands (S_A/σ_A = S_B/σ_B). Sources with colors larger (lower, respectively) than the sensitivity ratio will be detected with a higher S/N in band A (band B, respectively). If most of the sources colors are above (below, respectively) the sensitivity ratio, the survey will be more efficient in band A (band B, respectively).

In Fig. 4, we added the N2CLS sensitivity ratio between 2 mm and 1.2 mm (left panels, dot-dash lines) to enable a comparison with the colors of N2CLS sources. We also computed the sensitivity ratio of ALMA using the observing tool¹⁰. We estimated the expected sensitivity for 1.5 × 1.5 arcmin² continuum mosaic at 150 GHz (2 mm, band 4) and 250 GHz (1.2 mm, band 6) with similar total observing times (5.85 h, 0.02 mJy RMS at 2 mm and 0.044 mJy RMS at 1.2 mm). The chosen size and the total time used in our computation are not important for estimating the sensitivity ratio, if we neglect the discontinuities in the total survey time generated by the need for new calibration cycles when going deeper. In contrast, this method takes into account that more pointings are necessary at higher frequency to cover the same area due to the smaller primary beam. For NIKA2, the 1.2 mm is more efficient than the 2 mm for most of the sources across the full redshift range covered by the N2CLS (z≲ 6.5). For ALMA, the 1.2 mm is consistently more efficient. A similar comparison can be made between the 2 mm and 850 μm bands. In the case of ALMA, the 850 μm (353 GHz, band; 0.16 mJy RMS on 1.5 × 1.5 arcmin² observed in 5.85 h) observations are overall more efficient than those at 2 mm (band 4, 150 GHz). However, the weak trend of the color versus redshift (Sect. 3.3) suggests the possibility of a turnover in favor of the 2 mm band at high redshift. We discuss this possibility based on the SIDES simulation in Sect. 4.2.

The comparison of surveys that use different observational strategies (e.g., total observing time or survey area) and different instruments is less direct. For example, for a fixed total observing time, a narrower and deeper survey may appear more efficient. Therefore, comparisons between different instruments (e.g., SCUBA2 and NIKA2) should be interpreted with caution. N2CLS 2 mm could be considered more efficient than S2CLS 850 μm (Simpson et al. 2017) at all redshifts, although this comes with the caveat that S2CLS covers a larger area. In contrast, SUPER GOODS 850 μm (Cowie et al. 2017) is more efficient than the 2 mm band of the N2CLS at z≲ 4, but the trend is weak and the turnover redshift is highly uncertain. Overall, 1.2 mm emerges as the most effective wavelength for conducting a census of dusty galaxies in deep blank cosmological surveys.

4.2. Based on SIDES

We can also use the SIDES simulation to estimate the most efficient band across the redshift range. The 2 mm band becomes more efficient than the 1.2 mm band only at z≳ 7.5 for NIKA2, and never in the case of ALMA. At these very high redshifts, we do not expect a significant population of heavily dust-obscured sources to be formed (e.g., Michałowski 2015; Burgarella et al. 2020). Concerning the comparison between 2 mm and 850 μm, we confirm that the N2CLS 2 mm band is more efficient than the SUPER GOODS 850 μm data at z≳ 4, and than the S2COSMOS 850 μm data at all redshifts. For ALMA, SIDES predicts that the turnover in favor of the 2 mm occurs around z = 5.5.

Since SIDES is able to reproduce the current 2 mm observations (see also Bing et al. 2023 for the number counts), we can use this model to better understand the impact of a 2 mm selection. In Fig. 6, we show the predicted redshift distribution of galaxies selected using similar flux limits as N2CLS in GOODS-N (see Sect. 3.2). Contrary to previous observational redshift distributions, these distributions are not normalized, and indicate the total number of detections. This allows us to compare directly the number of expected detections in each redshift bin. The number of sources detected at 2 mm is consistently lower than at 1.2 mm across all bins covered by the simulation (0 < z < 7). As expected from the study of the color versus redshift, they start to be similar around z∼ 7, where the number of detectable galaxies per redshift interval is a factor of 25 and 13 lower than at the peak around cosmic noon (z ∼ 2.5) at 1.2 mm and 2 mm, respectively. At cosmic noon, the number of galaxies detected at 2 mm is approximately half of that detected at 1.2 mm.

Fig. 6. Redshift distribution of the sources predicted by the 117 deg2 SIDES simulation for flux cuts at 1.2 mm (blue filled histogram) and 2 mm (red open histogram) similar to N2CLS in GOODS-N.

The higher mean redshift observed in the 2 mm-selected sample is thus primarily due to the reduced number of galaxies detected around cosmic noon, compared to the 1.2 mm selection, rather than an increased sensitivity to higher-redshift sources. Casey et al. (2021) argued that 2-mm surveys are a way to find the needles (z> 4 objects) in the haystack (cosmic-noon sources). However, since most of the cosmic noon sources can be identified using ancillary data (Berta et al. 2025), a 1.2 mm selection from which cosmic noon sources are removed is unlikely to produce more higher-z candidates to follow up than a simple 2 mm flux selection. If we use only 2 mm data, we lose information about the dust content of galaxies at cosmic noon (Berta et al. 2025) without adding a significant number of new z≳ 4 targets.

Cooper et al. (2022) promoted an alternative approach in which sources are selected at shorter wavelengths (850 μm) from a single-dish survey and then followed up with ALMA at 2 mm to reduce the uncertainties in submm photometric redshifts from σ_z/(1 + z) = 0.3 to σ_z/(1 + z) = 0.2. However, this approach may be less efficient in the JWST era. For most sources, having the location of mm sources with a precision of ≲1 arcsec is sufficient to find their JWST counterparts and obtain their photometric redshifts. The typical galaxy flux is a factor of ∼4 (∼10, respectively) times higher at 1.2 mm (850 μm, respectively) than at 2 mm, while ALMA single pointings are only a factor of 1.7 (2.2, respectively) times more sensitive at 2 mm. We can thus follow up ∼6 (∼20, respectively) times more sources with the same amount of telescope time at 1.2 mm (850 μm, respectively)¹¹. Therefore, it may thus be more appropriate to follow up the sources directly at shorter wavelengths (850 μm or 1.2 mm).

Finally, we used SIDES to understand how the selections at various wavelengths are biased toward or against different type of SEDs. The SIDES SEDs are parametrized using the mean radiation field ⟨U⟩. High values of this parameter correspond to high dust temperatures. In Fig. 7, we show the redshift evolution of the mean ⟨U⟩ parameter together with the 16–84 percentile range corresponding to 1σ in the Gaussian case. In SIDES, the ⟨U⟩ distribution varies with redshift but not with the stellar mass. We thus use the full simulated sample including low-mass galaxies¹² to provide our unbiased reference (in grey in the figure). We compare this sample with our flux selections at 1.2 mm (blue) and 2 mm (red) using the N2CLS GOODS-N flux cut. The mean values are biased toward lower ⟨U⟩ values (colder dust temperature) by 1σ at z< 1 and ∼0.5σ at z≳ 4. The 2 mm selection is slightly more biased toward high-redshift sources, by approximately 10%. This is not surprising, since colder dust SEDs emit more in the mm at fixed total dust luminosity.

Fig. 7. Mean interstellar radiation field ⟨U⟩ as a function of redshift, as predicted by SIDES, illustrating the bias toward colder SEDs introduced by mm flux selections. The mean radiation field is used in SIDES to parametrize the SEDs, and a higher ⟨U⟩ corresponds to a higher dust temperature. On the right y-axis, we also indicates the corresponding peak wavelength in Sν units of the main-sequence SED templates. The solid lines are the mean evolutions as a function of redshift for the full star-forming sample (gray), and selections at 1.2 mm (blue) and 2 mm (red), while the shaded areas and the dashed lines correspond to the 16–84% interval.

Our modeling analysis is based on the SIDES model, which does not contain a numerous population of dusty galaxies with very cold dust or at very high redshift. In Fig. 8, we show the relation between color and redshift for both SIDES and several modified black body SED models. For SIDES, the 1.2 mm is systematically more sensitive than the 2 mm up to z ∼ 7. We might consider whether some very cold sources missing in SIDES could be detected at 2 mm, but not at 1.2 mm. To explore this hypothesis, we computed the NIKA2 color expected from modified black bodies (ν^β B_ν(T)) with various dust temperatures T_dust and emissivity indices β. Galaxies with T_dust = 25 K and β = 1.5 at z≳ 3.5 are easier to detect at 2 mm. Galaxies with T_dust= 35 K and β= 1.5 or T_dust= 25 K and β = 2 are detected first at 2 mm for z≳ 5. For other SEDs with a higher T_dust or β, the 1.2 mm is more efficient to detect them at least up to z∼ 7. If dusty populations at z≳ 4 would be mainly objects with very cold dust or a low β, we would have found large populations of galaxies detected at 2 mm or with high S_2 mm/S_1.2 mm colors. Their absence in N2CLS is thus compatible with warmer SEDs, in line with other studies on the evolution the dust temperature with redshift based on small targeted ALMA samples (e.g., Faisst et al. 2020; Sommovigo et al. 2022) or stacking (e.g., Béthermin et al. 2015a; Schreiber et al. 2018; Viero et al. 2022), which found an increasing T_dust with redshift.

Fig. 8. Relation between the N2CLS color and the redshift for various SED models. The grey solid line and the associated contours represent the median color of SIDES galaxies and the 16–84% dispersion. The solid and dashed colored lines (red for 25 K, green for 35 K, and blue for 45 K) correspond to modified black bodies with β= 1.5 and β= 2, respectively. The horizontal lines show the sensitivity ratio in COSMOS (dot-dash line) and GOODS-N (dotted line).

4.3. Exploring the most efficient approach for mm high-z surveys

The 2 mm surveys as GISMO (Staguhn et al. 2014; Magnelli et al. 2019), N2CLS (Bing et al. 2023; Ponthieu et al. 2026), and (Ex-)MORA (Casey et al. 2021; Long et al. 2024) have opened the 2 mm window for deep and wide high-z galaxy surveys. Unfortunately, these surveys did not reveal any new exotic population. However, they were very important to clarify which dusty sources could be found at z ≳ 4. With all these data in hands, we can now conclude that the optimal wavelength range to probe the z ≳ 4 dusty galaxy population is around 1.2 mm, and future surveys should target it in priority. Coverage of the same sky area at 850 μm and 2 mm can help identify some high-redshift candidates; however, this identification remains uncertain due to large color measurement errors caused by limited signal-to-noise, map-making, and source extraction.

Having access to both interferometric and single-dish surveys provides an opportunity to better understand the impact of the angular resolution on the measurements (see Sect. 3.5) and figure out the best observational strategy for future programs. Since high-z dusty sources have a low surface density with ≪1 source per ALMA or NOEMA pointing, probing a large volume could be more efficient using a single-dish mapping and a follow-up of the detections with an interferometer to unambiguously identify their counterparts (e.g., JWST).

We used Ex-MORA and N2CLS as examples to evaluate the merit of a direct interferometric survey against the single-dish approach followed up by interferometric observations. The number of detections per hour is comparable for (Ex-)MORA (37 detections in ∼44 h, combining MORA and Ex-MORA) and N2CLS in COSMOS (90 detections in 195 h). However, the carbon footprint (see Appendix A) and economical cost per hour of observation are dramatically lower for N2CLS at the IRAM 30 m telescope (∼2 /h, private communication, ∼0.1 tCO₂e/h) than (Ex-)MORA at ALMA (∼60 /h, ∼5 tCO₂e/h). The follow up of all robust N2CLS sources in COSMOS would require only 90 ALMA or NOEMA pointings, while (Ex-)MORA is based on 4851 individual ALMA pointings. Thus, these follow-up observations are shorter by several orders of magnitude than a blind interferometric survey. The sum of the cost of a blind single-dish survey and its interferometric follow-up is much lower both on the ecological and economical point of view. This should not be viewed as a shortcoming of (Ex-)MORA, since its final results were published before N2CLS. Rather, this highlights the importance of carefully considering survey strategies when designing future programs. Our discussion applies only to the case of deep surveys. Since a wide range of scientific research can only be conducted with ALMA, efficient use of its observing time is crucial for the scientific community.

In the future, the performance of ALMA will improve (e.g., larger bandpass increasing the continuum sensitivity, potential future multi-beam array). Meanwhile, new generations of single-dish instruments will emerge such as the upgraded version of NIKA2 or TolTEC (Wilson et al. 2020), as well as future single-dish telescopes (e.g., AtLAST, Mroczkowski et al. 2025) will emerge. Therefore, it will be essential to regularly reassess the optimal survey strategy in light of the evolving instrumental landscape and the performance of both interferometric and single-dish facilities.

5. Conclusion

The N2CLS survey produced an unprecedented sample of 115 robust source detections selected at 2 mm (25 in GOODS-N and 90 in COSMOS). This enabled us to constrain their statistical properties, test our models, and gain deeper insights into the types of sources selected by such surveys. Our main results are as follows:

• We measured the mean S_2 mm/S_1.2 mm color of sources detected at both 2 mm and 1.2 mm in the COSMOS field and found 0.222 ± 0.008. Our E2E simulations based on the SIDES model and the N2CLS data themselves show that the full observational process produces biases of ∼10%. We also obtained a scatter of 0.070 ± 0.010 from the observed fluxes. Our simulation shows that the observational process can be responsible for most of this scatter and the intrinsic value could be significantly lower.

• We studied the redshift distribution of the 2 mm sources in both fields. In GOODS-N, we find a mean redshift of 3.6 ± 0.3, which is marginally higher than the prediction from SIDES (2.9 ± 0.2). This high value in GOODS-N is driven by the N2CLS overdensity at z ∼ 5.2 (Lagache et al. 2026). In the COSMOS field, we found a mean redshift of 3.2 ± 0.2, which is identical to that obtained in SIDES.

• We found that the observed S_2 mm/S_1.2 mm and S_2 mm/S₈₅₀ colors have a large scatter and barely evolves from z = 0 to z = 6. This suggests that this color cannot be used as a reliable way to select high-z dusty galaxies.

• In COSMOS, there are eight sources detected at 2 mm without counterpart at 1.2 mm. Two of these sources are radio galaxies at z < 1. The six other sources are compatible with the expected number of spurious detections. In GOODS-N, only one 2 mm source has no 1.2 mm counterpart, but it is firmly identified as a z ∼ 2 galaxy thanks to NOEMA follow-up observations.

• We compared the N2CLS 2 mm flux measurements with the Ex-MORA ALMA survey. We observed a NIKA2 flux excess relative to ALMA, not only for sources with multiple components but also for those with a single component, with an average excess factor of 1.3. This is explained by a combination of multiple factors such as the different bandpass of the two surveys, source blending, and peak flux photometry used by Ex-MORA.

• In this paper, we discuss the relevance of 2 mm surveys using both N2CLS and SIDES. We find that the slightly higher mean redshift of 2 mm surveys compared to 1.2 mm is mainly caused by a massive loss of cosmic noon sources, without any significant gain even at z ≳ 5. The 2 mm sample is also slightly more biased toward colder dust temperatures. Finally, we show that 2 mm single-dish surveys are more efficient than ALMA surveys.

Although the deep 2 mm N2CLS data could presumably have revealed exotic populations of dusty galaxies with very cold dust temperatures or at extremely high redshifts, we found no evidence of such sources. Ruling this out is important for optimizing the wavelength selection of future surveys targeting dusty star-forming galaxies. Our results show that 1.2 mm surveys are more efficient for building large samples from cosmic noon (z ∼ 2) to the reionization era (z ∼ 8).

Data availability

The data can be accessed at https://data.lam.fr/n2cls/

Acknowledgments

We thank Arianna Long for providing various Ex-MORA products allowing us to compare our results. We thank Jürgen Knödlseder for providing us the carbon footprint estimates of ALMA and the 30 m telescopes. 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 programme (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). We would like to thank the IRAM staff for their support during the NIKA and NIKA2 campaigns. 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. 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 benefited from the support of the European Research Council Advanced Grant ORISTARS under the European Union’s Seventh Framework Programme (Grant Agreement no. 291294). F.R. acknowledges financial supports provided by NASA through SAO Award Number SV2-82023 issued by the Chandra X-Ray Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of NASA under contract NAS8-03060. The NIKA2 data were processed using the Pointing and Imaging In Continuum software (PIIC, Zylka 2013; Berta & Zylka 2019), developed by Robert Zylka at the Institut de Radioastronomie Millimetrique (IRAM) and distributed by IRAM via the GILDAS pages. PIIC is the extension of the MOPSIC data reduction software to the case of NIKA2 data. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. 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). 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.

References

Appendix A: Estimate of the carbon footprint of ALMA and the IRAM 30 m telescope

Knödlseder et al. (2024) estimated the carbon footprint of worldwide astronomical research facilities. They provided us their best estimate of the yearly carbon footprint of ALMA and the IRAM 30 m telescope on which the NIKA2 camera is installed. For ALMA, they estimated 256 000 tCO₂e for the construction and 21 000 tCO₂e/yr for the operations. For the 30 m, the construction produced 4700 tCO₂e and the operations 300 tCO₂e. If we distribute the construction impact over an assumed 50 yr of operation, we obtain 26 000 tCO₂e/yr for ALMA and 400 tCO₂e/yr. for the 30 m. Thus, there is a ratio of 65 between the two carbon footprints. Both observatories declare ∼4000-5000 h of useful science data per year, although the 30 m tends to observe at lower frequency. Even when trying to correct for the better weather at ALMA for 2 mm observations, it is rather safe to assume that the hourly carbon footprint ratio between the two facilites is at least 1.5 orders of magnitude higher.

Appendix B: Tables

The N2CLS robust (95 % purity) 2 mm catalogs are presented in Table B.1 for GOODS-N and Table B.2 and B.3 for COSMOS. All the data used in this paper (1.2 and 2 mm fluxes from N2CLS, 850 μm fluxes from SCUBA2, redshifts, and simulations) can be accessed at https://data.lam.fr/n2cls/.

All Tables

All Figures

Fig. 1. Distribution of the 2–1.2 mm color for the true N2CLS COSMOS catalog (red filled histogram, only S/N ≥4.6 sources at both wavelengths and various simulated catalogs based on SIDES. The solid black histogram corresponds to the simulated SIDES galaxy catalog after applying a flux cut similar to COSMOS data (see Sect. 3.1). The yellow histogram is based on a similar selection, but applied to the blob catalog extracted from the noiseless simulated map. The grey histogram is the distribution obtained using the E2E simulation described in Sect. 2.5. The error bars on the N2CLS histogram are computed assuming a Poisson law, and those from the E2E are the standard deviation between the various simulated fields. The histograms are normalized to have a unity area (∫n(c) dc = 1, where c is the color). The lower limits on the color of 1.2 mm non-detections are shown as right-pointing triangles. Their y-axis position is arbitrary and was chosen to reduce overcrowding in the figure.

In the text

Fig. 2. Redshift distribution of the GOODS-N high-reliability 2 mm sources. The red histogram is the distribution measured in the N2CLS. When sources are multiple, the redshift of the main component is used (see Sect. 3.4). The uncertainties correspond to a Poisson law. The grey histogram is derived from the SIDES galaxy catalog after applying a flux cut similar to the GOODS-N field (see Sect. 3.2). The uncertainties correspond to the field-to-field variance of a 15.9 × 10 arcmin2 field. The yellow histogram is based on a similar selection, but applied to the SIDES blob catalog and using the redshift of the brightest galaxy in the beam. The histograms are normalized to have a unity area (∫n(z) dz = 1).

In the text

	Fig. 3. Same as Fig. 2, but for the COSMOS field. The blue dashed histogram represents the Ex-MORA survey (Long et al. 2024).
In the text

Fig. 4. Color as a function of the redshift in the COSMOS (top) and GOODS-N (bottom) fields. The left plots show the internal NIKA2 color (S2 mm/S1.2 mm) and the right ones show the ratio between the N2CLS 2 mm and the SCUBA2 850 μm fluxes from S2COSMOS (Simpson et al. 2017) and SUPER GOODS (Cowie et al. 2017). The filled blue circle are the N2CLS sources with spectroscopic redshifts, while the open red circles corresponds to photometric redshifts. The grey dots are from the SIDES 2 deg2 simulated catalog after applying a 2 mm flux selection similar to SIDES (see Sect. 3.2). The yellow solid line corresponds to the median of the SIDES blob catalog, using the redshift of the brightest galaxy in the beam, and the colored area correspond to the 1σ confidence region. The dashed line is the mean color as a function of redshift from SIDES 2 deg2 without applying any selection. We also show the sensitivity ratio (σ2 mm/σ1.2 mm where σ2 mm and σ1.2 mm are the RMS of the noise at 2 mm and 1.2 mm, respectively) for ALMA (dotted line) and N2CLS (dot-dashed line). See the discussion in Sect. 4.1 for more details. The right-side panel of each plot show the color distribution of sources with a photometric redshift (red), a spectroscopic redshift (blue), and no redshift (green).

In the text

	Fig. 5. Comparison between the exMORA (x axis, Long et al. 2024) and N2CLS (y axis) 2 mm flux. The N2CLS sources with a single counterpart are in blue and the multiple sources are in red. The black dashed line indicates the one-to-one relation and the dotted line includes the various corrections (bandpass, source blending, and source sizes) discussed in Sect. 3.5.
In the text

	Fig. 6. Redshift distribution of the sources predicted by the 117 deg2 SIDES simulation for flux cuts at 1.2 mm (blue filled histogram) and 2 mm (red open histogram) similar to N2CLS in GOODS-N.
In the text

Fig. 7. Mean interstellar radiation field ⟨U⟩ as a function of redshift, as predicted by SIDES, illustrating the bias toward colder SEDs introduced by mm flux selections. The mean radiation field is used in SIDES to parametrize the SEDs, and a higher ⟨U⟩ corresponds to a higher dust temperature. On the right y-axis, we also indicates the corresponding peak wavelength in Sν units of the main-sequence SED templates. The solid lines are the mean evolutions as a function of redshift for the full star-forming sample (gray), and selections at 1.2 mm (blue) and 2 mm (red), while the shaded areas and the dashed lines correspond to the 16–84% interval.

In the text

Fig. 8. Relation between the N2CLS color and the redshift for various SED models. The grey solid line and the associated contours represent the median color of SIDES galaxies and the 16–84% dispersion. The solid and dashed colored lines (red for 25 K, green for 35 K, and blue for 45 K) correspond to modified black bodies with β= 1.5 and β= 2, respectively. The horizontal lines show the sensitivity ratio in COSMOS (dot-dash line) and GOODS-N (dotted line).

In the text