A Full Resolution of the 450 μm Extragalactic Background Light

The extragalactic background light (EBL) is the cumulative radiation outside the Milky Way. The determination of its corresponding primary emitting sources as well as its total energy level across the entire electromagnetic spectrum has profound implications for both cosmology and galaxy formation. However, the detailed origin of the EBL at far-infrared wavelengths, particularly those close to the peak of the cosmic infrared background, remains unclear. Here we report the results of our ongoing SCUBA-2 450 μm survey of 10 massive galaxy cluster fields. By exploiting the strong gravitational lensing offered by these clusters, we obtain significant counts down to an unprecedented depth of ∼0.1 mJy at this wavelength, about 10 times deeper than that reached by any other previous survey. The cumulative energy density based on the counts is 138.1 −19.3+23.9 Jy deg−2 or 0.45 −0.06+0.08 MJy sr−1. Comparing our measurements to those made by the COBE and Planck satellites, we find that at this flux density level, the 450 μm EBL is entirely resolved by our SCUBA-2 observations. Thus, we find for the first time that discrete sources produce fully to the 450 μm EBL, and that about half of it comes from sources with sub-mJy flux densities. Our deep number counts provide strong constraints on galaxy formation models.


Introduction
The cosmic energy budget, encompassing the total amount of energy radiated throughout the Universe, is a fundamental aspect of modern astrophysics.Central to this budget is the extragalactic background light (EBL), which represents the integrated emission across the electromagnetic spectrum from astrophysical sources outside the Milky Way (see, e.g., the review by Cooray 2016).Observational measurements of the EBL at different wavelengths allow one to understand its energy distribution, and, in principle, can provide insights into the dominant contributors to the EBL (see, e.g., the discussion in Hill et al. 2018).
One common way to put constraints on the EBL is to perform imaging surveys using ground-based facilities, which allow one to construct source number counts and calculate the integrated energy densities (e.g., the early results in the farinfrared (FIR)/submillimeter by Smail et al. 1997;Barger et al. 1999;Cowie et al. 2002).By comparing with satellite measurements with lower spatial resolution, which, in principle, take into account any diffuse emission, one can estimate the energy contributions to the EBL from galaxies and understand whether galaxies are the dominant contributors.
Surveys have found that the optical and near-infrared EBL primarily originate from the directly observed star formation in galaxies, while the FIR/submillimeter EBL mainly comes from the thermal emission from interstellar dust reradiated starlight (Puget et al. 1996;Fixsen et al. 1998;Dole et al. 2006).These dusty galaxies, sometimes called submillimeter galaxies (SMGs), are typically characterized by intense dust emission indicating high rates of star formation (see, e.g., review articles by Blain et al. 2002;Casey et al. 2014).
The advent of space-based telescopes, like the Herschel Space Observatory, and ground-based facilities, like the Atacama Large Millimeter/submillimeter Array (ALMA), has revolutionized our ability to perform either wide-field or deep FIR surveys, enabling unprecedented studies of galaxy number counts (Oliver et al. 2010;Berta et al. 2011;Gómez-Guijarro et al. 2022;Chen et al. 2023;Cowie et al. 2023).These observations, coupled with theoretical models, provide a comprehensive view of the FIR Universe and its connection to other astrophysical phenomena (e.g., Lacey et al. 2016;Lagos et al. 2019;Hayward et al. 2021).However, due to the confusion limits of Herschel and the inefficient survey capability of ALMA with its small field-of-view, the number counts at wavelengths close to the peak of the FIR EBL (∼200-300 μm; e.g., Odegard et al. 2019) remain limited to the brightest end (e.g., Oliver et al. 2010).
SCUBA-2, a state-of-the-art submillimeter camera, offers unmatched sensitivity and high angular resolution at 450 μm.Its optimized design for deep submillimeter observations makes it an excellent tool for studying dust-rich galaxies (Holland et al. 2013).By harnessing the sensitivity of SCUBA-2 at submillimeter wavelengths, we can detect and characterize faint sources at 450 μm, enabling comprehensive studies of galaxy number counts close to the peak of the FIR EBL and a deeper exploration of the submillimeter Universe (Casey et al. 2013;Geach et al. 2013;Wang et al. 2017;Zavala et al. 2017;Lim et al. 2020;Barger et al. 2022;Gao et al. 2024).
Additionally, the technique of gravitational lensing provides a unique means for studying galaxy number counts at the fainter end (Smail et al. 1997;Cowie et al. 2002Cowie et al. , 2022;; Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.Knudsen et al. 2008;Johansson et al. 2011;Chen et al. 2013;Hsu et al. 2016).Gravitational lensing occurs when the gravitational field of a massive object, such as a rich cluster of galaxies, bends and magnifies the light from distant galaxies.By exploiting the lensing effect, we can effectively boost the observed flux of background galaxies, allowing us to probe deeper into the Universe and uncover fainter sources that would otherwise remain undetected in blank-field surveys.
In this study, we utilize the deep data obtained with SCUBA-2 on 10 massive galaxy cluster fields to derive robust number counts at 450 μm.Our analysis properly accounts for selection biases, completeness, and uncertainties.By constructing the deepest 450 μm number counts ever, we aim to unravel the relative contributions of SMGs to the 450 μm EBL.In Section 2, we describe our SCUBA-2 data and data reduction.In Section 3, we describe our methodology, including source extraction, Monte Carlo simulations, and number count calculations.In Section 4, we present our number counts and the integrated energy density.We summarize our results in Section 5.

Data and Data Reduction
We retrieved the SCUBA-2 450 μm data from the CADC archive.We used the data taken under weather band 1 and band 2 conditions (τ 225 GHz < 0.08) between 2011 October and 2022 July (Cowie et al. 2022).The scan pattern used for these 10 cluster fields was CV DAISY, which has a roughly circular field size of ¢ 6  in radius.We summarize the data in Table 1.We reduced the data following Chen et al. (2013).We used the Dynamic Iterative Map-Maker (DIMM) method in the SMURF package contained in the STARLINK software (Chapin et al. 2013).This method models individual components that make up the time series recorded by the bolometer to produce science maps.We adopted the "blank field" configuration file, which is suitable for detecting faint point sources in extragalactic surveys.For each cluster field, we produced scan maps with a pixel scale of 1″ and then applied the recommended flux conversion factors (FCFs) from Mairs et al. (2021) to convert the pixel unit from picowatts to Jansky per beam.After calibration, we used the MOSAIC_JCMT_I-MAGES recipe from the Pipeline for Combing and Analyzing Reduced Data (PICARD) to coadd and mosaic calibrated scans for each field.
To improve source detection, we applied a matched filter using the SCUBA2_MATCHED_FILTER recipe in PICARD.This recipe first convolves the map with a Gaussian to estimate the low spatial frequency noise and then subtracts it from the original map.We adopted the default 20″ FWHM value for the Gaussian profile.We verified the flux recovery capability of SCUBA2_MATCHED_FILTER following Lim et al. (2020), and we adopted a mean upward correction of 5.3% for the flux loss to our 450 μm data.

Source Extraction
Before source extraction, for each field, we generated SCUBA-2 point-spread function (PSF) models by stacking the 10-20 highest signal-to-noise ratio (S/N) source images without neighboring sources.This method inherently assumes that the typical source size is much smaller than the beam size of SCUBA-2 in the submillimeter (∼7 5 at 450 μm), which is supported by recent ALMA observations that found typical sizes of subarcseconds (e.g., Simpson et al. 2015;Hodge et al. 2016;Fujimoto et al. 2017;Gullberg et al. 2018;Tadaki et al. 2020).We then fitted a double Gaussian profile to model the PSFs, and we used the best-fit model PSFs for source extraction.
We performed the source extractions as in Hsu et al. (2016).We searched for the maximum S/N pixel in the central circular region (6′ radius).(The pointing centers of the maps are normally close to the cluster centroids.)We recorded the location and flux density of the pixel, then subtracted a rescaled PSF centered at this pixel and searched for the next maximum S/N pixel.Following Hsu et al. (2016), we used a 3σ threshold, which allows us to obtain a better S/N in the number counts at the faint end.We repeated the extraction process until we reached this threshold.In Figure 1, we show one of our cluster fields as an example for the source extraction.
We determined the magnification factors μ i of the point sources from the magnification maps.The demagnified flux density of each source can be obtained from where S demag,i and S obs,i are demagnified and observed flux densities, respectively.We calculated the effective area A eff,i of each source on the source plane.We summed over the pixels whose S/Ns were greater than 3σ and then converted pixels to square degrees.We then calculated the delensed raw counts at the j-th flux bin as Here, X i represents the number density contribution of each source within that flux bin.We based the error calculation on Poisson statistics.

Simulations
Corrections for flux boosting, false detections, and incompleteness are needed in order to obtain the intrinsic number counts.To do this, we ran Monte Carlo simulations to find the underlying models for our fields.We used the Schechter function form as our number counts model: We generated artificial sources whose flux densities were assigned according to the underlying models for each cluster field.We then randomly distributed these sources in the source plane.Next, we projected the simulated sources onto the image plane using LENSTOOL.The outputs of LENSTOOL contain the new flux densities and positions of the simulated sources in the image plane.We convolved the simulated sources with the PSF and added them into the jackknife maps to produce mock observation maps.We produced jackknife maps following Hsu et al. (2016) by coadding even and odd scans separately.We then subtracted these two coadded maps and rescaled the value of each pixel by a factor of t t even odd /(t even + t odd ), where t even and t odd represent the integration times of each pixel from the even and odd coadded maps, respectively.We also applied matched filtering to the jackknife maps, as we did for the real data images.
In order to estimate the true number counts, we adopted an iterative procedure in our simulations.We generated 15-20 mock maps in each iteration step.We then performed source extraction and computed the averaged recovered counts of the mock maps.Next, we corrected the raw counts by using the ratio between the averaged recovered counts and the raw counts.Finally, we did a χ 2 fit to the corrected counts by using the Schechter function.This fit will be the new counts model for the next iteration.We iterated until the recovered counts converged with the raw counts to within the 1σ uncertainties for all the flux bins.Once we obtained the intrinsic number count models, we used these to produce 500 mock images and then performed the source extractions to create the source catalogs.Following Gao et al. (2024), we cross-matched the input and output catalogs within a half-beam FWHM as our search area to find the brightest counterparts.We considered S output /S input 3 a match in this analysis, which is similar to what has been adopted previously (e.g., Geach et al. 2017).Following Gao et al. (2024), we estimated the boosting factors, false detection rates, and completeness of the point sources using a two-dimensional binning method in our cross-matched catalogs.

Delensed Corrected Number Counts
To derive delensed corrected number counts, we first deboosted the fluxes of the sources.We then delensed their fluxes and estimated their effective areas in the source plane.where p false,i is the false detection rate, and C i is the completeness.We confirmed that the delensed corrected number counts agree with the intrinsic count models.We estimated the statistical uncertainties from the source positions and took these uncertainties into account in the analyses.(2024).The horizontal red dashed line is the EBL measured by the COBE +Planck satellites (Odegard et al. 2019).The red shaded region shows the range from COBE estimates (e.g., Puget et al. 1996;Fixsen et al. 1998;Gispert et al. 2000).

Systematic Uncertainties
Systematic uncertainties related to the modeling of gravitational lensing need to be taken into account for a proper assessment of the number counts error budget.In the following, we give our results with estimates of the systematic uncertainties caused by the redshift distribution of the background-lensed SMGs, the lens models, and the clustering of the source distributions.
In our methodology, we calculated corrected number counts by assuming a median redshift of 1.5.However, it is expected that the SMGs have a redshift distribution, which could affect the magnification estimates.This leads to additional uncertainties in the number counts.To address this, we randomly assigned redshifts to our sources using the redshift distribution from the STUDIES survey, which is the deepest 450 μm blankfield survey (Wang et al. 2017;Dudzevičiūtė et al. 2021;Gao et al. 2024).We calculated the corrected number counts again and compared the standard deviations of the number counts in the randomized redshift sample with the Poisson noise.We found that the uncertainties caused by the assumptions concerning the redshift are subdominant, and they are, on average, about 25% of the Poisson uncertainties.Nevertheless, we include this error budget in the total error budget.
To estimate the systematic uncertainties caused by the lens modeling, we utilized the various lens models provided for the five Hubble Frontier Fields (A370, A2744, MACS J0416.1-2403,MACS J0717.5+3745, and MACS J1149.5+2223).We ran the same Monte Carlo simulations as those done on the real data but using different lens models.We then estimated the systematic uncertainties by calculating the standard deviations of the corrected counts obtained using the different lens models.We found the systematic uncertainties to be subdominant, again about 25% of the Poisson uncertainties.We include this error budget in the total error budget.
We note that in our Monte Carlo simulations, we randomly distributed the positions of the artificial sources.This might bias our counts since we do not take clustering effects into account when we calculate the boosting factors, false detection rates, and completeness from the mock maps.To test whether neglecting clustering could significantly alter our count results, we used the empirical catalogs produced by the SIDES simulation (Béthermin et al. 2017), which inherently include clustering effects, since the simulation builds upon dark matter lightcones.
Before assessing the clustering effects, we first validated our methodology for estimating intrinsic counts by performing the following test.We clipped the original SIDES 2 deg 2 simulated map into a set of 25 smaller cutout maps with a size similar to our SCUBA-2 footprint, and we treated them as different cluster fields in the source planes.We then lensed these simulated maps to the image planes using LENSTOOL and used our methodology to find the corrected counts.We compared the corrected counts to the true counts provided by SIDES.In this test, we adopted the three-dimensional source positions from the SIDES catalog instead of having the sources randomly distributed on the sky with an assigned redshift.We show our results in Figure 2, where the averaged corrected number counts over the 25 cutout maps are consistent with the true number counts.
After validating our methodology, we moved on with a similar test.This time we randomly distributed the source positions to calculate the corrected number counts using the smaller SIDES cutout maps.We show our results in Figure 2.While not taking clustering into account can lead to larger uncertainties, evidently with a larger scatter in the averaged counts, on average, there is no significant difference between the SIDES counts and the corrected counts.Our results are therefore consistent with previous studies-either from other SCUBA-2 450 μm surveys (Wang et al. 2017) or from Herschel studies with similar beam sizes but at slightly shorter wavelengths (Béthermin et al. 2015)-where no significant impact of clustering on the number counts was found.

Results and Discussion
We present our corrected differential number counts for the 10 lensing cluster fields in the left panel of Figure 3.The solid black curve is the best-fit Schechter function for our corrected counts, which can be parameterized by Equation (4) with the following parameters: N 0 = 4437.5 ± 1399.5, S 0 = 10.4 ± 1.7, and α = −1.9± 0.1.The uncertainties on each data point include Poisson noise and the uncertainties from the redshift distribution and the lens model.Thanks to the powerful effects of strong gravitational lensing, we have pushed the detection limit down to ∼0.1 mJy at 450 μm, a factor of >10 improvement over the deepest blank-field counts (Wang et al. 2017;Gao et al. 2024).
We also calculate the weighted average counts by combining all 10 fields.We show the results in the right panel of Figure 3, and we provide the corresponding values in Table 2. Our results are consistent with the previous measurements shown in Figure 3.However, the various physical or empirical models (dotted and dashed curves in the right panel) tend to overpredict source densities at the bright end at 450 μm (1 mJy) by about 10%-30%.On the other hand, at the faint end (1 mJy), our counts are slightly higher, although not significantly.The physical reasons for the disagreement between the measurements and the models at the bright end need to be investigated.Gao et al. (2024) proposed that it could be due to the mismatch in halo masses between those inferred from clustering measurements (Lim et al. 2020) and those in the models.Studies of the physical properties of the 450 μm sources, such as stellar mass, could potentially shed more light on this issue.
To estimate the contributions of the detected sources to the 450 μm EBL, we integrated the best-fit Schechter model, which we show as a function of flux in Figure 4. Above 0.1 mJy where significant counts are obtained, we find the total energy density at 450 μm to be 138. of the total 450 μm EBL reported by Odegard et al. (2019).
Past works have suggested that a broken power law could be a viable alternative to the Schechter function for the underlying count models.It can be described as We ran further counts analyses based on this model form and found consistent results compared to those obtained based on Schechter functions.The best-fit parameters of the broken power law are N 0 = 151.0± 14.6, S 0 = 9.6 ± 0.4, α = 2.0 ± 0.1, and β = 5.6 ± 1.1.We plot the cumulative energy density of the best-fit broken power law in Figure 4.The cumulative energy density based on the broken power law is -+ 126.9 41.6 41.6 Jy deg 2 .Our work demonstrates for the first time that discrete sources are the dominant contributors to the 450 μm EBL.Interestingly, about half of the contribution comes from sources that are fainter than ∼1 mJy, below the typical confusion limit of SCUBA-2 450 μm images (Gao et al. 2024).Noticeably, as shown in Figure 4, the integrated energy density does not converge when integrating down to 0.14 mJy, suggesting that deeper data are needed in order to put a tighter constraint on the faint end counts and thus obtain a converged constraint on the EBL contributions from discrete sources.

Summary
In summary, our research represents a significant step forward in the study of the 450 μm EBL.Through the innovative combination of gravitational lensing and the unparalleled capability of SCUBA-2, we have achieved a complete resolution of the 450 μm EBL and established the dominance of discrete SMGs as its main contributors.These findings provide a broader understanding of the FIR/ submillimeter regime, the cosmic energy distribution, and the interplay between galaxies and the diffuse background radiation.Our measurements could also be helpful for the design of the next-generation submillimeter facilities, such as Large Submillimeter Telescope (Kohno et al. 2020) and AtLAST (Klaassen et al. 2020), which aim to obtain wide blank-field images to a depth that is similar to what has been reached by this work.

Figure 1 .
Figure 1.Left: an example SCUBA-2 450 μm flux density map of the cluster MACS J1149.5+2223 with a ¢ 6 radius circular footprint.The red circles on the 450 μm map show the extracted >3σ sources.The orange boxes mark the sources below 1 mJy after delensing.The white contours represent the CATS lens model for z = 1.5 at magnification values of 1.5, 2.0, and 2.5 (moving inwards).Right: histogram of the S/N values based on the map of MACS J1149.5+2223.The orange region shows the S/N distribution in the jackknife map.The blue represents detections in the data map.

Figure 2 .
Figure 2. Ratios between the corrected number counts and the model counts, where the model counts are based on the SIDES simulation source catalog, and the corrected counts were obtained by applying our methodology to the mock images that were made based on the SIDES source catalog.The top panel shows the results for the case where the source positions are adopted from SIDES, which includes clustering.The bottom panel shows the results for the case where the sky positions of the sources are randomized.

Figure 3 .
Figure 3. Left: corrected differential number counts for all 10 cluster fields at 450 μm.The solid black curve shows the best-fit Schechter function, while the shaded gray region is the 68% confidence interval.The uncertainties include Poisson noise, as well as those caused by the underlying redshift distributions and lens models.Right: weighted average counts of all 10 cluster fields in filled black circles, along with previous 450 μm surveys in open symbols (Chen et al. 2013; Hsu et al. 2016; Wang et al. 2017; Zavala et al. 2017; Gao et al. 2024).Model predictions from SIDES (Béthermin et al. 2017), GALFORM (Cowley et al. 2019), and SHARK (Lagos et al. 2020) are shown as dotted and dashed curves.
magnification.b Total number of sources in each flux bin.

Figure 4 .
Figure 4. Cumulative EBL as a function of flux density at 450 μm.The black solid curve was calculated based on our best-fit Schechter function with the 1σ uncertainties in gray shading.The blue solid curve was calculated based on a broken power-law function.The black dotted curve was calculated by combining our results with the deepest blank-field counts from Gao et al.(2024).The horizontal red dashed line is the EBL measured by the COBE +Planck satellites(Odegard et al. 2019).The red shaded region shows the range from COBE estimates (e.g.,Puget et al. 1996;Fixsen et al. 1998;Gispert et al. 2000).

Figure 5 .
Figure 5. Flux density maps for the remaining nine fields.The format follows that adopted in the left panel of Figure 1.

Table 2
Combined Differential Number Counts at 450 μm