3.1. Simulated maps

In preparation for testing the model set-up, we simulate SPIRE observations using the full galaxy catalogues from SIDES. The SIDES simulation is a $2 \, \mathrm{deg}^2$ mock galaxy catalogue at FIR to submillimetre wavelengths over the redshift range $0\lt z\lt7$ . Galaxies in SIDES are clustered through abundance matching of a dark matter simulation. While more recent simulations incorporating gas dynamics and dust modelling are now available, SIDES remains well-suited for our purposes, since the key requirement here is realistic clustering from dark matter halos rather than detailed baryonic physics. Figure 12 in Béthermin et al. (Reference Béthermin2017) illustrates this clustering effect in SIDES, showing the difference between the pixel histogram of the map with clustered galaxy positions and one with galaxy positions randomly distributed. In that paper it is thus proposed that spatial clustering of galaxies can bias the P(D) results of galaxy counts. The effect should be most significant in the SPIRE 500- $\mu$ m observations, given the larger beam at the longer wavelengths.

To illustrate the direct impact on the counts fit from clustering, we generate simulated Herschel-SPIRE 500- $\mu$ m maps in a similar manner to Béthermin et al. (Reference Béthermin2017). We obtain the full galaxy catalogues with simulated Herschel-SPIRE 500- $\mu$ m flux densities from SIDES to generate the clustered map. We then randomise the positions in the catalogue and generate another map, hereafter called a ‘randomised map’. We adopt a pixel size of 12 $\, \mathrm{arcsec}$ in the SPIRE 500- $\mu$ m maps. Both maps are convolved with a Gaussian beam of 36.3 arcsec FWHM, which corresponds to the SPIRE beam at 500 $\mu$ m. We check that the confusion noise ( $6.8$ mJy at this step) is consistent with the measurements of Nguyen et al. (Reference Nguyen2010). We then assume Gaussian instrumental noise with $\sigma = 1.03\, \mathrm{mJy}\, \mathrm{beam^{-1}}$ ,Footnote ^b which is the typical noise value in the central region of the deep GOODS-N field mapped by Herschel-SPIRE at 500 $\,\unicode{x03BC}$ m (see Section 5 for details about GOODS-N). We finally subtract the mean, since Herschel is not sensitive to the absolute value in its observations. The image edges are slightly trimmed by 5 pixels to avoid edge effects in the convolution step.

3.2. Model parameterisation

To enable Bayesian inference through P(D) to fit underlying galaxy number counts, we parameterise the galaxy number counts using a node-based model with interpolation of multiple power laws; in other words we have a function that is composed of straight-line segments on a log-log diagram. See Vernstrom et al. (Reference Vernstrom2014) for a discussion of the advantages of this approach compared with other parameterisations.

The positions of the nodes are selected to smoothly model the galaxy number counts. Number counts outside the node range are assumed to be zero. As a result, the faint limit of the nodes is particularly important and needs to be treated on a case-by-case basis, depending on the beam and noise of the data of interest. The first node is often treated as an upper limit in previous studies, since at this flux-density level galaxies do not change the shape of P(D) anymore, but only move the mean of P(D). While the mean can be calculated in the model set-up stage, we note later in this work that the mean (or equivalently the peak position of the histogram) should be fit in order to properly apply the clustering correction (see Section 4.1), and hence the mean is not computed in our model and the first node remains an upper limit.

In practice, we choose the first node using the full galaxy catalogue in the simulation and simulate maps of the beam and the noise level for SPIRE (see Section 3.4). We then calculate P(D) with the catalogue counts sampled in nodes and compare P(D) with the histogram of the simulated maps. We determine the first node where the map histogram is consistent with the P(D) calculated from the truncated counts, including the mean of P(D), i.e. the data are sensitive down to this flux-density limit. The second node is chosen to be at around a quarter of the instrumental noise level, following Vernstrom et al. (Reference Vernstrom2014). The brightest node is at around $0.1$ Jy, above which the number counts become well constrained by directly counting the resolved sources in SPIRE maps (e.g. Oliver et al. Reference Oliver2010; Clements et al. Reference Clements2010; Béthermin et al. Reference Béthermin2012; Valiante et al. Reference Valiante2016). The nodes between the second node and the brightest node are spaced logarithmically. The number of nodes is a trade-off between a full characterisation of the galaxy number counts model shape and computational speed, with too many nodes yielding highly-correlated results. After computing tests on different numbers of nodes on simulated SPIRE maps (see Section 3.4), we find that eight or nine nodes are sufficient to characterise the counts in SPIRE maps.

We apply some conditions on the model to ensure a physical parameterisation, as follows. (1) The integrated source count flux density $\int S\frac{dN}{dS} dS$ should be consistent with the upper limit of the total estimated cosmic infrared background (CIB) at 250, 350 and 500 $\,\mathrm{\unicode{x03BC} m}$ from Fixsen et al. (Reference Fixsen, Dwek, Mather, Bennett and Shafer1998). (2) The slope should not be steeper than -6 between adjacent nodes on a log-log scale. (3) The difference between two adjacent slopes should not be greater than 3.5, to maintain relatively smooth interpolation. (4) The magnitude of the faintest node should be between $10^9$ and $10^{16}$ Jy $^{-1}$ sr $^{-1}$ , a weak constraint based on observed galaxy number counts at these wavelengths (Marsden et al. Reference Marsden2009; Glenn et al. Reference Glenn2010; Béthermin et al. Reference Béthermin2012).

Figure 1.

Standard deviation by bin in the histograms of 1 000 randomised maps as a ratio of the Poisson error estimate. Randomised maps are simulated at the original pixel size and also with pixels sampled at $1\times$ FWHM. A grey band spanning 10% of the Poisson error is illustrated for visualisation. The sampled map is better represented by Poisson error, while the original map has additional scatter, mostly from pixel correlations.

We use the MCMC ensemble sampler emcee (Goodman & Weare Reference Goodman and Weare2010; Foreman-Mackey et al. Reference Foreman-Mackey, Hogg, Lang and Goodman2013) to fit the multinode model. For each step in the MCMC sampling, the model is evaluated by generating its P(D) prediction in Equation (3), and then comparing with the histogram of the input map to define the likelihood. We choose to estimate the log-likelihood as in Vernstrom et al. (Reference Vernstrom2014):

(4) \begin{equation}\log \mathcal{L}=-\frac{1}{2} \sum_{i} \frac{\left(D_{i}-M_{i}\right)^{2}}{\sigma_{i}^{2}},\end{equation}

where $D_i$ is the intensity in the ith bin, $\sigma_i$ is the error on $D_i$ , and $M_i$ is the P(D)-modelled histogram in the ith pixel for a given counts model. We assume the uncertainty to be Poisson distributed. We discuss reducing the uncertainty from pixel correlations tested on simulated maps in Section 3.3, and from the discretisation process of binning in Section 3.4.

3.3. Pixel correlations

The likelihood we adopt in Equation (4) assumes Poisson errors. However, the SPIRE beams are larger than the pixel area, where the pixel width is roughly a third of the FWHM of the beam by default in map making (Oliver et al. Reference Oliver2012). The unresolved sources can therefore contribute to multiple pixels, generating correlations between pixels. Ideally, the error in the likelihood should be evaluated using the full covariance matrix of pixel intensities. However, the exact form of the covariance matrix is challenging to determine in real data. Although it is feasible to obtain the covariance matrix through simulations (e.g. Glenn et al. Reference Glenn2010), the underlying counts model adopted in the simulation will have an effect (Vernstrom et al. Reference Vernstrom2014), as well as non-idealities in the data. Following Vernstrom et al. (Reference Vernstrom2014), we avoid introducing potential biases from the covariance matrix calculated from the simulations, and instead test the performance of sampling maps at larger pixel sizes.

We first explore the effect of pixel correlations contributing to the histogram bin error in simulations using the full galaxy catalogues from SIDES at 500 $\unicode{x03BC}$ m, where the beam is the largest. We use the 1 000 simulated maps with positions of galaxies randomised, as discussed in Section 3.1, with the pixel width the same as in the observed 500- $\unicode{x03BC}$ m map (12 arcsec $\,\simeq\,$ FWHM/3), and obtain an average histogram using the binning mentioned in Section 3.4. The error in each histogram bin can be described by the standard deviation of the number of pixels in the bin in 1 000 simulated maps, $\sigma_{\text{sim}}$ . On the other hand, we can calculate the Poisson error by taking the square-root of the number of pixels in the ith bin, $\sigma_{\text{Poisson}} = \sqrt{N_i}$ . A comparison between the two should tell us how well the error is described by the Poisson distribution. We repeat the same process for simulated maps sampled at different pixel sizes. In Figure 1, we specifically plot $\sigma_{\text{sim}}/\sigma_{\text{Poisson}}$ in each bin for the original pixel size (FWHM/3) and for a pixel size equal to the FWHM. The reference line is $\sigma_{\text{sim}}/\sigma_{\text{Poisson}}=1$ , where the Poisson error characterises well the error seen in the simulations. In the case of the original pixel width, we find that the overall simulated error is larger than the Poisson error, suggesting that there are underlying correlations contributing on top of the Poisson errors, especially at the bright and faint ends, while the larger sampled pixel size behaves consistently around $\sigma_{\text{sim}}/\sigma_{\text{Poisson}} \simeq 1$ . We therefore use $1\times$ FWHM as the pixel size throughout the work, for all three frequencies of Herschel-SPIRE data.

3.4. P(D) fitting

We implement the P(D) fit of the node-based model (Section 3.2) on the histograms of both the clustered and randomised maps. The numerical integration of the P(D) prediction in Equation (2) is calculated in $2^{18}$ flux density bins in the map flux range. In Equation (4), the likelihood is evaluated at each histogram bin, so that the choice of bins is important to trace the statistical significance in the comparison between data and model. We use the simulated clustered map as a reference to create the binning. The map is first sampled at the FWHM pixel width to reduce pixel correlations (Section 3.3). We make extra fine ( $2^{10}$ – $2^{12}$ ) linearly uniform bins across the flux density range of the sampled map. We then adaptively merge the fine bins to ensure more than 20 pixels in each bin for Poisson statistics to be reliable. When the number of pixels in a bin is smaller than a fraction of the number of pixels in the most populated bin (peak fraction $1/70$ – $1/50$ ), we switch to log-spaced bins. In practice we cut off the histogram where the adaptive bins become log-spaced, since we find that wider bins are not well characterised by their central fluxes and may bias the likelihood.

To ensure our results are robust to the choice of histogram binning, we perform a convergence test on the bin resolution. We vary the adaptive binning parameters (increasing the resolution of the initial linear grid and reducing the peak-fraction threshold) while maintaining the requirement of $N\gt20$ counts per bin to preserve Gaussian statistics. We found that beyond a certain resolution, the inferred model parameters and uncertainties stabilised and became independent of the binning scheme. We adopt the finer binning configuration as the optimal setup for our analysis to minimise any discretisation errors, while ensuring the binning is limited only by the map noise. We note that this adaptive binning is also treated case-by-case. In the following sections when we deal with maps of different size, beam, or noise, we rerun the adaptive binning process to find the optimal binning to reliably evaluate the model likelihood. Finally, we fix the brightest node to the simulation value, which is obtained through binning the simulation catalogue into fine flux-density bins and interpolating at that node.

Figure 2.

Left panel: Observed histograms of the clustered map (blue) and the randomised map (red), both generated from the simulated Herschel-SPIRE map at 500 $\,\unicode{x03BC}$ m. Right panel: P(D) galaxy number counts fit to the clustered map (blue) and the randomised map (red). The galaxy number counts are Euclidean normalised here. The first node is an upper limit. The shaded region is the 68% confidence region for each fit. The black dots are the real galaxy counts from the SIDES catalogue with Poisson errors in a finer binning for illustration.

Since the maps are not sensitive to the mean brightness in the field, the observed histograms are shifted compared to the P(D) predictions. We simply calculate the mean in each P(D) prediction and shift the P(D) histogram by this mean. We will show later that the clustering bias effectively changes the mean flux in the map, therefore fixing the mean flux to the value predicted by P(D) would not be desirable for clustering corrections (Section 4.1). For a bin-to-bin comparison, the P(D) prediction is interpolated onto the flux-density bins of the observed histogram defined in the last paragraph. We then evaluate the bin-to-bin difference using the likelihood function (Equation 4).

For the MCMC sampling process, the initial state of the positions of the walkers uses the interpolated nodes on the real galaxy counts in the simulation catalogue with the addition of a Gaussian scatter to generate some randomness in the initial states for different walkers that satisfy the conditions mentioned in Section 3. The problem becomes multimodal for different proposed models, and so we adopt a weighted move, with the differential evolution proposal described in Nelson et al. (Reference Nelson, Ford and Payne2014) of 80% probability and that with the ‘snooker updater’ in Ter Braak & Vrugt (Reference Ter Braak and Vrugt2008) of 20% probability.

The results of the fits on both the clustered map and the randomised map are shown in Figure 2. The left panel shows the histograms, and the right panel shows the corresponding galaxy counts models that generate the best fit of the histograms in the left panel. The 68% confidence region at each node is derived by obtaining the 16th and 84th percentiles of the samples after burn-in. The error bars in regions between the nodes are derived by interpolating 1 000 uniform bins on the samples of the nine fit nodes to the whole flux range in each step. We thereby confirm that the galaxy counts fit is biased when directly applying P(D) fluctuation analysis to the clustered map.

4.1. Generating maps of different clustering strength

Galaxy clustering can be characterised by the 2-point correlation function or the power spectrum of the map. Throughout this section, we discuss the effect of galaxy clustering only from 2-point statistics, since it is the main contributor to the change of the flux histogram compared to any higher-order correlations. In 2-dimensional source catalogues the 2-point correlation function $w(\theta)$ is defined as the excess probability of finding a second galaxy at an angular separation $\theta$ from a given galaxy. At small angular scales, it is empirically approximated by a power law of the form

(5) \begin{equation}w(\theta)=\left(\theta / \theta_{0}\right)^{-\alpha},\end{equation}

where $\theta_0$ is the characteristic clustering scale and thus the parameterisation of the clustering strength. We assume that the power-law index is a constant, $\alpha = 0.8$ , which has been found for a large number of optical galaxy surveys starting with Davis & Peebles (Reference Davis and Peebles1983), and is also found for radio galaxies (Blake & Wall Reference Blake and Wall2002a,Reference Blake and Wallb). We expect modest deviations from this assumption not to significantly affect the correction form and we note that this is a simplification of actual galaxy clustering, and we will discuss the potential effect of this simplification in Section 6.6.

We generate modified maps from the SIDES galaxy catalogue with different characteristic clustering scales $\theta_0$ , uniformly ranging from $0.036$ to $1.38\, \mathrm{arcsec}$ , using a method that we describe below. The range of $\theta_0$ is determined so that it is approximately centred at the actual clustering of the 2 $\, \mathrm{deg}^2$ simulated map from SIDES in the same configuration. This clustering scale is found by comparing the histograms of the generated clustered maps and the actual SIDES map. For simplicity, we assume $\theta_0$ to be the same for all galaxies, regardless of their brightness.

Figure 3.

Examples of generated maps with different clustering strengths or characteristic clustering scales $\theta_0$ (Equation 5), with exaggerated values for visualisation of the clustering effect. The map size is 2 $\text{deg}^2$ , the same as the size of the SIDES simulation.

We first generate, for different clustering strengths, the positions of all objects in the SIDES catalogue. There are several approaches that could be used to generate a density field that follows a designated 1D power spectrum. Here we assume a log-normal density field realisation, which is a good approximation for real cosmological density variations and a useful analytical form for simulation algorithms (e.g. Coles & Jones Reference Coles and Jones1991; Beutler et al. Reference Beutler2011). For a given $\theta_0$ , $w(\theta)$ can be converted to the angular power spectrum through the intrinsic Fourier relation (see Section 46 of Peebles Reference Peebles1980). For a sky area of 2 $\, \mathrm{deg}^2$ of interest here, we can adopt a flat-sky assumption and simplify the angular power spectrum to a 2D spatial Fourier transform of the map, following the convention used for CIB power spectra (e.g. Viero et al. Reference Viero2013). For $w(\theta)=\left(\theta / \theta_0\right)^{-\alpha} = A\theta^{-\alpha}$ , where A is the amplitude. In this relation we define $\theta$ to be in units of deg and A in deg $^{\alpha}$ . The resulting clustering power spectrum is the Fourier transform of $w(\theta)$ ,

(6) \begin{equation} \mathcal{P}(k_\theta) = \frac{2\pi A }{k_\theta^{2-\alpha}}\frac{2^{1-\alpha} \Gamma((2-\alpha) / 2)}{\Gamma(\alpha / 2)},\end{equation}

where $k_\theta$ labels the Fourier modes in units of radians $^{-1}$ . $ \mathcal{P}_\delta(k_\theta)$ is then in units of sr $^{-1}$ . We sample k-space with 1 000 modes and transform $ \mathcal{P}(k_\theta)$ to $ \mathcal{P}^\prime(k_\theta)$ using the relation of its covariance $w^\prime(\theta)=\ln(1+w(\theta))$ , so that this lognormal realisation of $ \mathcal{ P}^{\prime}(k_{\theta})$ is the same as the Gaussian realisation of $ \mathcal{P}(k_\theta)$ (see Coles & Jones Reference Coles and Jones1991). We generate the lognormal density field for a $1\,200 \times 1\,200$ grid on a $2\, \mathrm{deg}^2$ map through the inverse Fourier transform of $ \mathcal{P}^\prime(k_\theta)$ . We then draw Poisson samples from the lognormal density field for each cell plus a small random variation of the cell separation.

Next, we match flux densities of galaxies to the newly-generated positions by shuffling the flux densities to assure a uniform clustering of galaxies without any luminosity dependence. We assign galaxies with the new positions to a map with a pixel width of $12~\mathrm{arcsec}$ (as in the actual Herschel-SPIRE maps at 500 $\,\unicode{x03BC}$ m). We convolve the maps with the beam and add Gaussian noise with $\sigma=1.03\ \mathrm{mJy}~\mathrm{beam^{-1}}$ , chosen for the same reason as mentioned in Section 3.1. In Figure 3, we show examples of the generated clustered maps for five different clustering strengths to illustrate the effectiveness of the clustering algorithm. We adopt the adaptive binning set-up as in Section 3.4. We generate 100 clustered maps of each clustering strength by shuffling the flux densities and matching them with the generated positions. We find the histograms of 100 clustered maps of each clustering strength and stack them together to find the average histogram for each clustering strength. For comparison, we run 1 000 simulations, randomly assigning positions to galaxies in the SIDES catalogue, to produce maps in the same way, and obtain their stacked histograms and levels of scatter in each bin.

The map histograms are then transformed into consistent PDFs following

(7) \begin{equation}p_i = \frac{n_i}{N_{\mathrm{pix}}\Delta x_i},\end{equation}

where $p_i$ is the PDF in the ith bin, $n_i$ is the number of pixels in the ith bin, $N_{\mathrm{pix}}$ is the total number of pixels in the map, and $\Delta x_i$ is the bin width.

Figure 4.

Mean PDFs of 100 generated maps for each clustering strength or characteristic clustering scale $\theta_0$ (Equation 5) in coloured dashed lines. The PDF is computed from the average histogram by Equation (7). The red dashed line and the shaded region represent the mean PDF and standard deviation of 1 000 randomised maps. As the clustering increases, the peaks of the histograms shift to fainter flux densities and the tails grow.

We show the resulting mean PDFs of different clustering strengths in the left of Figure 4. The peaks of the histograms move to fainter flux densities as the clustering strength grows larger.Footnote ^c We perform a simple spline fit in the peak regions and align the peaks of these PDFs. In Figure 5, we notice a clear trend that as the clustering becomes stronger, that is, the characteristic clustering scale $\theta_0$ increases, the peak of the histogram is suppressed and the tails are boosted. In other words, as galaxies get more clustered, there is an increasing number of brighter pixels (sources merging into pixels) and fainter pixels (because underdense regions are also created by clustering).

Figure 5.

Top Panel: Same colour style as in Figure 4, but with the peaks of the histograms now aligned. Bottom Panel: The difference between the mean PDFs of the clustered maps and the randomised map, $\bar{p}_{\text{rand}} - \bar{p}_{\theta}$ , for each clustering strength. The red dashed line and the shaded region are the zero-residual and scatter from 1 000 randomised maps, respectively.

The effect seen in Figure 4 looks like a smoothly varying set of changes. We will find it convenient to first shift the means (as described in the next section), which is shown in the upper panel of Figure 5. In the bottom panel of Figure 5, we show the difference, $\bar{p}_{\text{rand}} - \bar{p}_{\theta}$ , between the mean PDFs of randomised maps, $\bar{p}_{\text{rand}}$ and those of clustered maps, $\bar{p}_{\theta}$ for each given $\theta_0$ .

We notice that the PDFs change in a systematic way, giving higher amplitude in the tails and a lower peak, representing an approximately single-parameter family of functions (and with the negative tail changing less significantly than the positive tail). Our goal is to create an asymmetric ‘W-shaped’ correction function, whose amplitude we expect to depend on the clustering strength, or $\theta_0$ . Then measuring the amplitude of the clustering should enable us to estimate how to correct the clustered histogram to provide the unclustered one, from which we can determine unbiased counts.

4.2. Correction function fit on the stacked histogram

We parameterise the correction function F(x) by first aligning the peaks of the histograms found by simple spline fits and fitting $\bar{p}_{\text{rand}} - \bar{p}_{\theta}$ . We do this using a function combining a second-order polynomial and skewed Gaussian on each side of the peak,

(8) \begin{equation}F(x)=\begin{cases}A_1\sigma_2 Q(x,\unicode{x03BC}_1,\sigma_1)\, \\\qquad \times \mathrm{skew}\!\left(x,\sigma_2,\unicode{x03BC}_2,a_1\right),& x\lt x_\mathrm{peak},\\[3pt]A_2\sigma_4 Q(x,\unicode{x03BC}_3,\sigma_3)\, \\\qquad \times \mathrm{skew}\!\left(x,\sigma_4,\unicode{x03BC}_4,a_2\right)+k,& x\gt x_\mathrm{peak}.\end{cases}\end{equation}

Here $Q(x, \unicode{x03BC},\sigma) = 1-(\frac{x-\unicode{x03BC}}{\sigma_{}})^{2}$ and $\mathrm{skew}\left(x, \sigma, \unicode{x03BC}, a\right)$ represents a skewed normal distribution with expectation value $\unicode{x03BC}$ , variance $\sigma^2$ , and skewness a.

After aligning the histogram peaks to the average histogram of randomised maps, we then adopt the Levenberg-Marquardt least-squares method (Section 15.5 in Press et al. Reference Press, Teukolsky, Vetterling and Flannery2007) to fit $\bar{p}_{\text{rand}} - \bar{p}_{\theta}$ with Equation (8), using equally-weighted bins to better recover the shape of the correction function. We explore the set of free parameters introduced in Equation (8), and find that some of them can be fixed to their mean values, i.e. independent of $\theta_0$ , while the other parameters can be fit by simple linear or log-linear relations with the characteristic clustering scale $\theta_0$ . The best-fit relations for all parameters are shown in Table 1. We note that the amplitudes of $\bar{p}_{\text{rand}} - \bar{p}_{\theta}$ vary with the beam size and noise level of the map, so the parameters in the correction function are treated case by case at this stage. In this section, we have adopted the SPIRE 500- $\unicode{x03BC}$ m beam and the GOODS-N noise level in this band, so the set of parameters in Table 1 is readily applicable for this particular field and band and is used later in Section 5.

Table 1.

Parameters in the correction function (Equation 8), either fixed or fit from a simple relation to the characteristic clustering scale $\theta_0$ in units of degree (Equation 5). We note that since we simulated SPIRE 500- $\unicode{x03BC}$ m maps with noise levels matching existing GOODS-N observations, the parameters in this table are applicable to this particular data set (see Section 5).

Table 1. Long description

A table with 11 rows and 3 columns. The columns are labeled Parameters, Fixed value, and Fit relation. The table presents parameters in the correction function, either fixed or fit from a simple relation to the characteristic clustering scale. Row 1: Parameters, A1; Fixed value; Fit relation, log A1 = (-0.089 ± 0.02) log² θ0 + (-0.02 ± 0.20) log θ0 - (5.21 ± 0.44). Row 2: Parameters, μ1; Fixed value, 3.48; Fit relation. Row 3: Parameters, σ1; Fixed value, 0.80; Fit relation. Row 4: Parameters, A2; Fixed value; Fit relation, log A2 = (0.78 ± 0.02) log θ0 + (-3.38 ± 0.07). Row 5: Parameters, μ2; Fixed value; Fit relation, log μ2 = (0.78 ± 0.03) log θ0 + (0.32 ± 0.10). Row 6: Parameters, σ2; Fixed value, 0.74; Fit relation. Row 7: Parameters, σ3; Fixed value; Fit relation, σ3 = (553.26 ± 22.59)θ0 + (1.12 ± 0.01). Row 8: Parameters, μ3; Fixed value; Fit relation, μ3 = (3247.44 ± 156.37)θ0 + (1.96 ± 0.04). Row 9: Parameters, α; Fixed value, 1.03; Fit relation. Row 10: Parameters, k; Fixed value; Fit relation, b = (-3560.90 ± 87.87)θ0 + (-0.01 ± 0.01).

In Figure 6 we show the correction function overlaid with $\bar{p}_{\text{rand}} - \bar{p}_{\theta}$ after fixing and fitting the simplified relation to $\theta_0$ . The correction function follows the shape of the residuals well and we show later in Section 4.4 that it sufficiently characterises the key features in the change of the histogram shapes, so that it almost completely debiases the P(D) galaxy counts from the effects of galaxy clustering.

Figure 6.

Examples of $\bar{p}_{\text{rand}} - \bar{p}_{\theta}$ at different clustering strengths or $\theta_0$ (dots), with their respective fits from Equation (8) in solid lines. The vertical dashed line represents the peak $x_{\text{peak}}$ in the equation. The shaded region represents the scatter in $\bar{p}_{\text{rand}} - \bar{p}_{\theta} $ in 100 simulated map pairs.

4.3. Two-point correlation function and power spectra

To determine the level of correction needed in the clustered map, we need to obtain the characteristic clustering scale, $\theta_0$ . The pointwise 2-point correlation function $w(\theta)$ is useful when there is a catalogue of galaxies in the observed field; however, the sources in a catalogue are usually flux-limited. For our application, we care about the clustering of relatively faint sources, since these are what bias the P(D) prediction as in Section 3.4. This means that the map-based power spectrum is a better place to focus our attention.

As mentioned in Section 4.1, the clustering power spectrum, $ \mathcal{P} (k)$ , is the Fourier transform of $w(\theta)$ , which describes the angular clustering of the source distribution. The angular power spectrum of an observed or simulated flux-intensity map is instead the 2D Fourier transform of the map intensity fluctuations in azimuthally-averaged Fourier component bins k:

(9) \begin{align} P(k)=\big\langle|\delta I(\vec{k})|^{2}\big\rangle,\end{align}

where $\delta I(\vec{k})$ is the intensity fluctuation field in Fourier space. Under the assumption that source fluxes are uncorrelated with source positions, the map power spectrum can be written as

(10) \begin{equation}P(k)=P_\mathrm{shot} + \bar I^2 \mathcal{P}(k),\end{equation}

where $\bar I$ is the mean map intensity, $\bar I = \Sigma_i{S_i}/\Omega$ for a source catalogue with galaxy fluxes $S_i$ over a total angular area $\Omega$ and $P_\mathrm{shot}$ is the contribution from Poisson shot noise.

We use the simulated clustered maps, without convolving with the beam or adding noise, so that we can calculate power spectra directly of the maps without any corrections for masking, beam functions, transfer functions, etc. We calculate the power spectra by taking the 2D Fourier transform at a sampling rate of the pixel width. We divide k-space into $2^8$ bins, sum over the power in each k bin, and normalise the power spectrum by the number of pixels in each k bin.

We calculate the power spectra for 100 maps of each clustering strength and stack the power spectra in each k-bin. We then remove the constant shot-noise component ( $P_\mathrm{shot}$ Footnote ^d ) from the power spectra by fitting a constant in the range  $1.93\, \text{arcmin}^{-1} \leq k_\theta \lt 3.00\, \text{arcmin}^{-1}$ , where $1.93 \,$ arcmin $^{-1}$ corresponds to the scale of the beam FWHM, and $3 \,$ arcmin $^{-1}$ is set by the map size. The resulting power spectra can be written as the power-law relation

(11) \begin{equation}\log \left(P(k_\theta) \right)=-1.2\log(k_\theta)+C,\end{equation}

where $P(k_\theta)$ is the amplitude of the power spectrum in units of Jy $^2$ sr $^{-1}$ and $k_\theta$ the Fourier modes in units of $\text{arcmin}^{-1}$ . The slope $-1.2$ is fixed from the slope of the 2-point correlation function $w(\theta)$ by Equation (6) and C is determined by the relation

(12) \begin{equation}\begin{aligned}C &= 0.8 \log(\theta_0)+ C_0, \\\end{aligned}\end{equation}

where $\theta_0$ is in units of deg. The slope 0.8 is similarly set by the power-law index of $w(\theta)$ . The intercept $C_0$ is set by the catalogue and can be analytically calculated by substituting $A\,{=}\,\theta_0^{\alpha}$ into Equation (6) and noting that $P_\mathrm{shot}$ in P(k) has already been subtracted from Equation (11). We find $C_0\,{=}\,$ 6.02, 5.82 and 5.39 at 250, 350, and 500 $\unicode{x03BC}$ m, respectively. We show the example power spectra and the fits at 500 $\unicode{x03BC}$ m in Figure 7.

Figure 7.

Examples of power spectra for generated maps at 500 $\unicode{x03BC}$ m with different clustering strengths or characteristic clustering scales $\theta_0$ (Equation 5). The constant shot noise term has been subtracted from these power spectra. The solid line is given by Equation (11).

Figure 8.

Flow chart illustrating the process for correcting the clustering bias in galaxy number counts fitting through P(D) fluctuation analysis. Since P(D) fits the flux histogram of the map, the process focuses on removing the impact of clustering (2-point statistics) on the flux histogram (1-point statistics). The specific form used in the correction is derived from simulations and the amplitude of clustering from the data themselves.

By inverting Equation (12), we can estimate the characteristic clustering scale $\theta_0$ from the map power spectrum and find the appropriate form of the correction function (Equation 8) from the parameters in Table 1. We note that in this section we focus on testing the correction at 500 $\unicode{x03BC}$ m, while the different clustering relations for 250 and 350 $\unicode{x03BC}$ m shown in Equation (12) are quoted later in Section 5. We now have the 2-point information on spatial clustering, and we can propagate its impact on the flux histograms (i.e. 1-point statistics) via the correction function, and then recover the unbiased galaxy number counts using the P(D) fluctuation analysis. In this way, we simultaneously obtain the 1-point (flux histograms or PDFs) and 2-point statistics (power spectrum or 2PCF) contained in a map, which together comprise nearly all the information available in an unresolved galaxy imaging survey.

4.4. Results for corrected galaxy number counts

We summarise the correction process in a flow chart (Figure 8). The usual P(D) framework is built exclusively on fitting the flux histogram, while the impact of clustering is not considered, leading to biases in the estimated galaxy number counts. However, we can determine the clustering strength from the data themselves and then correct the flux density histogram to find what it would be for unclustered sources. P(D) analysis conducted on this corrected histogram should therefore result in unbiased galaxy number counts.

To test our correction method, we use a simulated clustered map with a similar clustering strength as in the SIDES simulation, $\theta_0 = 0.33^{\prime\prime} \pm 0.003^{\prime\prime}$ calculated from Equation (11) after fitting the SIDES simulated map power spectrum by Equation (11). We follow the correction processes as summarised in Figure 8. We estimate the power spectrum using Equation (9) and fit it with Equation (11) after subtracting $P_{\text{shot}}$ , which is fit as a constant in the flat (high k) part of P(k). The uncertainty in the power spectrum is the dispersion of the power in each k bin. We then find the characteristic clustering scale $\theta_0$ by inverting Equation (12). We insert $\theta_{0, \text{fit}}$ into the relations in Table 1, where we find all the parameters needed to calculate the correction function in Equation (8). We add the correction function to the PDF of the map histogram, with the binning created adaptively as in Section 4.1. The corrected histogram is obtained by converting from the corrected PDF following Equation (7).

For the uncertainty in the characteristic clustering scale, $\theta_0$ , we draw samples from a normal distribution within its $1\sigma$ error, estimated from Equation (12) after fitting the power spectrum by Equation (11). We then find the error bar of the correction function by drawing 20 000 samples of each variable coefficient in Table 1 for each $\theta_0$ draw within its covariance range, and taking the 16th and 84th percentiles of the resulting correction function. To ensure numerical stability in the positive tail, we apply a mild winsorisation to the ensemble of correction function realisations, clipping extreme outliers at the 1st and 99th percentiles before estimating the uncertainty envelope. We assume the error bar to be symmetric, so $\sigma_{\text{corr, sample}} = (F_{\text{corr, sample}, 84} - F_{\text{corr, sample}, 16})/2$ . We then add in quadrature a 10% uncertainty ( $\sigma_{\text{sys}} = 10\%~\sigma_{\text{corr}}$ ) to account for systemetics such as small beam variations (see Section 6.3). The total uncertainty in the correction function is then $\sigma_{\text{corr}} = (\sigma_{\text{corr, sample}}^2 + \sigma_{\text{sys}}^2)^{1/2}$ . This is included in the likelihood function (Equation 4), where $\sigma_i$ becomes

(13) \begin{equation}\sigma_i =\sqrt{\sigma_{\text{Poisson}}^2+\sigma_{\text{corr}}^2} = \sqrt{n_i + \sigma_{\text{corr}}^2},\end{equation}

while for the standard process without using a correction function, included later for comparison, we use $\sigma_i = \sigma_{\text{Poisson}}$ .

Following this procedure, we show in Figure 9 the flux histograms before and after correction, and the corresponding final galaxy counts fit in the same fashion as in Section 3.4. The corrected counts clearly recover the true underlying counts much better than the fit to the uncorrected counts.

Figure 9.

Left panel: Observed histograms of the clustered map before (blue) and after (red) correction in solid lines, after aligning their peaks. The grey dashed lines are the flux cutoffs beyond which the data are not used in the fit. Right panel: P(D) galaxy number counts models before (blue) and after (red) correction. The galaxy number counts are Euclidean normalised here. The first node is an upper limit. The shaded band is the 68% confidence region for each fit. The black dots are the real galaxy counts from the SIDES catalogue with Poisson errors in an arbitrary binning for illustration.

To evaluate the goodness of fit of the galaxy counts in comparison to the real catalogue, we first bin the galaxies in the catalogue according to their brightness, with logarithmic bins at the brighter end to ensure there are more than 10 galaxies in each bin. The bins are extra fine here to fully trace the shape of the underlying counts, especially at the faint end. We then transform the number of galaxies in each bin to the galaxy number counts by dividing it by the bin size and the total area of the map. The uncertainty in the resulting galaxy counts from the catalogue is evaluated by taking the errors in each bin to be $\sqrt{n_i}$ , assuming a Poisson distribution of the number of galaxies $n_i$ in the ith bin, and propagating the error into the number counts, $\sigma_{\text{cat}}$ . The binned counts then represent the true galaxy counts in SIDES.

We can evaluate the goodness of fit of the galaxy counts by calculating $\chi^2$ . The fit nodes of the galaxy counts are interpolated by power laws onto the central bin values. We then calculate

(14) \begin{equation} \chi^{2}=\sum_{i} \frac{\left(D_{i}-M_{i}\right)^{2}}{\sigma_{i}^{2}},\end{equation}

where $D_i$ is the galaxy count in the ith bin calculated from the simulation catalogue, $\sigma_i$ is the error on $D_i$ , $M_i$ is the interpolated counts in the ith bin from the fit results, and there are 805 degrees of freedom in the galaxy counts. We find that $\chi^2$ drops from 10 581 to 915 after the correction, a dramatic improvement. This demonstrates that the uncorrected fit is systematically biased, and the proposed clustering correction dramatically improves the P(D) galaxy number counts fit. We have therefore shown that the proposed correction method is effective on simulated Herschel-SPIRE 500- $\unicode{x03BC}$ m maps, with noise levels typical of GOODS-N deep observations.

5.1. GOODS-N field

5.1.1. Data

The confusion noise for Herschel-SPIRE is 5.8, 6.3, and 6.8 mJy beam $^{-1}$ at 250, 350 and 500 $\,\unicode{x03BC}$ m (Nguyen et al. Reference Nguyen2010) and the GOODS-N maps are well into the confusion-limited regime. Pixel sizes are 6, 8.33 and 12 arcsec at 250, 350, and 500 $\unicode{x03BC}$ m, respectively. The original observations span around $0.64\,$ deg $^{2}$ , but we cropped the map to retain the central area of approximately $0.3\,$ deg $^2$ where the variation of the noise level is within 10%. The resulting noise levels of the cropped maps are shown in Table 2. The very deep SPIRE data in the GOODS-N field makes it ideal for studying the faint galaxy number counts. Although larger fields would effectively reduce field-to-field variance, their noise behaviour would be more complicated, requiring multi-zone P(D) analyses (e.g. Vernstrom et al. Reference Vernstrom2014), and the shallower depths of the other fields adds little constraining power at the faintest flux densities.

Table 2.

Specifications of the Herschel-SPIRE GOODS-N maps (Shirley et al. Reference Shirley2021), including the beam size FWHM, mean instrumental noise $\sigma_{\text{inst}}$ , mean instrumental noise of the cropped map used for this analysis $\sigma_{\text{inst, crop}}$ , and the confusion noise $\sigma_{\text{conf}}$ .

Table 2. Long description

A table comparing specifications of Herschel-SPIRE GOODS-N maps at different wavelengths. The table has four columns and three rows. The columns are labeled as SPIRE band (μm), FWHM (arcsec), σ_inst (mJy), σ_inst, crop (mJy), and σ_conf (mJy). The rows are labeled with the following values: Row 1: SPIRE band, 250; FWHM, 18.15; σ_inst, 1.85; σ_inst, crop, 0.88; σ_conf, 5.8. Row 2: SPIRE band, 350; FWHM, 25.15; σ_inst, 1.81; σ_inst, crop, 0.87; σ_conf, 6.3. Row 3: SPIRE band, 500; FWHM, 36.3; σ_inst, 2.15; σ_inst, crop, 1.03; σ_conf, 6.8.

5.1.2. Simulations and correction function

Similarly to Section 4, we simulate 1 000 randomised maps and 100 clustered maps of clustering strengths from 2.06 arcsec to 63.61 arcsec, excluding unrealistically strong clustering. The generated maps are 2 deg $^2$ , larger than GOODS-N to reduce the error after averaging the histogram. We then fit $\bar{p}_{\text{rand}} - \bar{p}_{\theta}$ with the correction function form in Equation (8) in a similar manner as in Section 4.2, and show the parameters of the correction function for 250 and 350 $\unicode{x03BC}$ m in Tables 3 and 4, respectively.

Table 3.

Same parameters as in Table 1 for the GOODS-N map at 250 $\unicode{x03BC}$ m (Equation 8).

Table 3. Long description

A table with 11 rows and 3 columns. The columns are labeled Parameters, Fixed value, and Fit relation. The table presents various parameters along with their fixed values and fit relations. The parameters include A1, μ1, σ1, μ2, σ2, σ3, A2, μ3, σ3, μ4, σ4, and k. Each parameter has a corresponding fixed value and a fit relation, which are detailed in the table. The fit relations are given in logarithmic forms with specific equations and values. The table provides a comprehensive overview of the parameters and their relationships for the GOODS-N map at 250 m.

Table 4.

Same parameters as in Table 1 for the GOODS-N map at 350 $\unicode{x03BC}$ m (Equation 8).

Table 4. Long description

A table with two columns: Parameters and Fit relation. The Parameters column lists various parameters with some having fixed values and others having fit relations. The Fit relation column provides mathematical expressions for these parameters. The table has 14 rows, each corresponding to a different parameter or fixed value. For example, the first row shows the parameter A1 with a fit relation involving log transformations and other parameters. The second row lists the fixed value for parameter mu1. The table continues with similar entries for other parameters and their respective fit relations or fixed values.

We take the auto-power spectrum at 250, 350 and 500 $\unicode{x03BC}$ m from Viero et al. (Reference Viero2013), choosing the deepest flux cut at $50~\mathrm{mJy}$ for the power spectrum. We then fit Equation (11), finding the intercepts $C= 2.64$ , $2.54$ , and $2.18$ , and the corresponding $\theta_\mathrm{0} = 0.26^{\prime\prime} \pm 0.002^{\prime\prime}$ , $0.37^{\prime\prime} \pm 0.004^{\prime\prime}$ , and $0.43^{\prime\prime} \pm 0.004^{\prime\prime}$ for 250, 350, and 500 $\unicode{x03BC}$ m, respectively, through Equation (12). We note that these $\theta_0$ values are derived from an empirical relationship, assuming $w(\theta)$ follows a single power law.

We then insert $\theta_0$ values into the relations in Tables 1, 3, and 4 to determine the parameters in the correction functions (Equation 8) at 250, 350, and 500 $\unicode{x03BC}$ m. Adding the corresponding correction functions to the original map PDFs and converting to flux histograms, we obtain the clustering-corrected histograms to feed into the P(D) framework. In light of the beam sizes and noise levels at 250 $\unicode{x03BC}$ m and at 350 $\unicode{x03BC}$ m (different than at 500 $\unicode{x03BC}$ m), the adaptive binning and flux cutoffs are carefully chosen following Section 3.4. The error bars in the correction functions are estimated and propagated into the likelihood in the same way as in Section 3.4.

5.1.3. Models and fits

We choose the node-based model for each band, (see Section 3). Each band has been tested on a simulated map to ensure the efficacy of the model. The brightest node is fixed the simulation when the maximum flux in the map is relatively small, thus the bright galaxy number counts are not well constrained by the map (see Section 3.2).

We perform the fitting on the corrected GOODS-N histograms using the same emcee configuration as in Section 3.4, and show the results of galaxy number count fits before and after correction in Figure 10 and in Table 5. At 10 mJy, we find the counts are corrected by a factor of $1.12 \pm 0.14$ , $1.25\pm 0.38$ and $1.63\pm 0.27$ at 250, 350, and 500 $\unicode{x03BC}$ m. We notice that the clustering bias is indeed the most significant at 500 $\unicode{x03BC}$ m, where the beam is the largest.

Figure 10.

Galaxy number count fits before (blue) and after (red) correcting the histograms at 250, 350 and 500 $\unicode{x03BC}$ m (top to bottom). Colours and symbols follow Figure 9.

6.1. Comparison to previous studies

We show the corrected galaxy number counts from this work compared with other measurements at the same wavelength in Figure 11. For all three wavebands, we include measurements from Glenn et al. (Reference Glenn2010), Béthermin et al. (Reference Béthermin2012), and Wang et al. (Reference Wang2019). At 500 $\unicode{x03BC}$ m, we also compare with deep 450 $\unicode{x03BC}$ m measurements using the Submillimetre Common-User Bolometer Array-2 (SCUBA-2; Holland et al. Reference Holland2013) camera on the James Clerk Maxwell Telescope (JCMT). We do not convert 450- $\unicode{x03BC}$ m results to 500 $\unicode{x03BC}$ m to avoid assuming a specific shape of the spectral energy distribution for all galaxies.

Table 5.

Corrected galaxy number count fits on GOODS-N maps at 250, 350, and 500 $\unicode{x03BC}$ m. The lower and upper limits are from the 68% confidence region (1 $\sigma$ ) of the MCMC fits, except for the first node which is an upper limit due to the assumption of the counts being zero outside the flux range of the fit nodes.

Table 5. Long description

A table comparing galaxy number counts at different wavelengths. The table has three main sections, each corresponding to a different wavelength: 250 micrometers, 350 micrometers, and 500 micrometers. Each section contains three columns: Node (in mJy), log10 of differential number counts (dN/dS) in units of (sr^-1 Jy^-1), and the lower and upper limits from the 68% confidence region of the MCMC fits. The table lists values for different nodes at each wavelength. For 250 micrometers, the nodes range from 0.01 to 0.21 mJy, with corresponding log10 (dN/dS) values and confidence limits. For 350 micrometers, the nodes range from 0.07 to 100.00 mJy, with corresponding log10 (dN/dS) values and confidence limits. For 500 micrometers, the nodes range from 0.11 to 100.00 mJy, with corresponding log10 (dN/dS) values and confidence limits. The table provides detailed data on the galaxy number counts and their uncertainties at different wavelengths.

Figure 11.

Corrected galaxy number counts at 250, 350 and 500 $\unicode{x03BC}$ m from top to bottom. The red line and the shaded region are labelled the same as in Figure 10. The first node in each panel is an upper limit. Data points are previous measurements at 250, 350, 500 $\unicode{x03BC}$ m (Glenn et al. Reference Glenn2010; Béthermin et al. Reference Béthermin2012; Wang et al. Reference Wang2019). For 500 $\unicode{x03BC}$ m measurements, deep 450 $\unicode{x03BC}$ m counts are also included from Chen et al. (Reference Chen2013), Wang et al. (Reference Wang2017), Hsu et al. (Reference Hsu, Cowie, Chen and Barger2024), and Gao et al. (Reference Gao2024).

Glenn et al. (Reference Glenn2010) implemented a P(D) fluctuation analysis on three fields of Herschel-SPIRE SDP data. GOODS-N is one of the three fields and the deepest one, dominating the measurement below 10 mJy. We find general agreement with Glenn et al. (Reference Glenn2010) in all three bands, the later with no clustering correction included. This is likely due to the noisy data in SDP that contributes to large error bars, as seen in Figure 10. Still, the trend of the correction can be seen from the median of the data points, especially at 500 $\unicode{x03BC}$ m, where Glenn et al. (Reference Glenn2010) lies slightly above our determination at $S\simeq10$ mJy, but falls below at $S\lesssim10$ mJy, in agreement with the trend of clustering bias shown in Figure 10. This indicates that our clustering correction is effective, while the depth of the full SPIRE data set probes more robustly down to a fainter limit.

Béthermin et al. (Reference Béthermin2012) used full-depth Herschel-SPIRE data and stacked on deep 24 $\unicode{x03BC}$ m catalogues in redshift slices. Empirical methods of clustering correction were applied to the stacked flux, based primarily on the clustering of sources in the 24 $\unicode{x03BC}$ m catalogue. We find our result in slight tension with the stacking result from Béthermin et al. (Reference Béthermin2012), with their counts lying slightly below ours at $S \lesssim 10$ mJy and above ours at $S \gtrsim 10$ mJy. This again mirrors the trend of clustering-biased counts in simulations through P(D) analysis, as shown in Figure 2. Although the methods differ in nature and careful examination of their implementations would be needed, we suspect the empirical clustering correction in Béthermin et al. (Reference Béthermin2012) based on the 24 $\unicode{x03BC}$ m catalogue is less accurate for the clustering in the SPIRE bands, which already varies across the three SPIRE bands, as we found in Section 5.1.2.

At all three SPIRE wavebands, our results are inconsistent with Wang et al. (Reference Wang2019) at both the fainter and brighter ends. At $S \lesssim 10$ mJy, the counts from Wang et al. (Reference Wang2019) show an excess compared to our work, while a deficit is seen at $S \gtrsim 10$ mJy. A similar discrepancy is noted in the 500- $\unicode{x03BC}$ m comparison with 450- $\unicode{x03BC}$ m results. Wang et al. (Reference Wang2019) performed a multiwavelength deblending based on compiled shorter-wavelength catalogue of the COSMOS field. Because this deblending method inherently relies on positional priors, its accuracy is heavily dependent on the completeness of the shorter-wavelength catalogue, which is known to miss highly obscured submillimetre sources. While a detailed analysis of their algorithm is beyond the scope of this paper, we suggest that this prior-dependency drives the discrepancy. Catalogue incompleteness could result in the simultaneous artificial deficit of bright sources and excess of faint sources seen in the figure. This could be improved with the James Webb Space Telescope (JWST)’s excellent recovery of submillimetre galaxies (e.g. Hill et al. Reference Hill2024; Liu et al. Reference Liu2026).

The 450- $\unicode{x03BC}$ m results from SCUBA-2 in general show excellent consistency with our 500- $\unicode{x03BC}$ m results. Wang et al. (Reference Wang2017) presented the preliminary results of the SCUBA-2 Ultra Deep Imaging EAO Survey (STUDIES) compared to Gao et al. (Reference Gao2024), which estimated direct counts in the full survey depth. Because of this, we include the full counts from Gao et al. (Reference Gao2024) and only the non-direct counts from Wang et al. (Reference Wang2017), the latter derived through P(D) fluctuation analysis with a single power-law model without clustering correction. We find that the direct counts from Gao et al. (Reference Gao2024) agree well with our results down to ${\simeq}$ $2\,$ mJy, while the single power law P(D) results from Wang et al. (Reference Wang2019) are largely consistent with our results at ${\simeq}$ $1\,$ mJy. Chen et al. (Reference Chen2013) and Hsu et al. (Reference Hsu, Cowie, Chen and Barger2024) measured combined counts on multiple fields with SCUBA-2, Chen et al. (Reference Chen2013) including two lensing-cluster fields and one blank field, while Hsu et al. (Reference Hsu, Cowie, Chen and Barger2024) extended this to four lensing-cluster fields and the same blank field. In general, we find excellent agreement with (Hsu et al. Reference Hsu, Cowie, Chen and Barger2024) down to ${\simeq}$ $0.1\,$ mJy, and we are consistent with (Chen et al. Reference Chen2013) except for their faintest counts around ${\simeq}$ $1.3\,$ mJy. Overall, we find great consistency with the direct, non-direct, and multi-field counts through deep SCUBA-2 observations at 450 $\unicode{x03BC}$ m, demonstrating the constraining power of P(D) with the clustering correction proposed in this work.

Finally, we note that at 250 $\unicode{x03BC}$ m where the beam is the smallest, our results at the faint end show a trend towards flattening of the Euclidean-normalised counts, which is expected as the counts approach the faint-end turnover. We anticipate that applying this method on deeper observations will fully resolve the extragalactic background light (EBL), and hence provide tighter constraints on models of galaxy formation.

6.2. Comparison to simulations

In Figure 11, we also show the counts from SIDES and GALFORM (Cowley et al. Reference Cowley2019). The SIDES counts show excellent consistency with our results for all three wavebands, despite a slight overestimation around $2\,$ mJy at 250 and 350 $\unicode{x03BC}$ m and around $4\,$ mJy at 500 $\unicode{x03BC}$ m. Semi-analytic models (SAMs) such as GALFORM historically struggle to reproduce observed populations of submillimetre galaxies. This tension is largely attributed to the complex treatment of baryonic feedback and dust properties within the simulations (e.g. Lovell et al. Reference Lovell, Geach, Davé, Narayanan and Li2021). We find that the GALFORM results reproduce the observed counts well at 250 $\unicode{x03BC}$ m, while an excess is found at 1 mJy $\lesssim S\lesssim10$ mJy for 350 and 500 $\unicode{x03BC}$ m. Resolving this discrepancy will require rigorous scrutiny of the underlying physical prescriptions. Consequently, our results emphasise that submillimetre number counts remain a stringent and highly valuable diagnostic for testing theories of galaxy formation.

Encouragingly, more recent SAMs and hydrodynamic simulations have made progress in reconciling simulated populations of dusty star-forming galaxies with certain observational constraints (e.g. Lovell et al. Reference Lovell, Geach, Davé, Narayanan and Li2021; Araya-Araya et al. Reference Araya-Araya2025). A direct comparison of our clustering-corrected Herschel-SPIRE counts against predictions from these updated simulations will therefore provide interesting constraints on galaxy evolution models in future.

6.3. Systematic effects

The application of P(D) fluctuation analysis requires any variation of the beam to be small and that the noise is nearly Gaussian, which are generally true for Herschel-SPIRE blank-field observations. Glenn et al. (Reference Glenn2010) tested the effect of these assumptions on P(D) analysis and found that a constant Gaussian beam was a good approximation for a small field like GOODS-N. Furthermore, the beam profile in the final SPIRE data is highly stable and well-characterised (Griffin et al. Reference Griffin2010), minimising its systematic impact. Similarly, (Glenn et al. Reference Glenn2010) simulated noise variations in SPIRE observations and found it to be negligible for P(D) results.

In this work, we introduced a correction process using the SIDES simulations. As mentioned in Section 4.2, we incorporated an additional 10% systematic uncertainty to account for potential clustering effects not captured by the simulation, such as the assumption of a single power-law model for the angular correlation function, $w(\theta)$ .

Finally, the roughly 5% SPIRE calibration uncertainty (Bendo et al. Reference Bendo2013; Hopwood et al. Reference Hopwood2015) has been neglected in our P(D) analysis, since it acts as a global scaling of the flux density and does not meaningfully alter the shape of the pixel histogram. Additionally, the propagated effect of this calibration error on the final number counts is entirely subdominant compared to the existing uncertainties.

6.4. Field-to-field variance

Accounting for field-to-field variance is important for measuring galaxy counts due not only to the variation of the counts in different fields, but also the variation of the clustering. In this work, we address this by measuring the counts in the data themselves, but using an estimate of the clustering scale, $\theta_0$ , from averaging over a much larger field. We tested the sensitivity of the correction function by shifting $\theta_0$ by $\pm 2\,\sigma$ from its average value mapped from the power spectrum (see Section 5.1.2). We find the variations in the corresponding correction function are smaller than its intrinsic uncertainty (shown in Tables 1, 3, and 4).

6.5. Gravitational lensing

Gravitational lensing can visibly alter the observed pixel histogram by boosting faint sources into brighter flux bins. At the extremely bright end, strong lensing results in a significant overestimation of the bright galaxy counts (Lima et al. Reference Lima, Jain and Devlin2010; Paciga et al. Reference Paciga, Scott and Chapin2009; Negrello et al. Reference Negrello2017). However, this strong lensing regime mostly biases the brighter counts, where the number-count slope is steep (see e.g. Bartelmann & Schneider Reference Bartelmann and Schneider2001), dominating above 100 mJy at 350 and 500 $\unicode{x03BC}$ m (Béthermin et al. Reference Béthermin2017). This is higher than the brightest source in GOODS-N and significantly above the bright-end histogram cut-off ( ${\simeq}$ 10 mJy; see Sections 3.4 and 5.1.2) we adopted in our P(D) analysis of the simulated 2 deg $^2$ maps or the observed GOODS-N maps.

In the faint regime, weak lensing introduces a magnification bias that acts as an extra term for $w(\theta)$ (e.g. Hildebrandt et al. Reference Hildebrandt2013). The effect of this extra term mimics that of clustering by widening the P(D) distribution. However, the spatial clustering of galaxies dominates over the weak lensing signal in these submillimetre wavelengths by roughly two orders of magnitude (e.g. Wang et al. Reference Wang2011).

6.6. More complicated clustering functions

We have demonstrated that our correction method performs well for realistic clustering typical of blank fields ( $\theta_0$ up to around 1 arcsec. Within this range, our proposed functional form of the clustering correction (Equation 8) well characterises the clustering bias induced in the shape of the P(D) histogram (Figure 6). However, we find that the correction function breaks down when clustering becomes even stronger, suggesting that a more accurate correction form would be required for extremely clustered environments, such as galaxy clusters.

Throughout this work, we have adopted a single power-law functional form for the 2PCF, $w(\theta)$ , to describe galaxy clustering in blank fields. While this simple model is well-motivated, in practice the clustering strength will depend on galaxy properties, i.e., more massive (or brighter) galaxies tend to cluster more strongly. Furthermore, when extending to larger scales the simple power-law approximation will break down. This all suggests that modelling realistic clustering with multi-component power-law forms may yield further improvements to the correction.

6.7. Other future improvements

All of the issues discussed in the above subsections are areas where improvements could be made in future. In addition, the shot-noise part of the power spectrum contains some useful information, which we have been neglecting here. We confirmed that the level is in rough agreement with what is expected from our counts, but in future it should be possible to include the shot-noise level in a self-consistent analysis of the 1- and 2-point statistics. Another direction for future study is the direct constraints on galaxy evolution models, where one could make use of the joint P(D) distribution for the three SPIRE wavebands.