An optimal ALMA image of the Hubble Ultra Deep Field in the era of JWST: obscured star formation and the cosmic far-infrared background

We focus on ALMA Band 6, which spans a wavelength (frequency) range of 1.1–1.4 mm (210–270 GHz). This band currently contains the most extensive ALMA observations of the HUDF, and so this is where we expect to produce the deepest archival map. To begin, we queried the ALMA archive¹¹1https://almascience.nrao.edu/asax for all Band-6 observations centred at 03^h32^m39.0^s $-$27^∘47^′29.1^′′ and overlapping within a radius of 1.5 arcmin. There are a total of 12 unique programmes that satisfy this criterion, and we selected all of these for our combined map. For some programmes, only a fraction of the total time was spent on the HUDF (e.g. for individual pointings of targeted galaxies or for large surveys extending beyond the HUDF), and for these cases we only selected the data that overlap with our region of interest. For more details see Table 1, where we summarize all of the data used here to produce the combined map. We also downloaded the raw data from the ALMA archive centred at the same position, but extending out to 3.5 arcmin, which we used to make a shallower but larger map; there are a total of seven additional programmes satisfying this criterion, and these are listed in Table 4.

For each of the observations, we downloaded the raw $uv$ data and calibrated it using the provided ScriptForPI and the CASA (McMullin et al., 2007) version used by the observatory at the time of the observations. We split the science targets from the calibration targets, then time-averaged the data by 30 s and averaged the frequency channels by a factor of 4 in order to reduce the volume of data. These tasks were carried out using the Canadian Advanced Network for Astronomy Research (CANFAR) platform (Major et al., 2019), which provides easy access to all versions of CASA, as well as ample storage space and memory.

In addition to downloading the $uv$ data for each programme, we also obtained the imaging products made available by the observatory. These data are used for understanding the transmission function of the combined data (see Section 4) and for testing an alternative data combination method (see Appendix C). For every tuning of every target given in Table 1, the observatory provides a single image made using the multi-frequency synthesis (MFS) mode, where visibilities in each channel are mapped to a single $uv$-plane, and therefore represent the mean value of the sky at a characteristic frequency (usually the central frequency of the channels) weighted by a spectral function. For programmes carried out in earlier ALMA cycles, we use the Additional Representative Images for Legacy (ARI-L; Massardi et al., 2021), which are reduced in a way more similar to later cycles. These images are primary-beam-corrected, but the primary-beam image is also available for download. The final number of HUDF images downloaded from the ALMA archive is 61, along with their corresponding 61 primary-beam images.

Our main goal is to produce the deepest possible map of the central HUDF region, which essentially corresponds to the footprint of ASPECS (see González-López et al., 2020). We do so by combining all of the available data within this region. Our secondary goal is to produce the deepest possible map of the wider HUDF region, corresponding essentially to the ASAGAO footprint (see Hatsukade et al., 2018). To do this, we carried out the following procedure.

All of the time-averaged and frequency-averaged $uv$ data were imaged using the standard CASA task tclean. For the data presented in Table 1, the pixel size was set to 0.2 arcsec. We ran tclean in MFS mode, which, as discussed above, scales $uv$ points observed at different frequencies to a single characteristic frequency, which in our case is the average frequency of the data being combined, 243 GHz (or 1230 $\mu$m). In principle spectral variations in the observed sources can lead to imaging artefacts; however, for frequency ranges of about 30 per cent these artefacts are expected to be small (Sault & Conway, 1999). The frequencies that are combined to produce our map range from about 211 to 275 GHz, or roughly 30 per cent, thus we do not attempt to correct for potential artefacts. We also chose natural weighting (each $uv$ point is weighted by its instrumental noise), but with a 250 k$\lambda$ $uv$ taper. We set the cutoff of the map to 0.2 times the primary beam, which is where the footprint of our map roughly matches the footprint of the ASPECS map. We also produced a map with no $uv$ taper, setting the pixel size to 0.18 arcsec (because the synthesized beam is somewhat smaller). Lastly, we cleaned the image down to 20 $\mu$Jy beam^-1, similar to the cleaning level chosen by ASPECS. The final synthesized beamsize of the $uv$-tapered map is $1.49\,{\rm arcsec}\,{\times}\,1.07\,{\rm arcsec}$, while for the map with no $uv$ taper the synthesized beamsize is $1.32\,{\rm arcsec}\,{\times}\,0.92\,{\rm arcsec}$. The pixel rms in the tapered image (which we use for further analysis in this paper) has a minimum value of 4.6 $\mu$Jy beam^-1 in the deepest region, while the rms is less than 9.0 $\mu$Jy beam^-1 over an area of 1.5 arcmin². Lastly, we tested tapers of 150 k$\lambda$ and 50 k$\lambda$, but we found that the trade-off in sensitivity for the larger synthesized beams led to fewer detected sources. We note that we found similar source flux densities across the different tapers, meaning that at 250 k$\lambda$ we are not resolving out much flux. On the other hand, the map with no $uv$ taper resolves a larger number of sources, making photometry more challenging.

For the shallower, larger map, we exclude all of the observations targeting the deep central region (this includes all of the ASPECS data, the data from Dunlop et al. 2017, and two deep pointings towards the eastern corner of the HUDF); this is to ensure that the synthesized beam of the larger map is not dominated by data from the central region (see Table 4 for details). This map is therefore mostly a combination of the ASAGAO survey data and the GOODS-ALMA survey data (both high and low resolution). However, we include all of the ASAGAO follow-up observations even if they lie within the deep central region (programme ID 2018.1.00567.S), since they roughly match the array configurations and frequencies of the ASAGAO survey and the GOODS-ALMA survey. We set the pixel size to 0.06 arcsec for the map with no $uv$ tapering, and 0.12 arcsec for the map with 250 k$\lambda$ $uv$ tapering, and set the cutoff of the map to 0.4 times the primary beam, which produced a map with a footprint roughly similar to the ASAGAO footprint. All other tclean parameters were kept the same as described above, except for the cleaning threshold, which was increased to 100 $\mu$Jy beam^-1 (similar to the level chosen by ASAGAO). The final synthesized beamsize of the $uv$-tapered map is $0.87\,{\rm arcsec}\,{\times}\,0.64\,{\rm arcsec}$, while for the map with no $uv$ taper the synthesized beamsize is $0.36\,{\rm arcsec}\,{\times}\,0.27\,{\rm arcsec}$. The pixel rms in the tapered image (which we also use for further analysis in this paper) has a typical value of 60 $\mu$Jy beam^-1 outside of individual pointings, and 20 $\mu$Jy beam^-1 within individual pointings. Similar to the deep map, we tested tapers of 150 k$\lambda$ and 50 k$\lambda$, but we found that the 250 k$\lambda$ taper recovers the most sources without resolving out too much flux.

Lastly, we produced maps with uniform noise properties by excluding programmes with individual pointings targeting known galaxies, and instead focusing on large survey programmes. Specifically, for our deep central mosaic we combined data from the programmes 2012.1.00173.S, 2015.1.00098.S, 2015.1.00543.S, 2016.1.00324.L and 2017.1.00755.S, while for the shallower, larger map we only used data from the programmes 2015.1.00098.S, 2015.1.00543.S and 2017.1.00755.S (see Tables 1 and 4, respectively). We used the same tclean parameters as above, generating versions with no $uv$ tapering and with 250 k$\lambda$ $uv$ tapering. The final synthesized beamsizes were effectively unchanged compared to the maps that included individual pointings.

For the remaining analysis in this paper we focus on the tapered maps with all of the available data combined. The additional maps, with no tapering and uniform noise properties, were used to check the reliability of our flux densities, and we found no systematic differences between the measurements.

To calculate the corresponding noise maps for the combined images, we first created a mask using catalogues of previously-detected galaxies. In particular we took the catalogues from Dunlop et al. (2017), ASAGAO (Hatsukade et al., 2018), ASPECS (González-López et al., 2020) and GOODS-ALMA (Gómez-Guijarro et al., 2022). For each galaxy we masked a region 3 times larger than the beamsize, and then multiplied this by the signal map. We next made cutouts of 5 times the beamsize around each pixel and calculated the rms. The resulting noise map was then smoothed with a Gaussian kernel with the same size as the cutout regions to remove artefacts.

All of these data products are made publicly available, including the primary-beam-corrected maps, the noise maps, the primary beam maps and the synthesized beam maps.²²2https://doi.org/10.5683/SP3/YWBVWH The final signal map, noise map and signal-to-noise ratio (S/N) map for the deep central mosaic with a 250 k$\lambda$ taper, including individual pointings, is shown in Fig. 1, while the same maps for the larger, shallower mosaic are shown in Fig. 13. The deep central mosaic covers an area of 4.2 arcmin² out to our primary beam threshold of 0.2, while the large shallow mosaic covers an area of 25.4 arcmin² out to our primary beam threshold of 0.4.

The HUDF has been the target of extensive multiwavelength observations that we can use to complement our 1.23-mm image. The Cosmic Assembly Near-IR Deep Extragalactic Legacy Survey (CANDELS) catalogue of galaxies in the GOODS-S field³³3https://archive.stsci.edu/prepds/candels (Guo et al., 2013) contains a summary of much of the deep optical-to-near-IR imaging in the HUDF, and it is expected that most of the ALMA-detected galaxies should have a counterpart in this catalogue. Detections in the CANDELS catalogue were obtained from HST’s WFC3 instrument in the F160W band, so it is a 1.5-$\mu$m-selected catalogue. The entire catalogue covers a much larger region than our deep map, but the deepest region of the catalogue covers roughly the same area.

The first public data release from the the JWST Advanced Deep Extragalactic Survey (JADES) NIRCam survey of the HUDF is also now available⁴⁴4https://archive.stsci.edu/hlsp/jades (Rieke & the JADES Collaboration, 2023), in addition to the JWST Extragalactic Medium-band Survey (JEMS, Williams et al., 2023). While the JADES survey is still ongoing and will eventually be deeper, it already contains considerably more galaxies per unit area than the CANDELS catalogue. JADES images range from 0.9 $\mu$m to 4.4 $\mu$m, and we chose to select our catalogue at $3.5\,\mu$m, as a compromise between the longest possible wavelength and the depth. Our entire deep map is covered by the JADES survey at approximately uniform sensitivity. Throughout this paper we primarily make use of the JADES imaging, with the CANDELS catalogue used to demonstrate new improvements thanks to JWST.

The central deep map shown in Fig. 1 contains all existing ALMA Band-6 data in this region and thus should be the best map currently achievable; for reference, the deepest part goes down to 4.6 $\mu$Jy beam^-1, which is 50 per cent deeper than any previous map. Over an area of 1.5 arcmin² our new map has an rms below 9.0 $\mu$Jy beam^-1, which is 5 per cent better than the deepest previously available map of this size (from ASPECS), and it has a noise level below 35 $\mu$Jy beam^-1 over the full $4.2\,{\rm arcmin}^{2}$. It is therefore of interest to see if any new galaxies are detected in this new map.

We ran the simple peak-finding algorithm find_peaks, available in the photutils python module, on both the S/N map contained within the region of interest, and on the negative of the same S/N map. One approach to setting the detection threshold is to use the most significant negative peak to set the level above which we might expect all positive peaks to be real galaxies; for reference, the most negative peak was found to have a S/N of $-4.1$, and there are 35 positive peaks with a S/N higher than $+4.1$. However, we can also lower the threshold to include more real sources at the expense of being less confident about the reality of each one.

As an alternative means of setting the threshold, we looked at how the ratio of the total number of positive to negative peaks greater than a given S/N thresholds varies. We can define the ‘purity’ to be $1\,{-}\,{\cal N}_{\rm neg}/{\cal N}_{\rm pos}$, such that a value of 0 means there are as many negative as positive peaks and 1 is reached when there are no more significant negative peaks. In Fig. 2 we plot this purity as a function of S/N threshold. The choice of a threshold is of course a trade off (a smaller threshold will result in more false positives), and to include more candidate sources (with the understanding that not all may be real) we choose a purity of 0.7, which corresponds to a peak ${\rm S/N}\,{=}\,3.6$. There are 45 sources above this threshold. The catalogue from ASPECS used a ‘fidelity’ threshold (the differential ratio of galaxies in S/N bins) that resulted in a least significant source with ${\rm S/N}\,{=}\,3.3$, while Dunlop et al. (2017) used a S/N threshold of 3.5 but removed sources with no HST counterpart. Our choice of threshold is therefore similar to (although slightly more conservative than) what was used in previous studies, and we can also use our JADES catalogue to investigate possible false positives.

In Fig. 3 we show the distribution of pixel values (in units of $\mu$Jy beam^-1) in our primary beam-uncorrected map. Here, instead of using the primary beam provided by CASA, we scale our own noise map by the minimum noise value to calculate the primary beam map, since this is the map used to search for sources. We can see that the distribution appears quite Gaussian at the negative end, while real sources heavily skew the positive distribution. Our source extraction relies on the Gaussianity of the noise, and so we fit a Gaussian to the pixel distribution, after masking all positive pixels brighter than the most negative pixel (here $-19\,\mu$Jy beam^-1). We find a mean value of $-0.08\,\mu$Jy beam^-1 (very close to zero, as expected) and a width of $4.60\,\mu$Jy beam^-1. Our expectation is that the width should be equal to the minimum noise value used to create our primary beam estimate, which here is $4.58\,\mu$Jy beam^-1, so we find good agreement. Lastly, we calculate the relative residual between the Gaussian fit and the actual pixel distribution (defined as data ${/}\,$model $-$ 1), and show this in the bottom panel of Fig. 3. We find that the Gaussian model deviates from the actual data by no more than 20 per cent for pixel values less than $20\,\mu$Jy beam^-1, after which real sources dominate and the distribution is no longer Gaussian. Of course, we still expect real sources to contribute to the map below $20\,\mu$Jy beam^-1 in this ‘P(D)’ analysis (see e.g. Vernstrom et al., 2014, and references therein), which could be the cause of the deviation.

As a check, we compare the 1-mm flux densities extracted here to the flux densities given in the four previously published surveys. For simplicity, in this comparison we simply match published sources to our new catalogue using a search radius of 1 arcsec. The mean frequency of the map from Dunlop et al. (2017) is 221 GHz (see Table 1), so we correct their flux densities to the mean frequency of our map (243 GHz) assuming a modified blackbody SED with spectral index $\beta\,{=}\,1.8$, a dust temperature of 25 K, and using each galaxy’s spectroscopic/photometric redshift when available, otherwise a redshift of 1.5. Similarly, the mean frequency of the map from GOODS-ALMA (Gómez-Guijarro et al., 2022) is 265 GHz, so we follow the same procedure and apply a correction factor of 0.8. The mean frequency of the remaining maps are effectively identical to ours, so we do not apply any further corrections.

The flux densities of matched peaks are shown in Fig. 4, and the cross-matched IDs are given in Table 3. We find generally good agreement with all of the previously-published flux densities after applying the above corrections. The flux densities from Dunlop et al. (2017) appear to have more scatter compared to later studies; however, those data were taken in ALMA’s early observing cycles when calibration was more uncertain. We also see that the brightest source in our map (ALMA ID 1) is somewhat fainter than what was reported in the ASPECS and GOODS-ALMA surveys, but this source is right at the edge of the ASPECS primary beam (and thus our primary beam as well), so we expect that there may be additional systematic uncertainties that are not being taken into account for this source.

In Fig. 5 we show our detected sources compared to those found by Dunlop et al. (2017), the ASAGAO survey (Hatsukade et al., 2018) and the ASPECS survey (González-López et al., 2020), as well as a few sources from the wider GOODS-ALMA survey (Gómez-Guijarro et al., 2022). It appears that most of the detections in our combined map coincide with published sources, but there are 13 new DSFGs. Of these 13 new DSFGs, six have a peak S/N$\,{>}\,4.1$, meaning that they have been detected at a purity level of 100 per cent. There are also a few published detections that do not appear here; typically, these are low-significance sources that could have been positive noise excursions.

In Appendix A we perform the same source extraction procedure on the larger but shallower map covering the ASAGAO footprint (excluding the central deep region, where we have already detected all of the sources in this shallow map), now fixing the detection threshold to a peal S/N$\,{>}\,4.5$ (as was done to make the ASAGAO catalogue). We find that the purity at this threshold is 0.9 (which is therefore fairly conservative) and we find a total of 27 additional galaxies, nine of which are new. We use a similar algorithm to measure flux densities, and find similar flux densities compared to those published by ASAGAO (Hatsukade et al., 2018) and GOODS-ALMA (Gómez-Guijarro et al., 2022). As a test, we also extract peaks from the central region of the shallow map and measured their flux densities, and find good agreement with the flux densities measured in our deep map.

The expected uncertainty in our ALMA positions is $\delta\mathrm{RA}=\delta\mathrm{Dec}=0.6\,{\times}\,{\rm FWHM}\,{\div}\,{\rm(S/N)}$ (Ivison et al., 2007), so we expect the $1\,\sigma$ radial uncertainty to be $\delta r\,{\approx}\,$0.76 arcsec/(S/N) (using the geometric mean of the elliptical ALMA beam). Since the probability density of finding a source a distance $r$ from its true position is proportional to $r\mathrm{e}^{-r^{2}/2\delta r^{2}}$, one must go out to a distance of 2.5$\delta r$ in order to find a correct match with 95 per cent certainty. For a given ALMA detection, we thus search the JADES catalogue out to a distance of 1.9 arcsec/(S/N), where S/N here is the peak S/N of each source found above our detection threshold. For the most significant ALMA sources this search radius is unphysically small (much less than 1 pixel), so we apply a minimal search radius of 0.3 arcsec.

We preform a similar counterpart search with the CANDELS catalogue, including the deep infrared data from the HUDF09 survey (Bouwens et al., 2011). Here, we simply apply a uniform search radius of 0.6 arcsec, as we found that the HST F160W morphologies of our DSFGs were often clumpy, leading to unrealistic offsets with respect to our ALMA centroids. We checked by eye that all identified CANDELS counterparts were indeed the same galaxy as the JADES counterparts. It should be noted that Dunlop et al. (2017) found a systematic offset between the CANDELS and ALMA astrometry, which was resolved by applying an offset of about 0.25 arcsec south to the CANDELS positions; the same offset was applied here. In Table 3 we provide the JADES and CANDELS IDs for these matches, and in Appendix B we show the positions of matched sources in our ALMA cutouts overlaid over JWST F356W imaging.

At the position of our ALMA source ID 45 we found that there were two blended galaxies in the longer-wavelength JWST images, yet at shorter JWST wavelengths one of these galaxies was undetected. In the JADES catalogue there was only one ID for these two sources (ID 207277) that had a spectroscopic redshift of 0.332. This spectroscopic redshift likely corresponds to the galaxy brighter at short wavelengths, while our ALMA source is probably the second galaxy that is only detected at longer wavelengths. We therefore use the longer wavelength photometry (where this galaxy is detected) to fit an SED, which results in a photometric redshift of 2.42. Throughout the rest of the paper we report results derived from this fit.

Lastly, our nominal search radius did not return a counterpart to ALMA ID 38, yet this source is located slightly less than 0.5 arcsec from a $z\,{=}\,1.998$, $M_{\ast}\,{\simeq}\,10^{10}\,$M_⊙ JADES galaxy (the JADES ID is 207221). Moreover, the spectroscopic redshift of this source was independently determined to be $z\,{=}\,1.999$ by Boogaard et al. (2023), nearly equal to that reported by JADES, so we assign JADES ID 207221 as the counterpart. There is another ALMA source (ID 44) located less than 2 arcsec south of source 38, whose closest JWST counterpart is in the JADES catalogue and not in our F356W-selected catalogue (the JADES ID is 267661) since it is very close to JADES ID 207221 and suffered blending issues. The reported JADES photometric redshift for ID 44 is 2.02, very close to the spectroscopic redshift of the JADES galaxy found near ID 38, thus it is likely that they are at the same redshift and are undergoing a merger, resulting in complicated morphologies that could easily lead to large offsets between our ALMA centroids and the corresponding JADES positions. We therefore classify ALMA ID 44 as having a JADES counterparts, although we do not fit an SED to it due to the blending issues.

Sources found to have counterparts within their given search radii in both the CANDELS catalogue and the JADES catalogue are marked in Fig. 5 with blue crosses. There are 37 ALMA-detected galaxies with counterparts in both catalogues, seven of which are ALMA sources that have not been previously reported. There are an additional two ALMA-detected galaxies with a counterpart only in the JADES catalogue, both of which are ALMA sources that have not been previously reported. We do not find any CANDELS counterparts with no corresponding JADES counterpart, which highlights one of the benefits of the deeper catalogue of the HUDF made possible with JWST (although the improvement is not dramatic, since CANDELS gets close to identifying all of our mm sources).

This leaves six ALMA-detected sources with no optical or near-infrared counterpart; one of these has previously-published ALMA identifications, the remaining five do not. It is worth pointing out that two of the sources with no counterparts are close to the edge of the map where the noise is the largest (IDs 23 and 25), but the other four are found near the centre of the map where the data are much deeper (IDs 34, 37, 41 and 43).

The number of beams in our ALMA map can be estimated as the area of the map (4.2 arcmin²) divided by the beam area ($2\pi\theta_{\rm maj}\theta_{\rm min}/8\ln 2$), which results in about 8500. From Gaussian statistics, we would expect around 1.5 beams to randomly have a significance greater than $3.6\,\sigma$; however, there are really twice as many statistically independent noise samples due to correlations with the beam (e.g., Condon, 1997; Condon et al., 1998; Dunlop et al., 2017), so really we would expect three false positive detections. The measured number of negative peaks more significant than 3.6 is 14, which is larger than the expectation from a Gaussian distribution – the negative peaks might not be perfectly Gaussian distributed, but this is not surprising given the non-uniform $uv$ coverage. On the other hand, the number of ALMA-detected sources with no optical or near-infrared counterpart (six) is more comparable to the number of negative peaks, and so we cannot rule out the possibility that some of them are false positives.

In Fig. 6 we show the distributions of redshifts, magnitudes, and colours from our JADES catalogue. For the redshifts and colours, we filter out two galaxies which have a S/N$\,{<}\,$5 in the F356W band and ALMA ID 44 due to blending issues as the photometry is not reliable. For the remaining galaxies, the median redshift uncertainty is $\Delta z/(1\,{+}\,z)\,{\simeq}\,0.03$, so we expect them to be accurate. We display the distributions of all 39/45 ALMA galaxies with JWST detections and also highlight the distributions of the eight ALMA galaxies that were not found in previous surveys. The first (top left) panel in Fig. 6 shows the distribution of redshifts (spectroscopic where available, otherwise photometric). We see that the redshift distribution is flat around $z\,{\simeq}\,2$, and our new galaxies appear to follow this trend. Next we show the distributions of the magnitudes of the sources in the NIRCam images (F277W, F335M, F356W, F410M and F444W). We see that most ALMA sources have magnitudes ranging from 19 to 25, with the two faintest sources extending out to ${>}\,$30 (these are the two sources that we do not fit SEDs to). Lastly we show distributions for three NIRCam colours; most ALMA galaxies are fairly red (i.e. brighter at longer wavelengths), and our new, fainter ALMA galaxies tend to follow the same distribution.

As described in Section 2.4.1, the photometric redshifts and stellar masses for most JADES galaxies have been estimated using the extensive multiwavelength imaging available for the HUDF. It is therefore of interest to see how the photometric redshifts and stellar masses of our mm-selected galaxies compare to the typical galaxies found in this field. For our list of sources with a JADES counterpart, we filter out all sources with S/N$\,{<}\,$5 in the F356W band, since these sources will not have reliable SED fits, and ALMA ID 44, which is heavily blended with ALMA ID 38. For the remaining sources, we check to see if there is a spectroscopic redshift, and otherwise use the photometric redshift.

In Fig. 7 we plot stellar mass versus redshift for all JADES galaxies (McLeod et al. in prep.), and highlight our ALMA Band-6-selected sources with good SED fits in red. We find that most of our ALMA-selected sources are high-stellar-mass galaxies between $z\,{=}\,1$ and 3, similar to what was found in earlier ALMA surveys of the HUDF (e.g., Dunlop et al., 2017; McLure et al., 2018; Aravena et al., 2020). In particular, there are 53 galaxies with $M_{\ast}\,{>}\,10^{10}\,$M_⊙, 22 of which have been selected in our ALMA image. Also of note is that at 1.23 mm we are sensitive primarily to $z\,{>}\,1$ objects, since in our sample of 36 galaxies with spectroscopic or photometric redshifts, only three are at $z\,{<}\,1$. In a similar vein, it may also be worth noting that all of the $\log(M_{\ast}/{\rm M}_{\odot})\,{>}\,10.3$ galaxies at $2\,{<}\,z\,{<}\,3$ are detected in our ALMA image.

In order to investigate the difference between galaxies with $M_{\ast}\,{>}\,10^{10}\,$M_⊙ that we have detected with ALMA versus those that we have not detected with ALMA, in Fig. 8 we show the distributions of three select NIRCam magnitudes and three NIRCam colours for the subsamples of 22 ALMA-detected galaxies and 31 ALMA-undetected galaxies. We find that there is no discernible difference in magnitudes between the two samples; however, mm-bright ALMA sources tend to have redder NIRCam colours than the galaxies with similarly high stellar masses that we have not detected with ALMA.

One major benefit of the deep ancillary catalogues available in the HUDF is the ability to stack on the positions of undetected galaxies in different stellar mass and redshift bins, thus estimating the statistical properties of fainter mm-emitting sources that are not individually detectable in our image. Such an analysis was done using the original ASPECS map and HST catalogue (Magnelli et al., 2020), but here the major improvements are a deeper 1-mm map and a deeper catalogue of infrared-selected galaxies from JWST.

To do this, we follow a procedure similar to Simstack (Viero et al., 2013b). However, since we adapt this for a non-circular beam and non-uniform noise distribution, it is worth describing what we do in a little detail. First, we mask out all galaxies that we have detected in our ALMA map (see Table 2), with a source mask set to be $3\times{\rm FWHM}$ in diameter. We also restrict this stacking analysis to the deep central map shown in Fig. 5. Then we subtract the weighted mean of the image – this is crucial, since the ‘stack’ is really the covariance between a map and catalogue (see section 3 of Marsden et al., 2009) and will give a biased result unless the map has zero mean.

Next, we define a grid of four stellar mass bins, logarithmically spaced between $\log\left(M_{\ast}/\mathrm{M}_{\odot}\right)\,{=}\,7.4$ and $\log\left(M_{\ast}/\mathrm{M}_{\odot}\right)\,{=}\,11.4$, with a width of $\Delta\log\left(M_{\ast}/\mathrm{M}_{\odot}\right)\,{=}\,1$. We also define a grid of four redshift bins between 0 and 8, with a width of 2. For the stacking catalogue, we again turn to JADES, using stellar masses and redshifts given by McLeod et al. (in prep.); here redshifts are spectroscopic if available, otherwise photometric. These bins are chosen because they are expected to contain the galaxies comprising the vast majority of the CIB (see Section 4.3 for more details).

For each redshift and stellar mass bin we produce a ‘hits’ map, which is simply a copy of our ALMA map with pixel values set to the number of JADES galaxies within the bin contained within each pixel. We convolve each hits map with the ALMA beam, thereby producing a model image of the sky where the only free parameter of each map is its amplitude, which in this case can be interpreted as the best-fit flux density of the galaxies within the given stellar mass and redshift bin. As was done with the data, we also subtract the weighted mean of the convolved maps.

To estimate the uncertainties in $\hat{S}_{\alpha}$, for each redshift and stellar mass bin we generate 1000 catalogues with the same number of sources in each bin but with random positions, and evaluate Eq. 3 for each one. We find the distribution of values to be well-described by a Gaussian, so we take the standard deviation of the random catalogue flux densities to be the $1\,\sigma$ error in $\hat{S}_{\alpha}$. In order to visually show our results, we also evaluate Eq. 3 after shifting the hits map $N_{\alpha}$ relative to the data map $D$ within boxes of 10 arcsec; this is effectively the cross-correlation between the two maps, where the central pixel is the zero-lag cross-correlation (or the covariance), which is the value we are most interested in.

It is worth noting that a simpler stacking technique (just summing pixel values at the positions of JWST galaxies; see Marsden et al. 2009) produces similar results, merely with slightly less significance. Furthermore, we also tried stacking directly on the dirty map, but we found that the resulting flux densities changed by ${<}1$ per cent.

As a completely separate null test, we also stacked our ALMA map at the positions of 21 sources from the JADES catalogue flagged as stars that happen to fall within our ALMA map. This stack is shown in Fig. 10 and we can see that it is consistent with noise.

In Fig. 9 we see that many high stellar mass bins are blank; this is because we have either detected all of the galaxies within these bins, or there were no galaxies within these bins to begin with. We also see that there are stacked peaks ${\gtrsim}\,3\,\sigma$ in all bins with stellar masses between $10^{8.4}\,{\rm M}_{\odot}$ and $10^{10.4}\,{\rm M}_{\odot}$ and with redshifts between 0 and 4. Across the lowest stellar mass the stacks tend to be positive but with error bars overlapping 0, meaning that galaxies with stellar masses ${<}\,10^{8.4}\,$M_⊙ barely contribute to our ALMA map (a topic that is explored further in Section 4.3). Finally, we note that the $z\,{=}\,$0–2, $M_{\ast}\,{=}\,10^{9.4}$–$10^{10.4}$M_⊙ bin has a very high S/N, thus we could extract more information from this bin by splitting it into smaller sub-bins: for $z\,{=}\,$0–1, $M_{\ast}\,{=}\,10^{9.9}$–$10^{10.4}$M_⊙, we find $S_{\nu}\,{=}\,(22.1\pm 7.9)\,\mu$Jy and $N_{\rm eff}\,{=}\,4$; for $z\,{=}\,$0–1, $M_{\ast}\,{=}\,10^{8.9}$–$10^{9.4}$M_⊙, we find $S_{\nu}\,{=}\,(15.3\pm 3.6)\,\mu$Jy and $N_{\rm eff}\,{=}\,15$; for $z\,{=}\,$1–2, $M_{\ast}\,{=}\,10^{9.9}$–$10^{10.4}$M_⊙ we find $S_{\nu}\,{=}\,(18.2\pm 3.8)\,\mu$Jy and $N_{\rm eff}\,{=}\,11$; and for $z\,{=}\,$1–2, $M_{\ast}\,{=}\,10^{8.9}$–$10^{9.4}$M_⊙, we find $S_{\nu}\,{=}\,(12.1\pm 2.9)\,\mu$Jy and $N_{\rm eff}\,{=}\,24$. These values are consistent with higher stellar mass and lower redshift galaxies being brighter at 1 mm.

In terms of stellar mass, the main conclusions of these stacking results are that: (1) there are about 60 galaxies with stellar masses around $10^{10}\,{\rm M}_{\odot}$ lying just below our $3.6\,\sigma$ ALMA threshold, which have flux densities around 15 $\mu$Jy; (2) there are more than 300 galaxies with stellar masses around $10^{9}\,{\rm M}_{\odot}$ that can also be statistically detected in the ALMA map, with individual flux densities of a few $\mu$Jy; and (3) galaxies with stellar masses around $10^{8}\,{\rm M}_{\odot}$ or below have flux densities that are too low to be detected in the ALMA map, even statistically.

An important question is whether the intensity of the CIB can be recovered by summing the contribution from known galaxies, or if there exists an additional population of sources or a genuinely diffuse component of the CIB. This question has been addressed by many previous studies at wavelengths around 1 mm (e.g. Penner et al., 2011; Viero et al., 2013b; Fujimoto et al., 2016; Dunlop et al., 2017; Hatsukade et al., 2018; González-López et al., 2020; Gómez-Guijarro et al., 2022; Chen et al., 2023), with results either in the tens of per cent range, or requiring a model fit to the number counts and extrapolated below the flux densities of the faintest galaxies detected. Here we explore this question with our new ALMA map at 1.23 mm and our new JADES catalogue, to see if we can recover the total CIB intensity solely with directly-detected galaxies. To start with, the sum of our detected source flux densities (${>}\,3.6\,\sigma$) is (7.33$\,{\pm}\,$0.10) mJy, and these sources are detected within an area of $1.15\,{\times}\,10^{-3}\,{\rm deg}^{2}$. Therefore we have resolved a total intensity of $(6.35\,{\pm}\,0.09)\,{\rm Jy}\,{\rm deg}^{-2}$ in individually-detected DSFGs.

In addition to this, our stacking analysis (Fig. 9) demonstrates that our map is also statistically sensitive to fainter galaxies. Multiplying $\hat{S}_{\alpha}$ by $N_{\rm eff}$ in each bin shown in Fig. 9, and summing, yields a flux density of $(2.6\,{\pm}\,0.5)$ mJy, or an intensity of $(2.4\,{\pm}\,0.5)\,{\rm Jy}\,{\rm deg}^{-2}$ (where the area used in our stack is $1.11\,{\times}\,10^{-3}\,{\rm deg}^{2}$, slightly less than the full map due to our source mask); this stacking result is larger than has previously been possible, due to the combination of a deeper ALMA map and larger catalogue from JWST. Lastly, noting that there are no galaxies with $S_{1230}\,{>}\,0.85$ mJy present in our map, we could also add a contribution from brighter sources. The GOODS-ALMA survey (Gómez-Guijarro et al., 2022) covered a much larger area with ALMA at 1 mm (about 72 arcmin²) and found 22 galaxies with $S_{1230}\,{>}\,0.85$ mJy, corresponding to an intensity of $(1.31\,{\pm}\,0.02)\,{\rm Jy}\,{\rm deg}^{-2}$ (including the factor of 0.8 to convert their flux densities from 268 GHz to 243 GHz) and so we can add this to our resolved CIB contribution as well.

In total, we directly or statistically detect a total CIB intensity of $(8.7\,{\pm}\,0.5)\,{\rm Jy}\,{\rm deg}^{-2}$ in our map, and by adding in a contribution from brighter sources, we estimate that the true value of the CIB is $(10.0\,{\pm}\,0.5)\,{\rm Jy}\,{\rm deg}^{-2}$. We note that if we instead perform our stacking analysis on the full map (without masking bright detected sources) we obtain a CIB estimate of $(8.7\,{\pm}\,0.8)\,{\rm Jy}\,{\rm deg}^{-2}$, consistent with our approach of detecting bright objects then masking them to add the stacking result. Now we have to determine what fraction of the background we have accounted for. While we have not performed any deboosting corrections for the individual flux densities measured in our catalogue, the fact that we find a nearly identical result when stacking on the full map with no masking shows that any bias from the flux density boosting of weaker S/N sources has very little effect on our conclusions about the total CIB intensity.

Estimating the absolute level of the CIB at 1 mm is subject to larger uncertainties than our estimate of the amount contributed by sources. The spectrum of the CIB was measured by COBE-FIRAS (Fixsen et al., 1998), but with fairly large systematics-dominated uncertainties. More recently Odegard et al. (2019) used Planck in combination with FIRAS to estimate the total intensity of the CIB in the Planck-High-Frequency Instrument channels, including at 217 and 353 GHz (the closest Planck frequencies to the mean frequency of our ALMA image, 243 GHz). The best-fit CIB spectral shape from Fixsen et al. (1998) is a modified blackbody function of the form

where $A$ is a constant, $\nu_{0}$ is a fiducial frequency, $\beta\,{=}\,0.65$, $B_{\nu}$ is the Planck function and $T_{\rm d}\,{=}\,18.5\,$K. The uncertainties in the three fit parameters $A$, $\beta$ and $T_{\rm d}$ are significantly correlated (the correlation coefficients reported in the paper are larger than 95 per cent), meaning that there is effectively only one free parameter in the fit. Lacking the data needed to do the fit ourselves, we simply fix $\beta$ and $T_{\rm d}$ (thus fixing the shape of the CIB spectrum), but take into account the uncertainties in the amplitude throughout our calculations.

With $T(\nu)$ in hand, we determine the effective frequency of the CIB in our image using the modified blackbody parameters fit to the CIB spectrum from FIRAS (Fixsen et al., 1998). Specifically, we calculate

with $I(\nu)$ from Eq. 6 and the constants $A$ and $\nu_{0}$ cancel out in the ratio of integrals. We find a value of $\nu_{\rm eff}\,{=}\,241.4\,$GHz, which is close to the ALMA data combined central frequency.

We finally estimate the amplitude of the CIB at $\nu_{\rm eff}$ from Eq. 6, using two different approaches to obtain the amplitude $A$. The first approach is to use $A\,{=}\,(1.3\pm 0.4)\,{\times}\,10^{-5}$, which is a direct result from Fixsen et al. (1998), assuming that the shape is fixed. The second approach is to fit a power law between the CIB intensities from Odegard et al. (2019) at 217 and 353 GHz (there is an exact solution because there are two measurements and two free parameters), then use the fit to interpolate the intensity at 241.4 GHz. To propagate uncertainties, we draw 10,000 Gaussian random numbers from the measured amplitudes and calculate the corresponding CIB value at 241.4 GHz.

We find absolute CIB values of $19_{-5}^{+6}\,{\rm Jy}\,{\rm deg}^{-2}$ using the fit from Fixsen et al. (1998) and $22_{-6}^{+7}\,{\rm Jy}\,{\rm deg}^{-2}$ using the interpolation from Odegard et al. (2019), where the error bars are the 68 per cent confidence intervals of the posterior distributions and the central values are the means of the distributions.

We show these estimates graphically in Fig. 11. Specifically we plot with a blue band the 68 per cent range of the absolute value of the CIB estimated using FIRAS data alone (Fixsen et al., 1998) and with a pink band we plot the estimate using FIRAS data in combination with Planck for foreground removal (Odegard et al., 2019). Although the Odegard et al. (2019) results substantially shrink some uncertainties compared to the Fixsen et al. (1998) results, that is really only the case at shorter wavelengths; at 1.2 mm the uncertainties are similar, but the background is actually a little higher in the more recent paper. The difference between the two estimates indicates that the background at these wavelengths is still quite uncertain.

In Fig. 11 we also present the contribution of sources to the CIB, by plotting the cumulative intensity of resolved galaxies from this work as a function of their flux density, including the estimated contribution from $S_{1230}\,{>}\,0.85$ mJy galaxies and also the stacking-estimate contribution from galaxies in the JADES catalogue. For reference we show the same analysis using the results from Dunlop et al. (2017) and ASPECS (González-López et al., 2020), where in both works an estimate of the contribution to the CIB from stacking was performed. Our new results are similar to those from previous studies, with our new analysis detecting the faintest sources and finding the highest total background value in the HUDF region.

Clearly the absolute value measurements of the CIB are larger than what we find by summing the flux densities of known galaxies (detected by both ALMA and JWST). It is therefore important to consider whether or not the two kinds of measurement are consistent.

There are fluctuations in the CIB (e.g., Planck Collaboration XXX, 2014) whose amplitude depends on the area observed. Studies of the CIB at submm wavelengths found that $\delta I/I\,{\simeq}\,15\,$per cent on scales around 10 arcmin (Viero et al., 2009, 2013a), which is larger than the area of our deep ALMA map. To find a better estimate of the expected amplitude of these fluctuations for maps the size of our HUDF image, we use the Simulated Infrared Dusty Extragalactic Sky (SIDES) mock catalogue (Béthermin et al., 2017; Gkogkou et al., 2023). Briefly, SIDES uses dark matter halos in a simulated light cone to obtain clustered positions on the projected sky, then attaches stellar masses and far-infrared SEDs to the halos using a two star-formation-mode model of galaxy evolution. The total simulated area of the SIDES simulation is 117 deg², but here we only use a 1 deg² tile from the full simulation.

Since each simulated galaxy has an SED, we first estimate their 243-GHz flux densities by integrating their SEDs through the transmission function derived above. We next take 100 random patches from the 1 deg² tile, each with the area of our image of the HUDF (here 4.2 arcmin²). The CIB intensity of each random patch is then just the sum of the flux densities of the galaxies in the patch divided by the area. However, since we find no galaxies with $S_{1230}\,{>}\,0.85$ mJy in our real ALMA map, then we subtract these from our random patches as well; this assumes that there are no additional fluctuations from the population of $S_{1230}\,{>}\,0.85$ mJy galaxies (i.e. we simply add a constant for the bright part of the background). Additionally, we neglect the contribution from galaxies at stellar masses where we were unable to detect any statistical signal in the real data.

Following this procedure, we find that the mean CIB value within 4.2 arcmin² patches of the SIDES simulation (after subtracting the contribution from $S_{1230}\,{>}\,0.85$ mJy galaxies) is 11.1 Jy deg^-2, comparable to the value of $(8.6\,{\pm}\,0.5)\,{\rm Jy}\,{\rm deg}^{-2}$ that we have measured. The standard deviation of the SIDES simulation patches is 1.9 Jy deg^-2, thus fluctuations are $\delta I/I\,{=}\,17\,$per cent. We must therefore take this into account when comparing the mean CIB value from FIRAS averaged over nearly the whole sky to the single 4.2 arcmin² patch we have observed.

We now estimate the probability of measuring a CIB value of $(10.0\,{\pm}\,0.5)\,{\rm Jy}\,{\rm deg}^{-2}$ assuming that the true value is either $19_{-5}^{+6}\,{\rm Jy}\,{\rm deg}^{-2}$ (Fixsen et al., 1998) or $23_{-8}^{+6}\,{\rm Jy}\,{\rm deg}^{-2}$ (Odegard et al., 2019) and that fluctuations are 17 per cent. We take our 10,000 CIB absolute value realizations discussed above and additionally draw 10,000 Gaussian-distributed values for a fluctuation amplitude, then multiply these together. We draw 10,000 Gaussian-distributed numbers for our measured CIB value, and take the difference between these and the possible absolute values. Finally, our statistic is the fraction of the area of the resulting posterior distribution that is less than zero, which can be interpreted as the probability of obtaining our actual measured CIB value or less while taking into account both the measurement uncertainties and intrinsic CIB fluctuations. We find that the probability of measuring a CIB value of $(10.0\,{\pm}\,0.5)\,{\rm Jy}\,{\rm deg}^{-2}$, given the absolute CIB measurement from Fixsen et al. (1998) is 8.9 per cent, or 5.3 per cent given the absolute CIB measurement from Odegard et al. (2019).

Assuming that the absolute value measurements from FIRAS are correct and that our measurement is also correct, we can calculate the required level of statistical excursion (in units of $\sigma\,{=}\,0.17$) corresponding to the HUDF. Using the same 10,000 random values for the absolute FIRAS values and the measured values, we find that the CIB fluctuation at the position of the HUDF must be $-2.4_{-1.4}^{+0.4}\sigma$ using the CIB value from Fixsen et al. (1998), or $-2.9_{-1.2}^{+0.3}\sigma$ from Odegard et al. (2019) (here the central values are the means of the posterior distributions and the error bars are 68 per cent confidence intervals). What this means is that the variance in HUDF fields is large enough that our results can explain the whole of the CIB provided that the HUDF happens to be a relatively mild (${\simeq}\,2\,\sigma$) underdense direction on the sky.⁵⁵5The HUDF was not selected entirely randomly; however, there is no particular reason to believe that the criteria used for its selection would make it likely to have a lower than average background (Beckwith et al., 2006).

As a final check, we use the simulated catalogue of galaxies from SIDES to estimate the fraction of the CIB emitted by galaxies with stellar masses between 10^7.4 M_⊙ and 10^11.4 M_⊙ and with redshifts between 0 and 8 (namely the parameter space over which we stacked on undetected JADES galaxies). For each of the 100 random patches described above, we also calculate the sum of the flux densities from all galaxies within our stacking range divided by the sum of all of the galaxies in the HUDF-sized region. We find that the average ratio is 97 per cent, with a standard deviation of 1 per cent. Thus if we accept that SIDES provides a reasonable model for counts at these wavelengths, we do not expect that we are missing a significant contribution to our ALMA measurement of the CIB from even fainter and lower stellar mass galaxies.

Turning back to the data, if we include sources detected by ALMA, the contribution to the CIB from 10^10.4–10^11.4 M_⊙ galaxies is $(3.5\,{\pm}\,0.1)\,{\rm Jy}\,{\rm deg}^{-2}$, from 10^9.4–10^10.4 M_⊙ galaxies is $(3.2\,{\pm}\,0.1)\,{\rm Jy}\,{\rm deg}^{-2}$, from 10^8.4–10^9.4 M_⊙ galaxies is $(1.1\,{\pm}\,0.2)\,{\rm Jy}\,{\rm deg}^{-2}$ and from 10^7.4–10^8.4 M_⊙ galaxies is $(0.5\,{\pm}\,0.4)\,{\rm Jy}\,{\rm deg}^{-2}$. If we stack on the next lowest mass bin (10^6.4–10^7.4 M_⊙) over the full redshift range we obtain a CIB contribution of $(-0.5\,{\pm}\,0.5)\,{\rm Jy}\,{\rm deg}^{-2}$. This is consistent with the CIB having essentially converged over the stellar mass range that we have probed.

A potential issue here is the completeness of the JADES catalogue within our stacking range. For our 1000 random SIDES patches we also keep track of the total number of galaxies with stellar masses and redshifts within our stacking region. We find that the average number of galaxies is 1890 (with a standard deviation of 150). In the JADES catalogue there are 1856 galaxies in our ALMA image with stellar masses and redshifts within our stacking range, which is consistent with the total number of galaxies expected from the SIDES simulation. For reference, there are 1561 galaxies from the CANDELS catalogue that are in our ALMA image footprint within the same stacking range. Adding JADES has helped us to find more of the CIB within our ALMA map, and it seems that going even deeper in the optical/near-infrared will not add significantly to the source-derived CIB estimate.

What these numbers indicate is that our estimate of the 1.23-mm CIB (from individually-detected galaxies, together with statistical stacking results) appears to contain essentially all of the possible galaxies that would contribute to the CIB, and that it is genuinely lower than what the mean value is estimated to be. However, the chances that our small patch of the extragalactic sky falls on a negative fluctuation of the CIB are not small enough to rule out the hypothesis that we have indeed recovered essentially the entire CIB from these known galaxies.