A Theory for the Variation of Dust Attenuation Laws in Galaxies

Desika Narayanan; Charlie Conroy; Romeel Davé; Benjamin D. Johnson; Gergö Popping

doi:10.3847/1538-4357/aaed25

1. Introduction

Astrophysical dust is pervasive in the interstellar medium (ISM) of most galaxies. This dust can both absorb stellar radiation and scatter light into and out of the line of sight (Witt & Gordon 2000). Physical properties of galaxies that derive from optical and ultraviolet (UV) radiation, such as the star formation rate (SFR), stellar mass (M_*), and stellar ages, must therefore account for these dust-driven effects. Unfortunately, how exactly to do this is complicated (see reviews by Walcher et al. 2011; Conroy 2013). Beyond the aforementioned effects of extinction and scattering, the spatial distribution between stars of different ages and the dusty ISM (hereafter simply referred to as "the geometry") can dramatically impact dust attenuation⁹ laws in galaxies. These effects are typically described via an attenuation curve, which quantifies the optical depth (τ) as either a function of wavelength (λ) or x ≡ 1/λ.

Decades of constraints of observed extinction attenuation curves have evidenced a number of clear features (see reviews by Calzetti 2001; Draine 2003; Galliano et al. 2018). First, there is a steep rise toward the UV due to absorption by small grains such that shorter wavelength photons are preferentially removed from the line of sight. In the near-UV (NUV), there is (sometimes) an absorption bump near 2175 Å, first reported by Stecher (1965), that is potentially associated with polycyclic aromatic hydrocarbons (e.g., Weingartner & Draine 2001), though some models suggest a graphite-based origin (e.g., Hoyle & Wickramasinghe 1962; Stecher & Donn 1965; Draine & Malhotra 1993) or other species (e.g., Wada et al. 1999). This 2175 Å UV bump is known to vary dramatically in strength, a topic we return to later. There is some evidence for an optical knee due to scattering by large grains (Galliano et al. 2018), and the extinction curve takes a power-law shape at longer NIR wavelengths. In Figure 1, we show a compilation of a number of literature extinction and attenuation curves that will serve as a useful reference throughout this paper.

**Figure 1.** Compilation of selected extinction and attenuation curve constraints taken from the literature. The extinction curves are the Cardelli et al. (1989) and Pei (1992) curves, while the remainder are attenuation curves. The solid lines show reference curves from the literature, while the green dashed line shows the intrinsic extinction curve used in our cosmological galaxy formation simulations in this paper (Weingartner & Draine 2001).
Download figure:
Standard image High-resolution image

These observations have revealed a diversity in derived extinction curves in the Milky Way and nearby galaxies. For example, when parameterizing the inverse of the slope of the extinction curve in the optical as ${R}_{{\rm{V}}}\equiv {A}_{{\rm{V}}}/({A}_{{\rm{B}}}-{A}_{{\rm{V}}})$ , observations have ranged from values of R_V as low as ∼2 (Welty & Fowler 1992) to nearly R_V ∼ 6 (Cardelli et al. 1989; Fitzpatrick 1999). Similarly, Fitzpatrick & Massa (1990) demonstrated a broad range of extinction curves among different Milky Way sight lines via International Ultraviolet Explorer (IUE) observations. Cardelli et al. (1989) and Fitzpatrick & Massa (1990) showed that variations in both the attenuation curve shape and 2175 Å UV bump strength can be well parameterized by the total to selective extinction R_V. The sight-line average attenuation curve of the Milky Way has R_V ≈ 3.1 (e.g., Rieke & Lebofsky 1985). Outside of our galaxy, extinction curves have been measured using this method for the Magellanic clouds and Andromeda. Pei (1992) and Gordon et al. (2003) derived attenuation curves for the Magellanic clouds and showed that the Small Magellanic Cloud (SMC) bar, for example, has a steeper curve than the Large Magellanic Cloud average when normalizing at optical wavelengths. On average, both of the Magellanic clouds have lower R_V ratios than the Milky Way, which may be due to their lower metallicities (Misselt et al. 1999; Calzetti 2001). At the same time, Clayton et al. (2015) found that M31 has a relatively similar average curve as the Milky Way.

At larger distances, emission from unresolved stellar populations is measured in aggregate (as opposed to in individual stellar pairs), and the role of star-to-dust geometry can introduce significant complications. In this regime, constraints on extinction curves become impossible, and instead, measurements are made of attenuation curves. When measuring attenuation toward stars in unresolved galaxies (i.e., A_star), methods include utilizing the ratio of the infrared excess (itself a ratio between the infrared luminosity to a monochromatic UV luminosity) to the UV slope (i.e., the IRX-β relation; Meurer et al. 1999; Johnson et al. 2007b; Casey et al. 2014b), color–magnitude diagram fitting (e.g., Dalcanton et al. 2015), and spectral energy distribution (SED) fitting techniques (e.g., Papovich et al. 2001; Buat et al. 2011; Conroy 2013; Kriek & Conroy 2013; Leja et al. 2017; Salim et al. 2018).

Akin to extinction curves within the Milky Way, attenuation curves in galaxies near and far have shown a dramatic range in observed shapes. In a series of papers, Calzetti et al. (1994), Kinney et al. (1994), Calzetti (1997), Böker et al. (1999), and Calzetti et al. (2000) investigated attenuation laws toward both stars and H ii regions in local universe starburst galaxies. These researchers found that these galaxies, on average, have grayer (also described as "flatter" or "shallower" throughout this paper) optical/UV attenuation slopes (when normalized in the optical) than the Milky Way, with an average R_V = 4.05. Johnson et al. (2007a) derived a mean attenuation law in a sample of ∼1000 galaxies of A_λ ∼ λ^−0.7, with little variation when binned by stellar mass. Wild et al. (2011) analyzed ∼23,000 galaxies from the Sloan Digital Sky Survey and found that the slope of the attenuation curves varies strongly with galaxy axial ratio, though only weakly with their specific SFR. The slopes of these curves vary with stellar mass surface density, with higher stellar mass surface densities correlating with steeper curves. Meanwhile, Battisti et al. (2016) found, for ∼10,000 star-forming galaxies at z < 0.1, an average attenuation curve similar to that derived by Calzetti et al. (2000), with an observed range comparable to the factor of ∼2 dispersion in τ_FUV/τ_V found by Wild et al. (2011). Battisti et al. (2017) expanded on this study and found little correlation between these curves and a range of galaxy physical properties, including mean stellar age, specific SFR, stellar mass, and metallicity. The consensus from local galaxy observations is that, aside from potentially stellar mass surface density (e.g., Wild et al. 2011), there appear to be relatively few correlations between the slope of dust attenuation curves and galaxy physical property.

At high redshift, the most commonly assumed attenuation curve is a Calzetti et al. (2000) law. Owing to the sensitivity of derived physical parameters on the assumed attenuation law in SED modeling (e.g., Salim et al. 2007, 2016), a number of authors have attempted to constrain observed attenuation laws at high redshift to assess the validity of the typical assumption of a Calzetti et al. (2000) relation. For example, Reddy & Steidel (2004) compared X-ray and radio (i.e., relatively extinction-free) derived SFRs with that derived from Calzetti-corrected UV photometry of UV-selected star-forming galaxies in the range 1.5 ≲ z ≲ 3 in GOODS-N and found consistent results. This implies, on average, that an attenuation curve similar to that derived by Calzetti et al. (2000) is a reasonable descriptor of the attenuation properties of these galaxies. A number of other works, using complementary methods, have found reasonable consistency with a Calzetti et al. (2000) attenuation curve at high redshift (e.g., Scoville et al. 2015; Cullen et al. 2017, 2018). Much of this work has focused on the location of galaxies at high redshift on the IRX–β plane and their consistency with a Calzetti et al. (2000) attenuation curve (e.g., Seibert et al. 2002; Reddy et al. 2008, 2012; Pannella et al. 2009; Siana et al. 2009; Heinis et al. 2013; To et al. 2014; Bourne et al. 2017; McLure et al. 2018).

At the same time, Reddy et al. (2015) derived attenuation curves for ∼200 z ∼ 2 galaxies from the MOSFIRE Deep Evolution Field survey and found evidence for a composite attenuation curve that is similar to the Calzetti et al. (2000) form at short wavelengths and the SMC extinction curve at λ ≳ 2500 Å. Similarly, Lo Faro et al. (2017) found evidence for a composite attenuation curve in a sample of z ∼ 2 galaxies selected for their infrared luminosity. Salmon et al. (2016) found a diverse range of attenuation laws in z ∼ 1.5–3 galaxies from the CANDELS survey with little correlation with galaxy physical property, and Shivaei et al. (2015) required curves steeper than Calzetti et al. (2000) to match the observed IR/UV ratios in ∼250 z ∼ 2 star-forming galaxies. And, in the same vein, dramatic outliers from the Meurer et al. (1999) IRX–β relation at high redshift have pointed to underlying variations in attenuation laws. This said, other complicating factors may contribute to deviations from the locally calibrated Meurer et al. (1999) relation (e.g., Bell et al. 2002; Kong et al. 2004; Grasha et al. 2013; Koprowski et al. 2016; Popping et al. 2017a; Safarzadeh et al. 2017; Narayanan et al. 2018).

Similar to variations in the attenuation curve slope, galaxies in both the local universe and at high redshift exhibit a wide range of 2175 Å UV absorption bump strengths. Both the SMC and average attenuation curve for nearby starburst galaxies lack strong bump features (Calzetti et al. 2000; Gordon et al. 2000), and Gordon et al. (1999) found relatively small bump strengths in attenuation curves for a sample of high-redshift galaxies from the Hubble Deep Fields. Meanwhile, a number of studies at both low redshift (Conroy et al. 2010a; Wild et al. 2011; Battisti et al. 2016; Salim et al. 2018) and high z (Motta et al. 2002; Burgarella et al. 2005; York et al. 2006; Noll et al. 2007, 2009; Stratta et al. 2007; Elíasdóttir et al. 2009; Buat et al. 2011; Scoville et al. 2015) have found evidence of 2175 Å absorption bumps of varying strengths. A substantial step forward was made by Kriek & Conroy (2013), who demonstrated not only evidence of a bump in the attenuation curves of z ∼ 2 star-forming galaxies but also that the bump strength varies with the slope of the attenuation curve (such that the steepest curves have the strongest bumps). Further evidence for this trend at z ∼ 2 was given by Tress et al. (2018).

Given the strong dependence on inferred galaxy physical properties on attenuation curves, understanding the origin of their shape variations is critical. Here modeling can provide significant insight as to the physical drivers of variations in both attenuation law slopes and feature strengths. Models designed to understand attenuation laws in galaxies generally fall into one of three methodology camps. The core of each involves radiative transfer modeling of stellar light through a dusty ISM, and the primary difference lies in how the geometry of stars and the dust is modeled. The general trade-off is that increases in model complexity are typically associated with increased astrophysical realism but decreased ability to develop controlled numerical experiments.

The most simplified models typically involve illuminating a slab or shell-like geometry with a stellar radiation source. Early works include Witt & Gordon (1996, 2000), who used these models to demonstrate that increasingly mixed star-to-dust geometries result in flatter (grayer) attenuation curves. Seon & Draine (2016) expanded upon these significantly by developing a model for a turbulent ISM with a lognormal density profile and a clumpy medium to investigate the origin of attenuation law shape and bump strength variations. While these models can neatly isolate individual physical effects, they do not model galaxies as a whole, taking into account the radiative transfer from stars with a distribution of stellar ages and metallicities through a cosmologically evolved ISM.

The latter two modeling categories analyze attenuation laws for galaxies as a whole. These include hydrodynamic simulations of idealized galaxies evolving without a cosmological context and simulations that employ semi-analytic modeling techniques within a cosmological framework. In the former category, Hayward & Smith (2015), Hou et al. (2017), Jonsson (2006), Natale et al. (2015) and Rocha et al. (2008) coupled idealized models of isolated disk galaxies and galaxy mergers with 3D dust radiative transfer to model their dust attenuation properties. Such simulations can achieve very high resolution but do not include the hierarchical growth of galaxies that can result in more diverse (and realistic) morphologies (e.g., Abruzzo et al. 2018) leading to variations in dust attenuation. This is especially true at early epochs when dusty disks are thicker and clumpier and are responding dynamically to high rates of gas inflows and outflows. Meanwhile, semi-analytic models (SAMs) have been developed to model dust attenuation in a cosmological context (e.g., Granato et al. 2000; Fontanot et al. 2009; Fontanot & Somerville 2011; Wilkins et al. 2012; Gonzalez-Perez et al. 2013). Such models do not track baryonic growth directly but rather track dark matter growth and make physically motivated but simplistic assumptions about the resulting galaxy properties. The benefit is that, when coupled with spectrophotometric dust radiative transfer calculations, SAMs are able to efficiently produce cosmological volumes of galaxies with SED information. But since the baryonic matter is characterized through analytic expressions (e.g., Popescu et al. 2011), the geometry of the ISM and stars is not directly predicted but rather based on simplified assumptions.

What is missing thus far is a theoretical interpretation of dust attenuation laws in the context of models that directly track the hydrodynamic growth and evolution of galaxies from cosmological conditions with sufficient resolution to model the impact of complex star-to-dust geometries. Cosmological hydrodynamic simulations are attractive in their ability to directly track the hierarchical growth of galaxies and therefore provide relatively realistic star-to-dust geometries for stellar populations with a distribution of metallicities and formation ages. In this paper, we present the first ever theoretical model for dust attenuation curves in galaxies from cosmological hydrodynamic simulations.

To do this, we will employ the cosmological "zoom" technique, where we focus on individual halos (and their associated baryons) in cosmological simulations at very high resolution. We couple these simulations with three-dimensional dust radiative transfer to model how the intrinsic stellar light is absorbed and scattered and consequently develop a model for dust attenuation laws in galaxies. Critical to the interpretation of our results: in this model, we will assume an underlying dust extinction law (i.e., we will hold the properties of dust grains fixed) and ask how geometry and radiative transfer effects drive variations in the attenuation law slope and UV bump strength. In Section 2, we describe our simulation setup; in Section 3, we build physical insight via simplified population synthesis models; and in Section 4, we apply this insight to direct results from cosmological zoom galaxy formation simulations. Here we explore the origin of variations in the slope and bump strength in attenuation laws in galaxies. In Section 5, we provide discussion, and we summarize in Section 6. Throughout, we assume a cosmology of $({{\rm{\Omega }}}_{0},{{\rm{\Omega }}}_{{\rm{\Lambda }}},{{\rm{\Omega }}}_{{\rm{b}}},h$ ) = (0.3, 0.7, 0.048, 0.68).

2. Simulation Methodology

2.1. Overview

In order to model dust attenuation curves, we must know both the intrinsic stellar spectrum in model galaxies and the observed SED. To generate these, we first simulate a series of model galaxies using the cosmological zoom technique. We then determine the stellar SEDs for these model galaxies based on the individual metallicities and ages for the star particles via fsps calculations. These stellar SEDs are propagated through the dusty ISM via hyperion dust radiative transfer (wrapped in the powderday code package) and compared to the final observed attenuated UV–optical SED in order to determine the attenuation curve. The underlying extinction properties of the grains are fixed throughout, and we utilize the aforementioned simulations to understand the impact of dust geometry and radiative transfer effects on the observed attenuation curves.

2.2. Cosmological Zoom Galaxy Formation Simulations

Our model galaxy suite is generated from the mufasa zoom simulation series, which are zoomed galaxies from the mufasa cosmological simulation (Davé et al. 2016, 2017a, 2017b). This zoom simulation methodology has been described in detail in Olsen et al. (2017), Narayanan et al. (2018), Abruzzo et al. (2018), and Privon et al. (2018). We therefore refer the reader to those works for details and summarize the salient points here.

All of our simulations are conducted with a modified version of the gizmo hydrodynamic code (Hopkins 2015, 2017), which builds off of the code base in gadget-3 (Springel et al. 2005). We first simulate a coarse-resolution dark matter–only run from z = 249 down to z = 0 using initial conditions generated from music (Hahn & Abel 2011). The initial conditions for this dark matter box are the exact same as those used in the mufasa cosmological hydrodynamic simulation. This coarse-resolution run is done in a 50h⁻¹ Mpc volume and includes 512³ dark matter particles, resulting in a dark matter mass resolution of $7.8\times {10}^{8}{h}^{-1}\,{M}_{\odot }$ . It is from this dark matter–only simulation that we select our model halos to resimulate at much higher resolution (and with baryons included).

We resimulate eight halos over a broad mass range to final redshifts z_final, where z_final varies based on the halo mass. Five of these halos are described in Narayanan et al. (2018), and we have include three additional models in this paper. In Table 1, we describe their physical properties¹⁰ and present their location in stellar mass–halo mass and SFR–M_* space in Abruzzo et al. (2018). These simulations span roughly three decades in halo mass and represent galaxies ranging from Milky Way mass (at z = 0) through halo masses comparable to luminous dusty galaxies at z ∼ 2 (e.g., Casey et al. 2014a; Narayanan et al. 2015).

Table 1. Description of the Simulated Galaxies

Name	${M}_{* ,\mathrm{central}}$	${M}_{* ,\mathrm{halo}}$	M_DM	oct_min	${R}_{* ,\mathrm{half}}$	${R}_{\mathrm{gas},\mathrm{half}}$	z_final
	${M}_{\odot };z=2$	${M}_{\odot };z=2$	${M}_{\odot };z=2$	pc	kpc	kpc
mz0	8.4 × 10¹⁰	4.1 × 10¹¹	4.1 × 10¹³	6.1	5.4	1.3	2.15
mz5	6.9 × 10¹⁰	8.3 × 10¹¹	6.3 × 10¹³	6.1	3.8	2.9	2
mz45	1.3 × 10¹⁰	1.3 × 10¹¹	3.7 × 10¹³	6.1	6.0	5.2	2
mz10	6.8 × 10¹⁰	1.6 × 10¹¹	1.1 × 10¹³	6.1	9.4	14.0	2
z0mz352	2.7 × 10⁹	7.2 × 10⁹	9.2 × 10¹¹	12.2	8.4	9.4	0
z0mz401	2.4 × 10⁹	3.8 × 10⁹	5.8 × 10¹¹	6.1	13.9	20.5	0
z0mz287	1.0 × 10⁹	2.3 × 10⁹	2.9 × 10¹¹	12.2	6.7	2.0	0.65
z0mz374	1.0 × 10⁸	2.9 × 10⁸	1.8 × 10¹¹	12.2	4.8	4.7	0.4

Note. We present all masses at z = 2, regardless of what the final simulation redshift is, in order to facilitate comparison between the models. These are ordered by either their actual or expected z = 2 halo mass. Here oct_min refers to the minimum cell size in the constructed octree for the radiative transfer, while ${R}_{* ,\mathrm{half}}$ and ${R}_{\mathrm{gas},\mathrm{half}}$ refer to the stellar and gas half-mass radii measured over 4π sr.

Download table as: ASCII Typeset image

We identify halos using caesar (Thompson 2014). For each halo of interest for resimulation, we build an ellipsoidal mask around all particles encapsulating a radius 2.5× the distance of the farthest particle from halo center and define this as the Lagrangian region for resimulation. We then track these particles back to z = 249, split these to the desired resolution, and rerun the simulation to z = z_final with hydrodynamics turned on. We have zero low-resolution particles contaminating any halo presented here within three virial radii.

Our simulations use the suite of physics developed for the mufasa cosmological hydrodynamic simulations (Davé et al. 2016, 2017a, 2017b). Here stars form in dense molecular gas, and the H₂ fraction is calculated using the Krumholz et al. (2009) methodology, tying the molecular fraction to the gas surface density and metallicity (Thompson et al. 2014). We impose a minimum metallicity for star formation of Z = 10⁻³ Z_⊙. This star formation occurs at a rate following a volumetric Schmidt (1959) relation, with an imposed efficiency per freefall time of _* = 0.02, as motivated by observations (Kennicutt 1998; Narayanan et al. 2008, 2012; Kennicutt & Evans 2012; Hopkins et al. 2013).

Feedback from massive stars is modeled via the mufasa decoupled two-phase wind scheme. The wind physics are presented in detail in Davé et al. (2016), and we refer the reader to that work for a full description of the wind model. Briefly, the modeled stellar winds have a probability for ejection that is modeled as a fraction of the SFR probability. This fraction is derived from the best-fitting relation from the Feedback in Realistic Environments (Hopkins et al. 2014, 2018; Muratov et al. 2015) simulation suite. The dependence of the ejection velocity on the galaxy circular velocity also derives from the Muratov et al. (2015) high-resolution simulations. The circular velocities are determined on the fly using a fast friends-of-friends finder (Davé et al. 2016).

Feedback from longer-lived stars (e.g., asymptotic giant branch (AGB) stars and Type Ia supernovae (SNe Ia)) are included as well (following Bruzual & Charlot 2003 stellar evolution tracks with a Chabrier 2003 initial mass function). We track the evolution of 11 elements: H, He, C, N, O, Ne, Mg, Si, S, Ca, and Fe. The yields for SNe Ia are taken from Iwamoto et al. (1999), assuming 1.4 M_⊙ of returned mass per supernova event. The AGB yields are drawn from the Oppenheimer & Davé (2008) lookup tables. Type II supernova yields derive from Nomoto et al. (2006) parameterizations, though they are reduced by 50% owing to studies that find that these yields return galaxies with metallicities roughly a factor of ∼2 too large at a fixed stellar mass when compared to the observed mass–metallicity relation (Davé et al. 2012).

The hydrodynamic simulations all use gizmo in the meshless finite mass (MFM) mode, with a cubic spline of 64 neighbors in the MFM hydrodynamics. This kernel is used to define the volume partition between gas elements, and the faces therefore over which the hydrodynamics is solved with the Riemann solver. Our final particle masses are as follows. The dark matter particle masses are M_DM = 1 × 10⁶h⁻¹ M_⊙, and the baryon particle masses are M_b = 1.9 × 10⁵h⁻¹ M_⊙. We employ adaptive gravitational softening for all particles throughout the simulation with minimum force softening lengths of 12, 3, and 3 pc for dark matter, gas, and star particles, respectively.

2.3. Dust Radiative Transfer Models

We generate the stellar SEDs and subsequent dust radiative transfer with powderday, a code package that wraps fsps (Conroy et al. 2009, 2010b; Conroy & Gunn 2010), hyperion (Robitaille 2011), and yt (Turk et al. 2011). This process is performed on all snapshots at redshifts z < 10 in post-processing.

We generate the stellar SEDs using fsps, in particular, their python hooks python-fsps.¹¹ The SEDs are calculated for each star particle as a simple stellar population based on age and metallicity, which is taken directly from the hydrodynamic simulations. We assume a Kroupa (2002) stellar IMF for the stellar SED generation and the Padova stellar isochrones (Marigo & Girardi 2007; Marigo et al. 2008). The sum of these stellar SEDs comprises the intrinsic stellar SED for any given model galaxy snapshot.

We then calculate the attenuation these SEDs experience by performing dust radiative transfer calculations. Functionally, the radiative transfer must occur on a grid. We therefore project the metal mass from the hydrodynamic simulations onto a 200 kpc octree grid centered on the central galaxy's center of mass. Each cell is recursively subdivided until a maximum of 32 particles are in a cell. In Appendix A, we show surface density images of our model galaxies next to the octree deposition. The dust mass is assumed to be a constant 0.4× the metal mass within a cell, as motivated by observational constraints from both local-epoch galaxies and those at high redshift (Dwek 1998; Vladilo 1998; Watson 2011).

The radiative transfer from stellar sources occurs in three dimensions in a Monte Carlo fashion using the dust radiative transfer code hyperion (Robitaille 2011). hyperion uses the Lucy (1999) algorithm for determining the converged equilibrium dust temperature and radiation field. Here radiation is emitted from stellar sources and then absorbed, scattered, and reemitted from dust in each cell in the octree grid. This process is iterated upon until the dust temperature has converged; convergence is determined when the energy absorbed by 99% of the cells has changed by less than 1% between iterations. We compute the emergent intensity for nine isotropic viewing angles around our model galaxies.

For the dust-grain properties themselves, we utilize the carbonaceous-silicate grain model from Draine & Li (2007) with a size distribution from Weingartner & Draine (2001). This model uses the Draine (2003) renormalization relative to H, and we assume ${R}_{{\rm{V}}}\equiv {A}_{{\rm{V}}}/E(B-V)=3.15$ . This curve is shown as a reference in Figure 1.

3. Insight from Simplified Toy Models

We begin the analysis for this paper by first conducting a series of simplified stellar population synthesis numerical experiments. It is from these that we will build the physical intuition necessary to understand the results of significantly more complex systems (i.e., the cosmological zoom galaxy formation simulations).

The relevant parameters in the dust attenuation curve for SED modeling are its normalization and shape. The overall normalization is trivially a function of the total dust column density seen by stars. As a result, we do not explore this further and concentrate hereafter on the shape of the attenuation curve.

When considering the shape (i.e., slope) of the attenuation curve, we assert that three principal physical drivers are at play: (i) the shape of the underlying extinction curve (which is itself dictated by the dust-grain size distribution), (ii) the star-to-dust geometry, and (iii) the average stellar age. We demonstrate the first (grain size variations) effect in Figure 2, the second effect (star-to-dust geometry) in Figure 3, and the final effect (stellar ages) in Figures 4 and 5.

**Figure 2.** Model SEDs and dust attenuation curves resulting from fsps stellar population synthesis calculations in which we vary the underlying dust extinction properties. In Figures 2–5, we intend to demonstrate that there are three fundamental drivers of variations in normalized attenuation laws: variations in the dust extinction properties, star-to-dust geometry, and stellar age effects. Here we concentrate on extinction properties. The stellar populations shown here are modeled as having a constant star formation history (with SFR = 1 M_⊙ yr⁻¹) and a stellar age of 13 Gyr. The left panel shows the UV–optical SED of a stellar population hidden behind a dust screen with a Kriek & Conroy (2013) extinction curve and a varying power-law index for that curve (δ). The right panel shows the modeled attenuation curves with these varying extinction curve indices. As is straightforward, variations in the dust extinction properties drive commensurate variations in the observed attenuation curves.
Download figure:
Standard image High-resolution image

**Figure 3.** Model SEDs and dust attenuation curves resulting from fsps stellar population synthesis calculations in which we vary the underlying dust extinction properties. In Figures 2–5, we intend to demonstrate that there are three fundamental drivers of variations in normalized attenuation laws: variations in the dust extinction properties, star-to-dust geometry, and stellar age effects. Here we concentrate on the star-to-dust geometry and model the population similarly as in Figure 2. The left panel shows the model SEDs with a varied fraction of stars that are not attenuated by this dust (i.e., a fraction of zero corresponds to all stars residing behind the dust screen, and a fraction of 1 corresponds to no obscuration). The right panel shows the resultant attenuation curves. As the fraction of unattenuated starlight increases in galaxies, the attenuation curves become shallower. This effect will map to increasingly complex star-to-dust geometries in high-redshift galaxies, with reduced optical depth sight lines to stars.
Download figure:
Standard image High-resolution image

**Figure 4.** Results from a population synthesis experiment aimed at simulating a galaxy dominated by old stars, with a minor amount of residual star formation. In this scenario, we see that the UV is dominated by young stars (that live in attenuating birth clouds), while the older stellar population dominates the optical. When comparing to Figure 5, we see that this scenario drives very steep attenuation curves. In detail, we create a composite stellar population in which 99% of the stellar mass in this stellar population is made up of a 13 Gyr population, and the remaining 1% is made up of a 0.01 Gyr population. The young stellar population sees a Charlot & Fall (2000) dust screen, while the older population does not.
Download figure:
Standard image High-resolution image

**Figure 5.** Results of numerical stellar population synthesis experiments in which we consider the impact of varying star formation histories and stellar ages on the derived dust attenuation curve. In short: galaxies whose mass is dominated by older stellar populations, though have residual star formation have relatively steep normalized attenuation curves. This is because the optical light, dominated by the older populations, sees little dust, while the UV (coming from new stars) is attenuated by dust associated with birth clouds. We demonstrate this by creating stellar populations with a range of star formation histories: exponentially decreasing (top row), constant (middle row), and rising (bottom row). These star formation histories are shown explicitly in the last column. The stars are placed behind a dust screen, with an escape time of 10 Myr (meaning any star older than this age will not reside in birth clouds anymore). The left column shows the attenuation curves normalized at 3000 Å, and the second column shows the absolute attenuation curve (i.e., not normalized at 3000 Å). The third column shows the UV–NIR SEDs for both the youngest and oldest stellar age in the star formation history. Owing to a combination of the mixing of young and old stars with the ISM and the relative UV and optical flux emitted by young and old stars, the stellar age of a galaxy can serve as a proxy for the underlying effects that drive variations in attenuation law slopes. See text for details.
Download figure:
Standard image High-resolution image

In Figures 2 and 3, we create (as an extreme example) a stellar population with fsps (Conroy et al. 2009, 2010b; Conroy & Gunn 2010) with age 13 Gyr formed with a constant 1 M_⊙ yr⁻¹ star formation history. Stars reside behind a dust screen whose opacity varies as a function of stellar age. Following Charlot & Fall (2000), we take

$\begin{eqnarray}\tau =\left\{\begin{array}{ll}{\tau }_{1}{(\lambda /5500\mathring{\rm A} )}^{{\delta }_{\mathrm{CF}}} & t\leqslant {10}^{7}\,\mathrm{yr}\\ {\tau }_{2}{(\lambda /5500\mathring{\rm A} )}^{{\delta }_{\mathrm{CF}}} & t\gt {10}^{7}\,\mathrm{yr},\end{array}\right.\end{eqnarray} \tag{ 1 }$

such that the dust screen in front of the stellar population in this experiment has a variable age-dependent normalization. Here we set the normalization for young stars (t_age < 10⁷ yr) as τ₁ = 1, while older stellar populations see a normalization τ₂ = 0.3 (see Conroy et al. 2010b). Following Charlot & Fall (2000), we set δ_CF = −0.7. The strong attenuation for young stars is intended to emulate the attenuation of starlight from newly formed stars by their birth clouds, while the more modest attenuation for older populations serves as a proxy for diffuse cirrus dust in galaxies. The extinction curve follows the Kriek & Conroy (2013) derived curve.¹²

In Figure 2, we vary the slope of the imposed Kriek & Conroy (2013) extinction curve as a proxy for variations in the underlying dust-grain properties and hold the dust covering fraction fixed at unity (meaning all stars see the dust screen). It is evident that variations in the intrinsic extinction curve propagate to a wide range of observed attenuation laws. At the same time, absent a model for the physical evolution of dust-grain properties in cosmologically evolving galaxies (e.g., McKinnon et al. 2016; Hou et al. 2017; Popping et al. 2017b), we are forced to assume a fixed underlying extinction curve and do not consider variations that originate in the physics of dust-grain property variations further. We instead focus on the impact of the star-to-dust geometry and stellar ages on the shape of the normalized attenuation curve.

In Figure 3, we show a similar model as in Figure 2, though we fix the slope of the Kriek & Conroy (2013) extinction curve (δ = −0.7) and instead vary the fraction of stars that do not see any dust. This model is intended to serve as a proxy for increasingly complex star-to-dust geometries in galaxies. A simplified geometry where all stars are enshrouded by dust is represented by a fraction of 0.0, while a model in which the stellar light and dust clouds are wholly decoupled is represented by a fraction 1.0. As the geometry becomes increasingly complex and more of the starlight is decoupled from dust, the effective attenuation curve becomes flatter (grayer), a well-known result dating back to seminal works by Witt & Gordon (1996).

Note that a "complex geometry" can encapsulate a broad range of effects. For example, a system that has a significant young stellar population that dominates the UV and is enshrouded in dust, but also a significant unobscured old star population that dominates the optical, will have a maximally steep attenuation curve. On the other hand, a galaxy with a significant unobscured young star population (e.g., Casey et al. 2014b; Geach et al. 2016; Narayanan et al. 2018) but obscured old star population would have an extremely shallow (gray) attenuation curve. Intermediate cases naturally fall between these two limits.

Because of this, alongside the typical geometry, the median stellar age of a galaxy can impact observed dust attenuation curves as well. In short, galaxies with a young median stellar age have both their UV and optical light dominated by young stars. Hence, the shape of the attenuation curve in this situation is dictated primarily by the fraction of obscured young stars. Galaxies with an older median stellar age, however, represent a different situation. Here the optical light is typically dominated by older stars, which tend to be free from obscuring dust, while the UV light is still dominated by young stars. In this case, the UV regime in the attenuation curve is dictated by the fraction of obscured young stars, while the optical is determined by the fraction of obscured old stars.

To see this explicitly, we develop a simplified population synthesis experiment in Figure 4. Here we have constructed a composite stellar population (red line) comprised of a 0.01 Gyr population (green line) and a 13 Gyr population (blue line). The 0.01 Gyr population constitutes 1% of the total mass, and the 13 Gyr population constitutes the remaining 99%. The intention behind this contrived experiment is to mimic a galaxy comprised principally (to the 99% level) of old stars, with only residual (1% of the total mass) current star formation. As usual in our population synthesis numerical experiments, the young stars are hidden behind a Charlot & Fall (2000) screen of dust with the same parameters used thus far. The old stars do not suffer any dust attenuation.

Evident in Figure 4 is that the new stars (that suffer attenuation from a dust screen) dominate the SED at far-UV (FUV) wavelengths, while the older stellar population dominates the optical. Because the older stars do not see a dust screen, this will result in increased transparency at optical wavelengths but depressed emission at FUV wavelengths (thanks to the young stars hidden behind the Charlot & Fall 2000 screen). This scenario can therefore play an important role in determining the shape of normalized attenuation curves.

To demonstrate the impact of stellar ages on the normalized attenuation curves even further, in Figure 5, we increase the complexity of the numerical experiments developed in Figure 4 and model three different star formation histories in population synthesis models. The top row shows an exponentially declining star formation history, the middle row a constant star formation history, and the bottom row a rising star formation history; the star formation history for each is quantified in the last (third) column. As in all simplified population synthesis experiments developed thus far in this paper, in Figure 5, we utilize the Charlot & Fall (2000)–type dust model, where young stars reside behind a dust screen but evolved stars do not.

In the first and second columns, we show the normalized and absolute attenuation curves for a range of stellar ages. The third column shows the star formation history. Let us first consider the case of the exponentially declining star formation history (row 1) as an instructive example. When examining the first column of the first row (normalized attenuation curves), we see that the models with the oldest ages have the steepest normalized curves. The reason for this becomes evident in the second column (absolute attenuation curve). As the star formation history evolves in the top row of Figure 5, the obscuration of the sources that dominate the FUV (x ≳ 6) photons does not change (thanks to the assumption of Charlot & Fall 2000 birth clouds), but that of the optical photons does. Galaxies with more evolved stellar populations still have their UV fluxes dominated by young stars, though their optical fluxes derive instead from older stellar populations. What this means is that if young stars are generally more obscured than older stars (in the stellar population synthesis models explored here, this is manifestly enforced), then galaxies with young median stellar ages will see significant obscuration in both the UV and the optical. Galaxies that are more evolved will still see significant obscuration in the UV (as the young stars are still enshrouded in dust), though the older stars that dominate the optical have lower optical depths (second column of Figure 5) and therefore steeper normalized attenuation curves.

The same effects are visible in the second and third rows of Figure 5 (constant and rising star formation history, respectively). That said, the dispersion in attenuation curve slopes is more modest compared to the exponentially declining history owing to the mixed contribution to optical light by both young and older stellar populations in the most evolved (t_age = 13 Gyr) stellar age bin.

4. Dust Attenuation Curves in Cosmological Galaxy Formation Simulations

With the insight we have built from our simplified stellar population synthesis models in Section 3, we now turn to the modeled curves in our cosmological zoom galaxy formation simulations. As a reminder: going forward, we will hold the underlying dust extinction properties fixed with a Weingartner & Draine (2001) size distribution. Beyond this, because the inclusion of subresolution birth clouds depends on tunable free parameters that can impact the amount of attenuated light (e.g., Narayanan et al. 2009, 2010; Hayward et al. 2013), we will abandon the usage of any subresolution "birth cloud" model in the zoom simulations.¹³ All attenuation seen by stars will be on larger scales (≳10 pc) from the geometry of the galaxy itself.

4.1. Diversity in Dust Attenuation Curves

To set the stage, in Figure 6, we show the diversity of attenuation curves for all snapshots of all model zoom galaxy formation simulations examined in this paper. The attenuation curves are normalized by their 3000 Å optical depth (which, hereafter, we refer to as τ₃₀₀₀). We additionally show a number of observationally constrained literature attenuation curves for reference. It is clear that a wide range of attenuation curves emerge from these simulations, with curves both steeper and grayer (shallower) than the standard literature assumptions. The goal of the remainder of this section is to unpack these curves further and examine in detail the physical drivers behind this diversity in attenuation curves.

4.2. Fitting Parameterizations

In our analysis of the physical drivers of the shapes of attenuation curves, we begin by describing the formulae that we employ to fit our model attenuation curves. It is useful to go through this exercise at this point, as this will allow us to introduce variables that we can use to characterize various properties of the modeled attenuation curve (e.g., their 2175 Å UV bump strengths or the normalized slope of the attenuation curve).

Broadly, we follow a modified version of the parameterizations described in Conroy et al. (2010b), which themselves are a modified version of the Cardelli et al. (1989) fits to the average Milky Way extinction curve. For the 2175 Å UV bump, we utilize a Drude profile, and in the NUV, we follow Noll et al. (2007).

In more detail, we first define x as the inverse wavelength with units of 1/μm. We fit the near-infrared (0.3 ≤ x < 2.5) with

$\begin{eqnarray}&&a{(x)}_{\mathrm{IR}}={a}_{1}\times {x}^{{\gamma }_{\mathrm{IR}}},\end{eqnarray} \tag{ 2 }$

$\begin{eqnarray}&&b{(x)}_{\mathrm{IR}}={b}_{1}\times {x}^{{\gamma }_{\mathrm{IR}}},\end{eqnarray} \tag{ 3 }$

$\begin{eqnarray}&&{f}_{\mathrm{IR}}=(a{(x)}_{\mathrm{IR}}+b{(x)}_{\mathrm{IR}})/{R}_{{\rm{V}}}.\end{eqnarray} \tag{ 4 }$

In the optical (2.9 < x ≤ 2.5), we define

$\begin{eqnarray}&&y=x-1.82\end{eqnarray} \tag{ 5 }$

and then parameterize the optical via

$\begin{eqnarray}\begin{array}{rcl}a{(x)}_{\mathrm{opt}} & = & 1+{a}_{1}y-{a}_{2}{y}^{2}-{a}_{3}{y}^{3}+{a}_{4}{y}^{4}\\ & & +\,{a}_{5}{y}^{5}-{a}_{6}{y}^{6}+{a}_{7}{y}^{7},\end{array}\end{eqnarray} \tag{ 6 }$

$\begin{eqnarray}\begin{array}{rcl}b{(x)}_{\mathrm{opt}} & = & {b}_{1}y+{b}_{2}{y}^{2}+{b}_{3}{y}^{3}-{b}_{4}{y}^{4}\\ & & -\,{b}_{5}{y}^{5}+{b}_{6}{y}^{6}-{b}_{7}{y}^{7},\end{array}\end{eqnarray} \tag{ 7 }$

$\begin{eqnarray}&&{f}_{\mathrm{opt}}=(a{(x)}_{\mathrm{opt}}+b{(x)}_{\mathrm{opt}})/{R}_{{\rm{V}}}.\end{eqnarray} \tag{ 8 }$

The NUV and 2175 Å UV bump are modeled via a Calzetti et al. (2000) law with a Lorentzian-like Drude profile. The latter is given by

$\begin{eqnarray}&&{D}_{{\lambda }_{0},\gamma ,{E}_{\mathrm{bump}}}=\displaystyle \frac{{E}_{\mathrm{bump}}{\lambda }^{2}{\gamma }^{2}}{{({\lambda }^{2}-{\lambda }_{0}^{2})}^{2}+{\lambda }^{2}{\gamma }^{2}},\end{eqnarray} \tag{ 9 }$

where λ₀ = 2175 Å is the central wavelength of the bump, γ is the width, and E_bump is the amplitude (Fitzpatrick & Massa 2007; Noll et al. 2009). Formally, the fit for the NUV and bump take the form

$\begin{eqnarray}\begin{array}{rcl}{f}_{\mathrm{UV},\mathrm{bump}} & = & {A}_{\mathrm{NUV}\mathrm{norm}}\left[(k(\lambda )+{D}_{{\lambda }_{0},\gamma ,{E}_{\mathrm{bump}}})\displaystyle \frac{1-1.12{c}_{{\rm{R}}}}{{R}_{{\rm{V}}}}+1\right]\\ & & \times \,{\left(\displaystyle \frac{\lambda }{5500\mathring{\rm A} }\right)}^{{\delta }_{\mathrm{NUV}}}.\end{array}\end{eqnarray} \tag{ 10 }$

Here A_NUV norm is a normalization of the fit, and k(λ) is the Calzetti et al. (2000) law over the wavelengths of interest. The $(\lambda /5500{)}_{\mathrm{NUV}}^{\delta }$ term allows for an arbitrary tilt to the curve in this wavelength regime, and the $1\mbox{--}1.12{c}_{{\rm{R}}}/{R}_{{\rm{V}}}$ term compensates for how R_V will change in the original Calzetti law due to the introduction of a tilt (Noll et al. 2009).

Finally, we model the FUV in the range 5.9 ≤ x < 10 via

$\begin{eqnarray}&&{f}_{a}={a}_{1,\mathrm{FUV}}{(x-5.9)}^{2}-{a}_{2,\mathrm{FUV}}{(x-5.9)}^{3},\end{eqnarray} \tag{ 11 }$

$\begin{eqnarray}&&{f}_{b}={b}_{1,\mathrm{FUV}}{(x-5.9)}^{2}+{b}_{2,\mathrm{FUV}}{(x-5.9)}^{3},\end{eqnarray} \tag{ 12 }$

$\begin{eqnarray}&&a(x)=1.752-0.316\times x-\displaystyle \frac{0.014B}{{(x-4.67)}^{2}+0.341}+{f}_{a},\end{eqnarray} \tag{ 13 }$

$\begin{eqnarray}&&b(x)=-3.09+1.825\times x+\displaystyle \frac{1.206B}{{(x-4.62)}^{2}+0.263}+{f}_{b},\end{eqnarray} \tag{ 14 }$

$\begin{eqnarray}&&{f}_{\mathrm{FUV}}=\left(\displaystyle \frac{a(x)+b(x)}{{R}_{{\rm{V}}}}\right){A}_{\mathrm{FUV},\mathrm{norm}}.\end{eqnarray} \tag{ 15 }$

In other parameterizations of attenuation curves (e.g., Conroy et al. 2010b), B is a parameter describing the bump strength. Because we employ the Noll et al. (2009) model for the NUV and bump, here (in the FUV) B serves simply as another free parameter in the fit. In this wavelength regime, we linearly interpolate over the Lyα line to avoid complications during fitting.

We note that while a piecewise fit to attenuation curves akin to what is presented here has been used regularly in the literature (e.g., Cardelli et al. 1989; Conroy et al. 2010b), a number of studies in the literature make use of a modified Calzetti et al. (2000) curve over the entire x = 1–10/μm wavelength range,

$\begin{eqnarray}&&{A}_{\lambda }={\rm{\Gamma }}({k}_{\mathrm{cal}}^{{\prime} }+D(\lambda ))+{\left(\displaystyle \frac{\lambda }{{\lambda }_{{\rm{V}}}}\right)}^{{\delta }_{\mathrm{cal}}},\end{eqnarray} \tag{ 16 }$

where ${k}_{\mathrm{cal}}^{{\prime} }$ is the Calzetti et al. (2000) relation, δ_cal is the index for a power-law modification to this relation, $D(\lambda )$ is the normal Drude profile (to represent the 2175 Å bump), and Γ is a constant free parameter. Various forms of this relation exist in the literature, where the constant may be multiplied by all or only some of the terms in Equation (16). Here Γ is a constant that typically relates either R_V or A_V to the ratio of total to selective extinction for the Calzetti et al. (2000) relation (i.e., ${R}_{{\rm{V}},\mathrm{cal}}=4.05$ ), but the exact implementation of this also varies in the literature (e.g., Noll et al. 2009; Kriek & Conroy 2013; Salmon et al. 2016; Salim et al. 2018). Similarly, not all authors include the Drude representation of a bump. In order to compare against observations, we will find it useful at times in this paper to employ Equation (16) to fit over the entire wavelength regime of the attenuation curve. In this case, we will note the power-law indices used in this case as δ_cal instead of δ_NUV, which we shall reserve for our normal fitting procedure.

Beyond this, some authors fix the width (γ) of the Drude profile representing the bump strength. In this situation, we denote the normalization of the bump as ${E}_{\mathrm{bump},\mathrm{KC}}$ and the width of the bump as γ_const:

$\begin{eqnarray}&&{D}_{{\lambda }_{0},\gamma ,{E}_{\mathrm{bump}}}=\displaystyle \frac{{E}_{\mathrm{bump},\mathrm{KC}}{\lambda }^{2}{\gamma }_{\mathrm{const}}^{2}}{{({\lambda }^{2}-{\lambda }_{0}^{2})}^{2}+({\lambda }^{2}\times {({\gamma }_{\mathrm{const}})}^{2})}.\end{eqnarray} \tag{ 17 }$

For clarity, we collect the fitting variables from this section that we will employ to compare with observations in Table 2.

Table 2. Definitions of Fitting Variables that Will Be Used to Characterize the Shapes of Attenuation Curves in This Paper

Fitting Variable	Definition	Equation
δ_NUV	Power-law index in NUV wavelength range	Equation (10)
${\delta }_{\mathrm{cal}}$	Power-law index of modified Calzetti et al. (2000) curve	Equation (16)
${\int }_{\mathrm{NUV}}{D}_{{\lambda }_{0},\gamma ,{E}_{\mathrm{bump}}}$	Integral of Drude profile in NUV: bump strength	Equation (9)
E_bump	Normalization of UV bump, given a variable bump width γ	Equation (9)
${E}_{\mathrm{bump},\mathrm{KC}}$	Normalization of UV bump, assuming a constant bump width ${\gamma }_{\mathrm{const}}$	Equation (17)
γ_const	Value of constant bump width for ${E}_{\mathrm{bump},\mathrm{KC}}$	Equation (17)

Download table as: ASCII Typeset image

4.3. What Sets the Slope of Dust Attenuation Curves?

Building from our intuition in Section 3, we examine the physical drivers of the slopes of our model dust attenuation curves. Fundamentally, the principal driver of variations in the slope of attenuation curves is the geometry (assuming a fixed extinction curve). Recalling Section 3, the steepest attenuation curves arise from galaxies with a significant fraction of young stars (dominating the emitted UV flux) obscured by dust but old stars (dominating the optical flux) that are not obscured. A more mixed geometry (i.e., more old stars obscured by dust or more young stars decoupled from dust) will flatten the attenuation curve from this extreme limit. In Section 3, we showed this explicitly with simplified population synthesis models; we now explore these effects in bona fide cosmological galaxy formation simulations. We center this discussion around Figure 7, where we show various incarnations of the attenuation curves of an example model galaxy, model mz10.

In the top left panel of Figure 7, we show the normalized attenuation curves for model mz10, color-coded by their fraction of unobscured young stars. This quantity is determined by calculating the fraction of stellar mass in the form of young (t_age < 50 Myr) stars that have no dust within 250 pc. In other words, larger values of this fraction correspond with increasingly decoupled geometry between young stars and dust in galaxies. To minimize the complicating effects of older stellar populations (an effect we will return to shortly), we only plot the attenuation curves for galaxies with a median stellar age (by mass) <0.5 Gyr; this said, in light gray, we show the attenuation curve for all snapshots of this model (i.e., of all median stellar ages). While the dynamic range in the attenuation curve slopes is relatively small, it is clear from the top left panel of Figure 7 that a larger fraction of unobscured young stars results in a flatter dust attenuation curve.

Similarly, the geometry between old stars and dust also plays a role in the tilt of the observed attenuation curves. To demonstrate this, in the top right panel of Figure 7, we show the same attenuation curves as shown in the top left panel, though in this case we color-code the attenuation curves by their fraction of unobscured old stars. As in the top left panel, we define "unobscured" as having no dust within 250 pc and old stars as those with t_age > 50 Myr (where t_age is the median stellar age). To minimize contamination by young stellar populations, we only include galaxies with a median stellar age >0.5 Gyr, but we show the curves for all snapshots in light gray. As in the simplified stellar population synthesis models presented in Figure 5, galaxies that contain a significant amount of old stars that are unobscured by dust have steep attenuation curves, while those that have larger obscuration of old stars have flatter (grayer) attenuation curves.

Why are the curves in the top left panel of Figure 7 that are associated with young galaxies uniformly shallower (grayer) than the older stellar age curves in the top right panel of Figure 7? As shown in Figures 4 and 5, in galaxies where the median stellar age is relatively young, the UV and optical light are both dominated by newly formed stars. If the light from these stars is attenuated by dust, then both the optical and UV emission from the galaxy will be attenuated. As a result, these curves will not be as steep as a situation where young stars (that are obscured) dominate the UV emission but old stars are relatively unobscured (a situation described by the top right panel of Figure 5).

What this results in is a situation where galaxies with older stellar ages, on average, have steeper attenuation curves. In the bottom left panel of Figure 7, we show this by plotting the same attenuation curves as in the other two panels but this time color-coded by the median stellar age. As the galaxy age increases, so does the slope. This is due to an increasingly large fraction of the optical emission coming from older stars as the galaxy ages. Older star particles, on average, are less likely to be associated with dust than young stars. This same effect was noted by Charlot & Fall (2000), who noted steeper attenuation curves with increasing galaxy age.

4.4. Variations in the 2175 Å Bump

As discussed in Section 1, while one of the strongest features in the attenuation curve of the Milky Way is the bump feature observed at ∼2175 Å, a number of observations show dramatic variations in bump strength in different galaxies (e.g., Gordon et al. 1999, 2000; Calzetti et al. 2000; Motta et al. 2002; Burgarella et al. 2005; York et al. 2006; Noll et al. 2007; Stratta et al. 2007; Elíasdóttir et al. 2009; Conroy 2010; Buat et al. 2011; Wild et al. 2011; Kriek & Conroy 2013; Scoville et al. 2015; Battisti et al. 2016; Salim et al. 2018). In this section, we utilize the model that we have developed thus far to understand the origin of variations in the UV bump strength in galaxy attenuation curves.

We define the strength of the UV 2175 Å bump as the integral over the best-fit Drude profile in the NUV (see Equation (9)):

$\begin{eqnarray}&&{\rm{b}}{\rm{u}}{\rm{m}}{\rm{p}}\,{\rm{s}}{\rm{t}}{\rm{r}}{\rm{e}}{\rm{n}}{\rm{g}}{\rm{t}}{\rm{h}}\equiv {\int }_{{\rm{N}}{\rm{U}}{\rm{V}}}{D}_{{\lambda }_{0},\gamma ,{E}_{{\rm{b}}{\rm{u}}{\rm{m}}{\rm{p}}}}.\end{eqnarray} \tag{ 18 }$

Note that literature definitions of the bump strength vary, with some definitions characterizing the strength as E_bump (i.e., one numerator term in Equation (9)), as well as ${E}_{\mathrm{bump},\mathrm{KC}}$ , where the width of the bump is held fixed; see Table 2.

Because the UV bump is an absorption feature, for a fixed underlying extinction curve, reduced bump strength in the observed attenuation curve signifies extra radiation filling in the bump. In principle, there are two possible sources of this extra radiation: light being scattered into the line of sight and unobscured sources of UV radiation (i.e., unobscured young stars).

In Figure 8, we show how the bump strength varies in our models with a number of relevant quantities. The top left panel of Figure 8 shows the relationship between the 2175 Å UV bump strength and the fraction of scattered light contributing to the total 2175 Å flux. There is a relatively weak correlation between the two: light scattered into the line of sight contributes modestly to reduced bump strengths in some attenuation curves, but it clearly does not dominate.

In the top right and bottom left panels of Figure 8, we investigate the role of unobscured stars in contributing to reduced bump strengths. The top right panel shows the bump strengths as a function of the fraction of unobscured old (t_age > 50 Myr) stars, while the bottom left panel shows the bump strengths as a function of the fraction of unobscured young (t_age < 50 Myr) stars. As is clear, the old stars, which put out relatively little flux in the NUV, do not contribute significantly to filling in the 2175 Å absorption bump, even at relatively large unobscured fractions (and, in fact, trend in the opposite direction). In contrast, however, there is a clear relationship between the fraction of unobscured young stars and the UV bump strength. In galaxies with increasingly complex young geometries, as more young stars find low-τ UV sight lines, this radiation fills in the UV bump and reduces its strength. In short, galaxies with very complex young geometries have reduced 2175 Å UV bumps.

Drawing on what we learned in Section 4.3, the fraction of unobscured young stars in our simulations also correlates with the slope of the dust attenuation curve. It follows transitively, then, that the bump strength should vary inversely with the attenuation curve slope. We quantify this relationship in the bottom right panel of Figure 8, where we show the attenuation curve bump strength versus NUV slope (δ_NUV; see Equation (10)). The grayest attenuation curves represent galaxies with the most complex young geometry and therefore have the weakest bump strengths. While we will return to this issue further in Section 5, it is worth briefly noting that Kriek & Conroy (2013) demonstrated that for z ∼ 2 galaxies in the NEWFIRM survey, the dust attenuation slope varies inversely with the measured bump strength. This may provide some tentative evidence for our interpretation.

4.5. Impact of Inclination Angle

A natural consequence of the impact of the geometry on attenuation curves is that, for a relatively disklike system, the inclination angle can play a role in driving the observed tilt. This said, especially at high redshift, the geometry of massive galaxies can often depart from well-ordered disks, thus minimizing the role of inclination angle on the observed attenuation curve slope.

As an estimate of the impact of inclination angle, in Figure 9, we show the attenuation curves for 25 randomly selected galaxies from our entire model sample. In each panel, we show the curves from nine isotropic viewing angles. Maximum dispersions in the FUV are of order ∼2–3, though the typical is much less than that. This is due primarily to the relatively small numbers of well-ordered disks in our simulation sample (which are predominantly the z ∼ 0 snapshots of models mz352 and mz401).

5. Discussion

5.1. Does a Single Attenuation Prescription Apply?

Thus far, we have explored the physical underpinnings of variations in dust attenuation curves in galaxies. We now ask the slightly more practical question: what range of attenuation curves can one expect for galaxies at a given redshift?

To answer this, we employ the (25/h Mpc)³ mufasa cosmological hydrodynamic galaxy formation simulation (Davé et al. 2016). This simulation has identical physics to our zoom models, save for the inclusion of the standard mufasa heuristic quenching model. The simulation is run with 512³ particles but in a volume half the size of the one our zooms are selected from, resulting in an effective mass resolution a factor of 8 worse (i.e., larger particle masses) than our zoom models. Achieving the full resolution of our zoom simulations in a cosmological volume requires a computational effort outside the scope of the current investigation, though in Appendix B, we quantify the level of uncertainty introduced by these lower-resolution simulations.

In Figures 10 and 11, we show heat maps of the attenuation curves for all galaxies¹⁴ in the 25 Mpc³ volume at redshifts z = 0, 2, 4, and 6. We additionally show the median attenuation curve (computed by calculating the median τ/τ₃₀₀₀ at every wavelength x = 1/λ) at each redshift via the solid pink line and the 1σ standard deviation in the dashed pink lines. In order to best compare with observations, which typically only select star-forming galaxies, we restrict our analysis here to galaxies with SFR ≥ 1 M_⊙ yr⁻¹. Figure 10 shows the heat map and median attenuation curves, and Figure 11 shows the median curves in comparison to literature references.

**Figure 10.** Heat map of attenuation curves at redshifts z = 0, 2, 4, and 6 derived from a mufasa (25/h Mpc)³ cosmological simulation. The physics of this simulation is identical to those in our zooms. The solid pink line is the median attenuation curve in all panels, while the dashed pink lines show the 1σ standard deviation. The median attenuation curves become grayer with redshift as geometries become more complex and median stellar ages become more uniform. There is significant dispersion at all redshifts, though the median is typically bounded by the family of curves that describe local galaxies (see Figure 11). We publish the best fit to our median curves at all integer redshifts between z = 0 and 6 in a public repository; see text for details.
Download figure:
Standard image High-resolution image

**Figure 11.** Median attenuation curves as a function of redshift as taken from Figure 10, with standard literature curves presented for comparison. The median curves are comparable to many locally calibrated curves through z = 4, though they become significantly grayer toward z ∼ 6. See text for details.
Download figure:
Standard image High-resolution image

At each redshift, there is significant dispersion about the median, as evidenced by the heat map. In the normalized attenuation curves, a wide range of curve tilts and bump strengths are exhibited. The significant dispersion seen in Figure 10 is expected from observational constraints that show a diverse range of slopes in attenuation curves (e.g., Wild et al. 2011; Battisti et al. 2017; Salim et al. 2018). This dispersion in curves decreases with increasing redshift. While galaxy geometry remains complex at these epochs, the metal content (and hence, manifestly in these simulations, dust content) decreases, thus reducing the impact of geometry on the attenuation curve shape. Beyond this, there is a much narrower distribution in median stellar ages with increasing redshift (see Section 4.3).

Despite the strong dispersion seen at nearly all redshifts, the medians in the broad distribution of modeled attenuation curves are always bounded by the standard literature attenuation curves. This is demonstrated in Figure 11, where we compare against standard literature curves. For convenience, we have published the best-fitting median curves in Figure 10 at integer redshifts from z = 0 to 6 on a publicly accessible website.¹⁵ It is important to note that while the median values are appropriate for ensemble averages, these median values (or any assumption of a locally calibrated curve, for that matter) may grossly mischaracterize the underlying attenuation on a case-by-case basis.

Finally, it is worth noting that the most common attenuation curves in any of our modeled redshift bins all have prominent bump features. It is indeed possible to generate attenuation curves with relatively small bump contributions: this is demonstrated explicitly in Figures 6 and 8. These curves with minimal bumps owe their origin solely to geometry and radiative transfer effects (i.e., no modification of the underlying dust properties is necessary). This said, this is not typical in our models: the median curve within any redshift bin in Figure 10 displays a prominent bump.

5.2. Comparison to Other Models and Constraints

Our work builds on a deep theoretical literature in this area. In this section, we aim to place the results from our work in this context. We painted the landscape for theoretical work in this field in Section 1. In this section, we aggregate some key results from these papers and compare them to our own investigation.

5.2.1. On the Role of Geometry

In Figures 3 and 7, we demonstrate via both population synthesis experiments and direct cosmological simulation that the role of geometry is paramount in driving the tilt of normalized attenuation curves. For young galaxies, an increasing fraction of unobscured young stars flattens normalized attenuation curves, while for galaxies dominated by older stellar populations, an increasing fraction of unobscured old stars steepens normalized curves.

In general, at least the former point (and, more broadly, the idea that geometry plays an important role in setting the shape of attenuation curves) is already well appreciated in the theoretical literature. Indeed, a broad range of modeling techniques arrive at the same conclusion. For example, Witt & Gordon (1996) examined radiative transfer models in a two-phase clumpy medium and found that as a general rule, attenuation curves became grayer (flatter) as the obscuration inhomogeneity increased. Ferrara et al. (1999) and Bianchi et al. (2000) developed analytic models of the geometry in galaxies and found that the inferred attenuation curve depends on the magnitude of clumping in the ISM, as well as the fraction of starlight emitted within these clumps. Seon & Draine (2016) expanded upon these earlier works by studying the radiative transfer in a turbulent medium and arrived at a similar conclusion, while Natale et al. (2015) utilized a coupling of 3D dust radiative transfer calculations (as in this paper) with idealized hydrodynamic models of disks in evolution to demonstrate this principle. Fischera & Dopita (2011) additionally employed nonhomogeneous ISM models to explain gray attenuation curves. Finally, Trayford et al. (2017) applied a similar analysis to EAGLE cosmological galaxy formation simulations and highlighted the role of geometry in driving attenuation curve variations at z ∼ 0.1.

It is important to note that the concept of increased complexity in the geometry driving grayer attenuation curves only applies to galaxies whose light is dominated by young stellar populations (see Figure 7). In our model, galaxies whose luminosity is dominated by older stellar populations have normalized attenuation curves that become steeper as more old stars are decoupled from dust.

5.2.2. Bump Strengths

When assuming an attenuation curve for the purposes of SED fitting in low-metallicity galaxies (especially at high redshift), a common assumption is an SMC-like attenuation curve. This is motivated by the assumption that the bumpless curve from the SMC may owe its origin to different grain compositions (presumably driven by the galaxy's low metallicity) than in the Milky Way. Similar logic oftentimes motivates the usage of a bumpless Calzetti et al. (2000) curve for observations of heavily star-forming galaxies.

Throughout this work, we have demonstrated that even without changing the underlying grain properties (i.e., our extinction curve in every model is identical), we are able to generate extinction curves that have a diverse range in bump strengths (see, e.g., Figure 10 for one such example). In other words, curves with very small bump strengths in our models result simply from the geometry of the system. This conclusion is not necessarily shared among other theoretical models.

For example, Witt & Gordon (2000) expanded on the radiative transfer models of Witt & Gordon (1996) and found that they required a change to the intrinsic dust curve in order to make a bumpless curve. Specifically, only by employing an extinction curve without a 2175 Å bump were they able to reproduce observed bump-free curves. Others (e.g., Fischera & Dopita 2011) have suggested models in which the UV bump carriers are destroyed at threshold column densities in order to reproduce the Calzetti et al. (2000) curve.

Hou et al. (2017) implemented a full dust formation and destruction model into cosmological hydrodynamic simulations and tracked the grain size distribution and extinction curve. While they did not model the radiative transfer associated with these simulations, they found in the extinction curve that 2175 Å bump strengths are naturally correlated with dust growth dominated by accretion of metals onto dust grains. While this study does not explore the geometry effects that drive attenuation curves, it is certainly conceivable that intrinsic dust properties can impact the strength of the UV bump.

At the same time, other models have suggested, like this work, that it may be possible to achieve bump-free attenuation curves without modifying the underlying extinction curve. For example, Granato et al. (2000) developed SAMs to demonstrate that the bump-free Calzetti et al. (2000) curve may be a result of a "birth cloud" model, wherein the radiation from young stars is heavily attenuated by their birth clouds, but older stars are only attenuated by diffuse ISM. In this picture, the observed UV is dominated by an older stellar population that sees relatively little diffuse dust; hence, the emission from these objects can fill in the 2175 Å absorption feature. This SAM was expanded upon by Panuzzo et al. (2007), who found similar results.

While the Granato et al. (2000) and Panuzzo et al. (2007) studies share with our model the idea that an attenuation curve with very small bump strengths is achievable without changing the underlying extinction properties of grains, the models share a key difference. Granato et al. (2000) and Panuzzo et al. (2007) found that, as young stellar populations become more obscured, 2175 Å bump sizes decrease in attenuation curves. In contrast, in Figure 8, we demonstrate that the opposite is true for our models: as a larger fraction of young stars becomes unobscured, observed bump strengths decrease, owing to the UV radiation from these young stars filling in the 2175 Å absorption feature. An important corollary to our model is that the same displacement between O and B stars and sites of dust obscuration results in grayer (flatter) attenuation curves, leading to a natural relationship where shallower attenuation curves have smaller bump features in our model (bottom right panel of Figure 8).

Seon & Draine (2016) similarly derived a model for dust attenuation curves that exhibits a relationship between the bump strength and slope of the attenuation curve in the same direction as both our models and the Kriek & Conroy (2013) observations. This too can be attributed to geometry effects, wherein as the clumping or size of the source distribution increase, attenuation curves become grayer and exhibit smaller 2175 Å bumps. We note, however, that a Calzetti et al. (2000) curve is only attainable using the Weingartner & Draine (2001) Milky Way dust model when the intrinsic bump is removed or suppressed in the Seon & Draine (2016) model. As a result, their model may be viewed as intermediate between the results of Granato et al. (2000), Panuzzo et al. (2007), and us.

5.3. Comparison with Observations

We now turn to a comparison of our models with relatively recent observational results in this area. We remind the reader of the discussion surrounding Equations (16) and (17) and Table 2. Specifically, while we find a piecewise fit to provide the best fits to our model attenuation curves as outlined in Section 4.2, many observational studies employ a modified Calzetti et al. (2000) relation, where both a Drude-like profile for the UV bump and a power-law modification are employed (e.g., Noll et al. 2009). As we clarify in Table 2, we distinguish the power-law index derived from this method of fitting (δ_cal) from the power law we typically employ just in the NUV bands (δ_NUV).

We compare to the observational results of Kriek & Conroy (2013), Salmon et al. (2016), and Salim et al. (2018) in Figure 12. Salmon et al. (2016) and Salim et al. (2018) derive attenuation laws for redshift z ∼ 2 and z ∼ 0 galaxies, respectively, via SED fitting techniques. Evident from both of these studies is a relationship between the V-band optical depth and the slope of the attenuation curve, δ_cal. This is similar to the power-law relationship modeled by Chevallard et al. (2013) and Leja et al. (2017) between the optical depth of diffuse dust and the power-law slope of the attenuation curve. In the left panel of Figure 12, we show a comparison between our models and this observed trend. We include every galaxy in our sample of zooms, though we note that the Salmon et al. (2016) and Salim et al. (2018) studies of course both employ individual selection techniques within particular redshift ranges. Given this, as well as the relative uncertainties involved in deriving attenuation curves from SED fitting, the trend of increasing δ_cal with V-band optical depth in the model galaxies is encouraging.

Kriek & Conroy (2013) employed SED fitting techniques to observations of z ∼ 2 galaxies to derive a relationship between the UV bump strength and slope power-law index, δ_cal. In Figure 8, we demonstrated a similar relationship, though we characterized this in terms of the integrated Drude profile, $\int {D}_{{\lambda }_{0},\gamma ,{E}_{\mathrm{bump}}}$ , and the NUV slope, δ_NUV. Kriek & Conroy (2013) fixed the width of the bump to γ_const = 350 Å and therefore characterized the strength of the bump by its normalization, ${E}_{\mathrm{bump},\mathrm{KC}}$ (see Table 2). In order to best compare to this study, we have reperformed our fits by fixing our bump widths similarly, and we report in the right side of Figure 12 our modeled relationship between E_bump and δ_cal. We show our model points in blue and compare these to the Kriek & Conroy (2013) data in orange. The best-fit lines for both are shown. Our model galaxies show a similar trend as the Kriek & Conroy (2013) observations in that steeper curves tend to have more prominent bump strengths. In our model, this is due primarily to unobscured young stars reducing bump strengths. Our modeled best-fit relation is

$\begin{eqnarray}&&{E}_{\mathrm{bump},\mathrm{KC}}=-0.46\times {\delta }_{\mathrm{cal}}+0.69.\end{eqnarray} \tag{ 19 }$

Our model galaxies exhibit a much shallower gradient than what is observed. This may be due to a number of causes. First, our modeled bump widths, γ (see Equation (10)), span a broad range of values, with some widths exceeding twice the assumed γ_const = 350 Å employed for the construction of Figure 12. Forcing a bump strength, therefore, may result in poorly performing fits, and therefore E_bump values that do not reflect the true dynamic range of bump strengths. As an example, examination of the bottom right panel of Figure 8 demonstrates that when considering the integral of the Drude profile of the bump strength, we see a dynamic range spanning an order of magnitude in bump strength, unlike the factor of ∼2–4 when characterizing the bump strength by E_bump alone. Beyond this, a true apples-to-apples comparison between our model points and the observed data would involve our fitting our model SEDs (assuming a given star formation history and IMF) and recovering the inferred attenuation law accordingly. Exploration of the differences resulting in inferred attenuation curves from SED modeling from those directly modeled will be presented in future work. We note that Seon & Draine (2016) derived relationships between ${E}_{\mathrm{bump},\mathrm{KC}}$ and δ_cal that were typically steeper than the Kriek & Conroy (2013) relation. More observational and theoretical work in this area is warranted.

5.4. On the Role of the Underlying Galaxy Formation Model

In recent years, there have been a number of major cosmological galaxy formation simulation campaigns by various groups that couple hydrodynamic simulations of galaxy evolution with 3D dust radiative transfer (e.g., Snyder et al. 2015a, 2015b; Torrey et al. 2015; Camps et al. 2016, 2018; Snyder et al. 2017; Liang et al. 2018). These simulation efforts span a wide range of physics implemented into the galaxy formation simulations, resolution scales, and assumptions regarding subresolution obscuration. In principle, it is not obvious that two individual galaxy formation models starting with the same initial conditions would result in the same calculated attenuation curve. While wind models are generally tuned to reproduce bulk galaxy characteristics (e.g., stellar mass functions), the individual star-to-dust geometries may vary substantially from model to model. To truly characterize the theoretical expectation for attenuation curve variations in galaxy formation simulations, therefore, it would be ideal to employ galaxy formation simulations with varying subresolution physical models but otherwise identical initial conditions, physical resolution, and dust radiative transfer algorithms. While such an endeavor is outside the scope of this work, it is in principle possible using either powderday or skirt3D (Camps & Baes 2015).

6. Conclusions and Summary

We have developed a model for the origin of variations in dust attenuation curve shapes and bump strengths. We accomplish this by combining high-resolution cosmological galaxy formation simulations with 3D dust radiative transfer calculations. Critically, for these radiative transfer models, we hold the underlying extinction curve fixed and ask how geometry and radiative transfer effects impact the resultant attenuation curves. Our main results are as follows.

1.
Despite the usage of a constant extinction curve in our underlying radiative transfer calculations, we find dramatic variations in the derived attenuation laws. These variations depend primarily on complexities in the star-to-dust geometry (Figures 3, 6, and 7). The details are as follows.
- (a)
  Increasing fractions of unobscured young stars result in flatter (grayer) attenuation curves as the galaxy becomes more transparent to UV radiation (Figures 3 and 7).
- (b)
  Increasing fractions of unobscured evolved stars result in steeper attenuation curves (Figures 4, 5 and 7).
- (c)
  These results taken together drive a trend where galaxies with highly obscured sight lines toward young stars but a significant (unobscured) evolved stellar population will have the steepest normalized attenuation curves (Figure 7).
2.
The 2175 Å UV bump strengths vary dramatically, despite our usage of an extinction curve with a UV bump present. Unobscured O and B stars result in reduced bump strengths in our model, with scattered light only having a secondary effect on the feature (Figure 8).
3.
The combined effect of unobscured young stars both flattening attenuation curve slopes and reducing the bump strength results in a natural relationship wherein the slope of the attenuation curve is related to the bump strength: flatter attenuation curves tend to have smaller bump strengths (Figure 8).
4.
We apply these results to a (25/h Mpc)³ cosmological volume to derive the median curve and expected dispersion at integer redshifts from z = 0 to 6. While the median curve at a given redshift is typically bounded by standard literature curves, the dispersion is significant. The average dispersion decreases with increasing redshift, and the median curves become grayer. This is due to reduced dispersion in star-to-dust geometry, as well as a narrower distribution in median stellar ages with redshift (Figure 10). We publish these median curves on a public-facing website.

DN is grateful to Andrew Battisti, Daniela Calzetti, Ignacio Ferreras, Rob Kennicutt, Maciej Koprowski, Danielle Pierni, Karin Sandstrom, Samir Salim, Brett Salmon, James Trayford, Monica Tress, George Privon, and Adolf Witt for valuable conversations during this study. We additionally thank Mariska Kriek, Samir Salim, and Brett Salmon for providing data to us for our comparisons with observations. After posting to the arXiv, George Privon, Samir Salim, and Mike Shull individually alerted us to typographical errors, which we are grateful for. The simulations published here were run on the University of Florida HiPerGator supercomputing facility, and the authors acknowledge the University of Florida Research Computing for providing computational resources and support that have contributed to the research results reported in this publication. This study was funded in part by NSF AST-1715206 and HST AR-15043.0001.

Appendix A: Octree Deposition

In Figure 13, we show gas surface density images (top row) of an example model galaxy (mz401 at $z\approx 0$ ) and the octree deposition for each snapshot.

Appendix B: Resolution Tests

In order to test the level of uncertainty introduced in our utilization of lower-resolution simulations in Figure 10, we have rerun simulation mz0 with one lower level of refinement in the generation of the initial conditions. This results in a zoom with the same resolution as the large-box cosmological simulations used in Figure 10. We present the normalized attenuation curve for both this low-resolution run and our fiducial run in Figure 14 for a sample of randomly selected snapshots.

**Figure 14.** Resolution test comparing 16 randomly selected snapshots from a lower-resolution run of model mz0 to our fiducial resolution. By and large, the attenuation curves are similar, with peak differences less than a factor of ≲2.
Download figure:
Standard image High-resolution image

For the most part, the resolution-related effects on the shape of the attenuation curve are relatively minor, with a few outliers showing differences of order a factor of ∼2. These outliers are relatively uncommon, however, and we can use these as an upper limit as to the level of uncertainty in our analysis in Section 5.1.

A Theory for the Variation of Dust Attenuation Laws in Galaxies

Article metrics

Permissions

Author e-mails

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction

2. Simulation Methodology

2.1. Overview

2.2. Cosmological Zoom Galaxy Formation Simulations

2.3. Dust Radiative Transfer Models

3. Insight from Simplified Toy Models

4. Dust Attenuation Curves in Cosmological Galaxy Formation Simulations

4.1. Diversity in Dust Attenuation Curves

4.2. Fitting Parameterizations

4.3. What Sets the Slope of Dust Attenuation Curves?

4.4. Variations in the 2175 Å Bump

4.5. Impact of Inclination Angle

5. Discussion

5.1. Does a Single Attenuation Prescription Apply?

5.2. Comparison to Other Models and Constraints

5.2.1. On the Role of Geometry

5.2.2. Bump Strengths

5.3. Comparison with Observations

5.4. On the Role of the Underlying Galaxy Formation Model

6. Conclusions and Summary

Appendix A: Octree Deposition

Appendix B: Resolution Tests

Footnotes

A Theory for the Variation of Dust Attenuation Laws in Galaxies

Article metrics

Permissions

Share this article

Author e-mails

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction

2. Simulation Methodology

2.1. Overview

2.2. Cosmological Zoom Galaxy Formation Simulations

2.3. Dust Radiative Transfer Models

3. Insight from Simplified Toy Models

4. Dust Attenuation Curves in Cosmological Galaxy Formation Simulations

4.1. Diversity in Dust Attenuation Curves

4.2. Fitting Parameterizations

4.3. What Sets the Slope of Dust Attenuation Curves?

4.4. Variations in the 2175 Å Bump

4.5. Impact of Inclination Angle

5. Discussion

5.1. Does a Single Attenuation Prescription Apply?

5.2. Comparison to Other Models and Constraints

5.2.1. On the Role of Geometry

5.2.2. Bump Strengths

5.3. Comparison with Observations

5.4. On the Role of the Underlying Galaxy Formation Model

6. Conclusions and Summary

Appendix A: Octree Deposition

Appendix B: Resolution Tests

Footnotes