# A dust attenuation model inspired by the NIHAO-SKIRT-Catalog

Nicholas Faucher  Center for Cosmology and Particle Physics, Department of
Physics, New York University, 726 Broadway, New York, New York 10003, USA
Michael R. Blanton  Center for Cosmology and Particle Physics, Department of
Physics, New York University, 726 Broadway, New York, New York 10003, USA

###### Abstract

We use simulated attenuation curves from the NIHAO-SKIRT-Catalog to test the
flexibility of commonly used dust attenuation models in the face of the
variations expected from realistic star-dust geometries. Motivated by lack of
flexibility in these existing models, we propose a novel dust attenuation
model with three free parameters that can accurately recover the simulated
attenuation curves as well as the best-fitting curves from the commonly used
models. This new model is fully analytic and treats all starlight equally, in
contrast to two-component dust attenuation models. We use the parametrization
to investigate the relationship between the overall attenuation law shape and
the strength of the 2175 Å bump. Our results indicate variation in star-dust
geometry leads these features to correlate tightly, with grayer attenuation
curves having weaker bumps.

##  1 Introduction

Untangling the history of galaxies over cosmic time allows us to better
understand the physical process affecting the formation of stars and the
growth of galaxies. One path to investigating this history is spectral energy
distribution (SED) modeling of photometry or spectra of low redshift galaxies
to infer their star formation and chemical evolution histories (Conroy & Gunn,
2010; Conroy, 2013). An important uncertainty in SED modeling is dust, which
reddens and attenuates the light in a manner that depends on geometry, viewing
angle, dust mass, and dust composition (Hayward & Smith, 2015; Trčka et al.,
2020; Lower et al., 2020; Faucher et al., 2023).

This paper explores how well various commonly used dust attenuation models can
reproduce the attenuation curves predicted from radiative transfer post-
processing of high-resolution hydrodynamical galaxy simulations from the
Numerical Investigation of Hundred Astrophysical Objects (NIHAO) project (Wang
et al., 2015; Faucher et al., 2023). Our goals are to understand the potential
uncertainties in commonly used models and to develop techniques to reduce
these uncertainties. Although our exploration here of dust attenuation models
is based on simulations run to redshift zero, astronomers apply the dust
attenuation models we test at all redshifts, and therefore our results may be
applicable to high redshift as well.

We expect that SED modeling inferences of star formation rates (SFRs) will be
sensitive to the assumed attenuation curve. For example, the slope of the
ultraviolet SED yields information about the attenuation if the attenuation
curve is steep, but might not if it is relatively gray; assuming the wrong
attenuation curve models will yield wrong total attenuation inferences and
therefore wrong SFRs (Conroy et al., 2010). Similar issues can arise in the
analysis of optical photometry and spectra, affecting not only SFRs but also
stellar masses and other aspects of the star formation history. (Greener et
al., 2020). Salim & Narayanan (2020) show that the inferred physical
parameters of simulated galaxies can vary dramatically when changing the
assumed dust attenuation curve.

Since attenuation curves are not directly observable, many models have been
proposed to explain the observed photometry of galaxies in the local Universe.
In this work, we focus on the dust attenuation models of Calzetti et al.
(2000), Charlot & Fall (2000), Kriek & Conroy (2013), and Cardelli et al.
(1989), which we refer to as Calzetti, Power Law, Kriek and Conroy, and
Cardelli, respectively. Of these four dust attenuation models, only the
Calzetti and Power Law attenuation models do not contain a UV-bump centered
around 2175 Å. We also note that all of these models except for Calzetti are
two-component models, in which the light from young stellar populations
receives additional extinction from a power law curve with its own free
parameters.

To improve our understanding of the accuracy of SED modeling inferences, in
this paper we test the flexibility of existing dust attenuation models by
performing fits to attenuation curves in mock observations from the NIHAO-
SKIRT-Catalog (Faucher et al., 2023), which take into consideration the three-
dimensional distribution of stars and dust of NIHAO simulated galaxies,
realistic dust grain size distributions and chemical compositions following
The Heterogeneous dust Evolution Model for Interstellar Solids (THEMIS) dust
model (Jones et al., 2017), and scattering by dust into and out of the line of
sight using the Stellar Kinematics Including Radiative Transfer (SKIRT) (Camps
& Baes, 2020) Monte Carlo radiative transfer code.

Section 2 describes how the simulated attenuation curves were calculated, the
models that we compare to, and our fitting procedures. Section 3 presents the
results. Section 4 discusses the implications of the results, and we conclude
in Section 5.

##  2 Methods

In this section, we describe our sample of simulated attenuation curves, the
dust attenuation models we use, and our fitting procedure. We have released
Python implementations of all of the dust attenuation models used in this
work111https://github.com/ntf229/dustModels.

###  2.1 NIHAO-SKIRT-Catalog attenuation curves

In order to test the flexibility of dust attenuation models, we utilize mock
photometry from the NIHAO-SKIRT-Catalog (Faucher et al., 2023) which includes
the effects of realistic star-dust geometries from the NIHAO galaxy
simulations (Wang et al., 2015). In addition to performing full 3D radiative
transfer, this mock catalog utilizes subgrid post-processing recipes that
mitigate limitations in the temporal and spatial resolution of the simulations
in order to produce far-ultraviolet (FUV) to far-infrared (FIR) photometry
that statistically reproduce the colors of galaxies in the nearby Universe.
The unattenuated stellar spectra in both the NIHAO-SKIRT-Catalog and those
used in this work are generated from FSPS (Conroy et al., 2009; Conroy & Gunn,
2010; Foreman-Mackey et al., 2014) using the MIST isochrones (Paxton et al.,
2011, 2013, 2015; Choi et al., 2016; Dotter, 2016), the MILES spectral library
(Sánchez-Blázquez et al., 2006), and assuming a Chabrier stellar initial mass
function (IMF) (Chabrier, 2003).

###  2.2 Calzetti

The Calzetti (Calzetti et al., 2000) dust extinction curve is the simplest
model we explore in this work. In the most flexible case, it is controlled by
only two parameters: The overall scaling, which we set by the attenuation in
the V-band, $A_{\rm V}\in[0,10]$, and a variable fraction of unattenuated
stellar light, $f_{\rm no-dust}\in[0,1]$, where the indicated ranges are those
allowed in the fitting procedure. Allowing some of the stellar light to reach
the observer unattenuated leads to more gray attenuation curves. This
extinction law has been used in in many influential works, e.g. investigating
the environment-dependent star formation rate and mass evolution (Peng et al.,
2010), associating galaxy growth history with dark matter halo assembly
history (Behroozi et al., 2019), and estimating photometric redshifts for
James Webb Space Telescope catalogs (Castellano et al., 2022; Treu et al.,
2022; Adams et al., 2023).

###  2.3 Power Law

The Power Law extinction curve is a two-component model due to (Charlot &
Fall, 2000), meaning that the young (age$\,<\,10$ Myr) stellar populations are
attenuated by an additional power law extinction curve with its own free
parameters, before being added to and attenuated along with the rest of the
stellar populations. These two components are meant to represent the
attenuation from molecular clouds surrounding young stellar populations and
the attenuation from the diffuse interstellar medium (ISM). This attenuation
model has 6 free parameters: The overall scaling of the diffuse attenuation,
$A_{\rm V}\in[0,10]$, the power law index of the diffuse attenuation,
$p\in[-5,0]$, the overall scaling of the young attenuation, $A_{\rm
V,\,young}\in[0,10]$, the power law index of the young attenuation, $p_{\rm
young}\in[-5,0]$, and a variable fraction of unattenuated stellar light for
both the diffuse and young attenuation, $f_{\rm no-dust}\in[0,1]$, and $f_{\rm
no-dust,\,young}\in[0,1]$, respectively. This attenuation law was assumed in
the heavily used MPA-JHU catalog (Kauffmann et al. 2003; Tremonti et al. 2004;
Brinchmann et al. 2004) including studies of AGN host galaxies, the origin of
the mass-metallicity relation, the mass dependence of the star formation rates
of galaxies, and many others.

###  2.4 Kriek and Conroy

The Kriek and Conroy (Kriek & Conroy, 2013) extinction curve is also a two-
component model. The young stellar populations are attenuated by a power law
extinction curve with its own free parameters. Then, an additional layer of
“diffuse” attenuation is applied uniformly to the entire stellar population.
This model has 6 free parameters: The overall scaling of the diffuse
attenuation, $A_{\rm V}\in[0,10]$, the power law index of the diffuse
attenuation $p\in[-5,0]$, the overall scaling of the young attenuation,
$A_{\rm V,\,young}\in[0,10]$, the power law index of the young attenuation,
$p_{\rm young}\in[-5,0]$, and a variable fraction of unattenuated stellar
light for both the diffuse and young attenuation, $f_{\rm no-dust}\in[0,1]$,
and $f_{\rm no-dust,\,young}\in[0,1]$, respectively. The diffuse attenuation
has a 2175 Å UV bump with a strength set deterministically by the power law
slope, $p$, with steeper slopes corresponding to stronger UV-bumps. However,
the effective power law slope of this model is not given by $p$, but rather
$\sim(p-0.8)$ due to the power law term being multiplied by the Calzetti
curve, which also causes the resulting attenuation curves to deviate from the
power law shape even when excluding the UV bump term:

$A_{\lambda}=\,\frac{A_{V}}{4.05}\,\left(k\left(\lambda\right)+D\left(\lambda\right)\right)\left(\frac{\lambda}{5500\,\mathrm{\AA}}\right)^{p}$
(1)

where $k\left(\lambda\right)$ is the Calzetti cruve, $\lambda$ has units of Å,

$D\left(\lambda\right)\,=\,\frac{E_{b}\left(\lambda\,\Delta\lambda\right)^{2}}{\left(\lambda^{2}-\lambda_{0}^{2}\right)^{2}+\left(\lambda\,\Delta\lambda\right)^{2}}$
(2)

with $E_{b}\,=\,0.85-1.9\,p$, $\Delta\lambda=350\,\mathrm{\AA}$, and
$\lambda_{0}=2175\,\mathrm{\AA}$.

This model has been used to investigate the progenitors of ultra-massive
galaxies (Marchesini et al., 2014), the evolution of the mass-size relation
(Suess et al., 2019), and the growth of metals and stellar mass at high
redshift (Tacchella et al., 2022).

###  2.5 Cardelli

The Cardelli (Cardelli et al., 1989) extinction curve, also known as the Milky
Way extinction or CCM (Cardelli, Clayton, and Mathis) law, is also a two-
component model of the same type as above. In this case the young stellar
population attenuation is again a power law, and the diffuse attenuation is a
Milky Way extinction curve. In contrast to the original model used in Cardelli
et al. (1989), we use the model as incorporated in FSPS which includes a
variable UV-bump term. The model has 7 free parameters: The overall scaling of
the diffuse attenuation, $A_{\rm V}\in[0,10]$, the ratio of total to selective
absorption which characterizes the Milky Way extinction curve, $R_{\rm
V}\in[0,10]$, the strength of the 2175 Å UV-bump, $b_{\rm UV}\in[0,10]$, the
overall scaling of the young attenuation, $A_{\rm V,\,young}\in[0,10]$, the
power law index of the young attenuation, $p_{\rm young}\in[-5,0]$, and a
variable fraction of unattenuated stellar light for both the diffuse and young
attenuation, $f_{\rm no-dust}\in[0,1]$, and $f_{\rm no-dust,\,young}\in[0,1]$,
respectively. Although this model has been widely used to account for the
foreground extinction curve of light by the Milky Way’s ISM in the context of
extragalactic surveys (Schlegel et al. 1998 among literally thousands of
others), it has also been used to account for the dust attenuation of host
galaxies in SED modeling such as in the study of the progenitors of ultra-
massive galaxies (Marchesini et al., 2014) and in order to find correlations
between galactic spectral type and the best-fitting dust attenuation law
(Kriek & Conroy, 2013).

###  2.6 TEA Model

Our proposed model, which we call the THEMIS Effective Attenuation (TEA)
model, has the following functional form with 3 free parameters $A_{\rm
V}\in[0,10]$, $p\in[-5,0]$, and $b_{\rm UV}\in[0,10]$:

$A_{\lambda}=\,A_{\rm
V}\,\left(\frac{\lambda}{5500\,\mathrm{\AA}}\right)^{p}+\,b_{\rm
UV}\left\\{\left[1+{\rm erf}\left(s_{\rm
UV}\,f_{1}\left(\lambda\right)\right)\right]\left(1+\left(6f_{1}(\lambda)\right)^{2}\right)^{-3/2}\,-\frac{1}{4}\left(\frac{f_{2}(\lambda)^{5}}{f_{2}(\lambda)^{4}+1}+f_{2}(\lambda)\right)\right\\}$
(3)

where $\lambda$ has units of Å, $s_{UV}=6.95$,

$\displaystyle f_{1}(\lambda)$ $\displaystyle=$
$\displaystyle\log_{10}(\lambda)-3.312,$ (4) $\displaystyle f_{2}(\lambda)$
$\displaystyle=$ $\displaystyle 16\left(\log_{10}(\lambda)-3.285\right),$ (5)

and ${\rm erf}()$ is the error function. The first term is a power law with
variable strength ($A_{\rm V}$) and power-law slope ($p$). The second term is
a skewed bump function where the bump strength is controlled by the parameter
$b_{\rm UV}$ and the skewness is set by $s_{\rm UV}$. The last term scales
with the bump strength $b_{\rm UV}$ and adds attenuation to the bump while
subtracting attenuation from wavelengths slightly shorter than the bump. We
can see from Figure 1, which shows attenuation curves from the TEA model over
a range of parameters, that the value of $A_{\rm V}$ only corresponds to the
actual attenuation in the V-band ($\log_{10}(\lambda_{\rm
eff}\,/\,\mathrm{\AA})=3.74$) when $b_{\rm UV}=0$, since $A_{\rm V}$ is the
attenuation in the V-band for the power law component only and the UV-bump
extends to large wavelengths.

Equation (3) is not a dust extinction curve, it is a dust attenuation curve
that was designed to encapsulate the effects of the star-dust geometry from
the NIHAO simulated galaxies and the subgrid attenuation from star-forming
regions as modeled by the photo-ionization code MAPPINGS-III using Starburst99
SEDs (Leitherer et al., 1999) as implemented in SKIRT.

Figure 1:  Visualization of TEA model attenuation curves in magnitudes over a
range of parameters.

###  2.7 Fitting Procedure

We test the flexibility of these dust attenuation models by fitting them to
attenuation curves from the NIHAO-SKIRT-Catalog. We fit each dust model to
these simulated attenuation curves by minimizing the root mean squared error
(RMSE) of $A_{\lambda}$ in magnitudes on a logarithmically spaced grid of 75
wavelengths between $10^{3}$ and $10^{4}$ Å using the dual annealing global
minimization function (Xiang et al., 1997) implemented in scipy.optimize
(Virtanen et al., 2020), set to 2000 maximum iterations. In our fits, we will
test both with and without freedom in $f_{\rm no-dust}$ and $f_{\rm no-
dust,\,young}$. We also note that we are allowing more flexibility in these
dust models than is typically done in SED modeling. For instance, both Leja et
al. (2017) and Hahn & Melchior (2022) implement the Kriek and Conroy dust
attenuation model within the SED modeling frameworks Prospector and SEDFLOW,
respectively, but only allow the $A_{\rm V}$, $p$, and $A_{\rm V,\,young}$
attenuation parameters to vary.

The two-component dust models require a model for the star formation history
to define an attenuation curve. Here we use the most favorable choice to these
models, which is the actual star formation history in these NIHAO galaxies.
For this star formation history we use the chemical evolution histories (CEHs)
published in the NIHAO-SKIRT-Catalog, which give the star formation rate as a
function of age and metallicity. We separate the populations into those
greater and less than 10 Myr in age. Then we apply each dust attenuation model
component to the dust-free spectra of the young and old populations, and
construct an attenuated spectrum $f_{A}(\lambda)$. The attenuation curve is
then $-2.5\log_{10}f_{A}(\lambda)/f(\lambda)$, where $f(\lambda)$ is the
unattenuated spectrum.

##  3 Results

In Section 3.1, we explore the ability of each dust model to recover the
attenuation curves from dusty galaxies from the NIHAO-SKIRT-Catalog. In
Section 3.2, we see how well the newly proposed TEA model can recover the best
fits of the other dust models. Finally, we use the NIHAO-SKIRT-Catalog fits
from the TEA model to explore correlations between its parameters in Section
3.3.

###  3.1 Fitting attenuation models to the NIHAO-SKIRT-Catalog

Figure 2 shows the fits of the face-on orientation of the simulated galaxy
g5.02e11 for the the five dust models being considered in this work. In this
case, we allow all parameters to be free. In total, the Calzetti attenuation
model has 2 free parameters ($A_{\rm V}$ and $f_{\rm no-dust}$), the Power Law
attenuation model has 6 free parameters ($A_{\rm V}$, $p$, $A_{\rm
V,\,young}$, $p_{\rm young}$, $f_{\rm no-dust}$, and $f_{\rm no-
dust,\,young}$), the Kriek and Conroy attenuation model has 6 free parameters
($A_{\rm V}$, $p$, $A_{\rm V,\,young}$, $p_{\rm young}$, $f_{\rm no-dust}$,
and $f_{\rm no-dust,\,young}$), the Cardelli attenuation model has 7 free
parameters ($A_{\rm V}$, $R_{\rm V}$, $b_{\rm UV}$, $A_{\rm V,\,young}$,
$p_{\rm young}$, $f_{\rm no-dust}$, and $f_{\rm no-dust,\,young}$), and the
TEA model has 3 free parameters ($A_{\rm V}$, $p$, and $b_{\rm UV}$).

Figure 2:  Fits of dust attenuation models to face-on orientation of the
attenuation curve in magnitudes from the simulated galaxy g5.02e11, allowing
all parameters to be free in the fits. In this case, Calzetti has 2 free
parameters, Power Law has 6 free parameters, Kriek and Conroy has 6 free
parameters, Cardelli has 7 free parameters, and the TEA model has 3 free
parameters.

Figure 3 shows the root mean squared error (RMSE) in magnitudes as a function
of galaxy stellar mass between the dust model fits and the simulated
attenuation curves. In these fits we treat as free parameters $f_{\rm no-
dust}$ and $f_{\rm no-dust,\,young}$, representing the fractions of starlight
allowed to reach the observer unattenuated. The TEA model outperforms all
other models for nearly all galaxies and orientations. Note that these RMSEs
are calculated from the unnormalized attenuation curves, which causes the
general trend of increasing RMSEs with increasing stellar mass, since more
massive galaxies tend to have stronger dust attenuation. Figure 4 shows the
RMSE, offset in log-space by the RMSE of the TEA model, as a function of
galaxy stellar mass. This plot makes it easier to see how the different
models’ RMSEs compare to each other by subtracting away the scatter from the
differences in RMSEs between galaxies.

Figure 3:  Root mean squared errors (RMSE) of the dust attenuation model fits
for the face-on (left) and edge-on (right) orientations of the simulated
galaxies’ attenuation curves in magnitudes, allowing all parameters to be
free.

Figure 4:  Similar to Figure 3, now with the RMSEs shifted by the TEA model’s
RMSE for each galaxy.

Figure 5 shows the same fits as in Figure 2, but with the parameters $f_{\rm
no-dust}$ and $f_{\rm no-dust,\,young}$ fixed to zero, meaning all starlight
gets attenuated. In total, the Calzetti attenuation model has 1 free parameter
($A_{\rm V}$), the Power Law attenuation model has 4 free parameters ($A_{\rm
V}$, $p$, $A_{\rm V,\,young}$, $p_{\rm young}$), the Kriek and Conroy
attenuation model has 4 free parameters ($A_{\rm V}$, $p$, $A_{\rm
V,\,young}$, $p_{\rm young}$), the Cardelli attenuation model has 5 free
parameters ($A_{\rm V}$, $R_{\rm V}$, $b_{\rm UV}$, $A_{\rm V,\,young}$,
$p_{\rm young}$), and the TEA model has 3 free parameters ($A_{\rm V}$, $p$,
and $b_{\rm UV}$).

Figure 5:  Similar to Figure 2, now not allowing any fraction of the stellar
light to reach the observer unattenuated. In this case, Calzetti has 1 free
parameter, Power Law has 4 free parameters, Kriek and Conroy has 4 free
parameters, Cardelli has 5 free parameters, and the TEA model has 3 free
parameters.

Figure 6 shows the root mean squared error (RMSE) between the dust model fits
and the simulated attenuation curves without allowing a fraction of starlight
to reach the observer unattenuated. We can see that in this case, the TEA
model outperforms nearly all other models for all galaxies and orientations.
Figure 7 shows the RMSE, offset in log-space by the RMSE of the TEA model, as
a function of galaxy stellar mass. This plot makes it easier to see how the
different models’ RMSEs compare to each other by subtracting away the scatter
from the differences in RMSEs between galaxies.

Figure 6:  Similar to Figure 3, now not allowing any fraction of the stellar
light to reach the observer unattenuated.

Figure 7:  Similar to Figure 6, now with the RMSEs shifted by the TEA model’s
RMSE for each galaxy.

###  3.2 Fitting the TEA model to best-fitting common model attenuation curves

The TEA model is also flexible enough to reproduce the behavior of the other
model attenuation curves. To demonstrate this fact we fit the TEA model to the
other models, using for the latter the parameters from all of the best fits to
the NIHAO-SKIRT-Catalog.

Figure 8 shows the result of this process applied to the face-on orientation
of the simulated galaxy g5.02e11, using the models that allow for freedom in
$f_{\rm no-dust}$ and $f_{\rm no-dust,\,young}$. It may appear counter-
intuitive at first glance that the TEA model fit to the power law curve has a
strong UV-bump, since if we set $b_{\rm UV}$ to zero the TEA model gives
exactly a power law. However, we must keep in mind that the power law model
isn’t a simple power law, its shape is complicated by its two-component nature
as well as the flexibility to only apply attenuation to a fraction of the
young and old stellar populations. Given that the TEA model has neither of
these features, it uses the UV-bump to match the power law model’s curve as
well as possible. This same line of reasoning also applies to the TEA model
fits to the Kriek and Conroy curves.

Figure 8: Fits of the TEA model to commonly used dust model fits of the face-
on orientation of attenuation curve in magnitudes from the simulated galaxy
g5.02e11, allowing all parameters to be free in the fits. In this case,
Calzetti has 2 free parameters, Power Law has 6 free parameters, Kriek and
Conroy has 6 free parameters, Cardelli has 7 free parameters, and the TEA
model has 3 free parameters.

Figure 9 shows the root mean squared error (RMSE) between the TEA model fits
and the best fits from the commonly used attenuation models of the simulated
attenuation curves, allowing for a fraction of the stellar light to reach the
observer unattenuated in the commonly used attenuation models.

Figure 9:  Root mean squared errors (RMSE) of the TEA model fit to commonly
used model fits for the face-on (left) and edge-on (right) orientations of the
simulated galaxies’ attenuation curves in magnitudes.

Figure 10 shows the fits of the TEA model to the best fitting curves of the
other dust models to the face-on orientation of the simulated galaxy g5.02e11,
fixing $f_{\rm no-dust}=f_{\rm no-dust,\,young}=0$, i.e. that all stellar
light is attenuated in the commonly used attenuation models. Comparing the TEA
model fit to the power law curves between Figures 8 and 10, we can see that
taking away the power law model’s ability to only apply attenuation to a
fraction of the young and old stellar light yields TEA model fits with weaker
UV-bumps. We can understand this by recalling that the TEA model is using its
UV-bump to compensate for its inability to separately attenuate the light from
young and old stellar light (both figures), as well as its inability to only
apply attenuation to a fraction of the young and old stellar light (Figure 8
only).

Figure 10:  Similar to Figure 8, now not allowing any fraction of the stellar
light to reach the observer unattenuated. In this case, Calzetti has 1 free
parameter, Power Law has 4 free parameters, Kriek and Conroy has 4 free
parameters, Cardelli has 5 free parameters, and TEA model has 3 free
parameters.

Figure 11 shows the root mean squared error (RMSE) in magnitudes between the
TEA model fits and the best fits from the commonly used attenuation models of
the simulated attenuation curves, fixing $f_{\rm no-dust}=f_{\rm no-
dust,\,young}=0$, i.e. that all stellar light is attenuated in the commonly
used attenuation models.

Comparing Figures 9 and 11 to Figures 3 and 6, we can see that the TEA model
is able to reproduce the best fitting attenuation curves of the common models
better than the common models are able to reproduce the attenuation curves
from the NIHAO-SKIRT-Catalog.

Figure 11:  Similar to Figure 9, now not allowing any fraction of the stellar
light to reach the observer unattenuated.

###  3.3 Correlation Between TEA model parameters

In this section we explore the TEA model’s best-fitting parameters to the
NIHAO-SKIRT-Catalog attenuation curves. We first consider the correlation
between overall attenuation strength ($A_{\rm V}$) and the power law slope
($p$) shown in Figure 12. While there is a substantial amount of scatter in
the relation, there is a clear trend for more heavily attenuated galaxies to
have grayer (less steep) attenuation curves.

Figure 12:  Distribution of best-fitting $p$ and $\log_{10}(A_{\rm V})$
parameters of the TEA model to the NIHAO-SKIRT-Catalog. Here we only include
galaxies which have at least one orientation with an $A_{\rm V}$ value greater
than 0.05 magnitudes. We also include the equation for the best-fitting line
in the legend.

Motivated by the imposed relation between power law index and bump strength in
the Kriek and Conroy model (Kriek & Conroy, 2013), we explore how power law
index in the TEA model, $p$, relates to bump strength, $b_{\rm UV}$. Figure 13
shows the distribution of best-fitting $\log_{10}(b_{\rm UV}\,/\,A_{\rm V})$
and $p$ parameters from the NIHAO-SKIRT-Catalog. To ensure we’re limiting our
analysis to sufficiently dusty galaxies, we only include galaxies which have
at least one orientation with an $A_{\rm V}$ value greater than 0.2. We can
see that the fits exhibit a strong correlation between these parameters,
qualitatively consistent with the findings of Kriek & Conroy (2013). We note
that this correlation is primarily due to the variation in the star-dust
geometries of the simulated galaxies, since only one dust extinction curve is
implemented in the radiative transfer calculations of the diffuse interstellar
medium (ISM) in Faucher et al. (2023).

To the extent that this correlation holds true for galaxies in the real
universe, the best fitting line included in the legend of Figure 13 may be
used to reduce the dimensionality of the TEA model from 3 to 2 without
significant loss of practical flexibility by parameterizing the bump strength
in terms of $p$ and $A_{\rm V}$. Alternatively, one could rewrite the TEA
model’s bump term by replacing the free parameter $b_{\rm UV}$ with the offset
from this relation, effectively reducing the size of the model’s parameter
space while retaining the flexibility to vary the bump strength and power law
index independently.

Figure 13:  Distribution of best-fitting $\log_{10}(b_{\rm UV}\,/\,A_{\rm V})$
and $p$ parameters of the TEA model to the NIHAO-SKIRT-Catalog. Here we only
include galaxies which have at least one orientation with an $A_{\rm V}$ value
greater than 0.2. We also include the equation for the best-fitting line in
the legend.

Due to the differences in the UV-bump shapes and functional forms of the Kriek
and Conroy model (Equations 1 and 2) and the TEA model (Equation 3), it is
difficult to directly compare the power law index to bump strength
correlations between the two models. To circumvent this difficulty, we choose
to explore this relation through the maximum UV-bump height relative to the
overall attenuation, which we calculate by finding the maximum value of
$A_{\lambda}/A_{\rm V}$ within the bump regions for the two models
($3.23<\log_{10}(\lambda)<3.6$ for the TEA model and
$3.3<\log_{10}(\lambda)<3.4$ for the Kriek and Conroy model).

As noted in Section 2.4, the effective power law index of the Kriek and Conroy
model is offset from the value of its $p$ parameter by $\sim-0.8$ due to its
Calzetti component. To account for this offset, we perform linear fits of the
attenuation curves on a log-log scale, with their UV-bumps masked out, to find
their effective power law indices.

Figure 14 shows how bump height correlates with the power law index model
parameters and with the effective power law indices, for the TEA model fit to
the NIHAO-SKIRT-Catalog galaxies and the Kriek and Conroy model, over a grid
of power law index model parameters. Although there is some discrepancy in the
slope and vertical offset of the relation, the TEA model fits follow a similar
trend to the one imposed by the Kriek and Conroy model, with steeper curves
correlating with stronger UV-bumps.

Figure 14:  Comparison of the power law index to UV-bump height correlation
between the TEA model and the Kriek and Conroy model. For both models, we
calculate the effective power law index of the attenuation curves by
performing a linear fit on log-log scale with the UV-bump regions masked out.

##  4 Discussion

Although the figures in the previous section demonstrate that the TEA model is
able to accurately reproduce the attenuation curves from the NIHAO-SKIRT-
Catalog, it is important to remember that the catalog attenuation curves are
specific to the NIHAO simulations, the THEMIS dust model, and the sub-grid
modeling of photo-disassociation regions (PDRs) around young stellar
populations.

In particular, all of these attenuation curves include a bump around 2175 Å,
the UV-bump. A number of authors have concluded that the attenuation curves of
some galaxies in the real universe such as the Small Magellanic Clouds (SMC)
contain an extremely weak or absent UV-bump (Prevot et al., 1984; Gordon &
Clayton, 1998; Misselt et al., 1999; Gordon et al., 2003). Some of these
conclusions have been challenged; for instance Calzetti et al. (1994) find
that M33 has no UV-bump but recent work by Wang et al. (2022) suggests that it
does.

Regardless, we do know that, theoretically, different galaxies should exhibit
different UV-bump strengths. We have designed the TEA model to have a variable
UV-bump strength, enabling it to accurately reproduce attenuation curves with
strong bumps, weak bumps, or no bump at all. This requirement was the primary
motivation against using symbolic regression on the NIHAO-SKIRT-Catalog
attenuation curves to find the optimal functional form, since based on that
training data that process would have tended to only select for strong UV-
bumps.

Although we have decided to use only a single component for the TEA model,
this choice necessarily limits the resulting attenuation curves to be smooth
and fully described by Equation 3. This contrasts to with two-component
models, which can yield attenuation curves with small scale features due to
the differences between the underlying young and old spectra. In particular,
emission lines that tend to be stronger in young stellar populations will
receive extra attenuation in two-component models that cannot be reproduced by
single-component models. This limitation is unlikely to affect the flexibility
when fitting only to broadband photometry, but two-component models may have
an advantage when fitting moderate or high resolution spectra. Furthermore,
despite being single-component, the TEA model is able to reproduce the
attenuation curves of the NIHAO-SKIRT-Catalog, which does include additional
extinction for young stellar populations, better than the existing two-
component models explored in this work. This indicates that the flexibility to
have small scale features in the attenuation curves of two-component models is
less important that the large scale flexibility, at least relative to the loss
function used in this work (the RMSE described in Section 2.7).

In addition to being able to accurately reproduce attenuation curves that
account for complex star-dust geometries and birth cloud attenuation, the TEA
model has the practical advantage of being fully analytic. This fact allows
for convenient calculation of analytic gradients with respect to model
parameters which can be used in optimization methods such as stochastic
gradient descent (SGD). While this is also technically possible for the other
dust models explored in this work, there are extra steps involved due to
having two separate extinction curves for the young and old stellar
populations (Power Law, Kriek and Conroy, and Cardelli), as well as having
extinction curves with wavelength dependent functional forms (Calzetti, Kriek
and Conroy, and Cardelli). As the need for computationally efficient SED
modeling methods grows with the increase in astronomical data from large
surveys, this ability to perform automatic differentiation will likely be
extremely useful.

The analysis performed in this work is not one that can be done for galaxies
in the real universe, since attenuation curves are not a direct observable,
but can only be inferred in the context of full SED modeling. In the context
of SED modeling, attenuation curves are applied to stellar population spectra
in order to match the models to the observed data, which could be in the form
of broadband photometry and/or spectra. The parameter spaces of these models
are often high-dimensional and lead to to loss-landscapes with many local
minima, which provides further motivation for a dust attenuation model with as
few free parameters as possible, but with sufficient flexibility.

##  5 Conclusion

In this work, we have tested the flexibility of commonly used dust attenuation
models by using them to fit the attenuation curves from the NIHAO-SKIRT-
Catalog. While the catalog attenuation curves are limited by the underlying
simulations and dust models, they do include a physically motivated
distribution of grain sizes and chemical compositions, nonlinear effects
resulting from complex star-dust geometries, and sub-grid modeling of PDRs
around young stellar populations.

Motivated by the lack of flexibility of these commonly used attenuation
models, we have proposed a new dust attenuation model that is able to
accurately reproduce the catalog attenuation curves with only three free
parameters. The new model, which we call the TEA model, is fully analytic and
single-component, meaning its shape does not depend on the stellar spectra
that it is applied to. This contrasts with two-component models that first
attenuate the light coming from young stellar populations with a power law
extinction curve before attenuating the composite spectra. The resulting shape
of the global attenuation curve from two-component models will depend on the
shape of the young stellar spectra relative to the old stellar spectra. These
two characteristics make the TEA model both computationally efficient to model
and easily compatible with automatic differentiation. In the era of big data
in astronomy, these features make the TEA model particularly appealing in the
context of SED modeling optimization.

Using our parameterization, we have studied the relationship between the
overall shape of the attenuation curve and the strength of the 2175 Å bump. We
find that these two features to correlate well, in a sense that agrees
qualitatively (though not quantitatively) with the work of Kriek & Conroy
(2013)— grayer attenuation curves have smaller bump strengths. Because the
extinction law is constant within our models, we conclude that the physical
cause of this close correlation is due to the star-dust geometry—patchier
extinction leads to both grayer attenuation and smaller bump strengths.

We have made the TEA model, as well as Python versions of the commonly used
dust attenuation models explored in this work, publicly
available.222https://github.com/ntf229/dustModels

## References

  * Adams et al. (2023) Adams, N. J., Conselice, C. J., Ferreira, L., et al. 2023, MNRAS, 518, 4755, doi: 10.1093/mnras/stac3347
  * Behroozi et al. (2019) Behroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. 2019, MNRAS, 488, 3143, doi: 10.1093/mnras/stz1182
  * Brinchmann et al. (2004) Brinchmann, J., Charlot, S., White, S. D. M., et al. 2004, MNRAS, 351, 1151
  * Calzetti et al. (2000) Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682, doi: 10.1086/308692
  * Calzetti et al. (1994) Calzetti, D., Kinney, A. L., & Storchi-Bergmann, T. 1994, ApJ, 429, 582, doi: 10.1086/174346
  * Camps & Baes (2020) Camps, P., & Baes, M. 2020, Astronomy and Computing, 31, 100381, doi: 10.1016/j.ascom.2020.100381
  * Cardelli et al. (1989) Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245, doi: 10.1086/167900
  * Castellano et al. (2022) Castellano, M., Fontana, A., Treu, T., et al. 2022, The Astrophysical Journal Letters, 938, L15, doi: 10.3847/2041-8213/ac94d0
  * Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763, doi: 10.1086/376392
  * Charlot & Fall (2000) Charlot, S., & Fall, S. M. 2000, ApJ, 539, 718, doi: 10.1086/309250
  * Choi et al. (2016) Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102, doi: 10.3847/0004-637X/823/2/102
  * Conroy (2013) Conroy, C. 2013, ARA&A, 51, 393, doi: 10.1146/annurev-astro-082812-141017
  * Conroy & Gunn (2010) Conroy, C., & Gunn, J. E. 2010, ApJ, 712, 833, doi: 10.1088/0004-637X/712/2/833
  * Conroy et al. (2009) Conroy, C., Gunn, J. E., & White, M. 2009, ApJ, 699, 486, doi: 10.1088/0004-637X/699/1/486
  * Conroy et al. (2010) Conroy, C., Schiminovich, D., & Blanton, M. R. 2010, ApJ, 718, 184, doi: 10.1088/0004-637X/718/1/184
  * Dotter (2016) Dotter, A. 2016, ApJS, 222, 8, doi: 10.3847/0067-0049/222/1/8
  * Faucher et al. (2023) Faucher, N., Blanton, M. R., & Macciò, A. V. 2023, ApJ, 957, 7, doi: 10.3847/1538-4357/acf9f0
  * Foreman-Mackey et al. (2014) Foreman-Mackey, D., Sick, J., & Johnson, B. 2014, python-fsps: Python bindings to FSPS (v0.1.1), v0.1.1, Zenodo, Zenodo, doi: 10.5281/zenodo.12157
  * Gordon & Clayton (1998) Gordon, K. D., & Clayton, G. C. 1998, ApJ, 500, 816, doi: 10.1086/305774
  * Gordon et al. (2003) Gordon, K. D., Clayton, G. C., Misselt, K. A., Landolt, A. U., & Wolff, M. J. 2003, ApJ, 594, 279, doi: 10.1086/376774
  * Greener et al. (2020) Greener, M. J., Aragón-Salamanca, A., Merrifield, M. R., et al. 2020, MNRAS, 495, 2305, doi: 10.1093/mnras/staa1300
  * Hahn & Melchior (2022) Hahn, C., & Melchior, P. 2022, ApJ, 938, 11, doi: 10.3847/1538-4357/ac7b84
  * Hayward & Smith (2015) Hayward, C. C., & Smith, D. J. B. 2015, MNRAS, 446, 1512, doi: 10.1093/mnras/stu2195
  * Jones et al. (2017) Jones, A. P., Köhler, M., Ysard, N., Bocchio, M., & Verstraete, L. 2017, A&A, 602, A46, doi: 10.1051/0004-6361/201630225
  * Kauffmann et al. (2003) Kauffmann, G., Heckman, T. M., Tremonti, C., et al. 2003, MNRAS, 346, 1055, doi: 10.1111/j.1365-2966.2003.07154.x
  * Kriek & Conroy (2013) Kriek, M., & Conroy, C. 2013, ApJ, 775, L16, doi: 10.1088/2041-8205/775/1/L16
  * Leitherer et al. (1999) Leitherer, C., Schaerer, D., Goldader, J. D., et al. 1999, ApJS, 123, 3, doi: 10.1086/313233
  * Leja et al. (2017) Leja, J., Johnson, B. D., Conroy, C., van Dokkum, P. G., & Byler, N. 2017, ApJ, 837, 170, doi: 10.3847/1538-4357/aa5ffe
  * Lower et al. (2020) Lower, S., Narayanan, D., Leja, J., et al. 2020, ApJ, 904, 33, doi: 10.3847/1538-4357/abbfa7
  * Marchesini et al. (2014) Marchesini, D., Muzzin, A., Stefanon, M., et al. 2014, ApJ, 794, 65, doi: 10.1088/0004-637X/794/1/65
  * Misselt et al. (1999) Misselt, K. A., Clayton, G. C., & Gordon, K. D. 1999, ApJ, 515, 128, doi: 10.1086/307010
  * Paxton et al. (2011) Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3, doi: 10.1088/0067-0049/192/1/3
  * Paxton et al. (2013) Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4, doi: 10.1088/0067-0049/208/1/4
  * Paxton et al. (2015) Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15, doi: 10.1088/0067-0049/220/1/15
  * Peng et al. (2010) Peng, Y.-j., Lilly, S. J., Kovač, K., et al. 2010, ApJ, 721, 193, doi: 10.1088/0004-637X/721/1/193
  * Prevot et al. (1984) Prevot, M. L., Lequeux, J., Maurice, E., Prevot, L., & Rocca-Volmerange, B. 1984, A&A, 132, 389
  * Salim & Narayanan (2020) Salim, S., & Narayanan, D. 2020, ARA&A, 58, 529, doi: 10.1146/annurev-astro-032620-021933
  * Sánchez-Blázquez et al. (2006) Sánchez-Blázquez, P., Peletier, R. F., Jiménez-Vicente, J., et al. 2006, MNRAS, 371, 703, doi: 10.1111/j.1365-2966.2006.10699.x
  * Schlegel et al. (1998) Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525, doi: 10.1086/305772
  * Suess et al. (2019) Suess, K. A., Kriek, M., Price, S. H., & Barro, G. 2019, ApJ, 877, 103, doi: 10.3847/1538-4357/ab1bda
  * Tacchella et al. (2022) Tacchella, S., Finkelstein, S. L., Bagley, M., et al. 2022, ApJ, 927, 170, doi: 10.3847/1538-4357/ac4cad
  * Tremonti et al. (2004) Tremonti, C. A., Heckman, T. M., Kauffmann, G., et al. 2004, ApJ, 613, 898, doi: 10.1086/423264
  * Treu et al. (2022) Treu, T., Roberts-Borsani, G., Bradac, M., et al. 2022, ApJ, 935, 110, doi: 10.3847/1538-4357/ac8158
  * Trčka et al. (2020) Trčka, A., Baes, M., Camps, P., et al. 2020, MNRAS, 494, 2823, doi: 10.1093/mnras/staa857
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Wang et al. (2015) Wang, L., Dutton, A. A., Stinson, G. S., et al. 2015, MNRAS, 454, 83, doi: 10.1093/mnras/stv1937
  * Wang et al. (2022) Wang, Y., Gao, J., Ren, Y., & Chen, B. 2022, ApJS, 260, 41, doi: 10.3847/1538-4365/ac63c1
  * Xiang et al. (1997) Xiang, Y., Sun, D., Fan, W., & Gong, X. 1997, Physics Letters A, 233, 216