Fade to grey: systematic variation of the galaxy attenuation curves with galaxy properties in EAGLE

We present a simple model for galaxy attenuation by distilling SKIRT radiative transfer calculations for ~100,000 EAGLE galaxies at redshifts z=2-0. Our model adapts the two component screen model of Charlot&Fall (2000), parametrising the optical depth and slope of the ISM screen using the average dust surface density, $\Sigma_{\rm dust}$. We recover relatively tight relations between these parameters for the EAGLE sample, but also provide the scatter in these parameter owing to the morphological variation and orientation of galaxies. We also find that these relations are nearly independent of redshift in the EAGLE model. By pairing our model with an empirical prescription for birth clouds below the resolution scale of the simulation, we reproduce the observed relation between attenuation slope and optical depth for the first time in a cosmological simulation. We demonstrate that this result is remarkably independent of the attenuation properties assumed for birth cloud screen, merely requiring a boosted attenuation for infant stars. We present this model with a view to interpreting observations, as well as processing semi-analytic models and other hydrodynamic simulations.


INTRODUCTION
Despite constituting only a small portion of the total mass in galaxies, cosmic dust plays an outsized role in our understanding of galaxy properties. The interaction between light and dust complicates the interpretation of observations, and the derivation of fundamental properties such as galaxy stellar masses and star formation rates (e.g. Hopkins et al. 2001;Sullivan et al. 2001;Zibetti et al. 2010;Taylor et al. 2011;Leja et al. 2019).
Dust obscuration tends to reduce the emergent light from galaxies in UV to NIR wavelengths through scattering and absorption (e.g. Witt & Gordon 2000), the net effect of which is termed attenuation. A distinction is often made between attenuation and extinction, where extinction properties are intrinsic to the dust itself, i.e. purely derived from the absorption and scattering cross-sections of grains. Attenuation, on the other hand, depends on how this dust is E-mail: trayford@strw.leidenuniv.nl (JWT) situated relative to stars, and the aspect from which a galaxy is observed (e.g. Calzetti 2013).
With attenuation dependent on both the unique structure and dust composition of a galaxy, correcting for dust observationally is highly challenging (see e.g. Conroy 2013, for a review). The radiative transfer equations do not typically harbor analytic solutions for realistic configurations, and the existence of well documented degeneracies between dust effects and intrinsic properties of galaxies, such as stellar age and metallicity, further exacerbate the situation.
As a starting point, the case of a thin intervening dust screen between source and observer does have an analytic solution. In this particular configuration, the wavelength dependence of extinction and attenuation have the same form. Consequently, assuming this geometry and utilising point source pairs has allowed the extinction properties of the Milky Way dust to be inferred (e.g. Cardelli et al. 1989;Fitzpatrick 1999;Sasseen et al. 2002). In some cases this approach can be extended to derive the extinction properties of external galaxies, through occultation of one galaxy by another (e.g. Keel & van Soest 1992;Holwerda & Keel 2017).
However, for many scenarios the uniform screen geometry has drawbacks; it represents a maximally attenuating configuration for a fixed dust mass, and has peculiar features such as indefinitely increasing reddening with increased optical depth (Witt et al. 1992).
Still, at least one class of low-redshift galaxies is found to exhibit screen-like attenuation: UV-bright starbursts. This is evidenced by a tight correlation between their observed UV slope index (β) and excess Infrared radiation (IRX), demonstrating the characteristic increase in reddening with absorption (Calzetti et al. 1994;Meurer et al. 1999;Gordon et al. 1997). The IRX-β correlation has been found to hold for high redshift analogues (e.g. Reddy et al. 2010). Calzetti et al. (2000) derive a widely used screen model applicable for these starbursts. Shallower than the extinction measured for the Milky Way and local systems locally, the Calzetti et al. (2000) attenuation curve is theorised to represent a clumpy or turbulent shell geometry, reflecting the configuration around the nuclei of starbursting systems (Witt & Gordon 2000;Fischera et al. 2003). This is therefore an example of an effective screen; dust is applied as a screen, but the attenuation curve accounts for the influence of geometric effects.
A caveat for the Calzetti et al. (2000) attenuation model was that emission lines emanating from nebular regions tend to exhibit a factor ∼ 2 stronger attenuation than the stellar continuum. Charlot & Fall (2000) devise a two-component screen model to account for this, where stellar populations 10 Myr are obscured by stellar birth clouds, in addition to an overall diffuse ISM screen attenuating all stars. This model couples the wavelength dependent attenuation to the recent star formation history in galaxies, and was demonstrated to reproduce the SEDs of Sloane galaxies reasonably well (Bruzual & Charlot 2003). This model has been extended with different parametrisations for the two screens to represent different galaxy types or orientations (e.g. Wild et al. 2007;da Cunha et al. 2008). These two components provide some level of galaxy-galaxy variation in the attenuation curve even for screens of a fixed shape, but stronger systematic variations evidenced by observations (e.g. Salim et al. 2018) are harder to account for. To investigate how screen models could vary or break down with galaxy properties, a better understanding of dust propagation in realistic systems is needed to complement empirical studies (e.g. Wild et al. 2011;Kriek & Conroy 2013).
In order to move beyond the screen representation, radiative transfer (RT) codes have developed, allowing more complex geometries to be solved numerically (see e.g. Steinacker et al. 2013;Whitney 2011). These codes provide insight into the role of 'clumpiness' and different geometries in shaping empirical attenuation curves (e.g. Witt & Gordon 2000;Baes & Dejonghe 2001;Bianchi 2008), but the high dimensionality and computational expense of RT models make their direct application to observations (i.e. inverse modelling) difficult and their solutions potentially highly degenerate. By using highly constraining data and making some simplifying assumptions sophisticated RT models of local galaxies are now being produced (e.g. De Looze et al. 2014;Viaene et al. 2017;Verstocken et al. in prep), but these remain limited to an exclusive set of nearby systems.
Instead, insight may be gained from forward modelling galaxy formation simulations. These directly pre-scribe galaxy properties and structure, simplifying the radiative transfer modelling. This has been performed for semianalytic models of galaxy formation (e.g. SAMs, Fontanot et al. 2009;Gonzalez-Perez et al. 2013), and compared with observations statistically. However, SAMs can typically only provide the basic morphological components of galaxies (e.g. bulge, disc) that are again modelled using idealised geometries. Thus, they may miss the influence of multi-scale clumping and inhomogenity in the galaxy structure. Hydrodynamical simulations have the potential to produce structures that are more representative of real galaxies, modelling clumping and inhomogeneities self consistently above the resolution scale. Processing 'zoom' hydrodynamical simulations with RT has allowed representative attenuation properties to be produced and compared with observations, selfconsistently accounting for the complex geometric effects that arise (Jonsson et al. 2009;Wuyts et al. 2009;Saftly et al. 2015;Feldmann et al. 2017). In particular, Narayanan et al. (2018) demonstrate how the attenuation curve is shaped by geometric effects using the mufasa zoom simulations (Davé et al. 2016), presenting a model for how the relative obscuration of young and old stars change the attenuation curve slope and the strength of the 2175Å bump (Noll et al. 2009).
While these zooms may provide insight into important properties influencing the attenuation curve, they lack the statistical power and morphological diversity of cosmological samples. With modern computing power, hydrodynamical simulations of cosmological volumes can now realistically be processed using radiative transfer (Camps et al. 2016;Trayford et al. 2017;Rodriguez-Gomez et al. 2019). This allows us to investigate how the attenuation curve may vary statistically for a diverse cosmological sample. We present, to our knowledge, the first systematic analysis of RT attenuation curves generated for a cosmological sample of simulated galaxies, with the aim of parametrising a screen model for applications to observations, SAMs and other hydrodynamical simulations.
The paper is organised as follows. Section 2 first describes the EAGLE simulations and SKIRT radiative transfer modelling foundational for this work. Section 3 then outlines the approach we use to derive general dust attenuation properties from radiative transfer analysis of individual simulated galaxies (Camps et al. 2016;Trayford et al. 2017), particularly the separating the ISM and birth cloud contributions, and the computation of a model dust surface densities. This method is applied in section 4 to derive the strength and spectral shape of the ISM attenuation of EA-GLE galaxies as functions of dust surface density. Section 5 then discusses the potential incorporation of sub-resolution and birth-cloud attenuation, comparing to previous observational and theoretical studies on the shape of the attenuation curve. Finally, we discuss the conclusions of our study in section 6. Questions of numerical convergence, measuring dust surface densities and the quality of attenuation law fits are expanded upon in the appendices.

SIMULATIONS AND PROCESSING
Here, we briefly summarise two elements upon which this study relies; the eagle suite of cosmological simulations ( § 2.1, for full details see Schaye et al. 2015;Crain et al. 2015), and the calculation of dust attenuation using the 3D radiative transfer code, SKIRT (Baes et al. 2003(Baes et al. , 2011.

The EAGLE simulations
The eagle simulations are a suite of smoothed particle hydrodynamics (SPH) simulations, following galaxy and structure formation self-consistently within a series of cosmological volumes Crain et al. 2015). eagle is built upon a modified version of the gadget-3 code (an update of gadget-2, Springel 2005). In this study we focus on the largest simulation using the fiducial 'reference' physics model, Ref100, but make use of a number of smaller diagnostic simulations for the purposes of testing convergence properties. For brevity, we focus on the Ref100 simulation properties here, and defer details of the diagnostic simulations to appendix A.
The eagle suite assumes a ΛCDM cosmology, with cosmological parameters taken from the initial Planck data release (Planck Collaboration et al. 2014). The simulations are initialised from a Gaussian random field. The Ref100 run follows the evolution of a periodic, cubic volume of 100 3 cMpc on a side, at a particle mass resolution of m g 1.82×10 6 M (m DM = 9.7×10 6 M ) in gas (dark matter), with a Plummerequivalent gravitational softening length of 0.7 proper kpc at redshift z = 0.
The gadget-3 SPH is changed to use the pressureentropy implementation of Hopkins (2013), alongside a number of modifications to standard SPH collectively termed Anarchy (summarised in appendix A of Schaye et al. 2015). In order to account for key physical processes that are not resolved, a number of subgrid modules are implemented. These include schemes for star formation, mass loss and enrichment by stars, radiative cooling, photoheating, as well as energetic feedback associated with stellar evolution and supermassive black holes.
Radiative processes regulating the temperature of gas are included to complement the SPH hydrodynamical interaction. Radiative cooling and photoheating are included following Wiersma et al. (2009a), under the assumption of photoionisation equilibrium with the UV background of Haardt & Madau (2001). As cool, dense gas clouds (T < 10 4 K, n H > 0.1 cm −3 ) posesses a Jeans length, λ J , below the resolution scale of the simulations, they will fragment artificially. A pressure floor is thus imposed via a polytropic equation of state, artificially puffing up the ISM of galaxies such that λ J is always marginally resolved.
Non-zero star formation rates are assigned to gas surpassing a metallicity-dependent density threshold (Schaye 2004). Star particles are then formed from the wholesale conversion of gas particles, stochatically sampled at each timestep. The star particles inherit the metallicity and mass of their parent gas particle. The mass loss and enrichment of gas by these star particles is derived from stellar isochrones, assuming a Chabrier (2003) stellar IMF (Wiersma et al. 2009b), distributed over the local gas kernel. Eleven individual element abundances are tracked, alongside a total metallicity, Z. Gas particles store an individual abundance for each element, as well as this quantity smoothed over the local gas kernel. Smoothed abundances somewhat mitigate large particle-to-particle metallicity fluctuation in the absence of metal diffusion, and are used in this work unless stated otherwise.
Energetic feedback is another key aspect, regulating the gas content of galaxies. Star particles inject thermal energy into the ISM, alongside black hole particles which seed and grow in halos of mass M H > 10 10 h −1 M (Rosas- Guevara et al. 2015). Thermal feedback is implemented following Dalla Vecchia & Schaye (2012), with temperature increments of ∆T SF = 10 7.5 K and ∆T BH = 10 9 K (∆T BH = 10 9.5 K) for stellar and black hole feedback in the Ref (Rec) models, respectively. ∆T values are chosen to be high enough to mitigate numerical losses. These parameters, along with the behaviour of the energy fraction coupled to the gas, f th , are used to calibrate the simulations to reproduce the galaxy stellar mass function, size-mass relation and black hole to stellar mass relation.
A friends-of-friends algorithm is run on the fly in order to identify halos that form within eagle, with constituent substructures identified through post-processing the simulation outputs with the subfind algorithm (Springel 2005). For our purposes, eagle galaxies are taken to comprise the material within the central 30 pkpc of each subfind structure, with centres defined using a recursive shrinking spheres approach (see Trayford et al. 2017).

Radiative transfer & property maps
Virtual imaging, spectra and photometry are generated for galaxies that form within eagle via the SKIRT code 1 . SKIRT uses 3D Monte Carlo modelling to calculate dust effects for a given star-dust configuration and viewing angle. SKIRT can thus provide a dust calculation representative of the complex geometries and spatial correlations that naturally emerge within the simulated galaxies. The SKIRT modelling approach and survey strategy are fully delineated by Camps et al. (2016) and Trayford et al. (2017), with the resultant dataset described by Camps et al. (2018).
In summary, stellar populations older than 10 Myr emit light according to the GALAXEV population synthesis models Bruzual & Charlot (2003). Younger stars and their associated nebulae are represented using the MAPPINGS-III spectra for HII regions (Groves et al. 2008). The intrinsic extinction properties of dust grains are assumed to follow that of Zubko et al. (2004). Dust is assigned to sufficiently cold (T < 8000 K) gas by assuming a fixed fraction of metal mass resides in dust grains. The metal mass fraction, f dust = 0.3, as well as the PDR covering fraction used as input to MAPPINGS-III, f PDR = 0.1, were set by Camps et al. (2016) to best reproduce FIR properties of the nearby galaxies from the Herschel Reference Survey (HRS, Boselli et al. 2010) within observational bounds (e.g. Inoue & Kamaya 2004;Mattsson et al. 2014).
SKIRT uses a Gaussian smoothing kernel, with smoothing lengths taken from the simulation, to construct a 3D dust grid through which stochastically emitted monochromatic photon packets can propagate. These interact with the dust grid through scattering and absorption, and are compiled into spectra for a fixed wavelength grid, with and without dust effects. Not counting dust re-emission in the IR, 333 wavelength points are set between 280 nm and 2500 nm, providing sufficient resolution to recover photometry across the UV-NIR range.
For this work, we make use of photometry with and without dust radiative transfer for galaxies with log 10 (M /M ) > 10 (log 10 (M /M ) > 9) at reference (high) resolution. In addition to the standard photometry, we also obtain photometry for light emanating from MAPPINGS-III HII regions alone, as well as coeval GALAXEV spectra 2 . We also utilise maps of physical properties for each galaxy, described in Trayford & Schaye (2018), which were produced using the pySPHViewer code (Benítez-Llambay et al. 2018). These allow us to derive dust maps following our dust prescription, and thus compute dust surface densities.

MODELLING ATTENUATION WITH EAGLE
The attenuation of a galaxy's intrisic SED emerges from the highly complex radiative transfer of light through an attenuating medium. Attenuation depends simultaneously on the stellar-dust geometry, the redirection of light by scattering and the detailed properties of the dust grains. Dust attenuation in galaxies is also multi-scale; it depends on both macroscopic structure of galaxies on kpc scales, and the cold pcscale clouds in which stars are born. As such, it is important to consider the limitations of simulations such as EAGLE in reproducing these complexities, alongside the insights the simulation provides.
With this in mind, we describe the fundamental aspects of our model, which predicts the ISM attenuation curves of galaxies using their dust surface density. First, the choice of screen paradigm is discussed in § 3.1. The measurement of dust surface densities is then explored in § 3.2. Finally, we describe the fitting of attenuation curves to individual galaxies in § 3.3. Calzetti et al. (2000, C00) derive an effective attenuating screen for starburst galaxies, which Fischera et al. (2003) interpret to represent a turbulent dusty medium of MWlike dust composition, with a lognormal distribution of column densities and a well-mixed stellar population. While this model may be a good approximation for some galaxies, it is unlikely a good approximation in general, particularly if certain stellar populations are preferentially located in dustier regions.

Screen models
The attenuation law of C00 has been adapted to account for such variation. Following e.g. Noll et al. (2009), a simple parametrisation is where k SB is the wavelength dependent form of the attenuation curve for starburst galaxies given by C00, with δ representing a power law tilt relative to the standard relation.
The value of δ must then be chosen. Such a form has been adopted by numerous studies (e.g. Chevallard et al. 2013;Salmon et al. 2016;Narayanan et al. 2018).
Charlot & Fall (2000, CF00) introduced a twocomponent screen model to account for the lines-of-sight towards younger stellar populations having higher column densities of dust. In this model, generalised by Wild et al. (2007), older stars are attenuated by a dust screen of fixed optical depth representing the interstellar medium (ISM), while lines of sight towards young stars have a boosted optical depth due to an additional term, taken to represent the birth clouds (BCs). This gives an attenuation curve, τ(λ), for a stellar population with age t age , of Where τ ISM represents the ISM optical depth, τ BC is the optical depth of stellar birth clouds, t disp is the birth cloud dispersal time and η ISM and η BC are the spectral slopes for the two attenuation curves. In their fiducial parametrisation, CF00 take η ISM = η BC = −0.7, t disp = 10 7 yr and τ BC = 2τ ISM .
This model was adopted for EAGLE galaxies by Trayford et al. (2015), and subsequently compared to radiative transfer calculations using the SKIRT code by Trayford et al. (2017). In the comparison with the SKIRT results, the τ ISM term is taken to represent the dust tracing the resolved gas distribution in EAGLE galaxies (on scales 1 kpc), and the τ BC term is taken to be unresolved, and represented by the built in dust attenuation in the MAPPINGS-III spectral libraries for HII regions (Groves et al. 2008) used to represent young stars. We find that this analogy between the components generally works well, and fitting the SKIRT results with the model of Trayford et al. (2015) recovers close to the fiducial model parameters.

Dust surface density
The surface density of dust, Σ dust , is clearly a key parameter in the conception of a physically motivated model for dust attenuation. As dust is not tracked explicitly in the EAGLE simulation, we assume a simple scaling, with where the dust mass, m dust , is deemed to constitute a fixed fraction, f dust , of the metal mass for a gas particle of mass m gas and metallicity Z. In order to measure dust properties, 256x256 pixel mass maps are produced for each galaxy within a 60 2 kpc 2 field of view using the approach described in Trayford & Schaye (2018), using f dust = 0.3, consistent with the SKIRT data.
In the case of an idealised uniform screen with fixed dust composition, Σ dust is the sole parameter dictating the absorption of light from a source. In a realistic configuration, however, the attenuation will depend on the relative distribution of dust with respect to stars in galaxies. To account for this, any measure of Σ dust should be related to the stellar distribution. To test how the measurement method of Σ dust may influence our results, we try two different definitions. The first uses a circular aperture to define the area over which the surface density is measured, Where M dust (< r) is the mass in dust within projected radius r. The projected stellar half-mass radius, R e , can then be used to relate this to the stellar distribution. We consider both Σ dust (R e ) and Σ dust (2R e ). A second method is to calculate the dust surface density weighted by the projected stellar mass. To do this, we combine the the stellar and dust mass maps, Where for pixel i, M dust,i and M ,i represent the projected dust and stellar mass, and l pix is the side length of each pixel. For our standard maps, l pix = 60/256 kpc ≈ 0.235 kpc.
In summary, for all of the results in this work, we have tested for the influence of dust surface density definition using three definitions of dust surface density; Σ dust (R e ), Σ dust (2R e ) and Σ dust . We find that there is no one definition that clearly correlates best with attenuation properties (see appendix B for details), so we take Σ dust = Σ dust (R e ) by default.

Measuring ISM attenuation
In order to measure the attenuation by the diffuse ISM modelled with SKIRT, we use two sets of photometry for each galaxy. The first is the full SKIRT photometry presented in Trayford et al. (2017), the second follows the same procedure with the ISM removed. For a given photometric band b, the AB rest-frame magnitudes with and without ISM attenuation are notated as b and b respectively. The ISM attenuation in that band is then The unique star-dust geometry and composition of each galaxy yields a distinct ISM attenuation for each dustbearing EAGLE galaxy in each band. We fit a power law, representing the ISM term of Equation 2, τ ISM (λ) = τ ISM (λ/5500Å) −η ISM , to the ISM effective optical depths of each of the ugriz SDSS bands for individual galaxies, treating τ ISM and η ISM as fitting parameters. The band representative wavelength is taken to be the effective wavelength of each filter.
We use a simple least-squares fitting routine to fit each curve, assuming representative SDSS photometric errors to weight each band. These are taken to be 0.01 mag in the griz bands and 0.02 mag in the u-band (Padmanabhan et al. 2008). We note that a power law generally fits the bands well and, as such, using unweighted bands produces very similar results. The quality of these power law fits are demonstrated in more detail in appendix C. Compiling the data, we turn to the relationships between the individual band ISM optical depths, best fit power attenuation law parameters and dust surface densities in Section 4.

EAGLE ISM ATTENUATION CURVES
We now present the attenuation properties of the interstellar medium in EAGLE galaxies calculated using SKIRT, and The non-parametric A V and R V values of galaxies calculated for the Ref100 EAGLE simulation run, as a function of their average dust surface density, Σ dust . We show individual galaxies as points at discrete redshifts of z = 0.1, 0.5, 1 and 2, coloured according to the legend. Coloured bars similarly indicate the 16-84 percentile range of A ISM r in uniform log 10 Σ dust /(M kpc −2 ) bins at each redshift. A coloured contour enclosing 90% of galaxies for is also shown for each redshift. The scatter on the overall relation (combining the distinct redshift samples) is plotted in black, with a cubic spline fit to the A ISM r medians shown as a dashed black line. We see that the EAGLE galaxies follow clear, redshift-independent relationships in these parameters.
following the methodology described in section 3. For reasons of economy 3 , we consider galaxies at 4 specific redshift outputs of the simulation; z =0.1, 0.5, 1 and 2, covering ≈80% of the age of the universe.
In Fig. 1 the non-parametric attenuation properties A V and R V are shown as functions of the dust surface density, R V represents the local gradient of the attenuation curve.
We first consider the median A V as a function of Σ dust for all redshifts, represented by the black dashed line. This relationship appears close to linear for Σ dust 10 5.5 M kpc −2 , and asymptotically approaches A V = 0 for Σ dust 10 5.5 M kpc −2 . The relationship also appears quite tight, as indicated by the 16 − 84th percentile ranges, suggesting the Σ dust alone is a good predictor of the V-band attenuation in EAGLE galaxies. The contributions of EAGLE galaxies at fixed redshifts are represented by the appropriately coloured points, contours and errorbars, as indicated in the legend. Perhaps the most striking feature of this is a remarkable concordance between the relation at different redshifts. While the contours indicate that higher Σ dust EAGLE galaxies are more prevalent at high redshift, the A V distribution at fixed Σ dust is highly consistent. The median values for each redshift in the same bins as the overall relation are shown as coloured bars (staggered for visibility), allowing a more direct comparison to be made. The redshift independence implies that the evolution in the typical A V of galaxies in EAGLE reflects galaxies sampling different ranges of a near-static A V -Σ dust relation.
A further detail is the scatter in values for a fixed Σ dust , as indicated by the 16 − 84th percentile ranges. We use half of this range in magnitudes to represent the dispersion, σ. For considerable levels of attenuation (A V 0.25), Σ dust captures most of the variance in A V , becoming more dominant toward higher Σ dust as the scatter becomes a smaller fraction of the median attenuation. Still the scatter at a fixed Σ dust is considerable, representing the influence of more complex aspects of the star-dust geometry and orientation. At log 10 Σ dust /(M kpc −2 ) ≈ 6 we find a typical attenuation of A V ≈ 1 and dispersion is σ ≈ 0.2 mag. This dispersion increases with Σ dust , reaching σ ≈ 0.5 mag by log 10 Σ dust /(M kpc −2 ) ≈ 7. For galaxies at fixed Σ dust , the difference in the dispersion of A V with redshift is marginal.
We now turn to the trend between Σ dust and R V , shown in the bottom panel of Fig. 1. Only galaxies with A ISM V > 0.25 mag are included in this plot, as R V values become noisy at very low attenuations owing to shot noise in the SKIRT photometry. Once again, we focus first on the redshiftaveraged trend, shown in black. At log 10 Σ dust /(M kpc −2 ) ≈ 6, where A V ≈ 1, the slope is intermediate between the C00 curve and the dedicated ISM attenuation curve used by CF00. We find that the ISM attenuation curve tends to become greyer (or shallower, higher R V ) at higher Σ dust . Across the range log 10 Σ dust /(M kpc −2 ) = 5 − 7 the typical ISM attenuation curve shallows from essentially a C00 curve, to significantly greyer than the fiducial CF00 attenuation.
As with the A V -Σ dust trend, the galaxy samples at fixed redshift appear to sample different Σ dust values along a static median R V relation. Focusing on individual galaxies, we see a significant tail to high R V , with points outside the plotted range indicated using up-arrow markers. It is worth noting that the definition of R V implies that as the slope of the attenuation curve tends to flat, R V tends towards arbitrarily high values. For a handful of galaxies R V is negative, implying that dust actually makes the galaxy appear bluer in B−V than it is intrinsically. We note that these galaxies are all at highly negative R V values, implying slightly positive (but essentially flat) attenuation curve slopes. For the majority of these galaxies, this can be attributed to shot noise in the photometry. Still, as mentioned in Trayford et al. (2017), it is possible that the star-dust geometry conspires to make the EAGLE galaxy bluer in some rare scenarios 4 .
The scatter in the R V values is again indicated by the 16-84th percentile range. While the total scatter in R V increases with Σ dust , this is primarily due to the scatter to higher R V , while the 16-50th percentile range is relatively constant. This is indicative of the skew to high R V becoming stronger with Σ dust . Considering the individual redshift subsamples, we find a marginal decrease in the scatter with increasing redshift.
Altogether, the shape and normalisation of the ISM attenuation curves of EAGLE galaxies show strong systematic variation with Σ dust , largely independent of redshift 5 . This is a remarkable result; EAGLE displays strong evolution in gas metallicity (Guo et al. 2016), gas content (Lagos et al. 2015), galaxy size (Furlong et al. 2017) and star formation Furlong et al. 2015), all of which may influence the attenuation properties of galaxies.
Given these non-parametric results, we now consider the result of fitting power law attenuation curves (equivalent to the CF00 ISM term, described in § 3.2). In Fig. 2, the best-fit τ ISM 550 and η ISM values are again compiled for redshift 0.1, 0.5, 1 and 2 snapshots in the top and bottom panels of Fig. 2, respectively. Black errorbars represent the 16-84th percentile range for all selected galaxies, in contiguous bins of log Σ dust .
The black dashed line then shows a cubic spline fit to the median values in each bin. Focussing first on the top panel, the τ ISM 550 relation shows a very similar trend to that of A V in Fig. 1, both in terms of shape and scatter. This is unsurprising, but demonstrates that the V-band attenuation is reproduced well by fitting a power law to the ugriz bands. The quality of the fits are explored further in Appendix C. The line-of-sight measurements of Kreckel et al. (2013) are included for comparison. Point markers are indicative of the nearby galaxy from which they are measured. Here, green points represent continuum measurements and grey points are measured from emission lines by Kreckel et al. (2013). Given that emission line measurements are biased towards HII regions, and thus associated with significant stellar birth cloud absorption, we deem the green points more appropriate as a comparison to the pure ISM absorption plotted for EAGLE. The Σ dust values are typically measured on sub-kpc scales, and are obtained by fitting IR emission maps pixel-by-pixel with a model for the dust mass (see Aniano et al. 2012 for details). The relation for two idealised cases are also plotted; a uniform screen (blue) and a homogeneous slab (orange). These curves are derived using the C00 attenuation law, following Kreckel et al. (2013).
The EAGLE attenuation demonstrates a relation intermediate between the two idealised cases. A primary factor causing the τ ISM 550 values to be lower than predicted for a screen can be attributed to the mixture of stars and dust along the line of sight. As a result, stars are typically situated behind lower dust surface densities than obtained by pro-  (2000) Example Dust Map Example Star Map Figure 2. Demonstrating the attenuation law derived using EAGLE, and its dependence on the average dust surface density, Σ dust . The top panel shows the effective ISM optical depth as a function of Σ dust . Points indicate the median optical depth of Ref100 EAGLE galaxies, with error bars indicating the 16-84 percentile range. Dashed line is a cubic spline fit to the median points. Thick blue and orange lines indicate idealised screen and slab geometries respectively, as in Kreckel et al. (2013). The data from Kreckel et al. (2013) is also plotted, using different symbols to indicate independent lines of sight from nearby galaxies. Grey and green symbols indicate the attenuation measured from emission line diagnostics and stellar continuum fitting, respectively. Bottom panel shows the power law slope of the attenuation curve in the optical range (0.2-0.8 µm). Black dashed lines and points represent the Ref100 run as above. The thick orange line again represents an idealised slab, with the fiducial CF00 value shown as a dot-dashed line. Example stellar and dust surface density maps of a contributing galaxy are inset. We see that the extreme, idealised geometries bracket the τ ISM 550 -Σ dust relation, while EAGLE attenuation curve slope varies significantly with Σ dust . jecting over the entire galaxies. Indeed, a foreground screen represents a maximally absorbing uniform configuration for the stars and dust Witt et al. (1992). While the slab model accounts for the mixing of dust and stars, the inhomogeneity of the star-dust geometry distinguishes it from this idealised case. In particular, radial gradients in the star and gas distribution, and the physical association of young stars with denser ISM (not just birth clouds) contribute to differences in the trends.
Generally the EAGLE attenuation is closer to the slab case, and appears to exhibit a similar plateau at high Σ dust values. In the slab, this is a skin effect where the influence of adding more dust to the configuration diminishes as sources on the near edge of the attenuating medium become dominant. This is likely similar in the EAGLE case, where τ ISM 550 saturates as the less obscured stellar populations come to dominate the observed light. This plateau is only marginally resolved in the EAGLE case, visible only at the highest Σ dust values. The scatter at a fixed Σ dust again also indicates geometric variations in the galaxies.
Comparing to the Kreckel et al. (2013) observations, we see these are also generally intermediate between the two idealised cases. The continuum measurements tend to occupy similar values as found for EAGLE over the 5.5 ≤ log 10 Σ dust /(M kpc −2 ) < 6.5 range, though exhibit strong scatter. At lower Σ dust , the observational limitations and uncertainties may contribute to the flat trend observed for e.g. NGC3077 (Kreckel et al. 2013). Uncertainties in the observed Σ dust values depend on both systematics in the modelling and observational uncertainties. While typical uncertainties for individual sight lines are quoted at 0.1 dex, this varies between galaxies and could boost scatter in the x-axis significantly for some systems.
At higher Σ dust , the observations appear to exhibit higher values than we measure. This could point to a more screen-like attenuation for theses lines-of-sight, or more preferential attenuation of bright young stars than found in EA-GLE. This seems consistent with the finding of more homogeneous and 'puffed-up' ISM in EAGLE galaxies relative to observations, due to numerical effects (Trayford et al. 2017;Benítez-Llambay et al. 2018), and may influence the EA-GLE attenuation measurement such that they are closer to a homogeneous slab case. However, we note that these data points are not directly comparable to the EAGLE trend as they represent lines of sight and not galaxy integrated quantities, and that they may also be influenced by small scale birth cloud attenuation, particularly in dusty starbursts such as NGC2146 (Jackson & Ho 1988).
The bottom panel of Fig. 2 demonstrates how the shape of the attenuation curve depends on Σ dust , via the best-fitting power law index for the ISM attenuation, η ISM . Again, the dashed line indicates the cubic spline fit to the bin medians, while error bars show the 16th-84th percentile range. The EAGLE trend becomes significantly greyer (i.e. less steep) with increasing Σ dust . As was also demonstrated for R V in the bottom panel of Fig 1, the power law slope of the attenuation curve plateaus to a value close that of the intrinsic extinction measured for the Milky Way. This is perhaps unsurprising, given that the average Milky Way dust extinction curve of Zubko et al. (2004) is the input extinction curve for dust in galaxies. With increased Σ dust , and thus optical depth, geometric effects produce greyer curves, reaching the fiducial power law slope of CF00 by log 10 Σ dust /(M kpc −2 ) ≈ 6.2. At higher Σ dust , the slope shallows further but also shows evidence of a plateau.
This can be compared again to a slab model. The slab follows a similar qualitative trend, tending towards to shape of intrinsic curve (C00) below log 10 Σ dust /(M kpc −2 ) ≈ 5.5, crossing the CF00 threshold at log 10 Σ dust /(M kpc −2 ) ≈ 6.8. At log 10 Σ dust /(M kpc −2 ) 6 the slab case is a good approximation to the slope of the EAGLE attenation curve. The higher attenuation in EAGLE relative to the slab is ascribed the preferential reddening of young stars due to being embedded in preferentially dustier ISM regions, even aside from the effects of stellar birth cloud attenuation not accounted for here. The brighter stars in young stellar populations combined with this preferential association has the consequence of yielding stronger attenuation overall relative to the slab case. We note that the slope for the screen case is constant by definition.
With both the V-band ISM attenuation and the power law slope of the ISM attenuation demonstrating strong trends with the underlying dust surface density, the relationship between these two attenuation properties is shown directly in Fig. 3. As before, the dashed line shows a spline fit to the median η ISM values in contiguous 0.33 mag bins bin for bins with > 50 galaxies, as indicated in the colour bar. We see that attenuation curves become greyer (flatter) as attenuation increases, becoming shallower than the standard η ISM = −0.7 assumption (CF00) at A ISM V ≥ 0.75. For a given attenuation, more edge-on (higher inclination) galaxies also tend to have greyer attenuation curves.
of A ISM V , with error bars indicating the 16th-84th percentile scatter. The significant flattening of the attenuation curve with attenuation can be seen, with the median slope value approaching that of CF00 at A ISM V ≈ 0.8 mag. The underlying colour map indicates the residual trend between galaxy orientation and η ISM as functions of A ISM V , with shading indicating the median offset of galaxies from the median inclination angle at that A ISM V . At a fixed A ISM V , higher inclination (more edge-on) galaxies exhibit greyer attenuation curves. Taking the difference in this way removes the primary trend that higher A ISM V measurements tend to be associated with more edge-on galaxies. This effect has been previously noted for real galaxies (e.g. Wild et al. 2011), and attributed to the effect that the attenuation in face-on galaxies is governed by the physical association between ISM and young stars, while in edge-on systems dust lanes are less discriminating between differently aged populations.
Together, these results indicate how trends in ISM attenuation properties (computed via radiative transfer) are established through the internal structure and orientation of galaxies, using a cosmological sample of virtual galaxies taken from the EAGLE simulations. This complex emergent structure within simulated galaxies modifies the attenuation, even when using a fixed input extinction curve (Zubko et al. 2004). Such trends support previous theoretical work (e.g. Fontanot et al. 2009;Chevallard et al. 2013;Narayanan et al. 2018) and motivate a more sophisticated prescription for dust attenuation in galaxy models and semi-analytic models (SAMs), with the dust surface density largely dictating the wavelength dependence and strength of attenuation in galaxies. This could also be incorporated into SED fitting tools, given the important implications such trends between spectral slope and attenuation strength are when inferring properties of the underlying stellar populations.
As previously stated, a caveat for this modelling is that internal galaxy structures are limited by resolution and numerical effects in the simulations, particularly with galaxy discs not being thin enough (e.g. Trayford et al. 2017;Benítez-Llambay et al. 2018), and eagle not resolving compact structures. The following section 5 explores how to incorporate the birth cloud term to enable a better comparison with observations.
Despite these limitations, we present this ISM screen model as a more nuanced alternative to idealised geometric models typically used when interpreting observations, incorporating the diverse galaxy morphologies that emerge within the eagle simulations.

THE BIRTH CLOUD TERM
As the SKIRT data can only provide constraints on the influence of star-dust geometry on scales above the EAGLE resolution limit ( 700 pc), the influence of nebular structures, represented by the τ BC term of equation 2, remains unconstrained. Below we explore the treatment of the τ BC term, and incorporate τ BC for comparison to data.

The contribution of infant stellar populations
To gauge the possible impact of differing birth cloud treatments on the overall attenuation, we must first assess how much the affected stars contribute to the total emitted light at different wavelengths. We assume stellar populations with ages t age < 10 Myr are affected by this BC term, following the original CF00 implementation, and we term these 'infant stars'.
In Fig. 4 we show the fraction of the total flux emanating from infant stars as a function of wavelength in z = 0.1 EAGLE galaxies, binned by specific star formation rate. In order to obtain these curves, we first stack the star formation histories of z = 0.1 galaxies taken from the Ref100 simulation in each sSFR bin. The composite histories for all stars and the infant portion alone (t age < 10 Myr) are run separately through SKIRT (without dust transfer) to generate spectra across the UV-NIR wavelength range. The flux fraction in infant stars is then the ratio of the infant to total stellar spectra. The wavelength ranges of a number of bands are plotted beneath these curves, indicated by the labelled shaded regions. We note that this fraction is in the absence of dust, so when birth clouds or a boosted ISM attenuation are included for the infant stars, the emergent flux fraction will be lower.
We see that, unsurprisingly, the flux contribution by infant stars generally increases towards shorter wavelengths and is systematically higher in higher sSFR bins. The infant star contribution remains low across most of the SDSS ugriz photometric bands, reaching ≈ 10 % ( 30 %) for the highest sSFR bin in the g-band (u-band). This suggests that while non-negligible, the infant stars contribute a small fraction across the ugriz range used to fit attenuation curves, and NUV u g r i z log10(SSFR/yr 1 ) 10 10.5 log10(SSFR/yr 1 ) < 10 11 log10(SSFR/yr 1 ) < 10.5 12 log10(SSFR/yr 1 ) < 11 log10(SSFR/yr 1 ) < 12 Figure 4. The fraction of the emitted flux emanating from infant stellar populations (with t age < 10 Myr) as a function of wavelength, computed using stacked EAGLE star formation histories. Binned star formation and enrichment histories of z = 0.1 Ref100 EAGLE galaxies with log 10 (M /M ) > 9.5 are stacked for bins of SSFR (see text), and Bruzual & Charlot (2003) spectra are used to compute the fraction of intrinsic flux emitted by infant stars. The 10th to 90th wavelength percentile range of the SDSS and GALEX filter transmission curves are indicated for reference. We see that in SDSS bands, infant stars remain subdominant, with galaxies in the highest SSFR bin only reaching a 10% contribution on the blue side of the g-band.
as such will remain a subdominant contributor to the effective attenuation across this range, even if the light from infant stars is completely obscured. While the star formation histories are compiled using z = 0.1 galaxies, The bins are chosen to include more extreme galaxies representative of the sSFRs at higher redshifts.
While this reveals that the ISM term generally dominates the effective attenuation in ugriz, the choice of birth cloud treatment may still be important in some regimes. At shorter wavelengths, we see that the infant star fraction rises markedly, up to ≈ 50% in the NUV. The infant star contribution will also be high for certain atomic transitions that emanate from compact HII regions across the UV-NIR range. In what follows, we explore how including an additional birth cloud term may influence the effective attenuation in galaxies.

The attenuation slope relation
The modified C00 law of equation 1 provides a general fit to the attenuation in galaxies, accounting for the influence of ISM and small-scale structure simultaneously. The power law modifier, δ, and normalisation, A V , parametrise the shape and strength of attenuation, and the relationship between the two has been predicted in theoretical studies (Witt et al. 1992;Chevallard et al. 2013;Narayanan et al. 2018) and inferred observationally (e.g. Salmon et al. 2016;Salim et al. 2018;Decleir et al. 2019). A shallowing trend of greyer (lower-δ) attenuation with increasing A V is observed in general, with quantitative differences between studies.
We first consider the ISM-only attenuation curves (black points, Fig. 5). The EAGLE ISM-only relation shows systematically greyer attenuation curves, with a wavelength dependence λ 0.4 times weaker. This offset appears remarkably constant, and the shape of the A V -δ relation appears to reproduce the observations very well.
Steeper attenuation curves are expected from models that include the influence of small scale and birth cloud absorption; the young stellar populations that are associated with birth clouds are relatively blue, so their preferential attenuation leads to redder stellar spectra. The MAPPINGS-III spectral libraries include an implicit dust attenuation associated with the starlight emanating from HII regions, combined with nebular absorption and emission. Given the inextricability of this dust correction, we instead replace the resampled MAPPINGS-III SEDs with GALAXEV spectra. This enables full control over any 'subgrid' attenuation that is applied.
EAGLE cannot provide insight into the nature of true birth cloud attenuation, as these structures are well below the scale that can be resolved. Instead, some assumptions must be made about the nature of birth cloud attenuation.
The high optical depths of stellar birth clouds makes it difficult to determine the wavelength dependence of birth cloud attenuation empirically, as the enshrouded stars contribute little to the emergent flux (e.g. da Cunha et al. 2008). Theoretically, a shell-like geometry may be a better description of the birth clouds that form around infant populations, driven by the strong radiation pressure from the short-lived OB stars within. In such a configuration, a dust screen is actually a reasonable representation of the birth cloud attenuation, and the extinction and attenuation curves converge. Fitting a power law to the extinction properties of an average Milky Way dust composition yields a power law index of η BC = −1.3 (Wild et al. 2007;da Cunha et al. 2008, see equation 2), significantly steeper than the η ISM values at high attenuation (Fig. 2) and the fiducial CF00 values of η BC = η ISM = −0.7. We denote this index value as η shell BC = −1.3.
The other primary uncertainty is in the normalisation of the optical depth assumed for the birth clouds. Typically, birth cloud optical depths are assumed to be higher than that of the ISM. A fixed ratio between τ BC and τ ISM is often assumed, with the default CF00 ratio of τ BC = 2τ ISM . Note that some studies fitting attenuation curves treat this ratio as a free parameter, but in lieu of a well motivated range for this parameter, we opt for a fixed ratio in our fiducial model.
Finally, in order to actually implement the birth cloud attenuation for eagle galaxies in post-processing, we use both the total SEDs of EAGLE galaxies, and the SED contributions of infant stars with ages t age < 10 Myr. These infant star SEDs were prepared using both the MAPPINGS-III prescription used in the fiducial SKIRT model for EAGLE (Camps et al. 2016;Trayford et al. 2017) and the GALAXEV (Bruzual & Charlot 2003) stellar population models, excluding nebular emission and in-built dust effects.
To obtain the new flux in a given band, f new , we sub-tract away the MAPPINGS-III contributions and add the GALAXEV SEDs for infant stars, corrected for the influence of stellar birth clouds, as where f total is the fiducial SKIRT flux, f HII is for MAPPINGS-III only and f BC03 is for the same stellar populations represented by the pure GALAXEV spectra. All fluxes used include the influence of ISM attenuation from SKIRT. f τ is the fixed ratio assumed for τ BC /τ ISM . The red points in Fig. 4 show the δ-A V relationship for the birth-cloud corrected EAGLE fluxes, assuming η BC = η shell ISM and f τ = 2. We see that δ shifts to systematically lower values at a fixed A V , indicating steeper (redder) attenuation curves. The shift decreases gradually with A V , largest in the low A V bin. This can be understood as a result of the higher contribution of infant stars to the total attenuation at low effective A V , due to the dominance of the τ BC term. At high A V infant stars are effectively hidden, and, given the limited fraction that they can possibly contribute (Fig. 4), the birth cloud attenuation no longer makes a significant difference to the emergent fluxes.
Comparing to the Salim et al. (2018) data, we see that including this birth cloud prescription improves the agreement with the observations. The BC-corrected fluxes agree within the scatter across the range, rather than only at the lowest attenuations (A V < 0.5). While the agreement of the ISM only relation is similar to that of Narayanan et al. (2018, their figure 11, left panel).
Despite this relative success, the BC-corrected fluxes do produce a noticeably steeper relationship than is recovered observationally, transitioning from underpredicting to overpredicting the observed δ. To see how this may be influenced by our choice of the BC prescription parameters, η BC and f τ , we also trial a number of different parametrisations, shown as thin red lines. We try the CF00 value of η BC = −0.7, as well as boosting f τ = 5 or uniformly sampling f τ over the range [2,10] for each galaxy. One might imagine a shallower relationship could be obtained using a shallower BC attenuation slope (η BC closer to 0), in combination with stronger attenuation (higher f τ ) to match the normalisation. However, this parametrisation (dashed), along with all the aforementioned variations, yield remarkably similar overall relations to the fiducial treatment. Clearly as f τ tends to zero the trend should approach the ISM only points. The dash-dotted line shows the result of using a lower f τ = 1 with η BC = −0.7, appearing marginally closer to the η BC = −1.3, f τ = 2 points. Together this suggests that, provided some reasonable attenuation boost is applied to infant stars ( f τ 2), the δ − A V relation is largely insensitive to the exact birth cloud prescription.
It is difficult, then, to ascertain exactly what factors lead to the slightly steeper relation we measure. Shortcomings in the realism of eagle galaxies could play a role here. In particular, we know that the pressure floor imposed in eagle leads to some level of artificial 'puffing-up' of the ISM, yielding galaxies that are somewhat more homogeneous. This could generally lead to more 'slab-like' attenuation properties (as seen in Fig. 2), where the stars and dust are more uniformly mixed than they is realistic. This is closely re-  Salim et al. (2018) for total attenuation is included. We see that the ISM-only attenuation is systematically greyer (lower-δ) than observed, but that the addition of birth cloud attenuation improves agreement.
lated to resolution effects, as there is missing structure on scales below ∼ 1 kpc and the birth clouds we account for are typically 10 − 100 pc in scale (see appendix A for further discussion). Also, because this is an effective screen model, the attenuation properties depend on the star formation histories (e.g. the burstiness of galaxies) which may not be representative.
Assumptions in the dust modelling could contribute too, e.g. a fixed dust composition and dust-to-metal ratio is assumed throughout the SKIRT modelling, that may vary in concert with these relations in real galaxies. We may also be hampered by the limitations of the two-component screen model representation we use; the single screen for infant stars does not account for the gradual dispersal and partial covering of realistic birth clouds. The true properties of birth clouds are fundamentally mysterious, and may well evolve with redshift alongside the size mass and turbulence of giant molecular clouds. In particular, we do not explore variation in the dispersal time for birth clouds in different environments, which remains a caveat for any implementation of such a screen model.
Finally, we caution that inferring the δ and A V from real galaxy spectra is highly challenging and likely highly degenerate. This may introduce systematic effects that distort the shape of the observed relation from its true value. It should be possible the same inference procedure applied to observations could be applied to eagle spectra, yielding a fairer comparison, but this is left for future work.

SUMMARY & CONCLUSIONS
In this study we present a model for effective dust attenuation parametrised solely by the projected dust surface density, Σ dust , of galaxies along a given line of sight. The model develops the two component screen model of Charlot & Fall (2000) by distilling information gained from full radiative transfer of ∼ 100, 000 eagle galaxies in random orientation (via the SKIRT code Trayford et al. 2017;Camps et al. 2016). This is motivated by first identifying a strong dependence of the non-parametric A V and R V values on Σ dust . By then fitting functional forms to the ISM attenuation curves for individual galaxies, we compare the relations between Σ dust and the slope and strength of ISM attenuation to those provided by both idealised geometries and inferred from observation.
Some consideration is made for the treatment of birth clouds when computing overall attenuation, given that such structures are far below the resolution of the simulation. We model the previously studied relationship between the strength and slope of the attenuation curve, with and without birth cloud attenuation. Our model can then be compared with prior works from observational (Wild et al. 2011;Salmon et al. 2016;Salim et al. 2018) and theoretical (e.g. Fontanot et al. 2009;Chevallard et al. 2013;Narayanan et al. 2018) perspectives.
Our primary findings are as follows: • The relationship between Σ dust and A V is both remarkably tight and highly independent of redshift. The Σ dust -R V dependence is less tight, but remains strong and demonstrates a similar level of redshift independence (Fig. 1).
• The best fit power law attenuation curves show optical depths, τ ISM 550 , comparable to line-of-sight observations in local galaxies, and bracketed by the theoretical extreme cases of a screen and a perfectly-mixed slab. The attenuation curve slope, η ISM , is close to that of the slab case at high Σ dust ( 10 6 M kpc −2 ) and close to a screen (i.e. extinction curve) case at low Σ dust ( 10 5 M kpc −2 ).
• We demonstrate that the scatter in the η ISM versus A V relation (Fig. 3) has a strong residual trend with inclination, such that greyer attenuation curves are preferentially found for edge-on systems in accordance with observations (e.g. Wild et al. 2011).
• The A V -δ 6 relation for ISM-only attenuation in eagle is shown to reproduce the shape of the Salim et al. (2018) data well, with an offset to greyer attenaution curves similar to that shown in Narayanan et al. 2018 (Fig. 5). We find that by applying a birth cloud attenuation term, attenuation curves become systematically redder showing improved agreement with the observations.
• We find that this result is largely insensitive to the assumed attenuation properties of infant stars, for reasonable values of the birth cloud optical depth and attenuation slope (assuming a fixed birth cloud dspersal time of 10 Myr). Indeed, a number of different parametrisations yield near identical results. The birth cloud parametrisation is therefore not considered as major source of uncertainty in our model.
Given the limited scope of this study, several aspects of the SKIRT-EAGLE attenuation curves are not investigated here. For example, by using full spectra the behaviour of the 2000Å bump could also be investigated, and the insights into how bump strength correlates with attenuation slope made by Narayanan et al. (2018) could be explored further in comparison to data (e.g. Kriek & Conroy 2013;Tress et al. 2019). It would also be desirable to extend this model to include how the absorbed light is re-emitted in the IR, using the self consistent modelling of SKIRT. This would allow the model to predict the FIR SEDs of galaxies, for example. These aspects are left for future work.
The model presented in this work allows users to sample the typical attenuation curves for the ISM in galaxies where the dust surface density can be measured, along with representative scatter in the slope and strength of attenuation. It is hoped that this will find use when forward modelling galaxies where the physical gas mass, metallicity and size are known or assumed. In particular, this can provide better motivated attenuation corrections for semi-analytic models of galaxy formation (SAMs), where galaxy structure is unknown or idealised. While radiative transfer have previously been assumed for SAMS via idealised geometries (e.g. Fontanot et al. 2009), the complex emergent geometries of galaxies represented by the eagle simulation can be incorporated. Lagos et al. (2019) present a first application of this result, generating panchromatic luminosity functions for the shark SAM (Lagos et al. 2018), including stellar and dust emission.
It is important to consider the convergence of the dust attenuation properties derived in this work. As discussed in section 3, our modelling associates the ISM term of the two component CF00 model to the attenuation by superkpc scale structures calculated for fiducial EAGLE galaxies using the SKIRT code. The birth cloud term may then be attributed to unresolved structures on the scale of HII regions. To test if smaller scale structures present in higher resolution simulations have a significant effect on the attenuation properties, we can appeal to higher resolution diagnostic EAGLE simulations, listed in Table A1.
The Ref25 simulation uses the same resolution and physics model as Ref100, but realises a 25 3 cMpc volume. The Rec25 and RefHi25 simulations are initialised Table A1. Key parameters of the primary simulations used in this work. From left to right: simulation label, side length of cubic volume L in co-moving Mpc (cMpc), initial mass of gas particles m g , Plummer equivalent gravitational softening prop at redshift z = 0 in proper kpc (pkpc), and the study reference for each volume. for the same volume and initial conditions, but with enhanced particle mass resolution of m g 2.26 × 10 6 M (m DM = 1.21 × 10 6 M ). Rec25 differs rom RefHi25 in that the physics model is modified from fiducial, as feedback efficiencies are recalibrated to better match the galaxy stellar mass function at z = 0.1 (see below). Fig. A1 shows convergence properties for parameters of the best-fit power law ISM attenuation curve for each simulated galaxy, previously seen in Fig. 2. The left panel shows that the relatively tight relationship between Vband (550 nm) optical depth and dust surface density found for Ref100 (black) is preserved extremely well between the higher resolution RecHi25 (blue) and RefHi25 (green) simulation runs. Additionally, the concordance of the Ref25 trend suggests that the limited environmental variation sampled by these 25 cMpc volumes has little influence over the τ ISM 550log 10 Σ dust /(M kpc −2 ) relation, and thus does not compensate for resolution difference. The V-band attenuation is thus deemed well converged across the 4 log 10 Σ dust /(M kpc −2 ) 6.5 range sampled by the 25 cMpc simulations.
The power law slope η ISM , however shows more variation, indicative of a more pronounced difference between the behaviour of attenuation curve shapes between simulations. At higher Σ dust (log 10 Σ dust /(M kpc −2 ) 5), the Ref25 and RefHi25 simulations follow the Ref100 trend closely, matching the average values and scatter in η ISM . However the RecHi25 simulation produces greyer (shallower) attenuation curves, by approximately a factor of λ 0.2 relative to Ref100. The slightly greyer attenuation curves could be attributed to a clumpier medium in the RecHi25 simulations, an effect theorised to affect attenuation curves in general (e.g. Fischera et al. 2003). The fact that this effect is observed in the RefHi25 and not the RecHi25 simulation suggest that this is not a purely numerical effect of sampling smaller scales and higher densities in the ISM, but rather requires the slightly enhanced stellar feedback employed by the RecHi25 model. Small deviations are also seen at low values of Σ dust . However, at log 10 Σ dust /(M kpc −2 ) 5 the attenuation values are small enough that such variations would likely not be measurable.
We note that the higher resolution EAGLE simulations are able to indicate how the smaller spatial scales sampled by higher resolution runs may effect the integrated attenuation, but are limited by insufficient small scale physics, lacking molecular gas modelling and harboring an artificially  Figure A1. Convergence properties of the best-fit V -band ISM optical depth, τ ISM 550 , and power law slope, η ISM , as functions of Σ dust , originally shown in Fig. 2. Left and right-hand panels plot the Σ dust dependence of τ ISM 550 and η ISM , respectively. The primary Ref100 simulation used in this work is compared to the Ref25, RefHi25 and Recal25 volumes (detailed in Table A1) to the influence of numerical convergence and related effects. The cubic spline fit to binned medians and the 16th-84th percentile scatter are plotted for each simulations following Fig. 2, and colour coded as indicated in the legend. The agreement suggests good convergence in the strength of ISM attenuation, and moderate convergence in the shape of the attenuation curves. pressurised ISM 7 . Pairing improved resolution with molecular cooling in future simulations may lead to thinner discs, and would likely have some effect on the attenuation properties of galaxies. However, we pose our attenuation model as an advancement over idealised geometric models, which begin to show the trends induced by preferential attenuation and complex emergent geometries in a physically motivated model.
Altogether, the attenuation properties are moderately well converged showing a highly consistent normalisation and slope for the ugriz range. The treatment of attenuation associated with sub-resolution birth clouds is investigated in section 5.

APPENDIX B: MEASURING DUST SURFACE DENSITY
A key part of the modelling performed in this work is the choice of how the physical dust surface density, Σ dust , is measured. On the one hand the choice of Σ dust should be representative of the typical dust surface density seen by stars, while on the other it should be a simple enough quantity to compute that it is easily applicable to physical models for galaxies over a broad range of complexity.
Perhaps the most complicated Σ dust value we compute is 7 See Trayford et al. (2017) for further dicussion of the effects of the puffing-up and homogenisation of the EAGLE ISM on the SKIRT results one where we take the M weighted average Σ dust over each pixel of our property maps described in section 2, Σ dust . The M weighting is intended to yield Σ dust values representing the typical Σ dust towards stars in the galaxies, e.g. an unobscured bulge would significantly down-weight Σ dust . A simpler method is to use the property maps to measure the Σ dust within some circular radius, which we take to be multiples of the stellar half mass radius, Σ dust (r < R e ) and Σ dust (r < 2R e ). While this is no longer dependent on the relative alignment of dust and stars within the aperture, the choice of stellar half mass radius means aperture size depends on the stellar distributions, such that more compact galaxies will tend to higher Σ dust . Fig. B1 shows the median relations of Fig. 2 for these three choices of Σ dust measurement. We find that, perhaps unsurprisingly, the Σ dust value is typically highest on average, as density profiles for stars and gas tend to increase towards the centers of galaxies. The Σ dust (r < R e ) values are systematically lower on average, shifting the relations to ≈ 0.25 dex lower Σ dust . The larger aperture Σ dust (r < R e ) yields values that are lower still, shifting the relations by a further ≈ 0.35 dex in Σ dust .
The scatter in the relations is remarkably similar for each Σ dust measure. While we observe that individual galaxies can scatter by very large values in Σ dust between metrics, the average scatter is largely preserved. As a result, there is no obvious candidate for the Σ dust measurement most representative of the attenuation. As we deem the aperture measurements the simplest, we choose Σ dust (r < R e ) as our fiducial measure.  Figure B1. Different assumptions for Σ dust , demonstrated by the relations between Σ dust and best fitting power law attenuation profile, as in Fig. 2. Σ dust is measured either as the stellar mass surface density weighted average value, Σ dust , and the average Σ dust within once and twice the stellar half-mass radius, Σ dust (r < R e ) and Σ dust (r < 2R e ), respectively. We see that Σ dust produces the highest surface densities on average, followed by Σ dust (r < R e ) and Σ dust (r < 2R e ) producing lower surface densities, with an overall ∼0.6 dex shift between the median relations.  Figure C1. The main panel shows the histogram of goodness-of-fit, via a reduced chi squared (χ 2 red ) measure, for the power law fit to the ISM attenuation of EAGLE galaxies with A V > 0.25. We note that χ 2 red is measured assuming representative errors for SDSS photometry on each band, as opposed to the true photometric error of SKIRT (due to shot noise in the photon sampling), which is typically much smaller. Three regions of the χ 2 red distribution are highlighted in green, amber and red, representing three levels of decreasing fit quality. 2 attenuation curves are randomly sampled from each region and are plot in the inset panel (solid lines), along with the best fitting CF00 power-law (Eq. 2, fine dotted lines) and modified C00 law (Eq. 1, coarse dotted lines). We see that the functional forms generally give qualitatively good fits, with discrepancies typically arising from a sharper uptick in the measured attenuation from g to u band.