How Well Can We Measure Galaxy Dust Attenuation Curves? The Impact of the Assumed Star-dust Geometry Model in Spectral Energy Distribution Fitting

The most common way to model the impact of the star-dust geometry on a galaxy's dust attenuation curve is to allow a variable UV-optical slope:

$$\begin{eqnarray}\begin{array}{rcl}{\tau }_{\lambda } & = & -\mathrm{ln}(I/{I}_{0})\\ {S}_{\mathrm{UV}-\mathrm{opt}} & \equiv & {A}_{1500}/{A}_{V},\end{array}\end{eqnarray} \tag{ 1 }$$

where τ_λ is the opacity calculated via the ratio of I, the composite spectrum, to I₀, intrinsic stellar spectrum, and S_UV−opt is the ratio of attenuation at 1500 Å to the attenuation in the V band. τ_λ relates to the attenuation as A_λ ≈ 1.086 τ_λ .

For example, in studies such as those of Noll et al. (2009) and Salim et al. (2018) the authors allow the UV-optical attenuation curve slope to vary by multiplying the Calzetti et al. (2000) curve with a power-law term centered at the V band:

$$\begin{eqnarray}&&k{(\lambda )}_{\mathrm{variable}}\propto k{(\lambda )}_{\mathrm{Calzetti}}{\left(\displaystyle \frac{\lambda }{5500\,\mathring{\rm A} }\right)}^{\delta }.\end{eqnarray} \tag{ 2 }$$

where k(λ)_variable is the model attenuation curve, k(λ)_Calzetti is the base Calzetti et al. (2000) curve, and δ is the variable power-law slope (hereafter referred to as the UV-optical slope). Negative (positive) values of the slope term give attenuation curves that are steeper (shallower) than the Calzetti curve at λ < 5500 Å, allowing the attenuation curve to model galaxy-to-galaxy variations in attenuation curves manifested by, e.g., varying star-dust geometries. We show an example of this model in the top panels of Figure 1, where we plot SEDs and attenuation curves generated from simple stellar populations assuming a Calzetti et al. (2000) curve with a variable power-law slope. A modest range in the δ power-law slope parameter space results in a wide range of SEDs and attenuation curves even for fixed values of τ_V .

Zoom In Zoom Out Reset image size

We consider this implementation of the Calzetti curve model a uniform screen, since all stellar light is assumed to be attenuated equally by the same optical depth regardless of stellar age. While the variable UV-optical slope model enables us to approximate the complex star-dust geometries, it is limited in cases where the attenuation curve deviates from, e.g., a power-law function. We explore in the following section an alternative model that can accommodate more diverse attenuation curve shapes by removing the limitation of a power-law model.

As discussed above, the effect of star-dust geometry can be approximated by a variable attenuation curve slope, commonly parameterized as a power-law deviation from the fiducial Calzetti et al. (2000) curve. What are the limitations of such a model? To explore this question, we introduce a parameter, called f_unobscured, that allows deviations from a uniform screen model. This parameter controls the fraction of the composite stellar SED that is attenuated by the dust attenuation curve; we can think of (1 − f_unobscured) as a covering fraction in which nonzero values of f_unobscured result in a nonuniform screen model. The model composite spectrum with this parameter becomes

$$\begin{eqnarray}&&I={I}_{0}{k}_{\lambda }(1-{f}_{\mathrm{unobscured}})+{I}_{0}{f}_{\mathrm{unobscured}}.\end{eqnarray} \tag{ 3 }$$

To see the impact of this parameter on the attenuation curve model, we revisit Figure 1. In the bottom panels, we plot the SED and attenuation curves for a nonuniform screen model where we vary the fraction of unobscured stellar light. For simplicity, we plot the attenuation curves for two fixed values of δ, but in practice the slope would be allowed to vary as in the top panels. For a fixed δ slope, increasingly larger fractions of unattenuated light result in grayer attenuation curves that flatten out at small wavelengths.

For the uniform screen model parameterization adopted in this paper (Kriek & Conroy 2013), the shape of the attenuation curve is tied to the shape of the Calzetti et al. (2000) curve plus a 2175 Å bump feature. Tying the attenuation curve shape to the Calzetti curve represents a restriction on our model; the slope in one wavelength range of the attenuation curve is tied to the slope at 5500 Å. This is shown explicitly in Figure 2, where we plot the prior distributions for the UV slope and the UV-optical slope for the uniform screen and nonuniform screen models. The uniform screen model restricts the possible shapes of the attenuation curves. In contrast, the nonuniform screen allows for a range of UV slope values for fixed values of UV-optical slope. Variations in f_unobscured result in attenuation curves that are distinctly different than the curves generated by a uniform screen hinting that a more complex distribution of stars and dust changes the attenuation curve in ways that cannot easily be captured by varying a power-law slope alone.

Zoom In Zoom Out Reset image size

The addition of the f_unobscured parameter to a variable power-law slope model enables physically motivated flexibility without needing to make choices about the exact configuration of the star-dust geometry. In our analysis, we opt to remain agnostic about the specific configuration of the star-dust geometry by instead modeling the aggregate effect of the geometry on the attenuation curve. The goals of this paper are to show that the addition of f_unobscured introduces a necessary degree of freedom into the attenuation curve model such that we are better able to reproduce the attenuation curve of galaxies from a cosmological simulation.

In order to understand the current efficacy of dust attenuation modeling, we employ galaxies from the simba hydrodynamical cosmological simulation coupled with post-processing radiative transfer to produce synthetic SEDs that we model with prospector, a flexible SED modeling code. Below we describe the simba galaxy formation model, the 3D dust radiative transfer, and the prospector SED models used to fit the powderday synthetic SEDs.

We utilize the (25 h⁻¹ Mpc)³ simba ¹¹ volume with 2 × 512³ particles, resulting in a baryonic mass resolution of 1.4 ×10⁶ M_⊙. The simba galaxy formation model is a hydrodynamical cosmological simulation based on the gizmo gravity and hydrodynamics code (Hopkins 2015) and assumes a Planck cosmology. The fundamental models describing galaxy formation and evolution include (1) an H₂-based SFR; (2) a chemical enrichment model that tracks 11 elements from Type II supernovae (SNe), Type Ia SNe, and asymptotic giant branch (AGB) stars; (3) subresolution models for stellar feedback including contributions from Type II SNe, radiation pressure, and stellar winds; (4) subresolution models for feedback via active galactic nuclei (AGN) including a two-phase jet (high accretion rate ) and radiative (low accretion rate) model; and (5) a self-consistent dust model (Li et al. 2019), in which dust is produced by condensation of metals ejected from SNe and AGB stars and is allowed to grow by metal accretion, and be destroyed by thermal sputtering and astration processes.

We identify galaxies in the z = 0 simba snapshot using the caesar (Thompson 2014) 6D friends-of-friends galaxy finder based on the number of bound stellar particles in a system: a minimum of 32 stellar particles defines a galaxy. We then select galaxies with at least one star particle formed in the last 100 Myr to provide an SFR for comparison purposes. Our galaxy sample is shown in the top and middle panels of Figure 3, where we show the distribution of stellar mass, SFR, and attenuation curve properties.

Zoom In Zoom Out Reset image size

3.3.1. Generating Synthetic SEDs

3.3.2. Calculating Attenuation Curves

We calculate the attenuation curves for each galaxy following Equation (1). In the left panel of Figure 4, we show the distribution of simulated galaxy attenuation curves. The heatmap shows the density of curve parameter space and the red line denotes the median attenuation curve. Though the Weingartner & Draine (2001) extinction curve, shown in dark blue, is the same for all galaxies, the effective attenuation curves cover a wide range of parameter space. This wide range is primarily driven by differences in the star-dust geometry and the stellar age distributions.

Zoom In Zoom Out Reset image size

In the context of the discussion in Section 2, one way to parameterize the star-dust geometry is shown in the right panel of Figure 4, where we plot the distribution of f_unobscured, defined here as the fraction of stars in each galaxy that are unobscured by dust, for different τ_λ limits defining "unobscured." We calculate f_unobscured by computing the line-of-sight dust column density for each star particle in the simulated galaxies over one sight line. Though this is not exactly how f_unobscured is defined in Section 2, this paints an idea of how the stars and dust particles are coupled in our simulated galaxies.

In order to model the dust attenuation of the simba galaxies, we use the flexible state-of-the-art SED modeling code prospector (Leja et al. 2017; Johnson et al. 2021). Similar to powderday, prospector is based on the stellar population synthesis code fsps to generate simple stellar populations and convolved with the assumed star formation history (SFH) model to generate a composite stellar spectra. We fit the powderday SEDs sampled at 19 broadband filters shown in the bottom panel of Figure 3. Prospector relies on dynesty (Speagle 2020), a dynamic nested sampler that estimates Bayesian posteriors and evidences, to sample model parameter space. In our analysis, we closely follow the models and priors outlined in Lower et al. (2020), including the stellar mass-stellar metallicity relation of Gallazzi et al. (2005) and the six-component nonparametric SFH using a flexible modified Dirichlet prior (Betancourt & Girolami 2013; Leja et al. 2019) on the fractional stellar masses formed in each time bin.

Prospector includes several models for dust attenuation within the fsps framework. Notably, parameters specifying the attenuated fraction of stellar light—that is, the fraction of stars, both young (<10 Myr old) and old, that are unobscured by dust—can be varied. Meaning, prospector enables the dust attenuation model to be more diverse than the traditionally assumed uniform screen model. Because we do not model birth clouds in either the simba model (due to limited mass resolution) nor the powderday radiative transfer, we exclude any explicit age-based differential attenuation such as the, e.g., Charlot & Fall (2000) model. We instead focus only on the diffuse dust component. The three model variations we consider in our analysis are described below:

Fixed curve shape with uniform screen: Here, we assume the Calzetti et al. (2000) model. The slope of this curve is fixed and we implement this curve as a uniform screen, meaning all stellar light is impacted by this attenuation. The Calzetti et al. (2000) curve does not have the 2175 Å bump feature. We allow the normalization (V-band optical depth) to vary with a truncated Gaussian prior, centered at A_V = 0 with a width of 0.3 and truncated at A_V = 2.5.

Variable curve shape with uniform screen: Here, we use the attenuation curve parameterization of Kriek & Conroy (2013), which is based on the above Calzetti et al. (2000) relation, modulated by a power-law slope and includes a variable 2175 Å feature where the strength of the bump is tied to the slope deviation from the Calzetti et al. (2000) curve, in the sense that steeper curves have larger bump strengths. This slope allowed is to vary, ranging from −1.8 to 0.3 with a uniform prior relative to the Calzetti et al. (2000) curve slope as in Equation (2). The prior on V-band attenuation is the same as the fixed curve slope model.

Variable curve shape with nonuniform screen: This model is similar to the above variable slope model except we also allow the covering fraction of dust to vary. This parameter, f_unobscured, modulates the covering fraction of the diffuse dust and allows us to model nonuniformity in the star-dust geometry. Though the power-law slope and the fraction of unobscured stellar light are expected to be somewhat degenerate, the use of a Bayesian sampler means we can fully map these covariances and our uncertainties on the model parameters will be an accurate reflection of model degeneracies. The priors for the power-law slope and V-band normalization are the same as above while the prior on f_unobscured is uniform from 0 to 1.

As a first exercise, we fit the SEDs in the z = 0 simba volume using the three aforementioned dust attenuation models. In the top panels of Figure 5, we show the distribution of offsets between the inferred attenuation curves and the true attenuation curves across the UV and optical wavelength range probed by the broadband photometry. We plot these offsets as a function of inverse wavelength to focus on the UV portion of the curve, as the three models converge around 5500 Å to a median offset of ∼0. We generate the inferred attenuation curve from the marginalized parameter values. As we expect degeneracies between model parameters, namely metallicity, the attenuation power-law slope, and f_unobscured, we sample the full parameter posteriors to generate the attenuation curves instead of using the maximum likelihood values. The median offset from the true attenuation curve for 530 galaxies, represented by the black line in Figure 5, is smallest for the nonuniform screen model at all wavelengths greater than 1500 Å, where we have available UV data.

Zoom In Zoom Out Reset image size

We explore this further by examining which aspects of the curve see improvements in modeling in the bottom two panels of Figure 5, where we plot the magnitude of offsets for the UV-optical slope and attenuation at 1500 Å. As above, we calculate the offsets between the marginalized model parameters and the true values. The slope here is defined as the power-law deviation from the UV-optical slope of the Calzetti et al. (2000) curve. While the decrease in bias of the curve slope between the uniform screen and nonuniform screen models is only marginal, there is a significant decrease in offset between the inferred and true attenuation in the UV for the nonuniform screen model. This implies that while the variable power-law slope uniform screen model can infer the shape of the curve between 1500 Å and the optical regime, it significantly underestimates the attenuation in the UV with a median bias of $-{0.3}_{-0.30}^{+0.24}$ magnitudes compared to $-{0.17}_{-0.21}^{+0.15}$ for the nonuniform screen model. We interpret this as a consequence of the limited capability of a power-law transformation of the Calzetti et al. (2000) curve to match the diversity of the simba attenuation curves. Supporting this, in Figure 6, we plot the the UV slope inferred by the two variable slope models as a function of UV-optical slope (similar to Figure 2) which demonstrates the inability of the uniform screen model to cover the parameter space occupied by the simba galaxies (shown in gray). The tight relation between the two slopes for the uniform screen model is a consequence of the fact that the shape of the attenuation curve is allowed to vary only as a power-law deviation from the Calzetti et al. (2000) curve shape.

Zoom In Zoom Out Reset image size

Lastly, we note that all three attenuation curve models also appear to significantly underestimate the strength of the 2175 Å bump feature. As shown in the left panel of Figure 4, the average simulated attenuation curve features the 2175 Å bump, as does the underlying extinction curve used to generate the powderday SEDs. Because the attenuation curve model with a fixed slope is based on the Calzetti et al. (2000) curve, this deficit in the 2175 Å bump strength is to be expected as the Calzetti et al. (2000) curve does not include this feature. The two variable slope models are based on the Kriek & Conroy (2013) curve, where the strength of the bump feature is a function of the slope of the curve. Kriek & Conroy (2013) found that galaxies exhibiting steeper attenuation curves had smaller bump features. The bump strengths the authors measured are factors of 2–3 smaller than the Milky Way bump strength for sight lines measured by Valencic et al. (2004). Thus, compared to the powderday attenuation curves, the bump strengths of the Kriek & Conroy (2013) model at similar slopes are 1.5× smaller. Therefore, these deficits in the inferred 2175 Å bump feature are also expected.

A challenging aspect of SED modeling is the necessity to model each component (stars, metallicity, dust, etc.) simultaneously and untangle the numerous degeneracies associated with each aspect to infer the physical properties of a galaxy. From the previous section, we demonstrated that the nonuniform screen attenuation model successfully inferred the shape of the true dust attenuation curve for a majority of galaxies (with a median bias in the near-UV regime of 0.17 magnitudes).

How does this translate to the model's ability to infer the physical properties of these galaxies? In Figure 7 we plot the cumulative distribution of offsets for two inferred properties: the galaxy SFR and stellar metallicity for the three attenuation models. For both the SFR and galaxy metallicity, the increased accuracy of the modeled attenuation curves does translate to a modest increase in accuracy. This is especially true for the stellar metallicity that sees a decrease in bias of 0.18 dex on average. The SFRs are modestly improved with a decrease in bias of 0.12 dex. To relate these improvements in physical properties to the attenuation curve model improvements, in Figure 8 we plot the offset in UV attenuation as a function of SFR offset and stellar metallicity offset and for SFR and metallicity together. For galaxies fit with the uniform screen model, a large bias in inferred SFR or metallicity is matched by a similarly large bias in the UV attenuation. The spread in offsets for all three properties imply that the uniform screen model struggles to model these SED components simultaneously—at most, only two properties can be accurately inferred.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As we explored in Lower et al. (2020), the dominant source of uncertainty when inferring galaxy SFRs and stellar masses from SED modeling is the form of the assumed SFH, namely, that relatively simple models for the SFH can bias the inferred physical properties. In this paper, we vary only the dust attenuation model between SED fits and note only negligible differences in the inferred stellar masses between the three models, but modest improvements in the inferred SFR and significant improvements in the stellar metallicity, indicating that the ability to correctly infer the shape of the dust attenuation curve is vital in correctly measuring galaxy physical properties.

By adding a single parameter to a flexible slope dust attenuation curve model used in SED modeling, we have shown that we can reasonably reproduce the true attenuation curves of model galaxies, leading to improved estimates for galaxy physical properties such as the SFR. This parameter allows a nonzero fraction of stellar light to be unattenuated by dust in the diffuse ISM, resulting in attenuation curves that are functionally different from modulating the power-law slope alone. A nonuniform screen model accommodates attenuation curves that show deviations from a power-law slope, due to some stars seeing different optical depths than other stars. Coupled with a variable power-law slope, we are better equipped to model galaxy-to-galaxy variations in star-dust geometries that manifest as changes to the shape of the dust attenuation curve.

However, this nonuniform screen model is still only an approximation of the impact of star-dust geometry on a galaxy's attenuation curve. For instance, a drawback of the our current implementation of f_unobscured is that we neglect any variation of the obscured fraction as a function of stellar age (and effectively wavelength). The light-weighted fraction of unobscured stellar light will therefore be different than the mass-weighted fraction. Narayanan et al. (2018) demonstrated that allowing larger fractions of older stars to remain unobscured will result in steeper normalized dust attenuation curves while larger fractions of unobscured young stars will result in flatter attenuation curves. For the model presented in this study, we are unable to differentiate between the impact of unobscured old stellar populations versus the impact of unobscured young stellar populations on the model attenuation curve. This means that we cannot determine which stars are unobscured and where this decoupling happens within the galaxy. In essence, while f_unobscured is not explicitly the star-dust geometry, it aggregates the distribution of stellar obscuration fractions in a galaxy into an approximate bolometric covering fraction.

To overcome these model limits, future implementations of dust attenuation curves in SED models can benefit from a nonparametric treatment as is done for the SFH model in this work and developed in Leja et al. (2017, 2019). A nonparametric dust attenuation model would consist of several piecewise functions whose normalizations could vary independently, effectively accommodating any dust attenuation curve shape. Such a model, when properly constrained by full SED wavelength coverage and carefully chosen informative priors, would be less biased than the approaches taken in this paper since no assumptions regarding dust properties or star-dust geometries would be made. While the use of parametric models may make SED fitting results precise, as is also the case for SFH modeling, the limited nature of those models means the results will potentially be precise and wrong.

Though the results presented in Section 4 demonstrate an increased ability to constrain galaxy attenuation curves and physical properties with the nonuniform screen model, adding further complexity to an already multivariate and physically complex model runs the risk of unconstrained results and overfitting. This is compounded by our reliance on broadband photometry, which is relatively lower information than provided by, e.g., spectroscopy. The nonuniform model, or any dust attenuation model that is more complex than the typical variable slope model, could be more flexible than the data warrants and thus our results would be dominated by the priors and less by the likelihood.

The use of dynesty, and more generally Bayesian inference, means we can fully map the model parameter covariances and understand whether we suffer from overfitting or underfitting from the use of a certain model. One way to measure this is by measuring the evidence in favor of one attenuation model over another. Following the methodology of Salmon et al. (2016), we quantify this with the posterior odds, comparing the posterior marginalized over all model parameters of one dust attenuation model (M) to another given the available data (D):

$$\begin{eqnarray}&&\mathrm{posterior}\,\mathrm{odds}=\frac{p({M}_{1}| D)}{p({M}_{2}| D)}.\end{eqnarray} \tag{ 4 }$$

Given Bayes' theorem, p(M_i ∣D) can be broken down into

$$\begin{eqnarray}&&{\rm{posterior}}\,{\rm{odds}}=\displaystyle \frac{p(D| {M}_{1})p(\theta {M}_{1})}{p(D| {M}_{2})p(\theta {M}_{2})},\end{eqnarray} \tag{ 5 }$$

where the first ratio is the Bayes factor and represents the integral of the posterior density over the parameter space for each model and the second ratio is the prior odds for model parameters θ. The marginal likelihoods in the Bayes factor can be written as the integral over all parameter space for each model:

$$\begin{eqnarray}&&p(D| {M}_{i})={\int }_{\theta }p(D| \theta ,{M}_{i})p(\theta | {M}_{i})d\theta ,\end{eqnarray} \tag{ 6 }$$

which is also called the evidence. Thus, to calculate the posterior odds of two models, we take the ratio of their evidence and their prior probabilities, all of which are accessible with dynesty nested sampling.

Lastly, we use the criteria outlined in Kass & Raftery (1995) to denote the strength of the evidence toward one model or another. These criteria use the natural logarithm of the Bayes factor ($\mathrm{ln}({{\rm{B}}}_{12})$) and place the posterior odds in categories of strong, moderate, or weak preference for one model over another; $\mathrm{ln}({{\rm{B}}}_{12})\gt 10$ denotes a strong preference for model one and $\mathrm{ln}({{\rm{B}}}_{12})\lt 2$ denotes a weak preference for model one. A preference does not necessarily point to one model being correct and the other incorrect. Rather, it describes the evidence against the opposing model, supporting the null hypothesis that the other model is correct.

We calculate the posterior odds for the two variable slope attenuation curve models and evaluate the preference for one model over another for each galaxy using the criteria outlined in Kass & Raftery (1995). In Figure 9, we plot the fraction of galaxies that fall into each category (strong, moderate, weak preference) for either model. The fraction of galaxies whose evidence does not prefer either model are labeled inconclusive. A majority of galaxies (61%) have Bayesian evidence that points to a preference for the nonuniform screen model while less than 9% of galaxies have evidence that points toward the uniform screen model. Because Bayesian evidence is highly sensitive to the model prior space, we also fit SEDs with a more informative prior on f_unobscured based on the distribution in Figure 4. We find that the preference for the nonuniform screen model has fallen to 53% of galaxies with 14% preferring the uniform screen model. We interpret these results as (1) the more flexible model is preferred by a majority of galaxies with little preference (over inconclusive preference) for the less flexible model and (2) that while the prior space has moderate influence on the preference toward the nonuniform screen model there are sufficient constraints from the data for such a preference, alleviating some dangers of overfitting by a more complex model than warranted by the data.

Zoom In Zoom Out Reset image size

Finally, we consider the need for FIR constraints on our ability to accurately derive dust attenuation curves from SED fitting. Specifically, FIR data can break the reddening degeneracy between age and dust. Our results presented thus far have benefited from a fully sampled SED. We therefore investigate the efficacy of these techniques for observations that do not benefit from the availability of FIR data.

To address this, we fit our mock SEDs without the Spitzer and Herschel bands in the FIR and compare the inferred attenuation curves in Figure 10. We show the results from our fiducial model (nonuniform screen, fully sampled SED including FIR photometry) in green, while the results from fitting with no FIR data are shown in red. We note that while the average bias between the true and inferred attenuation curves does not change, there is a significant increase in the scatter of this offset at fixed wavelength, implying a decreased ability to correctly infer the attenuation for an individual galaxy.

Zoom In Zoom Out Reset image size

Because only the total infrared luminosity is needed to constrain the amount of dust attenuation in the UV and optical regime, one solution is to provide L_IR as an input to the prospector fit in a similar fashion as done by Salim et al. (2018) with the code cigale. In Figure 10, the fits performed with L_IR as a constraint are given in blue. The results show a significant decrease in scatter at a fixed wavelength: at 2300 Å the distribution widths are 0.64, 0.25, and 0.21 magnitudes for the no FIR, L_IR, and full photometry fits, respectively. From this, we conclude that while a fully sampled FIR SED provides the best constraints for modeling the dust attenuation in a galaxy, approximating this with only the infrared luminosity sufficiently decreases the scatter in the inferred attenuation curves from fits with no FIR constraints.