The cosmic spectral energy distribution in the EAGLE simulation

TL;DR

This study compares the cosmic spectral energy distribution from observations and the EAGLE simulation, finding good agreement at low redshift but discrepancies at higher redshift due to simulation limitations and dust evolution uncertainties.

Contribution

It demonstrates that combining EAGLE simulation with SKIRT dust modeling produces a realistic local Universe CSED, highlighting areas for improvement at higher redshifts.

Findings

Excellent match with observed CSED at z=0 across UV to submm wavelengths.

Underestimation of CSED at higher redshift, especially in FIR/submm.

Discrepancies attributed to simulation volume and dust evolution uncertainties.

Abstract

The cosmic spectral energy distribution (CSED) is the total emissivity as a function of wavelength of galaxies in a given cosmic volume. We compare the observed CSED from the UV to the submm to that computed from the EAGLE cosmological hydrodynamical simulation, post-processed with stellar population synthesis models and including dust radiative transfer using the SKIRT code. The agreement with the data is better than 0.15 dex over the entire wavelength range at redshift $z=0$, except at UV wavelengths where the EAGLE model overestimates the observed CSED by up to a factor 2. Global properties of the CSED as inferred from CIGALE fits, such as the stellar mass density, mean star formation density, and mean dust-to-stellar-mass ratio, agree to within better than 20 per cent. At higher redshift, EAGLE increasingly underestimates the CSED at optical-NIR wavelengths with the FIR/submm emissivity underestimated by more than a factor of 5 by redshift $z=1$. We believe that these differences are due to a combination of incompleteness of the EAGLE-SKIRT database, the small simulation volume and the consequent lack of luminous galaxies, and our lack of knowledge on the evolution of the characteristics of the interstellar dust in galaxies. The impressive agreement between the simulated and observed CSED at lower $z$ confirms that the combination of EAGLE and SKIRT dust processing yields a fairly realistic representation of the local Universe.

Figures

Tables

Table 1. Table 1: Main characteristics of the different EAGLE runs used in this paper. Columns from left to right: EAGLE model name; comoving volume size; dark matter particle mass; initial baryonic particle mass; the total number of galaxies with M ⋆ > 10 8.5 ​ M ⊙ subscript 𝑀 ⋆ superscript 10 8.5 subscript 𝑀 direct-product M_{\star}>10^{8.5}~{}M_{\odot} for all 29 snapshots combined; the fraction of this total number of galaxies with sufficient dust to be included in the EAGLE-SKIRT catalogue; notes about the subgrid physics recipes. The recalibrated high-resolution simulation, RecalL0025N0752, is the run on which the main analysis in this paper is based.

EAGLE run	$L$	$N_{\text{tot}}$	$m_{\text{DM}}$	$m_{\text{gas}}$	$N_{\text{gal}}$	$f_{\text{dusty}}$	subgrid calibration
	[Mpc]		[$M_{\odot }$]	[$M_{\odot }$]		[%]
RefL0100N1504	$100$	$2\times {1504}^3$	$9.70\times {10}^6$	$1.81\times {10}^6$	$371,728$	$63.6$	fiducial calibration
RefL0050N0752	$50$	$2\times {752}^3$	$9.70\times {10}^6$	$1.81\times {10}^6$	$48,261$	$65.1$	fiducial calibration
AGNdT9L0050N0752	$50$	$2\times {752}^3$	$9.70\times {10}^6$	$1.81\times {10}^6$	$48.278$	$64.7$	different AGN feedback
RefL0025N0376	$25$	$2\times {376}^3$	$9.70\times {10}^6$	$1.81\times {10}^6$	$5742$	$67.4$	fiducial calibration
RefL0025N0752	$25$	$2\times {752}^3$	$1.21\times {10}^6$	$2.26\times {10}^5$	$8279$	$94.4$	fiducial calibration
RecalL0025N0752	$25$	$2\times {752}^3$	$1.21\times {10}^6$	$2.26\times {10}^5$	$5954$	$95.7$	recalibrated subgrid parameters

Table 2. Table 2: Global characteristics of the local Universe as derived from SED model fits to the CSED. The third column gives the intrinsic SFR density and stellar mass contributing obtained by directly summing the SFRs and stellar masses of all galaxies in the EAGLE-SKIRT volume. The fourth column corresponds to the CIGALE fit to the EAGLE-SKIRT CSED as derived in this paper. The fifth column corresponds to our CIGALE fit to the observed GAMA CSED data points from Andrews et al. ( 2017b ) . The last column corresponds to the detailed study by Driver et al. ( 2018 ) , based on MAGPHYS SED fits for individual galaxies from the GAMA, G10-COSMOS and 3D-HST surveys.

quantity	unit	EAGLE-SKIRT	EAGLE-SKIRT	GAMA	GAMA
		intrinsic	CIGALE	CIGALE	MAGPHYS
$\log {\epsilon }_{\text{bol}}$	W Mpc-3	$\cdots$	$35.00$	$35.03$	$35.10$
${\rho }_{\text{SFR}}$	$M_{\odot }$ yr-1 Mpc-3	$0.011$	$0.017$	$0.014$	$0.011$
${\Omega }_{\star }$	$-$	$0.0012$	$0.0015$	$0.0018$	$0.0015$
${\Omega }_{\text{d}}$	$-$	$\cdots$	$6.5\times {10}^{-7}$	$7.8\times {10}^{-7}$	$1.1\times {10}^{-6}$
${\Omega }_{\text{d}}/{\Omega }_{\star }$	$-$	$\cdots$	$4.4\times {10}^{-4}$	$4.3\times {10}^{-4}$	$7.4\times {10}^{-4}$
$f_{\text{abs}}$	$-$	$\cdots$	$0.27$	$0.30$	$\cdots$
$A_{\text{FUV}}$	mag	$\cdots$	$0.91$	$1.49$	$1.52$

Equations

$$\epsilon \equiv \nu {ϵ}_{\nu }\equiv \lambda {ϵ}_{\lambda }$$

$${\epsilon }_{\text{bol}}=\int {\epsilon }_{\nu }\text{d}\nu ,$$

Full text

††pubyear: 2019††pagerange: The cosmic spectral energy distribution in the EAGLE simulation–References

1 Introduction

The cosmic spectral energy distribution (CSED) is a fundamental observational characteristics of the Universe. It represents the total electromagnetic power generated within a cosmological unit volume as a function of wavelength. In the wavelength range between the UV and the submm, the vast majority of the emission is emitted by stars, gas and dust within galaxies, and hence the CSED is a complex function of both the volume density of different galaxy types and the different processes that shape the SEDs of individual galaxies.

Until a few years ago, most efforts to measure the CSED were concentrated on limited regions of the electromagnetic spectrum: UV (Wilson et al., 2002; Budavári et al., 2005; Wyder et al., 2005; Robotham & Driver, 2011), optical (Norberg et al., 2002; Blanton et al., 2003; Montero-Dorta & Prada, 2009), near-infared (Cole et al., 2001; Kochanek et al., 2001; Smith et al., 2009), mid-infrared (Saunders et al., 1990; Babbedge et al., 2006) and far-infrared/submm (Takeuchi et al., 2006; Bourne et al., 2012; Marchetti et al., 2016). The full CSED can then be obtained by combining these different measurements. However, this approach often results in mutually inconsistent results, due to different source selection criteria, different levels of completeness, photometric measurement discrepancies, or cosmic variance (Cross & Driver, 2002; Driver & Robotham, 2010; Driver et al., 2012).

The optimal approach to measure the CSED is to use a single and spectroscopically complete volume-limited sample with multi-wavelength data covering the entire UV-submm range. Such surveys have lately become available, with Galaxy And Mass Assembly (GAMA: Driver et al., 2011; Liske et al., 2015) probably being the most complete one in the local Universe. Driver et al. (2012) measured the local CSED from the UV to the NIR directly from GAMA data and used SED templates to extrapolate it to the submm range. Kelvin et al. (2014) considered the contribution of different galaxy types to the UV–NIR CSED and showed that all types contribute significantly to the ambient inter-galactic radiation field. Driver et al. (2016) used the GAMA Panchromatic Data Release to derive the first self-consistent measurement of the entire UV-submm CSED. Binning their galaxy sample in three redshift bins, they find a clear signature for evolution in the CSED over the past 2.3 Gyr. Very recently, Andrews et al. (2017b) took this study a significant step further. They included newly reduced panchromatic data from the Cosmological Origins Survey (COSMOS: Scoville et al., 2007; Davies et al., 2015; Andrews et al., 2017a), which enabled them to study the evolution of the CSED out to $z=1$. They demonstrated that the bolometric UV-submm energy output of the Universe has decreased by about a factor four over the past 8 Gyr.

Obviously, the CSED is nothing but the joint contribution of the individual spectral energy distributions (SEDs) of all galaxies within a cosmological unit volume. Overall, the shape and normalisation of the CSED, and its evolution with redshift, are strong and purely observable constraints for models of galaxy formation and evolution (e.g., Domínguez et al., 2011).

Galaxy formation and evolution models come in two broad classes: semi-analytical models and cosmological hydrodynamical simulations. The former have been around for more than two decades (e.g., Kauffmann et al., 1993; Cole et al., 1994; Somerville & Primack, 1999). Very recently, different versions of the GALFORM semi-analytical model (Cole et al., 2000; Gonzalez-Perez et al., 2014; Lacey et al., 2016) were used to predict the CSED, generally showing good agreement between the models and the data at low redshifts (Andrews et al., 2018; Cowley et al., 2019).

Cosmological hydrodynamical simulations, on the other hand, have only fairly recently achieved sufficient realism and statistics to become a serious contender for the semi-analytical models. Modern cosmological hydrodynamical simulations such as Illustris (Vogelsberger et al., 2014), EAGLE (Schaye et al., 2015), MassiveBlack-II (Khandai et al., 2015), MUFASA (Davé et al., 2016), and IllustrisTNG (Pillepich et al., 2018) can reproduce many characteristics of the present-day galaxy population, including the stellar mass function, the size distribution, the bimodality in colours, and the supermassive black hole mass function. To the best of our knowledge, the CSED of cosmological hydrodynamical simulations has never been calculated in a self-consistent way over the entire UV-submm wavelength range.

In order to do so, one needs to calculate the panchromatic SED of each galaxy in the simulation. This does not only include the stellar emission by different stellar populations, but also the distorting effect of interstellar dust. Indeed, dust attenuates roughly 30% of all the starlight in normal star-forming galaxies, and re-emits it as thermal emission in the infrared (Buat & Xu, 1996; Popescu & Tuffs, 2002; Viaene et al., 2016). This means that detailed 3D dust radiative transfer calculations are required. In the past few years, 3D dust radiative transfer has seen a remarkable development (Steinacker et al., 2013), and such detailed self-consistent calculations in a realistic, galaxy-wide setting have become possible (Jonsson et al., 2010; Hayward et al., 2011; De Looze et al., 2014; Domínguez-Tenreiro et al., 2014; Natale et al., 2015; Saftly et al., 2015; Guidi et al., 2015; Goz et al., 2017). In particular, panchromatic dust radiative transfer post-processing can nowadays be applied to large numbers of simulated galaxies from cosmological hydrodynamical simulations (Camps et al., 2016; Camps et al., 2018; Rodriguez-Gomez et al., 2019).

The goal of this paper is to compare the CSED of the EAGLE suite of simulations to the observed CSED as measured from GAMA observations. In Section 2 we briefly describe the EAGLE simulations and the mock observations database we use for our study. In Section 3 we compare the EAGLE-SKIRT CSED of the local Universe to observational data from the GAMA survey, and we derive and discuss a number of important characteristics of the local Universe. In Section 4 we compare the local Universe CSED corresponding to different EAGLE simulations, and discuss strong and weak convergence, and a variation of the subgrid model parameters. In Section 5 we discuss the cosmic evolution of the EAGLE-SKIRT CSED out to $z\sim 1$ and compare it again to observational data from GAMA and G10/COSMOS. In Section 6 we summarise our results.

Throughout this paper, we adopt $H_{0}=67.77$ km s*-1* Mpc*-1*, the Planck cosmology value adopted by EAGLE (Planck Collaboration XVI, 2014).

2 The simulations

2.1 The EAGLE simulations

EAGLE (Evolution and Assembly of GaLaxies and their Environments: Schaye et al., 2015) is a suite of cosmological hydrodynamics simulations. The simulations were run using an updated version of the $N$-body/SPH code GADGET-3 (Springel, 2005). Different runs correspond to different volumes (cubic volumes ranging from 25 to 100 comoving megaparsecs on a side), different resolutions, and different physical prescriptions for baryonic processes including star formation, AGN feedback and cooling. The main characteristics of the most important EAGLE runs are listed in Table 1.

The EAGLE simulations track the cosmic evolution of dark matter, baryonic gas, stars and massive black holes. As most other cosmological simulations, EAGLE lacks the resolution and the detailed physical recipes to model the cold phase in the ISM. To prevent artificial fragmentation of star-forming gas, the EAGLE simulations impose a pressure floor, corresponding to a polytropic equation of state (Schaye & Dalla Vecchia, 2008). As a consequence, the ISM does not consist of resolved molecular clouds, but rather of a fairly smoothly distributed, pressurised gas. The most important criterion for star formation in this simulated ISM is a metallicity-dependent density threshold (Schaye, 2004). Star formation is implemented according to the observed Kennicutt-Schmidt law (Schmidt, 1959; Kennicutt, 1998), and with a Chabrier (2003) initial mass function. The stellar evolution and chemical enrichment prescriptions in EAGLE are taken from Wiersma et al. (2009).

The EAGLE simulations have been calibrated to reproduce the local Universe stellar mass function, the galaxy-central black hole mass relation, and the galaxy mass-size relation, as described by Crain et al. (2015). The simulations were shown to be in reasonable to excellent agreement with many other observables not considered in the calibration, including the H2 galaxy mass function (Lagos et al., 2015), the relation between stellar mass and angular momentum (Lagos et al., 2017), the supermassive black hole mass function (Rosas-Guevara et al., 2016), the atomic gas properties of galaxies (Crain et al., 2017), the galaxy size evolution (Furlong et al., 2017), and the evolution of the star formation rate function (Katsianis et al., 2017).

2.2 The EAGLE-SKIRT database

EAGLE doesn’t include dust as a separate species. Camps et al. (2016) and Trayford et al. (2017) introduced an advanced framework to add interstellar dust to the EAGLE galaxies, and to calculate mock observables that fully take into account the absorption, scattering and thermal emission by this dust. The recipe includes a resampling procedure for star forming particles, the use of MAPPINGS SED templates (Groves et al., 2008) to model dusty Hii regions, and the inclusion of a diffuse dust distribution based on the distribution of metals in the gas phase. The final step in the procedure is a full 3D dust radiative transfer simulation using the SKIRT code. SKIRT (Camps & Baes, 2015) is an open-source 3D Monte Carlo dust radiative transfer code, equipped with advanced grids for spatial discretisation (Saftly et al., 2013, 2014), a hybrid parallelisation scheme (Verstocken et al., 2017), a library of input models (Baes & Camps, 2015), and a suite of optimisation techniques (Baes et al., 2003, 2011; Baes et al., 2016; Steinacker et al., 2013).

Camps et al. (2016) used this framework to calculate mock infrared and submm observations for a limited set of EAGLE $z=0$ galaxies, selected to match a sample of nearby galaxies selected from the Herschel Reference Survey (HRS, Boselli et al., 2010). The parameters in the post-processing scheme were calibrated to reproduce the observed HRS submm colours and dust scaling relations (Boselli et al., 2012; Cortese et al., 2012, 2014). Subsequently, Trayford et al. (2017) used this calibrated recipe to calculate mock optical images, broadband fluxes, colours, and spectral indices for more than 30,000 local Universe EAGLE galaxies. One of their conclusions is that this radiative transfer recipe shows a marked improvement in the color versus stellar mass diagram over a simple dust-screen model.

Camps et al. (2018) have used a refined version of this postprocessing recipe to populate a database of mock observations for the EAGLE simulations. The resulting EAGLE-SKIRT database contains mock UV to submm flux densities and rest-frame luminosities for nearly half a million simulated galaxies, distributed over 23 redshift slices from $z=0$ to $z=6$. The mock observations were calculated for six different EAGLE runs, corresponding to different box sizes, mass resolutions and physical ingredients (see Table 1). All these data are available in the public EAGLE database (McAlpine et al., 2016).

An important caveat is that only galaxies with at least 250 dust particles111The number of dust particles in a simulated EAGLE galaxy is defined as $N_{\text{dust}}=\max(N_{\text{coldgas}},N_{\text{SFR}})$, where $N_{\text{SFR}}$ and $N_{\text{coldgas}}$ indicate the number of (sub)particles in the sets representing the star-forming regions and cold gas particles, respectively. These sets may contain original SPH particles extracted from the EAGLE snapshot and/or resampled sub-particles replacing star-forming region candidates. See Camps et al. (2018, §3.1) for details. are considered in the EAGLE-SKIRT catalogue. As discussed in detail by Camps et al. (2018), a minimum number of dust particles is required for the radiative transfer postprocessing to be meaningful, and 250 was found to be an appropriate number. This threshold leads to an incompleteness of the EAGLE-SKIRT catalogue, with a bias against red and dead early-type galaxies (with a high stellar mass, but little dust) and against low-mass galaxies (with low stellar masses and low dust masses). The level of incompleteness differs for the different EAGLE runs: in the high-resolution RecalL0025N0752 run, 95.7% of all galaxies with stellar masses above $10^{8.5}~{}M_{\odot}$ are included in the database, whereas this fraction drops to 63.6% for the largest-volume RefL0100N1504 run (Camps et al., 2018, Table 1).

3 The local Universe CSED

We calculated the CSED at $\langle z\rangle=0.05$ as derived from the EAGLE-SKIRT database. We base our main analysis on the recalibrated 25 Mpc volume simulation, RecalL0025N0752, as this simulation has the highest resolution. An analysis of the effect of resolution, simulation volume and subgrid physics recipes will be presented in Section 4. Since the two last EAGLE snapshots correspond to $z=0$ and $z=0.1$, we averaged over these two snapshots. At each snapshot and in each broadband filter, we calculated the CSED by summing the observed flux densities of every single galaxy in the EAGLE-SKIRT database, and subsequently normalising the sum based on the snapshot co-moving volume and luminosity distance. We calculated the CSED in 31 broadband filters covering the UV-submm wavelength regime.

In Figure 1 we compare the EAGLE-SKIRT CSED at $\langle z\rangle=0.05$ to the observed GAMA CSED from Andrews et al. (2017b) for the redshift bin $0.02<z<0.08$, the lowest redshift bin they consider. The CSED is plotted as

[TABLE]

and has units of total power per unit volume (W Mpc*-3*). While the overall agreement between the GAMA observations and the EAGLE-SKIRT data points is very satisfactory, there are some minor differences to be noted, as we discuss next.

3.1 The UV-optical-NIR CSED

First of all, the EAGLE-SKIRT results underestimate the GAMA CSED at red optical and near-infrared wavelengths ($\sim$0.6–4 $\mu$m). The difference is small (about 0.13 dex or about 30%), but systematic. A first possible contributor to this systematic difference is the incompleteness of the EAGLE-SKIRT catalogue, in particular regarding the spheroidal galaxy population. Kelvin et al. (2014) find that more than half of the observed optical/NIR energy budget is dominated by galaxies with a significant spheroidal component. Herschel observations have shown that many early-types have dust masses below $5\times 10^{5}~{}M_{\odot}$ (Smith et al., 2012b; di Serego Alighieri et al., 2013). Such galaxies, with a substantial contribution in the optical but almost no dust, could easily drop out of the EAGLE-SKIRT catalogue.

Secondly, our simple method to estimate the CSED might also contribute to the systematic underestimation. We have estimated the CSED by simply adding the flux of each galaxy in the snapshot volume. As such, we might miss a non-negligible contribution from galaxies at the high-luminosity side of the luminosity function that is not properly sampled by the EAGLE RecalL0025N0752 simulation. Trayford et al. (2015) have presented optical and NIR luminosity functions for the EAGLE simulations, based on mock fluxes that take dust absorption into account with a simple heuristic recipe. Their Figure 3 clearly demonstrates that the RecalL0025N0752 simulation misses the more luminous sources in the optical/NIR bands due to the relatively small simulation volume sampled. More specifically, the RecalL0025N0752 luminosity functions drop to zero around or even before the knee of the luminosity function. This insensitivity to the exponential cut-off at high luminosities, due to the small simulation volume, definitely contributes to the underestimation of the EAGLE-SKIRT CSED in the optical/NIR bands.

Finally, part of this discrepancy is due to the EAGLE simulation itself: Trayford et al. (2015) note that their luminosity functions are consistent with a slight underestimate in the masses of more massive EAGLE galaxies, as already seen in the mass function shown by Schaye et al. (2015).

In the SDSS g and u bands, the discrepancy between the EAGLE-SKIRT and the observed CSED is smaller than in the red and near-infrared filters, and at UV wavelengths, the EAGLE-SKIRT results even overestimate the observations (by 50% in the GALEX NUV band, and even a factor 2 in the FUV band). The UV emission mainly originates from star forming regions, which are below the resolution limit for the EAGLE simulations. This resolution issue is handled using a subgrid approach, first employed by Jonsson et al. (2010): the star-forming regions are represented using template SEDs from the MAPPINGS library (Groves et al., 2008). These spherically symmetric models are controlled by different parameters, and it is well possible that the particular choice of these parameters can be further optimised. In particular, our calibration was based on submm colours and global dust scaling relations, and did not particularly focus on the UV wavelength regime (Camps et al., 2016; Camps et al., 2018). In addition, the limited spatial resolution in the gas component in EAGLE results in a more homogeneous ISM distribution than we expect in actual galaxies, and EAGLE galaxies have systematically thicker discs, yielding a puffed up interstellar medium (Trayford et al., 2017). It is well-known that inhomogeneities and clumping can easily cause differences in the UV attenuation of an order of magnitude or more (e.g., Witt & Gordon, 2000; Indebetouw et al., 2006; Saftly et al., 2015).

In summary, we believe that the discrepancies between the GAMA and EAGLE-SKIRT CSED at UV to NIR wavelengths are primarily due to the incompleteness of the spheroidal galaxy population in the EAGLE-SKIRT catalogue, and an underestimation of the UV attenuation in the radiative transfer post-processing recipe.

3.2 The infrared-submm CSED

In the far-infrared region, the GAMA data points are again systematically higher than the EAGLE-SKIRT data points. At 100 $\mu$m, the difference is 0.18 dex, between 160 and 350 $\mu$m it is reduced to about 0.1 dex, but at 500 $\mu$m it increases again to 0.34 dex. However, Andrews et al. (2017b) indicate that their CSED is only poorly constrained or partially extrapolated due to lack of data beyond 24 $\mu$m. We therefore also show the independent measurements of the infrared CSED as obtained by Marchetti et al. (2016), based on a very detailed analysis of Herschel and Spitzer data from the Herschel Multi-tiered Extragalactic Survey (HerMES: Oliver et al., 2012). Apart from the MIPS 24 $\mu$m band, where the agreement with GAMA was better, the EAGLE-SKIRT data agree much better with these independent measurements. In particular at the longest wavelengths, nearly perfect agreement is found between EAGLE-SKIRT and the Marchetti et al. (2016) HerMES data.

At far-infrared wavelength, i.e. between 60 and 250 $\mu$m, a small offset below 0.1 dex is observed. Given that the CSED at UV wavelengths overestimates the observations (Section 3.1), it seems likely that an underestimation of the UV attenuation in the radiative transfer post-processing recipe is the main reason for this offset. An additional factor can be the insensitivity of the RecalL0025N0752 simulation to luminous infrared sources (due to the limited simulation volume).

3.3 Contribution of different populations

In Figure 2 we split the local EAGLE-SKIRT CSED into contributions of galaxies in different bins in stellar mass (top row), star formation rate (middle row), and specific star formation rate (bottom row).

The top row shows that the CSED over the entire UV–submm wavelength range is dominated by galaxies within the stellar mass bin with $10<\log(M_{\star}/M_{\odot})<10.5$). The more massive galaxy population with $\log(M_{\star}/M_{\odot})>10.5$ comes second in the optical–NIR range. Note that the galaxies in this most massive bin are underrepresented in the EAGLE-SKIRT database due to the threshold on the number of dust particles, and in general, are underrepresented in the RecalL0025N0752 simulation because of the small volume (Furlong et al., 2015; Trayford et al., 2015). As expected, these most massive galaxies have a more modest contribution in the FIR-submm range, and a particularly low contribution in the UV. The population of increasingly lower-mass galaxies has an increasingly smaller contribution to the CSED, in spite of their increasing numbers. The lowest stellar mass bin has the smallest contribution along the entire spectrum.

The middle row shows that regular star-forming galaxies with a modest SFR around 1 $M_{\odot}$ yr*-1* contribute the bulk of the CSED over the entire wavelength range. Especially in the optical–NIR range they are clearly the most dominant contributor. The small population of galaxies in the highest SFR bin, with ${\text{SFR}}>1.8~{}M_{\odot}\,{\text{yr}}^{-1}$ has an almost equal contribution in MIR–FIR range. The population of galaxies with ${\text{SFR}}<0.06~{}M_{\odot}\,{\text{yr}}^{-1}$ has an almost negligible contribution to the CSED, even at optical–NIR wavelengths (but a fraction of these passive galaxies is missing in the EAGLE-SKIRT database).

The bottom row splits the contributions by populations in different sSFR bins. In this case, the contribution by galaxies centred around ${\text{sSFR}}\sim 10^{-10}~{}{\text{yr}}^{-1}$ is strongly dominant over the entire spectrum. This bin also corresponds to the largest number of galaxies. Interestingly, galaxies with a higher sSFR, although much less numerous, are the second-most important contributor in the UV and infrared range, while more passive galaxies with ${\text{sSFR}}\sim 10^{-10.5}~{}{\text{yr}}^{-1}$ contribute more to the optical–NIR CSED. Similarly, the galaxies with the lowest sSFR have the smallest contribution to the CSED in the UV and infrared range, whereas they contribute more in the optical–NIR range than the population of slightly higher sSFR galaxies.

The bottom line is that main sequence galaxies with $M_{\star}\sim 10^{10.25}~{}M_{\odot}$ and ${\text{sSFR}}\sim 10^{-10}~{}{\text{yr}}^{-1}$ dominate the local EAGLE-SKIRT CSED per dex in stellar mass or sSFR.

3.4 Implications

Overall, besides some minor differences that can readily be interpreted and explained, we find that the EAGLE-SKIRT CSED reproduces both the shape and the normalisation of the GAMA CSED very well over the entire UV-submm range. On the one hand, this is remarkable, given that no specific fitting or finetuning has been applied at this stage. This agreement can be interpreted as a confirmation that the combination of the EAGLE simulation and the SKIRT radiative transfer post-processing provides a solid representation of the present-day Universe.

One factor that we have not yet taken into account is the sample variance. Andrews et al. (2017b) estimate the total uncertainty on their CSED to be of the order of 20 per cent, which is almost completely due to sample variance. We estimated the sample variance in our CSED estimate by splitting the 100 Mpc simulation volume into 64 individual 25 Mpc volumes, and calculating the variance based on the CSEDs of these 64 volumes. The level of sample variance obtained this way is around 65% over the entire wavelength range, with only minor differences between the different wavelength regimes. This level of sample variance is consistent with the estimate of 56% cosmic variance based on the empirical formulae from Driver & Robotham (2010) for a cubical 25 Mpc volume. Given these levels of uncertainty, the agreement between the GAMA and EAGLE-SKIRT CSEDs is impressive.

On the other hand, the close agreement in the local Universe might not be too surprising. The subgrid physics parameters in the EAGLE simulations were calibrated to reproduce a number of characteristics of the local Universe, including the local stellar mass function (Crain et al., 2015; Schaye et al., 2015). In addition, the few free parameters in the SKIRT post-processing procedure were calibrated to reproduce the submm color-color relations and dust scaling relations as observed in the local Universe (Camps et al., 2016). The fact that this combination now also reproduces the UV-submm CSED to a large degree is a nice confirmation that these calibrations are fairly solid.

3.5 Bolometric energy output of the local Universe

From the estimates of the CSED in individual broadband filter bands, the logical next step is the determination of the bolometric energy output of the Universe222The quantity $\varepsilon_{\text{bol}}$ has the dimension energy per unit volume. However, it does not represent the bolometric radiative energy density, because it only contains the radiative energy generated in a representative unit volume, and not the energy connected to photons passing through this volume. This subtle distinction is important when considering the extragalactic background light (e.g., Andrews et al., 2017b; Andrews et al., 2018; Cowley et al., 2019).,

[TABLE]

where the integral covers the entire UV to submm wavelength range. This requires some interpolation scheme, and because of the complicated nature of translating broadband photometry into monochromatic flux densities, this integration is best done via SED fitting techniques (for a discussion, see Brown et al., 2016).

The solid grey line in Figure 1 is a panchromatic energy balance template fit to the observed GAMA CSED data points, based on the MAGPHYS code (da Cunha et al., 2008; Driver et al., 2018), and provided in tabular format by Andrews et al. (2017b). The line fits all of the GAMA data points very well, except the SPIRE 500 $\mu$m data point that seems to be incompatible with the typical dust emission profile. From this fit, Andrews et al. (2017b) recover a value of $\varepsilon_{\text{bol}}=1.26\times 10^{35}$ W Mpc*-3*.

We repeated this analysis for our EAGLE-SKIRT CSED, but we used CIGALE, another popular SED fitting package (Noll et al., 2009; Boquien et al., 2019). We used the CIGALE version 2018.0, equipped with the same model ingredients and parameter settings as Nersesian et al. (2019). For the stellar populations, we have adopted the Bruzual & Charlot (2003) library of single stellar populations with a Salpeter (1955) IMF, and a delayed and truncated star formation history model (Ciesla et al., 2016). Nebular emission was included, based on CLOUDY templates, and under the assumption that all the ionising radiation is absorbed by gas (Ferland et al., 1998; Inoue, 2011). For the dust, we adopt the THEMIS dust model (Jones et al., 2017; Davies et al., 2017), and a modified version of the standard starburst-like attenuation (Calzetti et al., 2000; Boquien et al., 2016). In total, about $8\times 10^{7}$ different models were considered in the fitting.

The solid green line in Figure 1 is the best fitting model. Integrating the CSED over the entire UV-submm wavelength range, we recover a value of $\varepsilon_{\text{bol}}=9.94\times 10^{34}$ W Mpc*-3*. This value is 0.10 dex below the GAMA value as measured by Andrews et al. (2017b). This difference is mainly due to the difference in SED shape between the two models in the ill-covered region between 24 and 100 $\mu$m. The MAGPHYS fit to the GAMA CSED shows a remarkable “boxy” shape in this part of the spectrum. We believe that this shape is due to the setup of the MAGPHYS code, which models the dusty medium as a combination of four different components with adjustable temperatures (da Cunha et al., 2008). Such a behaviour has been noted in MAGPHYS SED fits to individual galaxies, where other panchromatic SED fitting codes do not show this feature (e.g., Pappalardo et al., 2016; Hunt et al., 2019).

The dashed grey line in Figure 1 is our own CIGALE fit to the GAMA CSED data points from Andrews et al. (2017b). This SED model represents an equally good fit to the observed broadband data points, and avoids the boxy shape between 24 and 100 $\mu$m. Integrating the CSED now results in $\varepsilon_{\text{bol}}=1.08\times 10^{35}$ W Mpc*-3*, a difference of just 0.03 dex with our EAGLE-SKIRT value. To put these numbers in context: the bolometric energy output of the Universe corresponds to a single 50 W light bulb within a sphere of radius 1 AU.

In conclusion, the total EAGLE-SKIRT energy output agrees very well with the GAMA observations, in spite of the differences in the CSED at UV and infrared wavelengths. This indicates that the underestimation of the EAGLE-SKIRT CSED at FIR wavelengths is mainly due to an underestimation of the UV attenuation, as suggested in Section 3.2.

3.6 Physical global characteristics of the local Universe

The SED modelling exercise presented in the previous subsection was necessary to determine the total energy output of the Universe. But rather than just using it as a sophisticated integrator, we can use it to estimate a number of fundamental physical characteristics of the local Universe. Indeed, each SED model in the CIGALE library is characterised by a large number of physical parameters, such as the star formation rate, the stellar mass and dust mass. In fact, the goal of SED fitting usually is the determination and analysis of these parameters for a sample of galaxies (e.g. Smith et al., 2012a; Chang et al., 2015; Boquien et al., 2016; Pappalardo et al., 2016; Driver et al., 2018). In a similar way, we can interpret the physical parameters of our best fitting SED model to the CSED data points to estimate the quantities as the cosmic star formation rate density, stellar mass density, and dust mass density of the Universe, and compare them to measurements directly based on the EAGLE particle data.

This approach is not completely rigorous. In order to properly estimate, for example, the stellar mass density of the Universe, one should, for each individual galaxy in a representative cosmic volume, estimate the stellar mass using an SED fitting code, and subsequently sum all of these contributions. This will not necessarily yield the same answer as the stellar mass corresponding to the best fitting SED model to the sum of the individual galaxy SEDs. However, our CIGALE fit provides us with an easy and convenient shortcut to estimate these quantities, and it is instructive to compare these measurements to the actual values known in EAGLE. The most important characteristics are listed in Table 2. The third and fourth columns correspond to our CIGALE fits to the EAGLE-SKIRT and the GAMA CSED, respectively. The fifth column lists the values obtained by Driver et al. (2018), based on MAGPHYS SED fits to individual galaxies from the combined GAMA, G10-COSMOS and 3D-HST surveys.

For the cosmic stellar mass density, we find a value $\rho_{\star}=1.9\times 10^{8}~{}M_{\odot}$ Mpc*-3*, or equivalently $\Omega_{\star}=0.0015$, in almost scarily good agreement with the results from Driver et al. (2018). The value we obtain from our CIGALE fit to the GAMA CSED data points is 20% higher. These numbers are also fully consistent with independent estimates of the local stellar mass density based on the same GAMA data (Moffett et al., 2016; Wright et al., 2017), and fall nicely within the range of estimates based on independent studies (Cole et al., 2001; Bell et al., 2003; Panter et al., 2004; Eke et al., 2005; Li & White, 2009; Moustakas et al., 2013). It is also reassuring that our $\Omega_{\star}$ value determined from the CSED is completely consistent with the value obtained directly from the EAGLE stellar mass function: at $z\sim 0.05$, Furlong et al. (2015) found a value $\Omega_{\star}=0.0016$. Given the completely different methodology, the different initial mass function333The EAGLE simulations use a Chabrier (2003) IMF, whereas a Salpeter (1955) IMF is used for the CIGALE SED fitting., and the fact that the uncertainties on derived stellar masses due to stellar evolution models, even at a fixed IMF, are typically $\sim$0.3 dex (Conroy et al., 2009; Pforr et al., 2012), this agreement is remarkable. Finally, when we estimate the "intrinsic" value for $\Omega_{\star}$ by simply summing the actual stellar masses of all galaxies in the 25 Mpc volume, we recover the value 0.0012, again in very good agreement.

For the cosmic star formation rate density, one of the most fundamental parameters in galaxy formation and evolution models, we find a value of 0.017 $M_{\odot}$ yr*-1* Mpc*-3*. This value is 50% higher than the GAMA value (0.011 $M_{\odot}$ yr*-1* Mpc*-3*) reported by Driver et al. (2018), which agrees perfectly with the intrinsic value obtained by summing the individual SFRs of all galaxies in the EAGLE 25 Mpc volume (Katsianis et al., 2017). These values bracket the most credible recent literature values (Robotham & Driver, 2011; Gunawardhana et al., 2013; Madau & Dickinson, 2014; Davies et al., 2016). Given the lack of rigor in our method, and the intrinsic uncertainties in estimating SFRs from SED modelling due to the sensitive dependence on the assumed star formation history (Hunt et al., 2019; Carnall et al., 2019; Leja et al., 2019), we find this agreement very satisfactory.

The most significant deviation between our results and the GAMA results corresponds to dust-related parameters. For the cosmic dust mass density we find a value $\rho_{\text{d}}=8.3\times 10^{4}~{}M_{\odot}~{}{\text{Mpc}}^{-3}$, or equivalently $\Omega_{\text{d}}=6.5\times 10^{-7}$. This is almost 75% lower than the GAMA value of $1.1\times 10^{-6}$ reported by Driver et al. (2018), and at the bottom end of independent estimates for $\Omega_{\text{d}}$ in the local Universe, which range between $7\times 10^{-7}$ and $2.7\times 10^{-6}$ (Dunne et al., 2011; Clemens et al., 2013; Clark et al., 2015; Beeston et al., 2018). The low value we find is partly due to the underestimation of the UV attenuation and the resulting dearth of infrared emission, and to the incompleteness of our EAGLE-SKIRT catalogue. An additional factor is that the CIGALE fit to the CSED, on which our value for $\Omega_{\text{d}}$ is based, systematically lies below the EAGLE-SKIRT data points from 350 $\mu$m onwards (Figure 1). This is a consequence of the infrared templates using in our CIGALE setup, which assume a power-law distribution in the radiation field strength (Dale et al., 2001; Galliano et al., 2011). This is suitable to fit the SEDs of individual galaxies, but apparently less ideal to fit the broader CSED that results from the sum of many individual SEDs with different cold dust temperatures.

The dust-to-stellar mass ratio $\Omega_{\text{d}}/\Omega_{\star}$ is a valuable tool to probe the evolution of galaxies, as it represents an observable measure of how much dust per unit stellar mass survives the various destruction processes in galaxies. Theoretical models outline the strong dependence of this quantity on the underlying star formation history (Santini et al., 2014; McKinnon et al., 2017; Calura et al., 2017). For the dust-to-stellar-mass ratio of the local Universe we find $4.4\times 10^{-4}$, in perfect agreement with the value obtained from our CIGALE fit to the GAMA data. The ratio based on the MAGPHYS analysis of GAMA galaxies from Driver et al. (2018) is, not surprisingly, almost 70% higher. Observed dust-to-stellar mass ratios in individual galaxies vary over a wide range of values: the typical values for spiral galaxies are usually found in the range between $5\times 10^{-4}$ and 0.01, whereas the values for early-type galaxies go down to $10^{-5}$ and below (Skibba et al., 2011; Smith et al., 2012b; Cortese et al., 2012; Davies et al., 2012).

The FUV attenuation is probably the most important difference between the EAGLE-SKIRT simulations and the GAMA observations. Both the CIGALE and MAGPHYS fits to the GAMA observations indicate $A_{\text{FUV}}\sim 1.5$, whereas our EAGLE results yield an FUV attenuation below one magnitude. As discussed in Section 3, the limited resolution of the EAGLE simulation, even in the higher-resolution 25 Mpc volume, combined with limitations in our EAGLE-SKIRT subgrid model for star-forming regions and the thickness of the ISM in the EAGLE subgrid physics, are probably largely responsible. In this context, it is important to mention that most SED fitting codes, including CIGALE and MAGPHYS, adopt a rather simple treatment of attenuation, which is essentially a one- or two-component absorbing screen model. Our SKIRT radiative transfer modelling of each EAGLE galaxy, on the other hand, computes the attenuation in a fully self-consistent way, including multiple anisotropic scattering and in a realistic 3D setting. On the level of individual galaxies, this can yield fairly different effective attenuation curves, as demonstrated using idealised models (Witt & Gordon, 2000; Baes & Dejonghe, 2001; Inoue, 2005), and using EAGLE galaxies (Trayford et al., 2015, 2017). Given that the stellar populations and ISM properties of galaxies vary systematically with mass, we believe that the FUV attenuation as derived from an SED fit to the global CSED should be taken with a grain of salt. A more detailed analysis of the FUV attenuation, based on SED fits to individual EAGLE galaxies, will be considered in future work.

An interesting quantity that can be calculated from the CIGALE fitting is the cosmic bolometric attenuation $f_{\text{abs}}$, defined as the fraction of the total energy output that is absorbed and re-emitted by dust. We find a value of 27% for our EAGLE-SKIRT simulation, in very good agreement with the CIGALE fit to the GAMA data. These values are also very similar to the typical bolometric attenuation values found for individual star-forming galaxies of the local Universe. Viaene et al. (2016) recently performed a detailed investigation of the bolometric attenuation of 239 late-type galaxies from the Herschel Reference Survey (Boselli et al., 2010), and found an average value of 32%. They noted that this number is remarkably similar to the value obtained from previous studies, even including studies that hardly had any access to UV or submm data (Soifer & Neugebauer, 1991; Xu & Buat, 1995; Popescu & Tuffs, 2002; Skibba et al., 2011). Bianchi et al. (2018) present a larger study of the bolometric attenuation for more than 800 galaxies of different morphological types from the DustPedia sample (Clark et al., 2018). They find a slightly lower average value ($\langle f_{\text{abs}}\rangle=19\%$), but note that their sample is missing high-luminosity and high-specific-star-formation-rate objects. When only considering late-type galaxies, they find a average value $\langle f_{\text{abs}}\rangle=25\%$, very close to the mean value we find in our analysis.

4 Comparison of different EAGLE simulations

All EAGLE-SKIRT results presented so far were based on the high-resolution RecalL0025N0752 simulation. In this section we compare these results to EAGLE-SKIRT CSEDs based on a number of other EAGLE simulations. In all cases, the CSEDs were calculated using exactly the same procedure; in particular, all galaxies from each run were post-processed using SKIRT with the same “standard” parameter values, as discussed by Camps et al. (2018).

4.1 Strong and weak convergence

The top-left panel of Figure 3 compares the CSEDs of the RecalL0025N0752 and RefL0100N1504 simulations. The latter simulation is the reference EAGLE simulation, which has a volume 64 times larger, but a mass resolution eight times worse than the high-resolution RecalL0025N0752 run. Since the subgrid physics parameters in both simulations were independently calibrated to reproduce the galaxy properties in the local Universe, the comparison of both models can be regarded as a weak convergence test (Schaye et al., 2015). Such weak convergence tests have been done for various aspects of the EAGLE simulations, including the global stellar mass and SFR density (Furlong et al., 2015), the relation between stellar mass and angular momentum (Lagos et al., 2017), and optical luminosity functions (Trayford et al., 2015).

It is immediately obvious that the RefL0100N1504 results systematically underestimate the RecalL0025N0752 results, and hence also the observed GAMA CSED. The difference between the two EAGLE CSEDs is on average about 0.2 dex, and reaches up to 0.4 dex at UV wavelengths. This can largely be explained by the incompleteness of the EAGLE-SKIRT database. For the RefL0100N1504 simulation, the threshold of 250 dust particles translates to a dust mass of about $4\times 10^{6}~{}M_{\odot}$. Both low-mass late-type galaxies (Grossi et al., 2010, 2015; Skibba et al., 2011) and massive elliptical galaxies (Smith et al., 2012b; di Serego Alighieri et al., 2013) in the local Universe often have dust masses below this threshold.

The top-left panel of Figure 4 compares the same CSEDs as the top-left panel of Figure 3, but now the curves are normalised to the stellar mass density of the RecalL0025N0752 model. To some extent, this eliminates the differences in incompleteness of the EAGLE-SKIRT catalogues, and can highlight additional effects. Not surprisingly, the CSEDs for both models now agree nearly perfectly in the near-infrared. The largest differences between the normalised CSEDs are seen in the UV, with the RecalL0025N0752 model still 0.3 dex higher than the RefL0100N1504 model. This can be understood as the former model has a higher specific star formation rate density than the former, and thus a larger intrinsic UV output. In general, due to resolution effects, the RefL0100N1504 simulation is characterised by a relative overabundance of red/passive galaxies at $M_{\star}<10^{9}~{}M_{\odot}$ (Furlong et al., 2015; Trayford et al., 2015). In spite of the higher specific star formation rate, the normalised FIR CSED of the RecalL0025N0752 model is not much higher than that of the RefL0100N1504 model. This might be due to the lower metallicity, and hence dust content, of RecalL0025N0752 galaxies: indeed, this run has a lower, and more realistic, mass-metallicity relation (Schaye et al., 2015).

The top-right panel of Figure 3 can be regarded as a strong convergence test. This plot compares the CSEDs for models with exactly the same physical subgrid parameters (the parameters of the largest EAGLE volume) and the same simulation volume (25 Mpc on a side), but with a different resolution. The RefL0025N0376 simulation has the same resolution as the largest EAGLE volume simulation (RefL0100N1504), whereas the mass and spatial resolution of the RefL0025N0752 differ by factors of eight and two, respectively. It is no surprise that we see more or less the same effect as in the top-left panel: the higher threshold for the dust mass for the lower-resolution simulation causes an underestimation of the CSED over the entire wavelength range. Note, however, that this is most probably not the only reason: EAGLE simulations (and other cosmological hydro simulations) do not score well on strong convergence tests, which was exactly the reason why resolution-dependent recalibration has been implemented (Crain et al., 2015; Schaye et al., 2015).

The top-right panel of Figure 4 shows the corresponding version of this plot, but now again normalised to the stellar mass density of the RecalL0025N0752 model. A similar effect is noted as in the top-left panel: the higher resolution of the RefL0025N0752 simulation leads to a higher specific star formation rate, because feedback is less efficient. The result is a higher normalised CSED compared to the lower-resolution RefL0025N0376 simulation at UV and FIR wavelengths.

4.2 Variation of the subgrid model parameters

The bottom-left panel of Figure 3 shows the effect of the resolution-dependent recalibration of the EAGLE subgrid physics parameters. This panel compares the CSEDs corresponding to two EAGLE simulations with exactly the same volume and resolution (the highest resolution of all EAGLE runs), but with different subgrid parameters. RefL0025L0752 uses the same subgrid physical parameters as the standard 100 Mpc simulation, whereas the RecalL0025L0752 simulation has different values for a number of subgrid physics parameters, including the characteristic density for star formation, and the temperature increase of the gas during AGN feedback (Schaye et al., 2015).

The effects of this recalibration are clearly visible: the RefL0025L0752 simulation systematically gives larger values for the CSED than the RecalL0025L0752 simulation, over the entire wavelength range. The difference between both CSEDs is roughly 0.2 dex, with the largest difference (0.25 dex) at submm wavelengths. This is no surprise, given that the intrinsic stellar mass density and the SFR density of the RefL0025L0752 simulation are respectively 0.15 dex and 0.20 dex higher than the corresponding values of the RecalL0025L0752 simulation (see also Furlong et al., 2015). When these CSEDs are normalised to the same stellar mass density (bottom-left panel of Figure 4), the differences become almost negligible. The normalised CSEDs of both models are almost identical over the entire UV–NIR wavelength range. The RecalL0025L0752 simulation has a slightly lower FIR-submm normalised CSED, because of the lower mean metallicity and hence dust content.

Finally, the bottom-right panel compares the CSEDs of two EAGLE runs with the same simulation volume and resolution, but with a different parameterisation of the AGN feedback. The RefL0050N0752 simulation uses the standard subgrid physics parameterisation, whereas the AGNdT9L0050N0752 model has adjusted AGN parameters in order to further improve the agreement with observations for high-mass galaxies. The main difference is an increased gas heating temperature increase for AGN feedback, which corresponds to more energetic but less frequent bursts, and more effective gas ejection. This improves the comparison to X-ray observations of the intracluster medium, at least on the scales of groups of galaxies (Schaye et al., 2015).

On the scales of individual galaxies, this change in the subgrid physics parameters does not have a strong effect. For example, the intrinsic stellar mass density of both models is identical, and the SFR density of the AGNdT9L0050N0752 model is only slightly higher (0.05 dex). As a result, both CSEDs are almost identical in the UV–NIR range, and the AGNdT9L0050N0752 has a slightly increased FIR CSED compared to the RefL0050N0752 CSED. The latter effect might be due to the fact that massive galaxies in the AGNdT9L0050N0752 run are slightly more compact (Schaye et al., 2015), and as a result, their dust will be slightly warmer on average. Due to the relatively poor resolution (a mass resolution eight times worse compared to the high-resolution simulation), the corresponding incompleteness of the EAGLE-SKIRT catalogue, and the lack of recalibration of these simulations, both CSEDs underestimate the observed GAMA CSED over the entire wavelength range.

5 Cosmic evolution of the CSED

The comparison of the local Universe CSED provides a powerful test for galaxy evolution models, but an even more powerful test is the evolution of the CSED with cosmic time. Figure 5 shows the CSED at four different redshifts between 0 and 1. As in Figure 1, the grey squares correspond to the observed GAMA CSED from Andrews et al. (2017b). For the redshift bins $0.02<z<0.08$ and $0.14<z<0.20$, the CSED is based on data from the standard equatorial GAMA fields (Driver et al., 2011, 2016; Liske et al., 2015). For the redshift bins $0.45<z<0.56$ and $0.82<z<0.92$, the observed CSED is based on the G10/COSMOS data (Davies et al., 2015; Andrews et al., 2017a). The solid grey lines are MAGPHYS fits to the observed CSED points, as provided by Andrews et al. (2017b).

The coloured bullets correspond to the EAGLE-SKIRT simulation, derived in the same way as described in Section 3. It is immediately obvious that the nice agreement between the EAGLE-SKIRT and GAMA results in the local Universe does not continue to higher redshifts. At $\langle z\rangle=0.17$ and $\langle z\rangle=0.50$ the underestimation of the CSED in the red and NIR domain is somewhat stronger than in the local Universe. The main issue, however, is the FIR-submm range, where the EAGLE-SKIRT results significantly underestimate the observed GAMA CSED. In the highest redshift bin corresponding to $\langle z\rangle=0.91$, the disagreement is even worse. In the IRAC 1 band (corresponding to a rest-frame wavelength of $\sim 2~{}\mu$m), the EAGLE-SKIRT estimate is a factor 5 lower than the observed data point, and the underestimation in the FIR/submm dust emission peak is even up to a factor 10.

We note that a similar trend was also found by Andrews et al. (2018): they find that the GALFORM semi-analytical model can reproduce the CSED for $z<0.3$ with a fairly good agreement, but the agreement becomes increasingly poor towards $z=1$.

The increasing disagreement between our EAGLE-SKIRT CSED and the GAMA observations with increasing redshift is most probably due to a combination of different factors. Firstly, there is the incompleteness of the EAGLE-SKIRT catalogue due to the threshold of at least 250 dust particles per galaxy, which becomes worse for higher $z$ because galaxies are less massive on average (Camps et al., 2018). A second aspect that comes into play, particularly at FIR/submm wavelengths, is that luminous infrared sources contribute increasingly more when moving to higher redshift. Different surveys have clearly shown that the infrared/submm luminosity function has a strong and rapid cosmic evolution out to $z\sim 1$ (e.g. Eales et al., 2009, 2010; Dye et al., 2010). Marchetti et al. (2016) demonstrated that the characteristic luminosity in the SPIRE 250 $\mu$m band, $L_{250}^{\star}$, evolves as $(1+z)^{5.3}$. For the total infrared luminosity, they found an even stronger evolution, with $L_{\text{IR}}^{\star}\propto(1+z)^{6.0}$. In order to properly incorporate the contribution of rare but luminous galaxies to the CSED, the use of the small EAGLE-SKIRT simulation box is not ideal.

It is, however, very unclear whether these two aspects can explain the differences, and other aspect might very well also contribute to, or even dominate, this disagreement. As already discussed, the subgrid physics in the EAGLE simulations were calibrated based on observed relations in the local Universe, and hence are not necessarily optimised for the Universe at higher redshift (Crain et al., 2015; Schaye et al., 2015). We do note, however, that the EAGLE simulations show a reasonable level of agreement with the evolution of the galaxy stellar mass function (Furlong et al., 2015), the mass-size relation (Furlong et al., 2017), and the cosmic star formation rate density (Katsianis et al., 2017).

A similar argument applies to the EAGLE-SKIRT post-processing radiative transfer procedure. The calibration of this procedure was based on submm properties of galaxies in the local Universe (Camps et al., 2016), and turned out to be successful in also reproducing the optical colours in the local Universe (Trayford et al., 2017). However, as indicated by Camps et al. (2018), this does not guarantee that these settings are also ideal for higher redshifts. In particular, for the construction of the EAGLE-SKIRT database, the same dust properties and fixed dust-to-metal ratio were adopted at all redshifts. It is uncertain, however, whether these assumptions are realistic, as different studies on this topic yield different conclusions. Zafar & Watson (2013) performed a detailed study of the dust-to-metal ratio across a range of redshifts and sources, and do not find any obvious dependence of the dust-to-metals ratio on column density, galaxy type, redshift, or metallicity. Several other studies, including Brinchmann et al. (2013) and De Vis et al. (2017, 2019), do reveal a systematic evolution of the dust-to-metal ratio of galaxies. It is hence possible that our assumptions on the dust properties of the EAGLE galaxies are flawed, in particular towards higher redshift.

Dunne et al. (2011) argued that the observed evolution of the dust mass function is difficult to explain using standard dust evolution models and requires a combination of various processes, including a possible evolution of dust grain properties. Beeston et al. (2018) compares the local dust mass function with the predictions of the semi-analytical model of Popping et al. (2017) and the cosmological hydrodynamical model of McKinnon et al. (2016, 2017), both of which incorporate recipes for dust production and destruction. Both sets of theoretical predictions have difficulties in reproducing the dust mass function at both the low and high dust mass regimes, indicating that fundamental properties as stardust condensation efficiencies and dust grain growth time-scales are poorly understood. Beeston et al. (2018) conclude that the current theoretical models for dust evolution need to be improved. We agree that more work is required in this field, and we particularly hope that a more integrated approach combining cosmological hydrodynamical simulations and dust evolution modelling (Bekki, 2015; Zhukovska et al., 2016; McKinnon et al., 2016, 2017, 2018; Aoyama et al., 2018) will move this field forward.

6 Summary and outlook

We have presented, for the first time, the local Universe cosmic spectral energy distribution (CSED) as derived from a cosmological hydrodynamical simulation (EAGLE RecalL0025L0752), combined with a dust radiative transfer post-processing algorithm (SKIRT). The CSED was obtained by simply summing the flux densities of individual galaxies in the EAGLE-SKIRT database recently presented by Camps et al. (2018). Overall, we find an excellent agreement between the EAGLE-SKIRT CSED and the observed CSED based on GAMA observations, with some relatively minor differences:

• •

In the UV regime, the EAGLE-SKIRT CSED overestimates the observations, or conversely, underestimates the attenuation. This is likely a result of the limited resolution of the EAGLE simulations and the limitations in the SKIRT subgrid treatment of star forming regions.

• •

At optical and near-infrared wavelengths, the EAGLE-SKIRT CSED slightly underestimates the GAMA observations by about 0.13 dex. This is due to a combination of the incompleteness of the EAGLE-SKIRT database, which excludes galaxies with dust masses below $\sim 5\times 10^{5}~{}M_{\odot}$, the small volume of the EAGLE-SKIRT simulation and the resulting insensitivity for luminous galaxies, and a small EAGLE underestimation of the stellar mass function for massive galaxies.

• •

In the far-infrared and submm regime, our EAGLE-SKIRT results significantly underestimate the GAMA data points, which are only poorly constrained. This discrepancy largely disappears when our results are compared to independent estimates of the FIR/submm CSED based on deep Spitzer and Herschel data. The remaining underestimation of $\sim 0.1$ dex is probably mainly due to the underestimation of the attenuation at UV wavelengths.

• •

Splitting the local EAGLE-SKIRT CSED into different stellar mass, SFR and sSFR populations, we find that main sequence galaxies with $M_{\star}\sim 10^{10.25}~{}M_{\odot}$ and ${\text{sSFR}}\sim 10^{-10}~{}{\text{yr}}^{-1}$ are the dominant contributors to the local CSED over the entire UV–submm wavelength range.

Based on a physically motivated SED fit with the CIGALE code, we derive a number of global characteristics of the local Universe.

• •

Our EAGLE-SKIRT estimate of the bolometric energy output of the local Universe is in excellent agreement with the value obtained from GAMA observations. It is equivalent to a single 50 W light bulb in a sphere with radius 1 AU.

• •

The cosmic star formation rate density derived from our CIGALE fit to the EAGLE-SKIRT CSED is about 50% higher than the GAMA value, but perfectly within the range of independent values quoted in the recent literature for $z\sim 0.05$.

• •

For the cosmic dust density and the cosmic dust-to-stellar-mass ratio, we find values that are some 70% lower than the GAMA based estimates. This is due to a combination of the underestimation of the UV attenuation, the incompleteness of the EAGLE-SKIRT database, and the difficulties CIGALE has to properly fit the CSED data points at the longest wavelengths.

• •

We obtain the result that 27% of all radiative energy emitted by stars in the local Universe is absorbed by dust and re-emitted as thermal radiation in the infrared regime. This number is in very good agreement with the typical numbers found for the bolometric attenuation in individual galaxies.

Putting this all together, we interpret this excellent agreement between the simulated and observed local CSED as a confirmation that the combination of the EAGLE simulation and the SKIRT radiative transfer post-processing provides a reliable mock representation of the present-day Universe.

Unfortunately, the nice agreement that we found in the local Universe does not hold at higher redshifts. Even at $\langle z\rangle=0.17$, we find significant deviations between the EAGLE-SKIRT results and the GAMA results, and the magnitude of these differences only increases out to $z\sim 1$. The most important deviations are an increasing underestimation of the optical-NIR CSED, and an even greater discrepancy at far-infrared wavelengths. We believe that this discrepancy can be attributed to a combination of factors, including the incompleteness of the EAGLE-SKIRT database, the limited volume of the EAGLE 25 Mpc simulation (and hence the poor coverage of the high-luminosity end of the luminosity function), and our lack of knowledge on the evolution of the characteristics of the interstellar dust in galaxies. Indeed, one important caveat in our methodology is the assumption that the dust properties and the dust-to-metal ratio do not vary with increasing redshifts, but this assumption might be far too simplistic.

As the CSED is nothing but the joint contribution of the individual spectral energy distributions of all galaxies within a cosmological unit volume, it is a strong constraint for models of galaxy formation and evolution. We have shown in this paper that the EAGLE cosmological simulation, when combined with realistic mock observations based on detailed radiative transfer modelling, successfully withstands the comparison to the observed CSED. On the other hand, the CSED is not the ultimate test. A more stringent test would be to compare the EAGLE simulation to the CSED split up in different luminosity bins, i.e. luminosity functions covering the entire UV to submm wavelength range. Ideally, the luminosity functions could take into account information of EAGLE simulations at different resolutions, where the higher-resolution simulations constrain the low-luminosity regime, and the large volume simulations the high-luminosity tail. This is beyond the scope of this paper, and will be considered for future work.

Acknowledgements

The authors thank the anonymous referee for constructive comments that improved this paper. M.B., A.T., A.N., and W.D. gratefully acknowledge support from the Flemish Fund for Scientific Research (FWO-Vlaanderen), and from DustPedia, a European Union’s Seventh Framework Programme (FP7) for research, technological development and demonstration under grant agreement no. 606874.

The EAGLE-SKIRT database on which this work was based used the DiRAC Data Centric system at Durham University, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (http://www.dirac.ac.uk). This equipment was funded by BIS National E-infrastructure capital grant ST/K00042X/1, STFC capital grants ST/H008519/1 and ST/K00087X/1, STFC DiRAC Operations grant ST/K003267/ 1, and Durham University. DiRAC is part of the National E-Infrastructure.