From the far-ultraviolet to the far-infrared -- galaxy emission at $0\le z \le 10$ in the Shark semi-analytic model

TL;DR

This study combines a semi-analytic galaxy formation model with spectral energy distribution tools to accurately simulate galaxy emission across a wide range of wavelengths and redshifts, reproducing observed luminosity functions and cosmic SEDs without adjusting fundamental parameters.

Contribution

It introduces a comprehensive physical model for galaxy attenuation over cosmic time, integrating empirical dust relations with radiative transfer data, and successfully reproduces a broad set of observational galaxy properties from z=0 to z=10.

Findings

Shark model reproduces observed luminosity functions across multiple wavelengths and redshifts.

The model predicts star formation and dust attenuation behaviors consistent with observations.

Galactic disk and bulge contributions to emission evolve with redshift, with disks dominating in FUV-optical and bulges in NIR.

Abstract

We combine the Shark semi-analytic model of galaxy formation with the ProSpect software tool for spectral energy distribution (SED) generation to study the multi-wavelength emission of galaxies from the far-ultraviolet (FUV) to the far-infrared (FIR) at $0\le z\le 10$. We produce a physical model for the attenuation of galaxies across cosmic time by combining a local Universe empirical relation to compute the dust mass of galaxies from their gas metallicity and mass, attenuation curves derived from radiative transfer calculations of galaxies in the EAGLE hydrodynamic simulation suite, and the properties of Shark galaxies. We are able to produce a wide range of galaxies, from the $z=8$ star-forming galaxies with almost no extinction, $z=2$ submillimeter galaxies, down to the normal star-forming and red sequence galaxies at $z=0$. Quantitatively, we find that Shark reproduces the observed (i) the $z=0$ FUV-to-FIR, (ii) $0\le z\le 3$ rest-frame $K$-band, and (iii) $0\le z\le 10$ rest-frame FUV luminosity functions, (iv) $z\le 8$ UV slopes, (v) the FUV-to-FIR number counts (including the widely disputed 850$\mu$m), (vi) redshift distribution of bright $850\mu$m galaxies and (vii) the integrated cosmic SED from $z=0$ to $z=1$ to an unprecedented level. This is achieved without the need to invoke changes in the stellar initial mass function, dust-to-metal mass ratio, or metal enrichment timescales. Our model predicts star formation in galaxy disks to dominate in the FUV-to-optical, while bulges dominate at the NIR at all redshifts. The FIR sees a strong evolution in which disks dominate at $z\le 1$ and starbursts (triggered by both galaxy mergers and disk instabilities, in an even mix) dominate at higher redshifts, even out to $z=10$.

Tables

Table 1. Table 1: Attenuation models tested. “CF00” refers to Charlot & Fall ( 2000 ) , and “RR14” refers to Rémy-Ruyer et al. ( 2014 ) . The dependence of the CF parameters on Σ dust subscript Σ dust \Sigma_{\rm dust} is taken from the parametrisation of EAGLE galaxies by Trayford et al. ( 2019 ) . Our default option is the model EAGLE - τ 𝜏 \tau RR14.

Name	Description
CF00	Adopts default Charlot & Fall (2000)
	parameters.
EAGLE-$\tau$ $f_{\mathrm{dust}}$ const	Adopts CF parameters depending on ${\Sigma }_{\mathrm{dust}}$,
	using a constant $f_{\mathrm{dust}}$
EAGLE-$\tau$ RR14	Adopts CF parameters depending on ${\Sigma }_{\mathrm{dust}}$,
(default)	using the RR14 best-fit $f_{\mathrm{dust}}-Z_{\mathrm{gas}}$relation.
EAGLE-$\tau$ RR14-steep	Adopts CF parameters depending on ${\Sigma }_{\mathrm{dust}}$,
	using the RR14 $f_{\mathrm{dust}}-Z_{\mathrm{gas}}$ relation with
	a steeper slope

Equations

$$ϵ=\frac{V_{\mathrm{circ}}}{\sqrt{1.68G{M}_{\mathrm{disk}}/r_{\mathrm{disk}}}},$$

$${\Sigma }_{\mathrm{dust},\mathrm{disk}}=\frac{0.5M_{\mathrm{dust},\mathrm{disk}}}{\pi r_{50,\mathrm{d}}{l}_{50}},$$

$${\Sigma }_{\mathrm{dust},\mathrm{bulge}}=\frac{0.5M_{\mathrm{dust},\mathrm{bulge}}}{\pi r_{50,\mathrm{b}}^2},$$

$${\tau }_{\mathrm{ISM}}$$

$${\tau }_{\mathrm{BC}}$$

$${\tau }_{\mathrm{BC}}={\tau }_{\mathrm{BC},0}\left[\frac{f_{\mathrm{dust}}{Z}_{\mathrm{gas}}{\Sigma }_{\mathrm{gas},\mathrm{cl}}}{f_{\mathrm{dust},\mathrm{MW}}{Z}_{\odot }{\Sigma }_{\mathrm{MW},\mathrm{cl}}}\right],$$

Full text

††pubyear: 2019††pagerange: From the far-ultraviolet to the far-infrared - galaxy emission at $0\leq z\leq 10$ in the Shark semi-analytic model–From the far-ultraviolet to the far-infrared - galaxy emission at $0\leq z\leq 10$ in the Shark semi-analytic model

1 Introduction

Galaxy formation and evolution is one of the most outstanding questions in astrophysics. Galaxies are thought to form in the centre of the gravitational potential of dark matter (DM)-dominated halos, and hence are significantly affected by the growth of structures in the Universe. They are also subject to highly non-linear, complex astrophysical processes, such as gas accretion, star formation, energetic events that change the thermodynamics of the gas, just to mention a few (see Somerville et al. 2015 for a review on the topic). The clues we get about how galaxies form and evolve come mostly from the electromagnetic spectrum produced by the integrated contribution of gas, dust and stars in galaxies.

This integrated electromagnetic spectrum, also called spectral energy distribution (SED), encodes information of a galaxy’s stellar populations, via the light emitted by stars, as well as its interstellar medium (ISM) (both in terms of content and composition) through the absorption of the far-ultraviolet (FUV)-to-optical light, the re-emission in the infrared (IR) and via emission lines in the optical, IR and radio. In addition to this, bright events, such as active galactic nuclei (AGN) can significantly affect the observed SEDs of galaxies (see Conroy 2013 for a review on galaxy SEDs).

Truly multi-wavelength surveys, such as GAMA (Driver et al., 2009) in the local Universe and COSMOS (Scoville et al., 2007), CANDELS (Koekemoer et al., 2011) and DEVILS (Davies et al., 2018) in the high-redshift Universe, are becoming more common, and attempt to get a full picture of galaxy properties across the electromagnetic spectrum and cosmic time. This has allowed a full reconstruction of how the stellar mass, star formation rate (SFR), ISM and dust masses evolve with time for the overall population of galaxies (Santini et al., 2014; Scoville et al., 2016; Driver et al., 2018), the integrated SEDs (referred to as cosmic SED, CSEDs) of galaxies as a function of time (Andrews et al., 2017), the size-luminosity correlation as a function of wavelength in the local Universe (Lange et al., 2015), the IR-UV correlation as a function of redshift (Capak et al., 2015), among many others. The multi-wavelength nature of these surveys can also unveil the contribution from different galaxy populations to the cosmic SFR density of the Universe: at $z=0$ most star formation takes place in galaxies that are bright in the UV-to-optical, while at $z\gtrsim 1$ IR-bright galaxies tend to dominate (e.g. Casey et al. 2012; Magnelli et al. 2013; Madau & Dickinson 2014). These observations require cosmological galaxy formation simulations to be able to reliably predict SEDs of galaxies in as much of the electromagnetic spectrum as possible in order to offer a physical framework in which to interpret these observations, and to truly exploit their constraining power.

Multi-wavelength predictions covering from the FUV to the FIR have been challenging to produce because of the associated computational cost and uncertainties in the modelling process. In semi-analytic models (SAMs) of galaxy formation, a tool used to follow the formation and evolution of galaxies in DM halo merger trees from cosmological $N$-body simulations, this has been notoriously difficult. Early on Baugh et al. (2005), using GALFORM, noticed that there was significant tension arising when attempting to reproduce the FUV-to-near IR (NIR) and the FIR emission of galaxies simultaneously, and suggested that allowing for deviations from a universal initial stellar mass function (IMF) of stars in the case of starbursts helped solve the tension. This was done using a full radiative transfer (RT) approach in SAM galaxies, assuming a two-phase dust model in idealised geometries and employing the code GRASIL (Granato et al., 2000). Lacey et al. (2016) confirmed this conclusion in an updated version of GALFORM by adopting a more simplified method to predicting the FIR emission of galaxies. Cowley et al. (2019) showed that this tension also impacted the CSED and extra-galactic background light predictions.

Other SAMs, such as that of Somerville et al. (2012) have also attempted to predict the full FUV-to-FIR SEDs of galaxies. They used a different approach to Baugh et al., in that they used an attenuation model similar to that of Charlot & Fall (2000), with attenuation parameters varying with galaxy properties, and used observed dust templates to inform their model on how to re-emit the light in the IR. Somerville et al. (2012) scaled the optical depth with sensible galaxy properties, such as gas metallicity, gas mass and galaxy size, but without a theoretical motivation for their exact scaling. Despite this uncertainty, they found their model to provide a good match to the FUV-to-NIR emission of galaxies, but systematically underpredicted the emission at the FIR, finding a similar tension to that reported by Baugh et al. (2005).

In cosmological hydrodynamical simulations of galaxy formation the situation is not less different. Trayford et al. (2017) presented a full RT treatment of galaxies in the EAGLE simulations, which allowed the authors to produce FUV-to-FIR SEDs for all their galaxies. Camps et al. (2016); Baes et al. (2019); Cowley et al. (2019) showed that EAGLE was capable of reproducing the FUV-to-NIR emission of galaxies, but under-predicted the FIR emission, possibly suggesting the need for changes in their physical model by e.g. invoking a varying IMF.

A clear difficulty in providing predictions over the full FUV-to-FIR SED is how to simultaneously model the attenuation of stellar light and re-emission in the mid-to-far IR. To avoid this difficulty, many other SAMs and cosmological hydrodynamical simulations of galaxy formation limit themselves to modelling only the optical-to-NIR emission by using a slab or Charlot & Fall (2000)-like attenuation curves (see e.g. De Lucia & Blaizot 2007; Croton et al. 2016; Henriques et al. 2015; Yung et al. 2019 for examples from SAMs and Trayford et al. 2015; Nelson et al. 2018; Vogelsberger et al. 2019 for hydrodynamical simulations). Although the latter may be a pragmatic approach to tackle traditional galaxy surveys (e.g. SDSS, HST-based), future surveys are likely to move towards a more panchromatic view of galaxies, not only at $z\lesssim 2$ (e.g. GAMA and DEVILS, COSMOS, CANDELS, WAVES), but also at high redshift using the unprecedented combination of HST, JWST and ALMA.

Here, we use the recently introduced SAM of galaxy formation Shark (Lagos et al., 2018) in combination with RT results from the EAGLE simulations of Trayford et al. (2017) to produce a physically-motivated model for the attenuation of light in galaxies from the FUV to the NIR, and adopt an energy-conserving approach combined with observational IR templates (Dale et al., 2014) to re-emit the light in the mid-to-far IR. Our aim is to understand to what extent our state-of-the-art model can reproduce the observed FUV-to-FIR emission of galaxies and whether fine tuning and/or changes in the physical model (such as invoking a varying IMF) are required. Our approach is similar to Somerville et al. (2012) in that we start by adopting the Charlot & Fall (2000) parametric attenuation form, but we instead use the RT-predicted attenuation curves of EAGLE to inform Shark on how to scale the attenuation parameters with galaxy properties.

The advantage of using EAGLE to inform Shark, is that in EAGLE there is no need for assumptions about the geometry of the gas in galaxies and hence the derived attenuation parameters should not be biased by those assumptions (e.g. axi-symmetry, exponential radial profiles), which is a major risk in the case of RT applied to SAMs. Although using EAGLE allows us to relax typical assumptions made in SAMs, there are still important limitations. Most notably is the ISM model, which is sub-grid in simulations of coarse resolution such as EAGLE, directly impacting how “clumpy” the ISM of galaxies can be. Other sub-grid physical processes, such as stellar and active galactic nuclei (AGN) feedback, also impact the distribution of gas in galaxies, affecting the predicted attenuation. Hence, we ought to continue testing the validity of the attenuation model adopted here as simulations of higher resolution and improved ISM physics become available. Note that we do not attempt to tune to observations and instead combine the EAGLE RT results with Shark and, when necessary, adopt standard attenuation parameters widely adopted in the literature. The Shark model and SEDs presented here will be used to create panchromatic lightcones for the upcoming surveys DEVILS, WAVES, among others.

This paper is organized as follows. $\S$ 2 introduces Shark, describing the main physical processes included in the model, highlighting some key features and successes. We also describe how dust masses are computed. $\S$ 3 describes how we generate SEDs and the models we use for extinction and re-emission in the FIR. $\S$ 4 presents a comprehensive study of the galaxy LF from the FUV to the FIR, and from $z=0$ to $z=10$. We compare with available observations and analyse the physical drivers behind the predicted LF evolution. $\S$ 5 presents an analysis of galaxy number counts from the NUV-to-FIR, and the cosmic SED, how it is affected by extinction, compare with observations when available, and break down the total light budget into the contribution from different galaxy components. Finally, in $\S$ 6 we discuss the implications of our main findings, and the main successes and limitations of our work.

2 The Shark semi-analytic model

Shark, introduced by Lagos et al. (2018), is an open source, flexible and highly modular SAM111https://github.com/ICRAR/shark. The model includes all the physical processes that we think shape the formation and evolution of galaxies. These are (i) the collapse and merging of DM halos; (ii) the accretion of gas onto halos, which is modulated by the DM accretion rate; (iii) the shock heating and radiative cooling of gas inside DM halos, leading to the formation of galactic disks via conservation of specific angular momentum of the cooling gas; (iv) star formation in galaxy disks; (v) stellar feedback from the evolving stellar populations; (vi) chemical enrichment of stars and gas; (vii) the growth via gas accretion and merging of supermassive black holes; (viii) heating by AGN; (ix) photoionization of the intergalactic medium; (x) galaxy mergers driven by dynamical friction within common DM halos which can trigger starbursts and the formation and/or growth of spheroids; (xi) collapse of globally unstable disks that also lead to starbursts and the formation and/or growth of bulges. Shark adopts a universal Chabrier (2003) IMF. Lagos et al. (2018) include several different models for gas cooling, AGN, stellar and photo-ionisation feedback, and star formation. Here, we adopt the default Shark model (see models and parameters adopted in Lagos et al. 2018; their Table $2$).

An important assumption in Shark and any SAM is that galaxies can be described as a disk plus bulge at any time. The main distinction between these two components is their origin, while disks form stars from gas that is accreted onto the galaxy from the halo, bulges are built by stars that are accreted from satellite galaxies and starbursts that are driven by galaxy mergers or disk instabilities. Both disks and bulges in Shark form stars based on the surface density of molecular hydrogen, with the only difference being that in the latter the efficiency of conversion into stars is $10$ higher than for star formation in disks. In our default Shark model, we use the pressure relation of Blitz & Rosolowsky (2006) to estimate the radial breakdown between atomic and molecular gas. The higher H2-stars conversion efficiency in starbursts is found to be key to reproduce the cosmic star formation rate density (CSFRD) at $z\gtrsim 1.5$ in Shark (Lagos et al., 2018). As mentioned above, bulges can grow via disk instabilities, which happen when self-gravity dominates over centrifugal forces. This is evaluated by a global Toomre’s instability parameter (Ostriker & Peebles, 1973; Efstathiou et al., 1982),

[TABLE]

where $V_{\rm circ}$ is the maximum circular velocity, $r_{\rm disk}$ is the half-baryon mass disk radius and $M_{\rm disk}$ is the total baryon disk mass. Here baryon corresponds to gas plus stars. The numerical factor $1.68$ converts the disk half-baryon mass radius into a scalelength, assuming an exponential profile. If $\epsilon<\epsilon_{\rm disk}$ the disk is considered to be unstable. In the default Shark model used here, $\epsilon_{\rm disk}=0.8$. Simple theoretical arguments suggest $\epsilon_{\rm disk}\approx 1$ (Efstathiou et al., 1982). However, because the process of bar creation and thickening of the disk can be a very complex phenomenon (Bournaud et al., 2011) that can easily lead to the gas and stars not having the same $\epsilon$ parameter (Romeo & Wiegert, 2011; Romeo & Mogotsi, 2018), in Shark we treat $\epsilon_{\rm disk}$ as a free parameter. Note that many other SAMs do not include the effect of disk instabilities (e.g. Henriques et al. 2015; Xie et al. 2017), though Fanidakis et al. (2012) and Griffin et al. (2018), using the GALFORM SAM (Cole et al., 2000; Lacey et al., 2016), argue that disk instabilities are a key physical processes required to obtain a realistic population of QSOs throughout cosmic time.

In Shark, we numerically solve the differential equations (DEs) of mass, metals and angular momentum exchange between the different baryon reservoirs (see Eqs. $49$-$64$ in Lagos et al. 2018), only setting an accuracy to which these equations are solved. The baryon reservoirs in the model are: gas outside halos, hot and cold gas inside halos but outside galaxies, ionised/atomic/molecular gas and stars in disks and bulges in galaxies, and super-massive black holes. This approach makes our model less sensitive to the time-stepping of the $N$-body simulation compared to other models, and also means that the star-formation histories (SFH) of galaxies can have as complex shape as required to solve the DEs.

The model parameters of our default Shark model were tuned to the $z=0,\,1,\,2$ stellar mass functions (SMFs), the $z=0$ the black hole-bulge mass relation and the mass-size relations. The model also reproduces very well observational results that are independent from those used for the tuning, such as the total neutral, atomic and molecular hydrogen-stellar mass scaling relations at z=0, the cosmic star formation rate (SFR) density evolution at $z\approx 0-4$, the cosmic density evolution of the atomic and molecular hydrogen at $z\lesssim 2$ or higher in the case of the latter, the mass-metallicity relations for the gas and stars, the contribution to the stellar mass by bulges and the SFR-stellar mass relation in the local Universe (see Lagos et al. 2018 for more details). In addition, Davies et al. (2019) show that Shark also reproduces the scatter around the main sequence of star formation in the SFR-stellar mass plane, Chauhan et al. (2019) show that Shark reproduces very well the HI mass and velocity width of galaxies observed in the ALFALFA survey and Amarantidis et al. (2019) show that the AGN LFs agree well with observations in the X-rays and radio wavelengths. These represent true successes of the model as none of these observations were used in the processes of tuning the free parameters.

With the aim of building the SEDs of galaxies, Shark produces an output file star*-formation-histories, which contain the amount of stars that formed and the metallicity with which they formed throughout all the epochs sampled by the snapshots of the simulation until the point in which the output is being written. This is done separately for stars that end up in the disk and the bulge by the time of the output. Bulges are separated into stars built up by galaxy mergers and by disk instabilities. If a galaxy has a bulge that was built up by these two processes, then both arrays will have non-zero inputs. This information is then used by Viperfish (described in $\S$ 3) to create the SEDs and consequently calculate the galaxies’ emission in a large range of bands going from the far-UV (FUV) to the far-IR (FIR). Because we solve the DEs numerically, the arrays in star-formation-*histories show the average SFR and metallicity from which stars formed in the $200$ snapshots of the $N$-body simulation (see details below).

2.1 The surfs simulations

The results presented in Lagos et al. (2018) were produced using the surfs suite of N-body, DM-only simulations (Elahi et al., 2018b), most of which have cubic volumes of $210\,\rm cMpc/h$ on a side, and span a range in particle number, currently up to $8.5$ billion particles using a $\Lambda$CDM Planck Collaboration et al. (2016) cosmology. These correspond to a total matter, baryon and $\Lambda$ densities of $\Omega_{\rm m}=0.3121$, $\Omega_{\rm b}=0.0491$ and $\Omega_{\ L}=0.6751$, respectively, with a Hubble parameter of $h=100\,\rm Mpc\,km\,s^{-1}$ with $h=0.6751$, scalar spectral index of $n_{\rm s}=0.9653$ and a power spectrum normalization of $\sigma_{\rm 8}=0.8150$. All simulations were run with a memory lean version of the gadget2 code on the Magnus supercomputer at the Pawsey Supercomputing Centre. In this paper, we use the L210N1536 simulation, which has a cosmological volume of $(210\,\rm cMpc/h)^{3}$, $1536^{3}$ DM particles with a mass of $2.21\times 10^{8}\,\rm h^{-1}\,M_{\odot}$ and a softening length of $4.5\,\rm h^{-1}\,ckpc$. Here, cMpc and ckpc denote comoving Mpc and kpc, respectively. surfs produces $200$ snapshots for each simulation, typically having a time span between snapshots in the range of $\approx 6-80$ Myr.

Merger trees and halo catalogs, which are the basis for Shark (and generally any SAM), were constructed using the phase-space finder VELOCIraptor222https://github.com/icrar/VELOCIraptor-STF/ (Elahi et al., 2019a; Cañas et al., 2019) and the halo merger tree code TreeFrog333https://github.com/pelahi/TreeFrog, developed to work on VELOCIraptor (Elahi et al., 2019b). Poulton et al. (2018) show that TreeFrog+VELOCIraptor lead to very well behaved merger trees, with orbits that are well reconstructed. Elahi et al. (2018a) also show that these orbits reproduce the velocity dispersion vs. halo mass inferred in observations. Cañas et al. (2019) show that the same code can be applied to hydrodynamical simulations to identify galaxies and that the performance of VELOCIraptor is superior to space-finders, even in complex merger cases. We refer to Lagos et al. (2018) for more details on how the merger trees and halo catalogs are constructed for Shark, and to Elahi et al. (2019a, b); Cañas et al. (2019); Poulton et al. (2018) for more details on the VELOCIraptor and TreeFrog software.

2.2 Calculation of dust masses

In this paper we consider three models to compute the dust mass from the mass in metals and the gas metallicity:

• •

A constant dust-to-metals mass ratio, set to the Milky-Way value $M_{\rm dust}=0.33\,M_{\rm Z}$ (Rémy-Ruyer et al., 2014) (referred to as $f_{\rm dust}$-const).

• •

The best fit $M_{\rm dust}/M_{\rm Z}-Z_{\rm gas}$ relation of Rémy-Ruyer et al. (2014) (see thick dotted line in Fig 1; see Table $1$ in Rémy-Ruyer et al. 2014; referred to as RR14).

• •

The case in which a steeper relation is assumed with a break at higher gas metallicities, following the thin, black dotted line of Fig 1. This is motivated by the slope of the best fit $M_{\rm dust}/M_{\rm Z}-Z_{\rm gas}$ relation in Rémy-Ruyer et al. (2014) being quite significant and the recent observations of De Vis et al. (2019) seemingly favouring a break in the dust-to-metal ratio at higher gas metallicities (referred to as RR14-steep).

The three different options above are shown in Fig 1, and are expected to make a difference only in galaxies with $Z_{\rm gas}/Z_{\odot}<0.25$. This means that in the local Universe, only dwarf galaxies are expected to deviate from the constant dust-to-metal mass ratio significantly, and high redshift galaxies, as most of them have lower metallicities, deviating from $M_{\rm dust}=0.33\,M_{\rm Z}$.

Below, we describe how we compute $\Sigma_{\rm dust}$ for disks and bulges in Shark.

• •

Disks. We compute an average $\Sigma_{\rm dust}$ for the disk as:

[TABLE]

where $M_{\rm dust,disk}$ is the dust mass in the disk, $r_{\rm 50,d}$ is the half-gas mass radius of the disk and in a projected image represents the major axis, and $l_{\rm 50}$ is the projected minor axis, which is calculated as $l_{\rm 50}={\rm sin}(i)*(r_{\rm 50,d}-r_{\rm 50,d}/7.3)+r_{\rm 50,d}/7.3$, where $i$ is the inclination. The latter is $=r_{\rm 50,d}$ if the galaxy is perfectly face-on, and $=r_{\rm 50,d}/7.3$ if the galaxy is perfectly edge-on. The value $7.3$ comes from the scaleheight to scalelength observed relation in local galaxy disks Kregel et al. (2002). The inclination of a galaxy comes from the host subhalo angular momentum vector, or in the case of orphan galaxies, it is randomly chosen (see Chauhan et al. 2019 for details).

• •

Bulges. We assume bulges to be spherically symmetric and hence the inclination is unimportant. We then compute the bulge dust surface density:

[TABLE]

where $M_{\rm dust,bulge}$ is the dust mass in the bulge, and $r_{\rm 50,b}$ is the half-gas mass radius of the bulge.

Fig. 2 shows the resulting dust surface density evolution for the disks and bulges, computed as in Eqs. 2 and 3, respectively, of Shark galaxies at $z=0$ to $z=8$, for the RR14-steep scaling. Bulges display a monotonic evolution, with $\Sigma_{\rm dust}$ increasing with increasing redshift at fixed mass over the whole redshift range analysed here. This is due to a combination of the gas surface density evolution, in which high-z galaxies have higher $\Sigma_{\rm gas}$, and the fact that for bulges there is little evolution of the stellar mass-gas-metallicity relation.

Galaxy disks on the other hand, display a more complex behavior. At $M_{\rm star}\gtrsim 10^{9.5}\,\rm M_{\odot}$, galaxies show a $\Sigma_{\rm dust}$ that increases from $z=0$ to $z\approx 2$, followed by a decrease towards higher redshift, at fixed stellar mass. At $10^{8}\,\rm M_{\odot}\lesssim M_{\star}\lesssim 10^{9.5}\,\rm M_{\odot}$, this reversal happens at higher redshift, $z\approx 4$. At lower stellar masses we see that the reversal moves to even higher redshift. However, those masses are below what we would consider as “resolved” in our simulation box. Lagos et al. (2018) showed that the box used here is reliable down to $M_{\star}\approx 10^{8}\,\rm M_{\odot}$, but below that the number density of galaxies artificially drops, deviating from the values obtained from a higher resolution box of the same cosmology and initial conditions. The evolution of $\Sigma_{\rm dust}$ for disks is driven by the competing effects of the gas metallicity and $\Sigma_{\rm gas}$ evolution. At fixed stellar mass, Shark galaxies exhibit a strong $Z_{\rm gas}$ evolution, with galaxies at $z=3$ being $0.6$ dex metal poorer than galaxies at $z=0$ at fixed stellar mass. However, in the same redshift range, $z=3$ galaxies have a $\Sigma_{\rm gas}$ that is $\approx 1.2$ dex larger than the $z=0$ counterparts of the same stellar mass. As a result, the evolution seen in $\Sigma_{\rm dust}$ is more modest than that obtained for $\Sigma_{\rm gas}$ and the reversal displayed is due to the metallicity evolution overcoming the increase in $\Sigma_{\rm gas}$.

3 Lighting Shark galaxies through Viperfish

To generate SEDs for Shark, two packages have been developed: ProSpect444https://github.com/asgr/ProSpect and for an interactive ProSpect web tool see http://prospect.icrar.org/, which is recommended as an education tool. and Viperfish555https://github.com/asgr/Viperfish. ProSpect (Robotham et al. in prep) is a low-level package that combines the popular GALEXev stellar synthesis libraries Bruzual & Charlot (2003) (BC03 from hereafter) and/or EMILES (Vazdekis et al., 2016) with a multi-component dust attenuation model (Charlot & Fall, 2000) and dust re-emission (Dale et al., 2014). On top of this sits Viperfish, which allows for simple extraction of Shark SFHs, metallicity histories (ZFH), and generation of the desired SED through target filters.

ProSpect is designed in a pragmatic manner that allows for user-side flexibility in controlling the key components that affect the galaxy SED produced. Many of the design decisions were influenced by successful spectral fitting codes (e.g. MAGPHYS, da Cunha et al. 2008, and Cigale; Noll et al. 2009) with the emphasis here on a code that works in a fully generative mode with the types of outputs available from SAMs. Other differences lie in the specific choice of dust modelling (in particular the re-emission templates) and the manner in which SFHs and ZFHs are incorporated (highly flexibly).

For the production of galaxy SEDs, the decision was made early on to focus efforts on the BC03 stellar population (SP) libraries using a Chabrier (2003) IMF since these are well understood in the community, have a broad spectral range that makes them useful for current and next generation multi-band surveys and are the default in Shark. ProSpect can accept almost any functional form for the SFH or ZFH, which includes non-parametric, parametric or discontinuous specifications (the latter being most like the type produced in a modern SAM). The functional SFH or ZFH can in practice be arbitrarily complex, with internal interpolation schemes used to map the provided form onto the discrete library of temporal evolution available. For the ZFH, the metallicities are interpolated in log-space, producing a few tenths of a mag uncertainty at worst within the range available ($0.0001\leq Z\leq 0.05$). If the time-steps in which the SFH and ZFH are stored are too coarse, this interpolation may lead to large uncertainties in the predicted emission, particularly in the UV. Fortunately, the time-steps of our surfs simulations are sufficiently fine so that the UV emission is accurately predicted. In the worst case scenario of an extreme recent starburst, the UV would still be converged to better than $30$%, but in more common cases we expect an accuracy of $5$% or better.

The generative nature of ProSpect means it can be used in a number of ways: either to fit real data using Bayesian modelling via optimisation of Markov-Chain Monte-Carlo (MCMC; see Bellstedt et al. in prep. and Davies et al. in prep.), or in a purely generative mode given a SFH and ZFH evolution of, e.g., a simulated galaxy. For producing lightcones with SEDs from SAMs, this generative mode is obviously of most interest. However, some sensible assumptions must be made regarding light attenuation due to dust, and its re-emission at longer wavelengths. How to do this in a fully physical sense, given the limited range of knowledge we have about any single SAM galaxy, is a matter of ongoing research, but for the current purposes of Shark SED generation we settle on a deliberately simplified fiducial model of dust processing.

Firstly, the dust is attenuated by the dust model of Charlot & Fall (2000), in which the dust is assumed to be in a two-phase medium (birth clouds, BC, and diffuse ISM) in both the disk and the bulge (in which starbursts take place). Two different optical depths at $5500$Å are assumed for these phases, $\hat{\tau}_{\rm BC}$ and $\hat{\tau}_{\rm ISM}$, respectively. The absorption curves for the BCs and diffuse ISM are then defined as:

[TABLE]

The values we adopt as Charlot & Fall (2000) default are $\hat{\tau}_{\rm BC}=1$, $\hat{\tau}_{\rm ISM}=0.3$, $\eta_{\rm BS}={\eta_{\rm ISM}}=-0.7$ (suggested to be within a “reasonable” range in that paper). Stellar populations younger than $10$ Myr are in birth clouds, and hence their light is affected by both the optical depth of Eq. 5, while older stars which are in the diffuse medium are attenuated by Eq. 4.

With this model, light generated at different ages is attenuated differently, giving a natural means to simulate the effect of BC attenuation for younger stars. This absorbed light must then be re-emitted in a sensible fashion at longer wavelengths. For this process we adopt the Dale et al. (2014) FIR dust templates, with an assumption of no significant AGN emission, and an assumed dust radiation field of $\alpha_{\rm SF}=3$ for the diffuse ISM and $\alpha_{\rm SF}=1$ for the birth clouds. Since this re-emission process only makes use of the absorbed luminosity in the UV-NIR, the scaling is chosen to ensures energy balance. The $\alpha_{\rm SF}$ exponents represent the local interstellar radiation field the dust is exposed to, $0.3<U<10^{5}$, with $U=1$ being the local interstellar radiation field of the solar neighborhood. A power-law combination of local curves mimics the global dust emission, with a fraction ${\rm d}M_{\rm dust}$ of dust mass being heated by $U^{-\alpha_{\rm SF}}{\rm d}U$. The values adopted here for the screen and birth cloud components roughly correspond to effective dust temperatures of $20-25$ K and $50-60$ K, respectively. Note that emission from AGN can be included when using ProSpect to fit the SEDs of observed galaxies; however, we do not use it in Shark as it requires significant additional modelling to scale the AGN SED templates with meaningful AGN properties. We leave this for future work.

Once the full generative spectrum has been created (by adding the attenuated stellar light and the dust emission together), we redshift to the observed frame using the full spectral resolution available. Finally, we pass the spectrum through a chosen number of available filters that span the FUV to FIR, giving our final reduced outputs. Storing the spectral information of all galaxies is impractical, so care must be taken that all filters of interest are specified at this stage. Only a subset of these filters are discussed in this work, and user defined filters can be added easily if required. We warn the reader that in this work we do not include nebular emission lines, which can make an impact on narrow bands. Hence, in this work we focus solely on broad band emission.

The highest level code Viperfish allows for a very simple interface between the HDF5 outputs created by Shark and ProSpect. It effectively reduces a few hundred lines of R code to a single call with the path to the relevant HDF5 file. This makes it trivial to post-process any Shark outputs at any time (it does not need to be run in parallel), and it is designed to scale naturally with the computing resources available, e.g. it can use multiple cores.

3.1 Optical depth and reddening calculation of Shark galaxies

3.1.1 Attenuation due to the diffuse ISM

Trayford et al. (2017) used the RT code SKIRT to compute the attenuation curve for each galaxy in the EAGLE hydrodynamical simulation suite. From these curves, Trayford et al. (2019) found that they can be parametrized using the Charlot & Fall (2000) model, with values for $\tau_{\rm ISM}$ and $\eta_{\rm ISM}$ varying with the dust column density in the line of sight (hence, considering the effects of inclination). Trayford et al. (2019) in fact find that such parametrization is independent of redshift. Hence the redshift evolution obtained for the average optical depth and power-law index of Eq. 4 of galaxies is due to their dust surface density evolving.

Trayford et al. computed the median and $1\sigma$ scatter relationship between $\tau_{\rm ISM}$, $\eta_{\rm ISM}$ and $\Sigma_{\rm dust}$, from which we sample. In Shark, we use each galaxy’s dust surface density, $\Sigma_{\rm dust}$, to compute $\tau_{\rm ISM}$ and $\eta_{\rm ISM}$, and perturb the values by sampling from a gaussian distribution with width $\sigma$, where $\sigma$ is the $16^{\rm th}-84^{\rm th}$ percentile ranges predicted by Trayford et al. (2019). We compute $\Sigma_{\rm dust}$ for disks and bulges following Eqs. 2 and 3.

3.1.2 Attenuation due to birth clouds

For the birth clouds we follow Lacey et al. (2016), who assume the birth cloud optical depth to scale with the gas metallicity and gas surface density of the cloud, but modify it to use the dust surface density of clouds rather than the metal surface density,

[TABLE]

$f_{\rm dust}=M_{\rm dust}/M_{\rm Z}$ is the dust-to-metal mass ratio, $\tau_{\rm BC,0}=1$, $\Sigma_{\rm MW,cl}=85\,\rm M_{\odot}\,pc^{-2}$, $Z_{\odot}=0.0189$, and $f_{\rm dust,MW}=0.33$, so that in typical spiral galaxies $\tau_{\rm BC}\approx\tau_{\rm BC,0}$ as determined by Charlot & Fall (2000) and Kreckel et al. (2013). We compute the cloud surface density as $\Sigma_{\rm gas,cl}={\rm max}[\Sigma_{\rm MW,cl},\Sigma_{\rm gas}]$, with $\Sigma_{\rm gas}$ being the diffuse medium gas surface density, which is calculated as Eqs. 2 and 3, but using the gas masses of the disk and bulge, respectively. The reasoning behind this is that in the local group, galaxies ranging from metal-poor dwarfs to molecule-rich spirals seem to have giant molecular clouds (GMCs) with a constant gas surface density close to the value $\Sigma_{\rm MW,cl}$, which is surprisingly independent of galactic environment (see e.g. Blitz et al. 2007; Bolatto et al. 2008, and Krumholz 2014 for a review). However, as the ambient ISM pressure increases, the GMC surface density must increase in order to maintain pressure balance with the surrounding ISM. Hence, it follows that $\Sigma_{\rm gas,cl}\approx\Sigma_{\rm gas}$ in those extreme environments (Krumholz et al., 2009), which are expected to be more common at high redshift. We also impose the physical limit of $\tau_{\rm BC}\geq\tau_{\rm ISM}$.

For birth clouds we do not have a well informed choice for $\eta_{\rm BC}$, as we do for the diffuse ISM, and hence we adopt the default Charlot & Fall (2000) $\eta=-0.7$. Some models in the literature assume a more negative value of $-1.3$ (e.g. da Cunha et al. 2008; Wild et al. 2011) due to the expected shell-like geometry of BCs. We find, however, that the use of a steeper $\eta_{\rm BC}$ does not affect our results in any significant manner.

3.1.3 Summary of Attenuation models

Table 1 shows all the attenuation models used here: (1) the simplest assumption, which corresponds to fixed Charlot & Fall (2000) parameters (which are therefore constant and do not depend on galaxy properties or inclination; referred to as CF00); (2) the EAGLE attenuation parametrization of the Charlot & Fall (2000) parameters, assuming a constant fraction of the metals are locked in dust (referred to as EAGLE-$\tau$ $f_{\rm dust}$ const); (3) as model (2) but assuming the empirical relation of Rémy-Ruyer et al. (2014) between $M_{\rm dust}-M_{\rm Z}-Z_{\rm gas}$ (see thick dashed line in Fig. 1; referred to as EAGLE-$\tau$ RR14); (4) as model (3) but using a steeper dependence of $M_{\rm dust}/M_{\rm Z}$ on $Z_{\rm gas}$ within the errors of the best fit relation in Rémy-Ruyer et al. (2014) (see thin, black dotted line in Fig. 1; referred to as EAGLE-$\tau$ RR14-steep). Model (3) is our default model throughout the paper but we make it explicit in every figure caption which the model is shown.

3.1.4 Stellar mass dependence and redshift evolution of the optical depth of Shark galaxies

Fig. 3 shows the effective $V$-band optical depth, $\tau$, as a function of stellar mass at several redshifts for the disks and bulges of Shark galaxies. Here we adopt the attenuation model EAGLE-$\tau$ RR14 (see Table 1).

For the diffuse ISM, we obtain a steep increase of $\tau$ of galaxy’s disks with stellar mass at $M_{\star}>10^{10}\,\rm M_{\odot}$ at $z=0$, below which $\tau\rightarrow 0$. This stellar mass threshold moves to lower stellar masses as redshift increases, up to $z\approx 1$ for disks and at all redshifts for bulges. In the latter, $\tau\lesssim 0.5$ for all galaxies at $z=0$ due to the gas fractions of bulges being very small. This changes at $z\gtrsim 1.5$ due to bulges hosting large gas reservoirs and undergoing starbursts. Although galaxies at high redshift are more metal poor, their gas surface density is increasing rapidly, causing the redshift evolution seen in Shark galaxies. For the BCs, we obtain a relatively sharp transition from small to large extinctions at $M_{\star}\approx 10^{10}\,\rm M_{\odot}$ in disks, which is dictated by the gas metallicity, and at the high mass galaxies, by the average gas surface density. This transition moves to lower stellar mass for bulges (which tend to be more compact and more metal rich than disks), and to progressively lower masses as the redshift increases, mostly driven by the evolution of the bulge gas surface density. Adopting instead the attenuation models EAGLE-$\tau$ RR14 or EAGLE-$\tau$ $f_{\rm dust}$ const results in a shift of the $y$-axis values in both Figs. 3 and 4 to higher $\tau$ values, overall producing more attenuation (not shown here).

3.2 Example SEDs and SFHs

Fig. 5 shows examples of SFHs of $10$ randomly selected Shark galaxies at $z=0$ that have stellar masses $>10^{9}\,\rm M_{\odot}$ and stellar-mass weighted ages at $\pm 0.3$ Gyr around the values labelled in each panel, which span from $11$ Gyr to $3$ Gyr. The SFHs of Shark galaxies look anything but the idealized exponentially decay or composite instantaneous-burst plus exponential decay, which are typically assumed in observations when performing SED fitting (Mitchell et al., 2013; da Cunha et al., 2008). Pacifici et al. (2012) used SFHs and ZFHs from SAMs as inputs for the SED fitting of observed galaxies. This makes an important difference in the recovered stellar mass and SFR of up to a factor of $0.6$ dex (see e.g. Pacifici et al. 2015). This shows that using complex SFHs is important in the recovery of galaxy parameters.

Many Shark galaxies experience early starbursts seen as short-lived peaks in the SFH (quite common at look-back times $\gtrsim 10$ Gyr). The latter are more common in galaxies that have older stellar populations by $z=0$ than younger ones. At look-back times $\lesssim 6$ Gyr, starbursts are much less common, mostly seen in galaxies that by $z=0$ are very young. Also note that old galaxies tend to show sharp cut-offs in their SFH associated to stripping of their hot gas as they become satellite galaxies. On the contrary, galaxies that are on average young by $z=0$, tend to have very extended SFHs, that in some cases continue to rise to $z=0$. Most central galaxies that by $z=0$ are old tend to have SFHs that drop towards $z=0$, but less sharply than for satellites (see for example solid lines vs. dashed lines in Fig. 5).

Fig. 6 shows the broadband SED in $27$ bands for $3$ randomly selected Shark galaxies of different stellar ages. The SFHs of these galaxies are shown in the insets in each panel. Both the intrinsic emission and after dust attenuation and re-radiation are shown. As expected, young galaxies tend to have much more significant emission in the UV, which suffers from large extinction. Galaxies with ages $\gtrsim 11$ Gyr have very little intrinsic emission in the UV and little gas content, both of which result in a small extinction. We show in Fig. 7 the SEDs of three starburst galaxies at $z=2$ in the same $27$ bands of Fig. 6. These galaxies have widely different star formation histories, with one of them having significant star formation over the last $300\,\rm Myr$ but little before that. These galaxies differ significantly from the $z=0$ examples in that most of their emission happens at the FIR, and represent nice examples of sub-millimeter galaxies (SMGs) in Shark.

4 Galaxy emission and the effects of dust extinction on the galaxy LF

In this section we analyse the Shark predictions for the FUV-to-FIR emission of galaxies at $z=0$ and how this is affected by our new attenuation models. We specially focus on the properties of the different structural components of galaxies and the connection to their stellar populations and epochs.

4.1 The z=0 UV and FIR luminosity functions

Fig. 8 shows the $z=0$ GALEX FUV and NUV luminosity functions (LFs) predicted by Shark before and after dust attenuation is applied. We show four attenuation models corresponding to those in Table 1. The top panels show the total LFs. We also show the observations of Driver et al. (2012).

Galaxies at $z=0$ emit several orders of magnitude more UV emission than is observed (thin, solid lines in Fig. 8), meaning that extinction must play a very important role, particularly beyond the break of the LF, $L^{*}$. Adopting the CF00 extinction parameters leads to FUV and NUV LFs that are too shallow at the faint end ($>-17.5$ AB mags). The attenuation models based on the EAGLE RT massively improve the predicted faint end of the LFs The attenuation models EAGLE-$\tau$ $f_{\rm const}$ and EAGLE-$\tau$ RR14 produce almost identical LFs, due to most galaxies contributing to the UV LFs having $Z_{\rm gas}/Z_{\odot}>0.25$, which is the gas metallicity threshold above which galaxies converge to a constant dust-to-metal mass ratio (see Fig. 1). The extinction model EAGLE-$\tau$ RR14-steep, on the other hand, predicts a slightly brighter break of the LF (by $\approx 0.2$ mag). This difference is due to galaxies in this variant deviating from the constant dust-to-metal ratio at $Z_{\rm gas}\approx 0.7Z_{\odot}$ (see Fig. 1). Note that all the extinction models miss the sharp bright-end of the UV LFs, which indicate that Shark galaxies are slightly too star-forming and/or the attenuation for the brightest UV galaxies is too small. The obvious improvement obtained where going from the default CF00 to the EAGLE-like extinction models justifies the need for the added complexity, and nicely confirms that our RT-motivated extinction models allow Shark to predict more realistic UV LFs. The latter becomes even clearer at higher redshifts ($\S$ 4.3).

The middle and bottom panels of Fig. 8 show the contribution from disks and bulges to the FUV and NUV LFs at $z=0$, respectively. Bulges are only important at the very bright end; these galaxies correspond to the few rare local starburst. Note that the attenuation models based on the EAGLE-RT results produce virtually the same bulge UV LF, which is due to bulges having gas metallicities typically above $0.7\,Z_{\odot}$. This means that bulges have the same dust-to-metal ratio in the three EAGLE-RT variants of Table 1. This is not the case for disks, which is why the three EAGLE-RT model variants produce different UV LFs. Because disks dominate over the whole magnitude range, we end up with visible differences in the total UV LFs.

The better match to the faint end of the UV LFs by the EAGLE-$\tau$ attenuation models is the dependence of the gas surface density on stellar mass (which produces a differential optical depth): $z=0$ Shark galaxies of $M_{\star}\approx 10^{8}\,\rm M_{\odot}$ have $\Sigma_{\rm gas}\approx 10^{6.5}\,\rm M_{\odot}/kpc^{2}$, while $M_{\star}\approx 10^{10}\,\rm M_{\odot}$ galaxies have $\Sigma_{\rm gas}\approx 10^{7.3}\,\rm M_{\odot}/kpc^{2}$.

The changes seen in the UV LF are expected to be also seen in FIR, as the light that is extincted by the dust is then re-radiated in the FIR. This is shown in Fig. 9 for the same $4$ attenuation models of Table 1, but here we only show the total LF as we later analyse the contribution from disks and bulges. Significant differences are seen at the faint end of the FIR LFs of up to $\approx 0.5$ dex in number density, but that regime unconstrained by observations. All models, however, predict a very similar bright-end, which agree very well with the observed LFs. We remind the reader that here we assume two effective dust temperatures for the diffuse ISM and BCs to re-emit the extincted light in the FIR when using the Dale et al. (2014) templates. The values we adopt are typical of the local Universe and hence the agreement with the observations is not necessarily surprising.

4.2 The $z=0$ FUV-to-FIR luminosity functions

Fig. 10 shows the UV and optical luminosity functions at $z=0$ compared to the measurements of Driver et al. (2012) using the Galaxy and Mass Assembly(GAMA) survey. The thin lines show the intrinsic emission, while the thick lines show the emission after dust extinction and reprocessing. As we discussed in $\S$ 4.1, the effect of the latter is very important in the UV bands, shifting the luminosity function by up to $2$ magnitudes at the bright-end and in the FUV band. The effect becomes a lot weaker in the optical. For example, in the $r$-band the effect is only $\approx 0.3$ mag.

The observations of Driver et al. (2012) correspond to the observed luminosity functions and they should be compared to the thick lines. The agreement with the observations is remarkable across all the bands, considering that we do not use this information to tune the free parameters of the model. The latter is less obvious at the near-IR bands, as this luminosity correlates strongly with stellar mass, and as explained in $\S$ 2, the $z=0$ SMF was used to tune the parameters. Thus, it is not necessarily surprising that the $z$-band LF agrees well with the observations.

As discussed in $\S$ 4.1, Shark tends to produce slightly too many UV bright galaxies; $\approx 0.5-0.7$ dex more galaxies than Driver et al. (2012) at a FUV $-19.3$ and NUV $-20$ magnitudes, due to the contribution of starbursts in Shark. This is seen in the dashed and dot-dashed lines in Fig. 10, which show the LFs of the bulges that formed predominantly by galaxy mergers and by disk instabilities in Shark, respectively. Both mechanisms of bulge formation contribute similarly to the number density of bright UV galaxies. Although this changes significantly as we move towards redder bands. In the $r$ to $z$ bands, bulges built by disk instabilities make a similar contribution as disks at the bright-end, which is much smaller that that of bulges built by galaxy mergers.

Stars in the disks of galaxies always dominate the faint-end of the luminosity functions, but their contribution beyond the break in the LF is a strong function of wavelength. The bluer the band, the higher the contribution from disks at the bright-end. In the extreme cases of the FUV and NUV luminosity functions, disks dominate the number density over all but the brightest luminosity bin, while in the $u$- and $g$-bands, they contribute about half of the luminosity above $L^{*}$. This contribution becomes negligible in the $z$-band, where the bright-end beyond $-21.5$ mag is primarily tracking the bulge content of galaxies. We later show that this trend reverses for the mid- and FIR bands at low redshifts (fig. 12).

Fig. 11 shows the $z=0$ luminosity functions for the $4$ UKIDSS bands, Y, J, H and K, an the IRAC $3.6\mu$m, $4.5\mu$m, $5.8\mu$m and $8\mu$m of Shark galaxies, compared to Driver et al. (2012) and Dai et al. (2009). The agreement between the model and the observations in the UKIDSS and IRAC $3.6\mu$m, $4.5\mu$m bands is excellent, except in the brightest luminosity bin. Again, this is not surprising as Shark is tuned to fit the SMF at $z=0$. The overabundance of very bright galaxies is similar to the conclusion of Lagos et al. (2018) that the SMF has a high-mass end slope a bit too shallow compared to the observations, leading to slightly too many galaxies with stellar masses $\approx 10^{12}\,\rm M_{\odot}$, though still within the observational uncertainties. Note that here we see a continuation of the trend of the contribution from disks at the bright-end decreasing as the wavelength becomes longer. At the K-band, disks have a negligible contribution over the whole magnitude range above $L*$. The reasonable agreement at the IRAC $5.8\mu$m and $8\mu$m bands is more surprising and shows that our attenuation plus dust-remission models have a realistic effect on the UV light and re-emission at the mid IR. However, Shark does not reproduce perfectly the IRAC $5.8\mu$m and $8\mu$m LFs, with most tension seen at the faint end, and the bright end in the $5.8\mu$m band. These bands are particularly difficult as most of the emission comes from unidentified infrared emission (UIE), which is a ubiquitous component of the IR emission in galaxies and typically associated to polycyclic aromatic hydrocarbons (Li & Draine, 2012).

The LF of bulges built by disk-instabilities peaks below $L*$, but with the peak moving to brighter luminosities relative to $L^{*}$ as the wavelength shortens. This agrees with the overall picture of the stellar mass budget build up described in Lagos et al. (2018), who showed that the stellar mass contribution from bulges built via disk instabilities peaks at stellar masses of $10^{10.3}-10^{10.8}\,\rm M_{\odot}$. Those galaxies contribute little to the UV luminosity functions, as $\approx 30-40$% of them are passive (i.e. specific SFRs $>10$ times below the main sequence of star formation), while their contribution increases in the NIR bands as their stellar mass is large. The bottom panels of Fig. 11 show the comparison with LFs measured in the IRAC bands at $z\approx 0$. The IRAC $3.6\mu$m and $4.5\mu$m, behave similarly to the UKIDDS bands, but the $5.8\mu$m band starts to show an increase in the contribution from disk emission, and the LF starts to be dominated by the re-emission of light by dust rather than the intrinsic stellar light. By the IRAC band $8\mu$m, disks are back to contributing most of the light, and to dominate even above $L^{*}$.

Fig 12 shows the $z=0$ luminosity functions in the $160\mu$m, $250\mu$m, $350\mu$m, $500\mu$m bands of the Herschel Space Observatory (Pilbratt et al., 2010), and the James Clerk Maxwell Telescospe (JCMT) $850\mu$m band. We show observational measurements as symbols. Some of these LFs (e.g. those of Marchetti et al. 2016) correspond to LFs measured in very wide redshift ranges ($z<0.5$); hence, we include the $z=0.25$ LF to show how much evolution is expected in that redshift window. Disks are the primary contributor over the whole magnitude range in the FIR bands, except in the brightest two bins, where starbursts either driven by galaxy mergers or disk instabilities, are significant. This is because at these wavelengths the re-emission of the UV light that was absorbed due to dust starts to become the most dominant source of light (see difference between the thin and thick lines in the bottom-right panel of Fig 11).

In the Herschel bands, we see that Shark’s predictions agree well with the observations within the systematic uncertainties of the data. At the $850\mu$m, the model produces a bright end that is slightly too bright, but we will see in $\S$ 5 that the total emission at this band agrees quite well with the observations, possibly indicating systematic effects are important.

To the knowledge of the authors, the agreement of Shark with the observed LFs in such a broad wavelength coverage is unprecedented and a success of the overall modelling included in Shark+ProSpect. This implies that galaxies have roughly correct SFRs, gas content and gas metallicities (which were shown in Lagos et al. 2018), as well as sizes, which together provide realistic dust surface densities. We also remind the reader that the adopted empirical scalings (e.g. the dust-to-metal ratio vs. gas metallicity of Rémy-Ruyer et al. 2014) or theory-inspired relations (e.g. the attenuation parameters of Trayford et al. (2019) are not tuned to get the LFs correct. Instead, quite naturally they allow Shark to provide realistic multi-wavelength properties of galaxies.

4.3 Redshift evolution of the UV and K-band LFs

We now focus on the evolution of the galaxy LF in two broadly studied bands: the rest-frame K- and FUV bands. Fig. 13 shows the $K$-band LF from $z=0.5$ to $z=3$ in Shark using the four extinction models of Table 1. As expected, extinction is mostly unimportant in the $K$-band, except at $z=3$ where most of the stars are very young. The agreement with the observations, shown as symbols, is excellent. This is not necessarily surprising as the free parameters in Shark are chosen to provide a good fit to the $z=0,\,1,\,2$ stellar mass functions, which are strongly correlated with the rest-frame $K$-band luminosity. The tension seen at $z=3$ can in part be due to the BC03 SPs having a small contribution from Asymptotic Giant Branch (AGB) stars. Other SP models, such as those of Maraston (2005), produce more $K$-band emission from AGB stars at $z\approx 3$ than BC03 (see Gonzalez-Perez et al. 2014 for a discussion).

Because all the attenuation models produce very similar $K$-band LFs, we show the contribution from disks, and bulges formed via galaxy mergers and disk instabilities only for the EAGLE$\tau$ RR14 attenuation model. Galaxy disks tend to dominate at the faint-end, with the luminosity below which they dominate becoming fainter as the redshift increases. Bulges driven by disk instabilities have a contribution to the $K$-band luminosity that increases strongly with time. At $z=3$, bulges built via disk instabilities make only a small contribution throughout the magnitude range studied here; as time passes by, they become more important, and by $z=0.5$ they play a significant role in shaping $L^{*}$. Bulges built via galaxy mergers on the other hand dominate the number density of galaxies over the whole magnitude range at $z=3$, but their dominance shifts to brighter luminosities at lower redshifts. Note, however, that they always play an important role, even at the faintest magnitudes, contributing $\approx 15-25$% of the observed K-band luminosity in galaxies with $-20<M_{\rm K,rest-frame}(AB)<-16$. The integrated rest-frame $K$-band luminosity of galaxies is dominated by bulges even out to $z=8$. We come back to this in $\S$ 5.

The left panels of Fig. 14 show the total rest-frame UV LF evolution from $z=3$ to $z=10$ in Shark using the $4$ extinction models of Table 1. We show both the intrinsic emission and the one after attenuation. The latter is the one that should be compared to observations. A general trend obtained for all models is that the attenuation in the brightest UV galaxies at $z=3$ and $z=4$ tends to be extremely large, reaching even $\approx 3-4$ mags in some cases, a lot higher than what the values in Shark at $z=0$ (see Figs. 8 and 10). This shows that the extinction of the most star-forming galaxies tends to increase from $z=0$ out to $z=3$ and decrease towards higher redshift. This evolution is driven by these galaxies at $z=2-3$ being on average more dusty than those at $z=0$: they have dust surface densities peaking at higher values than at $z=0$ at fixed stellar mass (see Fig. 2), and a tail of galaxies with extremely large dust surface densities, $\Sigma_{\rm dust}>10^{10}\rm M_{\odot}/kpc^{2}$.

Comparing the different attenuation models of Table 1, it is clear that the model EAGLE-$\tau$ RR14-steep provides the best agreement with the observations at all the redshifts of Fig. 14. This is because this model produces the smallest $\tau$ in galaxies with $Z_{\rm gas}<0.5$, which most Shark galaxies are at $z>3$. The largest differences between models is seen for the bright galaxies, those with UV magnitudes $\lesssim-20$ mag. These galaxies have on average $0.25<Z_{\rm gas}/Z_{\odot}<0.7$, which in the models EAGLE-$\tau$ RR14 and EAGLE-$\tau$ $f_{\rm dust}$-const have the Milky-Way dust-to-metals ratio, while in the EAGLE-$\tau$ RR14-steep model can have $>10$ times less dust per metals mass. Although a different dependence of the dust-to-metal ratio on gas metallicity could provide a better fit to the observations, we decide not to force the agreement and simply explore whether local Universe empirical relations allow Shark to provide a reasonable match. We caution the reader, however, that the effect of cosmic variance in the observations is very large, which for the area of the Hubble Deep Field ($2.6$ arcmin2) is $\approx 77$% at $z\approx 4$ according to the cosmic variance calculator of Driver & Robotham (2010). The latter is generally not included in the errorbars of the observations.

We remind the reader that we are assuming the dust-to-metal mass ratio to be invariant with time. Vijayan et al. (2019) included explicit dust formation and destruction in the SAM L-galaxies and predict the dust-to-metal ratio to evolve strongly, with $z=8$ and $10$ values being about $1.5$ dex smaller than $z=0$ values at fixed stellar mass, which agrees with the observational inferences of De Vis et al. (2019). This not necessarily unexpected, as some sources of dust formation, such as AGB stars and formation in molecular clouds require at least few $100$ Myr before they start to contribute. If we were to apply such an evolution, our fit to the UV LF would improve. However, other SAMs, e.g. (Popping et al., 2017), after implementing similar models of dust formation and destruction find little to no evolution of the dust-to-metal ratio. These contradictory results therefore merit caution in using these relations.

Other SAM results for the UV LF at high redshift (e.g. Qiu et al. 2019; Yung et al. 2019) provide better fits to the UV LFs than those in Fig. 14. However, they tend to be tuned to the UV LFs at $z>3$, and it is unclear whether these models reproduce the panchromatic SEDs of galaxies and the lower redshift Universe observations simultaneously.

In the middle and right panels of Fig. 14 we split the UV LF into the contributions from galaxy disks and bulges, respectively. It is clear that the largest differences at $3\leq z\leq 6$ between different attenuation models in the total UV LF mostly come from how they predict the extinction for disks, with variations of up to $1.5$ mags at $z=3$ and $2$ mags at $z=6$ between the EAGLE-$\tau$ RR14-steep and the other models. Note that at the faint end, magnitudes $>-17$, the EAGLE-$\tau$ RR14-steep and EAGLE-$\tau$ RR14 extinction models converge to the same answer, as these galaxies have $Z_{\rm gas}<0.25\,\rm Z_{\odot}$. By $z=8$ and $10$ the EAGLE-$\tau$ RR14-steep predicts almost no extinction in the case of disks, and hence there are only marginal differences between the intrinsic and attenuated UV LFs of disks in this model. The EAGLE-$\tau$ $f_{\rm dust}$-const model produces a disk UV LF that is similar to the one obtained with the default CF00 parameters.

We shift our focus now to bulges, which at these redshifts mostly correspond to central starbursts, and are the main channel of bulge formation. At $z=3,4$, all the EAGLE-$\tau$ extinction models produce more extinction than the default CF00 model, and in fact there are little differences between the three EAGLE-$\tau$ models. This is because these starbursts have on average $Z_{\rm gas}>0.7\rm\,Z_{\odot}$. At $z\geq 6$ there are some significant differences, with the attenuation model EAGLE-$\tau$ RR14-steep producing much smaller attenuation, due to these starbursts having $Z_{\rm gas}<0.7\,\rm Z_{\odot}$. Note, however, that even at $z=8$ and even at $z=10$, the extinction in starbursts galaxies is predicted to be significant, with typical values at the bright end of $\gtrsim 2$ mags.

In Fig 15 we compare the predicted UV slopes of galaxies with an AB rest-frame UV magnitude of $-19.5\pm 0.2$, which we measure by fitting the spectrum in the range $0.1\mu\rm m<\lambda_{\rm rest}<0.3\mu\rm m$ with the function $\epsilon_{\nu}\propto\lambda^{2+\beta_{\rm UV}}$, which is equivalent to the fitting performed in observations with the flux in the wavelength space $f_{\lambda}\propto\lambda^{\beta_{\rm UV}}$. The two attenuation models based on RR14 produce similar evolution but with a zero-point offset of $\approx 0.3$. The other two attenuation models, CF00 and EAGLE-$\tau$ $f_{\rm dust}$ const, produce weaker redshift evolution. We compare with the observations of Bouwens et al. (2014) and find that the attenuation model EAGLE-$\tau$ RR14-steep, which reproduces the UV LFs the best, also reproduces the observed UV slopes very well. This is very encouraging as it shows that an attenuation model based on local Universe dust-to-metal scaling relations is capable of reproducing the UV emission of galaxies even out to $z=10$.

5 Number counts and the Cosmic SED across cosmic times

5.1 Number counts

Galaxy number counts are the most direct observable of galaxies: how many galaxies are observed in a given apparent magnitude in a given band. Because galaxies of different masses and at different cosmic epochs contribute to this observable, they have been difficult to reproduce in galaxy formation simulations (see Somerville et al. 2012; Lacey et al. 2016 for a discussion). Another obvious difficulty is that constructing number counts necessarily requires to predict the galaxy population over the entire age of the universe and in a wavelength range as wide as possible.

With the SED models presented here we can test Shark against the observed number counts. To do this, we build a lightcone of area $107$ deg2 including all galaxies with a dummy magnitude, computed assuming a stellar mass-to-light ratio of $1$, $<32$ and at $0\leq z\leq 8$. We then use the method described in $\S$ 3 to build SEDs. We refer to Chauhan et al. (2019) for more details about our lightcone construction. Fig. 16 shows the predicted number counts from the NUV to the $850\mu$m of this lightcone for the $4$ attenuation models of Fig. 16, compare with the observations of Driver et al. (2016b) and Geach et al. (2017).

The agreement with the observations is excellent across the entire wavelength range shown here and for all the attenuation models tested. Some tension is identified in the Herschel SPIRE bands, in which Shark tends to predict too few (many) galaxies with AB magnitudes $14-16$ ($<10$) by a factor of $\approx 2$ compared to Driver et al. (2016a). Interestingly, these differences are similar to those reported in Lacey et al. (2016) for the GALFORM semi-analytic model. Recently, Wang et al. (2019) showed that the Herschel number counts we show here likely suffer from systematic errors due to blending and confusion, and hence the tension with Shark could be due to those systematics.

The truly unexpected result of Fig. 16 is that we are able to match the observed number counts in the UV-optical and FIR bands simultaneously without the need to invoke a varying IMF. Baugh et al. (2005); Lacey et al. (2016) showed that in GALFORM this was only achieved by invoking a top-heavy IMF during starbursts. In the case of a universal IMF the numbers of bright $850\mu$m galaxies in their work was consistently under-produced, and not only that, but they tended to lie at low redshift, in clear tension with the observations (which find a peak at $z\approx 2$). Shark assumes a universal Chabrier (2003) IMF and hence this shows that in a fully cosmological galaxy formation model, this is possible. In order to confirm this claim, we show in Fig. 17 the predicted redshift distribution of bright $850\mu$m galaxies, fluxes $>5$ mJy, for the $4$ attenuation models of Fig. 1, compared to the observations of Wardlow et al. (2011). The agreement is outstanding with all the models that use the EAGLE attenuation curves, while for the model adopting the default CF00 parameters, the redshift distribution is less peaked at $z\approx 2$ than observations suggest. In any case, Shark captures well the redshift peak of the brightest $850\mu$m sources, and the tail towards high redshifts. We remind the reader that in all cases we assume an invariant relation between dust mass-gas metallicity and gas content that is informed by local Universe observations.

The reasons why Shark is able to reproduce the observed number counts from the UV to the FIR with a universal IMF and other models have not is difficult to pinpoint due to the many aspects that enter in the calculation: dust masses, gas metallicities, galaxy sizes, attenuation curves and dust temperature. Hence we here discuss some possibilities but warn the reader that these are not conclusive. An important quantity is the dust mass, which is tied to the gas metallicities and gas content. Both Shark and GALFORM reproduce well the gas content of galaxies, however, GALFORM predicts gas metallicities that are consistently too low compared to observations at $M_{\star}\lesssim 10^{10.5}\,\rm M_{\odot}$ by up to $1$ dex (see Fig. 11 in Guo et al. 2016). Galaxy sizes may also be too large in GALFORM compared to observations (see Fig. 21 in Lacey et al. 2016). Both these effects contribute to lowering the dust surface density. Shark on the other hand predicts sizes that agree with observations (by construction), and gas metallicities that are closer to those observed (but not perfect; see Figs. 10 and  15 in Lagos et al. 2018). We are also assuming two constant dust temperatures for the BCs and diffuse dust, while in GALFORM this is computed self-consistently, which produces a dust temperature that weakly increases with redshift (Cowley et al., 2017). The latter makes the $850\mu$m emission weaker at fixed total FIR luminosity. A definitive conclusion though is that the answer to whether a varying IMF is needed to reproduce simultaneously the UV-optical and FIR emission of galaxies or not is model dependent.

Fig. 16 also shows the contribution from star formation in disks and in bulges, the latter separated by triggering mechanism: galaxy mergers and disk instabilities. We show this for the EAGLE-$\tau$ RR14 steep model only for the sake of clarity. As expected, the NUV is dominated by star formation in disks over the whole magnitude range, while the $r$-, $Y$ and IRAC $3.6\mu$m bands are dominated by bulges at bright magnitudes, transitioning to disks dominating at fainter magnitudes. The exact transition is wavelength dependent, moving from $\approx 18$ AB magnitudes in the $r$-band to $22$ in IRAC $3.6\mu$m. In the FIR the opposite trend takes place: going from the IRAC $8\mu$m to the $850\mu$m bands, we see the transition from bulge-dominated to disk-dominated emission moving to fainter magnitudes, with a transition of $12$ AB magnitudes in the IRAC $8\mu$m band to $18$ mags at $850\mu$m. In Shark, bright $850\mu$m galaxies (also referred to SMGs) are a mix of starbursts driven by galaxy mergers and by disk instabilities in almost equal numbers, with a slight dominance of galaxy mergers.

5.2 The Cosmic SED

The integrated spectrum of galaxies at a given redshift is termed the cosmic SED (CSED), and holds important information of the star formation activity of galaxies, the amount of light that is absorbed and reprocessed by dust, and the type of galaxies that contribute to the light at different wavelengths.

In this section we compare our predictions with the observations of Andrews et al. (2017), which are based on the GAMA survey (Driver et al., 2009), as well as the re-analysis of the G10/COSMOS photometry and spectroscopy (Davies et al., 2015). These measurements are for $0\leq z\leq 1$, and hence any higher redshift result can be considered a prediction of Shark.

We compute the predicted CSEDs of Shark by simply adding the light from all the galaxies at any given redshift. Truncating the integration to AB magnitudes $<-14$ does not have an effect on the predicted CSED, which shows that the integral is well converged for the resolution of our simulation.

5.2.1 The effect of extinction in shaping the CSED

Fig. 18 shows the predicted CSEDs of Shark at $z=0.25,1,3,8$ for the four attenuation models of Table 1. All the models predict a similar FIR CSED that at $z=0.25$ and $z=1$ agree reasonably well with the observations of Andrews et al. (2017). The models tend to produce too much UV by $\approx 30-50$% at $z=0.25$ compared to observations, though at $z=1$ the agreement is excellent. Given all the modelling that goes into predicting the UV, such as the adopted IMF, SP templates, SFH, ZFH and dust attenuation, and the effects in observations that are more difficult to include in the errorbars, such as cosmic variance, the UV LF faint-end slope and uncertain extrapolations, we consider this level of disagreement to be acceptable.

Some important differences among models are seen at high redshift; by $z=8$, there are differences of up to $0.4$ dex in the power output at fixed wavelength at $10^{5}<\lambda_{\rm rest}/$Å$<10^{6}$. This is due to the large differences in extinction predicted by our attenuation models in galaxies with gas metallicities $<0.75\,\rm Z_{\odot}$. The NIR is consistent among all the attenuation models at all redshifts. This is not surprising as the light at these wavelengths tends to trace stellar mass closely, which is the same for all models.

In the UV-end of the CSED, all models predict an important steepening of the UV slope with increasing redshift. Observations of individual high redshift galaxies show similar steepening of the UV compared to local galaxies (e.g. Dunlop et al. 2012). Although the overall trends are qualitatively the same for the four attenuation models studied here, in the detail there are some important differences. In order to quantify them, we measure the UV slope of the CSED at different redshifts and show them in Fig. 19. The EAGLE-$\tau$ RR14 extinction model produces the strongest evolution with a difference of $1.8$ in $\beta_{\rm UV}$ between $z=0$ and $z=10$. We show in Fig. 19 the observational constraint of Davies et al. (2013) from stacking of Ly-break galaxies, which seem consistent with the predictions of all the EAGLE-$\tau$ attenuation models. The default CF00 attenuation model produces the weakest evolution, and in fact the values of $\beta_{\rm UV}$ at $z>5$ in this model are too large compared to Davies et al. (2013).

Cowley et al. (2019) analysed the CSED predictions of the GALFORM semi-analytic model, and unlike Shark, they find little evolution of $\beta_{\rm UV}$, with values that throughout redshift are close to $-2$. This is in clear tension with the observations at low redshift, as seen in Figs. 18 and 20, but at high redshift they are consistent with those in Shark (albeit some of our attenuation models produce bluer spectra). Somerville et al. (2012) presented CSEDs using the Santa-Cruz SAM, and although they did not quantify the UV slope, their results seem to qualitatively support a strong redshift evolution of $\beta_{\rm UV}$.

5.2.2 Breaking down the Light budget in the CSED across cosmic time

Fig. 20 shows the predicted CSED of Shark with the default EAGLE-$\tau$ RR14 attenuation model at $0.25\leq z\leq 1$. Small diamonds show the intrinsic emitted light, while bigger diamonds show the predicted light after we include the effects of attenuation and re-emission in the IR. We find that Shark predicts a CSED that overall agrees very well with the observations through the whole wavelength range tested here, within the observational uncertainties. The level of agreement displayed by Shark is unprecedented to the knowledge of the authors. Cowley et al. (2019) showed for the GALFORM semi-analytic model that their model variant with a universal IMF struggled to simultaneously reproduce the FUV-to-optical and FIR parts of the CSED, and a top-heavy IMF was required. Because our Shark model assumes a universal IMF, it suggests that this may be model dependent. This agrees with the findings discussed in $\S$ 5.1. Baes et al. (2019) presented the CSEDs of the EAGLE hydrodynamical simulations and showed excellent agreement at $z\lesssim 0.5$, but towards $z\approx 1$ they found EAGLE to produce too little FIR emission. Hence, we consider the agreement seen in Fig. 20 to be a key success of Shark. Some areas of tension at the $0.1-0.15$ dex level, however, remain. At $z=0.25$, Shark produces too much FUV emission, and at $z=0.5$, Shark tends to produce $0.1$ dex too much emission in the optical-to-NIR bands.

Fig. 20 shows the contribution from disks and bulges of galaxies to the total CSED. Bulges tend to dominate in the optical-to-NIR wavelength range at $z\leq 1$, while disks dominate in the FUV-NUV and FIR ranges. The importance of bulges in the FIR emission, however, evolves strongly with redshift. This is because at $z=0$ we transition from bulges with no or little star formation, to centrally concentrated starbursts at $z\gtrsim 1$, which tend to be very dusty (see Fig. 3).

Fig. 21 shows the evolution of the CSED of Shark using the EAGLE-$\tau$ RR14 attenuation model at $2\leq z\leq 8$. At these redshifts the FIR makes a more significant contribution to the integrated light than the FUV-NIR, with the peak of the CSED being at $10^{5.5}\lesssim\lambda_{\rm rest}/$Å$\lesssim 10^{6.2}$. The slope of the CSED in the FUV-to-optical wavelength range becomes increasingly steeper with increasing redshift, due to both the very high star-formation activity in galaxies and their low metal and dust content (see Fig. 3).

At $z\gtrsim 2$, bulges make up most of the FIR emission, due to their starburst and dusty nature, and their contribution continues to increase with increasing redshift. Disks, on the other hand, dominate at the FUV-NUV over the whole redshift range, and by $z=8$ they also dominate in the rest-frame $u$ and $g$ bands. Note that at the NIR, bulges dominate throughout the whole redshift range analysed here $0\leq z\leq 8$.

6 Conclusions

We presented an exhaustive analysis of the spectral energy distribution (SED) predictions of the Shark semi-analytic model (Lagos et al., 2018) at $0\leq z\leq 10$. We first introduced the modelling of galaxy’s SEDs, which make use of the ProSpect software tool, which takes as input the SFH and ZFH of galaxies, and uses the Bruzual & Charlot (2003) SPs to produce the intrinsic emitted light. We then use the parametric attenuation curves of Charlot & Fall (2000) to compute the amount of extinction, and re-emit that in the IR following the templates of Dale et al. (2014) and energy conservations arguments. For the latter, we adopt an effective dust temperature for the diffuse ISM and birth clouds of $\approx 20-25$ K and $\approx 50-60$ K, which are fixed for the whole redshift range analysed in this paper.

To compute the appropriate Charlot & Fall (2000) extinction parameters of individual Shark galaxies, we make use of the predicted attenuation curves of the RT analysis of EAGLE by Trayford et al. (2019) and how these vary with the dust surface density of galaxies. We compute the dust content of Shark galaxies by applying the local Universe scaling relation between the dust mass, gas content and gas metallicity of Rémy-Ruyer et al. (2014), and assume this relation to hold out to $z=10$. This method allows us to apply a physical model for the attenuation of UV-to-optical light and re-emission in the IR that scales with galaxy properties. After generating the FUV-to-FIR emission of Shark galaxies, we compare to observations without re-tuning the model.

We summarize our findings below:

• •

Our model is capable of reproducing the wide diversity of observed galaxies, from galaxies that are almost metal free and have negligible attenuation, which tend to be abundant at high redshift and at low stellar masses, to SMGs, which are most prominent at around $1\leq z\leq 3$, but exist in the model out to $z=6$.

• •

We tested different models for the conversion of gas mass and gas metallicity to dust mass within the observational uncertainties and find that these tend to produce different FUV LFs with the largest differences appearing at $z\geq 4$. Differences in the optical-to-NIR are negligible throughout $0\leq z\leq 10$, and in the FIR they are only important at faint magnitudes, below the current observational limits.

• •

Shark is capable of reproducing well the observed $z\approx 0$ LFs of galaxies from the FUV (GALEX) to the FIR ($850\mu$m). We compare our model with observed LFs in $27$ bands and found reasonable agreement in all of them. In a future paper (Bravo et al. in prep) we show that optical colours are also very well reproduced even at intermediate redshifts.

• •

We analysed the rest-frame K-band and UV LFs out to $z=3$ and $z=10$, respectively, and found Shark to reproduce them reasonably well. We find that the rest-frame K-band LF above the knee is always dominated by bulges in galaxies while the rest-frame UV LF sees a strong evolution, from being dominated by star-forming galaxy disks throughout most of the magnitude range at $z\lesssim 4$ to a bigger contribution from low metallicity galaxy mergers-induced starbursts at the bright end at $z\gtrsim 4$. UV-bright galaxies display a strong evolution of their UV slope from $\approx-0.2$ at $z=0$ to $\approx-2.5$ at $z=10$, with some variations between the different adopted dust-to-gas mass scalings. We find the attenuation of UV-to-optical light to be maximal at $z\approx 1-2$.

• •

By building a deep, wide area lightcone of $107$ deg2 with Shark galaxies, we compare the predicted number counts from the NUV to the $850\mu$m with observations and find unprecedented agreement. To confirm our SMG population is realistic, we also study the redshift distribution of bright, $>5$ mJy, SMGs and found that it peaks at $z\approx 2$ with a tail that extends out to $z\approx 6$, in very good agreement with observations. This is achieved without the need of invoking a top heavy IMF in starbursts and/or a redshift-dependent dust-gas mass-gas metallicity scaling showing that a fully cosmological galaxy formation model is capable of reproducing simultaneously the emission in the UV-optical to the FIR with a universal IMF.

• •

We integrate the galaxy LFs at different redshifts to produce a CSED from $z=0$ out $z=10$ and find Shark to reproduce well the observed CSEDs at $z\leq 1$, while there are no available observations at higher redshifts. Shark predicts the FIR emission to be dominated by star-forming disks at $z\lesssim 1.5$, and by starbursts at higher redshifts, even out to $z=10$. These starbursts are triggered by both disk instabilities and galaxy mergers, and we find that they contribute similarly to the IR emission. The rest-frame UV and NIR are dominated by star-forming disks and bulges at all redshifts, respectively.

The success of our model makes it an ideal tool for future galaxy surveys from $z=0$ to $z=10$. Possible applications include understanding the galaxy populations different color-based selections isolate, how observationally-based environment metrics trace the underlying halo population, the bias of flux-selected galaxies in different bands, systematic effects in photometric redshift determinations, among many others. The interested reader is encouraged to contact the authors of this manuscript for access to the simulated lightcones.

One of the most surprising results in this manuscript is the fact that we can simultaneously reproduce the UV-to-NIR and FIR properties of galaxies, including number counts and redshift distributions, without the need of varying the IMF of galaxies, which is unprecedented. The reason why previous models struggled with this and Shark does not is difficult to pinpoint as these models are complex and commonly a combination of processes are responsible for the differences seen among simulations. However, we discussed several possibilities, which we plan to explore in depth in the future, including (i) differences in the predicted gas metallicities and sizes among models (both of which affect the dust surface density) (ii) differences in the SFR function, particularly at $1\leq z\leq 3$, (iii) differences in the dust temperature evolution, and (iv) different attenuation curves. Nonetheless, we can certainly assert that the answer to what physical processes are required to simultaneously reproduce the FUV-to-FIR emission of galaxies is model dependent.

Acknowledgements

We thank Cedric Lacey, Carlton Baugh and Desika Narayanan for useful discussions about the results in this paper, and the anonymous referee for their constructive report. CL has received funding from the ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013. CL also thanks the MERAC Foundation for a Postdoctoral Research Award. JT and CL also thank the University of Western Australia for a Research Collaboration Award which facilitated the face-to-face interaction that led to this work. This work was supported by resources provided by The Pawsey Supercomputing Centre with funding from the Australian Government and the Government of Western Australia. Cosmic Dawn Centre is funded by the Danish National Research Foundation.