Candidate high-redshift protoclusters and lensed galaxies in the Planck list of high-z sources overlapping with Herschel-SPIRE imaging

Planck mapped the whole sky between 30 and 857 GHz, including channels centred at 217, 353 and 545 GHz (Planck Collaboration I, 2014). Using these maps, the Planck Collaboration generated a list of high-$z$ candidates by identifying red peaks in the CIB (Planck Collaboration XXXIX, 2016). Specifically, maps of the cleanest 25.8 per cent of the sky (i.e., the least contaminated by Galactic emission) were used, and the Infrared Astronomical Satellite (IRAS) 3-THz map was included in the analysis. The CMB was removed from each channel map using the 217-GHz map as a template, and Galactic cirrus was removed using the IRAS 3-THz map. S/N$\,{>}\,$5 sources were identified in the 545-GHz ‘excess’ map, constructed by subtracting a linear interpolation between the 353- and 857-GHz maps from the 545-GHz map. In addition, sources were required to have S/N$\,{>}\,$3 at 353 and 857 GHz. Lastly, to avoid selecting cold Galactic clumps and radio galaxies, all sources with $S_{545}/S_{857}\,{<}\,0.5$ or $S_{353}/S_{545}\,{>}\,0.9$ were removed. The resulting catalogue consists of 2151 potential high-$z$ sources, spread across the northern and southern Galactic hemispheres (see Fig. 1). The focus on picking out the apparently coldest, potentially highest redshift objects makes the PHz catalogue complementary to the Planck compact source catalogues, with almost no overlap (Planck Collaboration XXXIX, 2016).

Across three proposal periods, a total of 228 Planck-identified sources were followed up with dedicated Herschel-SPIRE observations (PlanckXXVII): 204 sources were selected with an algorithm that resembled the final PHz selection process and the other 24 came from the PCCS. The 228-source sample was constructed based on preliminary and evolving Planck data, as new products were made available throughout the three SPIRE observation periods. Additionally, there was a bias for preferentially selecting bright sources with regular shapes to follow up with Herschel-SPIRE (Planck Collaboration XXXIX, 2016). As a result, many of the sources catalogued in PlanckXXVII did not end up in the final PHz list, and the sample was biased towards sources with high S/N, small angular sizes, small ellipticities and low extinctions when compared with the full PHz. Nonetheless, the selection process for the 228 sources was qualitatively similar to the PHz catalogue, requiring compact (at 5-arcmin resolution) source-detections at 545 GHz.

Specifically, of the 228 Planck-identified sources presented in PlanckXXVII, 83/228 (36 per cent) ended up in the final PHz catalogue. Due to these differences, the authors of Planck Collaboration XXXIX (2016) remarked that ‘it is hard to draw definitive conclusions about the nature of the PHz sources based on the Herschel analysis’. In what follows, we study a large and unbiased sample of PHz sources using recently homogenized archival Herschel-SPIRE observations. For comparison, we repeated all of our analyses on the 228 PlanckXXVII sources (although in practice, we were unable to recover the Herschel-SPIRE maps of nine of the sources). Hereafter, the 228 Planck-identified sources followed up in PlanckXXVII will be referred to as the ‘PlanckXXVII sources’.

The Herschel satellite carried out a number of wide-field extragalactic surveys using the SPIRE instrument (see Fig. 1 and Table 1). HELP combines the Herschel-SPIRE survey fields into a consistent single data release spanning a total of 1270 deg² on the sky (Shirley et al., 2021a). In particular, HELP DR1 provides homogeneous SPIRE 250-, 350- and 500-$\mu$m images for each of the HELP fields. The HELP fields avoid regions of bright dust emission from our own Galaxy, span a range of right ascension and include areas of sky near both ecliptic poles. As a result, there is significant overlap between the coverage of HELP fields and sources in the PHz catalogue, allowing us to select an unbiased sample of PHz sources at the SPIRE wavelengths.

A total of 194 sources from the PHz list fall within the HELP fields, six of which lie too close to the edge of the HELP field coverage and one falls within a region masked by HELP because of Galactic cirrus, resulting in a total of 187 usable sources. Hereafter, these 187 PHz sources that lie in the HELP fields will be referred to as the ‘PHz HELP sources’. By design of the follow up in PlanckXXVII, there is no overlap between the PHz HELP sources and the PlanckXXVII sources. However, the PHz HELP sample includes two sources that have previously been studied extensively, namely PHz G95.50$-$61.59 (Flores-Cacho et al., 2016) and PHz G237.0$+$42.5 (Koyama et al., 2021; Polletta et al., 2021).

In order to study the number counts of SPIRE-detected galaxies in the PHz catalogue, we define an ‘IN region’ in the same way as in PlanckXXVII.¹¹1To avoid confusion, we will reserve the word ‘source’ for the Planck-selected targets and refer to the SPIRE-detected objects as ‘galaxies’. We draw contours enclosing 50 per cent of the Planck flux density at 545 GHz (the frequency used to select PHz sources), and use the galaxies within this region for further analysis. Other follow-up studies have adopted a circular IN region centred on the Planck position (e.g., Greenslade et al., 2018). We find that using a contour on the Planck 545-GHz map generally encloses sources more effectively than a circular IN region. That is, when the galaxies making up a PHz source cover an extended region, Planck partially resolves the source and the contour is slightly larger than the Planck beam, whereas for individual bright galaxies, the contour is approximately the size of the beam (see Fig. 2). As a result, defining the IN region using a Planck contour misses fewer contributing galaxies than a circular IN region (although the differences are not dramatic). For comparison, the contours have a median area of 59.2 arcmin² on the sky, corresponding to a search radius of 4.34 arcmin, which is slightly smaller than the 4.63 arcmin radius adopted in Greenslade et al. (2018). We determine the IN region contour by smoothing the 545-GHz map with a 2.0 arcmin full width at half-maximum (FWHM) Gaussian, then find the contour on the smoothed 545-GHz map corresponding to 50 per cent of the peak flux density of the source. For sources where the contour does not close or extends far beyond the Herschel-SPIRE map coverage, a circular IN region centred on the Planck position with a radius of 4.34 arcmin is adopted, corresponding to the median contour area.

Before performing photometry, we measure the noise at each wavelength across all of the HELP fields in a similar fashion to PlanckXXVII (Table 1). Specifically, in each HELP field, for each SPIRE band, we construct histograms including all of the pixels in the map. The histograms have a mean of 0 mJy, and appear Gaussian in shape, except for a positive tail of bright pixels corresponding to detected galaxies. We fit a Gaussian to each histogram in order to calculate the standard deviation of the pixel value distribution, $\sigma$, which quantifies the noise of the map (i.e., the combined confusion and detector noise). Note that all pixels, including those in the positive tail, are included when performing the Gaussian fit, thereby accounting for confusion noise (if instead the fit is restricted to the negative half of the histogram, the RMS noise levels change only slightly). For consistency with previous studies (e.g., PlanckXXVII; Greenslade et al., 2018), we use the same flux density thresholds for galaxy detection across all HELP fields. Because the noise levels of the HELP fields vary (Table 1), we adopt 250, 350 and 500 $\mu$m flux density thresholds that are significantly larger than the noise level of all the fields (${\gtrsim}\,3\,\sigma$). Specifically, we choose flux density cutoffs of 4$\bar{\sigma}$ at 250, 350 and 500 $\mu$m, that is, 4 times the mean noise of all the HELP fields. The cutoff flux densities are 32.0, 34.6 and 38.5 mJy at 250, 350 and 500 $\mu$m, respectively. These thresholds are comparable to those adopted in PlanckXXVII, although they are higher than the levels used in Greenslade et al. (2018), where 25.4 mJy was used at all three wavebands. We find that our results are not sensitive to the choice of flux density thresholds.

We detect Herschel-SPIRE galaxies independently in the 250-, 350- and 500-$\mu$m maps before matching galaxies across the three bands. Photometry is performed on sky-subtracted mean-0 maps using Gaussian point-spread functions (PSFs) with FWHMs of 18.15, 25.15 and 36.30 arcsec for the 250-, 350- and 500-$\mu$m maps, respectively (Griffin et al., 2010). First, we convolve the signal map with the corresponding Gaussian PSF. Then, we search for peaks in the convolved signal map that fall within the IN region. All local peaks in the 250-, 350- and 500-$\mu$m convolved maps above the corresponding galaxy detection thresholds are independently identified as potential galaxies. The values and positions of peak pixels in the convolved maps are used as initial guesses for the galaxy fluxes and positions at the wavelength of detection. A Gaussian PSF fit is then performed on each detected galaxy in the non-convolved map, with a fixed FWHM (depending on the wavelength), and initial position/amplitude guesses from the convolved map. During PSF fitting, the PSF amplitude is allowed total freedom, but the positions are allowed to shift only slightly from the initial position ($0.5\,{\times}\,$FWHM in each direction) to avoid accidentally fitting the PSF to a separate nearby galaxy instead of the intended galaxy. Finally, the amplitudes and positions of the fitted Gaussian PSFs are added to the 250/350/500 $\mu$m galaxy catalogue.

When galaxy blending is present, which is especially common in the overdense sources discovered with Planck, a more careful fitting procedure is required. Whenever multiple galaxies are detected within two FWHMs of each other (corresponding to 6.0, 5.1 and 5.2 pixels at 250, 350 and 500 $\mu$m, respectively), galaxy flux densities and positions are measured with simultaneous PSF fitting. That is, Gaussian PSFs are simultaneously fit to each of the nearby galaxies in the non-convolved map, avoiding the risk of inadvertently including the fluxes from other nearby galaxies during PSF fitting. In practice, twice the FWHM is a conservative threshold; galaxies separated by ${\gtrsim}\,2\,{\times}\,$FWHM are too far away for significant blending to occur, in which case simultaneously fitting the PSFs is equivalent to independently fitting Gaussian PSFs to both galaxies. After simultaneous PSF fitting, the amplitudes and positions of the fitted Gaussian PSFs are added to the 250/350/500 $\mu$m galaxy catalogue. We find that our galaxy detection and photometry values are largely consistent with the results from other comparable codes, including codes developed specifically for Herschel-SPIRE observations (e.g., sussextractor; Pearson et al. 2014).

Galaxies are matched across bands in order to measure single-galaxy properties (e.g., colours), as well as to identify candidate strong gravitationally-lensed galaxies. Matching galaxies detected across the 250-, 350- and 500-$\mu$m maps is done by focusing on the galaxies detected at 350 $\mu$m (also see PlanckXXVII). The 350-$\mu$m map has the advantage of higher angular resolution and lower noise than the 500-$\mu$m map, while also avoiding many of the low-$z$ galaxies that appear in the 250-$\mu$m map. Our band-merging procedure proceeds as follows. For each galaxy detected in the 350-$\mu$m map, if there is a galaxy within $0.5\,{\times}\,{\rm FWHM}$ of the 350-$\mu$m galaxy position in either the 250-$\mu$m map (FWHM 18.15 arcsec) or the 500-$\mu$m map (FWHM 36.30 arcsec), then the independently-calculated flux density is adopted from the corresponding map. If a galaxy is not detected within $0.5\,{\times}\,{\rm FWHM}$ in the 250- or 500-$\mu$m maps, we perform a ‘forced’ Gaussian PSF fit in the 250- or 500-$\mu$m maps (for which there are no enforced thresholds). That is, we use the 350-$\mu$m galaxy position as the initial guess position and the value in the nearest pixel in the convolved 250- or 500-$\mu$m map as the initial guess flux in order to fit a Gaussian PSF (as in Section 3.2). If multiple galaxies are detected within $0.5\,{\times}\,{\rm FWHM}$ of the 350-$\mu$m position at either 250 or 500 $\mu$m, we instead perform simultaneous forced Gaussian PSF fitting at the positions of the close galaxies in the other two maps (i.e., if a galaxy detected at 350 $\mu$m breaks up into multiple galaxies at 250 $\mu$m, simultaneous PSF fitting is carried out at the positions of the galaxies in the 350- and 500-$\mu$m maps). To be consistent with independent galaxy detection, when forced PSF fitting is performed, we carry out simultaneous PSF fits on all galaxies found within $2\times{\rm FWHM}$ of the galaxy being matched across maps. For instance, if a galaxy at 350 $\mu$m has a single corresponding galaxy within $0.5\,{\times}\,{\rm FWHM}$ at 250 $\mu$m, but there is a second galaxy $1.5\,{\times}\,{\rm FWHM}$ away at 250 $\mu$m, simultaneous Gaussian PSF fitting on both of the galaxies is performed at 350 $\mu$m and 500 $\mu$m, to avoid biases due to galaxy blending. If a 350-$\mu$m galaxy is successfully matched to a single galaxy at 250 $\mu$m, or is revealed to break up into multiple galaxies, we record the positions of the galaxy at 250 $\mu$m because of the superior angular resolution of the 250-$\mu$m map (otherwise, the 350-$\mu$m position is used). With this procedure, we generate a catalogue of galaxies with flux densities measured at 250, 350 and 500 $\mu$m.

As described in Section 2.1, the PHz sources were selected using Planck flux densities, or more specifically, the ‘excess’ flux density at 545 GHz. PlanckXXVII followed up a sample of 228 bright sources selected in a similar fashion as the PHz catologue, of which 83 ended up in the final version of the PHz catalogue. To evaluate the statistical differences resulting from this selection process, in Fig. 3 we show the distributions of the Planck flux densities for the 187 PHz HELP sources, about 100 of the PlanckXXVII sources (the sources for which Planck fluxes are publicly available) and all 2151 PHz sources. PlanckXXVII sources that came from earlier, non-public, PHz catalogue-like selection processes are necessarily excluded. The two-sample Kolmogorov–Smirnov (KS) test rejects the null hypothesis that the PlanckXXVII sources are drawn from the same sample as the full PHz catalogue with $p$-values of $9.7\times 10^{-8}$, $1.6\times 10^{-4}$ and $2.7\times 10^{-10}$ at 857, 545 and 353 GHz, respectively. For the PHz HELP sources, on the other hand, the two-sample KS test gives a $p$-value of $0.87$, $0.69$ and $0.90$ at 857, 545 and 353 GHz, respectively. In other words, the 187 PHz HELP sources are drawn from the same distribution as the entire PHz sample, and are therefore representative of the overall PHz catalogue. These trends are similar across other PHz source properties, such as FWHMs and ellipticities. It is worth noting that differences between the PlanckXXVII source flux densities and the PHz source flux densities (as well as other properties) are not simply caused by the inclusion of PCCS sources in the PlanckXXVII sample. We see an equally significant difference when we restrict the comparison to the 83 PlanckXXVII sources that ended up in the final PHz catalogue, indicating that there were biases in the selection of PlanckXXVII sources (see Planck Collaboration XXXIX 2016 for a discussion of these selection effects).

In Fig. 4 we show Euclidean-normalized differential number counts of SPIRE-detected galaxies, $S^{2.5}\,dN/dS$, where $N$ is the total number of galaxies per steradian in the IN regions and $S$ is the flux density at 250, 350 or 500 $\mu$m. Number counts are calculated using the galaxies detected independently in each of the three bands and are intended to be compared relative to one another, so no corrections (such as for completeness) have been applied. As in PlanckXXVII, we use the Lockman-SWIRE field to compare between the PHz HELP and PlanckXXVII sources. To compare the number counts of IN region galaxies with the number counts of random locations, we compute the number counts at 10,000 random positions in Lockman-SWIRE within a circle of radius 4.34 arcmin (the median size of the IN regions).

Similarly to the PlanckXXVII sources, when compared with random locations on the sky, the PHz HELP sources have a significant excess in number counts below approximately 100 mJy across the 250-, 350- and 500-$\mu$m bands, with the most substantial differences seen at 500 $\mu$m. This confirms that typical PHz sources contain an exceptional number of galaxies, just as the preferentially-selected PlanckXXVII sources do. That being said, the PHz HELP sources have lower number counts at 250 $\mu$m than the PlanckXXVII sources. Additionally, above about 100 mJy, there are substantial differences between the PHz HELP and PlanckXXVII number counts. As found in PlanckXXVII, the PlanckXXVII sources contain significantly higher numbers of SPIRE galaxies at large flux densities compared to random positions, the brightest of which are now known to be strong gravitational lenses. Across all three bands, the HELP sources have lower numbers of SPIRE galaxies at large flux densities compared to the PlanckXXVII sample, potentially indicating that fewer of the PHz HELP sources contain strong gravitationally-lensed galaxies. This is investigated further in Section 5.

We next quantify the projected overdensities of our PHz HELP sources. In Fig. 5 we compare the number densities of IN region galaxies for the PHz HELP sources, the PlanckXXVII sources and 10,000 random positions in the Lockman-SWIRE field. Just as we found in Section 4.2, the PHz HELP sources are significantly overdense when compared with random positions, across all three bands. The PlanckXXVII sources appear more overdense at 250 $\mu$m compared to the PHz HELP sources, whereas galaxy densisites appear more comparable at 350 and 500 $\mu$m.

After matching galaxies across the 250-, 350- and 500-$\mu$m bands, colours are quantified by comparing flux density ratios. Figure 7 shows the distributions of $S_{350}/S_{250}$ and $S_{500}/S_{350}$ for all of the galaxies within in the PHz HELP sources, the PlanckXXVII sources and 10,000 random Lockman-SWIRE positions. The PHz HELP sources display a clear excess of $S_{350}/S_{250}$ and $S_{500}/S_{350}$ ratios compared to random positions, indicative of an excess of red galaxies. The PlanckXXVII sources display a more significant excess of $S_{500}/S_{350}$ ratios compared to random positions; however, the distribution of $S_{350}/S_{250}$ is more consistent with random positions (these distributions are comparable to Fig. 7 from PlanckXXVII). Quantitatively, when PHz HELP colours are compared to random sources, we find that the null hypothesis is rejected at 20.3$\,\sigma$ and 14.2$\,\sigma$ for the $S_{350}/S_{250}$ colours and the $S_{500}/S_{350}$ colours, respectively. These differences demonstrate that the average PHz source contains an excess of red (and potentially higher $z$) sources, as was previously concluded for the preferentially-selected PlanckXXVII sources. However, we also find significant differences between the PHz HELP and PlanckXXVII galaxy colours. The two-sample KS test rejects (at 16.7$\,\sigma$) the null hypothesis that the distribution of $S_{350}/S_{250}$ colours for the PHz HELP sources is drawn from the same distribution as the PlanckXXVII sources. Similarly, for the distribution of $S_{500}/S_{350}$ colours, the KS test rejects the null hypothesis at 12.6$\,\sigma$. Furthermore, the peak in PHz HELP $S_{500}/S_{350}$ colours better aligns with the peak for all PHz sources from Planck colours.

To further illustrate the differences in galaxy colours between the PHz HELP, PlanckXXVII and random Lockman-SWIRE positions, Fig. 8 shows 2D histograms of $S_{350}/S_{250}$ versus $S_{500}/S_{350}$. The PHz HELP sources have a clear tail of red galaxies extending to the right, which is not seen in the PlanckXXVII sources or the random sources. However, the PlanckXXVII colour distribution appears shifted upwards when compared to the other two distributions. Based on the distribution of random galaxy colours, we set the thresholds $S_{350}/S_{250}\,{=}\,1.2$ and $S_{500}/S_{350}\,{=}\,0.8$ to select the reddest galaxies, falling in the upper-right region of Fig. 8. Overlaying these thresholds on the 2D colour histograms of the PHz HELP and PlanckXXVII sources illustrates the differences between the colour distributions; when compared with the PlanckXXVII sources, significantly more PHz HELP sources lie above the $S_{350}/S_{250}\,{=}\,1.2$ threshold and noticeably fewer lie above the $S_{500}/S_{350}\,{=}\,0.8$ threshold. Density histograms of red galaxies that fall above these two cuts confirm these trends (see Fig. 9). The two-sample KS test rejects the null hypothesis that the distribution of red PHz HELP sources is the same as the distribution the red PlanckXXVII sources at 8.0$\,\sigma$ for the $S_{350}/S_{250}\,{>}\,1.2$ sources and 5.3$\,\sigma$ for the $S_{500}/S_{350}\,{>}\,0.8$ sources.

Based on the 2D colour distributions, we can also infer that fewer PHz HELP galaxies peak in flux density at 250 $\mu$m as compared with both random galaxies and PlanckXXVII galaxies. Quantitatively, of random galaxies detected in the Lockman-SWIRE field, a total of 38 per cent peak at 250 $\mu$m, and comparatively, 37 per cent of PlanckXXVII galaxies peak at 250 $\mu$m. For the PHz HELP sources, on the other hand, we find that only 21 per cent of galaxies peak at 250 $\mu$m. Instead, most PHz HELP galaxies (66 per cent) peak at 350 $\mu$m, whereas only 41 per cent of galaxies in the PlanckXXVII sources peak at 350 $\mu$m. In this way, the PHz HELP sources typically consist of redder galaxies than the PlanckXXVII sources do. This suggests that the distribution of galaxy redshifts in PHz HELP sources may differ from the distribution of redshifts in PlanckXXVII sources, a possibility that is further explored in Section 4.5; however, it is difficult to draw conclusions about galaxy redshifts with SPIRE photometry alone. It is also worth noting that PlanckXXVII reports that none of the brightest galaxies within the PlanckXXVII sources peak at 250 $\mu$m (i.e., within each PlanckXXVII source, the brightest galaxy across all three bands is brightest at 350/500 $\mu$m), while we find that a significant fraction of the brightest PlanckXXVII galaxies peak at 250 $\mu$m (56 per cent), and a smaller, but still significant, fraction of the brightest PHz HELP galaxies peak at 250 $\mu$m (37 per cent).

Owing to the fact that the rest-frame far-IR emission observed by SPIRE probes a relatively simple, thermal portion of a galaxy’s spectral energy distribution (SED), SPIRE colours can be used to provide crude estimates of photometric redshifts, luminosities and therefore SFRs, after applying various simplifying assumptions (e.g., Blain et al., 2003; Draine, 2003). Assuming that any AGN contribution is sub-dominant, the SPIRE bands are sensitive to thermal dust emission generated by star formation, characterized by a modified blackbody distribution with a single dust temperature. However, due to redshifting, we observe a dust temperature that is decreased by a factor of $1\,{+}\,z$, preventing us from determining both the intrinsic dust temperature and the redshift simultaneously. Nonetheless, if we assume a typical dust temperature motivated by studies of similar SPIRE-detected galaxies, we can obtain rough redshift estimates and SFRs, which should at least be useful for comparitive studies.

Next, we calculate far-IR luminosities ($L_{\rm FIR}$) by integrating our fits to Eq. 2 between 8 and 1000 $\mu$m (in the rest-frame), and convert these to SFRs using a conversion factor of $0.95\,{\times}\,10^{-10}\,{\rm M}_{\odot}\,{\rm yr}^{-1}\,{\rm L}_{\odot}^{-1}$, which assumes a Chabrier initial mass function (Chabrier, 2003). The distributions of $L_{\rm FIR}$ and SFR are shown in Fig. 10 (right panel). We again find similar trends when the assumed dust temperature is changed by $\pm 5$ K, giving an uncertainty of $\delta L_{\rm FIR}\approx 5\times 10^{12}\,{\rm L}_{\odot}$ and $\delta\mathrm{SFR}\approx 5\times 10^{2}\,{\rm M}_{\odot}\,{\rm yr}^{-1}$.

In agreement with PlanckXXVII, we find that the redshift distribution of galaxies in the PHz HELP sources peaks at $z\,{\approx}\,2$. The PHz HELP galaxy redshifts are largely comparable to the PlanckXXVII galaxy redshifts and both have a positive tail of high-$z$ galaxies when compared with galaxies found in random positions. The KS test confidently rejects the null hypothesis (at ${>}\,15\,\sigma$) that the PHz HELP redshifts are drawn from the same distribution as the random galaxy redshifts.

It remains unknown whether the PHz sources are predominantly protoclusters of physically-associated galaxies at the same redshift or line-of-sight alignments of multiple groups of galaxies (e.g., Negrello et al., 2017; Gouin et al., 2022), as discussed further in Sect. 6. By comparing the SPIRE photometric redshifts of multiple galaxies detected within the same PHz source, we can attempt to quantify the presence of line-of-sight alignments within the PHz HELP sources, but we find that the SPIRE photometric redshifts are too imprecise ($\delta z\approx 1$) to distinguish between overdensities at a single redshift from random alignments along the line-of-sight.

Both the PHz HELP sources and PlanckXXVII sources generally have higher far-IR luminosities/SFRs than random sources (the KS test puts the significance of the differences at ${>}\,15\,\sigma$). Furthermore, the far-IR luminosity/SFR distribution of the PlanckXXVII sources is shifted to the right of the distribution of the PHz HELP sources; the difference between the PHz HELP and PlanckXXVII sources is significant at ${>}\,12\,\sigma$ by the KS test. The larger far-IR luminosities/SFRs of galaxies within the PlanckXXVII sources agrees with the larger flux densities of the PlanckXXVII sources, which is seen in the galaxy number counts (Fig. 4). Nonetheless, in the more representative PHz HELP sources, we still find an excess of high far-IR/SFR galaxies when compared with random positions.

We next investigate the spatial properties of the PHz HELP sources through a stacking analysis. Because we are primarily interested in the spatial distributions of the overdense sources, we remove the sources that consist of a single bright galaxy (namely the Planck GEMS from the PlanckXXVII sample, and the bright gravitational lens G160 from the PHz HELP sample). We create mean image stacks by averaging the 250-, 350- and 500-$\mu$m maps centred on the position of the brightest 350-$\mu$m galaxy in each IN region.²²2We explored other stacking choices as well, including stacking on the flux-weighted centroid positions, and found similar trends. Before computing the average, we set the background of each image to 0 Jy beam^-1 by subtracting the mean intensity value in pixels that are ${>}\,9$ arcmin from the centre. In Fig. 11 we show the resulting stacks at 350 $\mu$m for the PHz HELP sources, PlanckXXVII sources and random positions in the Lockman-SWIRE field centred on the brightest galaxy in each 4.34 arcmin radius circular region. It is clear from the 350-$\mu$m stacks that the PHz HELP sources are significantly less bright than the PlanckXXVII sources. Nonetheless, similar to the PlanckXXVII stack, the PHz HELP stack has clustering extending out to about 5 arcmin from the centre that is not seen around random galaxies in the Lockman-SWIRE field. This confirms that average PHz sources display clustering, just as the preferentially-selected bright PlanckXXVII sources do.

To quantify the extended nature of the PHz HELP and PlanckXXVII stacks, we plot the radial surface brightness profiles of the PHz HELP, PlanckXXVII and random position stacks (see Fig. 12). Confirming the visual conclusions drawn from Fig. 11, the PHz HELP stacks are less bright than the PlanckXXVII stacks, but still clearly extended when compared with random sources across 250, 350 and 500 $\mu$m. In comparison with the random sources, the PHz HELP sources display the strongest clustering at 500 $\mu$m and the weakest at 250 $\mu$m. This is likely caused by the PHz selection process, which relied primarily on detecting excess emission at 500 $\mu$m. To directly compare the PHz HELP and PlanckXXVII radial profiles, we arbitrarily normalize both profiles to 1.0 Jy beam^-1 in the 1.0–2.0 arcmin bin (Fig. 13). The PHz HELP and PlanckXXVII radial profiles closely match outside the central 1.0 arcmin bin; however, the PHz HELP sources are brighter within this central bin, suggesting that the PHz HELP sources are (relatively) more centrally-concentrated.

To confirm that the PHz HELP sources are genuinely more centrally-concentrated than the PlanckXXVII sources, we compare the ratios of the flux density contained within 1.0 arcmin of the brightest galaxy to the total flux density of the IN region ($S_{\mathrm{1\,arcmin}}/S_{\mathrm{IN\,region}}$; see Fig. 14). For the purpose of this plot, we adopt a 4.34 arcmin radius circle for the IN region of each source so that differences in $S_{\mathrm{1\,arcmin}}/S_{\mathrm{IN\,region}}$ are not caused by variations in the IN region sizes. As alluded to by the normalized radial profiles, a greater fraction of the PHz HELP source flux densities come from the 1.0 arcmin region around the brightest galaxy when compared to the PlanckXXVII sources. According to the KS test, the differences between the PHz HELP and PlanckXXVII sources is significant at ${>}\,3\,\sigma$ across the 250-, 350- and 500-$\mu$m distributions. We thus conclude that the PHz HELP sources are less bright, but more centrally-concentrated than the PlanckXXVII sources.

Because the photometric redshift information offered by the three far-infrared SPIRE bands is too uncertain to identify redshift clustering about a single target, we turn to the spatial overdensity data to identify the most promising protocluster candidates. We then statistically search for clustering within the combined photometric redshift distribution of the best protocluster candidates.

In PlanckXXVII, Herschel-SPIRE follow-up observations of 228 Planck-selected sources led to the discovery of 11 of the brightest known gravitationally-lensed submm galaxies (‘Planck’s dusty GEMS’; Cañameras et al. 2015). Due to the exceptional brightness of the GEMS and the relative paucity of other SPIRE galaxies in these Planck sources, identification of such bright gravitational lenses turned out not to be very challenging. PlanckXXVII classified sources as either overdensities or candidate gravitationally-lensed galaxies using two quantities: the flux density of the brightest galaxy in the source (at any wavelength); and the wavelength at which the brightest galaxy peaks. Specifically, any sources containing a galaxy with flux density ${>}\,400$ mJy were identified as candidate lens systems and all other sources whose brightest galaxy peaks at 350 or 500 $\mu$m were identified as overdensities. This classification scheme identified seven of the PlanckXXVII sources as candidate lensed galaxies, and the authors selected an additional five candidate lensed galaxies by visual inspection (all of which have peak flux ${\gtrsim}\,300$ mJy and appear as the only bright galaxy in the source). These 12 sources consisted of one previously-discovered strong gravitationally-lensed galaxy and 11 new lensing candidates, which have since been confirmed as bright gravitationally-lensed galaxies making up the GEMS sample (Cañameras et al., 2015).

Only one PHz HELP source, G160, contains a single exceptionally bright SPIRE galaxy, like the GEMs, and turns out to be a previously-discovered strong gravitationally-lensed galaxy (Geach et al., 2015; Harrington et al., 2016). However, not all strong gravitationally-lensed galaxies are as bright as Planck’s dusty GEMS; the PHz sources could contain high-$z$ strong gravitationally-lensed galaxies of somewhat lower brightness, comparable to other known gravitational lenses (e.g., Negrello et al., 2010; Weiß et al., 2013). Furthermore, unlike the GEMS, overdense PHz sources may contain high-$z$ strong gravitationally-lensed galaxies, making them difficult to identify visually. Studies of Herschel wide-field surveys have used a simple detection threshold of 100 mJy at 500 $\mu$m to identify candidate lensed galaxies. Applying this criterion to the galaxies within the PHz HELP sources, we find a total of just four sufficiently bright galaxies. To search for more candidates, here we use a machine-learning (ML) approach to select gravitationally-lensed galaxies in the PHz HELP sources.

Now that we have investigated the resolved rest-frame far-IR properties of a representative subsample of PHz sources, we can discuss the overall statistical properties of the PHz catalogue, including some corrections to the conclusions drawn with the biased subsample analysed in PlanckXXVII.

In the preferentially-selected PlanckXXVII subsample, which had a bias towards high S/N sources with compact/regular shapes and included sources from the PCCS, 12 sources turned out to be ultra-bright strong gravitationally-lensed galaxies (5 per cent of their subsample). Based on our more representative sample, we now see that the fraction of bright strong gravitational lenses (peak flux ${\gtrsim}\,300$ mJy) in the PHz catalogue has been overestimated; we find only one such galaxy (a rate of 0.5 per cent), namely G160. The larger fraction of bright strong gravitationally-lensed galaxies found in the PlanckXXVII sample suggests that such bright gravitational lenses are more easily found in the Planck compact source catalogues (see Trombetti et al. 2021), rather than the PHz catalogue.

On the other hand, we confirm that most PHz sources correspond to overdensities in SPIRE maps. This is particular illustrated through the stacking analysis, in which average stacks of the PHz HELP sources, similarly to the PlanckXXVII sources, show significant extended emission out to about 5 arcmin, well beyond the SPIRE beam at all three observed wavelengths (with the most substantial differences at $500$ $\mu$m). Additionally, we found a total of 40 sources having ${>}\,3\,\sigma$ overdensities of SPIRE galaxies, corresponding to 21 per cent of the PHz HELP sources; this is comparable to the 25 per cent of ${>}\,3\,\sigma$ overdensities found in the PlanckXXVII sources and includes a larger fraction of sources that are overdense at 350 and 500 $\mu$m.

Despite the similarity in these results, we find significant differences between the PHz HELP and PlanckXXVII sources when we examine the SEDs of the galaxies found across the two subsamples. For our PHz HELP sources, we measure lower overdensities at 250 and 500 $\mu$m when compared to the PlanckXXVII sources, with a slightly greater positive tail at 350 $\mu$m (see Fig. 5). This is related to the differences in colour seen between the galaxies in these two subsamples, for which we find larger values of $S_{350}/S_{250}$ in the PHz HELP galaxies and larger values of $S_{500}/S_{350}$ in the PlanckXXVII galaxies. This may indicate that galaxies found in the PlanckXXVII sources are at systematically higher redshift or contain systematically cooler dust temperatures, although we cannot provide any definitive conclusions because of the large uncertainties associated with the SPIRE photometric redshifts. Indeed, the distributions of photometric redshifts appear nearly identical for the two subsamples (see Fig. 10).

We know that the PHz sources are places on the sky corresponding to peaks in the CIB that are red. Hence we know that these are special directions, containing galaxies with red colours at these wavelengths. However, we would like to know whether the observed projected overdensities are predominantly protoclusters of galaxies at the same redshift or line-of-sight alignments of two or more groups of galaxies. Unfortunately, our current photometric redshifts are too uncertain to distinguish between these two situations. Instead, to make progress, we turn to a comparison with simulations. Specifically, we repeated the PHz selection algorithm on the Simulated Infrared Dusty Extragalactic Sky (SIDES; Béthermin et al., 2017) simulation data. Briefly, SIDES simulated 2 deg² of the extragalactic sky at the wavelengths observed by the SPIRE instrument (i.e., 250, 350 and 500 $\mu$m), as well as 850 $\mu$m, using a two-scenario model of star formation, where galaxies are either actively forming stars or quenched. Clustering was included by populating dark matter halos from dark matter-only simulations using an abundance matching technique (see Béthermin et al. 2017 for more details). We smoothed the SIDES catalogue to the 5 arcmin resolution of Planck, and reproduced the 500-$\mu$m excess selection criteria used to generate the PHz catalogue (see Planck Collaboration XXXIX 2016). Upon investigating the 30 brightest peaks in the resulting excess map, we did not find any strong evidence of large-scale clustering or line-of-sight alignments; instead, the selection criteria predominantly picked out individual galaxies bright at 500 $\mu$m, with only mild clustering. This is likely caused by the small size of the SIDES simulation, preventing it from reproducing the rare objects that make up the PHz catalogue; the 2151 PHz sources were found across 10,700 deg², meaning that their density is about 0.2 deg^-2, and thus a simulation at least as large as $10\,{\rm deg}^{2}$ would be needed to find comparably rare objects.

Given the limitations of the SIDES simulation, we also perform tests using an analytical formalism as an alternative way to investigate potential line-of-sight alignments. Specifically, we employ the analytic model developed by Negrello et al. (2017), which gives the luminosity function (LF) of clumps (either an isolated galaxy or a physically bound group of galaxies) at 545 GHz based on galaxy evolution models and the infrared LFs of individual galaxies; this model is found to match observations well. Following the simple Monte Carlo simulations done by Negrello et al. (2017), we draw clumps from the LF of their model, and distribute them randomly over 5,000 deg² of a mock sky — we expect this area to be large enough to generate the rare peaks in the CIB found by Planck, and we note that the 40 protocluster candidates found here were drawn from 1,270 deg².

We smooth the mock sky with a Gaussian corresponding to the Planck beam of 5 arcmin. From this smoothed map we select peaks with S/N$\,{>}\,$5 to mimic the PHz selection, where the noise is the rms of the map, that is, purely the confusion noise. Finally, from this sample, we select sources with 3 $\sigma$ overdensities of individual galaxies brighter than ${\simeq}\,32$ mJy at 545 GHz, corresponding to our selection of protocluster candidates as described in Sect. 5.1. The 500-$\mu$m excess selection used in Planck Collaboration XXXIX (2016) to eliminate low-redshift sources and other potential contamination is not necessary here because the LF adopted for the test only includes clumps with redshifts between 2 and 4.

From this test we find about 220 simulated protocluster candidates (or line-of-sight alignments), selected in the same way as our protocluster candidate catalogue. This roughly matches the number of candidates projected for the given simulated area. Upon investigating these simulated protocluster candidates, we find that 96 per cent of them contain two or more clumps aligned along the line-of-sight that are not physically connected to one another. Similarly, 81 per cent are found to have more than three clumps aligned along the line-of-sight by chance. The average number of physically unbound clumps along the line-of-sight is about four, while that for random directions is about one.

This suggests that most of the PHz sources, and our protocluster candidates selected from that catalogue, may be line-of-sight combinations of star-forming galaxies at different redshifts, rather than single clumps of physically associated galaxies. This is similar to the conclusions drawn in Negrello et al. (2017) and Gouin et al. (2022), which found that the Planck-derived SFRs for the PHz sources can be explained only by a dominant contribution from background/foreground interlopers, rather than by single haloes with very high SFRs. This conclusion, that there are foreground and background interlopers, is of course not very surprising. The more interesting question is whether there is a dominant star-forming protocluster along these special lines of sight. From the same simulations we find that the ratios of flux densities between the brightest and second brightest clumps along the selected lines of sight are much higher than that in random directions, with the brightest clump being a factor of approximately 2 (3) times brighter than the second brightest clump in 76 (60) per cent of the cases. The average flux density ratio between the brightest and second brightest clumps is about 4.4, in comparison with 1.5 for random directions. These results suggest that while the PHz sources and our protocluster candidates do tend to have line-of-sight interlopers, they still tend to be dominated by a significantly bright clump of galaxies at a single redshift, motivating continued follow-up studies with protocluster science in mind.

To summarize the conclusions of this section, we find that the PHz catalogue consists primarily of angular overdensities of $z\,{\sim}\,2$–3 galaxies, with only a very small fraction of bright strong gravitational lenses. Simulations suggest that there are multiple structures along each line of sight, but that they are often dominated by a single structure. In order to determine whether these spatial overdensities contain genuine protoclusters, further spectroscopic follow-up observations are needed. Our set of 40 ${>}\,3\,\sigma$ overdensities provides an optimal list for such studies, and should help us better understand how efficient CMB surveys are at finding the population of star-forming protoclusters around redshifts of 2.

Current models of galaxy formation underestimate the SFRs observed in high-$z$ protoclusters (Lim et al., 2021), including the PHz source G237 (Polletta et al., 2021). Even though simulations suggest that chance line-of-sight alignments are quite common among PHz sources, if each line-of-sight is at least dominated by one bright clump (as suggested above), then our Herschel-SPIRE observations may still be useful for constraining star-formation in typical protoclusters at $z\,{\sim}\,2$. In particular, in Section 4.6 we stacked cutouts of the PHz HELP sources centered on the brightest galaxy detected at 350 $\mu$m in each source. We expect that these brightest galaxies are part of the dominant clumps, and therefore our stacks may average out the fainter clumps along each line-of-sight.

Currently, statistical samples of protoclusters at higher redshift do not exist, but at $z\,{\approx}\,1$ Alberts et al. (2021) studied 53 near-IR-selected galaxy clusters observed by Herschel-SPIRE at 250, 350 and 500 $\mu$m. In Fig. 19 we compare our stacked $z\,{\sim}\,2$ surface brightness profiles (as well as those obtained from PlanckXXVII for reference) to the $z\,{\approx}\,1$ stack from Alberts et al. (2021) at 250, 350 and 500 $\mu$m. Here we have followed Alberts et al. (2021), adopting the same bin sizes and arbitrarily normalizing the profiles to a value of 0.5 at 27.5 arcsec.

Interestingly, we find that the average far-infrared extent of the PHz HELP sources is much larger than that of the average $z\,{\approx}\,1$ cluster; it is worth noting that angular diameter distances only differ by about 5 per cent between redshift 1 and 2, so the observed differences in size are physical. We know that the average SFR density of the Universe peaked around $z\,{\approx}\,2$ (e.g., Behroozi et al., 2013; Gruppioni et al., 2013; Madau & Dickinson, 2014; Chiang et al., 2017; Koprowski et al., 2017), and therefore we might expect $z\,{\sim}\,2$ protoclusters to be more extended in star-formation as well. However, this is speculation, and the exact mechanisms driving any increase in SFR and size would need to be investigated in more detail.

Similarly, it is interesting to note that the PHz HELP sources are more centrally concentrated than the PlanckXXVII sources. Since our stacking procedure relied on centering each image on the brightest 350 $\mu$m-detected galaxy, the centre of each stack is effectively a point source with a flux density given by the average flux density of the brightest PHz HELP/PlanckXXVII galaxies. We found that across the PHz HELP sources, the brightest central galaxy is relatively brighter than its surroundings when compared with the PlanckXXVII sources. This could come from the bias towards brighter sources in selecting the PlanckXXVII sample, which not only picked brighter central galaxies, but also sources with bright galaxies distributed throughout Planck’s 5 arcmin beam. With a more representative sample of PHz targets, we find that typical sources are, in fact, more likely to contain bright central galaxies relative to the surrounding faint galaxies, as opposed to collections of relatively bright galaxies. This may indicate the presence of brightest cluster galaxies (BCGs) forming in protocluster cores, although further work is needed to test this speculation.