<figure><img src="image_1.png" /></figure>

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

Morgane Cousin, V. Buat, Guilaine Lagache, M. Bethermin

► To cite this version:

Morgane Cousin, V. Buat, Guilaine Lagache, M. Bethermin. G.A.S. II: Dust extinction in galaxies; Luminosity functions and InfraRed eXcess. Astronomy and Astrophysics - A&A, 2019, 10.1051/0004-6361/201834674. hal-02008601

HAL Id: hal-02008601

https://hal.science/hal-02008601

Submitted on 29 Nov 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
G.A.S.

II. Dust extinction in galaxies: Luminosity functions and infrared excess

M. Cousin, V. Buat, G. Lagache, and M. Bethermin

Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France
e-mail: morgane.cousin86@gmail.com

Received 16 November 2018 / Accepted 17 January 2019

ABSTRACT

*Context.* 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.

*Aims.* We present the complete prescription of dust attenuation implemented in the new semi-analytical model called 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. Paper I).

*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 ultraviolet (UV) and the infrared (IR) luminosity functions from $z = 9.0$ to $z = 0.1$. Our two-phase ISM prescription catches the evolution of UV and IR luminosity functions very well. 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 far-ultraviolet (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.

**Key words.** galaxies: formation – galaxies: evolution – galaxies: fundamental parameters – ultraviolet: galaxies – infrared: galaxies – methods: numerical

**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 ultraviolet (UV) and visible and near-infrared (nIR) light, dust, heated by stellar radiation, radiates from the nIR to the submillimeter (submm). Attenuation of the stellar light due to dust dramatically modifies the fluxes, colors and spectral energy distributions (SEDs) 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 IR and the submm domain. This process is particularly efficient in dusty star-forming galaxies (DSFGs), which 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 these 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 as follows: $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). Because of 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 the corrections for this attenuation accurately is thus crucial to derive SFRs. One of the most popular methods 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 exten-
sively 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 complex star
formation histories. Indeed, works by Cortese et al. (2006),
Boquien et al. (2009), or more recently Pannella et al. (2015)
have revealed a correlation between the dust attenuation and gas-
phase metallicity, which is itself correlated with stellar mass (e.g.
Mannucci et al. 2010; Lilly et al. 2013; Zahid et al. 2014).

Understanding dust attenuation in galaxies is still an impor-
tant challenge and therefore a fundamental aspect of SAMs. The
impact of dust has been studied in SAMs with various 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) leading to number counts in
UV and/or luminosity functions that differ by a factor ≈2. The
recent work of Narayanan et al. (2018), based on hydrodynamic
zoom-in simulations and radiative transfer and the analytical pre-
scription 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 \approx 5$ to $z \approx 2$, their limited
number of simulated galaxies cannot 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 stel-
lar light. Physical prescriptions associated with the gas physics
is described in Cousin et al. (2019).

We investigate in this paper the redshift evolution of the far-
ultraviolet (FUV) and IR 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 \approx 1.5$. We explore
the impact of the stellar age, SFR, the absolute attenuation, disk
inclination and gas-phase metallicity. We 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 com-
ponents. 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 attenua-
tion law predicted by our model as a function of 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 and UV slope. We con-
clude in Sect. 6.

## 2. G.A.S. model

Our work is the continuation of the G.A.S. semi-analytical
model (SAM) described in Cousin et al. (2019). G.A.S. is
applied to dark-matter merger trees extracted from a pure N-
body simulation. This simulation uses a WMAP-5 yr cosmol-
ogy ($\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 \times 10^8$ M$_\odot$. Dark-matter halos and substructures (satel-
lites) are identified via the HaloMaker code (Tweed et al. 2009).
In merger tree structures, we only take into account halos with at
least 20 dark-matter particles, leading to a minimal dark-matter
halo mass of $2.050 \times 10^9$ 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-matter 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. Paper I, Sect. 2). We assume that the warm con-
densed 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 frag-
mentation following a turbulent cascade and is finally converted
in star-forming gas. Galaxies modeled by G.A.S. host therefore
two distinct ISM phases: a diffuse warm phase and a fragmented
dense cold phase (G.A.S. Paper I, Sects. 3 and 4). Young stars
are born in dense, structured gas. We assume that they stay in
this dense phase during $5 \times 10^7$ yr. Older stars then evolve in the
diffuse warm phase of the disk.

Our model takes into account supernovae (SNs) and active
galactic nuclei (AGN) feedback. A part of the power pro-
duced by these two types of feedback is directly injected in
the ISM of the galaxy and is used to disrupt the fragmented
dense gas. The remaining power generates large-scale ejecta
and heats 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. Paper I,
Sect. 5).

In addition to these processes, we introduce a regulation
mechanism operating in the surrounding hot gas phase which
reduces gas cooling in massive halos. It is based on thermal
instabilities (G.A.S. Paper I, 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 \approx 6$ to $z \approx 0.5$.

### 2.2. 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 ele-
ments, $^1$H, $^4$He, $^{12}$C, $^{14}$N, $^{16}$O and $^{56}$Fe. Masses are tracked sep-
arately: on one side in the diffuse and the fragmented/dense gas
phase of the disk, on the other side in the diffuse gas phase of
the central bulge (if any). The production and the reinjection
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 dense/fragmented gas are strewn in the hot atmosphere surrounding the galaxy with large-scale ejecta. Owing to this metal enrichment, the efficiency of the radiative cooling process affecting the hot surrounding gas is modified and impacts 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,

$$Z_g = 12 + \log(\text{O/H}) = 12 + \log_{10} \left( \frac{M_O m_H}{M_H m_O} \right), \quad (1)$$

where $M_O$ and $M_H$ are the oxygen and the hydrogen mass contained in a given gas phase and $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

Similar to the gas metallicity, the G.A.S. model tracks the evolution of the 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 with 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:

$$\eta_m = \frac{\text{MIN}(M_{1/2,1}; M_{1/2,2})}{\text{MAX}(M_{1/2,1}; M_{1/2,2})}, \quad (2)$$

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 \ge 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 ($\le 5 \times 10^7$ yr) of the two progenitor disks are still embedded in giant molecular clouds (GMCs). Young stars are therefore kept into the remnant disk.

Conversely, old stars (>5 $\times 10^7$ yr) 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 on 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.

<table><caption>**Table 1.** Dust properties.</caption><thead><tr><th>Type</th><th>ρ (g cm<sup>-3</sup>)</th><th>Size range (cm)</th><th>Size law</th></tr></thead><tbody><tr><td>PAH0</td><td>2.24</td><td>3.1 × 10<sup>-8</sup>–1.2 × 10<sup>-7</sup></td><td>mix-log n</td></tr><tr><td>PAH1</td><td>2.24</td><td>3.1 × 10<sup>-8</sup>–1.2 × 10<sup>-7</sup></td><td>mix-log n</td></tr><tr><td>amCBE</td><td>1.81</td><td>6.0 × 10<sup>-8</sup>–2.0 × 10<sup>-6</sup></td><td>log n</td></tr><tr><td>amCBE</td><td>1.81</td><td>4.0 × 10<sup>-7</sup>–2.0 × 10<sup>-4</sup></td><td>plaw-ed</td></tr><tr><td>aSil</td><td>3.00</td><td>4.0 × 10<sup>-7</sup>–2.0 × 10<sup>-4</sup></td><td>plaw-ed</td></tr></tbody></table>

Notes. The dust type and size law refer to DustEM keywords.

# 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; and (iii) the relative geometry of the dust and stars.

## 3.1. Dust composition

As proposed by for example Draine & Li (2001, 2007), Zubko et al. (2004) and 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 ISM, gas and dust are closely linked and the dust composition (PAH: $f_{PAH}$, BG: $f_{BG}$, VSG: $f_{VSG}$) depends on the gas metallicity.

Figure 1 shows measurements of the PAH fractions as functions 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 using Eq. (5) in Rémy-Ruyer et al. (2015) written as

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

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

Equation (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
<figure><img src="image_3.png" /><figcaption>Fig. 1. Fraction of PAH as a function of the gas metallicity in the dense/fragmented gas phase. The blue shaded area shows the full G.A.S. distribution. Blue circles and light-green squares show Ciesla et al. (2014) and Draine et al. (2007) data points respectively. Green triangles show Rémy-Ruyer et al. (2015) sample shifted by 0.2 dex. The black solid line corresponds to our mean relation (Eq. (3)). The original mean relation of Rémy-Ruyer et al. (2015; their Eq. (5)) is shown with the solid gray line. Dashed black lines mark the 3σ limits around our mean relation.</figcaption></figure>

Appendix A.2). We plot in Fig. 1 the whole distribution of $f_{\text{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.

Big grains 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 two-thirds 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 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 disk and bulge and dense/fragmented gas in disk), 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_{\text{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 detail by G.A.S. Indeed, 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_{\text{PAH}}$ relation assumed in G.A.S. (Eq. (3)) follows 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 on the dust attenuation.

### 3.3. Extinction curves

Once the dust composition ($f_{\text{PAH}}$, $f_{\text{BG}}$, $f_{\text{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).

Figure 2a shows the absorption and the scattering curves associated with the different dust types.

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

$$ \tau_{\lambda} = \sum_{\text{gt}} f_{\text{gt}} \times (\sigma_{\text{abs,gt}} + \sigma_{\text{sca,gt}}). \qquad (4) $$

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

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

$$ \omega_{\lambda} = \frac{\sum_{\text{gt}} f_{\text{gt}} \times \sigma_{\text{sca,gt}}}{\tau_{\lambda}}. \qquad (5) $$

Figure 2b 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\text{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. 2b, both the extinction ($\tau_{\lambda}$) and 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})$ which allows us to build the effective attenuation

$$ A_{\lambda} = -2.5 \log_{10} \Phi_{\lambda}(\tau_{\lambda}, \omega_{\lambda}) \qquad (6) $$

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 with 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
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Panel a: absorption ($\sigma_{abs}$, upper panel) and scattering ($\sigma_{sca}$ cross sections, lower panel). In each panel, PAH, VSG and BG are shown in blue, light blue and green, respectively. The gray curve shows the total absorption. We assume $f_{PAH} = f_{PAH,0} = 4.57\%$. The mass fraction of BG and VSG (given in the left panel) are then deduced using the random-walk technique. Panel b: extinction ($\tau_{\lambda}$) curves normalized to the FUV band. We compare those extinction curves to observational measurements: Mathis et al. (1983) Milky Way (dashed black line) and Gordon et al. (2003) for SMC (gold solid line), LMC (orange solid line) and Milky Way (red solid line). Lower panel: Albedo ($\omega_{\lambda}$). In the two panels, the color of the curves is linked to the PAH fraction assumed, from $10^{-2} \times f_{PAH,0}$ (blue) to $10 \times f_{PAH,0}$ (green).</figcaption></figure>

(Guiderdoni & Rocca-Volmerange 1987). However, the analysis done in Boquien et al. (2013; their Fig. 2) indicates that attenuation, measured in the FUV band, is more reliable than in the V band. These authors show that the relative scatter of the attenuation in the FUV band is 30%–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 × 10⁷ yr) and N_GMC GMCs hosting the youngest stars.

Three different geometrical functions are associated with these three different environments. For the bulge component we adopt a Dwek geometry (Sect. 3.4.2). We use the classical S1ab (Sect. 3.4.3) geometry for the diffuse ISM. In addition to this

S1ab 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 only affect the stellar radiation of the component if the gas metallicity is sufficiently high at 12 + log(O/H) > 6.5; if that is not the case, we neglect the dust attenuation.

### 3.4.2. Bulge environments: 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 generates 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 et al. 1996; Devriendt et al. 1999)

$$ \Phi_{\lambda}^{\text{SPH}} = \frac{a_{\lambda}}{1 - \omega_{\lambda} + \omega_{\lambda} a_{\lambda}}, \qquad (7) $$

where

$$ 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], \qquad (8) $$
and

$$k_{\lambda} = \tau_{\lambda}^{\text{SPH}} = \tau_{V}^{\text{SPH}} \times \frac{\tau_{\lambda}}{\tau_{\lambda}(\lambda_{V})}. \qquad (9)$$

The reference depth in V band, $\tau_V^{\text{SPH}}$, varies with the metal mass fraction<fnref n="1" /> $\chi_Z$ and the hydrogen column density $N_H$ following (Guiderdoni & Rocca-Volmerange 1987; Devriendt et al. 1999):

$$\tau_V^{\text{SPH}} = 2.619 \left(\frac{\chi_Z}{\chi_{\odot}}\right)^{1.6} \left(\frac{N_H}{2.1 \times 10^{21} \text{ at cm}^{-2}}\right). \qquad (10)$$

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: S1ab geometry

To describe the dust attenuation in the diffuse warm ISM of each galaxy disk we use the standard S1ab 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

$$\Phi_{\lambda}^{\text{ISM}} = \frac{1 - \exp\left(-\sqrt{1 - \omega_{\lambda}^{\text{ISM}} \tau_{\lambda}^{\text{ISM}} / \cos i}\right)}{\sqrt{1 - \omega_{\lambda}^{\text{ISM}} \tau_{\lambda}^{\text{ISM}} / \cos i}} \qquad (11)$$

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

$$\tau_{\lambda}^{\text{ISM}} = \tau_{\text{FUV}}^{\text{ISM}} \times \frac{\tau_{\lambda}}{\tau_{\lambda}(\lambda_{\text{FUV}})}. \qquad (12)$$

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 on the diffuse gas-phase metallicity $Z_g^{\text{ISM}}$. The face on UV depth $\tau_{\text{FUV}}^{\text{ISM}}$ is calculated as in Eq. (13) from Boquien et al. (2013) as follows:

$$\tau_{\text{FUV}}^{\text{ISM}} = (1.926 + 0.051 \times \Sigma_{\text{H}}^{\text{ISM}}) \times 10^{[0.947 \times Z_{\text{g}}^{\text{ISM}} - 9]}. \qquad (13)$$

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

### 3.4.4. Additional attenuation from GMCs

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

$$\Phi_{\lambda}^{\text{GMC}} = \exp\left(-\sqrt{1 - \omega_{\lambda} \tau_{\lambda}^{\text{GMC}}}\right). \qquad (14)$$

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

$$\tau_{\lambda}^{\text{GMC}} = \tau_{\text{FUV}}^{\text{GMC}} \times \frac{\tau_{\lambda}}{\tau_{\lambda}(\lambda_{\text{FUV}})}. \qquad (15)$$

The normalization factor $\tau_{\text{FUV}}^{\text{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_{\text{GMC}}$ spherical clouds of radius $R_{\text{GMC}} = 3h/8$ (Cousin et al. 2019). Normalization also depends on the structured gas metallicity $Z_g^{\text{GMC}}$. Metallicity in the dense phase is slightly higher than in the diffuse medium.

### 3.5. Attenuated stellar spectra

A galaxy contains up to three distinct stellar populations associated with three different geometrical distributions. We build up the three SEDs associated with each stellar population:

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

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

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

The overall extinguished stellar spectrum is then given by

$$(\lambda I_{\lambda}) = \Phi_{\lambda}^{\text{ISM}} \Phi_{\lambda}^{\text{GMC}} (\lambda I_{\lambda})_{\text{GMC}} \qquad (16)$$

$$+ \Phi_{\lambda}^{\text{ISM}} (\lambda I_{\lambda})_{\text{ISM}} \qquad (17)$$

$$+ \Phi_{\lambda}^{\text{SPH}} (\lambda I_{\lambda})_{\text{SPH}}. \qquad (18)$$

### 3.6. Infrared energy budget

The bolometric luminosity produced by stars is:

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

where $(I_{\lambda})_{\star}$ is the unextinguished 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:

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

where $(I_{\lambda})$ is the extinguished stellar SED built through 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 as

$$\Phi_{\lambda}^{\text{eff}} = \frac{(\lambda I_{\lambda})}{(\lambda I_{\lambda})_{\star}} \quad \text{and} \quad A_{\lambda}^{\text{eff}} = -2.5 \log_{10} \Phi_{\lambda}^{\text{eff}}. \qquad (21)$$

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. Figure 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.1-1]$ µ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]$ µ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

<footnote n="1">¹ We assume a solar metal mass fraction $\chi_{\odot} = 0.02$.</footnote>
<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Median effective attenuation curves as functions of the disk inclination for a random star-forming galaxy sample extracted between $z = 4$ and $z = 1.5$. The galaxy sample has been divided in 9 bins of disk inclination from face-on (0°) to edge-on disks (90°). The color code indicates different disk inclinations. The median trend of the full sample is plotted with the solid gray line. We compare this trend with a set of standard attenuation curves: Calzetti et al. (2000): circles, Reddy et al. (2016): squares and Gordon et al. (2003; SMC): triangles.</figcaption></figure>

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 that resulting from the variation of inclination. To illustrate this scatter we select galaxies with inclinations between 30° and 50°. In this subsample, 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 subsample 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

<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Variation of the median effective attenuation curves with $A_{FUV}$ (color coded). We only select here galaxies with a disk inclination between 30° and 50°. The median trend of the full sample is plotted with the solid gray line. We compare this trend to a set of standard attenuation curves, Calzetti et al. (2000): circles, Reddy et al. (2016): squares and Gordon et al. (2003; SMC): triangles.</figcaption></figure>

functions are built for stellar masses greater than $10^{7} M_{\odot}$. This mass threshold corresponds to our stellar mass resolution. Close to and below this mass, our predictions are incomplete and suffer from numerical instabilities. Particularly, gas-phase metallicity cannot be computed and dust composition cannot be defined. According to 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 to the incompleteness of the data.

## 4.1. Ultraviolet luminosity functions

### 4.1.1. Overall evolution

Figure 5 shows the predictions for the UVLF at different redshifts. Ultraviolet magnitudes were computed by convoluting the rest-frame spectrum of each galaxy with the GALEX (Martin et al. 2005a) FUV transfer function<fnref n="2" />. We compare our predictions with a large set of observational data coming from Arnouts et al. (2005), Oesch et al. (2010, 2013), Cucciati et al. (2012), McLure 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 \approx 9.0$ and at $z = 3.0$ our predictions are slightly below the observations, especially in the bright end of the luminosity function at $z = 3.0$ and in the low-luminosity range at $z \approx 9.0$. At $z \approx 9.0$ the number of galaxies taken into account in the predicted UVLF is low (<500). At these very high redshifts, mea-

<footnote n="2">² Obviously, similar results are obtained by convoluting, for example, *i*, *r*, or *u* filter bands with the redshifted galaxy spectrum.</footnote>
<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Ultraviolet luminosity functions. The different panels correspond to different redshift ranges. Predictions of our model are compared to a set of observations: Arnouts et al. (2005): dark blue circles; Oesch et al. (2010): dark green up-triangles; Cucciati et al. (2012): blue circles; McLure et al. (2013): orange down-triangles; Oesch et al. (2013) yellow up-triangles; Alavi et al. (2014): green low-triangles; Parsa et al. (2016): dark green diamonds; Hagen et al. (2015): grey low-triangles; Bouwens et al. (2015) lightblue pentagons; Finkelstein et al. (2015): light green squares; Atek et al. (2015): light-green up-triangles; and Livermore et al. (2017): purple pentagons. In each panel $f_{\text{ext}}$ indicates the fraction of the intrinsic UV emission of stars, which is attenuated by dust. The tick marks on the *left side* of each panel allow to us link the comoving density to the effective number of galaxies taken into account in the luminosity function.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Infrared luminosity functions. The different panels are associated with different redshift ranges. G.A.S. model predictions are only plotted above the stellar mass completeness limit ($10^7 M_{\odot}$). We compare our results with measurements done by Caputi et al. (2007): green squares; Rodighiero et al. (2010): light green diamonds; Gruppioni et al. (2013): blue circles; Magnelli et al. (2013): light blue pentagons; Marchetti et al. (2016): purple squares; and Koprowski et al. (2017): dark blue triangles. The solid gray line gives the total IRLF (including diffuse disk ISM, GMCs and bulge). The gray dashed line gives the IRLF that takes into account only the GMC environment and therefore only the young stars contribution. The tick marks on the *left side* of each panel allow us to link the comoving density to the effective number of galaxies taken into account in the luminosity function.</figcaption></figure>
surements are based on only few galaxies and may suffer from strong uncertainties. Therefore UVLFs are affected by statistical noise. At $z \approx 3.0$, UVLFs are built with a larger number of objects: the discrepancy between observational measurements and prediction cannot only be explained by uncertainties. At $z \approx 3.0$ specifically, modeled 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 2019). The IRLF (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 \approx 3.0$ therefore cannot be currently explained.

At all redshifts and for the range of magnitudes considered in this work, 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 UV luminosity. At high redshifts ($z \ge 6.0$), the shape of the extinguished UVLF is still compatible with a power law but with a steeper slope than the unextinguished 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 much 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 separately sum the extinguished and intrinsic UV luminosities over all galaxies. We define the UV extinguished fraction as the ratio of these two quantities

$$f_{\text{ext}} = 1.0 - \frac{\sum L_{\text{FUV}}^{\text{ext}}}{\sum L_{\text{FUV}}^{\text{int}}}. \qquad (22)$$

This fraction is indicated in each panel of Fig. 5. Our model predicts an UV extinguished fraction of 40% at $z \approx 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 \approx 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 \approx 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 \le 2.0$), the fraction of UV radiation absorbed by dust still increase but much more slowly. At $z \approx 0.3$ the fraction of UV radiation absorbed by dust reaches a value of $\approx 86\%$.

By integrating a compilation of UV and IR luminosity functions, Burgarella et al. (2013) found a similar ratio of 85 ± 10% and 84 ± 13% at $z = 2.2$ and $z = 3.15$, respectively. These results agree 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 ($f_{\text{ext}} = 50\%$) is reached at $z \approx 6$. Several teams also tried using various modeling approaches to estimate the evolution of the ratio of obscured to 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 cosmic infrared background anisotropies from Maniyar et al. (2018). In contrast, this value disagrees with Koprowski et al. (2017) who predict a transition around $z \approx 3$ based on their empirical modeling approach.

## 4.2. Infrared luminosity functions

### 4.2.1. The evolution of the total IRLF

Figure 6 shows our predictions for the total IRLF 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 in this work, the total IRLF 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_{\text{IR}} \approx 10^{12} L_{\odot}$ at $z \approx 6.0$. The break progressively shift to lower luminosity. At $z \le 1.0$ the knee appears around $L_{\text{IR}} \approx 10^{11.5} L_{\odot}$.

For the redshifts $z \approx 2.8$ and $\approx 3.7$, even if the bright ends are in good agreement with the recent Koprowski et al. (2017) estimates, G.A.S. seems to overpredict the density of IR galaxies for $L_{\text{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 \approx 7$. Metal production and dust assembly are 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 from environments, Possible cold dust at low redshift?

In addition to the total IRLF, 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_{\text{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 $\approx 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.
<figure><img src="image_10.png" /><figcaption>**Fig. 7.** IRX–$A_{FUV}$ relation derived from a sample of star-forming galaxies extracted at $z = 1.5$. The gray solid line indicates our best fit. For comparison we also plot the Buat et al. (2005) best fit (black solid line). Colored points show a subsample that is uniformly selected in the IRX–$A_{FUV}$ plan. The color scale shows the luminosity-weighted age of the stellar population hosted by the galaxy.</figcaption></figure>

At $z \le 0.6$ the fraction of IR luminosity coming from the GMC environment drops below 50% and even reaches 39% at $z \approx 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 overdensity of IR galaxies. The IR emission is slightly too high for the full IR luminosity range explored. However, if we focus only on 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. Paper III (Cousin 2019). Infrared 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. Impact of the latest disruption episode on massive IR bright galaxies

At all redshifts explored in this work, the gap existing between the total and GMC-only IRLF seems to indicate that the contribution of the GMC to the total IR luminosity seems lower at the bright end than at the faint end of the IRLF. 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, 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 \le 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} \ge 10^{10.5} M_{\odot}$, $L_{IR} \ge 10^{10.5} L_{\odot}$) the gas content

<figure><img src="image_11.png" /><figcaption>**Fig. 8.** IRX–$Z_g$ relation for star-forming main sequence galaxies extracted at $z = 1.5$. The gas-phase metallicity is computed as the oxygen over hydrogen abundance $Z_g = 12 + \log(O/H)$. We focus on the dense/fragmented gas-phase metallicity. The blue density shows the complete star-forming galaxy sample. The solid gray line shows the median trend of this sample. The dashed gray lines indicate the 15% and 85% quantiles of the distribution. We compare the G.A.S. model with observational measurements coming from Pannella et al. (2015) and Cortese et al. (2006).</figcaption></figure>

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. Paper I). 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) than in the less massive galaxies.

## 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 such as the slope of the UV spectrum, stellar mass or SFR (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 SFR (sSFR) is greater$^3$ than $1/4 \times (SFR_{MS}/M_{\star})$.

<footnote>$^3$ We checked that our results are not significantly affected if we use a threshold of one-third or one-fifth instead.</footnote>
<figure><img src="image_12.png" /></figure>

**Fig. 9.** IRX-$M_\star$ relation for star-forming main sequence galaxies extracted at $z = 1.5$. The blue density shows the complete sample. The dashed gray lines show the median relation. The solid gray line shows the mean relation computed using $\text{IRX} = \langle L_{IR} \rangle/\langle L_{UV} \rangle$, to be consistent with what is measured by stacking. Around the median relation, the dotted gray lines indicate the 15% and 85% quantiles of the distribution. We compare our predictions with observational measurements from Pannella et al. (2015; orange pentagons), Heinis et al. (2014; yellow squares) and Whitaker et al. (2014; solid red line).

## 5.1. IRX-$A_{FUV}$ 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 using radiative transfer models (e.g. Witt et al. 2006).

Figure 7 shows the IRX-$A_{FUV}$ relation based on a sample star-forming galaxies sample extracted at $z = 1.5$. Our prediction 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 $A_{FUV}$. 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 $A_{FUV}$ is associated with an older stellar population. Similarly, at a given $A_{FUV}$, 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-$A_{FUV}$ relation with the redshift.

## 5.2. 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 Boquien et al. (2009), for example. In their sample, galaxies with higher metallicities have higher IRX.

<figure><img src="image_13.png" /></figure>

**Fig. 10.** IRX-$M_\star$ relation for main sequence galaxies at $z = 1.5$. We plot the 15% and 85% percentiles of the distribution (dashed gray lines). Colored points correspond to a random sample uniformly distributed in the IRX-$M_\star$ plan, with a color scale linked to the extinguished UV luminosity of the galaxies.

Figure 8 shows the IRX-$Z_g$ relation built from our $z = 1.5$ star-forming galaxy sample. As discussed 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 discussed 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$). As a consequence of 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. IRX-$M_\star$ relation

Our median IRX-$M_\star$ relation is plotted in 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.
<figure><img src="image_14.png" /></figure>

**Fig. 11.** Variation of the IRX-$M_\star$ relation for main sequence galaxies from $z = 1.5$ to $z \approx 6.0$. The colored solid lines show the mean (obtained using IRX = $\langle L_{IR} \rangle / \langle L_{UV} \rangle$) relation. Our predictions are compared to measurements coming from Whitaker et al. (2014) at $z = 1.75$ and $z = 2.25$. We also add average relations from Heinis et al. (2014), Pannella et al. (2015) and McLure et al. (2018) which are available from $1.5 < z < 4.0$. The redshift of each sample and the associated dispersion are indicated in the bottom right corner.

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, i.e., $\text{IRX} = \langle L_{IR} \rangle / \langle L_{UV} \rangle$ (stacking estimator). We compute this estimator in our various 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, the 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).

<figure><img src="image_15.png" /></figure>

**Fig. 12.** IRX-$\beta$ relation for main sequence galaxies at $z = 1.5$. The median of the full sample (in blue) is shown with the dashed gray line. The IRX mean is shown with the solid gray line. Observations from Meurer et al. (1999), Overzier et al. (2011), Takeuchi et al. (2012) and, Heinis et al. (2013) are shown in solid red, dot-dashed orange, dotted green and dashed yellow line, respectively. We also add contours coming from measurements performed in the COSMOS field by Casey 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. Owing to this homogeneous selection, the dispersion of these points is much larger than the scatter of the relation given above. This random subsample 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 Figs. 6, C.1 and C.2).

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 to 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

Figure 11 shows the evolution of the IRX-$M_\star$ relation with redshift (from $z = 1.5$ to $z = 5.9$) obtained using our stacking estimator. At all redshifts explored, the scatter is similar and close to ~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
<figure><img src="image_16.png" /><figcaption>**Fig. 13.** IRX-β relation for main sequence galaxies at z = 1.5. Colored points show a random sample of star-forming galaxy uniformly distributed into the IRX-β plan. The color code is linked to the luminosity-weighted stellar population age.</figcaption></figure>

shifts to lower mass when the redshift is increasing (from $\approx 10^{10.5}$ to $\approx 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) which has a trend similar to that seen in our model. As previously (Sects. 4.2, 5.2), we explain this turnover 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 measured<fnref n="4" /> 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 at z > 3.0. At lower redshift the large dispersion ($\approx 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 must also be reduced 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. IRX-β relation

Figure 12 shows the IRX-β relation based on our star-forming galaxies sample extracted at z = 1.5. To compute the UV slope ($\beta$) we 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 pho-

<figure><img src="image_17.png" /><figcaption>**Fig. 14.** IRX-β relation for main sequence galaxies at z = 1.5 with colored points marking the values of the slope of the attenuation curve. The right-hand axis indicates the normalized slope of the attenuation curve. The parameter $\Sigma$ is the ratio of the effective attenuation curve slope and the slope of the Calzetti et al. (2000) attenuation law (here -0.54 dex $\mu$m). In other words, $\log_{10} \Sigma > 0$ (resp. $\log_{10} \Sigma < 0$) indicates a steeper (resp. flatter) attenuation curve than Calzetti et al. (2000) attenuation curve. We compare with average IRX-β computed by Salmon et al. (2016) for attenuation curves estimated by Gordon et al. (2003; SMC), Calzetti et al. (2000) and Reddy et al. (2015). These three average relations are plotted with solid colored lines. As for the random galaxies sample, the color code is linked to the normalized slope of the attenuation curve.</figcaption></figure>

tometric points associated with filters 1, 3, 5, 7 and 9 listed in Calzetti et al. (1994; their Table 2). Similar to what we did for the IRX-$M_{\star}$ 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^{-1}$) 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. We note that this bias also affects 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 $\log_{10}$ IRX $\in [0.0, 2.0]$, some galaxies can have $\log_{10}$ IRX > 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 on the IRX-β relation has been originally studied by Kong et al. (2004) and then discussed and confirmed in 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

<footnote n="4">⁴ These studies use different observational data: Spitzer for Whitaker et al. (2014), Herschel for Heinis et al. (2014) and Pannella et al. (2015).</footnote>
<figure><img src="image_18.png" /><figcaption>**Fig. 15.** Variation of the IRX-β relation with redshift, from $z = 1.5$ to $z \approx 6.0$. For each redshift, medians trends are shown with thick colored lines. Observational trends coming from Heinis et al. (2013) and Meurer et al. (1999) are shown with dashed and solid black lines, respectively. The Colored points (according to the redshift color coding) indicate the recent measurements performed by Koprowski et al. (2018) with ALMA.</figcaption></figure>

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-β 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 subsample of galaxies from our model uniformly selected in the IRX-β 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-β relation.

### 5.4.2. Impact of the attenuation curve

Recent observations done with the Atacama Large Millimeter Array (ALMA) by, for example, Capak et al. 2015; Bouwens et al. 2016 and Pope et al. 2017 ($z > 4$) or combined HST and Herschel/PACS measurements Reddy et al. (2018; $z > 1.5$) have mainly populated the IRX-β 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 Milky-Way attenuation curve) and the low metallicity expected for high-redshift galaxies.

Figure 14 shows the IRX-β 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. (21). 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 µ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-β 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 in Fig. 4, the slopes of these attenuation curves are different. We found a slope of -0.54 dex µm for Calzetti et al. (2000) (our reference), -0.64 dex µm for Reddy et al. (2015) and -1.23 dex µm for Gordon et al. (2003).

In Fig. 14, slopes of attenuation curves associated with our simulated galaxies and trends of average IRX-β 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 such as 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 that 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-β plan is fully linked to the capability of the G.A.S. model to associate a dedicated extinction curve and therefore a dedicated effective attenuation with each dust content/metallicity of each galaxy component (e.g., GMC and diffuse ISM).

### 5.4.3. Evolution with the redshift

Figure 15 shows the IRX-β relations for six different redshifts, from $z = 1.5$ to $z \approx 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 15% and 85% percentiles, we compute the scatter of the relation that slightly decreases with the redshift, i.e., from 0.5 dex at $z \approx 6.0$ to $\approx 1.0$ dex at $z = 1.5$. By taking into account these scatters, the slopes and behaviors of our IRX-β 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. Conclusions

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) a 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_\odot$) 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 us 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 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 \approx 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 0.8$. 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 IRLFs 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 IRLF, 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 \approx 7$.

The analysis of the IRLF 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 IRX. 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. Consistent 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 ($>100 M_\odot \text{ yr}^{-1}$) can reach IRX values that are 100 times larger than those predicted by the Meurer et al. (1999) fiducial relation.

With our model, we confirm that galaxies hosting an "old" ($>50$ Myr) 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 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 program managed by the ANR.

## References

Alavi, A., Siana, B., Richard, J., et al. 2014, ApJ, 780, 143

Andrews, B. H., & Martini, P. 2013, ApJ, 765, 140

Arnouts, S., Schiminovich, D., Ilbert, O., et al. 2005, ApJ, 619, L43

Atek, H., Richard, J., Jauzac, M., et al. 2015, ApJ, 814, 69

Baugh, C. M., Lacey, C. G., Frenk, C. S., et al. 2005, MNRAS, 356, 1191

Bell, E. F. 2002, ApJ, 577, 150

Béthermin, M., Le Floc'h, E., Ilbert, O., et al. 2012, A&A, 542, A58

Béthermin, M., Wu, H.-Y., Lagache, G., et al. 2017, A&A, 607, A89

Blumenthal, G. R., Faber, S. M., Flores, R., & Primack, J. R. 1986, ApJ, 301, 27

Boissier, S., de Gil Paz, A., Boselli, A., et al. 2007, ApJS, 173, 524

Boquien, M., Calzetti, D., Kennicutt, R., et al. 2009, ApJ, 706, 553

Boquien, M., Buat, V., Boselli, A., et al. 2012, A&A, 539, A145

Boquien, M., Boselli, A., Buat, V., et al. 2013, A&A, 554, A14

Bothwell, M. S., Kennicutt, R. C., Johnson, B. D., et al. 2011, MNRAS, 415, 1815

Bouwens, R. J., Illingworth, G. D., Franx, M., et al. 2009, ApJ, 705, 936

Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2010, ApJ, 708, L69

Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2012, ApJ, 754, 83

Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2015, ApJ, 803, 34

Bouwens, R. J., Aravena, M., Decarli, R., et al. 2016, ApJ, 833, 72

Bower, R. G., Benson, A. J., Malbon, R., et al. 2006, MNRAS, 370, 645

Bowler, R. A. A., Bourne, N., Dunlop, J. S., McLure, R. M., McLeod, D. J., et al. 2018, MNRAS, 481, 1631

Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000

Buat, V., Donas, J., Milliard, B., & Xu, C. 1999, A&A, 352, 371

Buat, V., Iglesias-Páramo, J., Seibert, M., et al. 2005, ApJ, 619, L51

Buat, V., Giovannoli, E., Burgarella, D., et al. 2010, MNRAS, 409, L1

Buat, V., Noll, S., Burgarella, D., et al. 2012, A&A, 545, A141

Burgarella, D., Buat, V., Gruppioni, C., et al. 2013, A&A, 554, A70

Calzetti, D. 2001, PASP, 113, 1449

Calzetti, D., Kinney, A. L., & Storchi-Bergmann, T. 1994, ApJ, 429, 582

Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682

Calzetti, D., Kennicutt, Jr., R. C., Bianchi, L., et al. 2005, ApJ, 633, 871

Capak, P. L., Carilli, C., Jones, G., et al. 2015, Nature, 522, 455

Caputi, K. I., Lagache, G., Yan, L., et al. 2007, ApJ, 660, 97

Casey, C. M., Berta, S., Béthermin, M., et al. 2012, ApJ, 761, 140

Casey, C. M., Narayanan, D., & Cooray, A. 2014, Phys. Rep., 541, 45

Chabrier, G. 2003, PASP, 115, 763

Chastenet, J., Bot, C., Gordon, K. D., et al. 2017, A&A, 601, A55

Chevallard, J., Charlot, S., Wandelt, B., & Wild, V. 2013, MNRAS, 432, 2061

Ciesla, L., Boquien, M., Boselli, A., et al. 2014, A&A, 565, A128

Clay, S. J., Thomas, P. A., Wilkins, S. M., & Henriques, B. M. B. 2015, MNRAS, 451, 2692

Compiègne, M., Verstraete, L., Jones, A., et al. 2011, A&A, 525, A103

Cortese, L., Boselli, A., Buat, V., et al. 2006, ApJ, 637, 242

Cortese, L., Boselli, A., Franzetti, P., et al. 2008, MNRAS, 386, 1157

Cousin, M. 2019, A&A, submitted (Paper III)

Cousin, M., Buat, V., Boissier, S., et al. 2016, A&A, 589, A109

Cousin, M., Guillard, P., & Lehnert, M. D. 2019, A&A, 627, A131

Cucciati, O., Tresse, L., Ilbert, O., et al. 2012, A&A, 539, A31

Cullen, F., McLure, R. J., Khochfar, S., Dunlop, J. S., & Dalla Vecchia, C. 2017, MNRAS, 470, 3006

da Cunha, E., Walter, F., Smail, I. R., et al. 2015, ApJ, 806, 110

De Lucia, G., & Blaizot, J. 2007, MNRAS, 375, 2

Devriendt, J. E. G., Guiderdoni, B., & Sadat, R. 1999, A&A, 350, 381

Draine, B. T., & Li, A. 2001, ApJ, 551, 807
Draine, B. T., & Li, A. 2007, ApJ, 657, 810

Draine, B. T., Dale, D. A., Bendo, G., et al. 2007, ApJ, 663, 866

Dwek, E., & Városi, F. 1996, AIP Conf. Ser., 348, 237

Ellis, R. S., McLure, R. J., Dunlop, J. S., et al. 2013, ApJ, 763, L7

Finkelstein, S. L., Papovich, C., Salmon, B., et al. 2012, ApJ, 756, 164

Finkelstein, S. L., Ryan, Jr., R. E., Papovich, C., et al. 2015, ApJ, 810, 71

Fontanot, F., Somerville, R. S., Silva, L., Monaco, P., & Skibba, R. 2009, MNRAS, 392, 553

Forrest, B., Tran, K.-V. H., Tomczak, A. R., et al. 2016, ApJ, 818, L26

Gordon, K. D., Clayton, G. C., Witt, A. N., & Misselt, K. A. 2000, ApJ, 533, 236

Gordon, K. D., Clayton, G. C., Misselt, K. A., Landolt, A. U., & Wolff, M. J. 2003, ApJ, 594, 279

Granato, G. L., Lacey, C. G., Silva, L., et al. 2000, ApJ, 542, 710

Gruppioni, C., Pozzi, F., Rodighiero, G., et al. 2013, MNRAS, 432, 23

Guiderdoni, B., & Rocca-Volmerange, B. 1987, A&A, 186, 1

Guo, Q., White, S., Boylan-Kolchin, M., et al. 2011, MNRAS, 413, 101

Hagen, L. M. Z., Hoversten, E. A., Gronwall, C., et al. 2015, ApJ, 808, 178

Hao, C.-N., Kennicutt, R. C., Johnson, B. D., et al. 2011, ApJ, 741, 124

Heinis, S., Buat, V., Béthermin, M., et al. 2013, MNRAS, 429, 1113

Heinis, S., Buat, V., Béthermin, M., et al. 2014, MNRAS, 437, 1268

Hernquist, L. 1990, ApJ, 356, 359

Hilton, M., Conselice, C. J., Roseboom, I. G., et al. 2012, MNRAS, 425, 540

Howell, J. H., Armus, L., Mazzarella, J. M., et al. 2010, ApJ, 715, 572

Inoue, A. K., Buat, V., Burgarella, D., et al. 2006, MNRAS, 370, 380

Jones, A. P., Tielens, A. G. G. M., Hollenbach, D. J., & McKee, C. F. 1994, ApJ, 433, 797

Jones, A. P., Tielens, A. G. G. M., & Hollenbach, D. J. 1996, ApJ, 469, 740

Karakas, A. I. 2010, MNRAS, 403, 1413

Kewley, L. J., & Ellison, S. L. 2008, ApJ, 681, 1183

Komatsu, E., Dunkley, J., Nolta, M. R., et al. 2009, ApJS, 180, 330

Kong, X., Charlot, S., Brinchmann, J., & Fall, S. M. 2004, MNRAS, 349, 769

Koprowski, M. P., Dunlop, J. S., Michałowski, M. J., et al. 2017, MNRAS, 471, 4155

Koprowski, M. P., Coppin, K. E. K., Geach, J. E., et al. 2018, MNRAS, 479, 4355

Lacey, C. G., Baugh, C. M., Frenk, C. S., et al. 2016, MNRAS, 462, 3854

Lilly, S. J., Carollo, C. M., Pipino, A., Renzini, A., & Peng, Y. 2013, ApJ, 772, 119

Livermore, R. C., Finkelstein, S. L., & Lotz, J. M. 2017, ApJ, 835, 113

Lo Faro, B., Buat, V., Roehlly, Y., et al. 2017, MNRAS, 472, 1372

Lucy, L. B., Danziger, I. J., Gouiffes, C., & Bouchet, P. 1989, in IAU Colloq. 120: Structure and Dynamics of the Interstellar Medium, eds. G. Tenorio-Tagle, M. Moles, & J. Melnick (Berlin: Springer Verlag), Lect. Notes Phys., 350, 164

Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415

Magnelli, B., Popesso, P., Berta, S., et al. 2013, A&A, 553, A132

Maniyar, A. S., Béthermin, M., & Lagache, G. 2018, A&A, 614, A39

Mannucci, F., Cresci, G., Maiolino, R., Marconi, A., & Gnerucci, A. 2010, MNRAS, 408, 2115

Marchetti, L., Vaccari, M., Franceschini, A., et al. 2016, MNRAS, 456, 1999

Martin, D. C., Fanson, J., Schiminovich, D., et al. 2005a, ApJ, 619, L1

Martin, D. C., Seibert, M., Buat, V., et al. 2005b, ApJ, 619, L59

Mathis, J. S., Mezger, P. G., & Panagia, N. 1983, A&A, 128, 212

McLure, R. J., Dunlop, J. S., Bowler, R. A. A., et al. 2013, MNRAS, 432, 2696

McLure, R. J., Dunlop, J. S., Cullen, F., et al. 2018, MNRAS, 476, 3991

Meurer, G. R., Heckman, T. M., & Calzetti, D. 1999, ApJ, 521, 64

Mo, H. J., Mao, S., & White, S. D. M. 1998, MNRAS, 295, 319

Narayanan, D., Davé, R., Johnson, B. D., et al. 2018, MNRAS, 474, 1718

Oesch, P. A., Bouwens, R. J., Carollo, C. M., et al. 2010, ApJ, 725, L150

Oesch, P. A., Bouwens, R. J., Illingworth, G. D., et al. 2013, ApJ, 773, 75

Overzier, R. A., Heckman, T. M., Wang, J., et al. 2011, ApJ, 726, L7

Pannella, M., Carilli, C. L., Daddi, E., et al. 2009, ApJ, 698, L116

Pannella, M., Elbaz, D., Daddi, E., et al. 2015, ApJ, 807, 141

Parsa, S., Dunlop, J. S., McLure, R. J., & Mortlock, A. 2016, MNRAS, 456, 3194

Peebles, P. J. E. 1969, ApJ, 155, 393

Pilyugin, L. S., & Thuan, T. X. 2005, ApJ, 631, 231

Pope, A., Montaña, A., Battisti, A., et al. 2017, ApJ, 838, 137

Popping, G., Puglisi, A., & Norman, C. A. 2017, MNRAS, 472, 2315

Reddy, N. A., Steidel, C. C., Pettini, M., et al. 2008, ApJS, 175, 48

Reddy, N. A., Erb, D. K., Pettini, M., Steidel, C. C., & Shapley, A. E. 2010, ApJ, 712, 1070

Reddy, N., Dickinson, M., Elbaz, D., et al. 2012, ApJ, 744, 154

Reddy, N. A., Kriek, M., Shapley, A. E., et al. 2015, ApJ, 806, 259

Reddy, N. A., Steidel, C. C., Pettini, M., & Bogosavljević, M. 2016, ApJ, 828, 107

Reddy, N. A., Oesch, P. A., Bouwens, R. J., et al. 2018, ApJ, 853, 56

Relaño, M., De Looze, I., Kennicutt, R. C., et al. 2018, A&A, 613, A43

Rémy-Ruyer, A., Madden, S. C., Galliano, F., et al. 2015, A&A, 582, A121

Rodighiero, G., Vaccari, M., Franceschini, A., et al. 2010, A&A, 515, A8

Salim, S., Rich, R. M., Charlot, S., et al. 2007, ApJS, 173, 267

Salim, S., Dickinson, M., Michael Rich, R., et al. 2009, ApJ, 700, 161

Salmon, B., Papovich, C., Long, J., et al. 2016, ApJ, 827, 20

Schreiber, C., Pannella, M., Elbaz, D., et al. 2015, A&A, 575, A74

Seibert, M., Martin, D. C., Heckman, T. M., et al. 2005, ApJ, 619, L55

Seon, K.-I., & Draine, B. T. 2016, ApJ, 833, 201

Takeuchi, T. T., Buat, V., & Burgarella, D. 2005, A&A, 440, L17

Takeuchi, T. T., Buat, V., Heinis, S., et al. 2010, A&A, 514, A4

Takeuchi, T. T., Yuan, F.-T., Ikeyama, A., Murata, K. L., & Inoue, A. K. 2012, ApJ, 755, 144

Tremonti, C. A., Heckman, T. M., Kauffmann, G., et al. 2004, ApJ, 613, 898

Treyer, M., Schiminovich, D., Johnson, B., et al. 2007, ApJS, 173, 256

Tweed, D., Devriendt, J., Blaizot, J., Colombi, S., & Slyz, A. 2009, A&A, 506, 647

Vieira, J. D., Marrone, D. P., Chapman, S. C., et al. 2013, Nature, 495, 344

Watson, D., Christensen, L., Knudsen, K. K., et al. 2015, Nature, 519, 327

Whitaker, K. E., Franx, M., Leja, J., et al. 2014, ApJ, 795, 104

Witt, A. N., Gordon, K. D., Vijh, U. P., et al. 2006, ApJ, 636, 303

Zahid, H. J., Kashino, D., Silverman, J. D., et al. 2014, ApJ, 792, 75

Zubko, V. G., Mennella, V., Colangeli, L., & Bussoletti, E. 1996, MNRAS, 282, 1321

Zubko, V., Dwek, E., & Arendt, R. G. 2004, ApJS, 152, 211
# 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:

$$ \frac{f_{\text{PAH}}}{f_{\text{PAH},0}} = 10^{-17.1+2.0(Z_g-0.2)} \quad \text{steeper envelope} \qquad \text{(A.1)} $$

$$ \frac{f_{\text{PAH}}}{f_{\text{PAH},0}} = 10^{-1.6+0.2(Z_g-0.2)} \quad \text{flatter envelope} \qquad \text{(A.2)} $$

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

In Fig. A.1 we compare the two modified relations with the original relation. According to Eqs. (A.1) and (A.2) the amount of PAH is reduced and increased in the gas phases, respectively. 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. A.2 the effective attenuation (Eq. (21)) produced by both the original ($A_{Orig}$) and modified ($A_{Modif}$) relations. The analysis is performed on a subsample of galaxies, randomly selected at different redshifts. Attenuations are computed for both the FUV and 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 normalization. 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. 2b). 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. A.3. 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

<figure><img src="image_19.png" /><figcaption>Fig. A.1. Fraction of PAH as a function of the gas metallicity of the dense/fragmented gas phase. The solid blue and green lines indicate the modified relations based on the flatter and steeper relations, respectively. The black solid line shows the original relation used in G.A.S. (black solid line in Fig. 1). The blue and green shaded areas show the full G.A.S. distribution generated by G.A.S. around the two modified relations. Data points are similar to those used in Fig. 1.</figcaption></figure>

decreases(increases). As for the effective attenuation in the FUV and NUV band, the higher is the redshift, the 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, the fraction of PAH evolves according to a random walk between the previous value, $f_{PAH}^{n-1}$ and a target value, $f_{PAH}^{n}$, expected from the empirical PAH-metallicity relation (see Eq. (3))

$$ f_{\text{PAH}}^{n} = f_{\text{PAH}}^{n-1} + r \times (f_{\text{PAH}}^{n} - f_{\text{PAH}}^{n-1}), \qquad \text{(A.3)} $$

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}^{n} = \frac{2}{3}(1 - f_{PAH}^{n})$. Then we apply

$$ f_{\text{BG}}^{n} = f_{\text{BG}}^{n-1} + r \times (f_{\text{BG}}^{n} - f_{\text{BG}}^{n-1}). \qquad \text{(A.4)} $$

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

$$ f_{\text{VSG}}^{n} = 1 - f_{\text{PAH}}^{n} - f_{\text{BG}}^{n}. \qquad \text{(A.5)} $$
<figure><img src="image_20.png" /></figure>

**Fig. A.2.** Impact of the modified $Z_g - f_{PAH}$ relations on the effective FUV and NUV attenuation (Eq. (21)). The difference between the original and modified relation is plotted as a function of the intrinsic FUV/NUV magnitude that does not evolve between the two configurations. The *upper* and *lower panels* show the impacts produced by the flatter and steeper relations, respectively. The solid and the dashed colored lines are associated with the FUV and NUV band, respectively. The effective attenuation, $A_{\text{Modif}}$ and $A_{\text{Orig}}$, are computed according to the modified and original $Z_g - f_{PAH}$ relations, respectively. Each colored line indicates a specific redshift from $z = 9.3$ to $z = 2.1$. The explicit list is given in the *left side* of the panels. The two other columns list the average dispersion measured for the FUV and NUV band around the mean relation.

<figure><img src="image_21.png" /></figure>

**Fig. A.3.** Impact of the modified $Z_g - f_{PAH}$ relations on the IRX. The *upper* and *lower panels* show the impacts produced by the flatter and steeper relations, respectively. The difference between the original and modified relation is plotted as a function of the intrinsic FUV magnitude that does not evolve between the two relations. Each colored line shows a given redshift. The list of redshift and the dispersion around the mean are indicated in the *bottom left corner* of each panel.

## Appendix B: GALAKSIENN library

The GALAKSIENN library stores the main results produced by our new G.A.S. SAM, especially MOCK galaxy catalogs and sky maps described in Paper III (Cousin 2019). It is available at the CDS and to the ZENODO platform<fnref n="5" />. A complete description of the GALAKSIENN library and the access through the CDS will be given in Paper III. In association with this Paper II, we distribute the ASCII tables of the FUV and IRLFs presented in Figs. 5 and 6, respectively.

<footnote n="5">⁵ https://zenodo.org, DOI: 10.5281/zenodo.1451229</footnote>