© The Authors 2026

1. Introduction

In the far-infrared (FIR) to radio domains, single-dish surveys face a fundamental limitation in sensitivity and photometric accuracy due to fluctuations in the background surface brightness. In surveys of so called cosmological fields (i.e. free from Galactic contamination, the dominant source of background fluctuations is unresolved extragalactic point sources, which contribute to confusion noise in both the radio (Condon 1974) and IR domains (Kiss et al. 2001; Lagache et al. 2003; Dole et al. 2004). This confusion noise establishes the fundamental sensitivity limit for single-dish surveys, particularly for blind detection and accurate flux measurements of point sources.

The angular resolution of a single-dish telescope directly depends on the observed wavelength, limiting its ability to resolve individual sources at longer wavelengths. For example, the Spitzer Space Telescope’s 85 cm mirror (Werner et al. 2004) resolved over 80% of the galaxies contributing to the total emission at 24 μm (Papovich et al. 2004), but only 23% at 70 μm and 7% at 160 μm (Dole et al. 2004). Similarly, the Herschel Space Observatory (Pilbratt et al. 2010), with its larger 3.5 m mirror, was able to resolve 58% and 74% of the sources at 100 μm and 160 μm, respectively, but only 25% at 250 μm (Berta et al. 2011; Oliver et al. 2010). To overcome these resolution limitations, innovative techniques have been developed, such as stacking analyses of known galaxy populations (Dole et al. 2006; Béthermin et al. 2012; Viero et al. 2013a) and analyses based on the probability of deflection to study unresolved source fluctuations (Glenn et al. 2010; Berta et al. 2011; Béthermin et al. 2012; Viero et al. 2013b). Additionally, studying the anisotropies of the cosmic infrared background (CIB) offers valuable insights into large-scale distribution and clustering of dusty star-forming galaxies, as well as the properties of the CIB itself (Lagache et al. 2000; Amblard et al. 2011; Planck Collaboration XVIII 2011; Viero et al. 2013b; Planck Collaboration XXX 2014). This helps achieve an effective approach to constrain the evolution of galaxy populations and cosmic star formation processes.

Following the initial theoretical studies of confusion noise in the radio domain (Condon 1974), its first empirical assessment was made in a 4 GHz radio survey using a half-difference technique to separate instrumental noise from observed noise (Ledden et al. 1980). This methodology has been extended to the FIR and submillimetre (submm) range, where confusion levels attributed to extragalactic sources and galactic cirrus fluctuations have been predicted (Helou & Beichman 1990; Hacking & Soifer 1991; Gautier et al. 1992). In the FIR domain, confusion noise was measured with ISOPHOT data (Kiss et al. 2001), which revealed fluctuations from extragalactic sources. Confusion noise, which represents the ultimate sensitivity limit for source detection, has become a critical parameter for assessing the performance of extragalactic deep fields. For the Spitzer Space Telescope, studies of confusion noise were undertaken to optimise deep IR surveys, particularly in regions with high source densities (Kiss et al. 2005). Observations from the MIPS instrument at 24 μm quantified confusion limits imposed by the high density of faint unresolved sources (Dole et al. 2003), influencing the design of deep field surveys (Dole et al. 2004). Measurements at longer wavelengths, such as 70 μm and 160 μm, revealed the combined role of extragalactic fluctuations and Galactic cirrus in setting confusion limits, providing a means for reliable source detection amidst the noise background (Frayer et al. 2006; Papovich et al. 2004).

In the submm range, confusion noise also determines the depth of the survey. This was first addressed by Blain et al. (1998), who highlighted the impact of confusion noise from faint, dusty galaxies on single-dish observations. The problem was further explored with SCUBA at the JCMT, where confusion from unresolved sources set the practical sensitivity limits for wide-field surveys (Hughes et al. 1998; Barger et al. 1998; Eales et al. 1999). More recently, observations with SCUBA-2 refined the understanding of submm confusion, particularly at 850 μm. Geach et al. (2017) measured the confusion noise at 0.42 mJy, highlighting the critical role of clustering in shaping the background fluctuations and limiting source extraction. The Herschel Multi-tiered Extragalactic Survey (HerMES) provided a detailed assessment of confusion noise using SPIRE, quantifying detection limits at 250, 350, and 500 μm (Nguyen et al. 2010). Follow-up analyses of Herschel data deepened our understanding of confusion noise contributions from clustered dusty galaxies and Galactic cirrus, aiding in the optimisation of survey strategies and robust source extraction (Lagache et al. 2003; Negrello et al. 2004; Nguyen et al. 2010). These studies have proven instrumental in improving our ability to characterise the faint submm galaxy population, which dominates the cosmic IR background.

Different methods are commonly used to estimate the level of confusion noise, especially in the radio, FIR, and submm bands. One category of methods is based on galaxy number count models. A widely used approach calculates the number of sources per beam based on the source count distribution, dN/dS (Dole et al. 2003). These methods estimate the confusion noise as the flux density variance induced by unresolved sources within the telescope beam. The flux cut-off threshold can be defined by a photometric criterion, such as a signal-to-noise ratio (S/N) threshold, or a source density criterion that assumes a uniform spatial distribution of sources (Takeuchi et al. 2001). While effective in sparse source environments, this approach becomes less reliable in clustered fields where source correlations amplify background fluctuations, rendering the uniformity assumption invalid.

Another type of method is data-driven and derives the confusion noise directly from the map noise measurements. In this approach, the confusion noise (σ_c) is derived as the asymptotic limit of the map noise (σ_m) as a function of observation time. The relation σ_m² = σ_i² + σ_c² is used, where σ_i represents the instrumental noise, which decreases with time as t^−0.5 (Nguyen et al. 2010; Geach et al. 2017). This technique has the advantage that it does not rely on assumptions about the underlying source distribution. However, it is sensitive to systematic effects in the instrument or data reduction processes that may not diminish with longer observing times, potentially biasing the estimated confusion limit.

These methods highlight a trade-off between model-based and empirical approaches. While number count models make use of well-understood theoretical frameworks, they falter in complex fields with clustered sources. Conversely, empirical methods offer robustness to such clustering but are vulnerable to systematic errors, highlighting the need for careful cross-validation between methods in confusion-limited surveys.

NIKA2 is the current IRAM camera for measuring the millimetre continuum emission of astrophysical objects at the 30 m telescope (Pico Veleta, Spain). It is described in detail, together with its performance in Perotto et al. (2020). In short, it is a dual-band camera whose detectors are kinetic inductance detectors (KID) (Monfardini et al. 2010). The bands, which are observed simultaneously, are centred at 260 and 150 GHz (1.2 and 2 mm, respectively). There are three focal planes (arrays of detectors), two operating at 260 GHz and consisting of 1140 KIDs each, as well as one operating at 150 GHz with 616 KIDs. About 90 and 84% of these detectors are effectively used for observations. The photometric bandwidths are 49 and 40 GHz. The main beams have FWHMs of 11.1 and 17.6 arcsec (at 1.2 and 2 mm), and the field of view (FOV) is of 6.5 arcmin in diameter. This large field of view combined with point source sensitivities of 30 and 9 mJy.s^1/2 provides mapping speeds of 111 and 1388 arcmin² mJy⁻² h⁻¹ at 1.2 and 2 mm respectively, an order of magnitude better than the previous generation photometric instruments at the 30 m (Perotto et al. 2020). Although it is not used in this work, linear polarisation measurements can be made with the 1.2 mm channel. NIKA2 has been available to the community since 2017. In return for its efforts to design, fund, install, and characterise the instrument, the NIKA2 collaboration¹ has been granted 1300 h of guaranteed time of observation that have been divided into five large programmes. This paper is one of those presenting the results of the NIKA2 Cosmological Legacy Survey (N2CLS), which consists of a deep integration of two extragalactic fields with many available ancillary data: GOODS-North and COSMOS.

This paper is organised as follows. In Sect. 2, we present the GOODS-North field and our observations. In Sect. 3, we describe the data processing that has been devised to produce the scientific maps. Sect. 4 presents the confusion noise estimator that we have built, the obtained measurements and their validation. In Sect. 5, we discuss the physical interpretation of the measured confusion. We present our summary and conclusions in Sect. 6.

2. Observing the GOODS-North field

GOODS-North is one of the two fields of the NIKA2 Large programme “NIKA2 Cosmological Legacy Survey” (N2CLS). An introduction to this IRAM 30 m Large Programme and a first analysis of its observations in terms of deep cosmological number counts have been proposed in Bing et al. (2023). For the sake of clarity, we recall here the main features of our observations that are relevant to this work.

From October 2017 to March 2021, across 12 NIKA2 runs (i.e. one or two consecutive weeks of observations), we have observed GOODS-North 749 times. In the following, we will refer to each of these observations as a ‘scan’. The observed region is a rectangle of 12′×6.3′ centred on RA = 12:36:55.03 and Dec = 62:14:37.59. Each of these raster scans is a collection of ‘subscans’ that encompass a straight trajectory of the telescope, that are parallel to each other and that are separated by 20″. A scan has its subscans oriented at 40 or −50° in (RA, Dec) so that we cover exactly the same region each time, but with alternating orthogonal directions. This was decided in anticipation of robustness tests (see Sect. 4.3). The scanning strategy on GOODS-North was designed to reach sensitivities comparable to the predicted confusion at 1.2 and 2 mm, respectively, estimated from simulations to be about 140 and 40 μJy/beam (Béthermin et al. 2017), in less than 100 h of observations. Each scan lasted about 7 and 8 mn.

During these 4.5 years of observations, some scans have suffered from poorer atmospheric conditions or instrumental instabilities. Thus, we rejected the ones that proved to be more of a liability than an asset. Our main criterion was the final average noise of an individual scan map compared to that of all the other scans. For this work, to reduce the noise in the data, we used an improved version of the pipeline described in Perotto et al. (2020), which has been used in previous papers published by our collaboration² (see Sect. 3.1). Although it is similar in its basic principles as far as noise reduction is concerned, this method also differs from the one used in Bing et al. (2023). Indeed, the need for extensive end-to-end simulations dedicated to this work was more conveniently achieved with the NIKA2 collaboration pipeline. Our methods have different exclusion criteria when it comes to selecting scans. In this work, we ended up using 668 scans, both at 1.2 and 2 mm, thereby giving 78.2 h of observation on the field. As for the opacities, the mean values are 0.25 and 0.15, with minima of 0.1 and 0.07, and maxima of 0.5 and 0.3 at 1.2 and 2 mm, respectively.

3. Data reduction

In this section, we present our data reduction pipeline. We retained the Planck idiom and refer indifferently to timelines or time-ordered information (TOI) for the data streams recorded by the detectors over time.

3.1. From raw to calibrated data

Most of the low-level processing that goes from the raw data to the calibrated data is the same as that described in detail in Perotto et al. (2020). Here, we present only the modifications that are relevant to this work. In particular, we have used a novel method to monitor the KID resonance frequency shifts and convert them into a TOI that is proportional to the incident flux. This method, described in detail in Appendix C, improves on the previous one in terms of the numerical integration, which improves the residuals from the atmosphere and low frequency electronic noise. Our data model has also changed. In Perotto et al. (2020), for each detector, we estimated a single composite mode of atmosphere and electronic noise from the other detectors that are most correlated with it. This composite mode was computed while taking into account a bright source mask and both can be improved by an iterative scheme. While this method has proven effective for all the other science objectives of NIKA2, the level of integration is so much deeper in this work that it proved insufficient and led us to explicitly separate the atmosphere and several different ‘modes’ of the sky noise and the electronic noise. We explored several approaches based on more sophisticated instrument models. The main idea was to use the distribution of detectors per electronic box and their analogue sub-bandpasses to derive educated models of the noise modes. Although promising, this idea did not meet our expectations. The output maps were indeed very clean, but the transfer function (see Sect. 4.2) was smaller than 50%. This led to too much uncertainty in the actual blind detection of the bright sources to be masked and to a large and more uncertain correction in the final quoted values of the confusion.

In the end, we limited ourselves to one mode per acquisition electronics box. This mode, however, is not built up from the KID TOIs, but from measures of the same feedlines where there are no KID resonances. This is the analog of dark or blind detectors in other instruments and we refer to them as ‘tones’. There are approximately 140 such measures available for each of arrays 1 and 3 (both at 1.2 mm) and 18 for array 2 (at 2 mm). This relatively low number is a limitation for our 2 mm channel. The same circle method is used for these off-resonance tones as for the standard KIDs to derive a frequency variation. Empirically, we have found that the signal that is orthogonal to the circle frequency variation (see Appendix C) is not useful for the decorrelation, so we used only the tangential on-circle variations, in the same way as for KIDs. Therefore, we here modelled the TOI of detector k as

$$\begin{array}{cc}\hfill m_t^k& ={\gamma }^k{P}_{\mathit{tp}}^k{S}_p+{\alpha }_0^k{A}_t^{\mathrm{LF}}+{\alpha }_1^k{A}_t^{\mathrm{HF}}+{\alpha }_2^k\frac{\partial A_t}{\partial t}+{\alpha }_3^k{A}_t^2\hfill \\ \hfill & +{\alpha }_4^{k}e{l}_t+{\alpha }_5^k+{\alpha }_6^{k}t+{\displaystyle \sum _{b=1}^{N_b}}{\beta }_b^k{E}_t^b\hfill \\ \hfill & +{\displaystyle \sum _{i=1}^{N_T}}({\eta }_i^{k}cos(2\pi {\nu }_{s}it)+{\xi }_i^{k}sin(2\pi {\nu }_{s}it))+n_t^k.\hfill \end{array}$$(1)

In Eq. (1), t is time, P_tp is the pointing matrix, and S_p is the sky brightness map in pixel, p. Here, A_t is the atmosphere as seen by all the KIDs, decomposed into its low frequency part: A_t^LF below 0.5 Hz and its high frequency part: A_t^HF = A_t − A_t^LF. To the first order, all KIDs see the same atmospheric signal, since its origin is in the near field of the telescope. Then, el_t is the elevation of the telescope and E_t^b is the series of noise components derived from off-resonance tones, which are partially correlated with a range of different detectors depending on their common readout electronic box, b. Then, n_t^k is the individual KID Gaussian white noise. We also included the time derivative of the atmosphere, $\frac{\partial A_t}{\partial t}$, as a proxy for second-order variations of the atmospheric contribution across the FOV, for possible residual non-linearity with the A_t² term and a time linear drift, α₅^k + α₆^k t. Finally, on a subscan basis, we also defined harmonic modes to filter out low-frequency residual noise. To filter short and long subscans equally, we set the effective number of trigonometric modes per subscan, so that the highest subtracted frequency was 0.2 Hz. In practice, this added N_T = 4 modes (2 cosine and 2 sine functions) for each of the 10 s subscans (6 × 12 arcmin case) or 6 modes for each of the 17 s subscans (12 × 6 arcmin case). Here, γ^k is the absolute calibration, known from specific observations (see Perotto et al. 2020) on Uranus and skydips, while the α^k, β^k, η^k, and ξ^k coefficients are determined during the processing, as described in the next subsection. The actual implementation of this model goes along with the map making presented in the next section.

3.2. Noise reduction and map making

In all of this work, we built maps at a resolution of 4 arcsec, so that it is less than a third of the main beam FWHM of the instrument, as the standard proxy for Nyquist’s criterion. Our data reduction scheme is then:

1. For each time, t, the median of the simultaneous samples of all the KIDs in the same array is computed. The median ensures robustness against possible outliers due to possible bright sources or anomalies that are seen by only a few KIDs at the same time, compared to the ∼1000 KIDs of the array. We internally iterate once on the mode production by eliminating KIDs that do not correlate with this initial mode (e.g. due to a bad tuning) and take the average of all selected KID TOIs.

2. This mean mode is dominated by the atmosphere and can therefore be used as an estimate of A_t. Its derivative ∂A/∂t is derived after a 0.5 s smoothing.

3. One noise mode per acquisition box E_t^b is built with the average signal of the off-resonance tones of this box.

4. Each KID timeline, m_t^k, is simultaneously linearly regressed against the atmosphere modes, the noise modes, and the trigonometric functions to estimate all the coupling coefficients α^k, β_i^k, η_j^k, and ξ_j^k. This regression is performed on a subscan basis to better account for the variability of the KID coupling to the different electronic modes. We also note that all electronic box modes are used altogether in the fit for each KID, and not just the box mode of that KID. This is to cope with the mixing of the electronic modes when A_t and ∂A/∂t are generated using all KIDs of the same array.

5. A few percent of the KIDs show a poor fit to this model and are excluded at this stage.

6. For each remaining KID, the estimated atmosphere and noise contributions are subtracted, leaving only γ^kP_tp^kS_p + n_t^k in Eq. (1).

7. Once all the KID timelines from a scan have been processed in this way, they can be projected onto a map. To do this, we use a standard nearest grid point average and apply a weight to each data sample. This weight is the inverse of the variance of the KID timeline during the subscan to which the data sample belongs. For the record, we note a slight overall excess of variance during the first two seconds of each subscan and we down-weight the corresponding data samples accordingly.

8. Low-level residuals can be seen in maps in Nasmyth coordinates, with a striping oriented in the same way as the electronic boxes. We decided to correct for this pattern with one template per array and per subscan.

Each scan was reduced according to this scheme and we obtained a map per scan. We then co-added all these scan maps to produce the final map, using inverse noise per pixel weighting. This process is based on the implicit assumption that the signal is zero or at least negligible compared to the parasitic components. To improve the final map, we iterated this process. We performed a point source detection on the map and kept only the brightest ones that appear at a SNR greater than 10. We built a model map that is just the sum of these bright point sources modelled as constant-width 2D Gaussian. At the start of the next iteration, we de-projected this model map into point source TOIs, which we then subtracted from the KID TOIs. In this way, the final map was no longer affected by the filtering residuals of the bright sources. In practice, we proceeded more progressively. At each iteration, iter, we actually subtracted (1 − 0.2^iter) times the source model map built at the previous iteration to mitigate the effect of errors on this map at the very first iterations. This process was iterated until convergence was achieved, with residual peaks on the difference between successive iterations at less than 10 μJy. In practice, three processing iterations were sufficient, the first one assuming a vanishing point-source model map. Figure 1 presents our final maps, together with the contours of the region used to derive the confusion (cf. Sect. 4). The noise map was derived in two steps. As described in step 7 of our data processing, each sample was projected onto the final map with a known weight. Straightforward algebra allowed us to derive the total noise per pixel on the final map. However, we note that the distribution of the S/N per map pixel did not follow the expected normalised Gaussian, probably due to some residual correlation in the noise. We concluded that our TOI noise propagation underestimated the effective noise on the map. To correct for this, we needed to multiply our noise maps by typically 1.6 and 1.3 at 1.2 and 2 mm, respectively.

Fig. 1. Maps of GOODS-North used in this work (left 1.2 mm, right 2 mm). The black contours define the region (i.e. mask) that we use to compute the confusion.

4. Measuring the confusion

The confusion, σ_c, is defined as the square root of the variance of the measured brightness, S (i.e. including the point spread function; PSF or beam of the instrument) when bright sources above a certain flux are masked out and in the limit of zero instrumental noise or (equivalently) infinite integration time (Condon 1974; Nguyen et al. 2010; Béthermin et al. 2024) as follows:

$$\begin{array}{c}\hfill {\sigma }_c^2=\frac{1}{N}{\displaystyle \sum _{p=1}^N}{(S_p-\overline{S})}^2,\end{array}$$(2)

where S_p is the signal at pixel, p, and $\overline{S}$ is the average brightness of the map. With this definition, the confusion depends on both the angular resolution of the instrument and the arbitrary flux cut of bright sources. In practice, maps are affected by the data processing and contaminated by instrumental noise so that Eq. (2) cannot be readily implemented. In the next subsection, we present our estimator of σ_c. Our final results and a discussion on their robustness are summarised in Sect. 4.3.

4.1. The cross-variance estimator

To define the sky area over which we calculate the confusion, we have to make a trade-off between the size of the patch and the homogeneity of the noise. To this end, we chose to restrict to pixels with a number of samples greater than two times the median number of hits per pixel on the full map for both bands. We built on the results of Bing et al. (2023) and masked sources that are confidently detected above a minimum S/N of 4.6 and correspond to fluxes greater than ∼0.54 and ∼0.17 mJy at 1.2 and 2 mm, respectively. We took the locations of these sources and mask all pixels within a radius of 19 arcsec at 1.2 mm and 27.8 arcsec at 2 mm. This leaves us with 64.7 and 62.1 arcmin² at 1.2 mm and 2 mm (see Fig. 1). The masks were independently derived at 1.2 and 2 mm.

Implementing Eq. (2) directly on the data leads to an estimate of ${\widehat{{\sigma }_c}}^2$ that is biased. In fact, the measured value in each pixel, A_p, is the sum of the signal, S_p, and the noise a_p, so the two variances add up. Furthermore, the signal has been affected by the data processing. In our case, where we are interested in a single number, the effect of the processing (i.e. the transfer function) can be represented by a single factor, f. In the limit of zero noise on the maps, replacing S_p by A_p in Eq. (2) gives f²σ_c². The derivation of f is detailed in Sect. 4.2 and is based on simulations. We went on to determine how we could estimate the variance of the noise, σ_a². We considered several possibilities.

First, we could use a null map produced by combining the 668 scans with alternating positive and negative weights. This would cancel out the astronomical signal, while leaving the noise properties unchanged. Applying Eq. (2) to this map would directly give σ_a². In practice, the noise on the null map actually appears to be less than the noise on the map. Indeed, when we computed the histogram of the S/N, with the noise derived from the null map, its width is rather 1.2 than 1. We attribute this to residual noise in the data maps that is difficult to estimate precisely and that may bias our estimate of σ_c. Nevertheless, if we computed ${\sigma }_c=\sqrt{{\sigma }^2-{\sigma }_a^2}$ using this method, we would find σ_c ≃ 145 and 34 μJy/beam at 1.2 and 2 mm, respectively. These values are in good agreement with our final results.

Second, we could use an approach such as the one proposed in Nguyen et al. (2010) and monitor the measured noise as a function of the integration time. Unlike them, we did not have enough integration time to do this analysis pixel per pixel, but we could generalise their approach to a sub-area of our map, accumulating scans, calculating the variance of the map pixels as we went along, and then fit its decrease as a linear function of the inverse of the integration time. The noise variance would decrease as 1/t and the confusion would be the remaining constant term in the fit³. We found that this does not solve the previous problem and comes with other difficulties, such as the determination of the (strongly correlated) uncertainties of each measurement of σ²(t) and the determination of the out-of-atmosphere equivalent integration time. Both of these factors impact the derived value of the confusion and its associated uncertainty. The former can be addressed by Monte-Carlo simulations, but the latter is poorly defined in the context of a non-uniform survey and is affected by uncertainties in the opacity as well as the overall absolute calibration; this can vary between different runs, depending on the number and quality of absolute calibration scans that could be performed. In Fig. 2, we still show the decrease of the variance of the map versus the integration time, but this is for illustrative purpose only. As such, the solid black line follows t^−1/2, which was fitted for integration times smaller than 20 000 s and it does not take the correlation between the error bars into account. For reference, completing this exercise yields σ_c≃ 142 and 49 μJy/beam at 1.2 and 2 mm, respectively. Given the actual significance of our detection (see below), the observed excess is real and dominated by the confusion and the agreement between this method and our final one is not surprising. However, the caveats we mention here have prevented us from using this estimator to derive our final results.

Fig. 2. Variance of the N2CLS maps vs. the out of atmosphere integration time on a region in the center equivalent to one NIKA2 field of view (6.5 arcmin diameter). To determine these data points, we simply take the variance of the unmasked pixels (see Fig. 1) in the central area. For each NIKA2 band, the lower plot shows the residuals compared to the t−1/2 instrumental noise integration that is fitted here for times of integration smaller than 20 000 s. The error bars are the square roots of the diagonal elements of the covariance matrix of the measures (Fig. A.2), derived from 100 Monte-Carlo simulations of map-based realistic noise (per scan). There is a clear excess of signal at high integration time due to the confusion, but as explained in the main text, this plot is for illustrative purposes only and is not used to derive our final results.

To improve on all these aspects, we developed a third approach based on the general idea of cross-correlations (e.g. Tristram et al. 2005). Instead of combining all the scans into a single map, we combined them into a set of independent maps (1 for each of the 12 observation runs). This gives us 66 pairs of different maps of the same sky (excluding auto-correlations) and for each pair (A_i, B_i), we can compute the cross-variance, as defined by:

$$\begin{array}{c}\hfill {\rho }_i\equiv \frac{1}{f^2}\frac{1}{N}{\displaystyle \sum _p}(A_p^i-{\overline{A}}^i)(B_p^i-{\overline{B}}^i),\end{array}$$(3)

where i is the pair index. Since the maps Aⁱ and Bⁱ are independent, the statistical expectation of ρ_i for unmasked pixels is directly σ_c² (see Eqs. (A.4) and (A.11)), without any bias. Indeed, since maps Aⁱ and Bⁱ were not coadded before computing the variance of the sum, their noise components did not add up in quadrature as in the previous two estimators. Finally, ρ_i makes no assumption about the uniformity of the integration time per pixel and is therefore immune to inhomogeneities in the sky coverage. To make the best use of all the available data, we computed the cross-variance of all the possible pairs and determine their full covariance matrix, X (see Appendix A for details). With these two in hand, we were able to compute the likelihoods of the confusion at 1.2 mm, 2 mm, and even the cross-confusion 1.2 × 2 mm:

$$\begin{array}{c}\hfill \mathcal{L}\propto exp[-({\rho }_i-{\sigma }_c^2)X_{\mathit{ij}}^{-1}({\rho }_j-{\sigma }_c^2)].\end{array}$$(4)

We relied on these likelihoods to derive final values and their associated confidence intervals. Our final results are discussed in Sect 4.3.

4.2. Determining the transfer function

To determine the transfer function factor, f, we relied on the state-of-the-art SIDES-UCHUU simulations (Béthermin et al. 2017, 2022; Gkogkou et al. 2023). We applied the NIKA2 bandpass and point spread function to the SIDES-UCHUU cube to create a simulated sky, sim_p, as it would be observed by NIKA2. We then ran our exact scan strategy and KID selection over this map, for each scan, to build simulated signal timelines, sim_t^k, for each KID, k. These simulated timelines can then be added to our real data timelines. To be immune to the true astronomical signal in our data, before adding the simulated timeline, we multiplied our data timelines by −1 for every second pair of scans (to account for both scan orientations), so that the true signal is cancelled out in the final map, leaving only the simulated component. This is the best possible simulation to characterise our pipeline, as it deals with the true atmosphere and electronic noise rather than approximations. The magnitude of the simulated timelines is smaller than these two components by a factor of more than 1000, so it does not bias their derivation. We then reduced these composite timelines in exactly the same way as the data and project them onto a map, M_p = sim_p^processed + N_p, where N_p is the noise. Basically, to determine f, we need to measure the rms of sim_p^processed and divide it by the variance of the input noiseless simulated map, sim_p. A first way to proceed is to process the data timelines with exactly the same mask and alternating sign weights, but without adding sim_t^k. This gives N_p which can be subtracted from M_p to get sim_p^processed. A second way is to do the same whole process as on the data, namely making maps per runs, calculating their cross-variance and deriving the estimated final confusion as the value that maximises their likelihood. The ratio of this value to the standard deviation of the input map, sim_p, on the same area gives an estimate of f. As described in Sect. 2, our observations consist of two sets of scans with orthogonal directions and different subscan lengths. This results in slightly different effective filtering between these two types of scans, so we computed a transfer function, f, per direction and we corrected accordingly. We discuss this in more detail in Sect. 4.3. We checked the robustness of our estimates of f by varying the mask definition criteria (size and masked sources) and by scaling the SIDES simulations by factors ranging from 1 to 3. The scaling improves or degrades the S/N of the output signal, while remaining negligible compared to the true noise of the data. Finally, we used the f factors derived from simulations scaled by a factor 3, as they provide a better dispersion due to their higher S/N. We estimated f on 27⁴ simulations and the two methods agree on average by 1%. The statistical uncertainty on the average of f is also about 1%, so we take this as our final uncertainty on f.

The dispersion of σ_c as estimated from these 27 end-to-end simulations (including the detection of bright sources and masks) also provides an estimate of the cosmic variance associated with our measurement. We find about 12, 4, and 6 μJy/beam at 1.2, 2, and 1.2 × 2 mm, respectively (exact values are given in Table 1).

4.3. Results

The 12 observation runs give us 66 cross-variance measurements, which are given in Fig. A.1 in the particular case of 1.2 mm. The likelihoods of the confusion at both wavelengths are shown in Fig. 3 as estimated on different subsets of our data. ‘Half 1’ denotes the ensemble of 334 of our scans that are oriented at 40°, ‘Half 2’ is the 334 at −50°. To check for consistency, we compute the confusion either exclusively on each half of the scans, on ‘all’ the scans combined or exclusively by crossing one half with the other. As mentioned in the previous section, we computed an effective transfer function, f, and corrected for it independently in each case. Indeed, this f factor relates to the variance of the signal and is implicitly azimuthally symmetric, whereas our filtering is mostly per subscan and is therefore directional. This particular feature prevents us from deriving f based on the combined directions taken from f coming from each direction. As shown in Fig. 3, all estimates are compatible within 1σ. Using the values derived from all the scans (as they have the highest S/N), we obtained the confusion estimates shown in Table 1; namely, $139.1_{-19.2}^{+15.9}\pm 11.9$, $38.6_{-13.1}^{+9.6}\pm 3.7$ and $49.6_{-24.8}^{+15.9}\pm 6.4$ μJy/beam at 1.2, 2 and 1.2 × 2 mm, respectively. The first set of error bars refers to the instrumental noise and the second to the cosmic variance. These are clear measures at a significance of 7σ at 1.2 mm and 3σ at 2 mm, but only a 2σ evidence at 1.2 × 2 mm. An overall absolute calibration uncertainty has to be considered. It affects the quoted value and its statistical error bars, but not the cosmic variance. This uncertainty has three main components. Firstly, there is the statistical variability of the measurement under various observation conditions from the ground, which was determined to be 5% during the NIKA2 commissioning. Secondly, there is uncertainty regarding the exact flux of our calibration source Uranus. As reported in Perotto et al. (2020), this uncertainty is about 5%. Taking all these uncertainties into account and adopting a conservative approach, we estimated an overall uncertainty of 10% in our absolute calibration.

Fig. 3. Likelihoods (normalised to their maximum) of the confusion at 1 mm, 2 mm, and the cross-confusion 1 × 2 mm for different subsets of our data. The central points and their 1-σ equivalent error bars are shown for convenience and have arbitrary ordinates. ‘Half 1’ refers to scans oriented at 40°. ‘Half 2’ refers to scans oriented at −50°. ‘All’ refers to the total sum of the two.

To check the consistency of our final values, we used the SIDES simulations and the bright source masks derived from our end-to-end simulations and calculated the expected confusion levels. We find 156, 49, and 84 μJy/beam at 1.2, 2, and 1.2 × 2 mm. If we compare this to our measures taken at face value, at 1.2 and 2 mm, our results are compatible with SIDES at 1σ (stat), even without considering the uncertainty onabsolute calibration or the cosmic variance. Our measurement of the 1.2 × 2 mm confusion is 1.4σ (stat) lower than the value expected from SIDES. Given the 10% uncertainty on our absolute calibration, if we allow it to increase by this factor, our results at 1.2 and 2 mm become 153 ± 21 and 43 ± 14. This shows an even better agreement with the prediction of SIDES. Our measurement at 1.2 × 2 mm becomes 55 ± 28 μJy/beam, which is now 1.04σ (stat) of the 84 μJy/beam expected from SIDES. This last compatibility can be further improved if we assume that our field has a low intrinsic confusion within one cosmic standard deviation (Table 1). This would actually also improve the agreement at 2 mm, but slightly degrade that at 1.2 mm, although keeping it within 1σ (stat). Overall, this comparison shows a good agreement between our results and the predictions of SIDES at 1.2 and 2 mm, but a marginally lower than expected value at1.2 × 2 mm.

The consistency of our results can be assessed in another way. Using Eq. (6), we can see that σ_c2 = rσ_c1θ₂/θ_1.2 and σ_c1 × 2 = (2r)^1/2σ_c1θ₂/(θ₁² + θ₂²)^1/2, where θ is the instrumental beam FWHM and r is a constant that can be derived from the SIDES confusions and is 0.21. Applying this to the measured value of the confusion at 1.2 mm and considering only the statistical uncertainty, we would expect σ_c2 ≃ 44 ± 6 and σ_c1 × 2 = 75 ± 10 μJy/beam. Both of these values are compatible with our measures and their associated error bars, although, again, the measured value at 1.2 × 2 mm appears slightly lower than the predicted one.

Finally, we compare these estimates with the instrumental noise in the final maps. On the central region of our maps equivalent to one NIKA2 field of view (a disk of 6.5 arcmin diameter) and still masking bright sources, we find total variances of map pixels of about 2.8 × 10⁵ and 4.5 × 10⁴ (μJy/beam)². Subtracting our estimates of confusion (139 and 39 μJy/beam) according to Eq. (B.5), we can derive the residual noise contribution per beam to be 218 and 60 μJy/beam at 1.2 and 2 mm, respectively. It means that the measured confusions are 0.64 and 0.65 times the residual noise in our final maps. These noise estimates agree within 25% with those obtained by Bing et al. (2023) (170 and 48 μJy/beam) with a data reduction optimised for point source detection; thus, it filters out more of the undetected sources that generate the confusion.

This was achieved after about 80 hours of observation of about 64 arcmin², that is, approximately 11 hours per field of view of 33.2 arcmin². These numbers provide with another cross-check when we use them to derive an approximate value of the instrumental sensitivity, also known as the Noise Equivalent Flux Density (NEFD). A thorough assessment of this sensitivity and its uncertainty with these data is beyond the scope of this paper. Still we note that taking the average encountered opacities and elevations (Sect. 2), we obtain effective out-of-atmosphere NEFDs of 32 and 10 mJy/beam at 1.2 and 2 mm, respectively. Although these values are only approximate, they are in good agreement with those presented in Perotto et al. (2020) and comfort our results.

5. Sources responsible for the confusion

The confusion is caused by all galaxies with flux below the detection limit S_lim. However, depending on the exact shape of the number counts, the confusion noise can be dominated by sources with different flux regimes. In this section, we use the SIDES simulation (Béthermin et al. 2017) to predict the contribution of the different flux densities to the confusion, and compare it with the contribution to the CIB, expressed as

$$\begin{array}{cc}\hfill B& ={\displaystyle {\int }_0^{\infty }}S\frac{\mathit{dN}}{\mathit{dS}}dS,\hfill \\ \hfill & ={\displaystyle {\int }_0^{\infty }}{S}^2\frac{\mathit{dN}}{\mathit{dS}}ln(10)d{log}_{10}(S),\hfill \end{array}$$(5)

where $\frac{\mathit{dN}}{\mathit{dS}}$ is the differential number counts. The factor $S^2\frac{\mathit{dN}}{\mathit{dS}}$ thus directly gives the relative contribution of a logarithmic flux density interval to the background.

In the absence of clustering, the confusion noise is expressed as

$$\begin{array}{c}\hfill {\sigma }_c={\displaystyle \int \int b^2}d\mathrm{\Omega }{\displaystyle {\int }_0^{S_{\lim }}}{S}^3\frac{\mathit{dN}}{\mathit{dS}}ln(10)d{log}_{10}(S),\end{array}$$(6)

where b is the beam function. Similar to the background case, the term $S^3\frac{\mathit{dN}}{\mathit{dS}}$ gives the contribution of the logarithmic flux density intervals to the fluctuations. As shown in Béthermin et al. (2017), the clustering tends to broaden the pixel flux distribution and thus increases the confusion. However, even for Herschel/SPIRE at 500 μm with its 36 arcsec beam, the effect is of the order of 15%. We can therefore reasonably neglect it in the case of NIKA2, which has a beam area four times smaller.

In Fig. 4, we show the relative contribution to the background (blue) and the confusion (red) as a function of the flux density. Only sources below the detection limit contribute to the confusion and the grey area indicates the region where they are masked and therefore do not contribute. The curves are calculated from the number counts derived from the SIDES simulation. At both 1.2 mm and 2 mm, the maximum of the red curve is in the masked regime and decreases strongly with decreasing flux density. The sources just below the detection limit are therefore the ones that produce most of the confusion we measure. The contribution to the background peaks at fainter flux densities, which are below the detection limit. However, we note that at 1.2 mm, this peak is just below the detection limit. At this wavelength, the flux density regime that dominates both the background and the confusion is therefore the same.

Fig. 4. Relative contribution to the background (blue curve) and to the fluctuations causing the confusion (red curve) as a function of the flux density derived with the SIDES simulation at 1.2 mm (left) and 2 mm (right) (see Sect. 5). The same quantities derived from the source counts measurements in GOODS-N by Bing et al. (2023) are shown as blue circles and red squares, respectively. The difference between the prediction of galaxy counts from SIDES and the source counts from N2CLS is well understood and caused by the blending of galaxies into the NIKA2 beam (see the analysis in Bing et al. 2023). The grey area shows the flux density regime where the sources are masked and do not contribute to the confusion.

The agreement of the confusion noise from SIDES and N2CLS is an important validation of SIDES in the flux regime corresponding to the peak of the contribution to the background and, thus, to the obscured star formation budget. At 2 mm, the peak is below the detection limit by a factor of ∼3. Therefore, we are mainly probing galaxies that are slightly brighter than those that make up the bulk of the background.

6. Conclusion

From October 2017 to March 2021, we observed the GOODS-North field. We were able to accumulate 78.2 h of observations on about 240 arcmin². By masking the sources that we confidently detected in a previous work (Bing et al. 2023), whose flux is greater than 0.54 or 0.17 mJy at 1.2 and 2 mm, and by restricting ourselves to the 62 arcmin² of our best sampled observations, we were able to derive the first estimates of the confusion ever obtained at the IRAM 30 m telescope. We measured $139.1_{-19.2}^{+15.9}\pm 11.9$ μJy/beam at 1.2 mm, where the first set of error bars is due to noise and the second one due to the cosmic variance. At 2 mm, we measured $38.6_{-13.1}^{+9.6}\pm 3.7$ μJy/beam. To derive these estimates, we developed a new estimator, referred to as the cross-variance. This estimate is less sensitive to residual noise on the maps and also allows us to measure the contribution to the confusion that is correlated between 1.2 and 2 mm. In this case of cross-band confusion, we only found evidence at the level of a 2σ: $49.6_{-24.8}^{+15.9}\pm 6.4$ μJy/beam. All these results are close to those expected from the independently developed SIDES simulations (Béthermin et al. 2017, 2022; Gkogkou et al. 2023). This agreement consolidates both our estimates and the model in the flux regime below the detection limit of N2CLS. Since the confusion is caused by a wide range of fluxes, we cannot formally exclude another scenario where the number counts would have a steeper slope than SIDES near the detection limit, but with a higher number of sources with very faint flux densities. However, many such faint galaxies would be required and this would contradict findings that state that low-mass galaxies tend to have lower dust attenuation and, thus, lower FIR emissions (e.g. Heinis et al. 2014; Fudamoto et al. 2020).

Acknowledgments

We acknowledge financial support from the “Programme National de Cosmologie and Galaxies” (PNCG) funded by CNRS/INSU-IN2P3-INP, CEA and CNES, France, from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (project CONCERTO, grant agreement No 788212) and from the Excellence Initiative of Aix-Marseille University-A*Midex, a French “Investissements d’Avenir” programme. This work is based on observations carried out under project numbers 192-16 with the IRAM 30 m telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). We would like to thank the IRAM staff for their support during the observation campaigns. 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 been sup- ported by the GIS KIDs. This work has benefited from the support of the European Research Council Advanced Grant ORISTARS under the European Union’s Seventh Framework Programme (Grant Agreement no. 291294). A. R. acknowledges financial support from the Italian Ministry of University and Research – Project Proposal CIR01_00010. M. M. E. acknowledges the support of the French Agence Nationale de la Recherche (ANR), under grant ANR-22-CE31-0010. R. A. acknowledges support from the Programme National Cosmology et Galaxies (PNCG) of CNRS/INSU with INP and IN2P3, co-funded by CEA and CNES. R. A. was supported by the French government through the France 2030 investment plan managed by the National Research Agency (ANR), as part of the Initiative of Excellence of Université Côte d’Azur under reference number ANR-15-IDEX-01.

References

Appendix A: Estimation of the confusion via the cross-variance

To motivate the definition of our estimator, we start with the simple case where the signal projected onto maps is unaffected by the data processing and the noise in the map pixels is Gaussian, white, and independent per run. We then return to these hypotheses.

We consider two different runs of observations. The combination of all scans for the first gives the map A, the combination of all scans for the second gives the map B. In the following, capital letters denote map raw values and p and q denote pixel indices. Let S_p be the sky signal at pixel p, so for map A we have A_p = S_p + a_p, where a_p is the noise of map A. We assume that masking bright sources leaves us with N unmasked pixels. Finally, we use an over bar to denote the spatial average of the sky signal over all these pixels:

$$\begin{array}{c}\hfill \overline{S}=\frac{1}{N}{\displaystyle \sum _{p=1}^N}{S}_p,\end{array}$$

With such conventions, the confusion at 1 or 2 mm is expressed as

$$\begin{array}{c}\hfill {\sigma }_c^2=\frac{1}{N}{\displaystyle \sum _p}{(S_p-\overline{S})}^2.\end{array}$$(A.1)

We generalise this definition to the cross-band confusion when A is a 1 mm map and B a 2 mm map.

$$\begin{array}{c}\hfill {\sigma }_{c_{1\times 2}}^2=\frac{1}{N}{\displaystyle \sum _p}(S_p^{1\mathrm{mm}}-{\overline{S}}^{1\mathrm{mm}})(S_p^{2\mathrm{mm}}-{\overline{S}}^{2\mathrm{mm}}),\end{array}$$(A.2)

Under the hypotheses of pure signal and ideal noise, we define the ideal cross-variance estimator as

$$\begin{array}{c}\hfill {\stackrel{\sim }{\rho }}_{\mathit{AB}}=\frac{1}{N}{\displaystyle \sum _p}(A_p-\overline{A})(B_p-\overline{B}),\end{array}$$(A.3)

We show that its statistical expectancy is indeed σ_c² via

$$\begin{array}{cc}\hfill ⟨{\stackrel{\sim }{\rho }}_{\mathit{AB}}⟩& =\frac{1}{N}{\displaystyle \sum _p}⟨(A_p-\overline{A})(B_p-\overline{B})⟩\hfill \\ \hfill & =\frac{1}{N}{\displaystyle \sum _p}⟨(S_p^A-{\overline{S}}^A+a_p)(S_p^B-{\overline{S}}^B+b_p)⟩\hfill \\ \hfill & =\frac{1}{N}{\displaystyle \sum _p}(S_p^A-{\overline{S}}^A)(S_p^B-{\overline{S}}^B)\hfill \\ \hfill & ={\sigma }_c^2\mathrm{or}{\sigma }_{c_{1\times 2}}^2.\hfill \end{array}$$(A.4)

Next, we can compute the covariance of this estimator. By definition, taking four maps, A, B, C, and D, the covariance matrix of $\stackrel{\sim }{\rho }$ is expressed as

$$\begin{array}{cc}\hfill {\stackrel{\sim }{X}}_{AB,CD}& =\langle ({\stackrel{\sim }{\rho }}_{\mathit{AB}}-⟨{\stackrel{\sim }{\rho }}_{\mathit{AB}}⟩)({\stackrel{\sim }{\rho }}_{\mathit{CD}}-⟨{\stackrel{\sim }{\rho }}_{\mathit{CD}}⟩)\rangle \hfill \\ \hfill & =\langle {\stackrel{\sim }{\rho }}_{\mathit{AB}}{\stackrel{\sim }{\rho }}_{\mathit{CD}}-{\stackrel{\sim }{\rho }}_{\mathit{AB}}⟨{\stackrel{\sim }{\rho }}_{\mathit{CD}}⟩-⟨{\stackrel{\sim }{\rho }}_{\mathit{AB}}⟩{\stackrel{\sim }{\rho }}_{\mathit{CD}}+⟨{\stackrel{\sim }{\rho }}_{\mathit{AB}}⟩⟨{\stackrel{\sim }{\rho }}_{\mathit{CD}}⟩\rangle .\hfill \\ \hfill \end{array}$$(A.5)

We define $T_p=S_p-\overline{S}$ for convenience and use Kronecker’s symbol δ, which is equal to 1 when its indices are equal, zero otherwise. The first of the four terms in Eq. (A.5) is expressed as

$$\begin{array}{cc}\hfill ⟨{\stackrel{\sim }{\rho }}_{\mathit{AB}}{\stackrel{\sim }{\rho }}_{\mathit{CD}}⟩& =\frac{1}{N^2}\langle {\displaystyle \sum _{p,q}}(A_p-\overline{A})(B_p-\overline{B})(C_q-\overline{C})(D_q-\overline{D})\rangle \hfill \\ \hfill & =\frac{1}{N^2}\langle {\displaystyle \sum _{p,q}}(S_p^A-{\overline{S}}^A+a_p)(S_p^B-{\overline{S}}^B+b_p)\hfill \\ \hfill & (S_q^C-{\overline{S}}^C+c_q)(S_q^D-{\overline{S}}^D+d_q)\rangle \hfill \\ \hfill & =\frac{1}{N^2}\langle {\displaystyle \sum _{p,q}}(T_p^A+a_p)(T_p^B+b_p)(T_q^C+c_q)(T_q^D+d_q)\rangle \hfill \\ \hfill & =\frac{1}{N^2}\langle {\displaystyle \sum _{p,q}}(T_p^A{T}_p^B+T_p^A{b}_p+T_p^B{a}_p+a_p{b}_p)\hfill \\ \hfill & (T_q^C{T}_q^D+T_q^C{d}_q+T_q^D{c}_q+c_q{d}_q)\rangle \hfill \\ \hfill & =\frac{1}{N^2}{\displaystyle \sum _{p,q}}{\delta }_{\mathit{pq}}[T_p^A{T}_p^B{T}_q^C{T}_q^D+\hfill \\ \hfill & (T_p^A{T}_q^C{\delta }_{\mathit{BD}}+T_p^A{T}_q^D{\delta }_{\mathit{BC}}){\sigma }_{B_p}^2+\hfill \\ \hfill & (T_p^B{T}_q^C{\delta }_{\mathit{AD}}+T_p^B{T}_q^D{\delta }_{\mathit{AC}}){\sigma }_{A_p}^2+\hfill \\ \hfill & ({\delta }_{\mathit{AD}}{\delta }_{\mathit{BC}}+{\delta }_{\mathit{AC}}{\delta }_{\mathit{BD}}){\sigma }_{A_p}^2{\sigma }_{B_p}^2].\hfill \end{array}$$(A.6)

In the previous equation, σ_Ap² for instance, denotes the variance of the noise a_p in pixel p of map A. It is derived at the map making stage when all scans are co-added (Sect. 3.2) and will be responsible for the statistical uncertainty on the determination of the confusion.

Assuming for now that A, B, C, and D are all maps of the same band, then Eq. (A.6) becomes:

$$\begin{array}{cc}\hfill ⟨{\stackrel{\sim }{\rho }}_{\mathit{AB}}{\stackrel{\sim }{\rho }}_{\mathit{CD}}⟩& =\frac{1}{N^2}{\displaystyle \sum _{\mathit{pq}}}{\delta }_{\mathit{pq}}[T_p^2{T}_q^2\hfill \\ \hfill & +T_p^2[({\delta }_{\mathit{BD}}+{\delta }_{\mathit{BC}}){\sigma }_{B_p}^2+({\delta }_{\mathit{AD}}+{\delta }_{\mathit{AC}}){\sigma }_{A_p}^2]\hfill \\ \hfill & +({\delta }_{\mathit{AD}}{\delta }_{\mathit{BC}}+{\delta }_{\mathit{AC}}{\delta }_{\mathit{BD}}){\sigma }_{A_p}^2{\sigma }_{B_p}^2\hfill \end{array}$$(A.7)

The other terms involved in Eq. (A.5) are

$$\begin{array}{cc}\hfill \langle {\stackrel{\sim }{\rho }}_{\mathit{AB}}⟨{\stackrel{\sim }{\rho }}_{\mathit{CD}}⟩\rangle & =⟨{\stackrel{\sim }{\rho }}_{\mathit{AB}}⟩\frac{1}{N}{\displaystyle \sum _q}{T}_q^C{T}_q^D\hfill \\ \hfill & =\frac{1}{N}{\displaystyle \sum _p}\langle T_p^A{T}_p^B+T_p^A{b}_p+T_p^B{a}_p+a_p{b}_p\rangle \hfill \\ \hfill & \times \frac{1}{N}{\displaystyle \sum _q}{T}_q^C{T}_q^D\hfill \\ \hfill & =\frac{1}{N^2}{\displaystyle \sum _{\mathit{pq}}}{T}_p^A{T}_p^B{T}_q^C{T}_q^D.\hfill \end{array}$$(A.8)

The same goes for the last two terms:

$$\begin{array}{c}\hfill \langle ⟨{\stackrel{\sim }{\rho }}_{\mathit{AB}}⟩{\stackrel{\sim }{\rho }}_{\mathit{CD}}\rangle =\frac{1}{N^2}{\displaystyle \sum _{\mathit{pq}}}{T}_p^A{T}_p^B{T}_q^C{T}_q^D=⟨{\stackrel{\sim }{\rho }}_{\mathit{AB}}⟩⟨{\stackrel{\sim }{\rho }}_{\mathit{CD}}⟩,\end{array}$$

hence,

$$\begin{array}{cc}\hfill {\stackrel{\sim }{X}}_{AB,CD}^{1\times 1\mathrm{or}2\times 2}& =\frac{1}{N^2}{\displaystyle \sum _p}{T}_p^2[({\delta }_{\mathit{AC}}+{\delta }_{\mathit{AD}}){\sigma }_{A_p}^2+({\delta }_{\mathit{BC}}+{\delta }_{\mathit{BD}}){\sigma }_{B_p}^2]\hfill \\ \hfill & +({\delta }_{\mathit{AC}}{\delta }_{\mathit{BD}}+{\delta }_{\mathit{AD}}{\delta }_{\mathit{BC}}){\sigma }_{A_p}^2{\sigma }_{B_p}^2.\hfill \end{array}$$(A.9)

In the specific case of the cross-band confusion, because we take maps from different runs and bands, if A and C refer to maps at 1 mm and B and D to maps at 2 mm, then eq. (A.9) simplifies further into

$$\begin{array}{cc}\hfill {\stackrel{\sim }{X}}_{AB,CD}^{1\times 2}& =\frac{1}{N^2}{\displaystyle \sum _p}{T}_p^2(1\mathrm{mm}){\delta }_{\mathit{BD}}{\sigma }_{B_p}^2+T_p^2(2\mathrm{mm}){\delta }_{\mathit{AC}}{\sigma }_{A_p}^2\hfill \\ \hfill & +{\delta }_{\mathit{AC}}{\delta }_{\mathit{BD}}{\sigma }_{A_p}^2{\sigma }_{B_p}^2,\hfill \end{array}$$(A.10)

We note that in each case, the covariance matrix involves terms like T_p², that is to say the pure underlying signal, that we actually don’t know in real data. In practice, we replace this term by the expected confusion squared σ_c². Because it is in any case much smaller than the noise per pixel, the uncertainty on its exact value has a negligible impact on the derivation of the covariance matrix and subsequently that of the final error bar.

In a more realistic case, the data are inevitably affected by the data processing. It is common to loosely refer to the overall effect as "filtering". In this work, where we determine a single scalar term, we model the impact of the filtering by a multiplicative factor f on the confusion, smaller than 1, that needs to be determined. Thus, our cross-variance estimator and its covariance become:

$$\begin{array}{cc}\hfill {\rho }_{\mathit{AB}}& =\frac{1}{f^2}{\stackrel{\sim }{\rho }}_{\mathit{AB}}\hfill \end{array}$$(A.11)

$$\begin{array}{cc}\hfill X_{AB,CD}& =\frac{1}{f^4}{\stackrel{\sim }{X}}_{AB,CD}\hfill \end{array}$$(A.12)

In this work, the timeline processing proceeds per subscan, and they have different orientations on the sky (cf. Sect. 2). As discussed in Sect. 4.2, this leads to different values of f. Filtering also affects the noise. In general, filtering mostly affects low frequency components and therefore has less effect on the statistical noise than on the signal. The 1/f² correction in the definition of X_AB, CD is therefore somewhat conservative in terms of statistical uncertainty.

Figure A.1 shows the 66 cross-variance estimates at 1.2 mm, and Fig. A.2 shows their associated covariance matrix.

Fig. A.1. Estimates of the cross-variance at 1.2 mm for all pairs of independent observation runs in (μJy/beam)2. Error bars are the square root of the diagonal elements of the covariance matrix of ρi shown on Fig. A.2. The dashed line marks zero and the solid line is the combined cross-variance (see Eq. (A.14)).

Fig. A.2. Covariance matrix of the 66 cross-variance estimates. Each pixel of this image is an element of this 66 × 66 matrix, and its value is given by the color bar. The color scale is saturated on purpose to display the off-diagonal terms. The top plot shows the diagonal elements of the matrix. The blocks, dots and lines that appear in the matrix reflect the correlations when the same map is used in Eq. (A.9).

Although we prefer to use likelihoods to derive our estimates of confusion and their associated confidence intervals, we here provide a direct way to estimate the overall cross-variance for the sake of completeness. We have checked that in the high-S/N case of 1.2 mm, the following derivation matches the maximum likelihood value given in the text.

We can gather all our cross-variance estimates in a single vector {ρ_i}_i and relabel the covariance matrix X_ij terms accordingly. The optimal way to combine ρ_i derives from the minimisation of the log likelihood

$$\begin{array}{c}\hfill -2ln\mathcal{L}={(\rho -\widehat{\rho })}^T{X}^{-1}(\rho -\widehat{\rho }),\end{array}$$(A.13)

thus,

$$\begin{array}{ccc}\hfill \frac{\partial (-2ln\mathcal{L})}{\partial \widehat{\rho }}=0& \iff \hfill & \hfill \frac{\partial }{\partial \widehat{\rho }}{(\rho -\widehat{\rho })}^T{X}^{-1}(\rho -\widehat{\rho })=0\\ \hfill & \iff \hfill & \hfill \frac{\partial }{\partial \widehat{\rho }}{\displaystyle \sum _{\mathit{ij}}}{\rho }_i{X}_{\mathit{ij}}^{-1}{\rho }_j-{\rho }_i{X}_{\mathit{ij}}^{-1}\widehat{\rho }-\widehat{\rho }{X}_{\mathit{ij}}^{-1}{\rho }_j+{\widehat{\rho }}^2{X}_{\mathit{ij}}^{-1}=0\\ \hfill & \iff \hfill & \hfill -2{\displaystyle \sum _{\mathit{ij}}}{X}_{\mathit{ij}}^{-1}{\rho }_j+2\widehat{\rho }{\displaystyle \sum _{\mathit{ij}}}{X}_{\mathit{ij}}^{-1}=0\\ \hfill & \iff \hfill & \hfill \widehat{\rho }=\frac{{\sum }_{\mathit{ij}}{X}_{\mathit{ij}}^{-1}{\rho }_j}{{\sum }_{\mathit{ij}}{X}_{\mathit{ij}}^{-1}},\end{array}$$(A.14)

and the final uncertainty of $\widehat{\rho }$ reads

$$\begin{array}{c}\hfill {\sigma }_{\widehat{\rho }}^2=\frac{1}{{\sum }_{\mathit{ij}}{X}_{\mathit{ij}}^{-1}}.\end{array}$$(A.15)

Appendix B: Note on photometry and conventions

The NIKA2 maps used in this work are calibrated in mJy/beam, and the absolute calibration is performed following Perotto et al. (2020). In this convention, a point source of flux ϕ would show up on the map m_p as a Gaussian with an amplitude equal to ϕ, i.e. m_p = ϕg_p, with g_p normalised such that their maximum is 1. In the limit of no noise and if the source falls exactly at the center of a map pixel, one could just read the value of this pixel to determine the flux. In practice, a weighted average of the pixel values m_p is required, and a simple estimator of the flux of the source is:

$$\begin{array}{c}\hfill \widehat{\varphi }=\frac{1}{{\sum }_p{g}_p^2}{\displaystyle \sum _p}{g}_p{m}_p\end{array}$$(B.1)

Under the assumption that the noise per pixel σ_p is independent, the variance of this estimator is expressed as

$$\begin{array}{c}\hfill \mathrm{Var}(\widehat{\varphi })={\left(\frac{1}{{\sum }_p{g}_p^2}\right)}^2{\displaystyle \sum _p}{g}_p^2{\sigma }_p^2.\end{array}$$(B.2)

Thus, under the assumption of uniform noise per pixel σ, the uncertainty on the flux of the source is simply:

$$\begin{array}{c}\hfill {\sigma }_{\varphi }=\frac{\sigma }{\sqrt{{\displaystyle \sum _p}{g}_p^2}}.\end{array}$$(B.3)

We can now relate the instrumental noise per beam to the overall variance of the unmasked pixels of our maps. Let us use the same notation as in the text to denote the value of a pixel A_p, its signal content S_p and its noise a_p, then:

$$\begin{array}{cc}\hfill {\sigma }^2& =\frac{1}{N}{\displaystyle \sum _{p=1}^N}{(A_p-\overline{A})}^2\hfill \\ \hfill & =\frac{1}{N}{\displaystyle \sum _{p=1}^N}{(S_p-\overline{S}+a_p-\overline{a})}^2\hfill \\ \hfill & =\frac{1}{N}{\displaystyle \sum _{p=1}^N}{(S_p-\overline{S})}^2+\frac{1}{N}{\displaystyle \sum _{p=1}^N}{(a_p-\overline{a})}^2\hfill \\ \hfill & ={\sigma }_c^2+{\sigma }_a^2.\hfill \end{array}$$(B.4)

According to Eq. (B.3), the instrumental noise per beam is now expressed as

$$\begin{array}{c}\hfill {\sigma }_{N/beam}=\sqrt{\frac{({\sigma }^2-{\sigma }_c^2)}{{\displaystyle \sum _p}{g}_p^2}}\end{array}$$(B.5)

Computing σ² on our final maps and using the values of the confusion derived in Sect 4.3, i.e. 139 and 39 μJy/beam, we thus find σ_N/beam = 218 and 60 μJy/beam at 1.2 and 2 mm, respectively.

Appendix C: The KID Möbius circle calibration method

In a nutshell, a KID is monitored by its response to two excitation signals, one in phase I and the other one in quadrature Q. The KID tone, set at frequency f is modulated by a known Δf frequency shift, every half millisecond, that induces variations of these excitation signals by those measured quantities dI and dQ (Catalano et al. 2014). Calvo et al. (2013) have shown different methods to use these quantities to monitor the variations of the resonance frequency f₀ of a KID when it collects photons, and how to relate them to the incoming power. We here present yet another method, that proves to be more robust against numerical integration. It relies on 1D (instead of 2D) polynomial fit. It is also easier to implement in algorithmic terms. The idea is to find how to project I, Q and dI, dQ onto an axis in a way that is as linear as possible with frequency, itself shown to be linear with optical power (Monfardini et al. 2014). This is provided by the Möbius transform.

We study the expected physical Z = I + jQ complex dependency of a KID with frequency. Grabovskij et al. (2008) have shown that

$$\begin{array}{c}\hfill Z=\frac{Z_r{z}_1}{Z_r{z}_2+z_3},\end{array}$$(C.1)

where z₁, z₂ and z₃ are constant complexes, and $Z_r=\frac{1}{2}+j(f-f_0)/w$ where f₀ and w are real values in Hz. The location of Z is on a circle, at least approximately when near the resonance (see Fig. C.2.2 and the subsection below on Circle Fit). The inverse of Z is a circle, but we can transform it to have an infinite radius circle i.e. a line, which can be expected to be linearised with the KID frequency. To simplify the solution of this inversion, we proceed in two steps. The first step is to scale, translate, rotate and reverse the initial circle so that it is identical to the reference circle defined as a $\frac{1}{2}$ radius, and $(\frac{1}{2},0)$ centre, in the complex plane, as defined by:

$$Z_{\mathit{ref}}=\frac{1}{2}+\frac{1}{2}[cos\varphi +jsin\varphi ]=cos\frac{\varphi }{2}expj\frac{\varphi }{2}.$$

The second step consists in inverting that reference circle (the Möbius transform), as $\frac{1}{Z_{\mathit{ref}}}=Z_{\mathit{res}}$ with

$$Z_{\mathit{res}}=\frac{1}{cos\frac{\varphi }{2}}exp-j\frac{\varphi }{2}=1-jtan\frac{\varphi }{2}=x_3+j{y}_3,$$

and $f=f_0+\frac{w}{2}tan\frac{\varphi }{2}$. The imaginary part of the inversion of the reference circle is just linearly dependent on the frequency, to first order. To calibrate this dependency (because the linearity is only approximate) we must rely on the dI, and dQ measurements (see the subsection below on 1D polynomial fit).

Circle fit

The first step is best achieved with a circle fit to the data samples⁵. Then, if r is the radius of the circle and x_c,  y_c its centre and $\alpha =arctan(\frac{y_c}{x_c})$, the complex Z becomes Z_norm with $I_{\mathit{norm}}=-\frac{1}{2r}[(I-x_c)cos\alpha +(Q-y_c)sin\alpha ]+\frac{1}{2}$, and $Q_{\mathit{norm}}=\frac{1}{2r}[-(I-x_c)sin\alpha +(Q-y_c)cos\alpha ]$. This is shown in Fig. C.1.1 where the top-left plot shows Z and the top-right plot shows Z_norm. This figure is made with a simple simulation using Eq. (C.1) and including some noise. The bottom-left plot shows the inverse of Z_norm. The data points lie mostly along the imaginary line of abscissa unity and of ordinate y₃. The transform from mostly circular data on a plane to a one dimensional data stream is a particular case of Möbius transforms. They preserve angles, and this property is used in the following to propagate the calibration from the circle to the imaginary line.

Fig. C.1. Illustration of our use of a Möbius transform to derive the total power measured by a KID on a simulation. Top-left: black diamonds are individual I,  Q points on a circle (blue curve). Top-right: normalised circle transformation of the points into the reference circle coordinates. The resonance is reached at (1,0). Bottom-left: Complex inverse of the points in the previous complex plane (the Möbius transform). Bottom-right: 1D polynomial fit (blue curve) to the derivative data points corresponding to dI,  dQ normalised, inverted and selected on the imaginary axis.

1D polynomial fit

We need to calibrate y₃ which is nearly proportional to f − f₀. The modulation technique gives us this calibration via the derivative function. Let us call dI, and dQ, the measured variations of the signal with a modulation of frequency of Δf. The normalised variations are $d{I}_{\mathit{norm}}=-\frac{1}{2r}[dIcos\alpha +dQsin\alpha ]$, and $d{Q}_{\mathit{norm}}=\frac{1}{2r}[-dIsin\alpha +dQcos\alpha ]$. The derivative of the inverse complex is $d{Z}_{\mathit{res}}=-d{Z}_{\mathit{norm}}/Z_{\mathit{norm}}^2$. We then take the imaginary part of dy₃ and use it to calibrate y₃ in this way: we fit $\frac{\mathrm{\Delta }f}{d{y}_3}=R_n(y_3)$ where R_n is a polynomial of degree n (see bottom-right panel of Fig. C.1.1) which is easy to integrate into P_n + 1 (with Ṗ_n+1 = R_n). We then obtain the relative frequency of the KID from y₃ only with f − f₀ = P_n + 1(y₃) which is plotted in Fig. C.2.2 along with the residual of the fit.

Fig. C.2. Comparison between the KID simulation and the 1D polynomial retrieval. Relative KID frequency as a function of the phase ϕ in the circle, as deduced from the 1D polynomial fit (grey curve) and from the initial simulation (red dashed curve). The small wiggly horizontal black line represents the residual frequency (grey curve minus red curve) multiplied by a factor ten for visibility, as a function of the same phase.

The method is robust if the polynomial fit is weighted. Indeed, the noise on the derivative measurement ($\frac{\mathrm{\Delta }f}{d{y}_3}$) increases as the phase changes (see the bottom-right panel of Fig. C.1.1 where the flaring away from the zero phase is apparent). The noise is proportional to the noise in dy₃ divided by $d{y}_3^2$. The noise in dy₃ is not constant. It can be shown that it is inversely proportional to the distance of the current I, Q point to the point in the circle which is opposite to the resonance point. This noise model can also be used to show that the noise on the recovered frequency is flaring as we go away from the resonance. To a good approximation, the rms noise is proportional to $1+{tan}^2\frac{\varphi }{2}=1+{(\frac{2(f-f_0)}{w})}^2$. Hence, although a KID is a linear device, its noise depends on the distance of the measuring tone frequency to the resonance frequency.

This Möbius Circle calibration method is better at accounting for the quick variations of dI, dQ and as such, is even more linear than previous methods. Another benefit of the method is that the sky noise, mostly due to water vapour fluctuations in the first kilometres of the atmosphere, is calibrated more accurately between KIDs, on the long time scales, whereas other methods have a random drift $\sim \sqrt{t}$. Here, the global circle method ensures that the measured KID frequency is unique for a given point in the [I, Q] plane. The atmospheric common mode can thus be removed with less residual noise.

All Tables

All Figures

	Fig. 1. Maps of GOODS-North used in this work (left 1.2 mm, right 2 mm). The black contours define the region (i.e. mask) that we use to compute the confusion.
In the text

Fig. 2. Variance of the N2CLS maps vs. the out of atmosphere integration time on a region in the center equivalent to one NIKA2 field of view (6.5 arcmin diameter). To determine these data points, we simply take the variance of the unmasked pixels (see Fig. 1) in the central area. For each NIKA2 band, the lower plot shows the residuals compared to the t−1/2 instrumental noise integration that is fitted here for times of integration smaller than 20 000 s. The error bars are the square roots of the diagonal elements of the covariance matrix of the measures (Fig. A.2), derived from 100 Monte-Carlo simulations of map-based realistic noise (per scan). There is a clear excess of signal at high integration time due to the confusion, but as explained in the main text, this plot is for illustrative purposes only and is not used to derive our final results.

In the text

	Fig. 3. Likelihoods (normalised to their maximum) of the confusion at 1 mm, 2 mm, and the cross-confusion 1 × 2 mm for different subsets of our data. The central points and their 1-σ equivalent error bars are shown for convenience and have arbitrary ordinates. ‘Half 1’ refers to scans oriented at 40°. ‘Half 2’ refers to scans oriented at −50°. ‘All’ refers to the total sum of the two.
In the text

Fig. 4. Relative contribution to the background (blue curve) and to the fluctuations causing the confusion (red curve) as a function of the flux density derived with the SIDES simulation at 1.2 mm (left) and 2 mm (right) (see Sect. 5). The same quantities derived from the source counts measurements in GOODS-N by Bing et al. (2023) are shown as blue circles and red squares, respectively. The difference between the prediction of galaxy counts from SIDES and the source counts from N2CLS is well understood and caused by the blending of galaxies into the NIKA2 beam (see the analysis in Bing et al. 2023). The grey area shows the flux density regime where the sources are masked and do not contribute to the confusion.

In the text

	Fig. A.1. Estimates of the cross-variance at 1.2 mm for all pairs of independent observation runs in (μJy/beam)2. Error bars are the square root of the diagonal elements of the covariance matrix of ρi shown on Fig. A.2. The dashed line marks zero and the solid line is the combined cross-variance (see Eq. (A.14)).
In the text

Fig. A.2. Covariance matrix of the 66 cross-variance estimates. Each pixel of this image is an element of this 66 × 66 matrix, and its value is given by the color bar. The color scale is saturated on purpose to display the off-diagonal terms. The top plot shows the diagonal elements of the matrix. The blocks, dots and lines that appear in the matrix reflect the correlations when the same map is used in Eq. (A.9).

In the text

Fig. C.1. Illustration of our use of a Möbius transform to derive the total power measured by a KID on a simulation. Top-left: black diamonds are individual I,  Q points on a circle (blue curve). Top-right: normalised circle transformation of the points into the reference circle coordinates. The resonance is reached at (1,0). Bottom-left: Complex inverse of the points in the previous complex plane (the Möbius transform). Bottom-right: 1D polynomial fit (blue curve) to the derivative data points corresponding to dI,  dQ normalised, inverted and selected on the imaginary axis.

In the text

Fig. C.2. Comparison between the KID simulation and the 1D polynomial retrieval. Relative KID frequency as a function of the phase ϕ in the circle, as deduced from the 1D polynomial fit (grey curve) and from the initial simulation (red dashed curve). The small wiggly horizontal black line represents the residual frequency (grey curve minus red curve) multiplied by a factor ten for visibility, as a function of the same phase.

In the text