TY - JOUR
AU - Obreschkow, Danail
AB - ABSTRACT We present a simple model for galaxy attenuation by fitting skirt radiative transfer calculations for ∼10 000 eagle galaxies at redshifts z = 2 − 0. Our model adapts the two-component screen model of Charlot & Fall, parametrizing the optical depth and slope of the interstellar medium screen using the average dust surface density, Σdust. We recover relatively tight relations between these parameters for the eagle sample, but also provide the scatter in these parameters 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 the 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. dust, extinction, galaxies: formation, galaxies: ISM, galaxies: structure 1 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, Charlot & Rix 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 grain species, and the distribution of different grain sizes. Attenuation, on the other hand, depends on how this dust is 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 Walcher et al. 2011 or Conroy 2013 for recent reviews). The radiative transfer (RT) equations do not typically harbour 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 exacerbates 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 utilizing point source pairs have allowed the extinction properties of the Milky Way dust to be inferred (e.g. Cardelli, Clayton & Mathis 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, Thronson & Capuano 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 generally found between the observed UV slope index (β) and excess Infrared radiation (IRX), demonstrating the characteristic increase in reddening with absorption (Calzetti, Kinney & Storchi-Bergmann 1994; Gordon, Calzetti & Witt 1997; Meurer, Heckman & Calzetti 1999). IRX-β correlations have been found to hold for some high-redshift analogues (e.g. Reddy et al. 2010, though see Casey et al. 2014 and Capak et al. 2015 for examples otherwise). Calzetti et al. (2000, hereafter C00) derive a widely used screen model applicable for these starbursts. Shallower than the extinction measured for the Milky Way and local systems, the C00 attenuation curve is theorized to represent a clumpy or turbulent shell geometry, reflecting the configuration around the nuclei of starbursting systems (Witt & Gordon 2000; Fischera, Dopita & Sutherland 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 C00 attenuation model was that emission lines emanating from nebular regions tend to exhibit a factor approximately two stronger attenuation than the stellar continuum. Charlot & Fall (2000, hereafter CF00) devise a two-component screen model to account for this, where stellar populations ≲10 Myr are obscured by stellar birth clouds (BCs), in addition to an overall diffuse interstellar medium (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 Sloan galaxies reasonably well (Bruzual & Charlot 2003). This model has been extended with different parametrizations for the two screens to represent different galaxy types or orientations (e.g. Wild et al. 2007; da Cunha, Charlot & Elbaz 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, Boquien & Lee 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, RT codes have developed, allowing more complex geometries to be solved numerically (see e.g. Whitney 2011; Steinacker, Baes & Gordon 2013). 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 preparation), but these remain limited to an exclusive set of nearby systems. Instead, insight may be gained from forward modelling galaxy formation simulations. These directly prescribe galaxy properties and structure, simplifying the RT modelling. This has been performed for semi-analytic models (SAMs) of galaxy formation (e.g. 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 idealized geometries. Thus, they may miss the influence of multiscale clumping and inhomogeneity 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 hydrodynamical simulations of a single galaxy halo with RT has allowed representative attenuation properties to be produced and compared with observations, self-consistently accounting for the complex geometric effects that arise (Jonsson, Groves & Cox 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é, Thompson & Hopkins 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 RT (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 parametrizing a screen model for applications to observations, SAMs, and other hydrodynamical simulations. The paper is organized as follows. Section 2 first describes the eagle simulations and skirt RT modelling foundational for this work. Section 3 then outlines the approach we use to derive general dust attenuation properties from RT analysis of individual simulated galaxies (Camps et al. 2016; Trayford et al. 2017), particularly separating the ISM and BC 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 eagle galaxies as functions of dust surface density. Section 5 then discusses the potential incorporation of sub-resolution and BC 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. 2 SIMULATIONS AND PROCESSING Here, we briefly summarize two elements upon which this study relies: the eagle suite of cosmological simulations (Section 2.1, for full details see Crain et al. 2015; Schaye et al. 2015), and the calculation of dust attenuation using the 3D RT code, skirt (Baes et al. 2003, 2011; Camps & Baes 2015). 2.1 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; Schaye 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 use 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 XX 2014). The simulations are initialized from a Gaussian random field. The Ref100 run follows the evolution of a periodic, cubic volume of 100 cMpc on a side, at a particle mass resolution of mg ⪆ 1.82 × 106 M⊙ (mDM = 9.7 × 106 M⊙) in gas (dark matter), with a Plummer-equivalent gravitational softening length of 0.7 proper kpc at redshift z = 0. The gadget-3 SPH is changed to use the pressure–entropy implementation of Hopkins (2013), alongside a number of modifications to standard SPH collectively termed anarchy (summarized 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, Schaye & Smith (2009a), under the assumption of photoionization equilibrium with the UV background of Haardt & Madau (2001). As cool, dense gas clouds (T < 104 K, nH > 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, stochastically sampled at each time step. 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 initial mass function (Wiersma et al. 2009b), distributed over the local gas kernel. 11 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 haloes of mass MH > 1010h−1 M⊙ (Rosas-Guevara et al. 2015). Thermal feedback is implemented following Dalla Vecchia & Schaye (2012), with temperature increments of ΔTSF = 107.5 K and ΔTBH = 109 K (ΔTBH = 109.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, fth, 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 haloes 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). 2.2 Radiative transfer and property maps Virtual imaging, spectra, and photometry are generated for galaxies that form within eagle via the skirt code (Camps & Baes 2015).1skirt 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 data set 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 H ii regions (Groves et al. 2008). The intrinsic extinction properties of dust grains are assumed to follow that of Zubko, Dwek & Arendt (2004), which are demonstrated to reproduce the Milky Way extinction curve of Fitzpatrick (1999) well. 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, fdust = 0.3, as well as the photodissociation region (PDR) covering fraction used as input to MAPPINGS-III, fPDR = 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 and 2500 nm, providing sufficient resolution to recover photometry across the UV–NIR range. For this work, we use photometry with and without dust RT for galaxies with log10(M⋆/M⊙) > 10 (log10(M⋆/M⊙) > 9) at reference (high) resolution. In addition to the standard photometry, we also obtain photometry for light emanating from MAPPINGS-III H ii regions alone, as well as coeval galaxev spectra.2 We also utilize 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. 3 MODELLING ATTENUATION WITH eagle The attenuation of a galaxy’s intrisic SED emerges from the highly complex RT 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 multiscale; it depends on both macroscopic structure of galaxies on kpc scales, and the cold pc-scale 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 Section 3.1. The measurement of dust surface densities is then explored in Section 3.2. Finally, we describe the fitting of attenuation curves to individual galaxies in Section 3.3. 3.1 Screen models C00 derive an effective attenuating screen for starburst galaxies, which Fischera et al. (2003) interpret to represent a turbulent dusty medium of MW-like 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. The attenuation law of C00 has been adapted to account for such variation. Following e.g. Noll et al. (2009), a simple parametrization is $$\begin{eqnarray*} {\tau} {(\lambda)} = {0.4}\ln {(10)} \, {E(B-V)}\, k_{\rm SB}{(\lambda)} \, \left ( \frac{\lambda}{5500\,\mathring{\rm A}} \right )^{\delta }, \end{eqnarray*}$$(1) where kSB 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). CF00 introduced a two-component screen model to account for the lines of sight towards younger stellar populations having higher column densities of dust. In this model, generalized by Wild et al. (2007), older stars are attenuated by a dust screen of fixed optical depth representing the ISM, while lines of sight towards young stars have a boosted optical depth due to an additional term, taken to represent the BCs. This gives an attenuation curve, τ(λ), for a stellar population with age tage, of $$\begin{eqnarray*} \tau (\lambda) = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\tau _{\rm ISM}\left(\frac{\lambda }{5500 {\mathring{\rm A} }}\right)^{\eta _{\rm ISM}} & {\rm for} \,\, t_{\rm age} \gt t_{\rm disp}\\ \tau _{\rm ISM}\left(\frac{\lambda }{5500 {\mathring{\rm A} }}\right)^{\eta _{\rm ISM}} + \tau _{\rm BC}\left(\frac{\lambda }{5500 {\mathring{\rm A} }}\right)^{\eta _{\rm BC}} & {\rm for} \,\, t_{\rm age} \le t_{\rm disp} \end{array}\right., \nonumber\\ \end{eqnarray*}$$(2) where τISM represents the ISM optical depth, τBC is the optical depth of stellar BCs, tdisp is the BC dispersal time, and ηISM and ηBC are the spectral slopes for the two attenuation curves. In their fiducial parametrization, CF00 take ηISM = ηBC = −0.7, tdisp = 107 yr, and τBC = 2τISM. This model was adopted for eagle galaxies by Trayford et al. (2015), and subsequently compared to RT 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 additional τBC term is taken to correspond to the unresolved attenuation, represented by the built-in dust attenuation in the MAPPINGS-III spectral libraries for H ii 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. With this separation, it is also more straightforward to avoid double-counting attenuation when we introduce additional BC attenuation in Section 5. 3.2 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 $$\begin{eqnarray*} m_{\rm dust} = f_{\rm dust} Z m_{\rm gas}, \end{eqnarray*}$$(3) where the dust mass, mdust, is deemed to constitute a fixed fraction, fdust, of the metal mass for a gas particle of mass mgas and metallicity Z. In order to measure dust properties, 256 × 256 pixel mass maps are produced for each galaxy within a 602 kpc2 field of view using the approach described in Trayford & Schaye (2018), using fdust = 0.3, consistent with the skirt data. In the case of an idealized 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, $$\begin{eqnarray*} \Sigma _{\rm dust}(\lt r) = \frac{M_{\rm dust}(\lt r)}{2\pi r^{2}}, \end{eqnarray*}$$(4) where Mdust(<r) is the mass in dust within projected radius r. The projected stellar half-mass radius, Re, can then be used to relate this to the stellar distribution. We consider both Σdust(Re) and Σdust(2Re). A second method is to calculate the dust surface density weighted by the projected stellar mass. To do this, we combine the stellar and dust mass maps, $$\begin{eqnarray*} \langle \Sigma _{\rm dust}\rangle = \frac{1}{l_{\rm pix}^2}\frac{\sum _{i=0}^{n_{\rm pix}}{M_{{\rm dust}, i}M_{\star , i}}}{\sum _{i=0}^{n_{\rm pix}}{M_{\star , i}}}, \end{eqnarray*}$$(5) where for pixel i, Mdust,i and M⋆,i represent the projected dust and stellar mass, and lpix is the side length of each pixel. For our standard maps, lpix = 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(Re), Σdust(2Re), 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(Re) by default. 3.3 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 Ab,ISM = b − b′, with the effective optical depth defined as τISM,b = −0.4ln (10)Ab,ISM. The unique star-dust geometry and composition of each galaxy yields a distinct ISM attenuation for each dust-bearing eagle galaxy in each band. We fit a power law, representing the ISM term of equation (2), |$\tau _{\rm ISM}(\lambda) = \tau _{{\rm ISM}}(\lambda /5500{\mathring{\rm A} })^{-\eta _{\rm 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-fitting power attenuation law parameters, and dust surface densities in Section 4. 4 eagle ISM ATTENUATION CURVES We now present the attenuation properties of the ISM in eagle galaxies calculated using skirt, and following the methodology described in Section 3. For reasons of economy,3 we consider galaxies at four specific redshift outputs of the simulation; z = 0.1, 0.5, 1, and 2, covering ≈80 per cent of the age of the Universe. In Fig. 1, the non-parametric attenuation properties AV and RV are shown as functions of the dust surface density, Σdust. AV represents the normalization of the ISM attenuation in the V band with AV = V − V′, while RV measures the total-to-selective reddening, i.e. RV = AV(B − V)/(B′ − V′). RV represents the local gradient of the attenuation curve. Figure 1. Open in new tabDownload slide The non-parametric AV and RV 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. Thecoloured bars similarly indicate the 16th–84th percentile range of |$A^{\rm ISM}_r$| in uniform log10Σdust/(M⊙ kpc−2) bins at each redshift. A coloured contour enclosing 90 per cent of galaxies 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^{\rm ISM}_r$| medians shown as a dashed black line. We see that the eagle galaxies follow clear, redshift-independent relationships in these parameters. Figure 1. Open in new tabDownload slide The non-parametric AV and RV 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. Thecoloured bars similarly indicate the 16th–84th percentile range of |$A^{\rm ISM}_r$| in uniform log10Σdust/(M⊙ kpc−2) bins at each redshift. A coloured contour enclosing 90 per cent of galaxies 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^{\rm ISM}_r$| medians shown as a dashed black line. We see that the eagle galaxies follow clear, redshift-independent relationships in these parameters. We first consider the median AV as a function of Σdust for all redshifts, represented by the black dashed line. This relationship appears close to linear for Σdust ≳ 105.5 M⊙ kpc−2, and asymptotically approaches AV = 0 for Σdust ≲ 105.5 M⊙ kpc−2. The relationship also appears quite tight, as indicated by the 16th−84th percentile ranges, suggesting that Σ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 error bars, 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 Σdusteagle galaxies are more prevalent at high redshift, the AV 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 AV of galaxies in eagle reflects galaxies sampling different ranges of a near-static AV–Σdust relation. A further detail is the scatter in values for a fixed Σdust, as indicated by the 16th−84th percentile ranges. We use half of this range in magnitudes to represent the dispersion, σ. For considerable levels of attenuation (AV ≳ 0.25), Σdust captures most of the variance in AV, becoming more dominant towards 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 log10Σdust/(M⊙ kpc−2) ≈ 6, we find a typical attenuation of AV ≈ 1 and dispersion is σ ≈ 0.2 mag. This dispersion increases with Σdust, reaching σ ≈ 0.5 mag by log10Σdust/(M⊙ kpc−2) ≈ 7. For galaxies at fixed Σdust, the difference in the dispersion of AV with redshift is marginal. We now turn to the trend between Σdust and RV, shown in the bottom panel of Fig. 1. Only galaxies with |$A^{\rm ISM}_{V} \gt 0.25$| mag are included in this plot, as RV values become noisy at very low attenuations owing to shot noise in the skirt photometry. Once again, we focus first on the redshift-averaged trend, shown in black. At log10Σdust/(M⊙ kpc−2) ≈ 6, where AV ≈ 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 RV) at higher Σdust. Across the range log10Σ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 AV–Σdust trend, the galaxy samples at fixed redshift appear to sample different Σdust values along a static median RV relation. Focusing on individual galaxies, we see a significant tail to high RV, with points outside the plotted range indicated using up-arrow markers. It is worth noting that the definition of RV implies that as the slope of the attenuation curve tends to flat, RV tends towards arbitrarily high values. For a handful of galaxies RV 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 RV 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 RV values is again indicated by the 16th−84th percentile range. While the total scatter in RV increases with Σdust, this is primarily due to the scatter to higher RV, while the 16th–50th percentile range is relatively constant. This is indicative of the skew to high RV becoming stronger with Σdust. Considering the individual redshift subsamples, we find a marginal decrease in the scatter with increasing redshift. Altogether, the shape and normalization 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; Schaye 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 Section 3.2). In Fig. 2, the best-fitting |$\tau ^{\rm 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. The black error bars represent the 16th−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. Figure 2. Open in new tabDownload slide 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. The points indicate the median optical depth of Ref100 eagle galaxies, with error bars indicating the 16th–84th percentile range. The dashed line is a cubic spline fit to the median points. The thick blue and orange lines indicate idealized screen and homogeneously mixed slab geometries, respectively, as in Kreckel et al. (2013). The data from Kreckel et al. (2013) are also plotted, using different symbols to indicate independent lines of sight from nearby galaxies. The grey and green symbols indicate the attenuation measured from emission line diagnostics and stellar continuum fitting, respectively. The bottom panel shows the power-law slope of the attenuation curve in the optical range (0.2–0.8 μm). The black dashed lines and points represent the Ref100 run as above. The thick orange line again represents the homogeneously mixed case, 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, idealized geometries bracket the |$\tau _{550}^{\rm ISM}$|–〈Σdust〉 relation, while the eagle attenuation curve slope varies significantly with Σdust. Figure 2. Open in new tabDownload slide 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. The points indicate the median optical depth of Ref100 eagle galaxies, with error bars indicating the 16th–84th percentile range. The dashed line is a cubic spline fit to the median points. The thick blue and orange lines indicate idealized screen and homogeneously mixed slab geometries, respectively, as in Kreckel et al. (2013). The data from Kreckel et al. (2013) are also plotted, using different symbols to indicate independent lines of sight from nearby galaxies. The grey and green symbols indicate the attenuation measured from emission line diagnostics and stellar continuum fitting, respectively. The bottom panel shows the power-law slope of the attenuation curve in the optical range (0.2–0.8 μm). The black dashed lines and points represent the Ref100 run as above. The thick orange line again represents the homogeneously mixed case, 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, idealized geometries bracket the |$\tau _{550}^{\rm ISM}$|–〈Σdust〉 relation, while the eagle attenuation curve slope varies significantly with Σdust. Focusing first on the top panel, the |$\tau ^{\rm ISM}_{550}$| relation shows a very similar trend to that of AV 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, the green points represent continuum measurements and the grey points are measured from emission lines by Kreckel et al. (2013). Given that emission-line measurements are biased towards H ii regions, and thus associated with significant stellar BC 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 idealized cases are also plotted; a uniform screen (blue) and a slab where stars and dust are homogeneously mixed (orange). These cases are modelled as described by Kreckel et al. 2013 (their Section 3.3), with the foreground screen employing the Milky Way dust model of Draine et al. (2007), and the mixed slab model following the ‘mixed’ example of Calzetti et al. (1994) including isotropic scattering (C00). We emphasize that these idealized cases are independent of the star formation histories, dust content, and geometries of eagle galaxies. The eagle attenuation demonstrates a relation intermediate between the two idealized cases. A primary factor causing the |$\tau ^{\rm 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 projecting 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 mixed slab model accounts for the mixing of dust and stars, the inhomogeneity of the star-dust geometry distinguishes it from this idealized case. In particular, radial gradients in the star and gas distribution, and the physical association of young stars with denser ISM (not just BCs) contribute to differences in the trends. Generally the eagle attenuation is closer to the mixed slab case, and appears to exhibit a similar plateau at high Σdust values. In the mixed 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 |$\tau ^{\rm 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 idealized cases. The continuum measurements tend to occupy similar values as found for eagle over the 5.5 ≤ log10Σ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. NGC 3077 (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 sightlines 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 these lines of sight, or more preferential attenuation of bright young stars than found in eagle. 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 eagle attenuation measurement such that they are closer to a homogeneous mixed 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 BC attenuation, particularly in dusty starbursts such as NGC 2146 (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 the 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 RV in the bottom panel of Fig. 1, the power-law slope of the attenuation curve plateaus to a value close to 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 log10Σ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 mixed slab model. The mixed slab follows a similar qualitative trend, tending towards to shape of intrinsic curve (C00) below log10Σdust/(M⊙ kpc−2) ≈ 5.5, crossing the CF00 threshold at log10Σdust/(M⊙ kpc−2) ≈ 6.8. At log10Σdust/(M⊙ kpc−2) ≳ 6, the mixed slab case is a good approximation to the slope of the eagle attenuation curve. The higher attenuation in eagle relative to the mixed slab is ascribed to the preferential reddening of young stars due to being embedded in preferentially dustier ISM regions, even aside from the effects of stellar BC 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 mixed 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 of |$A^{\rm ISM}_{V}$|⁠, with the 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^{\rm ISM}_{V} \approx 0.8$| mag. The underlying colour map indicates the residual trend between galaxy orientation and ηISM as functions of |$A^{\rm ISM}_{V}$|⁠, with shading indicating the median offset of galaxies from the median inclination angle at that |$A^{\rm ISM}_{V}$|⁠. At a fixed |$A^{\rm 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^{\rm 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. Figure 3. Open in new tabDownload slide Plot of the power-law slope of the best-fitting attenuation curve, ηISM, as a function of V-band ISM attenuation strength, |$A^{\rm ISM}_V$|⁠. The black dashed line shows a cubic spline fit to binned medians, while the black error bars show the 16th–84th percentile range in each bin. The underlying colour map shows the residual inclination offset for galaxies in each |$A^{\rm ISM}_V$| 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^{\rm ISM}_V \ge 0.75$|⁠. For a given attenuation, more edge-on (higher inclination) galaxies also tend to have greyer attenuation curves. Figure 3. Open in new tabDownload slide Plot of the power-law slope of the best-fitting attenuation curve, ηISM, as a function of V-band ISM attenuation strength, |$A^{\rm ISM}_V$|⁠. The black dashed line shows a cubic spline fit to binned medians, while the black error bars show the 16th–84th percentile range in each bin. The underlying colour map shows the residual inclination offset for galaxies in each |$A^{\rm ISM}_V$| 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^{\rm ISM}_V \ge 0.75$|⁠. For a given attenuation, more edge-on (higher inclination) galaxies also tend to have greyer attenuation curves. Together, these results indicate how trends in ISM attenuation properties (computed via RT) 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 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 BC term to enable a better comparison with observations. Despite these limitations, we present this ISM screen model as a more nuanced alternative to idealized geometric models typically used when interpreting observations, incorporating the diverse galaxy morphologies that emerge within the eagle simulations. 5 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. Next, we explore the treatment of the τBC term, and incorporate τBC for comparison to data. 5.1 The contribution of infant stellar populations To gauge the possible impact of differing BC 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 tage < 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 (sSFR). 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 (tage < 10 Myr) are run separately through skirt (without dust transfer) to generate pure stellar 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 BCs or a boosted ISM attenuation are included for the infant stars, the emergent flux fraction will be lower. Figure 4. Open in new tabDownload slide The fraction of the emitted flux emanating from infant stellar populations (with tage < 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 log10(M⋆/M⊙) > 9.5 are stacked for bins of sSFR (see the 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 per cent contribution on the blue side of the g band. Figure 4. Open in new tabDownload slide The fraction of the emitted flux emanating from infant stellar populations (with tage < 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 log10(M⋆/M⊙) > 9.5 are stacked for bins of sSFR (see the 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 per cent contribution on the blue side of the g band. 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 per cent (⪅30 per cent) 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 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 BC treatment may still be important in some regimes. At shorter wavelengths, we see that the infant star fraction rises markedly, up to ≈50 per cent in the NUV. The infant star contribution will also be high for certain atomic transitions that emanate from compact H ii regions across the UV–NIR range. In what follows, we explore how including an additional BC term may influence the effective attenuation in galaxies. 5.2 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 normalization, AV, parametrize 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 AV 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 AV–δ relation appears to reproduce the observations very well. Figure 5. Open in new tabDownload slide The relationship between the best-fittingAV and δ parameters, fitting the modified C00 relation of equation (1). The black colours indicate the ISM-only attenuation, whereas the red points include an additional BC term (see the text for details). Alternate BC prescription are demonstrated by the red lines. Median values are indicated by the solid points, with error bars indicating the 16th−84th percentile scatter. The contours and outlying points of the ISM-only attenuation are plotted in faint grey, to indicate the underlying distribution. For comparison, the observationally inferred relation of 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 BC attenuation improves agreement. Figure 5. Open in new tabDownload slide The relationship between the best-fittingAV and δ parameters, fitting the modified C00 relation of equation (1). The black colours indicate the ISM-only attenuation, whereas the red points include an additional BC term (see the text for details). Alternate BC prescription are demonstrated by the red lines. Median values are indicated by the solid points, with error bars indicating the 16th−84th percentile scatter. The contours and outlying points of the ISM-only attenuation are plotted in faint grey, to indicate the underlying distribution. For comparison, the observationally inferred relation of 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 BC attenuation improves agreement. Steeper attenuation curves are expected from models that include the influence of small scale and BC absorption; the young stellar populations that are associated with BCs 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 H ii 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, and avoids the double-counting of attenuation. eagle cannot provide insight into the nature of true BC attenuation, as these structures are well below the scale that can be resolved. Instead, some assumptions must be made about the nature of BC attenuation. The high optical depths of stellar BCs makes it difficult to determine the wavelength dependence of BC 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 BCs 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 BC 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 |$\eta _{\rm BC}^{\rm shell} = -1.3$|⁠. The other primary uncertainty is in the normalization of the optical depth assumed for the BCs. Typically, BC 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 BC attenuation for eagle galaxies in post-processing without double-counting attenuation on subgrid scales, we use both the total SEDs of eagle galaxies, and the SED contributions of infant stars with ages tage < 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, fnew, we subtract away the MAPPINGS-III contributions and add the galaxev SEDs for infant stars, corrected for the influence of stellar BCs, as $$\begin{eqnarray*} f_{\rm new} = f_{\rm total} - f_{\rm H\,{\small II}} + f^\prime _{\rm BC03}\exp \left(-f_{\tau } \tau _{\rm ISM} \left[\frac{\lambda }{5500\,{\mathring{\rm A} }}\right]^{\eta _{\rm BC}^{\rm shell}}\right), \end{eqnarray*}$$(6) where ftotal is the fiducial skirt flux, fH ii is for MAPPINGS-III only, and |$f^\prime _{\rm 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 δ–AV relationship for the BC-corrected eagle fluxes, assuming |$\eta _{\rm BC}=\eta _{\rm ISM}^{\rm shell}$| and fτ = 2. We see that δ shifts to systematically lower values at a fixed AV, indicating steeper (redder) attenuation curves. The shift decreases gradually with AV, largest in the low AV bin. This can be understood as a result of the higher contribution of infant stars to the total attenuation at low effective AV, due to the dominance of the τBC term. At high AV infant stars are effectively hidden, and, given the limited fraction that they can possibly contribute (Fig. 4), the BC attenuation no longer makes a significant difference to the emergent fluxes. Comparing to the Salim et al. (2018) data, we see that including this BC 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 (AV < 0.5). While the agreement of the ISM-only relation is similar to that of Narayanan et al. (2018, their fig. 11, left-hand 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 parametrizations, 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 that 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 normalization. However, this parametrization (dashed), along with all the aforementioned variations, yields 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 δ−AV relation is largely insensitive to the exact BC 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 that is realistic. This is closely related to resolution effects, as there is missing structure on scales below ∼1 kpc and the BCs 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 BCs. The true properties of BCs 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 BCs in different environments, which remains a caveat for any implementation of such a screen model. Finally, we caution that inferring the δ and AV 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 that the same inference procedure applied to observations could be applied to eagle spectra, yielding a fairer comparison, but this is left for future work. 6 SUMMARY AND CONCLUSIONS In this study, we present a model for effective dust attenuation parametrized 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 CF00 by distilling information gained from full RT of ∼10 000 eagle galaxies in random orientation (via the skirt code Camps et al. 2016; Trayford et al. 2017). This is motivated by first identifying a strong dependence of the non-parametric AV and RV 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 idealized geometries and inferred from observation. Some consideration is made for the treatment of BCs 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 BC 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 AV is both remarkably tight and highly independent of redshift. The Σdust–RV dependence is less tight, but remains strong and demonstrates a similar level of redshift independence (Fig. 1). The best-fitting power-law attenuation curves show optical depths, |$\tau ^{\rm 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 (⪆106 M⊙ kpc−2) and close to a screen (i.e. extinction curve) case at low Σdust (⪅105 M⊙ kpc−2). We demonstrate that the scatter in the ηISM versus AV 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 AV–δ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 attenuation curves similar to that shown in Narayanan et al. 2018 (Fig. 5). We find that by applying a BC 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 BC optical depth and attenuation slope (assuming a fixed BC dispersal time of 10 Myr). Indeed, a number of different parametrizations yield near identical results. The BC parametrization 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 SAMs of galaxy formation, where galaxy structure is unknown or idealized. While RT has previously been assumed for SAMS via idealized 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. ACKNOWLEDGEMENTS JT thanks Joop Schaye for useful discussions during the development of this work, and acknowledges the use of the DiRAC Data Centric system at Durham University, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). JT and CL thank the University of Western Australia for a Research Collaboration Award 2018 which facilitated face-to-face interactions which contributed to this work. CL has received funding from the ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through project number CE170100013 and the Cosmic Dawn Centre, which is funded by the Danish National Research Foundation. Footnotes 1 http://www.skirt.ugent.be 2 i.e. tage < 10 Myr stars, without the influence of nebular emission and absorption built into the MAPPINGS-III library. 3 The property maps produced by Trayford & Schaye (2018) are memory heavy, but were already produced for the four specified redshifts. 4 This is typically ascribed to strong obscuration of the evolved central regions, allowing a higher relative contribution of young stars in outer regions. 5 A caveat is that the intrinsic dust extinction and dust-to-metal ratios are held fixed in the skirt modelling with redshift. 6 δ is an alternative measure of attenuation curve slope, obtained through fitting a modified C00 law (equation 2). 7 See Trayford et al. (2017) for further dicussion of the effects of the puffing-up and homogenization of the eagle ISM on the skirt results. REFERENCES Aniano G. et al. ., 2012 , ApJ , 756 , 138 10.1088/0004-637X/756/2/138 Crossref Search ADS Baes M. , Dejonghe H., 2001 , MNRAS , 326 , 733 10.1046/j.1365-8711.2001.04626.x Crossref Search ADS Baes M. et al. ., 2003 , MNRAS , 343 , 1081 10.1046/j.1365-8711.2003.06770.x Crossref Search ADS Baes M. , Verstappen J., De Looze I., Fritz J., Saftly W., Vidal Pérez E., Stalevski M., Valcke S., 2011 , ApJS , 196 , 22 10.1088/0067-0049/196/2/22 Crossref Search ADS Benítez-Llambay A. , Navarro J. F., Frenk C. S., Ludlow A. D., 2018 , MNRAS , 473 , 1019 10.1093/mnras/stx2420 Crossref Search ADS Bianchi S. , 2008 , A&A , 490 , 461 10.1051/0004-6361:200810027 Crossref Search ADS Boselli A. et al. ., 2010 , PASP , 122 , 261 10.1086/651535 Crossref Search ADS Bruzual G. , Charlot S., 2003 , MNRAS , 344 , 1000 10.1046/j.1365-8711.2003.06897.x Crossref Search ADS Calzetti D. , 2013 , Star Formation Rate Indicators , Cambridge University Press , Cambridge, UK . p. 419 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Calzetti D. , Kinney A. L., Storchi-Bergmann T., 1994 , ApJ , 429 , 582 10.1086/174346 Crossref Search ADS Calzetti D. , Armus L., Bohlin R. C., Kinney A. L., Koornneef J., Storchi-Bergmann T., 2000 , ApJ , 533 , 682 10.1086/308692 Crossref Search ADS Camps P. , Baes M., 2015 , Astron. Comput. , 9 , 20 10.1016/j.ascom.2014.10.004 Crossref Search ADS Camps P. , Trayford J. W., Baes M., Theuns T., Schaller M., Schaye J., 2016 , MNRAS , 462 , 1057 10.1093/mnras/stw1735 Crossref Search ADS Camps P. et al. ., 2018 , ApJS , 234 , 20 10.3847/1538-4365/aaa24c Crossref Search ADS Capak P. L. et al. ., 2015 , Nature , 522 , 455 10.1038/nature14500 Crossref Search ADS PubMed Cardelli J. A. , Clayton G. C., Mathis J. S., 1989 , ApJ , 345 , 245 10.1086/167900 Crossref Search ADS Casey C. M. et al. ., 2014 , ApJ , 796 , 95 10.1088/0004-637X/796/2/95 Crossref Search ADS Chabrier G. , 2003 , PASP , 115 , 763 10.1086/376392 Crossref Search ADS Charlot S. , Fall S. M., 2000 , ApJ , 539 , 718 10.1086/309250 Crossref Search ADS Chevallard J. , Charlot S., Wandelt B., Wild V., 2013 , MNRAS , 432 , 2061 10.1093/mnras/stt523 Crossref Search ADS Conroy C. , 2013 , ARA&A , 51 , 393 10.1146/annurev-astro-082812-141017 Crossref Search ADS Crain R. A. et al. ., 2015 , MNRAS , 450 , 1937 10.1093/mnras/stv725 Crossref Search ADS da Cunha E. , Charlot S., Elbaz D., 2008 , MNRAS , 388 , 1595 10.1111/j.1365-2966.2008.13535.x Crossref Search ADS Dalla Vecchia C. , Schaye J., 2012 , MNRAS , 426 , 140 10.1111/j.1365-2966.2012.21704.x Crossref Search ADS Davé R. , Thompson R., Hopkins P. F., 2016 , MNRAS , 462 , 3265 10.1093/mnras/stw1862 Crossref Search ADS De Looze I. et al. ., 2014 , A&A , 571 , A69 10.1051/0004-6361/201424747 Crossref Search ADS Decleir M. et al. ., 2019 , MNRAS , 486 , 743 10.1093/mnras/stz805 Crossref Search ADS Draine B. T. et al. ., 2007 , ApJ , 663 , 866 10.1086/518306 Crossref Search ADS Feldmann R. , Quataert E., Hopkins P. F., Faucher-Giguère C.-A., Kereš D., 2017 , MNRAS , 470 , 1050 10.1093/mnras/stx1120 Crossref Search ADS Fischera J. , Dopita M. A., Sutherland R. S., 2003 , ApJ , 599 , L21 10.1086/381190 Crossref Search ADS Fitzpatrick E. L. , 1999 , PASP , 111 , 63 10.1086/316293 Crossref Search ADS Fontanot F. , Somerville R. S., Silva L., Monaco P., Skibba R., 2009 , MNRAS , 392 , 553 10.1111/j.1365-2966.2008.14126.x Crossref Search ADS Furlong M. et al. ., 2015 , MNRAS , 450 , 4486 10.1093/mnras/stv852 Crossref Search ADS Furlong M. et al. ., 2017 , MNRAS , 465 , 722 10.1093/mnras/stw2740 Crossref Search ADS Gonzalez-Perez V. , Lacey C. G., Baugh C. M., Frenk C. S., Wilkins S. M., 2013 , MNRAS , 429 , 1609 10.1093/mnras/sts446 Crossref Search ADS Gordon K. D. , Calzetti D., Witt A. N., 1997 , ApJ , 487 , 625 10.1086/304654 Crossref Search ADS Groves B. , Dopita M. A., Sutherland R. S., Kewley L. J., Fischera J., Leitherer C., Brandl B., van Breugel W., 2008 , ApJS , 176 , 438 10.1086/528711 Crossref Search ADS Guo Q. et al. ., 2016 , MNRAS , 461 , 3457 10.1093/mnras/stw1525 Crossref Search ADS Haardt F. , Madau P., 2001 , in Neumann D. M., Tran J. T. V., eds, Clusters of Galaxies and the High Redshift Universe Observed in X-rays , Editions Frontieres , Paris . p. 64 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Holwerda B. W. , Keel W. C., 2017 , in Gil de Paz A., Knapen J. H., Lee J. C., eds, Proc. IAU Symp. 321 , Formation and Evolution of Galaxy Outskirts . Kluwer , Dordrecht , p. 248 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Hopkins P. F. , 2013 , MNRAS , 428 , 2840 10.1093/mnras/sts210 Crossref Search ADS Hopkins A. M. , Connolly A. J., Haarsma D. B., Cram L. E., 2001 , AJ , 122 , 288 10.1086/321113 Crossref Search ADS Inoue A. K. , Kamaya H., 2004 , MNRAS , 350 , 729 10.1111/j.1365-2966.2004.07686.x Crossref Search ADS Jackson J. M. , Ho P. T. P., 1988 , ApJ , 324 , L5 10.1086/185079 Crossref Search ADS Jonsson P. , Groves B., Cox T. J., 2010 , MNRAS , 403 , 17 Crossref Search ADS Keel W. C. , van Soest E. T. M., 1992 , A&AS , 94 , 553 Kreckel K. et al. ., 2013 , ApJ , 771 , 62 10.1088/0004-637X/771/1/62 Crossref Search ADS Kriek M. , Conroy C., 2013 , ApJ , 775 , L16 10.1088/2041-8205/775/1/L16 Crossref Search ADS Lagos C. d. P. et al. ., 2015 , MNRAS , 452 , 3815 10.1093/mnras/stv1488 Crossref Search ADS Lagos C. d. P. , Tobar R. J., Robotham A. S. G., Obreschkow D., Mitchell P. D., Power C., Elahi P. J., 2018 , MNRAS , 481 , 3573 10.1093/mnras/sty2440 Crossref Search ADS Lagos C. d. P. et al. ., 2019 , MNRAS , 489 , 4196 Crossref Search ADS Leja J. et al. ., 2019 , ApJ , 877 , 140 10.3847/1538-4357/ab1d5a Crossref Search ADS Mattsson L. , De Cia A., Andersen A. C., Zafar T., 2014 , MNRAS , 440 , 1562 10.1093/mnras/stu370 Crossref Search ADS Meurer G. R. , Heckman T. M., Calzetti D., 1999 , ApJ , 521 , 64 10.1086/307523 Crossref Search ADS Narayanan D. , Conroy C., Davé R., Johnson B. D., Popping G., 2018 , ApJ , 869 , 70 10.3847/1538-4357/aaed25 Crossref Search ADS Noll S. , Burgarella D., Giovannoli E., Buat V., Marcillac D., Muñoz-Mateos J. C., 2009 , A&A , 507 , 1793 10.1051/0004-6361/200912497 Crossref Search ADS Padmanabhan N. et al. ., 2008 , ApJ , 674 , 1217 10.1086/524677 Crossref Search ADS Planck Collaboration XX , 2014 , A&A , 571 , A20 10.1051/0004-6361/201321521 Crossref Search ADS Reddy N. A. , Erb D. K., Pettini M., Steidel C. C., Shapley A. E., 2010 , ApJ , 712 , 1070 10.1088/0004-637X/712/2/1070 Crossref Search ADS Rodriguez-Gomez V. et al. ., 2019 , MNRAS , 483 , 4140 10.1093/mnras/sty3345 Crossref Search ADS Rosas-Guevara Y. M. et al. ., 2015 , MNRAS , 454 , 1038 10.1093/mnras/stv2056 Crossref Search ADS Saftly W. , Baes M., De Geyter G., Camps P., Renaud F., Guedes J., De Looze I., 2015 , A&A , 576 , A31 10.1051/0004-6361/201425445 Crossref Search ADS Salim S. , Boquien M., Lee J. C., 2018 , ApJ , 859 , 11 10.3847/1538-4357/aabf3c Crossref Search ADS Salmon B. et al. ., 2016 , ApJ , 827 , 20 10.3847/0004-637X/827/1/20 Crossref Search ADS Sasseen T. P. , Hurwitz M., Dixon W. V., Airieau S., 2002 , ApJ , 566 , 267 10.1086/337955 Crossref Search ADS Schaye J. , 2004 , ApJ , 609 , 667 10.1086/421232 Crossref Search ADS Schaye J. et al. ., 2015 , MNRAS , 446 , 521 10.1093/mnras/stu2058 Crossref Search ADS Springel V. , 2005 , MNRAS , 364 , 1105 10.1111/j.1365-2966.2005.09655.x Crossref Search ADS Steinacker J. , Baes M., Gordon K. D., 2013 , ARA&A , 51 , 63 10.1146/annurev-astro-082812-141042 Crossref Search ADS Sullivan M. , Mobasher B., Chan B., Cram L., Ellis R., Treyer M., Hopkins A., 2001 , ApJ , 558 , 72 10.1086/322451 Crossref Search ADS Taylor E. N. et al. ., 2011 , MNRAS , 418 , 1587 10.1111/j.1365-2966.2011.19536.x Crossref Search ADS Trayford J. W. et al. ., 2015 , MNRAS , 452 , 2879 10.1093/mnras/stv1461 Crossref Search ADS Trayford J. W. et al. ., 2017 , MNRAS , 470 , 771 10.1093/mnras/stx1051 Crossref Search ADS Trayford J. W. , Schaye J., 2018 , MNRAS , 485 , 5715 Crossref Search ADS Tress M. , Ferreras I., Pérez-González P. G., Bressan A., Barro G., Domínguez-Sánchez H., Eliche-Moral C., 2019 , MNRAS , 488 , 1799 Crossref Search ADS Viaene S. et al. ., 2017 , A&A , 599 , A64 10.1051/0004-6361/201629251 Crossref Search ADS Walcher J. , Groves B., Budavári T., Dale D., 2011 , Ap&SS , 331 , 1 10.1007/s10509-010-0458-z Crossref Search ADS Whitney B. A. , 2011 , Bull. Astron. Soc. India , 39 , 101 Wiersma R. P. C. , Schaye J., Smith B. D., 2009a , MNRAS , 393 , 99 10.1111/j.1365-2966.2008.14191.x Crossref Search ADS Wiersma R. P. C. , Schaye J., Theuns T., Dalla Vecchia C., Tornatore L., 2009b , MNRAS , 399 , 574 10.1111/j.1365-2966.2009.15331.x Crossref Search ADS Wild V. , Kauffmann G., Heckman T., Charlot S., Lemson G., Brinchmann J., Reichard T., Pasquali A., 2007 , MNRAS , 381 , 543 10.1111/j.1365-2966.2007.12256.x Crossref Search ADS Wild V. , Charlot S., Brinchmann J., Heckman T., Vince O., Pacifici C., Chevallard J., 2011 , MNRAS , 417 , 1760 10.1111/j.1365-2966.2011.19367.x Crossref Search ADS Witt A. N. , Gordon K. D., 2000 , ApJ , 528 , 799 10.1086/308197 Crossref Search ADS Witt A. N. , Thronson Jr. H. A., Capuano Jr. J. M., 1992 , ApJ , 393 , 611 10.1086/171530 Crossref Search ADS Wuyts S. et al. ., 2009 , ApJ , 700 , 799 10.1088/0004-637X/700/1/799 Crossref Search ADS Zibetti S. , Charlot S., Rix H.-W., 2010 , in Bruzual G. R., Charlot S., eds, Proc. IAU Symp. 262 , Stellar Populations – Planning for the Next Decade , Kluwer , Dordrecht , p. 89 10.1017/S1743921310002589 Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Crossref Zubko V. , Dwek E., Arendt R. G., 2004 , ApJS , 152 , 211 10.1086/382351 Crossref Search ADS APPENDIX A: RESOLUTION & CONVERGENCE 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 super kpc scale structures calculated for fiducial eagle galaxies using the skirt code. The BC term may then be attributed to unresolved structures on the scale of H ii 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. 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 mg, Plummer equivalent gravitational softening ϵprop at redshift z = 0 in proper kpc (pkpc), and the study reference for each volume. Name . L . mg . ϵprop . Ref. . . cMpc . 105 M⊙ . pkpc . . RefL025N0376 (Ref25) 25 18.1 0.70 S15 RefL025N0752 (RefHi25) 25 2.26 0.35 S15 RecalL025N0752 (Recal25) 25 2.26 0.35 S15 RefL100N1504 (Ref100) 100 18.1 0.70 S15 Name . L . mg . ϵprop . Ref. . . cMpc . 105 M⊙ . pkpc . . RefL025N0376 (Ref25) 25 18.1 0.70 S15 RefL025N0752 (RefHi25) 25 2.26 0.35 S15 RecalL025N0752 (Recal25) 25 2.26 0.35 S15 RefL100N1504 (Ref100) 100 18.1 0.70 S15 Open in new tab 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 mg, Plummer equivalent gravitational softening ϵprop at redshift z = 0 in proper kpc (pkpc), and the study reference for each volume. Name . L . mg . ϵprop . Ref. . . cMpc . 105 M⊙ . pkpc . . RefL025N0376 (Ref25) 25 18.1 0.70 S15 RefL025N0752 (RefHi25) 25 2.26 0.35 S15 RecalL025N0752 (Recal25) 25 2.26 0.35 S15 RefL100N1504 (Ref100) 100 18.1 0.70 S15 Name . L . mg . ϵprop . Ref. . . cMpc . 105 M⊙ . pkpc . . RefL025N0376 (Ref25) 25 18.1 0.70 S15 RefL025N0752 (RefHi25) 25 2.26 0.35 S15 RecalL025N0752 (Recal25) 25 2.26 0.35 S15 RefL100N1504 (Ref100) 100 18.1 0.70 S15 Open in new tab The Ref25 simulation uses the same resolution and physics model as Ref100, but realizes a 253 cMpc volume. The Rec25 and RefHi25 simulations are initialized for the same volume and initial conditions, but with enhanced particle mass resolution of mg ⪆ 2.26 × 106 M⊙ (mDM = 1.21 × 106 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-fitting power-law ISM attenuation curve for each simulated galaxy, previously seen in Fig. 2. The left-hand panel shows that the relatively tight relationship between V-band (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 |$\tau ^{\rm ISM}_{550}$|–log10Σ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 ⪅ log10Σdust/(M⊙ kpc−2) ⪅6.5 range sampled by the 25 cMpc simulations. Figure A1. Open in new tabDownload slide Convergence properties of the best-fitting V-band ISM optical depth, |$\tau ^{\rm ISM}_{550}$|⁠, and power-law slope, ηISM, as functions of Σdust, originally shown in Fig. 2. The left- and right-hand panels plot the Σdust dependence of |$\tau ^{\rm 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. Figure A1. Open in new tabDownload slide Convergence properties of the best-fitting V-band ISM optical depth, |$\tau ^{\rm ISM}_{550}$|⁠, and power-law slope, ηISM, as functions of Σdust, originally shown in Fig. 2. The left- and right-hand panels plot the Σdust dependence of |$\tau ^{\rm 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. 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 (log10Σ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 theorized 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 log10Σ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 affect the integrated attenuation, but are limited by insufficient small-scale physics, lacking molecular gas modelling and harbouring an artificially pressurized 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 idealized geometric models, which begins 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 normalization and slope for the ugriz range. The treatment of attenuation associated with sub-resolution BCs 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 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 < Re) and Σdust(r < 2Re). 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 centres of galaxies. The Σdust(r < Re) values are systematically lower on average, shifting the relations to ≈0.25 dex lower Σdust. The larger aperture Σdust(r < Re) yields values that are lower still, shifting the relations by a further ≈0.35 dex in Σdust. Figure B1. Open in new tabDownload slide 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 < Re) and Σdust(r < 2Re), respectively. We see that 〈Σdust〉 produces the highest surface densities on average, followed by Σdust(r < Re) and Σdust(r < 2Re) producing lower surface densities, with an overall ∼0.6 dex shift between the median relations. Figure B1. Open in new tabDownload slide 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 < Re) and Σdust(r < 2Re), respectively. We see that 〈Σdust〉 produces the highest surface densities on average, followed by Σdust(r < Re) and Σdust(r < 2Re) producing lower surface densities, with an overall ∼0.6 dex shift between the median relations. 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 < Re) as our fiducial measure.   APPENDIX C: FIT QUALITY While non-parametric attenuation properties are presented in Fig. 1, most of our analysis is presented using best-fitting parameters to individual eagle attenuation curves, via equations (1) and (2). Fig. C1 demonstrates the quality of these fits via the reduced chi-squared parameter, |$\chi ^2_{\rm red}$|⁠. This is defined in the standard way as the sum of the uncertainty-scaled squared offset between the measured and fit functional forms, divided by the number of degrees of freedom in the fit. We note that makes use of the representative photometric errors for the SDSS survey (Padmanabhan et al. 2008) rather than the intrinsic errors in the skirt attenuation curves, which are typically much smaller for AV > 0.25. Figure C1. Open in new tabDownload slide The main panel shows the histogram of goodness of fit, via a reduced chi squared (⁠|$\chi ^2_{\rm red}$|⁠) measure, for the power-law fit to the ISM attenuation of eagle galaxies with AV > 0.25. We note that |$\chi ^2_{\rm 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 |$\chi ^2_{\rm red}$| distribution are highlighted in green, amber, and red, representing three levels of decreasing fit quality. Two attenuation curves are randomly sampled from each region and are plotted in the inset panel (solid lines), along with the best-fitting CF00 power law (equation 2, fine dotted lines) and modified C00 law (equation 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. Figure C1. Open in new tabDownload slide The main panel shows the histogram of goodness of fit, via a reduced chi squared (⁠|$\chi ^2_{\rm red}$|⁠) measure, for the power-law fit to the ISM attenuation of eagle galaxies with AV > 0.25. We note that |$\chi ^2_{\rm 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 |$\chi ^2_{\rm red}$| distribution are highlighted in green, amber, and red, representing three levels of decreasing fit quality. Two attenuation curves are randomly sampled from each region and are plotted in the inset panel (solid lines), along with the best-fitting CF00 power law (equation 2, fine dotted lines) and modified C00 law (equation 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. Generally, we see that the majority of galaxies display a |$\chi ^2_{\rm red} \lesssim 1$|⁠, with a median of |$\bar{\chi }^2_{\rm red} = 0.65$|⁠. The distribution has a pronounced tail to high |${\chi }^2_{\rm red}$| values, and is well represented by a lognormal distribution. Inset, we show two randomly selected attenuation curves from each of the three highlighted |${\chi }^2_{\rm red}$| ranges in the main panel. For each curve, the best-fitting CF00 power law and modified C00 law appear to capture the shape of the attenuation curves generally well. The higher (highest) |${\chi }^2_{\rm red}$| bin examples, shown in amber (red), typically display a marginally better fit for the single power-law shape, with discrepancies between measurements and fits mainly arising from a sharper uptick in the u-band attenuation relative to g band than either model can capture. This is ascribed to the u band being blueward of the 4000Å break, and as such displays a discontinuous increase in the contribution by flux from young stars relative to redder bands. We emphasize that this is not indicative of a sudden change in slope towards bluer wavelengths, but rather a discontinuity in the attenuation curves to slightly higher value. This effect can be seen explicitly in fig. 6 of Trayford et al. (2017). Together, this demonstrates that eagle attenuation curves are generally well fit by the relations given by equations (1) and (2), suggesting that the parameters used in Figs 2, 3, and 5 are representative of the eagle attenuation curves. © 2019 The Author(s). Published by Oxford University Press on behalf of The Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. © 2019 The Author(s). Published by Oxford University Press on behalf of The Royal Astronomical Society.
TI - Fade to grey: systematic variation of galaxy attenuation curves with galaxy properties in the eagle simulations
JF - Monthly Notices of the Royal Astronomical Society
DO - 10.1093/mnras/stz3234
DA - 2020-01-21
UR - https://www.deepdyve.com/lp/oxford-university-press/fade-to-grey-systematic-variation-of-galaxy-attenuation-curves-with-Nqa0lYyqDv
SP - 3937
EP - 3951
VL - 491
IS - 3
DP - DeepDyve
ER -