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 combine archival ALMA data targeting the Hubble Ultra Deep Field (HUDF) to produce the deepest currently attainable 1-mm maps of this key, extragalactic survey field. Combining all existing data in Band 6, our deepest map covers 4.2arcmin^2, with a beamsize of 1.49"x1.07"at an effective frequency of 243GHz (1.23mm). It reaches an rms of 4.6uJy/beam, with 1.5arcmin^2 below 9.0uJy/beam, an improvement of>5% over the best previously published map and 50% improvement in some regions. We also make a wider, but shallower map, covering 25.4arcmin^2. We detect 45 galaxies in the deep map down to 3.6sigma, including 10 more 1-mm sources than previously detected. 38 of these galaxies have a JWST ID from the JADES NIRCam imaging and the new sources are typically faint and red. A stacking analysis on the positions of ALMA-undetected JADES galaxies yields detections for z<4 and stellar masses from 10^(8.4) to 10^(10.4)Msun, extracting 10% of additional stacked signal from our map compared to previous analyses. Detected sources and stacking contribute (10.0+/-0.5)Jy/deg^2 of the cosmic infrared background (CIB) at 1.23mm. Although this is short of the (uncertain) background level of about 20Jy/deg^2, after taking into account intrinsic fluctuations in the CIB, our measurement is consistent with the background if the HUDF is a mild (~2sigma) negative fluctuation. This suggests that within the HUDF, JWST may have detected essentially all of the galaxies that contribute to the CIB. Our stacking analysis predicts that the field contains around 60 additional galaxies with 1.23mm flux densities averaging around 15uJy, and over 300 galaxies at the few uJy level. However, the contribution of these fainter more modestly-obscured objects to the background is small, and converging, as anticipated from the now well-established strong correlation between galaxy stellar mass and obscured star formation.


INTRODUCTION
The Hubble Ultra Deep Field (HUDF) is probably the most wellstudied extragalactic region of the sky, containing some of the deepest optical and near-IR exposures obtained to-date (e.g.Beckwith et al. 2006;Illingworth et al. 2013;Koekemoer et al. 2013;Eisenstein et al. 2023).Studies of the HUDF can tell us about how galaxies have evolved from the earliest times to the present day.We now know that a significant fraction of the cosmic star formation occurred within distant galaxies full of dust (e.g.Casey et al. 2015;Koprowski et al. 2017), making them effectively invisible at optical and near-IR wavelengths, for all but the deepest images.Yet these galaxies are ★ E-mail: ryleyhill@phas.ubc.cabright at millimetre (mm) and submillimetre (submm) wavelengths, and so in order to complete our understanding of galaxy evolution, we must also survey the HUDF in this waveband.
The best telescope at these wavelengths available today is the Atacama Large Millimeter/submillimeter Array (ALMA).ALMA has been used several times to survey the HUDF at wavelengths around 1 mm (Dunlop et al. 2017;ASAGAO, Hatsukade et al. 2018;ASPECS, González-López et al. 2020;GOODS-ALMA, Gómez-Guĳarro et al. 2022), and it has also been used to follow up interesting individual targets within the HUDF (e.g.Fujimoto et al. 2017;Cowie et al. 2018).Yet in contrast to the tens of thousands of galaxies detected in the optical, only a few dozen mm-bright galaxies have been found in these follow-up observations so far.Although often referred to as 'SMGs' (for 'submillimetre galaxies'), since we are working here at wavelengths slightly longer than 1 mm, we will refer to them with the more generic acronym DSFG (for 'dusty starforming galaxy').
In order to detect more galaxies around 1 mm with ALMA, we can combine all of the data into a single image.In a similar way, two decades ago, a 'super-map' of the GOODS-N field was constructed by combining SCUBA data sets at 850 m (Borys et al. 2003), and multiwavelength counterparts from the Hubble Space Telescope (HST) were identified and used to study the properties of the detected submm-luminous sources (Pope et al. 2005).With this same motivation we have undertaken an archival project to combine all of the ALMA observations of the HUDF taken around 1 mm, with the goal of finding new DSFGs, identifying their counterparts, assessing how much of the background light has been resolved and providing an image that can be used by others for stacking analyses.
In Section 2 we describe how we retrieve the data from the ALMA archive and combine it into single continuum images in  space.In Section 3 we outline an alternative approach to combining data subimages in real space, which provides a test of the reliability of our -space combination.In Section 4 we provide our new galaxy catalogue and compare it to previous catalogues.In Section 5 we describe how our results lead to a new estimate of the resolved fraction of the cosmic infrared background (CIB) at 1 mm.We discuss in Section 6 improvements in our data products and results compared to what was previously available at 1 mm in the HUDF and we conclude in Section 7. Appendix A presents a wider (and shallower) map and gives a supplementary list of sources, while Appendix B provides cutouts of our ALMA sources overlaid over JWST F356W imaging.

Obtaining archival ALMA data
We focus on ALMA Band 6, which spans a wavelength (frequency) range of 1. .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 archive1 for all Band-6 observations centred at 03 h 32 m 39.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 A1.
For each of the observations, we downloaded the raw  data and calibrated it using the provided ScriptForPI and the CASA (Mc-Mullin 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.

Obtaining archival ALMA images
In addition to downloading the  data for each programme, we also obtained the imaging products made available by the observatory.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 -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 primarybeam-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 primarybeam images.

Data combination in the uv plane
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  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 all of the  points to the average frequency of the data being combined; here, the average frequency is about 243 GHz (or 1230 m).We also chose natural weighting (each  point is weighted by its instrumental noise), but with a 250 k  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  taper, setting the pixel size to 0.18 arcsec (because the synthesized beam is somewhat smaller).Lastly, we cleaned the image down to 20 Jy beam −1 , similar to the cleaning level chosen by ASPECS.The final synthesized beamsize of the -tapered map is 1.49 arcsec × 1.07 arcsec, while for the map with no  taper the synthesized beamsize is 1.32 arcsec × 0.92 arcsec.The pixel rms in the tapered image (which we use for further analysis in this paper) has a minimum value of 4.6 Jy beam −1 in the deepest region, while the rms is less than 9.0 Jy beam −1 over an area of 1.5 arcmin 2 .
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 A1 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 a Approximate map rms estimated using the observatory MFS data products.For programmes with a single tuning, this is the rms of the primary-beam-uncorrected map after masking all known sources.For programmes with multiple tunings, we follow the same procedure and then estimate the rms of the weighted mean of the images as 1/ Σ  1/2  .The actual rms from combining images at different tunings in  space using the MFS mode will generally be less than this estimate.b Maximum synthesized beam across all frequencies of a given data set.The actual synthesized beam from combining data taken in different tunings using the MFS mode may differ slightly from this estimate.c The sky coverage overlapping with our data selection criteria, which may be less than the total sky coverage of the given programme.d Programmes used to make a mosaic with uniform noise.e The imaging products from this data set were not used to create the combined image in the image plane because the beamsize is much larger than for the other data sets.
survey and the GOODS-ALMA survey.We set the pixel size to 0.06 arcsec for the map with no  tapering, and 0.12 arcsec for the map with 250 k  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 Jy beam −1 (similar to the level chosen by ASAGAO).The final synthesized beamsize of the -tapered map is 0.87 arcsec × 0.64 arcsec, while for the map with no  taper the synthesized beamsize is 0.36 arcsec × 0.27 arcsec.The pixel rms in the tapered image (which we also use for further analysis in this paper) has a typical value of 60 Jy beam −1 outside of individual pointings, and 20 Jy beam −1 within individual pointings.
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 A1, respectively).We used the same tclean parameters as above, generating versions with no  tapering and with 250 k  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-Guĳarro 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. 2 The final signal map, noise map and signal-to-noise ratio (S/N) map for the deep central mosaic with a 250 k taper, including individual pointings, is shown in Fig. 1, while the same maps for the larger, shallower mosaic are shown in Fig. A1.The deep central mosaic covers an area of 4.2 arcmin 2 out to our primary beam threshold of 0.2, while the large shallow mosaic covers an area of 25.4 arcmin 2 out to our primary beam threshold of 0.4.

Optical and infrared galaxy catalogues
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 field3 (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-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 available4 (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 m to 4.4 m, and we chose to select our catalogue at 3.5 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 and photometry, with the CAN-DELS catalogue used to demonstrate new improvements thanks to JWST.

JADES photometric redshifts and stellar masses
In order to construct photometry catalogues, we ran SourceExtractor (Bertin & Arnouts 1996) in dual-image mode, with NIRCam F356W used as the detection image.All imaging was homogenized to match the point-source function (PSF) of the F444W imaging in order to minimise any colour systematics arising from differences in the PSF.We performed isophotal photometry on all available JADES and JEMS NIRCam bands, as well as ancillary HST ACS imaging from the Hubble Legacy Fields GOODS-South data release (see Illingworth et al. 2013;Whitaker et al. 2019) and ground-based VIMOS -band imaging over GOODS-South (Nonino et al. 2009).To extract robust photometry from this lower-resolution imaging, we utilised TPHOT (Merlin et al. 2015), using positional and surfacebrightness information from the higher-resolution NIRCam imaging as input.
To calculate photometric redshifts and hence stellar masses, we performed spectral-energy-distribution (SED) fitting using Le-Phare (Arnouts & Ilbert 2011), with a Bruzual-Charlot (Bruzual & Charlot 2003) template set, incorporating a Chabrier (Chabrier 2003) initial mass function.The template set includes a Calzetti (Calzetti et al. 2000) dust-attenuation law, with range of reddening   = [0.0,6.0, 0.2].We used two metallicities, 0.2 Z ⊙ and Z ⊙ , and included exponentially declining star-formation histories with  ranging between 0.1 and 15 Gyr.Where the object has a known spectroscopic redshift, we fixed to this redshift in obtaining the stellar mass.
The median uncertainty in our photometric redshifts for all of our HUDF galaxies is Δ/(1 + ) ≃ 0.06.For stellar masses, in this paper we only care about relative uncertainties, since we only use this parameter to sort galaxies into different bins.The absolute uncertainties in stellar masses are expected to be large, but the relative uncertainties should be small, so ignore them.

Data combination in the image plane
In addition to our -combined images, we also tried combining archival ALMA data-product images (each one being produced individually using tclean in MFS mode) in real space.This provides a consistency check for how well we are able to combine in  space the data spanning a wide frequency range and large array configuration range.For this test we focus on the deep central mosaic (following the ASPECS footprint), where the data sets included are listed in Table 1.
The first step in this process is to convolve each image to the same resolution.A common beamsize was selected as the largest beam across all of the data sets, and for each image the convolution kernel required to produce this common beam was found numerically using the create_matching_kernel function from the python photutils module (Bradley et al. 2022).Each image was then convolved with this kernel.We expect that convolution to a lower resolution will effectively decrease the sensitivity to point sources.The largest beam across all of our data comes from the ASPECS pilot programme (2013.1.00718.S, see Aravena et al. 2016), which has a beam of about 1.5 arcsec × 0.8 arcsec.The most sensitive map (across a significantly larger area) comes from the full ASPECS survey (2016.1.00324.L, see González-López et al. 2020), where the beam is about 1.2 arcsec × 0.8 arcsec, or the second-largest beam across all of our images.We found that degrading the resolution of the full ASPECS map to match the ASPECS pilot map ultimately produces a map with a higher rms (after combining all the data) compared to removing the ASPECS pilot map from the combination and convolving each image to match that of the full ASPECS map; our final combined image therefore does not contain data from the ASPECS pilot, although this amounts to only a small loss of data (see Table 1).
We next created a pixel grid for the image.The size of the pixels in the grid was chosen as the largest pixel size within the set of 61 images (in this case 0.2 arcsec), and the grid was set to span 5 arcmin, centred at 03 h 32 m 39.0 s −27 • 47 ′ 29.1 ′′ .Next, each of the 61 convolved images (minus those from the ASPECS pilot) were reprojected onto this grid using the python function reproject_interp, part of the reproject module (Robitaille et al. 2020).
The next step was to create a noise map for each reprojected and convolved image.To do this, we created primary-beam-uncorrected images (simply the primary beam image multiplied by the primary beam), and for each image we created a mask using the same method described above.We then calculated the rms of the primary-beamuncorrected map multiplied by the mask, and divided this by the primary beam.
These images should now be aligned to the same pixel grid, have the same resolution, and have associated noise maps.The final step is to combine them, weighting each pixel by its inverse variance.In order to remove boundary artefacts associated with combining images with very different noise levels, an apodizing mask was applied to the edges of each image in the combination.The apodizing function chosen was a Gaussian with a standard deviation of the beamsize, and the apodization was applied out to 3 times the beamsize from the edge of each map.We investigated various options for apodization (e.g. a cosine function or different apodization lengths), and found that this simple choice removed most of the obvious edge effects.
Lastly, a new noise map was calculated for the combined image using the same algorithm used to calculate a noise map for the combined image.For reference, the image-combined map can also be downloaded.5

Comparison between uv-plane combination and image-plane combination
The optimal way to combine interferometric images is to add the observations in the  plane, then image the entire data set.However,  1.This region is defined by the contour where the primary beam reaches 0.2, covering a total of 4.2 arcmin 2 .Middle: Corresponding noise map.Right: S/N map, from dividing the signal map by the noise map.
the different frequencies and  sampling of the various observations may lead to unwanted behaviour.Moreover, adding the individual pointings is much more straightforward in the image plane than in  space.Therefore we would like to compare the properties of the two images (-space combination and real-space combination) to check that they are consistent with one another, and to see which one performs better.
We focus on the region defined by our deep central map with a 250 k taper, which was made going out to 0.2 times the primary beam (see Fig. 1).We extracted all sources with S/N > 4 within this region from both maps using the same source-extraction procedure as described in Section 4; this threshold was simply chosen so that enough overlapping sources could be extracted for comparison.There are 36 peaks with this significance in the -space combined map and 29 peaks in the image-space combined map.In Fig. 2 we show all of the sources found with S/N > 4, from which it can be seen that all bright and obvious sources are in agreement.However, we do see more sources detected in the  combination compared to the image combination, especially around the edges of our selected region, where the  map does much better.
To quantify the difference between the two maps, in Fig. 3 we plot the ratio of the peak S/N from the  combination to the peak S/N from the image combination as a function of 1-mm flux density in the -space combined map.For the galaxies with S/N > 4 in the -combined map but not the image-combined map, we simply extract the value of the image-combined map at the location of the  detections in order to include them on this plot.We find a relatively uniform improvement in S/N from combining the  data as a function of source brightness.The mean S/N of all matching sources is higher by about 1.2.Of course there is a balance between the increased complexity and computing resources required to combine the data in the  plane versus combining in the image plane, but our results indicate that the improvement from combining the data in the -plane is worthwhile.

Cataloguing 1-mm sources
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 Jy beam −1 , which is 50 per cent deeper than any previous map.Over an area of 1.5 arcmin 2 our new map has an rms below 9.0 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 Jy beam −1 over the full 4.2 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 value of −4.1 , and there are 35 positive peaks brighter 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 − N neg /N 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. 4 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 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 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.The positions of sources extracted from this search are shown in Fig. 6, with positions and 1-mm flux densities given in Table 2. Appendix B contains 5 arcsec × 5 arcsec cutouts overlaid on JWST F356W images.For single-dish surveys, DSFGs are almost always unresolved, but this is not the case for ALMA data.Hence we need to decide how to quote brightness values when some DSFGs are resolved.For the flux densities, we follow a procedure similar to the ASAGAO survey; for each source we fit a 2-dimensional Gaussian profile, fixing the position to the peak pixel and the amplitude to the value of the peak pixel, but allowing the size, ellipticity and position angle to vary.We then calculate the number of beams contained within each source as where  and  are the best-fit major and minor FWHM values, and  maj and  min are the synthesized beam major and minor FWHM, respectively.If the number of beams is greater than 1, we calculate the integrated flux density at 243 GHz, or 1230 m, as otherwise the integrated flux density is simply the peak pixel value.
We also flag fits where the minor/major axis ratio is less than 0.5 as bad fits, and use peak pixel values for these sources.Uncertainties are taken from the noise map, and uncertainties from the fits are propagated to the integrated flux densities.All of our sources and 1-mm flux densities are listed in Table 2.
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  = 1.5, a dust temperature of 30 K, and a redshift of 1.5; the correction factor is 1.3.Similarly, the mean frequency of the map from GOODS-ALMA (Gómez-Guĳarro 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. 5, 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.
In Fig. 6 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-Guĳarro et al. 2022).It appears that most of the detections in our combined map coincide with published sources, but there are 13 new DSFGs.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 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-Guĳarro 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.

Matching mm-selected sources to near-IR-selected sources
The expected uncertainty in our ALMA positions is RA = Dec = 0.6 × FWHM ÷ (S/N) (Ivison et al. 2007), so we expect the 1  radial uncertainty to be  ≈ 0.76 arcsec/(S/N) (using the geometric mean of the elliptical ALMA beam).Since the probability density of finding a source a distance  from its true position is proportional to e − 2 /2  2 , one must go out to a distance of 2.5 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 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 1.5 pixels.
We preform the same counterpart search with the CANDELS catalogue, including the deep infrared data from the HUDF09 survey (Bouwens et al. 2011).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.
We note that one ALMA source (ID 44) did not have a counterpart in our F356W-selected catalogue yet had a counterpart in the JADES catalogue, so we include this ID here.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 galacxies 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.
Sources found to have counterparts within their given search radii in both the CANDELS catalogue and the JADES catalogue are marked in Fig. 6 with blue crosses.There are 35 ALMA-detected galaxies with counterparts in both catalogues, seven of which are ALMA sources that have not been previously reported.There are an shown as brown circles.The gold contour shows the footprint of the deepest HST data from the CANDELS/HUDF09 survey, while the entire region is covered by the JADES survey.Galaxies with counterparts in both the CANDELS catalogue and the JADES catalogue are indicated with a blue cross and galaxies with only a JADES counterpart are indicated with a red cross.There are no galaxies with a CANDELS/HUDF09 counterpart but no JADES counterpart.Numbers refer to the labels in Tables 2 and 3.
Table 2. Positions and flux densities ( 1230 ) for the 45 sources found in the ALMA 1-mm combined image at S/N > 3.6.Here flux densities are measured by fitting Gaussian profiles to the sources in order to estimate the number of beams per source, then scaling the peak pixel value by the number of beams; for sources consistent with 1 beam or less, we set the flux density to be the peak pixel value.For galaxies with matches in the JADES catalogue and sufficient photometry to fit SEDs (S/N > 5 in the F356W band) we provide stellar masses and photometric redshifts from McLeod et al. (in prep.).Galaxies with spectroscopic redshifts from Rieke & the JADES Collaboration ( 2023) are indicated by asterisks, which we provide instead of photometric redshifts.The median photometric redshift uncertainty of our ALMA sample is Δ/(1 + ) ≃ 0.03, while the relative uncertainties in stellar masses are small.Corresponding identifications in other published catalogues are given in Table 3. ALMA ID 45 is at the position of a low- and a high- galaxy.We assume that this ALMA source corresponds to the high- galaxy, the photometric redshift and stellar mass are derived using only long-wavelength JWST photometry where the high- galaxy is detected.additional four ALMA-detected galaxies with a counterpart only in the JADES catalogue, one of which is an ALMA source that has 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 seven ALMA-detected sources with no optical or nearinfrared counterpart; two of these have 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 five are found near the centre of the map where the data are much deeper (IDs 34,37,38,41 and 43).

Name
The number of beams in our ALMA map can be estimated as the area of the map (4.2 arcmin 2 ) divided by the beam area (2 maj  min /8 ln 2), which results in about 8500.From Gaussian statistics, we would expect around 1.5 beams to randomly have a sig-Table 3. Cross-matched identifications between our catalogue of sources and previous studies.Near-infrared matches are from JADES (Rieke & the JADES Collaboration 2023) and CANDELS (Guo et al. 2013).Matches with other ALMA Band-6 surveys are from Dunlop et al. (2017), ASAGAO (Hatsukade et al. 2018), GOODS-ALMA (Gómez-Guĳarro et al. 2022), andASPECS (González-López et al. 2020).ALMA ID 44 is not detected in our F356W-selected catalogue, but has a counterpart in the JADES catalogue, which we provide below.The JADES ID corresponding to ALMA ID 45 is a blend of a low- and a high- galaxy, and we assume that our ALMA source corresponds to the galaxy at high-.nificance greater than 3.6 ; 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  coverage.On the other hand, the number of ALMA-detected sources with no optical or near-infrared counterpart (seven) is more comparable to the number of negative peaks, and so we cannot rule out the possibility that some of them are false positives.

Name
In Fig. 7 we show the distributions of redshifts, magnitudes, and colours from our JADES catalogue.For the redshifts and colours, we filter out three galaxies which have a S/N < 5 in the F356W band as the photometry is not reliable.For the remaining galaxies, the median redshift uncertainty is Δ/(1 + ) ≃ 0.03, so we expect them to be accurate.We display the distributions of all 38/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. 7 shows the distribution of redshifts (spectroscopic where available, otherwise photometric).We see that the redshift distribution is flat around  ≃ 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 26 to 32, with the three faintest sources extending out to 40 (these are the three 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.

Stellar mass-redshift 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.For the remaining sources, we check to see if there is a spectroscopic redshift, and otherwise use the photometric redshift.
In Fig. 8 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 ALMAselected sources are high-stellar-mass galaxies between  = 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  * > 10 10 M ⊙ , 21 of which have been selected in our ALMA image.Also of note is that at 1.23 mm we are sensitive primarily to  > 1 objects, since in our sample of 35 galaxies with spectroscopic or photometric redshifts, only three are at  < 1.In a similar vein, it may also be worth noting that all of the log( * /M ⊙ ) > 10.3 galaxies at 2 <  < 3 are detected in our ALMA image.
In order to investigate the difference between galaxies with  * > 10 10 M ⊙ that we have detected with ALMA versus those that we have not detected with ALMA, in Fig. 9 we show the distributions of three select NIRCam magnitudes and three NIRCam colours for the subsamples of 21 ALMA-detected galaxies and 32 ALMAundetected 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.

Stacking on near-IR-selected positions
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.
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 × FWHM in diameter.We also restrict this stacking analysis to the deep central map shown in Fig. 6.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 ( * /M ⊙ ) = 7.4 and log ( * /M ⊙ ) = 11.4,with a width of Δ log ( * /M ⊙ ) = 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 5.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.
As described in Viero et al. (2013b), in general there are correlations between redshift bins simply due to the presence of large-scale structure, and to take these into account one would have to fit for all of the amplitudes in the defined bins simultaneously.In our case, redshift bins have a width of 2, so these correlations are expected to be negligible.We therefore fit for the amplitude in each bin independently, which reduces to solving a simple weighted linear regression between the ALMA map and maps of 'hits' in the JADES catalogue for each bin.The solution is where  labels each map pixel,  denotes the redshift/stellar mass bin,   is the hits map of bin  convolved with the ALMA beam with the weighted mean subtracted,  is the data map with the weighted mean subtracted and   = 1/ 2  are the inverse-variance weights.The sum is performed over all the pixels in the map.Equation 3 is effectively the covariance between the map and the catalogue, but with the ALMA beam taken into account.
To estimate the uncertainties in Ŝ , 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  error in Ŝ .In order to visually show our results, we also evaluate Eq. 3 after shifting the hits map   relative to the data map  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.
The top panel of Figure 10 shows our results from the crosscorrelation, while in the bottom panel we list the central pixel values, which are the best-fit flux densities of the galaxies within each bin.For discussing the CIB contributions we want to know the pixel sum of the best-fit average surface brightness for each bin (weighted by the inverse-variance), which we calculate as ).The first (top left) panel shows the distribution of redshifts (spectroscopic when available, otherwise photometric).The next five panels show distributions of NIRCam magnitudes for the ALMA sources in the JADES catalogue.The final three panels show distributions of three selected NIRCam colours.The distributions from all 38/45 ALMA galaxies with JWST detections are shown in blue, while the distributions from the subset of eight new ALMA galaxies with JWST detections are shown in red.In the redshift panel (top left) there are only 35 ALMA galaxies with available redshifts, and five new ALMA galaxies with redshifts.In the bottom row we also only plot those galaxies with redshifts because the colours are quite uncertain.
Here  maj and  min are the beam major and minor axes (in standard deviation units), respectively, the quantity  is the solid angle of the map and  ,eff is the effective total number of sources from catalogue  in the map, weighted by the noise, (5) This is just the average number of sources per beam (weighted by the noise), times the number of beams in the map.In Fig. 10 we provide the values of  ,eff along with Ŝ .Note that when we calculate the weighted mean surface brightness of the sky from our model, the weighted mean should not be subtracted from    .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.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. 11 and we can see that it is consistent with noise.
In Fig. 10 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 ≳ 3  in all bins with stellar masses between 10 8.4 M ⊙ and 10 10.4 M ⊙ 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 5.3).
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 M ⊙ lying just below our 3.6  ALMA threshold, which have flux densities around 15 Jy; (2) there are more than 300 galaxies with stellar masses around 10 9 M ⊙ that can also be statistically detected in the ALMA map, with individual flux densities of a few Jy; and (3) galaxies with stellar masses around 10 8 M ⊙ or below have flux densities that are too low to be detected in the ALMA map, even statistically.

THE COSMIC INFRARED BACKGROUND IN THE HUDF
The absolute intensity of extragalactic light has been studied at all wavelengths from -rays to the radio (Hill et al. 2018).Peaking at around 200 m, the CIB is the brightness of the extragalactic sky at infrared wavelengths, averaged over the whole sky, after subtracting all contributions originating from the Solar System and the Milky Way.This absolute value tells us about the history of star formation and has been measured using the FIRAS instrument onboard the COBE satellite (Fixsen et al. 1998) the long-wavelength side of the CIB and we can try to use our new map to estimate what fraction can be accounted for in DSFGs.

Resolved source contribution to the cosmic infrared background
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;Dunlop et al. 2017;Hatsukade et al. 2018;González-López et al. 2020;Gómez-Guĳarro et al. 2022;Chen et al. 2023), with results typically in the tens of per cent range, depending on the precise wavelength.Here we explore this question with our new ALMA map at 1.23 mm and our new JADES catalogue.To start with, the sum of our detected source flux densities (> 3.6 ) is (7.33 ± 0.10) mJy, and these sources are detected within an area of 1.15 × 10 −3 deg 2 .Therefore we have resolved a total intensity of (6.35 ± 0.09) Jy deg −2 in individually-detected DSFGs.
In addition to this, our stacking analysis (Fig. 10) demonstrates that our map is also statistically sensitive to fainter galaxies.Multiplying Ŝ by  eff in each bin shown in Fig. 10, and summing, yields a flux density of (2.6 ± 0.5) mJy, or an intensity of (2.4 ± 0.5) Jy deg −2 (where the area used in our stack is 1.11 × 10 −3 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  1230 > 0.85 mJy present in our map, we could also add a contribution from brighter sources.The GOODS-ALMA survey (Gómez-Guĳarro et al. 2022) covered a much larger area with ALMA at 1 mm (about 72 arcmin 2 ) and found 22 galaxies with  1230 > 0.85 mJy, corresponding to an intensity of (1.31 ± 0.02) Jy 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 ± 0.5) Jy 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 ± 0.5) Jy deg −2 .We note that if we instead perform our .Top row: Distributions of three select NIRCam magnitudes for all ALMA-detected galaxies with JWST counterparts and with stellar masses above 10 10 M ⊙ , compared to all JWST galaxies with stellar masses above 10 10 M ⊙ that are not detected by ALMA.Bottom row: Same as top row but for three select NIRCam colours.While we do not see a significant difference between the two samples in terms of magnitude, the mm-bright galaxies detected by ALMA with high stellar masses tend to have redder NIRCam colours compared to ALMA-undetected galaxies with equally large stellar masses.
stacking analysis on the full map (without masking bright detected sources) we obtain a CIB estimate of (8.7 ± 0.8) Jy 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.

The absolute value of the cosmic infrared background
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  is a constant,  0 is a fiducial frequency,  = 0.65,   is the Planck function and  d = 18.5 K.The uncertainties in the three fit parameters ,  and  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  and  d (thus fixing the shape of the CIB spectrum), but take into account the uncertainties in the amplitude throughout our calculations.
To interpolate the CIB intensity to the frequency of our ALMA map, we first estimate what the transmission function of our ALMA map is, and use that to estimate the effective frequency of the image (which may be different from the mean frequency of 243 GHz).Each ALMA observation consists of two side bands 4 GHz wide, whose ) bins, after masking all of our detected ALMA sources.The cutouts are 10 arcsec × 10 arcsec.The colour represents S/N, where the signal is the variance-weighted mean and the noise is the uncertainty in the variance-weighted mean.Bottom: Central pixel values and uncertainties of the top panel (i.e.our best estimate of the average 1230-m flux density of galaxies within each bin), with the number of JADES galaxies contributing to each bin also given.For a definition of the quantity  eff , see Eq. 5.   ).We find that the resulting transmission function is dominated by ASPECS and is largely flat across Band 6, with slightly more weight towards lower frequencies.
With  () 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  () from Eq. 6 and the constants  and  0 cancel out in the ratio of integrals.We find a value of  eff = 241.4GHz, which is close to the ALMA data combined central frequency.
We finally estimate the amplitude of the CIB at  eff from Eq. 6, using two different approaches to obtain the amplitude .The first approach is to use  = (1.3±0.4) × 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 We show these estimates graphically in Fig. 12. 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. 12 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  1230 > 0.85 mJy galaxies and also the stackingestimate 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.

Consistency between resolved sources and the absolute value of the cosmic infrared background
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 / ≃ 15 per cent on scales around 10 arcmin (Viero et al. 2009(Viero et al. , 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 starformation-mode model of galaxy evolution.The total simulated area of the SIDES simulation is 117 deg 2 , but here we only use a 1 deg 2 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 2 tile, each with the area of our image of the HUDF (here 4.2 arcmin 2 ).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  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  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 2 patches of the SIDES simulation (after subtracting the contribution from  1230 > 0.85 mJy galaxies) is 11.1 Jy deg −2 , comparable to the value of (8.6 ± 0.5) Jy deg −2 that we have measured.The standard deviation of the SIDES simulation patches is 1.9 Jy deg −2 , thus fluctuations are / = 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 2 patch we have observed.
We now estimate the probability of measuring a CIB value of (10.0 ± 0.5) Jy deg −2 assuming that the true value is either 19 +6 −5 Jy deg −2 (Fixsen et al. 1998) or 23 +6 −8 Jy 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 ± 0.5) Jy 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  = 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 +0.4  −1.4  using the CIB value from Fixsen et al. (1998), or −2.9 +0.3  −1.2  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 (≃ 2 ) underdense direction on the sky. 6s 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.
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-IR 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.

IMPROVEMENT OVER PREVIOUS STUDIES
The deepest previously-published ALMA survey of the HUDF at 1 mm is ASPECS and so a question of interest is how much of an improvement have we achieved by including all of the additional data.Qualitatively, looking at  2017), ASAGAO and GOODS-ALMA overlap with the entire ASPECS map and so by combining them with ASPECS the result must be deeper.Additionally, by including more individual pointings, we have been able to make some regions even deeper still.
To quantify the difference between ASPECS and our combined image, we downloaded the ASPECS map produced using CASA in the MFS mode and made public by the ASPECS team, 7 then ran the primary-beam-corrected image through our algorithm for generating the noise map (see Section 2.3).The pixel size of the ASPECS map is 0.2 arcsec, the same as our -combined map, and the beamsizes are very similar (1.49 arcsec × 1.07 arcsec versus 1.53 arcsec × 1.08 arcsec), so we simply computed the ratio of the two noise maps to assess the amount by which the noise improves with the additional data.Unsurprisingly we find that across the entire ASPECS region the noise (meaning here the rms after masking sources) in our map is smaller, ranging from about 5 per cent smaller in the central region to about 50 per cent smaller near the edges and around the deepest individual pointings.The ASPECS map was made going out to a primary beam value of 0.1 and covers a total area of 4.2 arcmin 2 , 3.3 arcmin 2 of which has a noise level less than 35 Jy beam −1 , whereas our map was made going out to a primary beam value of 0.2 and all 4.2 arcmin 2 has a noise level less than 35 Jy beam −1 .Thus by combining the different data sets we not only reduce the noise across the majority of the map by 5 per cent, but we also expand the area where the map is at its most sensitive.González-López et al. (2020) searched the ASPECS map for sources down to a fidelity of 0.5, and their lowest significance source had a S/N of 3.3.Their catalogue contains a total of 35 sources (four of which are not detected in our map), whereas we find 45 sources; it is worth noting that all 45 of our sources are contained within the area mapped by ASPECS.Thus the extra depth we have been able to achieve by combining archival ALMA data with ASPECS has led to a significant increase in detected sources.
In addition to detecting more sources, we are able to recover more of the CIB through stacking thanks to our deeper JADES catalogue.If we stack on the 1661 galaxies from the CANDELS catalogue with stellar masses between 10 7.4 M ⊙ and 10 11.4 M ⊙ and with redshifts between 0 and 8 (after masking detected ALMA sources) we obtain (1.5 ± 0.4) Jy deg  2017) and the first GOODS-ALMA survey (Franco et al. 2018); their final map has a mean rms of about 75 Jy beam −1 , with a beamsize of 0.59 arcsec × 0.53 arcsec and a pixel size of 0.1 arcsec after applying a 250 k taper.Our shallow and wide combined map (presented in Appendix A) contains additional GOODS-ALMA data that were not available when the ASAGAO map was constructed, as well as more individual pointings, while maintaining a similar beamsize (0.87 arcsec × 0.64 arcsec) and pixel size (0.12 arcsec).We thus also downloaded the ASAGAO 250 k-tapered map8 to quantitatively check the improvement.After running the ASAGAO map through our noise algorithm, we find that the new GOODS-ALMA data reduces the noise by 10-20 per cent throughout, while some individual pointings go about 50 per cent deeper.

CONCLUSIONS
We have produced a series of 1-mm maps of the HUDF by combining all of the previously-published survey data in the  plane in various ways, reducing the noise compared to previous studies.We specifically constructed a deep map covering 4.2 arcmin 2 and a shallower map covering 25.4 arcmin 2 .Our deep map has a pixel rms that ranges from 5 to 50 per cent lower than in the best previous study, with an area of about 1.5 arcmin reaching below 9 Jy beam −1 and a minimum of about 4.6 Jy beam −1 reached in some regions.Our shallow map has a pixel rms that ranges from 10 to 50 per cent lower than in the best previous wider and shallower study.We make all of our maps publicly available 9We searched our deep map for sources down to a signal-to-noise threshold of 3.6, finding a total of 45 peaks in the S/N map, 13 of which are new.Nearly all (38/45) of these ALMA sources have near-IR counterparts detected by JWST and HST.We additionally find 27 sources in our wider map, nine of which are new.The JADES data enable stellar masses and photometric redshifts to be estimated for the ALMA source counterparts, and we find that they are all relatively high  * galaxies.Compared with ALMA-undetected galaxies at similar  * , the ALMA-detected galaxies typically have redder JWST colours.With our larger sample of mm-selected sources in the HUDF, other studies investigating the statistical properties of the faintest starforming galaxies could be carried out.For example, Aravena et al. (2020) studied the SFR-stellar mass relation of ASPECS galaxies and Boogaard et al. (2023) looked at the morphologies of ASPECS galaxies in JWST's MIRI filters; expanding to a larger sample size and using improved SED fits from JWST photometry could lead to more robust conclusions.
Since the vast majority of near-IR-selected galaxies are not directly detected in our 1-mm map, we performed a stacking analysis on their positions.We found significant average signals for all galaxies in the range  = 0 to 3 and with stellar masses between 10 9.4 M ⊙ and 10 10.4 M ⊙ , as well as a roughly 3  signal from 10 8.4 M ⊙ to 10 9.4 M ⊙ stellar mass galaxies between  = 0 and 2. The evidence is that  * ∼ 10 10 M ⊙ ALMA-undetected galaxies have a 1-mm flux density around 15 Jy and would be individually detected in even deeper integrations.The stacking results also show the value of our new data products for performing similar statistical analyses on sets of galaxies detected in other wavebands.
We used our galaxy detections, as well as our stacking analysis, to estimate the level of the CIB at 1.23 mm that we have resolved.We thus account for a background level of (10.0 ± 0.5) Jy deg −2 , with an expectation that even fainter galaxies will hardly change this number.There is still a large uncertainty in the total background estimate, and we also stress that the variance expected in a region as small as the HUDF is around 20 per cent.We do not recover all of the background, and there are a number of possible resolutions: (1) HUDF is a 2-3  negative fluctuation in the CIB; (2) the absolute level of the CIB has been overestimated; (3) there is a new population of galaxies that contribute at 1 mm, but are not detected by JWST; or (4) a genuinely diffuse component of the background exists.Given that the simple explanation (1) cannot be excluded, then the suggestion is that deep ALMA+JWST observations may account for all of the mm-wave background.
The HUDF is probably the best-studied extragalactic field, and it is crucial to probe this field at longer (> 500 m) wavelengths in order to understand how the earliest galaxies began forming their stars, complementing the data obtained by HST, JWST and other telescopes at shorter wavelengths.The 4.2 arcmin 2 map presented here is the deepest such image of the HUDF available to-date, yet it is clear that there are still many more galaxies left to detect.Deeper ALMA observations of the HUDF will inevitably uncover these galaxies, providing a more complete understanding of galaxy evolution at early times.a Approximate map rms estimated using the observatory MFS data products.For programmes with a single tuning, this is the rms of the primary-beam-uncorrected map after masking all known sources.For programmes with multiple tunings, we follow the same procedure and then estimate the rms of the weighted mean of the images as 1/ Σ  1/ 2  .The actual rms from combining images at different tunings in  space using the MFS mode will generally be less than this estimate.b Maximum synthesized beam across all frequencies of a given data set.The actual synthesized beam from combining data taken in different tunings using the MFS mode may differ slightly from this estimate.c The sky coverage overlapping with our data selection criteria, which may be less than the total sky coverage of the given programme.d Programmes used to make a mosaic with uniform noise.Table A2.Positions and flux densities ( 1230 ) for the 27 sources found in the ALMA 1-mm combined image.Here flux densities are measured by fitting Gaussian profiles to the sources in order to estimate the number of beams per source, then scaling the peak pixel value by the number of beams.For sources consistent with 1 beam or less and for sources with peak pixels > 1000 Jy we set the flux density to be the peak pixel value.For the extended galaxy ALMA-HUDF-69 (marked with a dagger), we use an aperture with a radius of 2 arcsec to calculate a flux density.For galaxies with matches in the JADES catalogue and sufficient photometry to fit SEDs we provide stellar masses and photometric redshifts from McLeod et al. (in prep.).Galaxies with spectroscopic redshifts from Rieke & the JADES Collaboration ( 2023) are indicated by asterisks, which we provide instead of photometric redshifts.The median photometric redshift uncertainty of our ALMA sample is Δ/(1 + ) ≃ 0.03, while the relative uncertainties in stellar masses are small.

Figure 1 .
Figure1.Left: Signal map of the deepest region of the HUDF after combining all of the available archival ALMA Band-6 data given in Table1.This region is defined by the contour where the primary beam reaches 0.2, covering a total of 4.2 arcmin 2 .Middle: Corresponding noise map.Right: S/N map, from dividing the signal map by the noise map.

Figure 2 .Figure 3 .
Figure 2. Results from combining all of the data in Table 1 in the  plane before imaging, and imaging the same data individually before combining.Signal maps are shown in the two top left panels, noise maps are shown in the two top right panels and S/N maps are shown in the bottom two panels.In each case, the  combination is on the left and the real-space image combination is on the right.Sources with S/N > 4 are marked with green squares in the S/N maps.

Figure 4 .
Figure 4. Purity, i.e. 1 minus the the total number of positive to negative peaks greater than a given S/N map, plotted as a function of S/N.A purity of 0.7 corresponds to a S/N threshold of 3.6, which we use to extract sources from our deep map (Fig.1).

Figure 6 .
Figure 6.Signal-to-noise ratio map of the deep central region, covering 4.2 arcmin 2 , with peaks > 3.6  indicated as green boxes.Galaxies found by Dunlop et al. (2017) are shown as yellow circles, galaxies from the ASPECS survey (González-López et al. 2020) are orange circles and galaxies from the ASAGAO survey (Hatsukade et al. 2018) are red circles.A few of the brightest galaxies are also detected in the GOODS-ALMA survey (Gómez-Guĳarro et al. 2022) andshown as brown circles.The gold contour shows the footprint of the deepest HST data from the CANDELS/HUDF09 survey, while the entire region is covered by the JADES survey.Galaxies with counterparts in both the CANDELS catalogue and the JADES catalogue are indicated with a blue cross and galaxies with only a JADES counterpart are indicated with a red cross.There are no galaxies with a CANDELS/HUDF09 counterpart but no JADES counterpart.Numbers refer to the labels in Tables2 and 3.

Figure 7 .
Figure7.Distribution of ALMA source properties matched to our JADES catalogue(Rieke & the JADES Collaboration 2023, McLeod et al. in prep.).The first (top left) panel shows the distribution of redshifts (spectroscopic when available, otherwise photometric).The next five panels show distributions of NIRCam magnitudes for the ALMA sources in the JADES catalogue.The final three panels show distributions of three selected NIRCam colours.The distributions from all 38/45 ALMA galaxies with JWST detections are shown in blue, while the distributions from the subset of eight new ALMA galaxies with JWST detections are shown in red.In the redshift panel (top left) there are only 35 ALMA galaxies with available redshifts, and five new ALMA galaxies with redshifts.In the bottom row we also only plot those galaxies with redshifts because the colours are quite uncertain.
Figure9.Top row: Distributions of three select NIRCam magnitudes for all ALMA-detected galaxies with JWST counterparts and with stellar masses above 10 10 M ⊙ , compared to all JWST galaxies with stellar masses above 10 10 M ⊙ that are not detected by ALMA.Bottom row: Same as top row but for three select NIRCam colours.While we do not see a significant difference between the two samples in terms of magnitude, the mm-bright galaxies detected by ALMA with high stellar masses tend to have redder NIRCam colours compared to ALMA-undetected galaxies with equally large stellar masses.

Figure 10 .
Figure10.Top: Stacks at the positions of JADES galaxies in redshift and stellar mass(from McLeod et al. in prep.)bins, after masking all of our detected ALMA sources.The cutouts are 10 arcsec × 10 arcsec.The colour represents S/N, where the signal is the variance-weighted mean and the noise is the uncertainty in the variance-weighted mean.Bottom: Central pixel values and uncertainties of the top panel (i.e.our best estimate of the average 1230-m flux density of galaxies within each bin), with the number of JADES galaxies contributing to each bin also given.For a definition of the quantity  eff , see Eq. 5.

Figure 11 .SFigure 12 .
Figure 11.Stack at the positions of 21 JADES sources flagged as stars that happen to lie within our ALMA 1-mm map; as expected, we do not find a significant signal.
central frequencies are separated by 12 GHz, and for each observation we have already computed the central rms of the observatoryproduced MFS image (  ; see Section 3).Assuming each observation's individual transmission function is flat (in ) across the two sidebands, we can calculate the mean transmission function of all of the relevant observations used to produce our final map, weighted by each observation's inverse variance.To each weight we must also include a factor of the ratio of the given observation's beam area (set by  maj, and  min, ) to the final image's beam area (set by  maj and  min ), thus the transmission function weight for observation  is  = 1  2   maj,  min,  maj  min .(7)Toobtain a transmission function,  (), appropriate for the average of the entire map, we only consider the observations from surveys that uniformly covered the HUDF, namely programmes 2012.1.00173.S (from Dunlop et al. 2017), 2015.1.00098.S (ASAGAO), 2015.1.00543.S and 2017.1.00755.S (GOODS-ALMA) and 2016.1.00324.L (ASPECS find absolute CIB values of 19 +6 −5 Jy deg −2 using the fit from Fixsen et al. (1998) and 22 +7 −6 Jy 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.

Figure A1 .Figure A3 .
Figure A1.Left: Signal map of the shallower, larger region of the HUDF after combining all of the available archival ALMA Band-6 data given in TableA1.This region is defined by the contour where the primary beam reaches 0.4, covering a total of 25.4 arcmin 2 .Middle: Corresponding noise map.Right: S/N map, made by dividing the signal map by the noise map.In all panels, the cyan contour shows the footprint of our combined deeper central map, which we exclude from our source extraction procedure in the shallower region.

Figure
Figure B1. 5 arcsec × 5 arcsec cutouts.The background shows JWST F356W imaging, while the contours are from our deep ALMA map in steps of 2  , with  = 1, 2, 3, 4, 5 and 6.ALMA centroids are indicated by magenta circles and the positions of JWST counterparts are indicated by green crosses.ALMA IDs 28, 39, and 44 all have a F356W flux density signal-to-noise ratios less than 5, and we do not use their photometry to fit SEDs.

Table 1 .
ALMA projects downloaded from the ALMA archive and used to create the final combined 1-mm image.
. The 1-mm background lies on Figure 8. Stellar mass versus redshift for JADES galaxies (McLeod et al. in prep.), with JADES galaxies detected in our  1230 image highlighted in red (circles indicate photometric redshifts, stars indicate spectroscopic redshifts).The apparent vertical features are a result of the grid used for the photometric redshift fitting.
Table 1 we can see that ASPECS is by far the deepest of the individual surveys.The maps from Dunlop et al. (

Table A1 .
Archival ALMA projects downloaded from the ALMA archive and used to create the larger and shallower combined 1-mm image.