GAMA/DEVILS: Constraining the cosmic star-formation history from improved measurements of the 0.3-2.2 micron Extragalactic Background Light
Soheil Koushan, Simon P. Driver, Sabine Bellstedt, Luke J. Davies, Aaron S. G. Robotham, Claudia del P Lagos, Abdolhosein Hashemizadeh, Danail Obreschkow, Jessica E. Thorne, Malcolm Bremer, B.W. Holwerda, Andrew M. Hopkins, Matt J. Jarvis, Malgorzata Siudek, Rogier A. Windhorst
MMNRAS , 1– ?? (2021) Preprint 25 February 2021 Compiled using MNRAS L A TEX style file v3.0
GAMA/DEVILS: Constraining the cosmic star-formationhistory from improved measurements of the 0.3-2.2 µ mExtragalactic Background Light Soheil Koushan (cid:63) , Simon P. Driver , Sabine Bellstedt , Luke J. Davies ,Aaron S. G. Robotham , Claudia del P Lagos , , Abdolhosein Hashemizadeh ,Danail Obreschkow , Jessica E. Thorne , Malcolm Bremer , B.W. Holwerda ,Andrew M. Hopkins , Matt J. Jarvis , Malgorzata Siudek , , andRogier A. Windhorst ICRAR, The University of Western Australia, 35 Stirling Highway, Crawley WA 6009, Australia ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D) Astrophysics Group, School of Physics, University of Bristol, Tyndall Avenue, Bristol, BS8 1TL, UK Department of Physics and Astronomy, 102 Natural Science Building, University of Louisville, Louisville KY 40292, USA Australian Astronomical Optics, Macquarie University, 105 Delhi Rd, North Ryde, NSW 2113, Australia Oxford Astrophysics, Denys Wilkinson Building, Keble Road, Oxford OX1 3RH, UK Institut de Fisica d’Altes Energies (IFAE), The Barcelona Institute of Science and Technology, 08193 Bellaterra (Barcelona), Spain National Centre for Nuclear Research, ul. Hoza 69, 00-681 Warsaw, Poland School of Earth & Space Exploration, Arizona State University, Tempe, AZ85287-1404, USA
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
We present a revised measurement of the optical extragalactic background light (EBL),based on the contribution of resolved galaxies to the integrated galaxy light (IGL).The cosmic optical background radiation (COB), encodes the light generated by star-formation, and provides a wealth of information about the cosmic star formationhistory (CSFH). We combine wide and deep galaxy number counts from the GalaxyAnd Mass Assembly survey (GAMA) and Deep Extragalactic VIsible Legacy Survey(DEVILS), along with the Hubble Space Telescope (HST) archive and other deepsurvey datasets, in 9 multi-wavelength filters to measure the COB in the range from0.35 µ m to 2.2 µ m. We derive the luminosity density in each band independently andshow good agreement with recent and complementary estimates of the optical-EBLfrom very high-energy (VHE) experiments. Our error analysis suggests that the IGLand γ -ray measurements are now fully consistent to within ∼ . − .
076 M (cid:74) yr − Mpc − at a lookback time of 9.1 to 10.9 Gyrs. Key words: cosmology: cosmic background radiation - cosmological parameters -diffuse radiation – galaxies: statistics - evolution – methods: data analysis (cid:63)
E-mail: [email protected]
The extragalactic background light (EBL) is defined asthe total flux received today (nWm − steradian − ) fromall sources of photon production since the epoch of recom- © a r X i v : . [ a s t r o - ph . C O ] F e b S. Koushan et. al. bination plus the Cosmic Microwave Background (CMB).The CMB component is somewhat distinct, as it representsthe radiation field from the epoch of recombination asthe temperature, pressure, and density of the Universedropped to the point where nuclei and electrons combinedto form neutral atoms (Hauser & Dwek 2001; Kashlinsky2006; Cooray et al. 2009; Grazian et al. 2009; Dom´ınguezet al. 2011; Inoue et al. 2013; Driver et al. 2016b). Hence,the EBL represents the integrated luminosity cone ofthe universe from the CMB surface to the present day.Understanding this energy is essential as it encodes allphoton production pathways. This photon-energy pro-duction pervades the entire electromagnetic spectrum,but is typically broken into distinct wavelength ‘windows’(see Hill et al. 2018), aligned with the detection technologies.Overall, the CMB, the redshifted thermal remnantof the hot early universe (Penzias & Wilson 1965), is thedominant component of the EBL. The subsequent photonproduction represents the remainder of the EBL, and hasbeen produced by stars, gas, accreting black holes (BHs)and dust emission — and other processes — arising overthe entire path length of the Universe (i.e. since ∼ ∼
40 timeshigher ( ∼
960 nWm − sr − ) than either the infrared oroptical backgrounds (Scott 2000; Wright 2004; Dole et al.2006; Driver et al. 2016b; Driver 2020).The cosmic optical background (COB) light is the sec-ond most important part of the EBL in terms of photonproduction, and is tied with the cosmic infrared background(CIB) in terms of final photon density and integrated pho-ton energy (NB: most CIB photon energy originates fromphotons produced in the COB, which are absorbed and rera-diated to form the CIB, hence the COB dominates in termsof photon production). The COB is one of the least well-known components of the EBL spectra, primarily due tothe uncertainties in subtracting the intervening foregrounds.Formally, the COB is defined to be the emission emitted inthe wavelength range 0.1 µ m to 8 µ m (i.e. the UV to mid-infrared, MIR), and arises primarily from the formation andevolution of stellar mass as galaxies form and evolve (Hauser& Dwek 2001). Hence, the energy density of the COB isclosely related to the growth of stellar mass over all time(Bernstein 1999). There is also a non-negligible contribu-tion from Active Galactic Nuclei (AGN) whereby light, dueto the accretion onto a supermassive BH, escapes from thehost galaxy and radiates through the inter-galactic medium(e.g. Matsuoka et al. 2011; Driver et al. 2016a).The COB also encodes information about the entirecosmic star formation history, as well as the origin and fateof the very first stars responsible for reionisation whichconstitutes 0.1 to 1% of the total EBL (Cooray et al.2009; Windhorst et al. 2018). The reasons above make theCOB one of the most highly studied portions of the back-ground light with tremendous potential to provide strongastrophysical constraints on galaxy formation and evolution.At longer wavelengths the cosmic infrared background(CIB) light is defined as emission at wavelengths between8 µ m and 1000 µ m (Puget et al. 1996; Lagache et al. 1999; Hauser & Dwek 2001; Kashlinsky 2006). As men-tioned, a significant amount of the starlight produced atUV and optical wavelengths, is absorbed by dust in thehost galaxies and later re-radiated in MIR to far-infrared(FIR), contributing to the CIB. The CIB spectral energydistribution (SED) has a peak of the emission at 100 µ mequating to thermal dust emission with temperatures inthe range 20-60K. As detected from the Earth, both theCOB and CIB have approximately the same contribution of24 and 26 nWm − sr − , respectively (Driver et al. 2016b).Thus, the combination of the IR and optical backgroundlight equates, energy-wise, to about ∼
5% of the CMB.Unfortunately, despite significant efforts, the currentmeasurements of the COB and CIB exhibit significantinconsistencies most likely associated with the measure-ment methods (see reviews by Hill et al. 2018; Mattila& V¨ais¨anen 2019; Driver 2020). Conventionally, theseefforts are classified as direct (background sky brightness;Matsumoto et al. 2005; Cooray et al. 2012; Matsuoka et al.2011; Zemcov et al. 2013; Zemcov et al. 2017; Kawaraet al. 2017), and indirect (based on counting galaxies;Madau & Pozzetti 2000; Totani et al. 2001; Xu et al. 2005;Grazian et al. 2009; Keenan et al. 2010; Driver et al. 2016b).Direct measurement is the traditional technique appliedto estimate the COB, including both discrete and diffuseemission. This is based on detecting the absolute flux of thenight sky from space platforms and subtracting the con-tributing foregrounds. The dominant foreground componentoutside of the Earth’s atmosphere is the zodiacal dust glow.This foreground is caused by minuscule dust fragmentsreflecting sunlight and forming the so-called zodiacal light(ZL; Leinert et al. 1981; Bernstein et al. 2002; Aharonianet al. 2006). In addition to the ZL, the diffuse sky emissionmay contain residual airglow (unless the observing platformis entirely beyond Earth’s orbit), and diffuse galactic light(DGL). The optimal technique for separating the DGLand ZL+Airglow remains controversial (Kim et al. 2019).The intensity of the DGL is 5 to 10 times larger than theEBL, and the ZL+Airglow varies from 30 to 100 timesthe EBL (Mattila et al. 2012). These foregrounds are alsotime/day variable (i.e. ZL+Airglow), and so they need tobe estimated for each observation. Therefore, searching foran optimised and well-established model to remove theforegrounds efficiently is of high priority. One suggestedtechnique designed to overcome these backgrounds is theDark Cloud method (Mattila et al. 2017a), in which theeffect of the foregrounds on the background light is mea-sured by estimating the differential measurement betweena dark nebula at high galactic latitude and a surroundingclear region. This results in an upper limit to the UV andblue portion of the EBL.The indirect method, alternatively, represents theintegrated galaxy counts from deep observations anddefines the entire photon flux from discrete sources (i.e.galaxies, quasars) only. In this work, we aim to measurethe integrated galaxy light from a compendium of wideand deep galaxy number count data from the optical to thenear-infrared (NIR). By having sufficiently wide and deepphotometric data, the EBL obtained from the galaxy count
MNRAS , 1– ?? (2021) onstraining the CSFH from the improved EBL method should converge to that from direct measurements.The difference in the two methods can, in due course, pro-vide constraints on any diffuse sources (see the discussion inDriver et al. 2016b). Over the past decades, there has beena rapid increase of new data provided by e.g. HST (HSTCOSMOS; Scoville et al. 2007, HST GOODs; Giavaliscoet al. 2004, HST CANDLES; Grogin et al. 2011),
Spitzer (Sanders et al. 2007, Frayer et al. 2009),
IRIS (De Pontieuet al. 2014),
Gaia (Gaia Collaboration et al. 2016),
WISE (Wright et al. 2010), and
Herschel (Pilbratt et al. 2010).Along with the unprecedented number of new data fromspace-based facilities, there has been a major improvementin deep and wide surveys with ground-based telescopes. Inparticular:
ESO-VLT (Vernet et al. 2009),
Pan STARRS (Kaiser et al. 2002),
Subaru (Taniguchi et al. 2007), and
CFHT (Capak et al. 2007). Nevertheless, the discrepancybetween direct and indirect methods is far larger thanexpected, and the subject of much debate (Mattila &V¨ais¨anen 2019). At present the discrepancy in the CIBestimates are ∼ γ -rays of distant blazars(e.g. MAGIC Collaboration et al. 2008; VERITAS Collab-oration et al. 2011; Dom´ınguez et al. 2011; Biteau 2013;HESS Collaboration et al. 2018; Fermi-LAT Collaborationet al. 2018). This technique relies on the interaction of high-energy TeV photons from blazars with the COB through γ - γ absorption (Desai et al. 2017). Accordingly, this TeV tomicron photon coupling allows the use of γ -rays to probethe optical portion of the EBL over the redshift rangecovered by the blazar distribution. The advantage of thisapproach, in comparison with direct estimation, is that it isunaffected by the foregrounds (ZL and Galactic emission),although arguably susceptible to intergalactic magneticfields. At the present time, the integrated galaxy light (IGL)measurements (EBL estimates in the absence of any diffuselight) agree with the very high-energy (VHE), Dark Cloud,and are mostly consistent with the recent direct estimatesfrom deep space probes ( New Horizon and
Pioneer , whichobserve from beyond the ZL and airglow). This consensusmotivates us to improve our IGL analysis further by lookingto reduce the inherent uncertainty from our previous IGLmeasurements of ∼
20% as reported by Driver et al. (2016b).In Driver et al. (2016b), published data from a variety ofsources was assembled spanning from far-UV (FUV) to FIRwavelengths, i.e., covering both the COB and CIB. Thesedata were used to provide compilations of galaxy number-counts from bright to faint magnitudes and integrated toobtain FUV to FIR IGL estimates. The error analysis ofthat work indicated an uncertainty in the measurements ofaround 20% stemming from several sources. In this paper,we look to make two improvements. The first is to replacethe bright and intermediate data with a more homogeneousanalysis where particular care has been taken to overcomefragmentation and false source detection. The second is animproved error-analysis, in particular an improvement in thesample variance handling and alignment of photometry ontoa common filter set. Through these improvements we look to improve theoptical/NIR IGL uncertainty from 20% to below 10%,as well as identify the key areas for future improvementto eventually bring the uncertainties down to a few percent. This paper is structured as follows: In Section 2, wesummarise the newly processed data from various wide anddeep surveys followed by the colour transformations to acommon filter system. Section 3 shows the revised galaxynumber-counts and the fitting process. In Section 4, wepresent our analysis for estimating the integrated galacticlight of the optical EBL with an improved focus on erroranalysis. Comparison to VHE analysis and implications formissing diffuse light are discussed in Section 5. Finally, weconstrain the CSFH using the EBL estimates described inSection 6.Throughout this paper, all magnitudes are in theAB system and we assume a cosmology with H =70 kms − Mpc − , Ω Λ = 0 .
7, and Ω M = 0 . In Driver et al. (2016b) data for the optical/NIR IGL es-timates were taken primarily from the GAMA, COSMOS,and HST surveys (HST/ACS/WFC3 + LBT/LBC + ESOVLT Hawk-I). Recently improvements have been made inthe analysis of both the GAMA and COSMOS photometry,using the new
ProFound source finding code (Robothamet al. 2018; Bellstedt et al. 2020a; Davies et al., in prep).Here we exclusively focus on three contributing datasets:GAMA, DEVILS-COSMOS, and HST, as briefly describedbelow.
Galaxy And Mass Assembly (GAMA; Driver et al. 2011;Liske et al. 2015) is a spectroscopic (98% completeness)multi-wavelength campaign across ∼
230 deg in the North-ern and Southern Galactic caps, centered on R.A.=2h, 9h,12h, 14.5h, and 23h, with the AAOmega fiber spectrographon the 3.9 meter Anglo Australian Telescope (AAT) (Hop-kins et al. 2013). GAMA includes a panchromatic surveyin 20 filters extending from the UV to the FIR (Driveret al. 2016a), composed of observations from various wideand deep ground-based (VST KiDS; VISTA VIKING) andspace-based ( GALEX ; WISE ; Herschel ) facilities as sum-marised in Driver et al. (2016a). GAMA’s main science goalis to quantify the evolution of mass, energy, and structureover the history of the cosmos (Driver et al. 2018; Moffettet al. 2018; Driver & GAMA Team 2016). GAMA targetsmore than 200,000 galaxies to a depth of r < < ∼ . Previously, GAMA optical pho-tometry as used in Driver et al. (2016b), was based on SEx-tractor analysis of the SDSS data (see Driver et al. 2016b),whereas the reanalysis is based on
ProFound (Robotham
MNRAS , 1– ????
MNRAS , 1– ???? (2021) S. Koushan et. al. et al. 2018). Improvements include the use of dilated seg-ments, watershed deblending, higher-resolution data (KiDS– SDSS), improved star-galaxy separation, extensive visualinspection, and manual fixing of fragmented galaxies (seeBellstedt et al. 2020a for full details). These enhancementsimprove the integrity of the counts at both very bright andintermediate magnitudes, as well as extending the data tofainter limits than was possible with the SDSS data, and inparticular to the depths that dominate the IGL measure-ment.
The Deep Extragalactic VIsible Legacy Survey (DEVILS;Davies et al. 2018), is a spectroscopic survey on the Anglo-Australian Telescope (AAT), sampling more than ∼ Y < . of the sky aiming for more than 95% completeness inthree deep fields (COSMOS; Scoville et al. 2007, ECDFS;Virani et al. 2006, XMM-LSS; Pierre et al. 2004). Here, weadopted the up-to-date photometric catalogue for the D10(COSMOS) and D02 (CMM-LSS) regions which cover ∼ ( Y JHK s ) of the UltraVISTA and VIDEO datasets(McCracken et al. 2012; Jarvis et al. 2013) centered at RA10h and DEC 2.2 ◦ for D10 and RA 2.2h and DEC -4.5 ◦ for D02 (Davies et al. in prep.). The imaging datasets in ugrizYJHKs filters come from the following sources: the u band originates from the CFHT, the griz bands from the Hy-per Suprime Camera Strategic Program (HSC; Aihara et al.2018), and the YJHKs bands from either UltraVISTA-DR3or VIDEO XMM-LSS. For more detail on the photometriccatalogue, see the DEVILS photometry description paper(Davies et al., in prep). In the same manner as the updatedGAMA photometry, the photometric catalogue has also re-cently been reassembled using the
ProFound source find-ing analysis package from the FUV to FIR (Davies et al, inprep).Improvements over previous photometry include thosehighlighted for GAMA through the use of
ProFound ,bringing consistency of methodology, with comparable levelsof visual checks and manual fixes.
For the deepest data, we adopt the deep number countsfrom Driver et al. (2016b), which were mostly obtainedfrom the HST, Wide Field Camera 3 (WFC3) Early ReleaseScience (ERS) archive (WFC3-ERS; Windhorst et al. 2011)in combination with the HST Ultra-Violet Ultra-Deep fields(UVUDF; Teplitz et al. 2013). The data reduction of thepanchromatic WFC3 ERS filters are discussed in Windhorstet al. (2011). They derived deep field galaxy number countsacross the WFC3 ERS filters (F336W, F435W, F606W,F775W, F850LP, F098M, F125W, and F160W) from uband to Y band. The panchromatic photometric countsfor the same filters from NUV to NIR in the Hubble UltraDeep Field (HUDF) are described in Rafelski et al. (2015).We note that other HST fields exist, however, wehave chosen to exclusively use the ERS and UV-UDFdata as they have independent detections in each band, thereby circumventing the colour bias that is introducedthrough forced aperture photometry. We note that a fullreanalysis of the entire ACS WFC3 archive is currentlyunderway, as part of the SkySURF HST Legacy Program(PI: Windhorst), and the inclusion of additional HST fieldsin IGL EBL analysis will be the subject of a future paper.In the u band we retain the LBT data from Grazianet al. (2009) and in the K s band we also retain the deepdata from Fontana et al. (2014) using ESO VLT Hawk-Iwith an exposure time of 13 hours, and reaching a depth of27.3 AB magnitude (1 σ limit).Table 1 summarises the wavelength coverage and depthof data used in this work. The datasets described above come from a variety of facil-ities with aligned but slightly differing filters and facilitythroughput transmissions. To derive galaxy number countsacross consistent bands, we need to implement colour trans-formations to a common filter set. Here we elect to use theVST ugri and VISTA
ZY JHK s filters as our EBL refer-ence filters. The transformations we need include the DEV-ILS optical data (CFHT u ∗ and HSC griz filters), all HSTdata (F336W, F435W, F606W, F775W, F850LP, F098W,F105W, F125W and F160W), the LBT data ( U
360 filter),and the ESO Hawk-I ( K s filter). For the DEVILS data, wehave full access to the flux measurements for each individualgalaxy through our segment-matched ProFound photom-etry. Hence, we can convert robustly from native to KiDSfilters, by including a colour term to fold in the second-order dependence on the spectral shape. For the HST+datasets we only have access to the number-count data, which have mostly been determined independently in eachband , and hence can only determine a mean first order fil-ter offset given some reasonable assumption as to the likelymean spectral shape and redshift distribution of the underly-ing galaxy population (see also the discussion in Windhorstet al. 2011).To determine our filter transformations, we use the
ProSpect package (Robotham et al. 2020) to generate sev-eral thousand galaxy spectra with a uniform redshift dis-tribution from z=0 to z=2, and using the
ProSpect de-fault metallicity, dust properties, a Chabrier IMF (Chabrier2003), and a randomly selected star formation history(SFH). We then determine the fluxes for these spectra asobserved through our reference filters and through the filtersfor each of the datasets we wish to transform. We fit for thedifference between the measured fluxes in the original andtarget filters, with a linear dependence on the colour termfor DEVILS, and without a colour term for the HST+ data.As part of the fitting process, we determine an estimateof the intrinsic uncertainty in the transformation process,which will of course be higher for the HST+ datasets (be-cause we lack the colour terms). We note that in the futurewe are looking to reprocess the entire HST archive as partof the
SkySurf
HST Archival program (PI: R. Windhorst),which will allow us, in due course, to factor in the properHST colour terms. Table 2 shows the resulting transforma-tions for the filters used in this work which we now apply.
MNRAS , 1– ?? (2021) onstraining the CSFH from the improved EBL Table 1.
A summary of multi-wavelength data used in this analysis. The reference surveys for the analysis are as follows: GAMA(Bellstedt et al. 2020a); DEVILS (Davies et al, 2020 - in prep); HST ERS Data 3 and UVUDF (Windhorst et al. 2011, Rafelski et al.2015, Grazian et al. 2009, and Fontana et al. 2014). The pivot wavelength refers to the GAMA data.EBL Survey EBL Pivot Depth (5 σ AB)Filter WavelengthGAMA DEVILS HST+ ( µ m) GAMA DEVILS u VLT KiDS CFHT(u ∗ ) F336W, LBT a U360 b g VLT KiDS HSC c g F435W 0.4744 25.1 27.3 r VLT KiDS HSC r F606W 0.6312 24.9 26.9 i VLT KiDS HSC i F775W 0.7584 23.8 26.7 Z VISTA VIKING HSC z F850LP 0.8833 23.1 26.3 Y VISTA VIKING UltraVISTA F098W, F105W 1.0224 22.3 24.7 J VISTA VIKING UltraVISTA F125W 1.2546 22.1 24.5 H VISTA VIKING UltraVISTA F160W 1.6477 21.5 24.1 K s VISTA VIKING UltraVISTA ESO K d Table 2.
Filter transformations used to place all data onto the
VST/VISTA filter systems (i.e., ugriZYJHKs) using our onlinetool at: https://transformcalc.icrar.org .Filter transformation IntrinsicUncertainty u VST = u ∗ CFHT + 0 . u ∗ CFHT − g HSC ) − . g VST = g HSC + 0 . g HSC − r HSC ) + 0 . r VST = r HSC − . g HSC − r HSC ) + 0 . i VST = i HSC + 0 . r HSC − i HSC ) + 0 . Z VISTA = z HSC + 0 . i HSC − z HSC ) + 0 . u VST = U LBT + 0 . u VST = F W HST − .
113 0.089 g VST = F W HST − .
159 0.112 r VST = F W HST − .
114 0.070 i VST = F W HST + 0 .
031 0.019 Z VISTA = F LP HST + 0 .
057 0.049 Y VISTA = F W HST − .
076 0.053 Y VISTA = F W HST + 0 .
066 0.078 J VISTA = F W HST − .
015 0.022 H VISTA = F W HST − .
066 0.024 Ks VISTA = Ks Hawk − I + 0 .
000 0.002
Figure 1 shows our revised galaxy number count compila-tions in ugriZY JHK s bands as derived from the deep andwide data (GAMA, DEVILS, and HST+). Figure 2 showsthe same data, but with a linear term of d ( log ( N )) /dm =0 . ∝ (cid:90) N m × − . · m dm (1)Where N is the number of galaxies per magnitude bin be-tween m and m + δ m .Each dataset is truncated at the faint end at thepoint at which it becomes inconsistent with deeper data(see Fig. 1). As can be seen, the overlap between thedatasets is smooth, with no obvious discontinuities. As aresult, the galaxy number counts indicate good agreementin all bands with the largest discrepancies seen in the u and g bands where more robust Galactic extinctioncorrections may become important (in terms of bothspatial distribution and reddening). Driver et al. (2016b)discussed possible sources of both systematic and statisti-cal error. These will be discussed in more detail in Section 4.Figure 2 shows the same plot but now in terms of thedifferential luminosity density contribution from each mag-nitude interval. Hence it is the integral under this curve thatprovides the final EBL measurement in that band. Figure 3shows the cumulative version of the same plot. In order to estimate the total EBL from discrete sources ineach band, we fit a smooth spline with 10 degrees of free-dom (df=10) to the combined luminosity density data, in-verse weighted by the variance. We then integrate the splinefit from -100 to +100 AB magnitude. Assuming no unex-pected behaviour of the spline outside the data range, Fig.2, shows the contribution of each magnitude bin to the EBLin all filters (as indicated) in the unit of WHz − m − per 0.5magnitude bin, and normalised by effective area. As oughtto be expected, the contribution to the IGL is limited atboth very faint and very bright magnitudes with most ofthe EBL contribution coming from intermediate magnitudes MNRAS , 1– ????
Figure 1 shows our revised galaxy number count compila-tions in ugriZY JHK s bands as derived from the deep andwide data (GAMA, DEVILS, and HST+). Figure 2 showsthe same data, but with a linear term of d ( log ( N )) /dm =0 . ∝ (cid:90) N m × − . · m dm (1)Where N is the number of galaxies per magnitude bin be-tween m and m + δ m .Each dataset is truncated at the faint end at thepoint at which it becomes inconsistent with deeper data(see Fig. 1). As can be seen, the overlap between thedatasets is smooth, with no obvious discontinuities. As aresult, the galaxy number counts indicate good agreementin all bands with the largest discrepancies seen in the u and g bands where more robust Galactic extinctioncorrections may become important (in terms of bothspatial distribution and reddening). Driver et al. (2016b)discussed possible sources of both systematic and statisti-cal error. These will be discussed in more detail in Section 4.Figure 2 shows the same plot but now in terms of thedifferential luminosity density contribution from each mag-nitude interval. Hence it is the integral under this curve thatprovides the final EBL measurement in that band. Figure 3shows the cumulative version of the same plot. In order to estimate the total EBL from discrete sources ineach band, we fit a smooth spline with 10 degrees of free-dom (df=10) to the combined luminosity density data, in-verse weighted by the variance. We then integrate the splinefit from -100 to +100 AB magnitude. Assuming no unex-pected behaviour of the spline outside the data range, Fig.2, shows the contribution of each magnitude bin to the EBLin all filters (as indicated) in the unit of WHz − m − per 0.5magnitude bin, and normalised by effective area. As oughtto be expected, the contribution to the IGL is limited atboth very faint and very bright magnitudes with most ofthe EBL contribution coming from intermediate magnitudes MNRAS , 1– ???? (2021) S. Koushan et. al.
Figure 1.
The galaxy number count in each filter for a combination of datasets derived from the GAMA, DEVILS, and the HST surveysin blue, green and orange respectively. We fitted a smooth spline in a red dashed line through all datasets. For a sanity check, we compareto published data (Driver et al. 2016b). The black solid line shows a fit to their data. The error bars indicate the Poisson error. (and redshift). Following integration we convert to the stan-dard EBL units using Eq. 2.
EBL = u × . × × ( c/λ ) (2)where u is the energy density, c is the speed of light in m/s , and λ is the wavelength in meters, and the constantconverts the measurements into units of nWm − sr − . Notethat the fraction of the IGL in the extrapolated of the spline fit is typically less than a few per cent, see Fig. 3.Figure 2 shows the luminosity density with the splinefits overlaid (red dashed line) for each band. We compare ourbest-fitted splines in red dashed lines to the compendium ofDriver et al. (2016b) shown as solid black lines. Our data areclearly bounded at both ends in all filters and show a welldefined peak at intermediate magnitude levels in all bands. MNRAS , 1– ?? (2021) onstraining the CSFH from the improved EBL Figure 2.
The contribution of each magnitude bin to the energy density. The count has been normalised and scaled for each survey indepth and area.
In this section, we explore error estimation taking into con-sideration the systematic and random errors for each band.Our systematic errors break down into cosmic variance (CV)and zero-point (ZP) uncertainty, while random errors aredue to Poisson error, as well as the uncertainty inherent inthe fitting process. The final uncertainty is then the combi-nation of these four separate errors, which we explore in afull Monte Carlo (MC) analysis. Below we discuss each of these individually to address their impact on our final EBLCOB measurements.
To estimate the impact of Poisson uncertainty alone, weperformed a Monte Carlo analysis with 1001 iterations.Within each iteration, we perturb the galaxy number-countsby their associated density uncertainty error, and repeatour full analysis. For the perturbation of each data point
MNRAS , 1– ????
MNRAS , 1– ???? (2021) S. Koushan et. al.
Figure 3.
The cumulative contribution function ( ρ L ) to each band. we sample a Gaussian distribution with mean zero and astandard deviation equivalent to the quoted number-counterror for that magnitude bin. The top row of Fig. 4 showsa range of the MC fits (left panel), and the recovered IGLmeasurement (right) for 1001 realisations in the Z-band.From this we determine a mean value of 10.658 and astandard deviation of 0.035 i.e., <
1% error. Table 3 showsthis estimate due to Poisson error from u to K s and ingeneral this error is fairly negligible throughout, as expectedgiven the large sample size. Errors in the absolute zero-point calibration will also af-fect the number-counts in a systematic manner for a givendataset, but randomly between the datasets, and randomlybetween filters. Here we estimate the zero-point errors inthe GAMA/KiDS data by comparing the KiDS/VIKINGphotometry to the SDSS/2MASS photometry, incorporatingappropriate filter conversions as provided by Kuijken et al.(2019) for KiDS (only relevant for u and r band, see theirfigure 5 and eq 5), and Gonz´alez-Fern´andez et al. (2018) MNRAS , 1– ?? (2021) onstraining the CSFH from the improved EBL Table 3.
Integrated EBL results considering individual errors. Column 2 shows our IGL measurement in each band. Columns 3 to 6show the mean offset in % for the contribution of each uncertainty. Finally, column 7 presents the total error (added in quadrature) asa percentage of the IGL.Filter EBL Poisson Error Cosmic Variance Zero-Point Fitting Error Total( nW/m sr − ) (%) Error (%) Error (%) (%) (% of IGL)u 4.13 0.38 2.54 6.36 0.35 6.87g 5.76 0.49 3.29 2.07 0.92 4.02r 8.11 0.37 3.55 1.44 1.34 4.08i 9.94 0.35 3.98 1.67 0.99 4.44Z 10.71 0.30 3.98 3.16 1.06 5.20Y 11.58 0.32 4.11 1.62 0.88 4.52J 11.22 0.29 4.25 2.35 0.75 4.92H 11.17 0.26 4.38 1.63 0.65 4.73Ks 9.42 0.19 4.12 3.09 0.74 5.21 for VISTA (in ZY JHK s bands, see their eqs 5—9 and alsoAppendices A and C). In effect, we are attempting to pro-vide an independent albeit basic calibration check of theKiDS and VIKING data, compared to the more sophisti-cated calibration conducted by the KiDS and VISTA CASU(Cambridge Astronomy Survey Unit) teams. The zero-pointverification analysis is done on a square-degree by squaredegree (i.e., TILE ) basis, and we exclusively use stars in arestricted magnitude range, where bright stars are not satu-rated in KiDS/VIKING, and faint stars are not dominatedby sky noise for SDSS/2MASS. Left panel in Fig. 5 showsan example comparison for the ugriZY JHK s filters for asingle TILE centred at RA=8.5h and Dec= − . o , the yel-low data points highlight the magnitude range used, whilethe grey data points showing the full range of the compar-ison. From the yellow data points, we determine a simpleoffset and a linear fit for each filter and for each TILE , asindicated by the cyan and mauve lines respectively. In gen-eral the fits are good, showing relatively little magnitudedependence, and hence we use the simple offsets as our indi-cator of the zero-point offset for each
TILE . Figure 5 (right)shows the histogram of these zero-point offsets for the 220
TILES that make up the GAMA KiDS/VIKING datasetreflecting a distribution of zero-point uncertainty. For eachfilter, we determine the median absolute zero-point offset,and the 1- σ zero-point uncertainties from a simple Gaussianfit, and the empirical standard deviations (as indicated inthe Figure panels). We also note that the KiDS team iden-tified a Galactic Latitude dependency of –0.02 to 0.1 mag inthe u band, potentially explaining the much broader spreadthat we see here. We also note that the VISTA zero-pointsevolved from VISTA reduction 1.3 to 1.5 with implied zero-point changes of –0.03, 0.018, –0.0200, 0.0067 and 0.0106,for the five VISTA bandpasses ( ZY JHK s ) respectively (seeGonz´alez-Fern´andez et al. 2018, Appendix D). We thereforeelect to adopt worst case scenario for the zero-point errorsby combing all of the uncertainties in quadrature (i.e., lit-erature offset + median offset + (1- σ ) spread + intrinsicuncertainty from Table 2) to obtain ultra-conservative zero-point uncertainties of: 0.11, 0.022, 0.011, 0.018, 0.041, 0.021,0.032, 0.022 and 0.037 for the GAMA ugriZY JHK s respec-tively.For the DEVILS and HST+ datasets we adopt the zero-point uncertainties as the maximum of either those foundfor GAMA, or an error floor of 0.03 mag for DEVILS (i.e.,that adopted in Driver et al. 2018) or HST+. In addition we also fold into the DEVILS and HST+ zero-point errors theintrinsic uncertainty from the filter transformations shownin Table 2. Technically, this overestimates the likely colourtransformation error (as it will be random not systematic),but this is a convenient place to incorporate it that errs onthe side of caution. Table 4 shows the final adopted zero-point errors for each dataset and each band.With conservative zero-point errors in hand, we can nowMonte Carlo simulate the impact of the zero-point error byperturbing the magnitudes for each of our three datasets(GAMA, DEVILS, and HST) in a systematic manner, butindependently for each dataset and for each filter. To dothis, we draw a random value for each dataset and filterfrom a Gaussian distribution with a mean of zero and astandard-deviation defined by the relevant zero-point un-certainty (Tab 4), and systematically adjust the magnitudebins of each survey by this amount for that filter. Finally, werefit our splines and re-derive our optical/NIR IGL measure-ments 1001 times to build up a distribution of the impact ofthe zero-point errors. The results are shown in Table 3 andFig. 4. Pragmatically, cosmic variance can be considered as a sys-tematic uncertainty inherent in observational estimates ofthe volume of the extragalactic sky due to large-scale densityfluctuations. This error is quite significant in most galaxysurveys, and is often the dominant error in measuring theIGL (especially in deep galaxy surveys).For the GAMA survey, we obtained a mock catalogueusing the
SHARK semi-analytic model of galaxy formationand evolution (Lagos et al. 2019; Bravo et al. 2020) used toconstruct light-cones for the GAMA regions. To estimate ourCV error, we can determine the effective volume from the16% and 84% quantiles of the predicted redshift distribu-tion, n ( z ), for the magnitude slice in question, and evaluateequation 3 from Driver & Robotham (2010). We use an ef-fective sky area of 60 deg , redshift ranges of the comovingcone, and set the number of separated regions to 4 for eachsurvey to obtain the CV error. Top panel in Fig. 6 showsan example redshift distribution for the magnitude range of20.5 ± r band. Fig. 6 (bottom panel) also high-lights the variation of the cosmic variance as a percentage(solid red line) compared to the zero-point (solid green line)and Poisson (solid blue line) uncertainties. MNRAS , 1– ????
SHARK semi-analytic model of galaxy formationand evolution (Lagos et al. 2019; Bravo et al. 2020) used toconstruct light-cones for the GAMA regions. To estimate ourCV error, we can determine the effective volume from the16% and 84% quantiles of the predicted redshift distribu-tion, n ( z ), for the magnitude slice in question, and evaluateequation 3 from Driver & Robotham (2010). We use an ef-fective sky area of 60 deg , redshift ranges of the comovingcone, and set the number of separated regions to 4 for eachsurvey to obtain the CV error. Top panel in Fig. 6 showsan example redshift distribution for the magnitude range of20.5 ± r band. Fig. 6 (bottom panel) also high-lights the variation of the cosmic variance as a percentage(solid red line) compared to the zero-point (solid green line)and Poisson (solid blue line) uncertainties. MNRAS , 1– ???? (2021) S. Koushan et. al.
Table 4.
Adopted zero-point uncertainties for each dataset and each filter.Dataset u g r i Z Y J H K s GAMA 0.11 0.02 0.01 0.02 0.04 0.02 0.03 0.02 0.04DEVILS 0.11 0.03 0.03 0.03 0.04 0.03 0.03 0.03 0.04HST+ 0.11 0.08 0.13 0.06 0.04 0.04 0.07 0.03 0.04
Figure 4.
An example of a luminosity density plot including allpossible sources of error including the Poisson, CV, ZP, and fittingin blue, orange, green, and red for the Z band, respectively. Leftcolumn shows the contribution of each uncertainty to the EBL.Right column shows the distribution of our IGL for each source oferror. The mean and standard deviation (1 σ ) are shown in blackand gray solid lines respectively. For the DEVILS/D10+D02 and HST surveys, we alsoadopted the cosmic variance uncertainty following eq 3 ofDriver & Robotham (2010) for the appropriate sky area forthe particular survey, and adopting an approximate redshiftrange. To determine the approximate redshift ranges to usefor the DEVILS/D10 + D02 and HST surveys, we sliced amock galaxy catalogue generated by the SHARK team intothe 16% and 84% quantile ranges independently for eachmagnitude bin and each bandpass. For our three datasets, wethen calculated CV fluctuations for each magnitude bin, andperturb the galaxy counts in the magnitude direction usingthis value for each filter and independently for each dataset.We perform a Monte Carlo simulation with 1001 iterationsand, consequently, rederived the IGL EBL. Throughout theCV analysis, the procedure was forced to start from the samerandom seed to ensure that the CV error is correlated acrossthe filters. Table 3 shows the error from the CV.
To calculate the IGL, we performed a spline fit using anarbitrarily selected 10 degrees of freedom. Figure 4 (lower)shows the shift in the Z band IGL value if we use a spline oforder 8, 9, 10, 11 or 12. On the whole, the shift is marginal,as compared to that from the Poisson error (see top panel inFig. 4), and far less than the effect of either cosmic varianceor the zero-point errors. Table 3 summarizes the integrated EBL results showing theindividual errors. Figure 7 shows the contribution of eacherror as a percentage of the optical/NIR IGL in a singlegraph. We show the photometric, spline fitting, Poissonand cosmic variance error in green, red, blue and orangelines, respectively. As can be seen, the systematic cosmicvariance and zero-point errors dominates for all filters, whilethe contribution of random errors (fitting and Poisson)is negligible. The solid line also highlights the total errorbudget for all optical and near-IR filters ( ugriZY JHK s )comparing to the previous result from Driver et al. (2016b);shown as a black line.We finally combine the errors in quadrature, to produceour final IGL error estimates. However we note, that theerrors are a mixture of both random and systematic errors.Table 3 (last column) shows our final error estimate asa percentage of the IGL value. As a consequence of therevised data, colour transformations and revised erroranalysis, the IGL error, compared to the earlier results ofDriver et al. (2016b), is reduced from ∼
20% to the meanerror of ∼ MNRAS , 1– ?? (2021) onstraining the CSFH from the improved EBL Figure 5.
Left panel: A direct comparison of the GAMA
ProFound photometry to either SDSS ( ugriZ ) or 2MASS (
Y JHK s ) data,using objects classified as stars. The region we use for determining offsets for each TILE and for each filter are shown in yellow, with theoffsets and line fits for each sample shown in blue or mauve respectively. Right panel: the compendium of zero-point offsets for all 220
TILES with the σ of the Gaussian fit and the standard deviation shown on each panel. The top panel in Fig. 8 shows the comparison between ournew IGL measurements and the results from Driver et al.(2016b), who used essentially the same methodology, butwhich we have improved by the replacement of SDSS datawith VST KiDS data, a full reanalysis of the core GAMAand COSMOS/DEVILS number-counts data, as well as animproved error analysis. We also include the very recentdata from the New Horizons LORRI (Lauer et al. 2020)using the Zemcov et al. (2017) (green triangle) estimateswhich were obtained at an even greater distance of > µ m at high galactic latitude after sub-tracting estimates of the outer Kuiper belt light and DGL.The Lauer et al. (2020) IGL estimate is measured as 7.37 ± ± − sr − , which looks consistent with an interpolationof our IGL points at 0.6 microns.The ratio of our revised measurement of the IGL for each fil-ter, to that from Driver et al. (2016b) is shown in the lowerpanel of Fig. 8, and highlights the changes resulting fromthe current analysis. As the areas covered by the contribut-ing surveys are essentially identical the reduction in error comes mainly from the improvement in our modeling of thecosmic variance, which now relies on the SHARK galaxyformation simulations (see Lagos et al. 2019), and improvedcolour corrections (see Robotham et al. 2020). In all wave-bands, we see modest increase in the recovered IGL valueof approximately 5-15 per cent (see lower panel of Fig. 8).While in all bands the change is within the errors, it is sys-tematic across the filters. We believe this is due to a combi-nation of ProFound recovering more accurate total fluxesthan the previous source detection algorithms (see Fig. 14in Robotham et al. 2018), the improved depth of KiDS overSDSS, improved star-galaxy separation (see Bellstedt et al.2020a), and to a lesser extent more accurate accounting ofthe masked area, i.e., that lost near bright foreground stars.We note that this increase in the IGL measurements bringsinto line the IGL and the VHE data, available at that time,in which Driver et al. (2016b) noted a 10-20% discrepancy.At the time this was tentatively attributed to the possibil-ity of stripped mass (e.g., ICL and Intra-Halo Light) and/ormissing populations. See also all discussion in Ashcraft et al.(2018). They found ≤ ∼
32 mag.
MNRAS , 1– ????
MNRAS , 1– ???? (2021) S. Koushan et. al.
Figure 6.
Top panel: Cosmic variance density distribution for asample of a magnitude bin in the r band. The dashed black linedenotes the mean redshift. Upper/lower quantiles are displayedin dashed red lines. Bottom panel: The cosmic variance comparedto the zero-point and Poisson uncertainties. As indicated in the introduction, VHE studies estimated theEBL intensity through the production of electron-positronpairs as a by-product of the interaction of high energyphotons, emanating from the distant blazars, with micronwavelength photons in the EBL. This interaction resultsin the decrement of the expected power-law distribution ofhigh energy photons with 100GeV blazar photons preferen-tially interacting with micron wavelength EBL photons. Atpresent, most VHE applications require an input EBL modelcurve and constrain the normalisation of the curve. Figure9 presents the most recent results from Fermi-LAT Collab-oration et al. (2018), who studied 759 active galaxies as ob-served via the Fermi Large Area Telescope (LAT). In Fig. 9we compare our new IGL measurements to published VHECOB measurements (Biteau & Williams 2015; Ahnen et al.2016; HESS Collaboration et al. 2018; Fermi-LAT Collabo-ration et al. 2018). This latest result shows full consistencywith our optical/NIR IGL data, however, we note that animproved comparison could be made by a re-analysis of theVHE data which adopts an EBL model fitted to our IGLdata (see Sec. 6).
Figure 7.
A contribution of each systematic and random uncer-tainties to the IGL for each filter. The contribution of the cosmicvariance to the IGL is dominant across the optical and NIR fil-ters. Random errors have an insignificant contribution to the IGLestimation. Our total error (in gray line) has been compared toDriver et al. (2016b) marked with black line.
Current attempts, within the VHE community, arenow ongoing to not only constrain the EBL normalisation,but also the shape of the EBL COB SED (see for exampleAbeysekara et al. 2019). Recently, the VERITAS team(Biteau & Williams 2015) demonstrated this by detectingand measuring the spectra of 14 blazars at very high energy( >
100 Gev) in order to enable a full reconstruction of theSED of the EBL. However, at this stage the analysis errorsare too large to be useful.The shaded area obtained by the VHE COB measure-ment also depicts the range of the allowed EBL intensityin the optical/near-IR filters (from 0.24 µ m to 4.25 µ m forMAGIC and extending to 10.4 µ m for H.E.S.S). At this stageit can only really be stated that the VHE data are fullyconsistent with our IGL measurements in the u − K s rangewithout the need to include any significant additional sourceof diffuse light.As the errors in both the optical/NIR IGL and VHECOB data improve to the few percent level, the comparisonbetween these complementary methods should provide in-teresting constraints on the presence, or lack of non-boundoptical/NIR emission. Figure 9 shows a comprehensive compilation of nearby Solarsystem measurements of the COB. These mainly stem fromdirect background measurements of HST data, and are of-ten shown as either upper limits (Bernstein 2007) or some-times shown as tentative measurements (Hill et al. 2018).The open orange diamonds and gray square show direct es-timates from the Pioneer 10/11 (Matsuoka et al. 2011) andNew Horizons (Zemcov et al. 2017) spacecraft at distances > MNRAS , 1– ?? (2021) onstraining the CSFH from the improved EBL Figure 8.
Top panel: The IGL as a function of wavelength esti-mated in Driver et al. (2016b) and this work marked with cyandiamonds, and red circles, respectively. Shaded areas correspondto the total error. For comparison, we show the IGL estimate fromLauer et al. (2020) as a lower limit. As can be seen, their estimateis consistent with our measurement. Bottom panel: The ratio ofIGL measured in this work to the IGL measured by Driver et al.(2016b). The red dot-dashed line corresponds to the equality linebetween the two measurements. the reanalysed EBL beyond 4.5 a.u. as 7 . − . . nWm − sr − and 7 . − . . nWm − sr − at 0.44 and 0.64 µ m respectively.New Horizon also reported the EBL as 4 . − . . nWm − sr − at 0.66 µ m. Both of these measurements are sufficiently farfrom the Sun to be considered near free of Zodiacal lightcontamination.As clearly shown in Fig. 9, our optical/NIR IGLmeasurements are inconsistent with the interpretation ofthe HST -based direct estimates as tentative measurements,but consistent as upper limits. As outlined in Bernstein(2007), the key issue is robust subtraction of bright fore-grounds such as the diffuse Galactic light and the Zodiacalemission. The discrepancy is less prominent between thedirect measurements from deep-space probes such as thePioneer 10/11 or the New Horizon mission but to a lesserextent. It would therefore seem that the likely issue, at leastin comparison to the HST data, is with the subtraction ofthe Zodiacal light, which drops substantially as one movesoutward in the Solar system, or potentially due to addi- tional airglow contamination from the upper atmosphere.For completeness, we also include recent result from thedark cloud method (purple diamond; Mattila et al. 2017b),which excluded bright galaxies < ∼
22 mag and resulted in11 . − . . nWm − sr − at 400 nm as an upper limit.Hence some small tension does still remain between thedeep space probe data, the dark cloud data, and the IGLdata which as argued by Lauer et al. (2020) may leave someroom for a diffuse COB component of unknown origin. Oneargument against this, is the strong agreement between theIGL data which is sensitive to the discrete COB, and theVHE data (discussed earlier) which is sensitive to both thediscrete and diffuse COB. Hence we are in the intriguingsituation where the IGL and VHE are fully consistent, butat low-level tension with the deep space platform and darkcloud data. An obvious solution is that the additional ra-diation being detected by the space-platform data may beassociated with either the outer Solar System, additionalDiffuse Galactic Light, or the outer Milky-Way dark-matterHalo, i.e., local. Resolving this discrepancy will likely requireimprovements in all three methods.Finally, in Fig. 9 we also show three recent modelsalong with our data points. SHARK (Lagos et al. 2019;dashed green line) is a semi-analytical model of galaxyformation and evolution (SAM) extending from FUV toFIR, the model of (Khaire & Srianand 2019, solid orangeline) is a synthesis model of the EBL from FIR to highenergy γ -rays, and Andrews et al. (2018) (solid purple line)presents a phenomenological model of the cosmic spectralenergy distribution (CSED) based on simple forwardmodelling from an input cosmic star formation history(CSFH). While all the models show good consistency theymostly lie outside of the error ranges of our data heraldingthe prospect of using our data to now test the models morestringently as our errors continue to improve.In this work, we have shown that the optical/NIR IGLuncertainty has been improved (see Fig.7) in the opticalregime from 20% (Driver et al. 2016b) to below 10%. Figure6 highlights the dominant error across the magnitude rangein the r band. This shows that in order to further improveIGL measurements we need wider area data at brightermagnitudes as well as improved calibration, including filterconversions, at fainter magnitudes. Both of these will besatisfied through the WAVES and HST SkySurf projects.The former will reprocess the entire VISTA VIKING/VSTKiDS region increasing the wide area coverage from 230deg to 1300 deg , while the SkySurf project will includea complete reanalysis of HST data using
ProFound withthe opportunity to perform more robust colour correctionson a galaxy-by-galaxy basis. Both of these improvementshave the potential to reduce the IGL uncertainty to < < MNRAS , 1– ????
ProFound withthe opportunity to perform more robust colour correctionson a galaxy-by-galaxy basis. Both of these improvementshave the potential to reduce the IGL uncertainty to < < MNRAS , 1– ???? (2021) S. Koushan et. al.
Figure 9.
Our measurement of the IGL based on galaxy number counts from wide and deep surveys (GAMA+DEVILS+HST) in ninemulti-wavelengths from u to K s band (red circles). We compared our IGL to literature measurements taken from various techniques. Directestimates: (Bernstein 2007) in blue triangles; New Horizons (Zemcov et al. 2017) as an upper limit in gray; Pioneer 10 /
11 (Matsuoka et al.2011) in orange diamonds; recent estimates from the New Horizons LORRI (Lauer et al. 2020) in green triangle. VHE estimates: (MAGICCollaboration et al. 2008) in yellow shades; (HESS Collaboration et al. 2018) in pink shades; (Fermi-LAT Collaboration et al. 2018) ingray shades; γ -ray observation from (Biteau & Williams 2015) in violet circles. Dark Cloud: (Mattila et al. 2017a) and upper limits inpurple. We also compared the EBL estimation obtained from the models: SHARK (Lagos et al. 2019) green dashed line; (Khaire &Srianand 2019) orange line; GAMA (Andrews et al. 2018) purple line. Finally, we added a data assembled by Hill et al. (2018) from radioto γ -ray wavelengths as an upper and lower limit following by the red and blue arrows respectively. As outlined in the introduction, the COB (and CIB) isthe radiation by-product of cosmic star-formation over alltime, and as such is directly predictable given a cosmicstar formation history (CSFH), and what we know aboutstellar evolution and dust attenuation. Recently, Driveret al. (2012) and Andrews et al. (2018), built models topredict the EBL (Driver et al. 2016b), and its subdivisionby redshift into CSED slices. These models require as inputa CSFH parametrisation, stellar libraries, e.g.,
Pegase2 orBC03 (Bruzual & Charlot 2003), an initial mass function(IMF), and simple axioms around metallicity evolution(typically linked to star-formation), dust attenuation anddust emission, and, in the Andrews et al. (2018) case,the ability to incorporate obscured and unobscured AGNcomponents. Here, we essentially repeat this process, butwe now adopt the comprehensive package
ProSpect (Robotham et al. 2020). This package has been developedto both model and generate individual galaxy SEDs, withsignificant robustness, flexibility and functionality. It alsoincludes a sub-module to provide an EBL prediction from UV to FIR, given the kind of assumptions described above.
ProSpect therefore provides us with the opportunity touse the COB to constrain model inputs to the EBL code,and in particular the CSFH normalisation (see Robothamet al. 2020; Bellstedt et al. 2020b).The CSFH is a fundamental empirical measurement de-scribing the evolution of the galaxy population from thepast to the present time. Empirically the CSFH has beenmeasured via a number of complementary methods i.e.broadband and spectral line measurements (the compila-tion of measurements by Madau & Dickinson 2014, hereafterMD14); The core sampling technique results from
MAG-PHYS (Driver et al. 2018, D18); VHE constraints (Fermi-LAT Collaboration et al. 2018, FL18); SED forensic recon-struction (Bellstedt et al. 2020b, B20) and IFU forensic re-construction (S´anchez et al. 2019, S19; L´opez Fern´andezet al. 2018, L18). On the whole the distinct methods for re-constructing the CSFH agree on the overall trend and shape,i.e., a rapid rise to cosmic noon at z ≤
2, and thereafter agentle continuous decline. In detail, the range of CSFH mea-
MNRAS , 1– ?? (2021) onstraining the CSFH from the improved EBL surements exhibit an uncertainty of almost a factor of 3 atcosmic noon (see Fig. 10). In this work we attempt to re-duce the uncertainty by including our COB measurementsto the modeling as a new constraint. To predict the COB using the
ProSpect
EBL model,we represent each measured CSFH as a simple analyticfunction for computational purposes. We elect to use theskewed normal function (
Snorm ), defined in
ProSpect as massfunc snorm function given by Eq 3, which can bedefined by four free parameters:SFR(age) = mSFR × e − X(age)22 (3)where X is age dependent and defined as: X (age) = (cid:18) age − mPeakmPeriod (cid:19) (cid:0) e mSkew (cid:1) asinh ( age − mPeakmPeriod ) (4)where “age” is in lookback time (Gyrs) and the fourfree parameters are defined as: • mSFR : Star formation rate normalisation in units ofM (cid:12) yr − ; • mPeak : The lookback age of the star formation history peakin Gyr; • mPeriod : Time period of the star formation i.e. the stan-dard deviation of the star formation history period; and • mSkew : The skew of the star formation history (0 for nor-mal, < > SFR is the peak SFR value andmodifying its value scales the entire plot up and downwithout changing the shape and hence it also acts as anormalisation parameter.We fit each of our adopted CSFHs with the
Snorm function and show the fitted values in Table 5 and also plot-ted as lines in Fig. 10. The solid gold line shows the rawCSFH values from Bellstedt et al. (2020b) where the yellowshaded region shows the measurement uncertainty from thesampling of the MCMC chains that determine the SFH. Ingeneral the
Snorm parametrisation is a reasonable fit to thevarious data and places all datasets onto a simple compara-ble footing, with perhaps the exception of the S´anchez et al.(2019) dataset at very low redshift. Note: the S´anchez et al.(2019) study used a Salpeter IMF, and the values have beenrescaled to a Chabrier IMF by multiplying by × With a simple representation of the CSFH in place, wecan now make our COB predictions using the
ProSpect
EBL routine (see EBL vignette on github ). In brief the ProSpect prediction adopts the Chabrier IMF (Chabrier2003) and defines the COB/CIB by condensing the full https://github.com/asgr/ProSpect galaxy population down to a canonical representative galaxyintended to reflect a mean mass-weighted SED. This canon-ical representative SED, includes the following elements: • an invariant Chabrier IMF (Chabrier 2003). • closed box stellar evolution using BC03 (Bruzual &Charlot 2003) simple stellar population libraries. • two-components dust model to represent the birthclouds and the inter-stellar medium (ISM) • a FUV-IR attenuation prescription following the freeform variant of the Charlot & Fall (2000) prescription (in-dependent for each dust component) • energy-conserving flux redistribution to the far-IR usingtemplates by Dale et al. (2014) (independent for each dustcomponent) • evolving gas-phase metallicity that tracks star-formation, i.e., grows with the locked-up stellar mass (Driveret al. 2013)In running our model for our canonical galaxy we adopt theopacities ( τ BC = 1, τ ISM = 0 . α BC = 1 . α ISM = 3 . ProSpect defaults, exceptfor the α BC value which is increased from 1.0 to 1.75 im-plying slightly cooler dust, and was chosen to shift the pre-dicted FIR dust peak to longer wavelengths to better alignwith the observed CIB measurements. With these parame-ters fixed throughout, we can now feed in the parameterisedCSFH into the ProSpect
EBL routine to produce a fullCOB/CIB prediction of the extragalactic background lightas observed at z = 0 .
0. The dashed lines in the left panelsof Fig. 11 show the resulting EBL predictions. Note that,for CSFH measurements that are normalised higher thanthe others (S´anchez et al. 2019 and Fermi-LAT Collabora-tion et al. 2018), the resulting EBL over-predicts our EBLmeasurements.For each of our six CSFH paramaterisations, we nowgenerate a grid of COB/CIB predictions, as viewed throughthe ugriZY JHK s filter-set, by fitting for a range of CSFHnormalisations (mSFRs), and final metallicities ( Z final ). Asdiscussed in detail by Bellstedt et al. (2020b), the SED ofa galaxy is influenced by not only the age distribution ofthe stellar population, but also the metallicity of the stellarpopulation. It is therefore important for us to simultane-ously fit for the modelled final gas-phase metallicity ( Z final ),in addition to the normalisation of the CSFH (mSFR). Us-ing the closed-box metallicity evolution prescription within ProSpect , the shape of the metallicity evolution is pre-scribed by the build-up of stellar mass, but the target finalmetallicity value is a free parameter. Inherently, this valuerepresents the light-weighted metallicity at which stars areforming at the present day, which could be interpreted asthe light-weighted ISM metallicity. We note that this is nota measurement of the overall “cosmic” gas-phase metallicity,as this approach is only sensitive to star-forming gas. Con-sequently, the Z final parameter can be regarded as a kindof nuisance parameter, required for fitting in order to accu-rately model the metallicities of the stellar populations inthe galaxy sample.In order to compute the revised mSFR, we simply scaleup and down the mSFR value until we get the best χ between the EBL model prediction and the EBL data in MNRAS , 1– ????
0. The dashed lines in the left panelsof Fig. 11 show the resulting EBL predictions. Note that,for CSFH measurements that are normalised higher thanthe others (S´anchez et al. 2019 and Fermi-LAT Collabora-tion et al. 2018), the resulting EBL over-predicts our EBLmeasurements.For each of our six CSFH paramaterisations, we nowgenerate a grid of COB/CIB predictions, as viewed throughthe ugriZY JHK s filter-set, by fitting for a range of CSFHnormalisations (mSFRs), and final metallicities ( Z final ). Asdiscussed in detail by Bellstedt et al. (2020b), the SED ofa galaxy is influenced by not only the age distribution ofthe stellar population, but also the metallicity of the stellarpopulation. It is therefore important for us to simultane-ously fit for the modelled final gas-phase metallicity ( Z final ),in addition to the normalisation of the CSFH (mSFR). Us-ing the closed-box metallicity evolution prescription within ProSpect , the shape of the metallicity evolution is pre-scribed by the build-up of stellar mass, but the target finalmetallicity value is a free parameter. Inherently, this valuerepresents the light-weighted metallicity at which stars areforming at the present day, which could be interpreted asthe light-weighted ISM metallicity. We note that this is nota measurement of the overall “cosmic” gas-phase metallicity,as this approach is only sensitive to star-forming gas. Con-sequently, the Z final parameter can be regarded as a kindof nuisance parameter, required for fitting in order to accu-rately model the metallicities of the stellar populations inthe galaxy sample.In order to compute the revised mSFR, we simply scaleup and down the mSFR value until we get the best χ between the EBL model prediction and the EBL data in MNRAS , 1– ???? (2021) S. Koushan et. al. the ugriZY JHK s bands. Each CIB/COB prediction is re-gressed against the COB data, the error-weighted (reduced) χ values determined for each prediction, and the final best-fit and error contours determined for each of our adoptedCSFH representations. The model results, showing the orig-inal and final COB/CIB predictions, for each of our adoptedCSFHs are shown in Fig. 11 (left panels), along with the 1 σ ,2 σ and 3 σ confidence levels (right panels).Table 5 summarises the result of the χ analysis, show-ing the four initial free parameters of the Snorm function,the final metallicity, the revised mSFR and best χ -valueafter applying the EBL constraint to each CSFH represen-tation. Note that the reduced- χ values are reasonably dis-tributed around unity, as would be expected if errors are rea-sonable and the model appropriately describes the data. Therevised Snorm values represent a range of adjustments withthe CSFHs of MD14 (Madau & Dickinson 2014) and B20(Bellstedt et al. 2020b) requiring the least re-scaling and es-sentially fully consistent with their original specified values.The remaining CSFH’s all required adjustment to match theCOB/CIB data. Figure 12 shows the set of CSFHs renor-malised, and the scatter is notably reduced around cosmicnoon from × < (cid:74) yr − Mpc − at z ∼ ∼ The top panel in Fig. 13 shows a zoom of the best fitEBL from u to K s band as obtained by ProSpect usingthe renormalised D18 CSFH. We also show the associatedprediction for the stellar mass density (SMD) based on thisCSFH (bottom panel).This SMD prediction, shown as the red curve inthe bottom panel of Fig. 13 compares well to data fromWilkins et al. (2008) (grey data points) and D18 (blue datapoints), where the latter represents the full GAMA/G10-COSMOS/3D-HST combined dataset, spanning almost theentire age of the Cosmos. In detail, we see that the modelonly just lies within the errorbars of the data which ismost likely a reflection of the fact that
ProSpect doesproduce slightly higher stellar masses, given the same data,than MAGPHYS (see Robotham et al. 2020 figure 33,lower panel) which was used to determine the datapointsin D18 (Driver et al. 2018). However, there are also twoother plausible physical explanations. One is that mass-lossis very slightly higher than assumed for a Chabrier IMF.The second and perhaps more interesting is that such adiscrepancy could potentially represent stripped or missinglocal mass, as the IGL is predominantly derived from light production at intermediate redshifts while the SMD ismeasured locally (see Ashcraft et al. 2018 and Ashcraft2018 PhD dissertation). Hence if some stellar mass isstripped from the progenitors of the present day populationin the latter half of the Universe, this could manifest asa modest discrepancy between the EBL predicted SMDand that measured. The obvious way to distinguish, is bybreaking the EBL into its constituent CSED time slices,where one might start to identify less energy than expectedif mass is being stripped and not accounted for. The reverseis also true, in that a good match between the CSFH,SMD and COB/CIB data suggests, mass stripping is a verysmall factor (few percent), otherwise the predicted SMDwould significantly over-predict the measured stellar-massdensity. To bring our predicted SMD perfectly in line withthe data, would require and offset of about 20 percent,which is consistent with the 15 per cent stellar massoffset between MAGPHYS and
ProSpect (see Robothamet al. 2020; Bellstedt et al. 2020b). The implication, isthat our results are consistent with mass stripping ofanywhere from 0-10 per cent. This constraint is likely to betightened as we improve both our stellar mass estimates ofthe nearby population, as well as our COB/CIB constraints.
Finally, Fig. 14 presents our best-fit EBL model based onthe modified D18 (Driver et al. 2018) CSFH, to the fullCOB and CIB data, shown as black data points. Notethat the new IGL data from our revised number counts( ugriZY JHK s ) are shown as closed black circles, andthis is now complemented by our previous IGL resultsDriver et al. (2016b) shown as open black circles to extendthe plot to the UV, MIR and FIR regimes. Overlaid arerecent model predictions from our earlier phenomenologicalmodel (purple line; Andrews et al. 2018), data from thenew SHARK semi-analytic model (green line; Lagos et al.2019), and the recent Khaire & Srianand phenomenologicalmodel (orange line; Khaire & Srianand 2019). All providereasonable descriptions of the observed distribution, withthe optimised ProSpect fit providing the closest matchover not just the u − K s regime but the entire FUV to FIRrange. We provide our optimised model in electronic Tableform, with Table 6 showing the first and last three lines ofthe data file, and advocate the use of these data as the bestcurrent representation of the COB and CIB. A related pointto consider is that the comparison between SHARK andthe observation at the FIR likely requires further analysis.The cyan dashed line shows the EBL calculation fromthe EAGLE data . We compute the predicted EBL fromthe EAGLE simulations (Schaye et al. 2015) by retrievingfrom the EAGLE public database (McAlpine et al. 2016)all galaxies from z=0 to z=9 that had an observer-framer-band flux > Data obtained from the public database (McAlpine et al. 2016)MNRAS , 1– ?? (2021) onstraining the CSFH from the improved EBL Figure 10.
The cosmic star-formation rate before our constraint using mSFR (the peak value) parameter (as described in text) shownfor the various methods: SED core sampling in sky-blue (Driver et al. 2018); VHE constraints in gray (Fermi-LAT Collaboration et al.2018); SED forensic reconstruction in gold (Bellstedt et al. 2020b); IFU forensic reconstruction in purple and green (L´opez Fern´andezet al. 2018 and S´anchez et al. 2019 respectively). We also show the curve from Madau & Dickinson (2014) as the dashed red line. Thesolid gold line shows the raw CSFH value from Bellstedt et al. (2020b).
Table 5.
Free parameters obtained from the
Snorm function and our best SFR and metallicity values for six different methods beforeand after rescaling. Method m SFR final m SFR initial m Peak m Period m Skew Z best chi best Driver et al. (2018); D18 0.069 ± ± ± ± ± ± ± ± ± ± ± ±
100 and 160 microns, SPIRE 250, 350 and 500 microns,and SCUBA2 450, and 850 microns. These fluxes werecomputed via postprocessing of EAGLE galaxies using theradiative transfer code SKIRT (see Trayford et al. 2017for details). For each redshift, we compute the impliedprojected sky area of a (100 cMpc) square assuming theEAGLE cosmology (Planck Collaboration et al. 2014) andthe energy contributed by that redshift, E ( z ), as the sumof the fluxes of all galaxies at that redshift divided by thesolid angle corresponding to the area above. We then inte-grate under the E(z) vs. z curve to obtain the predicted EBL. We note that the fluxes used to compute the pre-dicted EBL for EAGLE, only include galaxies that havea number of (sub-)particles representing the galaxy’s bodyof dust >
250 (see Camps et al. 2018 for details). Thisroughly ignores about ≈
30% of galaxies with stellar masses > M (cid:12) at z < .
4, which are expected to be bright inthe optical and NIR given their large stellar mass. We henceconsider the presented EBL here as a lower limit for EA-GLE, but it’s likely that the true EBL in the simulation isonly 20 −
30% higher in the optical-NIR, which is still ≈ MNRAS , 1– ????
30% higher in the optical-NIR, which is still ≈ MNRAS , 1– ???? (2021) S. Koushan et. al.
Figure 11.
The reduced χ analysis for each of our star-formation histories showing the initial and final EBL model predictions comparedto the data, and (right panels) the 1, 2 and 3 σ contours in the mSFR-Zfinal plane. The y-axis shows the final metallicities ( Z final ). Thisanalysis has been summarised in Table 5. MNRAS , 1– ?? (2021) onstraining the CSFH from the improved EBL Figure 12.
The cosmic star-formation rate after our constraint.
Table 6.
A sample of our EBL spectrum from
ProSpect . Thefull table in csv format ( filename : EBL model.csv ) is providedelectronically as supplementary material. The first column iden-tify the wavelength in angstrom for a range of FUV to FIRand the second column shows the derived flux in the unit of( nW/m sr − ).Wavelength (˚A) Flux ( nW/m sr − )1122.0185 0.23051148.1536 0.27531174.8976 0.3160... ...9549925.9 0.64879772372.2 0.604610000000.0 0.5636 hence the under-prediction here is in tension with observa-tions. Furthermore, there are reasons to believe the simula-tions are underestimates by 20-30 per cent. However even ifone incorporated this correction the simulations would stillunder-predict the EBL observations by 50 per cent. This isnot necessarily surprising given the discrepancy in the pre-dicted number counts of EAGLE and observations in theFIR bands (e.g. Cowley et al. 2019; Wang et al. 2019). This paper presents a measurement of the optical part ofthe extra-galactic background light via integrated galaxynumber-counts to significantly greater accuracy than pre-vious measurements, enabling us to use the COB data toplace constraints onto the CSFH. The main results of thiswork are summarized as: • We combined three complementary datasets and de-rived in relatively wide and deep fields the galaxy countsin the optical and NIR bands ( ugriZY JHK s ) from theGAMA, DEVILS and various deep ground-based and HSTdatasets. The compilation of wide and deep datasets allowus to measure the galaxy number counts in a wide magni-tude interval (spanning more than 15 magnitudes) acrossthe multi-wavelength filters. Where possible filter conver-sions have been applied to map our measurements onto theESO VST/VISTA ugriZY JHK s system. In general our re-vised number-counts agree well with other published datain all bands and with significantly improved errors over thecritical intermediate magnitude range that dominates theIGL signal. • From the number-count data we determine the contri-bution in each magnitude interval to the total luminositydensity. We fit a smooth spline and extrapolate over the en-tire magnitude range to derive the integrated flux density,
MNRAS , 1– ????
MNRAS , 1– ???? (2021) S. Koushan et. al.
Figure 13. (top) the best fit EBL model from
ProSpect ;(bottom) the derived stellar mass density compared to datafrom Wilkins et al. (2008) (gray) and Driver et al. (2018) (blue). or IGL, of discrete sources from a steradian of sky over thefull path-length of the universe and in each filter. Comparedto previous results we find a 5-15 per cent increase in ourIGL measurements from u − K s . We attribute this increaseto several factors which include: ProFound ’s propensity torecover closer to total fluxes, appropriate filter transforma-tions, and improved accounting for the area lost to brightstars and their associated halos (ghosting). • We compared our optical/NIR IGL measurements toVHE COB (i.e., those from Fermi-LAT, MAGIC, H.E.S.S.,and VERITAS), and direct estimates from space platforms (HST, Pioneer 10/11, New Horizons). Our derived opti-cal/NIR IGL results agree well within the errors with allVHE COB measurements, and those from the deep spacebased missions (i.e., Pioneer 10/11 and New Horizons). Theimplication is that significant uncertainty must exist in ei-ther or both our understanding of dust in the inner Solarsystem (zodiacal light), or upper atmosphere Earth-shine inthe HST dataset. • The new optical/NIR IGL measurements include an im-proved error-analysis to manage the uncertainty introducedby cosmic variance. We assessed both systematic (zero-point, filter uncertainty, and cosmic variance), and random(spline fitting and Poisson) errors. In doing so we have re-duced the uncertainty from ∼
20 per cent to below ∼ • To model the EBL, we introduced the
Snorm functionfitted to the data of Driver et al. (2018). This is a simplefour-parameters function that describes the adopted CSFHsextremely well. Using
ProSpect we provide a predictionof the EBL and regress against our data to determine theoptimal normalisation of the
Snorm function and the gas-phase metallicity evolution. The model provides an excel-lent fit from UV to far-IR and places strict constraints onthe normalisation of the cosmic star-formation history witha constraint of 0 . − . M (cid:74) yr − Mpc − at cosmicnoon. The constraint is in close agreement with that pro-vided by Madau & Dickinson. We conclude that both therevised Snorm function and the Madau & Dickinson (2014)result provide an excellent description of the CSFH whichis fully consistent with both the stellar mass growth andenergy production over all time.
ACKNOWLEDGEMENTS
MNRAS , 1– ?? (2021) onstraining the CSFH from the improved EBL Figure 14.
Final EBL model from Driver et al. (2018) in sky blue is compared to the other data from semi-analytical model (green),synthesis model (orange) and phenomenological model (purple). We overplotted the IGL data from u to K s obtained from our numbercounts in closed circles combined with the IGL data for UV, mid-IR and far-IR from Driver et al. (2016b) in open circles to present bothCIB and COB. EAGLE data is shown in cyan. tralian Research Council and the participating institutions.The DEVILS website is https://devilsurvey.org. The DEV-ILS data are hosted and provided by AAO Data Central .Parts of this research were conducted by the AustralianResearch Council Centre of Excellence for All Sky Astro-physics in 3 Dimensions (ASTRO 3D) through project num-ber CE170100013.This research is supported by the UPA awarded by the Uni-versity of Western Australia Scholarships Committee andthe awarded by the Astronomical Society of Australia.This work was supported by resources provided by thePawsey Supercomputing Centre with funding from the Aus-tralian Government and the Government of Western Aus-tralia.MS has been supported by the European Union’s Horizon2020 Research and Innovation programme under the MariaSklodowska-Curie grant agreement (No. 754510), the Pol-ish National Science Centre (UMO-2016/23/N/ST9/02963),and the Spanish Ministry of Science and Innovation throughthe Juan de la Cierva-formacion programme (FJC2018-038792-I). https://datacentral.org.au/ The data underlying this article are derived from threedistinct sources as described below:1- The GAMA catalogue (
GAMAKiDSVIKINGv01.fits
DEVILSProFoundPhotomv01 ) is currently available forinternal team use and will be presented in Davies at alin prep, and subsequently publicly releases through datacentral (https://datacentral.org.au/).3- The HST data is obtained from Driver et al. (2016b).All the data products that are used and/or derived in thiswork are available on request.
REFERENCES
Abeysekara A. U., et al., 2019, ApJ, 885, 150Aharonian F., et al., 2006, Nature, 440, 1018Ahnen M. L., et al., 2016, A&A, 590, A24Aihara H., et al., 2018, PASJ, 70, S4Andrews S. K., Driver S. P., Davies L. J. M., Lagos C. d. P.,Robotham A. S. G., 2018, MNRAS, 474, 898Ashcraft T. A., et al., 2018, PASP, 130, 064102MNRAS , 1– ????
Abeysekara A. U., et al., 2019, ApJ, 885, 150Aharonian F., et al., 2006, Nature, 440, 1018Ahnen M. L., et al., 2016, A&A, 590, A24Aihara H., et al., 2018, PASJ, 70, S4Andrews S. K., Driver S. P., Davies L. J. M., Lagos C. d. P.,Robotham A. S. G., 2018, MNRAS, 474, 898Ashcraft T. A., et al., 2018, PASP, 130, 064102MNRAS , 1– ???? (2021) S. Koushan et. al.
Bellstedt S., et al., 2020a, Galaxy And Mass Assembly(GAMA): Assimilation of KiDS into the GAMA database( arXiv:2005.11215 )Bellstedt S., et al., 2020b, arXiv e-prints, p. arXiv:2005.11917Bernstein R. A., 1999, HST/LCO Measurements of the OpticalExtragalactic Background Light. p. 487Bernstein R. A., 2007, The Astrophysical Journal, 666, 663Bernstein R. A., Freedman W. L., Madore B. F., 2002, ApJ, 571,85Biteau J., 2013, in Cambresy L., Martins F., Nuss E., PalaciosA., eds, SF2A-2013: Proceedings of the Annual meeting of theFrench Society of Astronomy and Astrophysics. pp 303–312( arXiv:1310.2989 )Biteau J., Williams D. A., 2015, ApJ, 812, 60Bravo M., Lagos C. d. P., Robotham A. S. G., Bellstedt S.,Obreschkow D., 2020, arXiv e-prints, p. arXiv:2003.11258Bruzual G., Charlot S., 2003, MNRAS, 344, 1000Camps P., et al., 2018, ApJS, 234, 20Capak P., et al., 2007, ApJS, 172, 99Chabrier G., 2003, PASP, 115, 763Charlot S., Fall S. M., 2000, ApJ, 539, 718Cooray A., et al., 2009, in astro2010: The Astronomy and Astro-physics Decadal Survey. ( arXiv:0902.2372 )Cooray A., et al., 2012, in Tuffs R. J., Popescu C. C.,eds, IAU Symposium Vol. 284, The Spectral EnergyDistribution of Galaxies - SED 2011. pp 482–488,doi:10.1017/S1743921312009672Cowley W. I., Lacey C. G., Baugh C. M., Cole S., Frenk C. S.,Lagos C. d. P., 2019, MNRAS, 487, 3082Dale D. A., Helou G., Magdis G. E., Armus L., D´ıaz-Santos T.,Shi Y., 2014, ApJ, 784, 83Davies L. J. M., et al., 2018, preprint, ( arXiv:1806.05808 )De Pontieu B., et al., 2014, Sol. Phys., 289, 2733Desai A., et al., 2017, ApJ, 850, 73Dole H., et al., 2006, A&A, 451, 417Dom´ınguez A., et al., 2011, MNRAS, 410, 2556Driver S., 2020, in The Build-Up of Galaxies through MultipleTracers and Facilities. p. 7, doi:10.5281/zenodo.3756434Driver S. P., GAMA Team 2016, Galaxy And Mass Assembly(GAMA): A Study of Energy, Mass, and Structure (1 kpc-1Mpc) at z<0.3. p. 269Driver S. P., Robotham A. S. G., 2010, MNRAS, 407, 2131Driver S. P., et al., 2011, MNRAS, 413, 971Driver S. P., et al., 2012, MNRAS, 427, 3244Driver S. P., Robotham A. S. G., Bland-Hawthorn J., Brown M.,Hopkins A., Liske J., Phillipps S., Wilkins S., 2013, MNRAS,430, 2622Driver S. P., et al., 2016a, MNRAS, 455, 3911Driver S. P., et al., 2016b, ApJ, 827, 108Driver S. P., et al., 2018, MNRAS, 475, 2891Edge A., Sutherland W., Kuijken K., Driver S., McMahon R.,Eales S., Emerson J. P., 2013, The Messenger, 154, 32Fermi-LAT Collaboration et al., 2018, Science, 362, 1031Fontana A., et al., 2014, A&A, 570, A11Frayer D. T., et al., 2009, AJ, 138, 1261Gaia Collaboration et al., 2016, A&A, 595, A2Giavalisco M., et al., 2004, ApJ, 600, L93Gonz´alez-Fern´andez C., et al., 2018, MNRAS, 474, 5459Grazian A., et al., 2009, A&A, 505, 1041Grogin N. A., et al., 2011, ApJS, 197, 35HESS Collaboration et al., 2018, MNRAS, 476, 4187Hauser M. G., Dwek E., 2001, ARA&A, 39, 249Hill R., Masui K. W., Scott D., 2018, Applied Spectroscopy, 72,663Hopkins A. M., et al., 2013, MNRAS, 430, 2047Inoue Y., Inoue S., Kobayashi M. A. R., Makiya R., Niino Y.,Totani T., 2013, ApJ, 768, 197Jarvis M. J., et al., 2013, MNRAS, 428, 1281 Kaiser N., et al., 2002, in Tyson J. A., Wolff S., eds,Proc. SPIEVol. 4836, Survey and Other Telescope Technolo-gies and Discoveries. pp 154–164, doi:10.1117/12.457365Kashlinsky A., 2006, New Astron. Rev., 50, 208Kawara K., Matsuoka Y., Sano K., Brand t T. D., Sameshima H.,Tsumura K., Oyabu S., Ienaka N., 2017, PASJ, 69, 31Keenan R. C., Barger A. J., Cowie L. L., Wang W. H., 2010, ApJ,723, 40Khaire V., Srianand R., 2019, MNRAS, 484, 4174Kim M. G., Matsumoto T., Lee H. M., Jeong W.-S., Tsumura K.,Seo H., Tanaka M., 2019, PASJ, 71, 82Kuijken K., et al., 2019, A&A, 625, A2Lagache G., Abergel A., Boulanger F., D´esert F. X., Puget J. L.,1999, A&A, 344, 322Lagos C. d. P., et al., 2019, MNRAS, 489, 4196Lauer T. R., et al., 2020, arXiv e-prints, p. arXiv:2011.03052Leinert C., Richter I., Pitz E., Planck B., 1981, A&A, 103, 177Liske J., et al., 2015, MNRAS, 452, 2087L´opez Fern´andez R., et al., 2018, A&A, 615, A27MAGIC Collaboration et al., 2008, Science, 320, 1752Madau P., Dickinson M., 2014, ARA&A, 52, 415Madau P., Pozzetti L., 2000, MNRAS, 312, L9Matsumoto T., et al., 2005, ApJ, 626, 31Matsuoka Y., Ienaka N., Kawara K., Oyabu S., 2011, ApJ, 736,119Mattila K., V¨ais¨anen P., 2019, Contemporary Physics, 60, 23Mattila K., Lehtinen K., V¨ais¨anen P., von Appen-Schnur G.,Leinert C., 2012, in Tuffs R. J., Popescu C. C., eds, IAUSymposium Vol. 284, The Spectral Energy Distributionof Galaxies - SED 2011. pp 429–436 ( arXiv:1111.6747 ),doi:10.1017/S174392131200957XMattila K., Lehtinen K., V¨ais¨anen P., von Appen-Schnur G., Lein-ert C., 2017a, MNRAS, 470, 2133Mattila K., V¨ais¨anen P., Lehtinen K., von Appen-Schnur G., Lein-ert C., 2017b, MNRAS, 470, 2152McAlpine S., et al., 2016, Astronomy and Computing, 15, 72McCracken H. J., et al., 2012, A&A, 544, A156Moffett A. J., et al., 2018, VizieR Online Data Catalog, p.J/MNRAS/462/4336Penzias A. A., Wilson R. W., 1965, ApJ, 142, 419Pierre M., et al., 2004, J. Cosmology Astropart. Phys., 2004, 011Pilbratt G. L., et al., 2010, A&A, 518, L1Planck Collaboration et al., 2014, A&A, 571, A11Puget J. L., Abergel A., Bernard J. P., Boulanger F., BurtonW. B., Desert F. X., Hartmann D., 1996, A&A, 308, L5Rafelski M., et al., 2015, AJ, 150, 31Robotham A. S. G., Davies L. J. M., Driver S. P., Koushan S.,Taranu D. S., Casura S., Liske J., 2018, MNRASRobotham A. S. G., Bellstedt S., Lagos C. d. P., Thorne J. E.,Davies L. J., Driver S. P., Bravo M., 2020, arXiv e-prints, p.arXiv:2002.06980S´anchez S. F., et al., 2019, MNRAS, 482, 1557Sanders D. B., et al., 2007, ApJS, 172, 86Schaye J., et al., 2015, MNRAS, 446, 521Scott D., 2000, Cosmic Glows: A CMB Review. p. 403Scoville N., et al., 2007, ApJS, 172, 1Taniguchi Y., et al., 2007, ApJS, 172, 9Teplitz H. I., et al., 2013, AJ, 146, 159Totani T., Yoshii Y., Iwamuro F., Maihara T., Motohara K., 2001,ApJ, 550, L137Trayford J. W., et al., 2017, MNRAS, 470, 771VERITAS Collaboration et al., 2011, Science, 334, 69Vernet J., D’Odorico S., Christensen L., Dekker H., Mason E.,Modigliani A., Moehler S., 2009, The Messenger, 138, 4Virani S. N., Treister E., Urry C. M., Gawiser E., 2006, AJ, 131,2373Wang L., Pearson W. J., Cowley W., Trayford J. W., B´etherminMNRAS , 1– ?? (2021) onstraining the CSFH from the improved EBL M., Gruppioni C., Hurley P., Micha(cid:32)lowski M. J., 2019, A&A,624, A98Wilkins S. M., Trentham N., Hopkins A. M., 2008, MNRAS, 385,687Windhorst R. A., et al., 2011, ApJS, 193, 27Windhorst R. A., et al., 2018, ApJS, 234, 41Wright E. L., 2004, New Astron. Rev., 48, 465Wright E. L., et al., 2010, AJ, 140, 1868Xu C. K., et al., 2005, ApJ, 619, L11Zemcov M., et al., 2013, ApJS, 207, 31Zemcov M., Immel P., Nguyen C., Cooray A., Lisse C. M., PoppeA. R., 2017, Nature Communications, 8, 15003MNRAS , 1– ????