The Herschel-SPIRE Dark Field – II. A P(D) fluctuation analysis of the deepest Herschel image of the submillimetre universe

ABSTRACT

The Herschel-SPIRE Dark Field is the deepest field produced by the SPIRE instrument pushing down below the galaxy confusion limit in each of the 250, 350, 500 |$\mu$|m bands. Standard source extraction techniques inevitably fail because of this, and we must turn to statistical methods. Here, we present a P(D)  – probability of deflection – analysis of a 12|$^{\prime }$| diameter region of uniform coverage at the centre of the Herschel-SPIRE Dark Field. Comparing the distribution of pixel fluxes from our observations to the distributions predicted by current literature models, we find that none of the most recent models can accurately recreate our observations. Using a P(D) analysis, we produce a fitted differential source count spline with a bump in the source counts at faint flux densities, followed by a turnover at fainter fluxes, required to fit the observations. This indicates a possible missing component from the current literature models that could be interpreted perhaps as a new population of galaxies, or a missing aspect of galaxy evolution. Taking our best-fitting results, we also calculate the contribution to the cosmic infrared background (CIB) in each of the bands, which all agree with the Planck CIB measurements in this field.

galaxies: evolution, cosmology: observations, submillimetre: diffuse background, submillimetre: galaxies

1 INTRODUCTION

In 1967, a background of light emitted from young galaxies was predicted by Partridge & Peebles (1967a, b). Subsequent surveys in the optical to identify the galaxies responsible for this background were unsuccessful (Davis & Wilkinson 1974). Ultimately, it was discovered that interstellar dust played a key role in obscuring the background’s origin, absorbing and reprocessing the ultraviolet (UV) light emitted from young stellar objects to far-infrared (FIR) wavelengths (Kaufman 1976), meaning this background could lie between 30 and 500 |$\mu$|m (Stecker, Puget & Fazio 1977). Detection of this predicted cosmic infrared background (CIB) was ultimately made with the Far Infrared Absolute Spectrophotometer (FIRAS) instrument (Mather, Fixsen & Shafer 1993; Puget et al. 1996; Fixsen et al. 1998) aboard the Cosmic Background Explorer (COBE) satellite (Boggess et al. 1992), and later confirmed by the Diffuse Infrared Background Experiment (DIRBE) instrument (Schlegel, Finkbeiner & Davis 1998). Infrared observations are important for understanding galaxy evolution; the energy contained within the CIB is comparable with the energy contained within the cosmic optical background (COB) emission, i.e. approximately half of the emission from galaxies is at FIR wavelengths. As such, it follows that to fully understand the evolution of galaxies we should also understand their emission at FIR wavelengths, not just their emission in optical bands (Ryter & Puget 1977; Soifer, Neugebauer & Houck 1987).

Launched in 2009, the Herschel Space Observatory provided a revolutionary step forward in exploring the cosmic history of dust-obscured galaxies and large-scale structure formation (Pilbratt et al. 2010). The SPIRE instrument on-board Herschel, in particular, was able to make simultaneous observations at wavelengths of 250, 350, and 500 |$\mu$|m (Griffin et al. 2010). The Herschel mission carried out many legacy programmes over its lifespan, including the HerMES (Oliver et al. 2012), H-ATLAS (Clements et al. 2010; Eales et al. 2010), and HLS (Egami et al. 2010) surveys. Of particular importance to this study is the Herschel Lensing Survey (HLS). However, all deep Herschel surveys (except the HLS survey that used gravitational lensing of clusters to reach the deepest fields) suffered from confusion noise due to the relatively low resolution of the Herschel-SPIRE instrument at FIR wavelengths, and the high number density of FIR sources on the sky. During Herschel’s 4-yr mission, the Herschel-SPIRE instrument routinely observed a patch of dark sky for calibration of the instrument (Pearson et al. 2014). By combining 141 of these observations, we have assembled the deepest Herschel-SPIRE image ever produced (Pearson et al. 2025): roughly 10 times the integration time of the HLS survey.

Pearson et al. (2025) presented source count results from conventional source extraction (SUSSEXtractor) and list-driven cross-identification (XID) methods (Hurley et al. 2017) on the Herschel-SPIRE Dark Field, revealing a good agreement with literature counts at fluxes |$\gtrsim$|20 mJy. With conventional source extraction methods (SUSSEXtractor), our results at fainter fluxes (⁠|$\lesssim$|20 mJy) were significantly affected by confusion noise. We, therefore, applied a list-driven XID source extraction technique – using 24 |$\mu$|m catalogues from the Spitzer-MIPS instrument as a prior – to begin to probe around and below the confusion-limited flux level.

However, to probe to even fainter fluxes (⁠|$\lesssim$|1 mJy), we must turn to statistical methods, such as P(D), to exploit the unique depth of the Herschel-SPIRE Dark Field (Scheuer 1957). Previous work on Herschel data has already shown this technique to be a valuable method of overcoming the significant confusion noise seen in Herschel-SPIRE observations (Glenn et al. 2010).

In Section 2, we describe the data and observations of the Herschel-SPIRE Dark Field used in this study. In Section 3, we present the P(D) method and results from this analysis of the Herschel-SPIRE Dark Field. In Section 4, we use the results from our P(D) analysis to calculate the contribution of these sources to the CIB and compare the results to Planck HFI observations. In Section 5, we discuss the potential interpretations of these new P(D) differential source count results. We finally conclude by proposing potential extensions and implications of these results in Section 6.

2 OBSERVATIONS

In this study, we make use of observations of the Herschel-SPIRE Dark Field (Pearson et al. 2025). The Herschel-SPIRE Dark Field is composed of 141 stacked, mean-subtracted, calibration observations in scan-map mode of a roughly circular patch of sky, 30|$^{\prime }$| in diameter, near the North Ecliptic Pole at RA = 17h40m12s, Dec. = +69d00m00s. It is the deepest field observed by the SPIRE instrument. These observations reach below the source confusion limits in every SPIRE band: PSW (250 |$\mu$|m), PMW (350 |$\mu$|m), and PLW (500 |$\mu$|m). Because of this, the Herschel-SPIRE Dark Field is highly confused, which makes it ideally suited for a statistical P(D) approach once standard source extraction techniques are exhausted.

Within the Herschel-SPIRE Dark Field, we define a smaller region at the centre of the maps which we use for our analysis. We refer to this 12|$^{\prime }$| diameter, circular region as the Deep Region. This smaller field is selected as it has an almost uniform coverage and is void of any contamination from significant diffuse emission or bright sources. When using a statistical technique like P(D) analysis, it is important to select a uniform field so as not to affect the statistical interpretation of the pixel counts. Non-uniformities in the coverage will affect the interpretation of the pixel brightnesses.

The instrumental noise for the (250, 350, 500) |$\mu$|m SPIRE band observations has been measured to scale as (9.0, 7.5, 10.8) mJy per |$\sqrt{N_{\rm reps}}$|⁠, where |$N_{\rm reps}$| is the total number of map repeats in the observation (Griffin et al. 2010). The SPIRE Dark Field is assembled from 141 maps, each with 4 repeats (for a total of 564 repeated observations). Table 1 presents the calculated instrument noise levels for the observations used in our study. In our maps, we expect scaled instrumental noise levels of (0.379, 0.316, 0.455) mJy in the (250, 350, 500) |$\mu$|m bands, respectively.

Table 1.

Instrument noise (⁠|$\sigma _{\text{inst}}$|⁠), limiting sensitivity (⁠|$\sigma _{\text{limit}}$|⁠), and confusion noise (⁠|$\sigma _{\text{conf}}$|⁠), for each SPIRE band. The confusion noise here is reproduced from Nguyen et al. (2010). The instrument noise for each band is from Griffin et al. (2010), scaled by the square root of the number of repetitions in each band in the Dark Field. All noises and limits are given in mJy.

SPIRE band	|$\sigma _{\text{inst}}$| (mJy)	|$\sigma _{\text{limit}}$| (mJy)	|$\sigma _{\text{conf}}$| (mJy)
250 |$\mu$|m (PSW)	0.379	0.00643	5.8 |$\pm$| 0.3
350 |$\mu$|m (PMW)	0.316	0.00714	6.3 |$\pm$| 0.4
500 |$\mu$|m (PLW)	0.455	0.0151	6.8 |$\pm$| 0.4

Generally speaking, for an image containing a total of |$N_{\rm pix}$| pixels, we might assume the fluxes at each pixel is an independent measurement. However, as Scheuer (1957) notes, this assumption is invalid for confused fields with low angular resolutions. Instead, since the SPIRE beam is larger than a single pixel, the emission from a point source will be spread over multiple pixels. As a result, there is a non-zero correlation between pixel intensities. So, our standard error (or limiting sensitivity of the P(D) technique) becomes |$\sigma _{\rm inst}/\sqrt{N_{\rm beams}}$|⁠. Since a P(D) analysis considers the statistical distribution of pixel intensities (see Section 3) rather than attempting to identify and extract sources, we are able to reach below the confusion limit, and the instrument noise, all the way down to the limiting sensitivity. Table 1 presents the calculated limiting sensitivity for each SPIRE band.

3 P(D) ANALYSIS

To extend previous galaxy survey work done on the Deep Field (Pearson et al. 2025), we have conducted a P(D) – probability of deflection – analysis on the region. A P(D) analysis considers the probability distribution of flux deflections above/below the mean flux level in an observation, due to an underlying population of faint-flux sources (Scheuer 1957). This distribution is built from a given source count model, and the point spread function of the band/instrument in question (Barcons & Fabian 1990). We can therefore vary this input source count model to match the observed distribution of flux deflections as seen in the Deep Field observations. This technique has already been used for Herschel-SPIRE observations by Glenn et al. (2010) to great effect, but on a significantly shallower field.

3.1 P(D) Method

Following the method of Glenn et al. (2010), we make use of the pofd_affine package¹ for our P(D) analyses. This software represents the source count models as splines. These splines are made up of knots, where each knot defines a Euclidean normalized differential source count value [Jy|$^{1.5}$| deg|$^{-2}$|] (the knot value) at a specified flux [mJy] (the knot flux). The spline interpolation smoothly connects these knots.

The pofd_affine package uses a Markov Chain Monte Carlo (MCMC) method to optimize the source count model to match a given observed flux deflection distribution. This flux deflection distribution is calculated by taking a histogram of the pixel intensity values from an observation image.

We initially constrain the values of any knots – where the knot fluxes |$\gtrsim 2$| mJy – using the Glenn et al. (2010) results. For knots with fainter fluxes |$\lesssim 2$| mJy, where there are no longer any observed number count results, the initial source count values of the knots were chosen manually to follow the trends predicted by the current literature models (Negrello et al. 2017; Pearson et al. 2025). These knot values were then allowed to evolve from these initial conditions, according to the MCMC optimization routine used by the pofd_affine package. On the order of |${\sim }100$| optimizations were run with these different input conditions to ensure that our optimization had found the global optimal solution, and not just local minima.

To assess the goodness of fit of our P(D) results, we calculate the reduced chi-squared (RCS) statistic for the predicted pixel flux distribution and the observed pixel flux distribution. Additionally, we compared our best-fitting results with the results predicted by current literature models Béthermin et al. (2017), Negrello et al. (2017), and Pearson et al. (2025). These comparisons are presented in Fig. 1 and Table 2.

Figure 1.

P(D) results. Top row: Comparison between observed pixel flux distribution (black histogram) with flux distribution predicted by the (Negrello et al. 2017) model (red curve). Bottom row: Comparison with the best-fitting splines presented in this paper.

Table 2.

RCS statistics are presented for three literature models for each of the three SPIRE bands. None of the current literature models are able to accurately describe the Herschel-SPIRE Dark Field P(D) results. Sorted in order of descending Avg. |$\chi _\nu ^2$|⁠.

	|$\chi _{\nu }^{2}$|
Model	250 |$\mu$|m	350 |$\mu$|m	500 |$\mu$|m	Avg.
Béthermin et al. (2017)	18.44	15.99	16.00	16.81
Negrello et al. (2017)	22.77	6.92	4.31	11.33
Pearson et al. (2025)	16.17	7.11	4.49	9.26
Our P(D) analysis	1.45	1.26	1.38	1.36

We estimate the uncertainty in our fitted source count spline result using the following procedure. Taking a copy of the best-fitting spline, randomized displacements in the knot value are applied to each knot. The RCS statistic is recalculated for this perturbed spline. If the new spline’s RCS value lies within |$1\sigma$| of the best-fitting spline’s RCS statistic, we record the knot values for each knot flux. Starting again from a new copy of the best-fitting spline, this process is repeated until we have accumulated enough accepted splines that we can take statistics of the distribution of the possible values for each knot. Our estimated |$1\sigma$| uncertainty region is then given by the mean and standard deviation of each distribution of knot fluxes for each knot and interpolated linearly between fluxes in log space since we treat each knot independently in this procedure.

3.2 P(D) Results

Fig. 1 presents comparisons between the observed flux distributions from the Deep Region of the Herschel-SPIRE Dark Field with the predicted flux deflection distributions from our P(D) analysis, and the predictions of the Negrello et al. (2017) galaxy source count model. In all three bands, the predicted distribution from the Negrello et al. (2017) source count model struggles to fit the observed fluctuations at faint fluxes. Therefore, to fit our observed pixel distributions these models will need modifying. Following the method outlined in Section 3.1, we produced best-fitting P(D) splines for the Dark Field. The predicted distributions from our best-fitting P(D) splines, as seen in Fig. 1, match the observed pixel distributions (unlike the literature models). These qualitative conclusions are supported by the relevant RCS values for these fits (seen alongside the plots in Fig. 1, and in Table 2).

Table 2 presents the RCS statistics (goodness of fit) for our best-fitting source count splines from our P(D) analysis, along with comparisons to three current source count models from the literature. A well-fitted model will have a RCS, |$\chi _{\nu }^2$|⁠, close to 1. To represent an overall goodness-of-fit considering all three SPIRE bands, we compare the average |$\chi _{\nu }^2$| for each model. As is clear from the reduced chi-square statistics and the visual comparison with the observed pixel flux distribution, none of the current literature source count models can accurately describe our observations of the Herschel-SPIRE Dark Field. We compare the following:

• Béthermin et al. (2017): A phenomenological model, based on earlier work, with dark-matter haloes included to better match observed clustering.

• Negrello et al. (2017): Extends an earlier model by Cai et al. (2013) (which includes only late (spirals) and early (spheroids) populations), by including a population of lensed spheroids.

• Pearson et al. (2025): A backwards evolution infrared-based model which includes quiescent, starburst, luminous, ultraluminous, and AGN populations.

We generated expected flux deflection histograms for each of the three literature models above. Comparing these histograms with the observed flux deflection histograms (see Fig. 1 for an example), we see that none of the models can properly describe our Deep Field observations. Quantitatively, we can see that the RCS statistics for these fits agree that the models cannot describe the data. Table 2 presents these goodness-of-fit results.

Fig. 2 presents our best-fitting source count splines and regions of 1|$\sigma$| and 2|$\sigma$| uncertainty (grey-shaded regions) resulting from our P(D) analysis. Table 3 lists the values of the plotted spline knots. In all three SPIRE bands, we see a good agreement with the previously observed source counts down to |${\sim }$|20 mJy (Clements et al. 2010; Oliver et al. 2010). We also continue to see a good agreement at fainter fluxes to the stacked counts from Béthermin et al. (2012), and also the P(D) counts from Glenn et al. (2010), down to |${\sim }$|2 mJy. At fluxes |$\gtrsim$|2 mJy, we also see good agreement with the source count models from the literature (Béthermin et al. 2017; Negrello et al. 2017; Pearson et al. 2025).

Figure 2.

Differential source counts for all three SPIRE bands (250, 350, 500) |$\mu$|m. Observed differential source count results (black points) show very good agreement at bright fluxes. (Clements et al. 2010; Oliver et al. 2010; Béthermin et al. 2012). Prior SPIRE P(D) results are also presented (Glenn et al. 2010, black squares), and demonstrate our field is capable of reaching much fainter fluxes. Our P(D) results (grey regions) show a bump-turnover feature going towards fainter fluxes. Our P(D) results deviate from current literature galaxy evolution models (solid, dashed, and dotted lines) at these fainter fluxes (Béthermin et al. 2017; Negrello et al. 2017; Pearson et al. 2019). We indicate the instrument noise (⁠|$\sigma _{\text{inst}}$|⁠) and limiting sensitivity (⁠|$\sigma _{\text{inst}}/\sqrt{N}$|⁠) fluxes for reference (grey dotted and dashed lines, respectively).

Table 3.

Euclidean normalized best-fitting and 1|$\sigma$| differential source count results from our P(D) analyses in all three Herschel-SPIRE bands. The best-fitting results (along with the 1|$\sigma$| uncertainty region bounded by the 1|$\sigma$| lower and upper values) are plotted in Fig. 2.

250 |$\mu$|m				350 |$\mu$|m				500 |$\mu$|m
	Best fit	|$1\sigma$|			Best fit	|$1\sigma$|			Best fit	|$1\sigma$|
Flux	|${}^{\text{dN}}/_{\text{dS}} \cdot S^{2.5}$|	Lower	Upper	Flux	|${}^{\text{dN}}/_{\text{dS}} \cdot S^{2.5}$|	Lower	Upper	Flux	|${}^{\text{dN}}/_{\text{dS}} \cdot S^{2.5}$|	Lower	Upper
|$[\log _{10} (\text{Jy})]$|	|$[\log _{10} (\text{Jy}^{1.5} \text{deg}^{-2})]$|			|$[\log _{10} (\text{Jy})]$|	|$[\log _{10} (\text{Jy}^{1.5} \text{deg}^{-2})]$|			|$[\log _{10} (\text{Jy})]$|	|$[\log _{10} (\text{Jy}^{1.5} \text{deg}^{-2})]$|
|$-$|4.30	|$-$|17.0	|$-$|20.7	|$-$|16.5	|$-$|4.00	|$-$|15.0	|$-$|26.8	|$-$|16.8	|$-$|4.82	|$-$|1.98	|$-$|6.90	|$-$|1.50
|$-$|4.00	|$-$|6.09	|$-$|8.22	|$-$|5.65	|$-$|3.80	|$-$|2.82	|$-$|8.13	|$-$|2.74	|$-$|4.50	0.110	|$-$|3.24	0.256
|$-$|3.67	|$-$|1.19	|$-$|1.87	|$-$|0.575	|$-$|3.55	1.30	|$-$|1.77	1.16	|$-$|4.00	0.330	|$-$|0.882	0.552
|$-$|3.42	0.884	0.115	0.772	|$-$|3.30	|$-$|0.0623	|$-$|1.13	0.0868	|$-$|3.50	0.260	|$-$|0.266	0.185
|$-$|3.29	1.17	0.547	0.962	|$-$|3.00	|$-$|0.900	|$-$|1.02	|$-$|0.483	|$-$|3.00	0.0500	|$-$|0.152	0.0976
|$-$|3.06	0.712	0.723	0.848	|$-$|2.70	0.0400	|$-$|0.216	0.0659	|$-$|2.70	0.0250	0.0408	0.268
|$-$|2.70	0.241	0.396	0.596	|$-$|2.30	0.567	0.805	1.01	|$-$|2.30	0.440	0.384	0.543
|$-$|2.30	0.584	0.483	0.684	|$-$|2.00	0.810	0.679	0.747	|$-$|2.00	0.470	0.480	0.564
|$-$|2.00	0.911	0.907	0.996	|$-$|1.70	0.809	0.616	0.807	|$-$|1.70	0.209	0.0396	0.254
|$-$|1.70	0.918	0.868	0.996	|$-$|1.35	0.311	0.311	0.311	|$-$|1.35	|$-$|0.353	|$-$|0.353	|$-$|0.353
|$-$|1.35	0.588	0.588	0.588	|$-$|1.00	|$-$|0.668	|$-$|0.668	|$-$|0.668	|$-$|1.00	−1.36	|$-$|1.36	|$-$|1.36
|$-$|1.00	0.0268	0.0268	0.0268	|$-$|0.699	|$-$|0.864	|$-$|0.864	|$-$|0.864	|$-$|0.699	|$-$|1.50	|$-$|1.50	|$-$|1.50
|$-$|0.699	|$-$|0.375	|$-$|0.375	|$-$|0.375	|$-$|0.347	|$-$|0.906	|$-$|0.906	|$-$|0.906	|$-$|0.347	|$-$|1.50	|$-$|1.50	|$-$|1.50
|$-$|0.347	|$-$|0.439	|$-$|0.439	|$-$|0.439	0.000	|$-$|0.950	|$-$|0.950	|$-$|0.950	0.000	|$-$|1.50	|$-$|1.50	|$-$|1.50
0.000	|$-$|0.290	|$-$|0.290	|$-$|0.290	−	−	−	−	−	−	−	−

However, at fluxes |$\lesssim$|2 mJy – below the limit of previously observed counts – our P(D) results begin to deviate from the trend predicted by the literature models. Going to fainter fluxes, as we approach the instrument noise limit in all three bands we see a significant increase in our P(D) source count results compared with the literature models. In the 350 |$\mu$|m band, we also see a significant dip in the source counts (at |${\sim }$|1 mJy) preceding this steep rise. In all three bands, the literature models do not predict this trend.

Below the instrument noise level (Section 2), we see a steep drop off in the Euclidean normalized source counts as we go to fainter fluxes in the 250- and 350-|$\mu$|m bands. The resulting ‘bump-turnover’ shape in the source counts peaks at |${\sim }0.5$| mJy in the 250 |$\mu$|m band and at |${\sim }0.3$| mJy in the 350 |$\mu$|m band. The turnover seen is also much more sudden than the gentle roll-off predicted by the literature models. In previous, complementary work on the CIB, Duivenvoorden et al. (2020) were similarly able to push below the instrument noise level using a stacking technique.

In the PLW band, we see a similar albeit slightly different trend below the instrument noise level. Instead of a sharp peak near the instrumental noise limit, we see a gentle rise above the predictions given by the literature models, leading into a roughly steady Euclidean normalized source count that does not change with decreasing flux. This trend peaks slightly at |${\sim }0.05$| mJy, before turning over and beginning a rapid decrease in source counts with decreasing flux. We are not able to measure the full turnover due to the cut-off imposed by the limiting sensitivity. This turnover in the counts in the PLW band happens at much lower (by around an order of magnitude) fluxes than in the PSW and PMW bands. The peaks in the source count models observed in each band occur at lower fluxes with increasing wavelength. We will discuss the potential interpretations of these results in Section 5.

In each band, we see our uncertainty regions follow the general trend of the best-fitting spline. The region also widens as we go to fainter fluxes. Strangely, however, our best-fitting spline deviates from the uncertainty region, most notably in the PSW band. We hypothesize this is due to how we have treated the spline knots when estimating the uncertainty region. The knots in the best-fitting spline are highly correlated and must represent a coherent spline whereas the points defining the uncertainty region are treated independently. These uncertainty regions are ranges in number counts where a spline could lie which does not contradict with the observations at that sigma level. For our best-fitting spline, this coincides with the point not quite matching in the histogram and suggests that there may be more structure in the number counts than we can reproduce. It is important to emphasize that our best-fitting spline is not a rigorous physical model, but a fit to the data.

Another possible source of uncertainty could be in the mean-subtraction performed on the Dark Field maps. If the mean of the map is incorrect, this will add an offset to all pixel flux values in the map, and could affect the differential counts. An uncertainty of 0.07 MJy sr⁻¹ at 250 microns corresponds to an offset of approximately 0.5 mJy beam⁻¹. Adjusting the mean by this value could potentially move the counts from the peak in the best-fitting spline at 0.4 mJy to the lowest flux knot. However, considering the pixel histograms in Fig. 1, a mean offset would equate to a lateral shift of the histogram in flux, S. As we can see, with the current literature models we still cannot explain our observations by simply shifting the distributions in S. We still require new models to produce a different histogram shape overall. This can be quantified by comparing the reduced chi-squared statistics for maps with/without an offset applied to the mean. For our best-fitting PSW spline, applying a shift of |$\pm 0.5$|mJy beam⁻¹ increases the reduced chi-squared statistic for our fit from |${\sim }1.5$| (with no mean shift) to |${\sim }6.5$| (with the shift applied). Table 4 summarizes the results for the Negrello et al. (2017) model in all three bands. In all cases, adding a shift to the mean produces a higher reduced chi-squared statistic, and as such a less good fit to the observations.

Table 4.

Quantitive goodness-of-fit tests on whether a shift applied to the mean could explain the poorly-fitting literature models. For each row, a constant offset – |$\pm 0.5$| mJy/beam – is applied to the maps before the histogram is calculated and compared to the literature model presented in Fig. 1 (Negrello et al. 2017).

	Reduced chi-squared, |$\chi _\nu ^2$|
Mean shift (mJy beam−1)	250 |$\mu$|m	350 |$\mu$|m	500 |$\mu$|m
|$-$|0.5	30.30	10.97	5.66
0.0	22.77	6.92	4.31
+0.5	23.78	7.77	6.70

4 CONTRIBUTION TO THE CIB

The emission from a large population of galaxies contributes to a cosmic (infrared) background. Using our best-fitting |${\rm d}N/{\rm d}S$| results from our P(D) analysis, we can integrate the |${\rm d}N/{\rm d}S$| curve to obtain a contribution to the extra-galactic background intensity from our population,

$$\begin{eqnarray} \int S \cdot \frac{{\rm d}N}{{\rm d}S} \cdot {\rm d}S, \end{eqnarray}$$

(1)

where |${\rm d}N/{\rm d}S$| are the differential (not Euclidean normalized) source counts. Integrating over the entire flux range finds the total emission from the entire population of faint galaxies. Comparing this with independent measurements of the CIB then enables us to verify the results of our P(D) analysis. With pofd_affine, it is possible to constrain the fit using a known integrated CIB value. However, we chose not to do this so that we could independently verify the result using observations of the same region from the Planck satellite.

4.1 Comparison with Planck Observations

The high-frequency instrument (HFI) aboard Planck captured all-sky maps² at 545 GHz (550 |$\mu$|m) and 857 GHz (350 |$\mu$|m; Ade et al. 2014, Planck Collaboration 2020). Fig. 3 shows the 857 GHz (350 |$\mu$|m) map with the locations of the Dark Field region and other SPIRE fields indicated. Using these maps, we cut out the region corresponding to the Deep Region within the Herschel-SPIRE Dark Field, and integrated the Planck measurements taken within this region. This gave us a total intensity in MJy sr⁻¹ for the region of sky in which our P(D) |${\rm d}N/{\rm d}S$| models were fitted.

Figure 3.

An all-sky map, in galactic coordinates, at 857 GHz (or 350 |$\mu$| m) from the Planck satellite. Locations of various survey fields (including the Dark Field) are indicated. Lockman-SWIRE (cyan square), GOODS-N (yellow triangle), and the Herschel-SPIRE Dark field (white circle).

In order to compare these observations with our Herschel-SPIRE data, we converted the 545 GHz (550 |$\mu$|m) maps to 500 |$\mu$|m by assuming the intensity ratio of the CIB as measured by Planck at these two wavelengths was the same as the ratio of the CIB as measured by the FIRAS instrument aboard COBE at these two wavelengths (Fixsen et al. 1998). This is in essence a colour correction. At a wavelength of |$\lambda$||$\mu$|m, we measure an intensity of |$I_{F,\lambda }$| and |$I_{P,\lambda }$| with COBE-FIRAS and Planck, respectively. So, we have

$$\begin{eqnarray} \frac{I_{F,550}}{I_{F,500}} = \frac{I_{P,550}}{I_{P,500}} \longrightarrow I_{P, 500} = \frac{I_{P, 550} \cdot I_{F, 500}}{I_{F, 550}} \end{eqnarray}$$

(2)

and using this we can calculate the intensity at 500 |$\mu$|m, as would be measured by the Planck and Herschel satellites. Similarly, we also colour-corrected the Planck results from 350 |$\mu$|m to 250 |$\mu$|m, to enable a direct comparison with the Herschel-SPIRE PSW results.

4.2 CIB Results

Table 5 presents the results from our integration of both the SPIRE P(D) fit and the Planck data. This includes the |$250{-}$| and |$500{-}\mu$|m colour-corrected Planck data, as described in Section 4.

Table 5.

Integrated CIB results for each SPIRE band: PSW (250 |$\mu$|m), PMW (350 |$\mu$|m), and PLW (500 |$\mu$|m). We have the integrated CIB from our best-fitting source count model for each band and for the upper and lower 1|$\sigma$| uncertainty regions. We take the bounding |${\rm d}N/{\rm d}S$| curve for this uncertainty region and integrate it to find the lower and upper bounds. We also present an integrated CIB result from the Planck satellite of the Deep Region. These independent measurements lie within our 1|$\sigma$| bounds for the SPIRE results, and we see the best agreement in the 350 |$\mu$|m band.

	Intensity (MJy sr−1) in SPIRE band
	250 |$\mu$|m	350 |$\mu$|m	500 |$\mu$|m
1|$\sigma$| Lower bound	1.56|$\pm$|0.05	0.53|$\pm$|0.02	0.54|$\pm$|0.02
Best fit	2.27|$\pm$|0.07	2.20|$\pm$|0.07	2.59|$\pm$|0.08
1|$\sigma$| Upper bound	2.73|$\pm$|0.08	2.63|$\pm$|0.08	4.0|$\pm$|0.1
Planck (Deep Region)	2.6|$\pm$|0.6	2.01|$\pm$|0.07	0.98 |$\pm$| 0.11

All three of the Planck measurements lie within the |$1\sigma$| region bounded by our SPIRE results. However, we get a closer match between the SPIRE best-fitting result (⁠|$2.20\pm 0.07$|⁠) and the Planck result (⁠|$2.011\pm 0.069$|⁠) in the 350 |$\mu$|m band. According to our SPIRE best-fitting results, the deep region is redder (more emission in the longer wavelength bands) than what was measured byPlanck.

5 DISCUSSION

Without pushing to fainter fluxes, our results agree well with the current literature models for fluxes |$\gtrsim 20$| mJy. However, it is the unique depth of the Herschel-SPIRE deep region – that we can exploit with statistical techniques like P(D) – that has enabled us to explore fainter fluxes than ever before with the Herschel-SPIRE instrument. At these fainter fluxes |$\lesssim 20$| mJy, the bump-then-turnover in the source counts models found via the P(D) analysis indicates our observations account for 100 per cent of the FIR galaxies at 250 and 350 |$\mu$|m. Duivenvoorden et al. (2020) applied a novel source extraction technique to the COSMOS field (Scoville et al. 2007) and found results consistent with recent literature galaxy count models and the Fixsen et al. (1998) CIB measurements. The CIB convergence and agreement with the Fixsen et al. (1998) result suggests that Duivenvoorden et al. (2020) have observed all the galaxies within the COSMOS field between 250 |$\mu$|m and 500 |$\mu$|m.

Our results are also – at faint (⁠|$\lesssim 20$| mJy) fluxes – distinct from the current literature models. The bump-turnover shape, which cannot be explained by any of the current models, indicates that something is missing from these models. This could be an additional contribution from a new population of galaxies, or a population that is already known but simply absent in the models. Perhaps it could be an evolution mechanism that is not included or poorly understood. With the data we have, we are unable to conclude exactly what this missing component is. Further observations or modelling is necessary to discover this, but this is beyond the scope of this paper.

We also notice that the peak of the bump in the source counts shifts to fainter fluxes as the wavelength increases. If this bump is due to an unaccounted-for population in the literature models, this shifting towards fainter fluxes with increasing wavelength suggests the population might be bluer in colour (brighter at shorter wavelengths).

As with the Duivenvoorden et al. (2020) results, the close agreement of our 350 |$\mu$|m CIB result with that from the Planck satellite is yet another indication (in addition to the sharp turnover in the source counts) that we have observed all the FIR galaxies at this wavelength within the Herschel-SPIRE Deep Region.

Although the Planck measurement of the CIB at 500 |$\mu$|m lies comfortably within the region bounded by the lower and upper |$1\sigma$| bounds, our SPIRE measurement from integrating the best-fitting source count spline predicts a value almost two and a half times larger. It is important to note that the uncertainty region for our 500 |$\mu$|m spline result (Fig. 2) is the largest of all three bands.

6 CONCLUSIONS

Using data from the Herschel-SPIRE Dark Field (constructed from stacked calibration maps), we identified a region of uniform coverage at the centre we call the Deep Region (Pearson et al. 2025). This region is highly confused due to the large beam size of the Herschel-SPIRE instrument, and as such, standard source extraction techniques are limited. Instead, we performed a statistical P(D) analysis of the region (Scheuer 1957; Glenn et al. 2010). A P(D) analysis attempts to fit a distribution of flux deflections due to an underlying population of faint-flux sources (galaxies).

Comparing current literature models with our best-fitting |${\rm d}N/{\rm d}S$| source counts from the P(D) analysis, we found that none of the current models were able to describe our data, particularly at faint (⁠|$\lesssim 20$| mJy) fluxes. In all three bands, our best-fit spline had a bump-turnover shape going to lower fluxes. A bump – above the trend in the source counts predicted by current literature models – was observed at |${\sim }0.5$|⁠, |${\sim }0.3$|⁠, and |${\sim }0.05$| mJy in the 250, 350, and 500 |$\mu$|m bands, respectively. Both the 250 and 350 |$\mu$|m bands exhibited a rapid turnover in the source counts, indicating that we had observed the entire population of galaxies in those two bands with our observations. In the 500 |$\mu$|m band, we still see a turnover, albeit at much fainter fluxes, and less steep than in the other two shorter-wavelength bands.

Integrating our best-fitting source counts splines, we also saw the CIB contribution agrees with independent Planck observations at 350 |$\mu$|m. At 500 |$\mu$|m, the Planck emission lies within the |$1\sigma$| bounds of our Herschel-SPIRE measurement, although it agrees less closely. Our 500 |$\mu$|m CIB measurements are less certain, likely due to most of the emission being at such faint fluxes (⁠|$\lesssim 0.1$| mJy) and therefore our uncertainty on these values is greater. However, the turnover and agreement with the CIB in the 350 |$\mu$|m band does suggest that we have detected the entire 350 |$\mu$|m population of galaxies in this region with our observations.

The Deep Region overlaps deep MIPS and IRAC observations in the Spitzer IRAC validation field (Krick et al. 2008), and while this already provides a wealth of ancillary data, extending the observations of the Deep Region to submillimetre wavelengths via ground-based facilities may help to further constrain the fainter population through higher resolutions. The Herschel-SPIRE Dark Field observations presented here are the deepest far-IR observations currently available, and they will remain so until the next generation of space-based far-IR telescopes start to operate. Such telescopes will be necessary for us to better understand the nature and origin of the sub-mJy bump in counts found by our P(D) analysis. The ideal instrument for this task should be both more sensitive than Herschel-SPIRE and also capable of operating at a higher angular resolution so that it is not as confusion-limited in the far-IR. Several far-IR mission concepts are currently being studied by NASA to fulfil the far-IR Probe recommendations of the US Decadal Review (National Academies of Sciences, Engineering, and Medicine 2023), most of which will be able to achieve some of these capabilities either directly, or through the use of ancillary data at other wavelengths to break confusion. Until such a future mission flies we are limited to using data from Herschel-SPIRE, the deepest of which, as discussed here, is found in the Herschel-SPIRE Dark Field.

ACKNOWLEDGEMENTS

SPIRE has been developed by a consortium of institutes led by Cardiff University (UK) and including University of Lethbridge (Canada); NAOC (China); CEA, LAM (France); IFSI, University of Padua (Italy); IAC (Spain); Stockholm Observatory (Sweden); Imperial College London, RAL, UCL-MSSL, UKATC, University of Sussex (UK); and Caltech, JPL, NHSC, University of Colorado (USA). This development has been supported by national funding agencies: CSA (Canada); NAOC (China); CEA, CNES, CNRS (France); ASI (Italy); MCINN (Spain); SNSB (Sweden); STFC, UKSA (UK); and NASA (USA). HIPE is a joint development by the Herschel Science Ground Segment Consortium, consisting of ESA, the NASA Herschel Science Center, and the HIFI, PACS and SPIRE consortia.

Our work also made use of observations from the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory (JPL), California Institute of Technology, under a contract with the National Aeronautics and Space Administration (NASA).

Our research made use of the NASA/IPAC Infrared Science Archive, which is funded by the National Aeronautics and Space Administration (NASA) and operated by the California Institute of Technology Funding for this work was provided, in part, by the Science and Technology Facilities Council (STFC).

We acknowledge computational resources and support provided by the Imperial College Research Computing Service (http://doi.org/10.14469/hpc/2232). We also acknowledge support via the Rutherford Appleton Laboratory (RAL) Space In House Research programme funded by the Science and Technology Facilities Council (STFC) of the UK Research and Innovation (UKRI) (award ST/M001083/1).

Our paper makes use of perceptually uniform, colour-vision-deficiency-friendly colour maps (Crameri 2023) to prevent the visual distortion of the data and promote accessibility.

The authors would also like to thank the anonymous referee for their comments and suggested improvements to the paper.

DATA AVAILABILITY

The processed Herschel-SPIRE Dark Field observations (including the Deep Region) are available on Zenodo at https://dx.doi.org/10.5281/zenodo.14846645, and are discussed in more detail in Pearson et al. (2025). P(D) and |${\rm d}N/{\rm d}S$| model files, along with plotting scripts (python jupyter notebooks) for the data analysed in this study are also available on Zenodo under the same directory. The best-fitting source count splines for all three SPIRE bands (250, 350, 500 |$\mu$|m) are presented in Tables 3, and are available in full online. Data sets for the CIB analysis were derived from sources in the public domain. The Planck all-sky maps are accessible athttps://irsa.ipac.caltech.edu/data/Planck/release_3/all-sky-maps/. The IRAS ISSA all-sky maps can be found at https://lambda.gsfc.nasa.gov/product/iras/iras_thirdp_prod_table.cfm

Footnotes

1

The pofd_affine package was developed by A. Conley and is available at https://github.com/aconley/pofd_affine.

2

The Planck all-sky maps can be found at https://irsa.ipac.caltech.edu/data/Planck/release_3/all-sky-maps/. See also: Planck Collaboration (2020).