Galaxy And Mass Assembly: the evolution of the cosmic spectral energy distribution from z = 1 to z = 0

We present the evolution of the Cosmic Spectral Energy Distribution (CSED) from $z = 1 - 0$. Our CSEDs originate from stacking individual spectral energy distribution fits based on panchromatic photometry from the Galaxy and Mass Assembly (GAMA) and COSMOS datasets in ten redshift intervals with completeness corrections applied. Below $z = 0.45$, we have credible SED fits from 100 nm to 1 mm. Due to the relatively low sensitivity of the far-infrared data, our far-infrared CSEDs contain a mix of predicted and measured fluxes above $z = 0.45$. Our results include appropriate errors to highlight the impact of these corrections. We show that the bolometric energy output of the Universe has declined by a factor of roughly four -- from $5.1 \pm 1.0$ at $z \sim 1$ to $1.3 \pm 0.3 \times 10^{35}~h_{70}$~W~Mpc$^{-3}$ at the current epoch. We show that this decrease is robust to cosmic variance, SED modelling and other various types of error. Our CSEDs are also consistent with an increase in the mean age of stellar populations. We also show that dust attenuation has decreased over the same period, with the photon escape fraction at 150~nm increasing from $16 \pm 3$ at $z \sim 1$ to $24 \pm 5$ per cent at the current epoch, equivalent to a decrease in $A_\mathrm{FUV}$ of 0.4~mag. Our CSEDs account for $68 \pm 12$ and $61 \pm 13$ per cent of the cosmic optical and infrared backgrounds respectively as defined from integrated galaxy counts and are consistent with previous estimates of the cosmic infrared background with redshift.


INTRODUCTION
The Cosmic Spectral Energy Distribution (CSED, e.g.Driver et al. 2008Driver et al. , 2016b;;Domínguez et al. 2011;Gilmore et al. 2012;Somerville et al. 2012), describes the total energy generated as a function of wavelength for a cosmologically representative volume at some specified time.This is different to the photon budget, which describes the photons passing through the same volume at that time, and the extragalactic background light (EBL, e.g.Dwek et al. 1998;Hauser & Dwek 2001), which describes the radiation received per unit solid angle.One can define the resolved CSED which arises from discrete sources, and the total CSED which includes both discrete and diffuse light.This work measures the resolved CSED, which we will refer to as simply "the CSED" unless indicated.
The CSED encodes statistical information about the ongoing processes of galaxy formation and evolution.This link becomes obvious when one considers the (resolved) CSED to be the sum of the individual spectral energy distributions (SEDs) of all galaxies in some cosmologically representative volume.The optical to near-infrared portion of a galaxy's SED encodes information about stellar mass, both gas and stellar phase metallicity, and the ages of stellar populations (e.g.Taylor et al. 2011;Madau & Dickinson 2014), while the ultraviolet and total infrared emission can be used to estimate a galaxy's current star formation rate (e.g.Kennicutt 1998;Madau & Dickinson 2014;Davies et al. 2016).About half the energy produced by stars is attenuated by dust within their galaxy and re-radiated in the mid-and far-infrared, with the shape and magnitude of a galaxy's far-infrared SED depending on dust mass, temperature, geometry and composition (e.g.da Cunha, et al. 2008;Dunne et al. 2011;Symeonidis et al. 2013).Finally, emission from an active galactic nucleus can also shape a galaxy's SED, and this contribution becomes increasingly significant at higher redshifts (see e.g.Richards et al. 2006).Analysing the total CSED allows for the extraction of the population average, weighted by energy and number density, of key quantities within the specified volume.Tracing the evolution of the CSED with lookback time therefore allows a reconstruction of the evolution of the energy output from stars, dust and AGN.
The integrated galactic light (IGL), which represents the resolved component of the EBL, can be determined from the redshifted CSED using a volume-weighted integral over all time.To derive the EBL from the IGL, one should consider additional contributions from diffuse radiation from the epoch of reionization (Cooray et al. 2012a) and faint intra-halo light (Cooray et al. 2012b;Zemcov et al. 2014).
Recent measurements (H.E.S.S. Collaboration 2013; Ahnen et al. 2016;Driver et al. 2016c) suggest the diffuse components of the EBL represent an excess of approximately 20 per cent over the IGL in the near-infrared.This excess is marginally significant given uncertainties in the IGL (dominated by cosmic sample variance) and EBL measurements.
The EBL and IGL have historically received more interest than the full, panchromatic CSED (e.g.Partridge & Peebles 1967a,b;Franceschini et al. 2008;Finke et al. 2010;Domínguez et al. 2011;Inoue et al. 2013;Driver et al. 2016c).Multiple groups have also examined how the cosmic infrared background at specific wavelengths builds up as a function of redshift, which if combined is equivalent to deriving the far-infrared portion of the CSED (e.g.Marsden et al. 2009;Jauzac et al. 2011;Berta et al. 2011;Béthermin et al. 2012b;Viero et al. 2013).Measurements of the CSED -an instantaneous quantity, as opposed to an integrated measurement over all of Universal history -at multiple epochs have a greater ability to constrain cosmological models of galaxy formation.
Multiple groups (e.g.Béthermin et al. 2010;Carniani et al. 2015;Fazio et al. 2004;Gardner et al. 2000;Madau & Pozzetti 2000;Xu et al. 2005) have measured the IGL, but these efforts are generally restricted to one region of the electromagnetic spectrum (Driver et al. 2016c).To study the CSED from the far-ultraviolet to the far-infrared requires a combination of sufficiently deep and wide multiwavelength imaging and spectroscopic datasets and a means of deriving consistent photometry across the wavelength range.With the aid of SED fitting tools, the CSED and IGL can be characterized over the wavelength range while simultaneously avoiding systematic errors induced by inhomogenous data reduction methods.
The Galaxy And Mass Assembly survey (Driver et al. 2011;Liske et al. 2015) is ideally suited to measuring the recent time evolution of the CSED.GAMA is an integrated multiwavelength imaging and spectroscopic campaign to examine the distribution of matter and energy on kpc to Mpc scales in the low-redshift (z < 0.25) Universe.GAMA covers 180 square degrees of the equatorial sky in three fields to a spectroscopic completeness of ∼98 per cent.An intermediate redshift (0.5 < z < 1) analogue to GAMA (in sample size and wavelength coverage) was assembled by Davies et al. (2015) and Andrews et al. (2017) using existing public data in the Cosmological Origins Survey (COSMOS; Scoville et al. 2007) field.This work aims to characterise the evolution in the CSED using a combination of the GAMA and COSMOS datasets.
Three previous measurements of the CSED at low redshifts (z < 0.2) have been made using the GAMA dataset.Driver et al. (2012) calculated the CSED using the GAMA far-ultraviolet to Ks luminosity functions, and extrapolated it to the far-infrared using the Dale & Helou (2002) models of dust attenuation.Kelvin et al. (2014) measured the CSED as a function of galaxy morphology for z < 0.06.More recently, Driver et al. (2016b) measured the CSED and its evolution over the redshift range 0 < z < 0.2 by fitting SEDs to a cosmologically representative sample of galaxies.
Here we improve on the latter CSED measurement in two respects: i) using an updated photometric catalogue (Wright et al. 2016a) with improved deblending and reduced gross photometric errors, and ii) expanding the redshift range using re-reduced publicly available data in the COSMOS field (Andrews et al. 2017).These CSEDs complement the recent IGL measurement by Driver et al. (2016c).In Section 2, we describe the multiwavelength data sets, techniques used to derive photometry and the SED modelling techniques we use to interpolate between these photometric measurements.Then in Section 3, we determine the CSED and examine its reliability.In Section 4, we present concluding remarks.

GAMA equatorial regions
The GAMA (Driver et al. 2011;Liske et al. 2015) spectroscopic campaign targeted 180 square degrees of the equatorial sky in three fields, centered on R.A.=9h, 12h and 14.5h, using the AAOmega spectrograph mounted on the 3.9 m Anglo-Australian Telescope.GAMA obtained redshifts for ∼200k galaxies to a depth of r < 19.8 mag (SDSS Petrosian).This is complemented by ultraviolet imagery from the Galaxy Evolution Explorer (GALEX ; Martin et al. 2005), optical data from the Sloan Digital Sky Survey (SDSS DR7; Abazajian et al. 2009), near-infrared data from VISTA (Visible and Infrared Telescope for Astronomy) Kilodegree Infrared Galaxy survey (VIKING; Edge et al. 2013), mid-infrared data from the Wide-field Infrared Survey Explorer (WISE ; Wright et al. 2010) and far-infrared data from Herschel -Atlas (Eales et al. 2010), as summarized by Driver et al. (2016b).The combined grasp of the GAMA filters is shown in Figure 1 (upper).Wright et al. (2016a) implemented a novel program, the Lambda-Adaptive MultiBand Deblending Algorithm in R (lambdar), which is capable of deriving (forced) aperture matched photometry from the far-ultraviolet to the farinfrared.lambdar is explicitly designed to correctly deblend objects in the far-infrared, where resolution and sensitivity are the lowest.Particular care was taken in building a set of robust aperture definitions in order to obtain an optimal photometric solution.Deblended fluxes are then obtained using lambdar across all 21 bands using these aperture definitions convolved with the point spread function.This approach minimises the potential for gross photometric inconsistencies, including table and aperture mismatches.The resulting catalogue, LDRSummaryPhotometryv01, contains consistent flux and error measurements for approximately 220,000 sources across the GAMA wavelength range.Wright et al. (2016a) demonstrate that using this catalogue leads to clear improvements in SED fits and star formation rate estimators over table matching.

G10/COSMOS
The G10 region (R.A. = 149.55• -150.65 • , Dec = +1.80-+2.73• ) is a subset of the broader COSMOS region chosen for its relatively high spectroscopic completeness (∼ 45% for extra-galactic sources with i + < 22 mag).It forms an intermediate-redshift comparison to GAMA.Davies et al. (2015) combined existing spectra from the PRIsm MUltiobject Survey (PRIMUS; Coil et al. 2011;Cool et al. 2013), the VIMOS-VLT Deep Survey (VVDS; Garilli et al. 2008) and SDSS DR10 (Ahn et al. 2014), re-reduced zCOSMOS (Lilly et al. 2007(Lilly et al. , 2009) ) spectra and photometric redshifts (Ilbert et al. 2009), obtaining reliable spectroscopic redshifts for over 22000 sources.This is complemented by publicly available ultraviolet images from GALEX (Zamojski et al. 2007), optical images from the Canada-France-Hawaii and Subaru telescopes (Capak et al. 2007;Taniguchi et al. 2007Taniguchi et al. , 2015) ) and near-infrared data from UltraVISTA (Mc-Cracken et al. 2012).In the mid-infrared, data is available from the Spitzer survey in COSMOS (S-COSMOS; Sanders et al. 2007) and the Spitzer Large Area Survey with Hy-perSuprimeCam (SPLASH; Capak et al. in prep), while farinfrared data has been published by PACS (Photodetector Array Camera and Spectrometer) Evolution Probe (PEP; Lutz et al. 2011) and the Herschel Multitier Extragalactic Survey (HerMES; Oliver et al. 2012).The combined filter curve for these datasets is shown in Figure 1 (lower).Andrews et al. (2017) used an equivalent procedure to Wright et al. (2016a) to derive consistent flux measurements in 38 bands and photometric redshifts from this dataset for approximately 186,000 sources brighter than i + < 25 mag.Andrews et al. (2017) demonstrated that using their catalogues results in goodness of SED fits being equivalent to or better than those derived from the recent COSMOS2015 catalogue (Laigle et al. 2016).et al. (2017) fitted SEDs to the GAMA and G10/COSMOS data using the spectral analysis code magphys (Multi-wavelength Analysis of Galaxy Physical Properties;da Cunha, et al. 2008).These models are computed using the Bruzual & Charlot (2003) stellar libraries, a Chabrier (2003) initial mass function, and the Charlot & Fall (2000) prescription for dust attenuation which consists of a two-component description of the interstellar medium (stellar birth clouds and diffuse interstellar medium).The dust emission is computed via energy balance -the energy absorbed by dust in the ultraviolet to near-infrared range is reemitted in the mid-to far-infrared range using empiricallycalibrated dust emission components: polycyclic aromatic hydrocarbons, warm dust emitting in the mid-infrared, and two components of dust in thermal equilibrium with varying temperatures.magphys was modified to use the upper limits in the Andrews et al. (2017) and Wright et al. (2016a) catalogues, increase the range of possible dust masses and increase the photometric error floor from 5 to 10 per cent.Driver et al. (2017) demonstrate that the SED fits are generally robust and obtain stellar and dust mass densities and cosmic star formation rates from z = 5 to z = 0 in line with literature estimates.However, the fitting procedure and underlying imaging data may leave systematic impacts on measurements of the CSED:

Driver
• The magphys templates do not, as yet, incorporate AGN emission.This affects predominantly the higher redshift end of the sample.
• Dust properties are poorly constrained or extrapolated for a significant portion of the combined sample due to the low depth and resolution of the far-infrared data.This is especially important for the high-redshift end of the sample, where many objects only have adoped upper limits in the far-infrared.
• The GAMA sample does not have 70 µm imaging.This results in a near-complete inability to constrain the warm dust properties of individual galaxies.
We address these potential systematic effects in Section 3.3.
From the GAMA catalogue, we select all galaxies with good spectroscopic redshifts (NQ > 2 from SpecObjv27 ) and coverage in FUV, NUV and 250 µm, which play an important role in constraining star formation and dust properties.This reduces the sample area to 129 0 0.2 0.4 0.6 0.8 1 10 When using a flux-limited sample to estimate the resolved CSED, one inevitably misses objects that fall below either the apparent magnitude or surface brightness limits of the survey.As a consequence, the raw CSEDs derived from simply stacking the SEDs are not comparable across epochs.Firstly, each redshift bin requires a 1/Vmax correction to correct for incompleteness in the yellow shaded areas of Figure 2. We compute these from the r or i + absolute magnitude of each galaxy assuming limiting apparent magnitudes of r = 19.8 and i + = 25.0 mag for GAMA and G10 respectively.We then compute a (luminosity) distance and corresponding volume Vmax for each object beyond which it should no longer be observable given these limits.Each object is then assigned a weight wi = (Vu − V l )/(Vmax − V l ) where Vu and V l correspond to the upper and lower edges of each redshift bin.These weights are capped at 10 to stop a single galaxy on the boundary of a particular redshift bin being upscaled to dominate emission in that bin.The (restframe) CSED (λ) can then be derived by simply computing: where SEDi(λ) is the rest frame best fit SED of a galaxy in λf λ units.This Vmax correction biases the contribution to the CSED from lower mass systems by overweighting observed galaxies in a given redshift bin to compensate for systems below the apparent magnitude limit.This effect is mitigated by the small duration in lookback time (approximately 750 Myr each) of our redshift bins.Furthermore, in a sufficiently deep sample these systems only represent a small contribution to the total luminosity.This effect can only be addressed with deeper surveys.Secondly, each redshift bin samples a different range of stellar masses -higher redshift bins are Malmquist biased towards high mass galaxies and will not include systems in the red shaded areas of Figure 2.This incompleteness affects both the shape and normalization of the CSED because the CSED shape is mass dependent.
To correct for the mass limit bias, we compute a total optical luminosity (Lopt) for each galaxy by integrating the respective SED fit between 100 nm < λ < 8 µm and construct a corresponding Lopt distribution (histogram) as shown in the top panels of Figure 3.We then compute the successive contributions of each bin to the CSED (bottom panel) by multiplying by Lopt.Finally we fit a 10 point spline (green curve), weighted by the inverse fractional error squared, to where the energy density is complete and reliable.The fitting limit is determined by eye (see dotted lines on Figure 3).The spline fit is then extrapolated outside this range.Manual fitting helps ensure that the gradient of the extrapolated Lopt function is reasonable, reducing the impact of small number statistics at the high-mass end where the gradient is rapidly changing (see e.g. the 0.08 < z < 0.14 and 0.45 < z < 0.56 bins).The ratio of the integral under the spline to the total CSED (the integral under the red/blue lines) gives a redshift-independent correction factor as depicted in Figure 3.While the spline fit is bound at both ends, the integration is performed from 10 34 to 10 39 W (reflecting the extent of the data and to reduce extrapolation) to reduce the impact of error in the faint/bright end slope on the CSED.In the lowest redshift bin, extending the correction to 10 33 W or 10 30 W yields only 1.4 and 2.2 per cent extra flux respectively.The small bump in the Lopt distribution for the G10 in combined redshift bins is likely to be a systematic effect in either the photometric redshifts or SED fitting -no such bump exists in the spectroscopic redshift only GAMA sample; we adjust the fitting range accordingly.Using a full bolometric luminosity (100 nm < λ < 1 mm) correction would be ideal for avoiding bias, however the quality of the far-infrared data does not permit this.
We also define three bins -0.20 < z < 0.28, 0.28 < z < 0.36 and 0.36 < z < 0.45 -where we combine the two samples (blue areas in Figure 2).The (Vmax corrected) samples are spliced such that the CSED is the sum of GAMA objects whose Lopt > 10 37 , 10 37.2 and 10 37.6 W (with increasing redshift) and G10 objects below these thresholds.The cutoff corresponds to the peak in the contribution to the total CSED of the GAMA sample.Both samples are consistent at high Lopt, barring discontinuities due to cosmic variance.Note there is a significant overdensity in G10/COSMOS at z ∼ 0.35 (Darvish et al. 2017).The combined sample is Lopt corrected in the manner described above.
We show the resulting attenuated and unattenuated CSEDs, rescaled using the correction factors shown in Figure 3,in Figures 4 and 5 and Tables 1, 2, 3 and 4. The dotted lines in the Figures show regions where CSEDs are potentially unreliable due to the underlying data either being missing or of too low sensitivity and resolution as described in Driver et al. (2017).

The CSED error budget
Using the prescription in Driver & Robotham (2010)1 , we derive the cosmic (sample) variance (CV) for each of our redshift bins based on areas of 129 deg 2 for GAMA and 0.915 deg 2 for the G10 region (see Table 5).For the combined bins, we compute the fraction of Lopt (55, 37 and 16 per cent for GAMA) arising from each of GAMA and G10/COSMOS; CV is the weighted average (by Lopt) of the respective CVs of GAMA and G10/COSMOS.The Driver et al. ( 2012 .Top panel: the total optical (100 nm < λ < 8 µm) luminosity function for each redshift interval in luminosity bins of log(Lopt) = 0.05 (GAMA: red, G10: blue).Bottom panel: the contribution to the total CSED of each luminosity bin (GAMA: red, G10: blue) fitted by a 10 point spline (green) to the completeness limit (black dashed vertical line).The correction to the CSED normalisation, computed from the ratio of integral of the spline to the histogram with error derived from 1000 Monte Carlo simulations (1 in 10 plotted), is also given.
absolute deviation of 0.007 against a sample of zCOSMOSbright spectroscopic redshifts.However, to be conservative, we add an extra 5 per cent to the G10/COSMOS CSED error budget in quadrature to account for unknown systematics.
We also perform a jackknife resampling by deleting 10 per cent of the sample (as determined by the last digit of the GAMA/G10 catalogue number) for 10 iterations to check whether any given CSED is being dominated by a small number of SED fits.The error is given by σ , where N = number of iterations, mean(f ) is the mean of f across 100 nm < λ < 24 µm (to avoid regions where the CSED is substantially extrapolated), xi is the ith resampled CSED and x the total CSED.This analysis shows our CSED stacks are robust, with errors not exceeding 2 per cent.
The Lopt correction (Figure 3) is also subject to uncertainty.We estimate this error in 1000 Monte Carlo simu- We then refit the spline and rederive the correction factors.The quoted error (Lopt in Table 5) is the standard deviation of correction factors.
The final error budget is shown in Table 5, with the entries combined in quadrature.As expected, CV is the dominant error term.However, all error terms will be reduced by next-generation spectroscopic surveys of galaxy evolution, such as the Wide Area VISTA Extragalactic Survey (WAVES; Driver et al. 2016a) -which will probe significantly larger areas to higher redshifts.

Biases, caveats and missing energy
magphys does not provide an error estimate on the interpolated (best fit) spectrum.We create an indicative estimate of the uncertainty in the CSED shape by recomputing the CSED within each volume through four distinct methods: (i) A strict lower bound to the CSED can be derived by summing the Vmax corrected catalogue's fluxes and errors (in quadrature) for each galaxy in each redshift bin (i.e.ignoring the Lopt correction).Non-detections and non-measurements are treated as zero flux with an error equal to the appropriate 1σ upper limit.This is consistent with the use of 1σ limits in both the underlying photometry and SED fitting.These values are shown in Tables 1 and 2 (first in each set of three measurements) and Figure 6 (pink points).As expected, photometric errors scale with resolution and sensitivity, with the high-resolution optical and near-infrared being associated with negligible error and the low-sensitivity far-infrared data in G10 having the largest error.
(ii) A strict upper bound can be derived from the sum of the Vmax and Lopt corrected fluxes and treating nondetections and non-measurements as having flux equal to the 1σ upper limit.We present this in Tables 1 and 2 (third in each set of three measurements) and as the upper pink error bar in Figure 6.
(iii) We compute the CSED for galaxies with non-zero flux in all photometric bands in each redshift bin and rescale to the full CSED at 1 µm.We depict these CSEDs and respective normalisation constants with blue curves on Figure 6.This sample is biased toward dusty systems, especially those with colder dust populations (due to relatively low sensitivity at 500 µm), and thus provides an alternative worst case scenario (upper limit) for the far-infrared CSED.This curve does not have the same properties in the ultraviolet part of the spectrum because some sources with measurable flux are excluded.
(2) CSED convolved with the filter curve in the observed frame.
(3) Sum of Vmax and Lopt corrected fluxes and upper limits.
restricted CSED represents a lower limit to the ultraviolet CSED.
We then combine these alternative estimates of the CSED to obtain an indicative error range.For the upper bound, we take the minimum of the rescaled no upper limits CSED (blue curve) and a spline interpolation of the sum of fluxes and upper limits (upper pink error bar).The lower bound is a spline interpolation of the photometric measurements minus the photometric error.We also impose a minimum 10 per cent error from the measured CSED (black curve) to account for SED modelling error.The result is the grey regions in Figure 6.The error terms described in Section 3.2 should be added to this range in quadrature.
The large errors present in the far-infrared portion of G10/COSMOS demonstrate the CSEDs in that region are mostly extrapolated and are depicted with dotted lines in Figure 4. Similarly, we depict the GAMA CSEDs with dotted lines beyond 23 µm; the underlying lack of 70 µm imaging prevents the CSED being constrained between 23 and 100 µm.The relatively low detection rate in the Herschel bands reduces the reliability of the far-infrared CSED. Figure 7 shows the fraction of the attenuated, filter convolved (in the observed frame) CSED accounted for by summing Vmax (but not Lopt) corrected fluxes in PACS 160 and SPIRE 250.20 -70 per cent of the attenuated, filter convoved CSED is accounted for by > 3σ detections depending on redshift.The majority of the CSED at 160 and 250 µm below z = 0.45 is accounted for by > 1σ excesses over the background, while above this redshift this fraction is about half.The remainder consists of interpolated and extrapolated flux.
To estimate the effect of incompleteness on the CSED shape, we refer to the ratio of ultraviolet emission between the mass-restricted (green) and unrestricted (black) CSEDs.Assuming that this ratio does not change with redshift, this suggests that the higher redshift GAMA bins are potentially missing significant (∼20 per cent) ultraviolet flux.We do not correct for this effect because it is likely to be redshift dependent.Some of this emission will be attenuated by dust and appear in the far-infrared.This should not significantly affect the overall energy output due to the relatively low amount of energy escaping into the intergalactic medium at these wavelengths.
We also consider the effect of AGN contamination on the CSED.magphys does not incorporate AGN emission into its template library.In G10, we sum fluxes for the objects excluded from the Driver et al. (2017) analysis as AGN using the Donley et al. (2012) criteria.Summing fluxes and errors in quadrature for these objects gives the orange points in Figure 6.When compared to the sum of fluxes for all objects that enter the CSED sample, AGN and their host galaxies consist of ∼4 per cent of the (observed frame) optical and near-infrared flux but ∼ 10 per cent in the midinfrared for z = 0.505 and 12-17 per cent for z = 0.905.As an aside, it is interesting to note that AGN at these redshifts appear to be entirely dominated by their quiescent host galaxies and/or are heavily obscured.For GAMA and (1) Sum of Vmax corrected fluxes.
(2) CSED convolved with the filter curve in the observed frame.
(3) Sum of Vmax and Lopt corrected fluxes and upper limits.the combined bins, we assume AGN contamination to be negligible due to low number density at z < 0.5 (in line with Driver et al. 2017).
At ultraviolet wavelengths, the Driver et al. ( 2012) and the Driver et al. (2016b) CSEDs are based on curve of growth photometry that captures emission from extended UV discs (Gil de Paz et al. 2005;Thilker et al. 2007).The Wright et al. (2016a) catalogue computes consistently derived aperture-matched photometry (defined in the r band) across all wavelengths.This approach trades potentially missing extended flux for consistent errors from the ultraviolet to the infrared, with the latter being more advantageous for SED fitting.
We do not see a significant impact on the near-infrared Table 4. Unattenuated, filter convolved CSED measurements convolved for given filters for G10 in the observed frame, as in Table 3.The unattenuated CSED is negligible in bands longwards of W2.All estimates are subject to the errors described in Table 5 CSED at low redshifts caused by uncertainty in modelling thermally-pulsating asymptotic giant branch (TP-AGB) stars (see e.g.Maraston et al. 2006;Bruzual 2007;Conroy & Gunn 2010;Bruzual et al. 2013;Noël et al. 2013;Capozzi et al. 2016) -our CSEDs lie within the photometric bounds for GAMA.It is unknown whether magphys compensates for imprecision in TP-AGB modelling by adjusting the other SED fitting parameters (such as stellar phase metallicity or ages) in order to achieve a mathematically good, but unphysical, fit.Using a SED fitting code based on the Maraston (2005) (which arguably overestimate the TP-AGB contribution) or other stellar library to determine the uncertainty range would naturally be part of a thorough evaluation of SED fitting, which is also out of scope of this paper.
The aperture photometry used within may systematically miss flux in highly concentrated systems, which are typically massive, quiescent and dust-free galaxies.One could determine an aperture correction by fitting a Sérsic profile to each source and computing the ratio of the inte-grated flux to the aperture flux (see e.g.Taylor et al. 2011).These come with their own set of biases and random errors, which may be wavelength dependent.No such analysis has been performed on G10/COSMOS, so for consistency across the redshift range we do not include it.
The Lopt correction only serves to adjust the normalisation of the CSED.This creates a bias against fainter systems if they have substantially different SEDs than the average population and, in aggregate, contribute significantly to the total CSED.This may be evident in e.g. the peak of the far-infrared CSED at 0.14 < z < 0.2 compared to 0.2 < z < 0.28, with the caveat that a greater portion of the far-infrared CSED arises from extrapolation in the higher redshift bin and the complete lack of 70 µm data in GAMA compared to at least an upper limit in G10/COSMOS.This is best addressed with deeper data.
The r and i + selections used to define the GAMA and G10/COSMOS samples are equivalent to a luminosity cut at successively shorter wavelengths in the rest frame with redshift.The most obvious resulting bias is against highly obscured or quiescent systems which lie below the respective magnitude limits.More subtle biases may occur depending on the SEDs of individual sub-populations.The largest effect occurs where r or i + is equivalent to u (and to a lesser extent, g) in the rest frame, namely the high-redshift ends of the GAMA and G10/COSMOS samples.This is not an issue in GAMA out to z = 0.10, as our low-redshift CSED is consistent with that constructed from luminosity functions in each band from FUV through K by Driver et al. (2012).The effect of these biases can only be quantified with a deep, near-infrared selected sample, which will also enable studies of the CSED at higher redshifts.

The energy output of the Universe
Integrating under our CSEDs gives the expected trend of declining energy output (Figure 8), consistent with star formation activity winding down since cosmic noon at z = 1 ) Scaling: 24.01 Scaling: 1.2 q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q 10 −7 10 −6 10 −5 10 −4 10 −3 ) Scaling: 39.36 Scaling: 1.22 q q qq q q q q qq q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q Driver+17 AGN 10 −7 10 −6 10 −5 10 −4 10 −3 ) Scaling: 48.17 Scaling: 1.21 q q qq q q q q qq q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q 10 −7 10 −6 10 −5 10 −4 10 −3 ) Scaling: 75.39 Scaling: 1.14 q q qq q q q q qq q q q q q q q q q q q q q q q q qq q q q q qq q q q q q q q q q q q q q q Table 6.Energy output as a function of redshift computed by integrating under the (attenuated) CSED over the full wavelength range and between 100 nm and 8 µm (optical only).The quoted errors are derived from Table 5.
One exception is the decrease in energy output at optical wavelengths between the 0.14 < z < 0.20 and 0.45 < z < q q q q q q q q q q 0 0.05 0.1 0.15 0.2 0.25 0.3 1 2 5 10 q q q q 0 0.1 0. 0.56 bins.This we measure as the total energy output over the rest frame wavelength range 100 nm < λ < 8 µm with increasing redshift.The decline at intermediate redshifts is not significant with respect to CV. Ultraviolet and far-infrared emission, where significant extrapolation is present, continue to increase at this epoch, as shown in Figure 4. Overall, the approximately factor four decline -from 5.1 ± 1.0 × 10 35 h70 W Mpc −3 at z = 0.91 to 1.3 ± 0.3 × 10 35 h70 W Mpc −3 at z = 0.06 is robust to all forms of error.The decline in the total integrated energy output (E = λ dλ) is well fit by a smooth power law: log(E) = (2.18 ± 0.12) log(1 + z) + (0.07 ± 0.02) + 35 (2) This decline occurs at a slightly slower rate than the cosmic star formation density (see Figure 8).This is expected, as the energy density reflects not only current star formation, but also remaining stellar populations from previous star formation.The size of this effect is small -resulting in a difference of a factor of ∼1.25 over the redshift range -indicating that the energy output of the Universe is still dominated by young stars.Note that the cosmic star formation density has been arbitrarily scaled to match the above CSED fit in the highest redshift bin.Beyond z = 1 (not shown on Figures 4, 5 and 8), the CSED unphysically declines rapidly due to incompleteness in the Andrews et al. (2017) catalogue.A deep, near-infrared selected catalogue is required to probe to higher redshifts.

The integrated photon escape fraction
Dividing the attenuated CSED by the unattenuated CSED gives the integrated photon escape fraction (IPEF), which we show in Figure 11.Also shown in Figure 11 are the Driver et al. (2012) values adopted from the Millenium Galaxy Catalogue data spanning 0 < z < 0.18 (Driver et al. 2008).The IPEF is a simple measure of the effect dust attenuation has on galaxy light.The IPEF is not subject to cosmic variance, jackknife, Lopt correction or photometric redshift errors as the uncertain normalization of the CSED at a given redshift is divided out.This leaves errors in the CSED shape only and dust attenuation as the dominant sources of uncertainty.These uncertainties are hard to estimate without a thorough evaluation of the SED modelling process, so we assign an indicative 20 per cent error to account for this.
The IPEFs shown here are consistent with increasing opacity with lookback time, which presumably is linked to an increasing dust mass density.This is consistent with the findings of Dunne et al. (2011) and Driver et al. (2017) and not surprising given the correlation between star formation and dust opacity and the cosmic star formation density rising with redshift.This is also consistent with the evolution in the ultraviolet and infrared luminosity functions from z > 1 to z = 0 as noted by Bernhard et al. (2014).The escape fraction at 150 nm decreases from 24 ± 5 per cent at z = 0.05 to 16 ± 3 per cent at z = 0.915.The escape fraction at 250 nm also declines from 40 ± 8 per cent to 26 ± 5 per cent over the same time intervals.This is in line with Driver et al. (2016b), with the caveat that the dust attenuation and re-emission curves may be unreliable as explained in Section 3.3.These attenuation curves are similar to that of the Milky Way (Cardelli et al. 1989) when normalised to an escape fraction of 0.6 at 549.5 nm (black dashed line), with the lower redshift curves being more transparent.It is clear that further reduction in SED modelling errors is required to obtain a Universal extinction curve.
Convolving the IPEF with a filter curve and converting to magnitudes gives an attenuation coefficient.We show the resulting FUV attenuation coefficients (AFUV) with the 20 per cent indicative uncertainty in the IPEF in Figure 12.Our estimates are in line with those derived from Cucciati et al. (2012) and Burgarella et al. (2013) using VVDS, PEP and HerMES data across the redshift range.

The integrated galactic light
The integrated galactic light (IGL) at z = 0 can be derived from the CSED as follows: where d l (z) is the luminosity distance and dV (z) is the differential comoving volume of each redshift bin subtending a solid angle of 1 sr. Figure 13  0.8 1 q q q q q q q q q q q q Cardelli+89 Driver+12 0.02 < z < 0.08 0.08 < z < 0.14 0.14 < z < 0.2 0.20 < z < 0.28 0.28 < z < 0.36 0.36 < z < 0.45 0.45 < z < 0.56 0.56 < z < 0.68 0.68 < z < 0.82 0.82 < z < 0.99 Rest frame wavelength (m) Integrated photon escape fraction are also shown as is the average Milky Way attenuation curve (Cardelli et al. 1989) assuming an escape fraction of 0.6 at 549.5 nm.Curves are subject to CSED shape errors as described in Section 3.3.0 0.2 0.4 0.6 0 0.5 1 1.5 2 2.5 q q q q q q q q q q q q q q q q q q q 0 0.25 0.5 1 1.5 2 2.5 3 4 This work Cucciati+12 Burgarella+13 Figure 12.A FUV as a function of redshift (computed from convolving the IPEF curves with the GALEX FUV filter, with a 20 per cent indicative error in the IPEF) compared to Cucciati et al. (2012) and Burgarella et al. (2013).λ < 1 mm) and IGL as a function of redshift shown in Table 7. Figure 14 shows the fraction of the cosmic optical and infrared backgrounds accounted for by emission prior to the redshift.As expected, low redshifts dominate FUV emission and high redshifts dominate in the mid-infrared and longwards of the cold dust peak.We do not recover the total IGL at any wavelength, indicating a significant portion was emitted before z = 1.Areas with relatively low level of IGL recovery should be easily explained by the high-redshift galaxy population.
Contributions to the cosmic optical background (see Table 7) are approximately constant as a function of redshift, while contributions to the cosmic infrared background increase.For comparison, Driver et al. (2016c) measured cosmic optical and infrared backgrounds of 24 ± 4 and 26 ± 5 nW m −2 sr −1 respectively, giving the shaded areas in Figure 14.Our measurements are in excellent agreement with contributions to the cosmic infrared background measured by Béthermin et al. (2012a) (green points/lines).
The combination of WAVES and a high-redshift sample based on James Webb Space Telescope (JWST ), Wide Field Infrared Survey Telescope (WFIRST ) and Euclid data will allow the determination of the optical CSED out to the epoch of reionisation and thus the full characterisation of the IGL.

CONCLUSION
In this work, we measured the CSED for 0 < z < 1 by stacking SED fits from Driver et al. (2017).We binned the SEDs into 10 different redshift intervals and measured the CSED and integrated photon escape fraction for each through stacking.We found that energy output declines as a function of redshift, from (5.1 ± 1.0) × 10 35 h70 W Mpc −3 to (1.3 ± 0.3) × 10 35 h70 W Mpc −3 (Figure 8).This decline is robust despite cosmic variance (CV) and other uncertainties and occurs at a rate slightly slower than the decline in cosmic star formation.Combined with the reddening of the unattenuated CSED, this is consistent with the mean age of stellar populations becoming older.AGN do not contribute significantly to the CSED at any z < 1.We also show 1 10 100 0.02 < z < 0.08 0.08 < z < 0.14 0.14 < z < 0.20 0.20 < z < 0.28 0.28 < z < 0.36 0.36 < z < 0.45 0.45 < z < 0.56 0.56 < z < 0.68 0.68 < z < 0.82 0.82 < z < 0.99 Total q q q q q q q q q q COB log (1 + z) q q q q q q q q q q q This work Driver+16 Bethermin+12

CIB
log (1 + z) that the photon escape fraction also declines with increasing redshift, equivalent to an increase in AFUV of 0.4 mag, consistent with an increase in the cosmic dust density.The CSEDs presented here complement the recent measurements of the integrated galactic light reported by Driver et al. (2016c).We will follow up the CSEDs presented in this work with a comparison to semi-analytic (e.g.Gilmore et al. 2012) and empirical models (e.g.Driver et al. 2013) of the CSED and integrated galactic light as a function of cosmic time (Andrews et al. in prep).We will also explore the physical properties obtained in our SED fits for various (sub-)populations of galaxies (e.g.Kelvin et al. 2014) and investigate the biases inherent in the CSEDs in more detail.Finally, we will compare our physical parameters with measurements made using independent methods derived from the GAMA database and use them to improve the SED fitting process.
The dominant source of directly quantifiable error in our CSED estimates is CV.Incompleteness remains a significant impediment towards fully characterising the CSED over cosmic time.While the optical (100 nm < λ < 8 µm) portion of the CSED is well constrained, the relatively low sensitivity and resolution of the far-infrared data prevents the full characterization of the CSED at all redshifts.Incompleteness also causes a systematic underestimation of the ultraviolet and far-infrared flux of 10-30 per cent, as low luminosity galaxies are preferentially star forming and dusty.Systematic errors also arise from SED modelling and the non-modelling of AGN.
Data from next-generation surveys of galaxy evolution, such as WAVES and the Maunakea Spectroscopic Explorer, will significantly reduce incompleteness and CV.These surveys aim to cover hundreds of square degrees of sky to 2-4 magnitudes fainter than GAMA with a high level of spectroscopic completeness.In addition, this will enable us to replace photometric redshifts with spectroscopic data for a significant portion of the sample.The large survey footprint will also reduce the uncertainty in the CSED normalization due to cosmic variance -for example the 26 per cent CV in the 0.915 deg 2 observable portion of the G10 field for 0.45 < z < 0.56 will be reduced to 7 per cent for an illustrative survey of two 25 deg 2 fields.The combination of WAVES with observations equivalent to COSMOS using JWST and WFIRST will, once combined using the procedure in this work, constrain the ultraviolet to near-infrared CSED out to the epoch of reionisation and thus how the cosmic optical background builds up over the history of the Universe.A. et al. 2007, ApJ, 172, 468 Zemcov M. et al. 2014, Science, 346, 732

Figure 2 .
Figure 2. Stellar mass versus redshift for the GAMA (black) and G10/COSMOS (grey) samples (1 in 10 plotted).Complete regions are denoted by green (as determined from the peak of the distribution of stellar masses in each redshift bin), incomplete regions are denoted by yellow and unobservable regions denoted by red.Redshift bins where a combined sample was used are denoted by blue.
Figure3.Top panel: the total optical (100 nm < λ < 8 µm) luminosity function for each redshift interval in luminosity bins of log(Lopt) = 0.05 (GAMA: red, G10: blue).Bottom panel: the contribution to the total CSED of each luminosity bin (GAMA: red, G10: blue) fitted by a 10 point spline (green) to the completeness limit (black dashed vertical line).The correction to the CSED normalisation, computed from the ratio of integral of the spline to the histogram with error derived from 1000 Monte Carlo simulations (1 in 10 plotted), is also given.

FigureFigure 7 .
Figure 6 -continued shows the fraction of theDriver et al. (2016c) IGL accounted for by redshifting and summing the CSEDs measured in this work.Our CSEDs constitute a roughly constant 40-80 per cent of the IGL across the entire wavelength range, with contributions to the cosmic optical background (the portion of the IGL between 100 nm and 8 µm), cosmic infrared background (8 µm <

Figure 11 .
Figure 11.The integrated photon escape fraction for the GAMA and G10 datasets.The IPEFs estimated by Driver et al. (2012) (grey line and black points) and Driver et al. (2016b) (dashed)are also shown as is the average Milky Way attenuation curve(Cardelli et al. 1989) assuming an escape fraction of 0.6 at 549.5 nm.Curves are subject to CSED shape errors as described in Section 3.3.

Figure 13 .
Figure 13.The fraction of the total IGL, measured by the fitting function of Driver et al. (2016c) (top: as a function of wavelength and redshift, bottom: per photometric filter from GALEX, SDSS, VISTA, IRAC, MIPS, PACS and SPIRE and summed over all redshifts), contributed by our CSEDs.Measurements are subject to the errors described in Table5and the CSED shape errors described in Section 3.3.

Figure 14 .
Figure 14.The portion of Driver et al. (2016c) cosmic optical (left) and infrared (right) backgrounds recovered as a function of redshift compared against Béthermin et al. (2012a).The shaded areas represent the uncertainty range of the Driver et al. (2016c) cosmic optical and infrared background measurements.

Table 1 .
The attenuated CSEDs for the GAMA and G10 datasets.Dashed lines indicate regions where the respective CSEDs are poorly constrained or partially extrapolated due to lack of data.The curves are subject to the normalisation errors described in Table5, with uncertainty in the shape discussed in Section 3.3.The unattenuated CSEDs for the GAMA and G10 datasets.The curves are subject to the normalisation errors described in Table5, with uncertainty in the shape discussed in Section 3.3.Attenuated, filter convolved CSED measurements for GAMA.At each epoch, three measurements are given -(1) the sum of the Vmax corrected fluxes, (2) the CSED convolved with the filter curve in the observed frame and (3) the sum of the Vmax and Lopt corrected fluxes and upper limits.These measurements are subject to errors described in Table5, photometric error is negligible.

Table 2 .
Attenuated, filter convolved CSED measurements for the intermediate redshift G10 region in the observed frame, as in Table1.All measurements are subject to the errors described in Table5.

Table 3 .
Unattenuated, filter convolved CSED measurements for the given filters for GAMA in the observed frame.The unattenuated CSED is negligible in bands longwards of W2.All estimates are subject to the errors described in Table5.

Table 5 .
. The CSED error terms arising from cosmic variance (CV), jackknife resampling, use of photometric redshifts and the Lopt correction.The total error in per cent for a given redshift bin is the sum of all columns in quadrature for that bin.
Madau & Dickinson (2014) Universe as measured from this work,Driver et al. (2012)andDriver et al. (2016b).Errors depicted are those described in Table5only.For comparison, we add the Driver et al. (2017) cosmic star formation rate density measurements and theMadau & Dickinson (2014)cosmic star formation rate density fitting function, arbitrarily scaled to the fit to the total energy output for the highest redshift bin.

Table 7 .
Contributions to the cosmic optical and infrared backgrounds (COB, CIB) and IGL as a function of redshift.Quoted uncertainties are derived from Table5.

Table 8 .
Extract from the lowest redshift CSED.The full tables for all redshifts are available as supporting information.