G.A.S. II: Dust extinction in galaxies; Luminosity functions and InfraRed eXcess

TL;DR

This paper introduces a detailed dust attenuation model in the G.A.S. semi-analytical galaxy formation framework, predicting UV and IR luminosity functions and exploring dust effects on galaxy observables across cosmic time.

Contribution

It presents a comprehensive dust model with three components and a two-phase ISM, improving predictions of galaxy luminosity functions and IRX relations over previous models.

Findings

Accurately reproduces UV and IR luminosity function evolution from z=9 to z=0.1.

Predicts large scatter in IRX-$eta$ relations, influenced by dust and stellar age.

Identifies the slope of the attenuation curve as more dependent on FUV attenuation than disk inclination.

Abstract

Dust is a crucial component of the interstellar medium of galaxies. The presence of dust strongly affects the light produced by stars within a galaxy. As these photons are our main information vector to explore the stellar mass assembly and therefore understand a galaxy's evolution, modeling the luminous properties of galaxies and taking into account the impact of the dust is a fundamental challenge for semi-analytical models.We present the complete prescription of dust attenuation implemented in the new semi-analytical model: G.A.S. This model is based on a two-phase medium originating from a physically motivated turbulent model of gas structuring (G.A.S. I paper). Dust impact is treated by taking into account three dust components: Polycyclic Aromatic Hydrocarbons, Very Small Grains, and Big Grains. All three components evolve in both a diffuse and a fragmented/dense gas phase. Each phase has its own stars, dust content and geometry. Dust content evolves according to the metallicity of it associated phase.The G.A.S. model is used to predict both the UV and the IR luminosity functions from $z=9.0$ to $z=0.1$. Our two-phase ISM prescription catches very well the evolution of UV and IR luminosity functions. We note a small overproduction of the IR luminosity at low redshift ($z<0.5$). We also focus on the Infrared-Excess (IRX) and explore its dependency with the stellar mass, UV slope, stellar age, metallicity and slope of the attenuation curves. Our model predicts large scatters for relations based on IRX, especially for the IRX-$\beta$ relation. Our analysis reveals that the slope of the attenuation curve is more driven by absolute attenuation in the FUV band than by disk inclination. We confirm that the age of the stellar population and the slope of the attenuation curve can both shift galaxies below the fiducial star-birth relation in the IRX-$\beta$ diagram.

Tables

Table 1. Table 1: Dust properties. The dust type and size law refer to DustEM keywords.

Type	$\rho$ [$g/c{m}^3$]	size range [cm]	size law
PAH0	$2.24$	$3.1\times {10}^{-8}$ - $1.2\times {10}^{-7}$	mix-logn
PAH1	$2.24$	$3.1\times {10}^{-8}$ - $1.2\times {10}^{-7}$	mix-logn
amCBE	$1.81$	$6.0\times {10}^{-8}$ - $2.0\times {10}^{-6}$	logn
amCBE	$1.81$	$4.0\times {10}^{-7}$ - $2.0\times {10}^{-4}$	plaw-ed
aSil	$3.00$	$4.0\times {10}^{-7}$ - $2.0\times {10}^{-4}$	plaw-ed

Equations

$$Z_g=12+log(O/H)=12+lo{g}_{10}\left(\frac{M_O}{M_H}\frac{m_H}{m_O}\right)$$

$${\eta }_m=\frac{MIN(M_{1/2,1}\text{leavevmode};\text{leavevmode}{M}_{1/2,2})}{MAX(M_{1/2,1}\text{leavevmode};\text{leavevmode}{M}_{1/2,2})},$$

$$\frac{f_{PAH}}{f_{PAH,0}}=10^{-11.0+1.30Z_g(-0.2)}$$

$${\tau }_{\lambda }=\sum _{gt}{f}_{gt}\times \left({\sigma }_{abs,gt}+{\sigma }_{sca,gt}\right)$$

$${\omega }_{\lambda }=\frac{{\sum }_{gt}{f}_{gt}\times {\sigma }_{sca,gt}}{{\tau }_{\lambda }}.$$

$$A_{\lambda }=-2.5\text{leavevmode}lo{g}_{10}{\Phi }_{\lambda }({\tau }_{\lambda },{\omega }_{\lambda })$$

$${\Phi }_{\lambda }^{SPH}=\frac{a_{\lambda }}{1-{\omega }_{\lambda }+{\omega }_{\lambda }{a}_{\lambda }},$$

$$a_{\lambda }=\frac{3}{4k_{\lambda }}\left[1-\frac{1}{2k_{\lambda }^2}+\left(\frac{1}{k_{\lambda }}+\frac{1}{2k_{\lambda }^2}\right)exp(-2k_{\lambda })\right],$$

$$k_{\lambda }={\tau }_{\lambda }^{SPH}={\tau }_V^{SPH}\times \frac{{\tau }_{\lambda }}{{\tau }_{\lambda }({\lambda }_V)}.$$

$${\tau }_V^{SPH}=2.619{\left(\frac{{\chi }_Z}{{\chi }_{\odot }}\right)}^{1.6}\left(\frac{N_H}{2.1\times 10^{21}at\cdot c{m}^{-2}}\right)$$

$${\Phi }_{\lambda }^{ISM}=\frac{1-exp\left(-\sqrt{1-{\omega }_{\lambda }^{ISM}}{\tau }_{\lambda }^{ISM}/cos\text{leavevmode}i\right)}{\sqrt{1-{\omega }_{\lambda }^{ISM}}{\tau }_{\lambda }^{ISM}/cos\text{leavevmode}i}$$

$${\tau }_{\lambda }^{ISM}={\tau }_{FUV}^{ISM}\times \frac{{\tau }_{\lambda }}{{\tau }_{\lambda }({\lambda }_{FUV})}$$

$${\tau }_{FUV}^{ISM}=\left(1.926+0.051\times {\Sigma }_H^{ISM}\right)\times 10^{[0.947\times Z_g^{ISM}-9]}$$

$${\Phi }_{\lambda }^{GMC}=exp\left(-\sqrt{1-{\omega }_{\lambda }}{\tau }_{\lambda }^{GMC}\right)$$

$${\tau }_{\lambda }^{GMC}={\tau }_{FUV}^{GMC}\times \frac{{\tau }_{\lambda }}{{\tau }_{\lambda }({\lambda }_{FUV})}$$

$$\begin{array}{rl}(\lambda I_{\lambda })& ={\Phi }_{\lambda }^{ISM}{\Phi }_{\lambda }^{GMC}(\lambda I_{\lambda }{)}_{GMC}\\ & +{\Phi }_{\lambda }^{ISM}(\lambda I_{\lambda }{)}_{ISM}\\ & +{\Phi }_{\lambda }^{SPH}(\lambda I_{\lambda }{)}_{SPH}\end{array}$$

$$L_{\star }=\int (I_{\lambda }{)}_{\star }d\lambda \text{leavevmode}\text{leavevmode}[L_{\odot }]$$

$$L_{IR}=L_{\star }-\int (I_{\lambda })d\lambda \text{leavevmode}\text{leavevmode}[L_{\odot }],$$

$${\Phi }_{\lambda }^{eff}=\frac{(\lambda I_{\lambda })}{(\lambda I_{\lambda }{)}_{\star }}\text{leavevmode}\text{leavevmode}\text{leavevmode}\mathrm{and}{A}_{\lambda }^{eff}=-2.5lo{g}_{10}{\Phi }_{\lambda }^{eff}.$$

$$f_{ext}=1.0-\frac{\sum L_{FUV}^{ext}}{\sum L_{FUV}^{int}}$$

$$\frac{f_{PAH}}{f_{PAH,0}}=10^{-17.1+2.0(Z_g-0.2)}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\mathrm{steeper}\text{leavevmode}\text{leavevmode}\mathrm{envelope}$$

$$\frac{f_{PAH}}{f_{PAH,0}}=10^{-1.6+0.2(Z_g-0.2)}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\mathrm{flatter}\text{leavevmode}\text{leavevmode}\mathrm{envelope}$$

$$f_{PAH}^n=f_{PAH}^{n-1}+r\times (f_{PAH}^{†}-f_{PAH}^{n-1}).$$

$$f_{BG}^n=f_{BG}^{n-1}+r\times (f_{BG}^{†}-f_{BG}^{n-1})$$

$$f_{VSG}^n=1-f_{PAH}^n-f_{BG}^n$$

Full text

*Aims. *We present the complete prescription of dust attenuation implemented in the new semi-analytical model (SAM): G.A.S. . This model is based on a two-phase medium originating from a physically motivated turbulent model of gas structuring (G.A.S. I paper).

*Methods. *Dust impact is treated by taking into account three dust components: Polycyclic Aromatic Hydrocarbons, Very Small Grains, and Big Grains. All three components evolve in both a diffuse and a fragmented/dense gas phase. Each phase has its own stars, dust content and geometry. Dust content evolves according to the metallicity of it associated phase.

*Results. *The G.A.S.  model is used to predict both the UV and the IR luminosity functions from z=9.0 to z=0.1. Our two-phase ISM prescription catches very well the evolution of UV and IR luminosity functions. We note a small overproduction of the IR luminosity at low redshift ($z<0.5$). We also focus on the Infrared-Excess (IRX) and explore its dependency with the stellar mass, UV slope, stellar age, metallicity and slope of the attenuation curves. Our model predicts large scatters for relations based on IRX, especially for the IRX-$\beta$ relation. Our analysis reveals that the slope of the attenuation curve is more driven by absolute attenuation in the FUV band than by disk inclination. We confirm that the age of the stellar population and the slope of the attenuation curve can both shift galaxies below the fiducial star-birth relation in the IRX-$\beta$ diagram. Main results presented in this paper (e.g. luminosity functions) and in the two other associated G.A.S.  papers are stored and available in the GALAKSIENN library through the ZENODO platform.

1 Introduction

In galaxies, stars, gas, and dust are mixed and strongly interact. Dust plays a major role in the chemistry and the physics of the interstellar medium (ISM). Whereas stars produce Ultra-Violet (UV) and visible and near-infrared (NIR) light, dust, heated by stellar radiation, radiates from the near-InfraRed (nIR) to the sub-millimeter (submm). Attenuation of the stellar light due to dust dramatically modifies fluxes, colors and spectral energy distributions of galaxies. Measurements of star formation rates (SFRs) and even stellar mass (e.g. Lo Faro et al. 2017) are therefore strongly affected. A large fraction of the UV radiation coming from newly formed O and B stars is obscured by dust (e.g. Casey et al. 2014; Burgarella et al. 2013; Takeuchi et al. 2005) and re-emitted mainly in the Mid- and Far-InfraRed (IR) and the sub-millimeter (submm) domain. This process is particularly efficient in dusty star forming galaxies (DSFGs), that can reach IR luminosities up to $10^{13}L_{\odot}$ and thus SFRs of thousands of solar masses per year (e.g. Casey et al. 2012; Vieira et al. 2013; da Cunha et al. 2015).

To take into account the overall star formation activity, both UV and IR contributions must be considered. By combining both measurements, Burgarella et al. (2013) inferred the redshift evolution (0 $<$ z $<$ 4) of the total cosmic star formation activity. Overall attenuation affecting the galaxies is often parametrized by the IR excess: $IRX=L_{IR}/L_{UV}$, where $L_{IR}$ is the integrated infrared luminosity between 8 and 1000 microns and $L_{UV}$ is the rest-frame 1500 Å  luminosity.

A large number of studies have measured the IRX. In the nearby Universe (e.g. Martin et al. 2005b; Bothwell et al. 2011), half of the total UV budget is obscured. At higher redshift, the fraction of the star formation activity directly traced by the UV strongly decreases, being less than 15% at $z>1$ (e.g., Takeuchi et al. 2005; Cucciati et al. 2012). Due to the decrease of the metallicity in very young galaxies, at $z>$3, it is expected that the UV flux traces a higher fraction of the star formation activity ($<$ 20%, Reddy et al. 2008, Bouwens et al. 2010). Thanks to recent observational progresses, building large UV-selected samples is feasible over most of the evolution of the Universe (e.g. Bouwens et al. 2012; Ellis et al. 2013; Reddy et al. 2012) and the majority of studies are based on (rest-frame) UV selected samples.

At high redshift (z $>$ 3) or in fields where IR/submm measurements are not available, only the stellar emission (UV-NIR) attenuated by dust is observed. Estimating accurately the corrections for this attenuation is thus crucial to derive SFRs. One of the most popular method to derive dust obscuration is based on the observed correlation between the slope of the UV continuum ($\beta$) and IRX. Initially measured and calibrated in nearby starbursts (e.g. Meurer et al. 1999), this method has been extensively used at all redshifts. However, it was found to be invalid for nearby secular star forming galaxies (e.g. Boissier et al. 2007; Hao et al. 2011). The relation between $\beta$ and IRX is also affected by selection effects (e.g. Buat et al. 2005) and could also be sensitive to dust properties (e.g. Calzetti 2001; Inoue et al. 2006), star formation histories (e.g. Kong et al. 2004; Boquien et al. 2012; Forrest et al. 2016) and shape/slope of the attenuation law (e.g. Salmon et al. 2016; Reddy et al. 2016; Bouwens et al. 2016).

In parallel to the IRX$-\beta$ relation another correlation has been measured between IRX and the stellar mass (e.g. Pannella et al. 2009; Buat et al. 2012; Finkelstein et al. 2012; Heinis et al. 2014; Whitaker et al. 2014). This correlation is currently not fully understood. In this context, stellar mass could be a crude tracer of the metal enrichment process through the complex star formation histories. Indeed, work by Cortese et al. (2006), Boquien et al. (2009), or more recently Pannella et al. (2015) has revealed a correlation between the dust attenuation and the gas-phase metallicity, which is itself correlated with the stellar mass (e.g. Mannucci et al. 2010; Lilly et al. 2013; Zahid et al. 2014).

Understanding dust attenuation in galaxies is still an important challenge and therefore a fundamental aspect of SAMs. The impact of dust has been studied in SAMs with different methods (e.g. Devriendt et al. 1999; Granato et al. 2000; Baugh et al. 2005; Bower et al. 2006; De Lucia & Blaizot 2007; Fontanot et al. 2009; Guo et al. 2011; Lacey et al. 2016) and yields to predicted number counts in UV and/or luminosity functions that differ by a factor $\simeq$2. The recent work of Narayanan et al. (2018), based on hydrodynamic zoom-in simulations and radiative transfer, and the analytical prescription proposed by Popping et al. (2017) highlight the complex effects of the average stellar population age, attenuation curves, and dust/stars mixing. However, the analytical model by Popping et al. (2017) cannot take into account complex star formation histories and, even if the hydrodynamic simulations from Narayanan et al. (2018) have been performed from z $\simeq$ 5 to z $\simeq$ 2, their limited number of simulated galaxies does not allow to explore the redshift evolution of the IRX-$\beta$ relation.

In this paper (paper II of the G.A.S. model presentation set), we present the set of prescriptions implemented in the G.A.S. model to describe the effects of the dust attenuation onto the stellar light. Physical prescriptions associated to the gas physics is describe in Cousin et al. (2018) (paper G.A.S. I).

We investigate in this paper the redshift evolution of the UltraViolet (UV) and InfraRed luminosity functions. We also explore the main drivers of both the IRX-$M_{\star}$ and $IRX$-$\beta$ relations in a realistic cosmological context from z $>$ 6 to z $\simeq$ 1.5. We explore the impact of the stellar age, the star formation rate, the absolute attenuation, the disk inclination and the the gas-phase metallicity. We will question the validity of the IRX-$\beta$ relation to correct UV observed flux from the dust attenuation.

The paper is organized as follows: in Sect. 2, we briefly describe our G.A.S. model. In Sect. 3 we describe dust components. We also present the recipes implemented to take into account dust attenuation and explore the impact of geometry and absolute attenuation. We present the evolution of the attenuation law predicted by our model as function of both the absolute attenuation in FUV band and disk inclination. We also compare our attenuation curves with standard prescriptions. In Sect. 4 we compare our predictions for the UV and IR luminosity functions to a large set of observational measurements. In Sect. 5 we focus on correlations between IRX, stellar mass, stellar population age, UV luminosity, gas phase metallicity, UV slope. We conclude in Sect. 6.

2 The G.A.S. model

Our work is the continuation of the G.A.S. semi-analytical model described in the G.A.S. I paper (Cousin et al. 2018). G.A.S. is applied to dark-matter merger trees extracted from a pure N-body simulation. This simulation uses a WMAP-5yr cosmology ($\Omega_{m}=0.28$, $\Omega_{\Lambda}=0.72$, $f_{b}=0.16$, $h=0.70$ Komatsu et al. (2009)) and describes a volume of $[100/h]^{3}Mpc$ containing $1024^{3}$ particles. Each of these particles has a mass of $m_{p}=1.025\leavevmode\nobreak\ 10^{8}\leavevmode\nobreak\ M_{\odot}$. Dark-matter halos and sub-structures (satellites) are identified by using the HaloMaker code (Tweed et al. 2009). In merger tree structures, we only take into account halos with at least twenty dark-matter particles leading to a minimal dark matter halo mass of $2.050\times 10^{9}\leavevmode\nobreak\ M_{\odot}$.

2.1 General description

Galaxy assembly starts from gas accretion. In G.A.S.halos are fed through a cosmological smooth metal-free baryonic accretion following the smooth dark-mater accretion. The total accreted gas is divided in two channels – a cold and a hot mode – depending on the dark-matter halo mass. While in the cold mode the gas directly feeds the galaxies, the shock-heated gas of the hot mode has to cool down and condensate before feeding the galaxy (G.A.S. I paper, Sect. 2). We assume that the warm condensed gas feeds a gaseous disk with an exponential profile. Gas acquires angular momentum during the transfer (Peebles 1969) and the disk is then supported by this angular momentum (e.g. Blumenthal et al. 1986; Mo et al. 1998).

In the disk the accreted gas is initially assumed to be mainly diffuse and non-star-forming. It goes through a progressive fragmentation following a turbulent cascade and is finally converted in star-forming gas. Galaxies modelled by G.A.S. host therefore two distinct ISM phases: a diffuse warm phase and a fragmented dense cold phase (G.A.S. I paper, Sect. 3 and 4). Young stars are born in the dense, structured gas. We assume that they stay in this dense phase during $5\times 10^{7}$ yrs. Older stars then evolve in the diffuse warm phase of the disk.

Our model takes into account supernovae (SNs) and active galaxy nuclei (AGN) feedback. A part of the power produced by these two types of feedback is directly injected in the interstellar medium of the galaxy and is used to disrupt the fragmented dense gas. The remaining power allows to generate large-scale ejecta and heat the gas. The gas ejected from the galaxy ends up in the hot gas atmosphere surrounding the galaxy. These feedback mechanisms mainly affect the evolution of low-mass and intermediate-mass galaxies ($<10^{9.5}M_{\odot}$) (G.A.S. I paper, Sect. 5).

In addition to these processes, we introduce a new regulation mechanism operating in the surrounding hot gas phase, that allows to reduce gas cooling in massive halos. It is based on thermal instabilities (G.A.S. I paper, Sect. 6).

Using these new gas-cycle and gas-regulation prescriptions, the G.A.S. model is able to reproduce the stellar mass function of galaxies from z $\simeq$ 6 to z$\simeq$0.5.

2.2 The chemodynamical model

A chemodynamical model (Cousin et al. 2016) is used to track the mass evolution (e.g. $M_{O}(t)$) of six of the main ISM elements: ${}^{1}H$, ${}^{4}He$, ${}^{12}C$, ${}^{14}N$, ${}^{16}O$, and ${}^{56}Fe$. Masses are tracked separately: on one side in the diffuse and the fragmented/dense gas phase of the disc, on the other side in the diffuse gas phase of the central bulge (if any). The production and the re-injection of these elements are taken into account for stars with initial mass between 0.1$M_{\odot}$ and 100$M_{\odot}$ and for metal-free to super-solar metal fractions. We assume that new stars are formed in the dense/fragmented gas phase following a Chabrier (2003) initial mass function (IMF). Metals are injected initially into the dense/fragmented gas phase. The kinetic energy coming from both SNs and AGNs partly disrupts the dense/fragmented gas phase and injects metals in the diffuse gas phase. Then metals evolving in both the diffuse and the dense/fragmented gas are strewn in the hot atmosphere surrounding the galaxy with the large scale ejecta. Due to this metal enrichment, the efficiency of the radiative cooling process affecting the hot surrounding gas is modified and impact the gas accretion.

Oxygen is the most abundant element formed into stars. It is commonly used as a tracer of the gas-phase metallicity. We define the gas metallicity $Z_{g}$ of a given phase as the number of oxygen atoms per hydrogen atoms with a logarithmic scale:

[TABLE]

$M_{O}$ and $M_{H}$ are the oxygen and the hydrogen mass contained in a given gas phase. $m_{O}$ and $m_{H}$ are the atomic mass of oxygen and hydrogen, respectively. In this formalism we adopt $Z_{\odot}=8.94$ (Karakas 2010).

2.3 Stellar spectra

Similarly to the gas metallicity, the G.A.S. model tracks the evolution of the spectral energy distribution (SED) of stellar populations. The star formation history of a galaxy component (disk or central bulge) is continuously followed according to the stellar age and the initial stellar metallicity. Thus, at each step of the evolution, we associate to a given stellar population its mass-weighted stellar SED: $(\lambda I_{\lambda})_{\star}$ based on Bruzual & Charlot (2003) SED libraries.

2.4 Merger events

We define the merger type (major or minor) of two galaxies 1 and 2 using the following mass ratio:

[TABLE]

where $M_{1/2}$ is the total mass (galaxy and dark-matter halo) inside the galaxy half-mass radius. We consider the merger event as minor when $\eta_{m}<1/3$.

In the case of a major merger ($\eta_{m}\geq 1/3$), the galaxy structure and the dynamics are strongly modified. All the gas (from the disk and the bulge) is attributed to the remnant disk.

We assume that young stars ($\leq 5\times 10^{7}$yrs) of the two progenitor disks are still embedded in GMCs. Young stars are therefore kept into the remnant disk.

Conversely, old stars ($>5\times 10^{7}$yrs) coming from the two progenitors are transferred to the central bulge. The stellar mass distribution of the central bulge is modeled by a Hernquist (1990) profile.

Between two merger events, all the gas produced by stellar winds from stars in the bulge is kept in the bulge. The mass that is generated is in average low ($M_{g}<10^{5}M_{\odot}$) but the metallicity can be relatively high ($Z_{g}>8.0$).

During a minor merger, we add separately the stellar population of the bulges and disks to form the remnant bulge/disk. All gas reservoirs are added into the remnant disk.

3 Dust in galaxies

In a galaxy, the overall extinction is the result of different factors: i) the amount/composition of dust in the different components of the galaxy (disk, bulge); ii) the relative contribution of the various stellar populations (in age and metallicity) in the various components; iii) the relative geometry of the dust and stars.

3.1 Dust composition

As proposed by for example Draine & Li (2001), Zubko et al. (2004),Draine & Li (2007), Compiègne et al. (2011) (and references therein), dust can be described by a set of three different grain types: Polycyclic Aromatic Hydrocarbons (PAH), Very Small Grains (VSG) and Big Grains (BG). In G.A.S.  we assume that these three dust types are formed in all galaxies at all redshifts. However, the relative mass fraction of this three dust types can vary from one environment to another (Sect. 3.2). The intrinsic properties of dust assumed in our model are listed in Table 1. These properties are fully compatible to those used in previous work (e.g. Zubko et al. 2004; Draine & Li 2007; Compiègne et al. 2011).

3.2 Dust content

In the interstellar medium, gas and dust are closely linked and the dust composition (PAH:$f_{PAH}$, Big Grains:$f_{BG}$, Very Small Grains:$f_{VSG}$) depends on the gas metallicity.

Fig. 1 shows measurements of the PAH fractions as a function of the gas metallicity performed by Draine et al. (2007), Ciesla et al. (2014) and Rémy-Ruyer et al. (2015).

These studies used different prescriptions for gas metallicity measurements. As shown by Kewley & Ellison (2008), the prescription proposed by Pilyugin & Thuan (2005) (used in Rémy-Ruyer et al. (2015)) leads to lower values than all other prescriptions. A systematic shift of +0.2 dex applied to Rémy-Ruyer et al. (2015) measurements leads to a good agreement between the studies.

In the two distinguished ISM phases hosted by our modeled galaxies, the fraction of PAH formed in a given gas-phase is settled by using Rémy-Ruyer et al. (2015) (their Eq. 5):

[TABLE]

where $f_{PAH,0}=4.57\%$ is the reference value for our Galaxy (Zubko et al. 2004).

Eq. 3 is used to generate a continuous and smooth evolution of the PAH fraction in relation to the gas phase metallicity. Around the mean relation, we introduce a scatter by applying between two time-steps a random walk between the current value and a target value given by Eq. 3 (for more details, see Appendix A.2). We plot in Fig. 1 the whole distribution of $f_{pah}$ obtained following this random walk procedure and measured in the dense/fragmented gas phase. The scatter obtained around the mean relation is fully compatible with observational measurements.

BG represent the main part of the dust mass (e.g. Zubko et al. 1996; Draine & Li 2001; Compiègne et al. 2011), therefore when the PAH fraction is settled, we then treat the BG mass fraction. This fraction is assumed to be around to $2/3$ of the residual mass. This fraction is typical of what is measured in our Galaxy (e.g. Zubko et al. 2004; Compiègne et al. 2011). As for PAH, we computed the fraction of BG following the random-walk technique. At each time step, the mass fraction of BG is randomly settled between the current value and the target value ($2/3$ of the residual mass). Finally, mass conservation requires that VSGs are made of the rest of the dust mass.

The dust composition of each gas phase (diffuse gas in disc and bulge, dense/fragmented gas in disc), is derived independently according to its metallicity, using Eq.3 and our random-walk technique (Eq. 1).

The empirical relation described by Eq. 3 cannot fully account for the complex dust processing in the ISM: fragmentation, destruction by UV photons, heating, or interstellar shocks (Jones et al. 1994, 1996; Draine & Li 2007; Relaño et al. 2018). When $f_{PAH}$ increases or decreases between two time-steps, according to Eq. 3 and the random-walk technique, the chosen ”effective evolution” hides the complexity, which cannot be treated in full details by G.A.S.İndeed, even if G.A.S.tracks the kinetic energy produced by SNs and AGN, this energy is only used to (re-)distribute the gas between the different gas phases and does not affects directly the dust contents/composition.

Moreover, G.A.S.assumes a single relation characterized by a smooth and continuous trend from low to high metallicities. As pointed by Rémy-Ruyer et al. (2015) or Draine & Li (2007), the PAH fraction is low in dwarf galaxies and can be larger in more evolved galaxies. The $Z_{g}$–$f_{PAH}$ relation assumed in G.A.S.(Eq. 3) follow the trend of the majority of current observational measurements. We show in Appendix A.1 that the variations of this relation have only a small impact onto the dust attenuation.

3.3 Extinction curves

Once the dust composition ($f_{PAH}$, $f_{BG}$, $f_{BG}$) is fixed in the two ISM phases, we build the associated extinction curve. For that, we take into account absorption and scattering properties of dust through the DustEM model (Compiègne et al. 2011).

Fig. 2(a) shows the absorption and the scattering curves associated to the different dust types.

We define the extinction curve $\tau_{\lambda}$ as:

[TABLE]

The sum is performed over the three grain types and $f_{gt}$ is the mass fraction associated to a given grain type (see Sect. 3.2).

Based on the scattering and extinction curves, we define the albedo $\omega(\lambda)$,

[TABLE]

Fig. 2(b) shows the extinction and albedo wavelength dependencies computed for different PAH fractions. The set of extinction curves produced by G.A.S. is compared to LMC, SMC and Milky Way extinction curves measured by Mathis et al. (1983) and Gordon et al. (2003). In the wavelength range $\lambda\in[0.1,2.0]\,\mu$m, according to the PAH abundance, our model allows to reproduce a large range of slopes. However, our assumption about the compositions of dust, PAH, VSG and BG leads especially to Milky-Way like extinction curves.

As highlighted in our Fig. 2(b), both the extinction ($\tau_{\lambda}$) and the albedo ($\omega(\lambda))$ curves are very sensitive to dust composition/properties. These effects are discussed in the literature (e.g. Seon & Draine 2016; Chastenet et al. 2017) and are strongly linked to the choices of dust content/composition. They have to be kept in mind.

3.4 Overall attenuation, dust geometrical distribution

In a galaxy, the overall attenuation affecting the stellar radiation is governed by the dust content but also by the geometrical distribution of the dust inside the galaxy. On the one hand, the impact of the dust content is quantified by the face-on FUV/optical depth. On the other hand, the geometrical distribution of stars and dust is taken into account through different geometrical functions $\Phi_{\lambda}(\tau_{\lambda},\omega_{\lambda})$. $\Phi_{\lambda}$ allows to build the effective attenuation:

[TABLE]

through the extinction curve $\tau_{\lambda}$ and the albedo curve $\omega_{\lambda}$. We define in the next sections different geometrical functions $\Phi_{\lambda}$ which are associated to the different environments/ISM-phases or galaxy components (disk/central bulge) of our modeled galaxies.

In Devriendt et al. (1999) or De Lucia & Blaizot (2007) the effective attenuation is calibrated according to the Visible (V) band, which is linked to both gas surface density and metallicity (Guiderdoni & Rocca-Volmerange (1987)). However, the analysis done in Boquien et al. (2013) (their Fig. 2) indicates that attenuation, measured in the Far Ultra-Violet (FUV) band, is more reliable than in the V-band. They show that the relative scatter of the attenuation in the FUV-band is 30% to 50% smaller than in the V-band. These results are only valid for disk galaxies and we thus applied a FUV-band normalization based on Boquien et al. (2013) to the disks only. We kept the V-band scaling for the bulges.

3.4.1 Contents and morphologies of galaxies; impact on the effective attenuation.

In a given galaxy, the gas/dust distribution is strongly linked to the morphology. Each galaxy is composed of a disk and a potential central bulge. In each disk, the G.A.S. model follows the evolution of two main gas phases: a diffuse warm gas hosting the oldest stars ($>5\times 10^{7}$ yrs) and $N_{GMC}$ Giant Molecular Clouds hosting the youngest stars.

Three different geometrical functions are associated to these three different environments. For the bulge component we adopt a Dwek geometry (Sect. 3.4.2). We use the classical Slab (Sect. 3.4.3) geometry for the diffuse ISM. In addition to this Slab geometry, young stars evolving in GMCs are described with a screen geometry (Sect. 3.4.4). In all these cases, we assume that the dust content can affect the stellar radiation of the component only if the gas metallicity is sufficiently high: 12 + log(O/H) $>$ 6.5, else we neglect the dust attenuation.

3.4.2 Bulge environments: the Dwek geometry

As explained in Sect. 2.4, during a major merger old stars are transferred into a central bulge. Some gas, with metals, is produced by stellar evolution into the bulge. dust formed in this gas phase generate an extinction of the stellar light. We assume that stars and dust are homogeneously mixed and we use the Dwek geometrical function (Lucy et al. 1989; Dwek & Városi 1996; Devriendt et al. 1999):

[TABLE]

where

[TABLE]

and

[TABLE]

$\tau^{SPH}_{V}$, the reference depth in (V-band), varies with the metal mass fraction111We assume a solar metal mass fraction: $\chi_{\odot}=0.02$ $\chi_{Z}$ and the hydrogen column density $N_{H}$ following (Guiderdoni & Rocca-Volmerange 1987; Devriendt et al. 1999):

[TABLE]

The hydrogen column density is computed within the half mass radius and we assume an homogeneous gas metallicity in the bulge.

3.4.3 disk diffuse ISM environments: the Slab geometry

To describe the dust attenuation in the diffuse warm ISM of each galaxy disk we use the standard Slab geometry. This distribution is defined as an infinite plane-parallel in which stars, gas and dust are perfectly mixed over the same scale. We use:

[TABLE]

where $i$ is the inclination angle of the disk with respect to the $z$ axis of the simulated box. $\tau_{\lambda}^{ISM}$ is defined as followed:

[TABLE]

The extinction curve $\tau_{\lambda}$ is built according to the dust composition (see Sect. 3.2, Eq. 4) of the diffuse ISM. This composition depends of the diffuse gas phase metallicity $Z_{g}^{ISM}$. The face on UV depth $\tau_{FUV}^{ISM}$ is calculated as in Boquien et al. (2013) (their Eq. 13):

[TABLE]

To compute the hydrogen surface density $\Sigma_{H}^{ISM}$ we assume that half of the hydrogen mass, is stored in the half mass radius of the galaxy. $Z_{g}^{ISM}$ is settled to the gas-phase metallicity specifically associated with the diffuse warm gas contained into the disk. In Eq. 11, the albedo is taken into account through the factor $\sqrt{1-\omega_{\lambda}}$.

3.4.4 Additional attenuation from Giant Molecular Clouds

New stars are embedded in some GMCs. These dense structures generate an additional attenuation that we model by applying a Screen geometrical function:

[TABLE]

Albedo $\omega_{\lambda}$ and extinction curve $\tau_{\lambda}^{GMC}$ are computed following the dust composition of the structured gas phase. $\tau_{\lambda}^{GMC}$ is normalized to the FUV band thought:

[TABLE]

The normalization factor $\tau_{FUV}^{GMC}$ is settled with Eq. 13. For that, all GMC formed in a galaxy are assumed to be similar. For a GMC, the average hydrogen surface density is computed by assuming that all the hydrogen mass is stored in $N_{GMC}$ spherical clouds of radius $R_{GMC}=3h/8$ (see Cousin et al. (2018)). Normalization also depends on the structured gas metallicity $Z_{g}^{GMC}$. Metallicity in the dense phase is slightly higher than in the diffuse medium.

3.5 The attenuated stellar spectra

A galaxy contains up to three distinct stellar populations associated with three different geometrical distributions. We build up the three spectral energy distributions (SED) associated with each stellar population:

• •

the young stellar population, in the disk: $(\lambda I_{\lambda})_{GMC}$,

• •

the old stellar population in the disk: $(\lambda I_{\lambda})_{ISM}$,

• •

the old stellar population in the central bulge: $(\lambda I_{\lambda})_{SPH}$.

The overall extinguished stellar spectrum is then given by

[TABLE]

3.6 Infrared energy budget

The bolometric luminosity produced by stars is:

[TABLE]

where $(I_{\lambda})_{\star}$ is the non-extinguished stellar SED (Sect. 2.3).

By assuming that dust in the diffuse ISM and in GMC absorbs a fraction of this radiation and re-emits it the IR range, the energy balance leads to a total IR luminosity given by:

[TABLE]

where $(I_{\lambda})$ is the extinguished stellar SED build throught Eq. 16.

3.7 Effective dust attenuation

The attenuation is a combination of the three distinct attenuations associated with the distinct stellar populations. For simplicity, and to allow a comparison with observations, we define the effective attenuation:

[TABLE]

Once dust/stars geometry and dust properties are fixed, the FUV depth (used as normalization) and the disk inclination are the two last key parameters for the effective attenuation. Fig. 3 shows the impact of disk inclination. As expected (e.g. Chevallard et al. 2013), face-on disks show a steeper attenuation curves than edge-on disks. Between these two extreme cases we observe a slow evolution of the slope of the attenuation curve. We compare these median attenuation curves to both Calzetti et al. (2000) and Reddy et al. (2016) laws. In the wavelength range $\lambda\in[0.5-1]\,\mu$m our median attenuation curves are in excellent agreement with Calzetti et al. (2000) measurements. Our median attenuation curves corresponding to the highest inclinations are above that of Calzetti et al. (2000). For $\lambda\in[0.1-0.5]\,\mu$m and for the majority of disk inclinations, our median attenuation curves are steeper than the Calzetti et al. (2000) law and closer to those measured by Reddy et al. (2016). However, our median attenuation curves are shallower than the attenuation law measured by Reddy et al. (2016). The sample of galaxies used to explore the evolution of the effective attenuation curve are randomly selected between z = 4.0 and z = 1.0 and for stellar mass $M_{\star}>10^{7}M_{\odot}$. It is important to note that for this large stellar mass window, our results do not vary with the redshift.

In a given disk inclination bin, the variation of the optical depth in FUV band leads to a large scatter around the median trend. This scatter can be larger than the one resulting from the variation of inclination. To illustrate this scatter we select galaxies with inclinations between $30^{\circ}$ and $50^{\circ}$. In this sub-sample, we show in Fig.4 the trend of the attenuation curve as a function of the absolute attenuation in the FUV band ($A_{FUV}$). As expected, a lower $A_{FUV}$ is associated with a steeper attenuation curve. In the wavelength range $\lambda\in[0.1-0.5]\,\mu$m, the median trend of this sub-sample of galaxies lies between the Calzetti et al. (2000) and the Reddy et al. (2016) attenuation curves. Galaxies experiencing a very low FUV attenuation can have a steeper attenuation curve than the measurements of Reddy et al. (2016). In contrast, galaxies with a strong attenuation in FUV band show a shallower attenuation curve than Calzetti et al. (2000) law. According to these trends, the absolute extinction in FUV band seems to be the main driver of the extinction even if it is difficult to extract a unique main contributor driving the attenuation curve slope.

The observed scatter of the distribution of our effective attenuation curves is fully consistent with recent models implemented in hydrodynamic simulations, especially the recent study performed at z = 5.0 by Cullen et al. (2017) (their Fig. 12) in the context of the First Billion Years projects.

4 Luminosity functions

We now compare our predictions with some luminosity function measurements. We focus on the UV Luminosity Function (UVLF) and the IR Luminosity Functions (IRLF). Luminosity functions are built for stellar masses greater than $10^{7}M_{\odot}$. This mass threshold corresponds to our stellar mass resolution. Close and below this mass, our predictions are incomplete and suffer of numerical instabilities. Particularly, gas-phase metallicity cannot be computed and dust composition cannot be defined. According to the star formation histories and the evolution of the stellar populations metallicity, this mass threshold corresponds to different UV magnitudes and IR luminosities limits and varies with z. In Figs. 5 and 6 we present luminosity function only above the turnover, which is specific of the incompleteness of the data.

4.1 UV luminosity functions

4.1.1 The overall evolution

Fig. 5 shows the predictions for the UVLF at different redshifts. UV Magnitudes have been computed by convoluting the rest-frame spectrum of each galaxy with the GALEX (Martin et al. 2005a) FUV transfer function222Obviously, similar results are obtained by convoluting for example i, r, or u filter bands with the redshifted galaxy spectrum.. We compare our predictions with a large set of observational data coming from Arnouts et al. (2005), Oesch et al. (2010), Cucciati et al. (2012), McLure et al. (2013), Oesch et al. (2013), Alavi et al. (2014), Parsa et al. (2016), Hagen et al. (2015), Bouwens et al. (2015), Finkelstein et al. (2015), Atek et al. (2015) and Livermore et al. (2017).

For all redshifts, we show the attenuated and intrinsic luminosity function. At all redshifts, the attenuated UVLFs and measurements are in very good agreement. However, at z$\simeq$9.0 and at z = 3.0 our predictions are slightly below the observations, especially in the bright end part of the luminosity function at z = 3.0 and in the low-luminosity range at z$\simeq$9.0. At z$\simeq$9.0 the number of galaxies taken into account in the predicted UVLF is low ($<$500). At these very high redshifts, measurements are based on only few galaxies and may suffer from strong uncertainties. UVLFs are therefore affected by statistical noise. At $z\simeq$3.0, UVLFs are built with a larger number of objects: the discrepancy between observational measurements and prediction cannot be only explained by uncertainties. At $z\simeq 3.0$ specifically, modelled galaxies seem to be affected by a too strong attenuation, especially in the high luminosity (and also mass) range. This trend is also confirmed by the behavior of the K-band luminosity function at this redshift (Cousin, in preparation). The IR luminosity function (Fig. 6) indicates that the IR energy absorbed by dust in this high luminosity regime is in agreement with observations. In addition, the stellar mass function presented in paper I also shows a good agreement with measurements. The discrepancy noted at $z\simeq 3.0$ cannot be therefore currently explained.

At all redshifts and for the range of magnitudes considered here, the shape of the intrinsic UVLF is compatible with a power law. At highest luminosities a break appears. It is only due to the exponential decrease of the number of galaxies present in the volume. On the contrary, the shape of the extinguished UVLF evolves with both the redshift and the UV luminosity. At high redshifts ($z\geq 6.0$), the shape of the extinguished UVLF is still compatible with a power law but with a steeper slope than the un-extinguished UVLF. For redshifts z$<$6.0, the extinguished UVLF clearly shows two different shapes. In the low luminosity range $M_{UV}>-20$, the shape is compatible with a power law. However, in the high luminosity range ($M_{UV}<-20$), we find an exponential decrease.

4.1.2 How many energy is absorbed by dust ?

By comparing the attenuated (ext) and the intrinsic (int) UV luminosities of all the galaxies produced by our model at a given redshift, it is possible to compute the fraction of UV radiation absorbed by dust. At a given redshift, we sum separately the extinguished and the intrinsic UV luminosities over all galaxies. We define the UV extinguished fraction as the ratio of these two quantities:

[TABLE]

This fraction is indicated in each panel of Fig. 5. Our model predicts an UV extinguished fraction of 40% at z$\simeq$9.0. At this redshift, only the most massive galaxies ($M_{\star}>10^{10.5}M_{\odot}$) are affected by dust extinction. For the galaxies with lower mass, the average metallicity is too low. The fraction of UV radiation absorbed by dust progressively increases when redshift decreases. At z$\simeq$8.2 our model predicts that roughly half of the UV radiation produced by young stars is absorbed by dust. This fraction reaches 60%, 70% and 80% at z $\simeq$ 7.0, 5.0 and 3.0. Obviously this progressive increase of the dust attenuation is linked to the progressive increase of the gas-phase metallicity in galaxies. At lower redshift (z$\leq$ 2.0), the fraction of UV radiation absorbed by dust still increase but much more slowly. At z$\simeq$0.3 the fraction of UV radiation absorbed by dust reaches a value of $\simeq 86$%.

By integrating a compilation of UV and IR luminosity functions, Burgarella et al. (2013) found a similar ratio of 85$\pm$10 % and 84$\pm$13 % at z=2.2 and z=3.15, respectively. These results agrees at less than 1$\sigma$ with our model, which predicts 82 and 79 %, respectively. At higher redshifts, only weak constraints are available in the far-infrared (Madau & Dickinson 2014), since the limited depth of the instruments does not allow us to constrain accurately the position of the knee of the luminosity function. Burgarella et al. (2013) used the attenuations estimated by Bouwens et al. (2009) to put constraints at higher redshifts and found that the equality between UV and IR (fext = 50 %) is reached at z$\simeq$6. Several teams also tried using various modeling approaches to estimate the evolution of the ratio between obscured and unobscured star formation density. A similar value is found by the empirical model of Béthermin et al. (2017) and in the measurement of the obscured star formation rate density from the CIB anisotropies from Maniyar et al. (2018). In contrast, it disagrees with Koprowski et al. (2017) who predicts a transition around z$\simeq$3 based on their empirical modeling approach.

4.2 IR luminosity functions

4.2.1 The evolution of the total IR luminosity function

Fig. 6 shows our predictions for the total IR luminosity function at different redshifts. We compare our results with measurements coming from Caputi et al. (2007), Rodighiero et al. (2010), Gruppioni et al. (2013), Magnelli et al. (2013), Marchetti et al. (2016) and Koprowski et al. (2017).

Our predictions are in general in very good agreement with observations. At all redshifts presented here, the total IR luminosity function predicted by our G.A.S. model is fully compatible with a double power-law. A knee between the two regimes (low and high luminosities) appears around $L_{IR}\simeq 10^{12}L_{\odot}$ at z$\simeq$6.0. The break progressively shift to lower luminosity. At $z\leq 1.0$ the knee appears around $L_{IR}\simeq 10^{11.5}L_{\odot}$.

For the redshifts z$\simeq$2.8 and $\simeq$3.7, even if the bright ends are in good agreement with the recent Koprowski et al. (2017) estimates, G.A.S. seems to over-predict the density of IR galaxies for $L_{IR}<10^{12}L_{\odot}$. However we also note strong tensions between different measurements at these high redshifts.

Even if measurements are not available at z $>$ 4.0, we show the G.A.S. predictions at z = 4.7, 6.1 and 6.9 in the three lower panels of Fig. 6. These IRLF highlight that IR luminosities higher than $10^{12}L_{\odot}$ can already be reached by galaxies at z$\simeq$7. Metal production and dust assembly is therefore already very active at this epoch. This prediction is in agreement with hydrodynamic simulations, other SAM taking into account dust obscuration (e.g. Clay et al. 2015) and also recent observations of high redshift (z$>$6) heavily-obscured galaxies (e.g. Bowler et al. 2018; Watson et al. 2015).

4.2.2 Contributions of the different environments, Is there an excess of cold dust at low redshift ?

In addition to the total IR luminosity function, we also plot in each panel of the Fig. 6 the luminosity function build only with the IR luminosity coming from GMC environments (young stars). The fraction $f_{GMC}$ of the IR luminosity coming from GMC environments is also indicated in each panel.

At high redshift ($z>3.5$) this fraction reaches very high values ($>$80, 95%). The large majority of the IR light comes from dense strongly attenuated GMCs. Between $z=1.5$ and $z=3.5$ this fraction is still larger than 60%. For $z<1.5$ the contribution of GMC progressively decreases and reaches $\simeq$ 53% at $z=0.9$. The contribution of the oldest stars evolving in the diffuse ISM appears progressively. However, the majority of the IR luminosity always comes from GMCs illuminated by young stars.

At $z\leq 0.6$ the fraction of IR luminosity coming from the GMC environment drops below 50% and even reaches 39% at z$\simeq$0.1. The contribution of old stars is therefore slightly dominant. In parallel, at $z=0.4$ and $z=0.1$, our G.A.S. model generates a slight over-density of IR galaxies. The IR emission is slightly too high for the full IR luminosity range explored. However, if we focus only onto the GMC contribution the agreement at these low redshift is very good. The IR emission coming from the cold dust component hosted by the diffuse ISM is probably too high. This trend is confirmed in the G.A.S.III paper. IR number counts and the comparison of the IR SEDs with recent measurements reveal an excess of IR light emitted by the cold dust component, mainly in the massive local ($z<0.5$) galaxy population.

4.2.3 The impact of the latest disruption episode on massive IR bright galaxies

At all redshifts explored here, the gap existing between the total and the GMC-only IR luminosity functions seems indicated that the contribution of the GMC to the total IR luminosity seems lower at the bright-end than at the faint-end of the IR luminosity function. To explore this trend we compute the contribution of GMCs for both the faint IR-galaxies ($L_{IR}<10^{11.5}L_{\odot}$) and the bright IR-galaxies ($L_{IR}>10^{11.5}L_{\odot}$). These median fractions are indicated in the bottom of each panel of Fig. 6.

At all redshifts shown here, the contribution of the GMCs ($f_{GMC}$) is found to be larger in luminous IR galaxies ($L_{IR}>10^{11.5}L_{\odot}$) than in faint IR galaxies ($L_{IR}<10^{11.5}L_{\odot}$). This trend is especially true at high redshifts. At z $\leq$ 1.5, the contributions of the two galaxy populations are closer and vary between 60% and 40%.

This trend is explained by the progressive evolution with the stellar mass (IR luminosity) of the dense/fragmented gas fraction. In the star-forming galaxy sample, massive/IR-luminous galaxies ($M_{\star}\geq 10^{10.5}M_{\odot}$, $L_{IR}\geq 10^{10.5}L_{\odot}$) the gas content is less structured/fragmented than in less massive and less IR-luminous galaxies. In massive galaxies, we observe a progressive decrease of the gas accretion. This leads to a decrease of the star formation activity. However, SN feedback still disrupts the dense/fragmented gas after the last star formation episode (G.A.S. I paper). The impact of the latest disruption episode impacting the GMCs leads to a more diffuse ISM and therefore to an attenuation slightly lower in massive (IR-luminous) galaxies than in the less massive ones.

5 Dust attenuation

The IRX can be used to measure the attenuation in FUV $A_{FUV}$. Many efforts have been made to link IRX (and therefore $A_{FUV}$ ) to galaxy properties like the slope of the UV spectrum, the stellar mass or the star formation rate (e.g. Meurer et al. 1999; Buat et al. 1999; Gordon et al. 2000; Bell 2002; Kong et al. 2004; Calzetti et al. 2005; Seibert et al. 2005; Cortese et al. 2006; Boissier et al. 2007; Salim et al. 2007, 2009; Treyer et al. 2007; Boquien et al. 2009; Buat et al. 2010; Takeuchi et al. 2010).

In this section we explore the variations of IRX with galaxy properties such as gas-phase metallicity and stellar population (age and mass). For the IRX-$\beta$ relation we also analyze its dependency with disk inclination and slope of the attenuation curve.

Our analysis is based on a star-forming galaxies sample. We only keep galaxies with a mass higher than $10^{8}M_{\odot}$. A star-forming galaxy is defined following the Main Sequence SFR–$M_{\star}$ relation calibrated by Schreiber et al. (2015). We consider that a galaxy is star forming if its specific star formation rate sSFR is greater333We checked that our results are not significantly affected if we use instead a threshold of 1/3 or 1/5. than $1/4\times(SFR_{MS}/M_{\star})$.

5.1 The IRX – AFUV relation

In star forming galaxies, the link between IRX and $A_{FUV}$ has been explored by Buat et al. (2005), Cortese et al. (2008) or Boquien et al. (2012) and also by using radiative transfer models (e.g. Witt et al. 2006).

Fig. 7 shows the IRX-AFUV relation based on a star forming galaxies sample extracted at z = 1.5. Our predictions agrees with the relation from Buat et al. (2005).

Even if the large majority of galaxies lies around the parametric relation, we see from Fig. 7 that a scatter progressively appears at high IRX and high AFUV. As highlighted on Fig.  7 this scatter is linked to the age of the stellar population hosted by galaxies. At a given IRX, a lower AFUV is associated with an older stellar population. Similarly, at a given AFUV, larger IRXs also corresponds to older stellar populations. We clearly see that the scatter of the relation increases with IRX. The variation with the luminosity-weighted stellar population age observed in our G.A.S model is in agreement with the analysis done by Cortese et al. (2008). They conclude that, at a given IRX, active star forming galaxies (with young stars) are more attenuated than passive galaxies (with old stars). In Fig.  7, our results are obtained at z=1.5, but similar trends and behaviors can also be extracted at other redshifts. We do not observe any variation of the IRX-AFUV relation with the redshift.

5.2 The IRX – gas metallicity relation

The production of dust in the ISM is closely linked to the metal enrichment of the gas phase and a correlation between IRX and gas phase metallicity is expected. Such a correlation has been found in nearby galaxies by e.g. Boquien et al. (2009). In their sample, galaxies with higher metallicities have higher IRX.

Fig. 8 shows the IRX–$Z_{g}$ relation built from our z=1.5 star forming galaxy sample. As discuss previously, the effective attenuation is mainly produced by the GMC environment associated with the dense/fragmented gas phase. The average gas metallicity used in Fig. 8 is therefore computed in the dense/fragmented gas phase according to Eq. 1.

In Fig. 8, the majority of our star forming galaxies lie at $Z_{g}<8.5$. In this low metallicity regime, the dispersion is high (close to 2 dex) and it is difficult to extract any trend between the IRX and the gas-phase metallicity.

However, at high metallicity, i.e. $Z_{g}>8.5$ we observe a clear correlation with the IRX increasing with the average gas phase metallicity. In this domain, the median trend is in good agreement with the measurements performed by Pannella et al. (2015). However, the data points obtained from star-forming galaxies evolving in nearby clusters by Cortese et al. (2006) mainly lie below our median trend. At a given metallicity, local star forming galaxies are therefore less attenuated. This behavior could be associated with a less structured and less dense ISM in low redshift galaxies. As discuss previously (Sect. 4.2), when the redshift decreases, we observed in our simulated galaxies a decrease of the structured gas fraction. In parallel we can also note an increase of the characteristic size (exponential radius) of disks. These trends are especially marked in massive galaxies (i.e with $Z_{g}>8.0$). Due to these two conjugated trends, the average gas density and the effective attenuation are therefore probably lower in low redshift galaxies ($z<1.0$) than at higher redshifts (here z = 1.5).

5.3 The IRX–M*⋆* relation

Our median IRX–$M_{\star}$ relation is plotted on Fig. 9. It agrees well with observational measurements (e.g. Heinis et al. 2014; Whitaker et al. 2014; Pannella et al. 2015). As observed, our model predicts a broad correlation between the stellar mass and the dust attenuation. In this context, the stellar mass has to be seen as a tracer of the galaxy activity, which is connected to the amount of metals produced by a galaxy since its formation (e.g. Cousin et al. 2016; Tremonti et al. 2004; Andrews & Martini 2013). Through the stellar mass, we roughly trace the dust enrichment of the ISM.

It is often difficult to get both UV and IR luminosities for individual objects, especially at high redshift. The IRX is therefore mainly measured using stacking techniques (e.g. Béthermin et al. 2012; Hilton et al. 2012): for a given stellar mass, IRX is defined as the ratio of the average IR luminosity to the average UV luminosity: $IRX=<L_{IR}>/<L_{UV}>$ (stacking estimator). We compute this estimator in our different stellar mass bins. As shown in Fig. 9, we do not observe strong differences between our median and this stacking estimator in the low mass range ($M_{\star}<10^{10}M_{\odot}$). At larger stellar masses $M_{\star}>10^{10.5}M_{\odot}$, our stacking estimator reaches lower IRX values than the median trend and therefore, in this high mass regime, values and shape of the stacking estimator are in better agreement with observational measurements performed by Heinis et al. (2014) and Pannella et al. (2015). In parallel, in this high mass regime, our median estimator shows a better agreement with Whitaker et al. (2014).

Around the median trend, we also plot in Fig. 9 the dispersion limits of the relation. These two limits correspond to the 15% and 85% percentiles of the distribution in a given stellar mass bin. For stellar masses $M_{\star}>10^{9.5}M_{\odot}$, we observe a large scatter around the median relation of about 1.2 dex in average (and therefore a average dispersion of 0.6 dex). This dispersion is in good agreement with the largest error bars of the measurements which are estimated to be 0.4-0.7 dex (Heinis et al. 2014; Pannella et al. 2015). We note from Fig. 9 that the dispersion mainly comes from the high IRX domain. Our G.A.S. model predicts IRX larger than $10^{3}$ for all the stellar mass domain explored. These very extinguished galaxies are probably lost by the UV-selected sample used by (Heinis et al. 2014).

5.3.1 Correlation with UV luminosity

In Fig. 10, we highlight the typical UV luminosity of galaxies depending on their position in the IRX–$M_{\star}$ diagram by selecting a sample uniformly in the IRX–$M_{\star}$ plan. Due to this homogeneous selection, the dispersion of these points is much larger than the scatter of the relation given above. This random sub-sample reveals a clear evolution of the extinguished UV luminosity with the IRX. At a given stellar mass, the observed UV luminosity decreases when the IRX increases in very good agreement with observational measurements of Heinis et al. (2014) (see their Fig. 6, C1 and C2).

As mentioned previously, Fig. 10 shows that our model predicts the existence of galaxies with a very strong attenuation (IRX $>10^{3}$). These galaxies are associated with $<10^{8}L_{\odot}$ observed UV luminosities, close from the lower limit of the Heinis et al. (2014) UV-selected sample. These very UV-faint galaxies can therefore explain the large scatter predicted by our model. To confirm this hypothesis, we limit our star-forming galaxies to galaxies with an extinguish UV luminosity higher than $L_{UV}>10^{8.5}L_{\odot}$. With this sample, we do not observe any variation of our mean or median IRX–$M_{\star}$ relation, which are dominated by luminous UV galaxies. However the scatter of the relation falls to 0.9 and therefore the dispersion to 0.45 dex which is fully compatible with dispersion measurements (Heinis et al. 2014).

5.3.2 Evolution with redshift

Fig. 11 shows the evolution of the IRX–M*⋆* relation with redshift (from z = 1.5 to z = 5.9) obtained using our stacking estimator. At all redshifts explored here, the scatter is similar and close to $\sim$1.2 dex.

At stellar mass above $10^{10.5}M_{\odot}$, we observe a decrease of the IRX with the stellar mass. The mass threshold of this decrease shifts to lower mass when the redshift is increasing (from $\simeq 10^{10.5}$ to $\simeq 10^{9.5}M_{\odot}$ from z = 1.5 to z = 6, respectively). This decrease at high stellar mass is also observed by Whitaker et al. (2014) with a trend similar to which is seen in our model. As previously (Sects. 4.2, 5.2), we explain this turn-over with the evolution of the contribution of the GMC environment. Massive galaxies are less-structured and therefore show a lower attenuation than in less massive and more structured galaxies.

For M${}_{\star}<10^{10.5}M_{\odot}$, at a given stellar mass, our model predicts an increase of the average IRX with the redshift. Even if we note some tension in the absolute IRX values measured444These studies use different observational data: Spitzer for Whitaker et al. (2014), Herschel for Heinis et al. (2014) and Pannella et al. (2015) by Whitaker et al. (2014) and Heinis et al. (2014), these two studies do not observe the increase predicted by our G.A.S. model. For example, the average relations given by Heinis et al. (2014) or Pannella et al. (2015) are valid from z = 4.0 to z = 1.5. The discrepancy is mainly visible a $z>3.0$. At lower redshift the large dispersion ($\simeq 1.1$ dex) can explain the discrepancy. At $z>4.0$, IR and therefore IRX measurements are not currently available on a sufficiently large sample of galaxies. It is therefore difficult to put strong constraints on the amount of energy effectively absorbed by dust. At these high redshifts our UVLF agrees well with the measurements. To reduce the effective IRX and keep the UVLF similar, the overall attenuation has to be slightly reduced and the amount of young stars also reduced such as to keep a constant apparent UV emission. However, the amounts of stars currently predicted by G.A.S.at these redshifts agree well with the most recent measurements. In that context, it is difficult to explain the possible discrepancy.

5.4 The IRX–$\beta$ relation

Fig. 12 shows the IRX–$\beta$ relation based on our star forming galaxies sample extracted at z = 1.5. To compute the UV slope ($\beta$) we have fitted the rest-frame continuum of each galaxy with a power law expressed as $F(\lambda)\propto\lambda^{\beta}$. We used a set of five photometric points associated to filters 1, 3, 5, 7 and 9 listed in Calzetti et al. (1994) (their Table 2). Similarly to what we did for the the IRX–M*⋆* relation, we derived both the median trend of the IRX and the stacking estimator, i.e. the ratio of the mean IR luminosity to the mean UV luminosity in a given bin of $\beta$.

Our two estimators show similar behaviors, but the median is systematically lower than the mean. Our sample is in fact dominated by low-SFR galaxies ($<1\,M_{\odot}/yr$) and therefore low-luminosity galaxies. The median estimator simply reflects this distribution, while the stacking mean estimator is more sensitive to high SFRs (luminosity) values. Note that this bias also affect observational measurements.

Our mean relation agrees well with measurements done by Heinis et al. (2013) and Takeuchi et al. (2012). It is important to note the large scatter of the distribution. Even if the majority of our simulated galaxies are located between $\beta\in$ [-1.5, 0.0] and log10IRX$\in$ [0.0, 2.0], some galaxies can have log10IRX $>$ 2.5. Our simulated galaxies located above the median relations are in a region that is known to be populated by IR bright galaxies at all redshifts (e.g. Casey et al. 2014; Howell et al. 2010; Reddy et al. 2012; Lo Faro et al. 2017).

5.4.1 Impact of the star formation history

The impact of the star formation history onto the IRX–$\beta$ relation has been originally studied by Kong et al. (2004) and then discussed and confirmed by several studies (Cortese et al. 2006, 2008; Boquien et al. 2009; Popping et al. 2017; Narayanan et al. 2018): older stellar populations lead to higher intrinsic values of $\beta$. Narayanan et al. (2018) also found that galaxies hosting older stellar populations are less obscured (the old stellar population evolves in an ISM, which is less dense than in GMC) leading to a lower IRX.

Star formation history is therefore a relevant parameter affecting the IRX–$\beta$ relation. However it is difficult to find a good tracer of the average star formation history (see Boquien et al. 2009 for a complete discussion about the D4000 and b parameters). In our case we use the luminosity-weighted age of the stellar population as tracer of the star formation history.

In Fig. 13, we highlight the luminosity-weighted age of a sub-sample of galaxies from our model uniformly selected in the IRX–$\beta$ plan. We confirm the trend found by the previous studies mentioned above: galaxies with older stellar populations preferentially lie below the canonical Meurer et al. (1999) IRX–$\beta$ relation.

5.4.2 Impact of the attenuation curve

Recent observations done with the Atacama Large Millimetre Array (ALMA) by e.g. Capak et al. 2015, Bouwens et al. 2016, and Pope et al. 2017 (z $>$ 4) or combined HST and Hershel/PACS measurements (Reddy et al. 2018) (z $>$ 1.5) have mainly populated the IRX-$\beta$ diagram below the canonical relation of Meurer et al. (1999). Those studies have increased the observed range of $\beta$ ($<$-1.0) at low IRX ($<$10). They concluded that such low IRX and $\beta$ values were consistent with a SMC attenuation curve (stepper than MW attenuation curve) and the low metallicity expected for the high-redshift galaxies.

Fig. 14 shows the IRX–$\beta$ relation at z = 1.5 with a color coding of the attenuation curve. The slope is computed using a linear regression in the log$\lambda$–log$(A_{\lambda}/A_{V})$ plan. The attenuation $A_{\lambda}/A_{V}$ is computed according to Eq. 19. The fitting procedure is applied with seven photometric points (FUV, NUV and filters 1, 3, 5, 7 and 9 of Calzetti et al. 1994). To allow for an easy comparison with the standard Calzetti et al. (2000) attenuation law, we show in Fig. 14 a value normalized to the slope of the Calzetti et al. (2000) attenuation law (-0.54 dex/$\mu$m). In other words, in Fig. 14 slope log${}_{10}\sigma>$ 0 (resp. log${}_{10}\sigma<$ 0) indicates a steeper (resp. flatter) attenuation curve than Calzetti et al. (2000) attenuation curve.

For comparison, in addition to our stacking estimator we also include in Fig. 14 the average IRX–$\beta$ relations computed by Salmon et al. (2016) and derived from attenuation curves estimated by i) Gordon et al. (2003) (SMC), ii) Calzetti et al. (2000) and Reddy et al. (2015). As shown on Fig. 4, slopes of these attenuation curves are different; Calzetti et al. (2000): -0.54dex/$\mu$m (our reference), Reddy et al. (2015): -0.64 dex/$\mu$m, Gordon et al. (2003): -1.23 dex/$\mu$m.

In Fig. 14, slopes of attenuation curves associated to our simulated galaxies and trends of average IRX–$\beta$ computed by Salmon et al. (2016) are in good agreement. We clearly observe that galaxies lying below the reference Meurer et al. (1999) relation have steeper attenuation curves (up to three times), consistent with the SMC-like attenuation curves. Our G.A.S. model is in perfect agreement with the recent studies performed by i.e. Salmon et al. (2016), Lo Faro et al. (2017), or Reddy et al. (2018) and with the simulations analyzed by Popping et al. (2017) or Narayanan et al. (2018).

On the contrary, galaxies lying above the reference relations have flatter attenuation curves (up to ten times, Lo Faro et al. 2017). As expected, for galaxies lying close to the reference relations, the slope of the attenuation curve is fully compatible with the one of the Calzetti et al. (2000) sample.

As shown in Figs. 3 and 4, the G.A.S. model allows us to produce a large set of effective attenuation curves agreeing with recent measurements. These effective attenuation curves can be flatter or steeper than the Calzetti et al. (2000) law. This diversity and the variation of the slope across the IRX–$\beta$ plan is fully linked to the capability of the G.A.S. model to associate with each dust contents/metallicities of each galaxy components (GMC, diffuse ISM) a dedicated extinction curve and therefore a dedicated effective attenuation.

5.4.3 Evolution with the redshift

Fig. 15 shows the IRX–$\beta$ relations for six different redshifts, from z = 1.5 to z $\simeq$ 6.0. We present relations built with the stacking estimator. At a given $\beta$-slope, we note an increase of the average attenuation with redshift (about 0.4 dex between z = 1.5 to z = 6.0). Using the 15% and 85% percentiles, we compute the scatter of the relation that slightly decreases with the redshift, from 0.5 dex at z $\simeq$ 6.0 to $\simeq$1.0 dex at z = 1.5. By taking into account these scatters, the slopes and behaviors of our IRX–$\beta$ relations are compatible with the measurements, especially with Heinis et al. (2013). While the measurements performed with ALMA by Koprowski et al. (2018) are above our mean relation at $z=3.35$, the measurements at $z=3.87$ and $z=4.79$ agree fairly well with our mean trend. As for the IRX–$M_{\star}$ relation, we see a progressive increase of the IRX with the redshift. This trend is difficult to validate/invalidate with current observational measurements.

6 Conclusion

We presented our implementation of dust extinction processes in the G.A.S. galaxy evolution model. The effective attenuation is assumed to be caused by three types of dust grain (PAH, VSG and BG). Each galaxy in our model can host up to three components: i) an homogeneous and diffuse ISM hosting old stars, ii) some GMCs composed of fragmented/dense gas in which the youngest stars are still embedded and iii) a central bulge component hosting old stars transferred during major mergers. The proportion of the different grains is driven by the PAH mass fraction which is linked to the metallicity of each gas phase. The metallicity is computed separately in each phase according to their metal content ($m_{O}$) and following the oxygen to hydrogen abundance. These different metallicities lead to different dust compositions in the different components. In each gas phase, the evolution of the dust composition with the gas metallicity allows to compute the attenuation curves that can be shallower or steeper (SMC-like) than the extensively used Calzetti et al. (2000) attenuation law. Even if the disk inclination affects the slope of the attenuation curve, the absolute attenuation in the FUV band appears to be the main driver.

For each simulated galaxy, we built the intrinsic and the attenuated stellar spectrum. In each galaxy environment, the intrinsic stellar spectrum is computed following the explicit star formation history and the metal enrichment. We then produced intrinsic and extinguished UV luminosity functions, from $z\simeq 9$ to $z=0.1$ and deduced the evolution of the absorbed fraction of UV emission. The predictions of our model are in very good agreement with measurements. We confirm that at least half of the UV radiation is absorbed and reprocessed by dust at z $\sim$ 8.0. Metal enrichment is therefore already very active at this epoch.

By assuming that all the energy absorbed in the UV-optical wavelength range is reprocessed in the IR, we computed the IR luminosity of each galaxy. The IR luminosity function predicted by our model are in good agreement with existing observational measurements (between z = 4.0 and z = 0.3). At z = 0.1 and z = 0.4 the model slightly overestimates the total IR luminosity function, especially in the bright-end. At higher z, our model predicts that the most massive and most active star-forming galaxies can reach IR luminosities close to $10^{12.5}$, even at $z\simeq 7$.

The analysis of the IR luminosity function reveals that the main part of the IR radiation comes from the strongly-attenuated GMCs, more than 80% at $z>3$ and at least 50% at $z=1.0$. The contribution coming from GMCs in the total IR emission is slightly higher in low mass/luminous IR galaxies, in which the ISM is more structured/fragmented, than in high mass/luminous IR galaxies.

We then focused on the InfraRed Excess. For z$<$3.0, the IRX–$M_{\star}$ relation predicted by our model agrees with the recent measurements (e.g. Pannella et al. 2015; Heinis et al. 2014; Reddy et al. 2010). This supports previous suggestions indicating that the stellar mass may be used as a rough proxy for dust attenuation. The choice of the $Z_{g}-f_{PAH}$ relation does not impact this result. At $z>3.0$, our model predicts an increase of the average IRX for a given stellar mass, although this trend is not observed (Pannella et al. (2015); Heinis et al. (2014); Whitaker et al. (2014)).

The IRX-$\beta$ correlation is one of the key tools to correct the UV emission of star forming galaxies from dust attenuation. This correlation is widely explored and discussed in the literature and large scatters and departures are found. Consistently with the observations, the dispersion of the IRX-$\beta$ relation predicted by our model is large. At a given UV slope, some galaxies hosting a strong star-formation activity ($>100M_{\odot/yr}$) can reach IRX values 100 times larger than those predicted by the Meurer et al. (1999) fiducial relation.

With our model, we confirm that galaxies hosting an “old” ($>$ 50Myr) stellar population lie preferentially below the Meurer et al. (1999) relation. An attenuation law steeper than the Calzetti et al. (2000) average attenuation law also shifts galaxies below the Meurer et al. (1999) relation. The deviations from the standard relation are essentially driven by these two factors.

This analysis confirms that $\beta$ is not a good tracer of dust attenuation, especially in galaxies hosting a strong star formation activity. The choice of the attenuation curve used to correct for attenuation is therefore crucial and has to be taken with caution.

Acknowledgements.

MC thanks Laurent Verstraete for very useful discussions about the DustEM model and Mederic Boquien for discussions about extinction prescriptions. MC thanks the Centre National d’Etude Spatial for its financial support. We acknowledge financial support from the “Programme National de Cosmologie and Galaxies” (PNCG) funded by CNRS/INSU-IN2P3-INP, CEA and CNES, France, from the ANR under the contract ANR-15-CE31-0017 and from the OCEVU Labex (ANR-11-LABX-0060) and the A*MIDEX project (ANR-11-IDEX-0001-02) funded by the “Investissements d’Avenir” French government programme managed by the ANR.

Appendix A PAHs

A.1 Impact of the $Z_{g}-f_{PAH}$ relation

The $Z_{g}-f_{PAH}$ relation is a strong assumption for the dust attenuation model implemented in G.A.S.. To illustrate its impact on the effective attenuation, we created two modified relations, mostly compatible with the lowest and the highest $f_{PAH}$ values associated with the observed sample Draine et al. (2007); Ciesla et al. (2014); Rémy-Ruyer et al. (2015). The two relations are defined following:

[TABLE]

[TABLE]

with $f_{PAH,0}=4.57\%$, which is the reference value for our Galaxy.

In Fig. 16 we compare the two modified relations with the original one. According to Eq. 21 and Eq. 22 the amount of PAH is respectively reduced and increased in the gas phases. In the two cases, the dust composition of low-metallicity galaxies ($Z_{g}<$ 8.0) is strongly affected.

To highlight the impact of these new $Z_{g}-f_{PAH}$ relations, we compare in Fig. 17 the effective attenuation (Eq. 19) produced by both the original ($A_{\rm{Orig}}$) and the modified ($A_{\rm{Modif}}$) relations. The analysis is performed on a sub-sample of galaxies, randomly selected at different redshifts. Attenuations are computed for both the FUV and the NUV bands. For the two bands and the two modified relations, the impact is small (only up to 0.15 Mag) and mainly affects bright galaxies. As expected, a smaller(larger) fraction of PAH leads to galaxies slightly more(less) luminous. We note that the impact is higher in the NUV than in the FUV band. This trend is fully expected. Indeed, as presented in Sects. 3.4.3 and 3.4.2 the effective attenuation in the FUV band is fixed according to both the hydrogen column density and the gas phase metallicity, following Boquien et al. (2013). The modification of the $Z_{g}-f_{PAH}$ does not affect this normalisation. The fraction of PAH only slightly affects the slope of the extinction curve in the FUV band. The measured variations ($<$0.05 Mag) are only due to this small slope change integrated through the FUV filter response. In the galaxy rest frame, the NUV band falls on the PAH bump. This bump can strongly modify the local slope of the extinction curves around the NUV band (Fig. 2(b)). The impact of the new $Z_{g}-f_{PAH}$ relations are therefore larger in the NUV than in the FUV band. We also note that the dispersion around the mean trend is also larger in the NUV.

The overall impact of these modified $Z_{g}-f_{PAH}$ relations can be measured through the IRX variation, as shown in Fig. 18. The variation is plotted as a function of the intrinsic FUV magnitude that does not evolve between the two relations. The impact is larger when the PAH fraction increases (up to 4%) than when the PAH fraction decreases (only -0.8%). By reducing(increasing) the effective amount of PAH, IRX slightly decreases(increases). As for the effective attenuation in the FUV and NUV band, higher is the redshift, smaller is the impact. In addition, UV bright galaxies are more affected than UV faint galaxies. These two trends are linked to the progressive increase of the average metallicity, from high to low redshift and from low- to high-stellar masses.

A.2 PAH fraction: random walk algorithm

To allow for a smooth and continuous evolution of the PAH, VSG and BG fractions, from one time-step to the next one, the fraction of PAH evolves according to a random walk between the previous value, $f_{PAH}^{n-1}$, and a target value, $f_{PAH}^{{\dagger}}$, expected from the empirical PAH-metallicity relation (see Eq. 3):

[TABLE]

where $r$ is a random number following a uniform distribution. For the BG component, we apply a similar algorithm. The target mass fraction is settled $f_{BG}^{{\dagger}}=\frac{2}{3}(1-f_{PAH}^{n})$. Then we apply:

[TABLE]

The mass conservation rule then requires that VSG are distributed following the residual mass.

[TABLE]

Appendix B The GALAKSIENN library

The GALAKSIENN library stores the main results produced by our new G.A.S.  semi-analytical model, especially MOCK galaxy catalogs and sky maps described in our paper III. It is available online through the ZENODO platform: https://zenodo.org/. A complete description of the GALAKSIENN library is given in paper III. In association with this paper II, we distribute the ASCII tables of the FUV and IR luminosity functions presented in Fig. 5 and Fig. 6 respectively.