# A Bayesian Population Model for the Observed Dust Attenuation in Galaxies

Gautam Nagaraj Department of Astronomy & Astrophysics, The Pennsylvania State
University, University Park, PA 16802, USA  Institute for Gravitation and the
Cosmos, The Pennsylvania State University, University Park, PA 16802, USA
Center for Computational Astrophysics, 162 Fifth Avenue, New York, NY, 10010,
USA  John C. Forbes Center for Computational Astrophysics, 162 Fifth Avenue,
New York, NY, 10010, USA  Joel Leja Department of Astronomy & Astrophysics,
The Pennsylvania State University, University Park, PA 16802, USA  Institute
for Computational & Data Sciences, The Pennsylvania State University,
University Park, PA, USA  Institute for Gravitation and the Cosmos, The
Pennsylvania State University, University Park, PA 16802, USA  Daniel Foreman-
Mackey Center for Computational Astrophysics, 162 Fifth Avenue, New York, NY,
10010, USA  Christopher C. Hayward Center for Computational Astrophysics, 162
Fifth Avenue, New York, NY, 10010, USA

###### Abstract

Dust plays a pivotal role in determining the observed spectral energy
distribution (SED) of galaxies. Yet our understanding of dust attenuation is
limited and our observations suffer from the dust-metallicity-age degeneracy
in SED fitting (single galaxies), large individual variances (ensemble
measurements), and the difficulty in properly dealing with uncertainties
(statistical considerations). In this study, we create a population Bayesian
model to rigorously account for correlated variables and non-Gaussian error
distributions and demonstrate the improvement over a simple Bayesian model. We
employ a flexible 5-D linear interpolation model for the parameters that
control dust attenuation curves as a function of stellar mass, star formation
rate (SFR), metallicity, redshift, and inclination. Our setup allows us to
determine the complex relationships between dust attenuation and these galaxy
properties simultaneously. Using Prospector fits of nearly 30,000 3D-HST
galaxies, we find that the attenuation slope ($n$) flattens with increasing
optical depth ($\tau$), though less so than in previous studies. $\tau$
increases strongly with SFR, though when $\log~{}{\rm SFR}\lesssim 0$, $\tau$
remains roughly constant over a wide range of stellar masses. Edge-on galaxies
tend to have larger $\tau$ than face-on galaxies, but only for
$\log~{}M_{*}\gtrsim 10$, reflecting the lack of triaxiality for low-mass
galaxies. Redshift evolution of dust attenuation is strongest for low-mass,
low-SFR galaxies, with higher optical depths but flatter curves at high
redshift. Finally, $n$ has a complex relationship with stellar mass,
highlighting the intricacies of the star-dust geometry. We have publicly
released software for users to access our population model.

Hierarchical models (1925), Galaxy evolution (594), Spectral energy
distribution (2129), High-redshift galaxies (734)

††facilities: HST (WFC3), Spitzer (MIPS)††software: pymc3 (Salvatier et al.,
2016), AstroPy (Astropy Collaboration et al., 2013, 2018), SciPy (Jones et
al., 2001–; Virtanen et al., 2020), FSPS (Conroy et al., 2009; Conroy & Gunn,
2010), Prospector (Leja et al., 2017, 2019), dynesty (Speagle2020)

##  1 Introduction

One of the least understood components of baryonic matter in galaxies is dust.
While constituting only a negligible fraction of the total mass (e.g., Driver
et al., 2018), dust is important in interstellar medium (ISM) evolution and
chemistry, including cooling and catalysis of molecule formation. It plays an
important role in star and planet formation and reprocesses up to 30% of the
bolometric luminosity of a galaxy (e.g., Draine, 2003a, b).

Perhaps the most influential challenge dust presents to astrophysical studies
is the significant amount of absorption and scattering of light it causes,
especially in the ultraviolet (UV)—in a galaxy setting we call this dust
attenuation. An estimate of these effects is needed in order to measure
several crucial physical parameters of galaxies, including star formation
rates (SFRs) In fact, studies of galaxies at high redshift find the effects of
dust to be one of the most important sources of uncertainty (e.g., Bouwens et
al., 2012; Finkelstein et al., 2012; Oesch et al., 2013).

Many observational studies at high redshift have employed the “Calzetti Law”
(Calzetti et al., 2000) as a model for the dust attenuation curve (see §3 for
description of attenuation laws) as a function of wavelength (e.g., Erb et
al., 2006; Daddi et al., 2007; Mannucci et al., 2010; Peng et al., 2010b;
Elbaz et al., 2011; Madau & Dickinson, 2014) and/or used empirical relations
between observed and intrinsic UV slope from Calzetti et al. (1997); Meurer et
al. (1999); Calzetti (2001) to infer dust attenuation from UV spectra (e.g.,
Shapley et al., 2001; Bouwens et al., 2015). The aforementioned flagship
studies provide simple methods to approximate dust attenuation based on ample
observations in the local universe, making them quite useful.

However, the non-universality of dust attenuation curves has been known for
decades. While Stecher (1965) found evidence for a 2175 Å feature in
extinction curves of the Milky Way and Large Magellanic Cloud (a feature not
found in the Calzetti Law), Prevot et al. (1984) found a steeper extinction
curve with no evidence for the bump in the Small Magellanic Cloud. Witt &
Gordon (2000) showed through radiative transfer modeling that clumpier dust
distributions lead to flatter attenuation curves. Observational studies have
found that galaxies exhibit a wide variety of dust attenuation curves (e.g.,
Burgarella et al., 2005; Noll et al., 2009; Wild et al., 2011; Buat et al.,
2012; Arnouts et al., 2013; Kriek & Conroy, 2013; Reddy et al., 2015; Salim et
al., 2016; Leja et al., 2017; Salim et al., 2018). Furthermore, the details of
dust attenuation for a given galaxy have been shown to greatly influence the
measurement of other properties, such as SFRs and histories (e.g., Kriek &
Conroy, 2013; Reddy et al., 2015; Shivaei et al., 2015; Salim et al., 2016;
Salim & Narayanan, 2020).

Several studies have looked into ensemble trends in attenuation curves, both
in the local universe (e.g., Burgarella et al., 2005; Salim et al., 2016; Leja
et al., 2017; Salim et al., 2018) and at high redshift (e.g., Arnouts et al.,
2013; Kriek & Conroy, 2013; Salmon et al., 2016; Tress et al., 2018). These
studies have greatly improved our empirical understanding of dust attenuation,
but there is still progress to be made on understanding the evolution of
attenuation and the complex interplay among variables (Salim & Narayanan,
2020, and references therein), which is something we begin to tackle in this
work.

The measurement of dust attenuation curves from photometry and/or spectra is
done through empirical or modeling approaches (see Salim & Narayanan 2020 for
a detailed summary of the literature). Empirical methods often use templates
(e.g., Calzetti et al., 1994, 1997, 2000; Johnson et al., 2007; Reddy et al.,
2015; Battisti et al., 2016, 2017), in which galaxy spectra are binned by
Balmer decrement, infrared excess, or other relevant quantities, and averaged
in each bin to create templates. The template corresponding to little or no
dust is used to calibrate attenuation curves. While template empirical methods
are useful for learning about broad trends in galaxy populations, they
struggle with the complex correlations between parameters that are present in
galaxies.

This leads us to the modeling approach, which is dominated by spectral energy
distribution (SED) analysis (see Walcher et al. 2011 and Conroy 2013 for
reviews on this subject). SED modeling integrates various simple stellar
populations—collections of stars created simultaneously from (single)
molecular clouds—over a star formation history (SFH) and computes the
integrated SED. Of course, to match the observed SED, dust attenuation must be
taken into account. The goal of the SED model is to find the dust attenuation
(and possibly emission if the mid- and/or far-IR is included), SFH, and
metallicity parameters that are consistent with the data.

While SED fitting is powerful in its ability to extract physical parameters of
a galaxy based on photometry and/or emission line strengths or spectra, it
faces many challenges. Relevant to this project; differing ages,
metallicities, and dust attenuation curves can produce similar effects in
spectra, causing degeneracies (e.g., Conroy, 2013; Santini et al., 2015).

Bayesian techniques, especially parameter space exploration methods like
Markov Chain Monte Carlo (MCMC) algorithms, can help highlight the covariances
among parameters. Such approaches have become increasingly common in SED
modeling efforts over the past two decades, including CIGALE (Noll et al.,
2009; Boquien et al., 2019), GALMC (Acquaviva et al., 2011), Prospector (Leja
et al., 2017; Johnson et al., 2021), MCSED (Bowman et al., 2020), and others.

While SED fitting is an effective way to derive properties of individual
galaxies, the process of determining trends in ensembles of galaxies based on
individual fits is often limited to binning properties by their maximum
likelihood solutions and inspecting the results (e.g., Salim et al., 2018;
Nagaraj et al., 2021b, a) due to the computational expenses of the
alternative. In reality, observed SEDs are often noisy, and also, fitted
parameters are highly correlated and subject to degeneracies, as mentioned
earlier. Including proper treatment of correlated errors will lead to more
accurate inferences and avoid potential methodological bias (see Foreman-
Mackey et al. 2014 for a detailed discussion). The method we adopt to account
for covariances is population modeling, as portrayed in Figure 1.

The base layer of the pyramid represents observations. Technically, converting
raw telescope data to photometry and spectroscopy would require another one or
two pyramid layers, but for the scope of this paper we begin at the reduced
data layer. In this work, the observations come from the 3D-HST photometric
catalogs (Skelton et al., 2014). The middle layer represents measurements of
physical properties of galaxies, which are directly input into our dust
attenuation models (the top layer). In this work, we use Prospector parameter
fits of 3D-HST galaxies.

As a technical clarification, in this work, we use the arrows in the upward
direction only, making our model strictly a population model. If we used the
constraints on the population “hyperparameters” to inform better fits of
individual galaxies in the sample, it would be a true hierarchical model.

Figure 1: Schematic showing the structure of a simple hierarchical model—the
combination of all arrows to link the top and bottom layer—in the context of
galaxy evolution. In this diagram, the ensemble property $X$ (e.g., a dust
attenuation parameter) requires information about stellar mass and SFR of all
of the individual galaxies (all arrows between the middle and top layers).
Meanwhile, the information about stellar mass and SFR for individual galaxies
is derived from the data, including spectroscopy and photometry, reflected by
all arrows between the bottom and middle layers.

Only a few studies have conducted hierarchical studies of dust involving
galaxy SEDs (Kelly et al., 2012; Galliano, 2018). These focus on just dust
emission. In our study, we present the first population Bayesian modeling of
dust attenuation. As described in §3, we model dust attenuation curves as a
function of physical parameters of galaxies using posterior samples from
Prospector fits of nearly 30,000 galaxies at a wide range of redshifts:
$0.5<z<3.0$. As our galaxy sample is mass-complete, we do not need selection
function modeling to make population-level inferences.

Through this effort, we highlight more complex trends in dust attenuation,
pointing to gaps in our theoretical understanding. Our models, derived through
statistically rigorous means, can be used by modelers who do not wish to incur
the expense or systematic uncertainties of dust radiative transfer modeling.
They can also be used to strengthen the prior for future SED studies of
individual galaxies, allowing for tighter posterior constraints.

In §2 we describe the data we are using in our population model. This includes
both observations (photometry and grism data) and measurements of physical
parameters of galaxies using SED fitting. In §3, we describe the hierarchical
Bayesian framework for our dust attenuation model and describe results from
the tests used to demonstrate its accuracy. In §4, we share our results when
applying the population Bayesian model to the Prospector posterior samples and
compare them to what we get when simply taking posterior medians for parameter
values. We discuss some of the scientific statements we can make from the
fitted models in §5. In §6, we describe the public package we have created for
users to adopt our models. Finally in §7 we summarize our findings and suggest
steps for the future.

##  2 Data

The 3D-HST Treasury program (Brammer et al., 2012; Momcheva et al., 2016) used
the G141 grism on the WFC3 camera on the Hubble Space Telescope (HST) to
measure spectra of galaxies in 625 arcmin2 of the CANDELS fields (Grogin et
al., 2011; Koekemoer et al., 2011), a set of five highly studied patches of
sky with a wealth of available observations.

Skelton et al. (2014) used point spread function (PSF) matching techniques to
properly combine photometry from 147 data sets at wavelengths $0.3-8.0$ $\mu$m
for $\sim 200,000$ galaxies over a wide variety of redshifts. Momcheva et al.
(2016) used state-of-the-art reduction techniques and grism-fitting software
to calculate emission line fluxes and redshifts for nearly 80,000 galaxies.
The most important lines included (at different redshift intervals) are [O II]
$\lambda 3727$, H$\beta$, [O III] $\lambda 5007$, and H$\alpha$. For example,
H$\alpha$ is in the grism wavelength range (1.08-1.67 $\mu$m) at redshifts
$0.64<z<1.54$ whereas [O III] $\lambda 5007$ is present at $1.90<z<3.48$.

The data from the 3D-HST catalogs represent the bottom layer of the pyramid in
Figure 1. The actual inputs (middle layer) to our population model are the
values of the physical parameters like stellar mass and SFR that are derived
from the data catalogs via SED fitting. By including all posterior samples of
the desired input variables from their SED fitting results, we are properly
accounting for the uncertainties in the photometry as well as the correlations
between SED-fitted parameters. The formalization of the process of creating a
population model is presented in §3.1.

For our project, we have used Prospector (Leja et al., 2017; Johnson et al.,
2021) fits of 3D-HST galaxies because the combination of ample broadband data
and sophisticated SED modeling helps yield more accurate parameter fits.
Prospector has undergone a number of posterior predictive checks to verify it
produces accurate results, such as predicting the observed H$\alpha$
luminosities using constraints only from the photometry, (Leja et al., 2017)
and removing the discrepancy between the growth of the stellar mass density
and the observed cosmic SFR density (Leja et al., 2019). By employing a non-
parametric SFH, it allows more flexibility; indeed, the accumulation of more
mass at earlier times, which would not be possible in standard parametric
approaches, helped bridge the aforementioned disconnect (Leja et al., 2019).
In addition, Prospector allows for custom or informative priors, permitting
greater inclusion of physics and empirical results in the Bayesian framework.

Leja et al. (2019, 2020) used Prospector to fit 63,427 objects from the 3D-HST
catalogs at $0.5<z<3.0$. Of those, $\sim 30,000$ belong to a mass-complete
regime in 3D-HST (based on the completeness estimation by Tal et al. 2014). In
this work, we use the posterior samples from the mass-complete galaxies’ fits.
Restricting our analysis to the regime of completeness in 3D-HST allows us to
make statements about a complete galaxy population without having to worry
about selection effects.

Specifically, the stellar mass, sSFR, metallicity, and dust attenuation
measurements that we use in our modeling process are posterior samples from
Prospector fits of the galaxies. Redshift measurements come from the 3D-HST
team’s analysis of grism spectra (Momcheva et al., 2016), though $\sim 70$% of
sources have only photometric redshifts available. These are not allowed to
vary; in other words, each galaxy has a fixed redshift. Finally, axis ratios
were derived using GALFIT (Peng et al., 2002, 2010a) and provided by van der
Wel et al. (2014). For these, posterior samples were generated using the
maximum likelihood solutions as the means and $1-\sigma$ errors as Gaussian
standard deviations.

The assumptions used in Prospector include the Chabrier (2003) IMF and a
flexible binned SFH. The dust attenuation law parameterization is detailed in
§3.2 (Equations 11, 12, 13, 14). Table 1 in Leja et al. (2019) lists the
priors on all free parameters of their SED models.

##  3 Methodology

In this section we derive the equations governing our algorithms and explain
how the models work. As there are many mathematical symbols defined throughout
this section, we have created a glossary of terms (Table 1) in Appendix §A.
Before we dive into the details, we define dust attenuation laws below.

The reduction of light intensity through material along a given line of sight
is called extinction. However, the overall effect dust has on light from an
unresolved galaxy is called attenuation. Attenuation is a complex quantity
that is a function of grain composition, dust-star-gas geometry, galaxy
orientation, and the total amount of dust.

One simple way to model dust attenuation is a screen model, in which the
attenuation is modeled with a single dust screen extinguishing the light from
the stars (and gas) of a galaxy. In this scenario, the complex effect of star-
dust geometry is approximated by changes in the dust attenuation law, an
equation for effective extinction vs wavelength. Calzetti et al. (1994, 2000)
used such a method to find that dust attenuation curves of local starbursts
have a relatively homogeneous form. As mentioned in §1, the Calzetti Law
(2000) has played an integral role in dust attenuation studies over the last
two decades.

As configured above, the screen model produces an effective dust attenuation
curve. However, the success of the Charlot & Fall (2000) two-component model
suggests the utility of separating dust into diffuse dust and birth cloud
dust. The two-component dust model is treated as a combination of two screens,
one that affects all stars (diffuse) and one that affects only stars under 10
million years old (birth cloud). We use both effective and two-component dust
screen models in our work. The details are provided in §3.2.

###  3.1 Probability Theory for Population Model

We introduce the mathematical basis for the Bayesian model that we use to
study dust attenuation as a function of various other properties, including
stellar mass, SFR, galaxy inclination, metallicity, and redshift. Throughout
this work, $p(\ldots)$ refers to the probability density of the quantity in
the parentheses, possibly conditional on quantities behind a “$|$”.

Our presentation largely follows that of Foreman-Mackey et al. (2014). We let
$\\{\mathbf{x}_{k}\\}$ represent the entire set of data (mainly photometry)
from the 3D-HST database used in the Prospector fits. Next,
$\\{\mathbf{w}_{k}\\}$ is the set of all the parameters of the galaxies
(including stellar mass, SFR, $n$, etc.). More precisely,
$\\{\mathbf{w}_{k}\\}$ represents the unknown “true” values of these
parameters. In this notation, braces signify a set of galaxies, with the index
($k$) referring to individual galaxies. Finally, bold-face indicates a vector
of parameters.

In reality, not all vector components of $\\{\mathbf{w}_{k}\\}$ factor into
every part of the computation process. For example, our population model does
not consider Active Galactic Nuclei (AGN) fraction, though it is still a
parameter in Prospector. To guarantee we are always considering the right
subset of parameters, we can break down $\\{\mathbf{w}_{k}\\}$ into four
components.

  * •

$\\{\mathbf{t}_{k}\\}$ represents components like stellar mass that we include
in our population Bayesian model that are not directly fitted in Prospector
but are rather derived using other parameters.

  * •

Next, we consider parameters that factor into our model and are directly
fitted in Prospector (such as dust index $n$), which we will call
$\\{\mathbf{u}_{k}\\}$.

  * •

The last two components include parameters that are fitted in Prospector but
are not directly considered in our model (such as the star formation history
parameters). Those that go into calculating $\\{\mathbf{t}_{k}\\}$ we call
$\\{\mathbf{v}_{k}\\}$.

  * •

Finally, we christen the “nuisance” parameters in Prospector that do not have
any bearing toward the population model $\\{\mathbf{y}_{k}\\}$. This includes
parameters such as AGN fraction.

For the majority of this section, we will be focusing on probability
distributions for individual galaxies. To keep the notation as uncluttered as
possible, we will not include the subscript $k$, but it is still implied as
the equations are valid for any given galaxy.

We must remember that $\mathbf{t}$ is a function of $\mathbf{u}$ and
$\mathbf{v}$ since they are simply derived parameters. In other words,
$\mathbf{t}=f(\mathbf{u},\mathbf{v})$, where $f$ is some deterministic
function. Therefore, we can represent $\mathbf{w}$ as
$\left(\mathbf{u},\mathbf{v},\mathbf{y}\right)$.

We let $\boldsymbol{\theta}$ be the parameters of our population model.
Finally, we let $\boldsymbol{\alpha}$ represent the interim priors on
$\mathbf{u},\mathbf{v},\mathbf{y}$ Prospector uses:
$p(\mathbf{u},\mathbf{v},\mathbf{y}|\boldsymbol{\alpha})$.

As mentioned in §2, we are using the posterior samples from Prospector for
nearly 30,000 galaxies to inform our population model. We can use the notation
defined above to specify the distribution from which these samples are taken.

$(\mathbf{u},\mathbf{v},\mathbf{y})\sim
p(\mathbf{u},\mathbf{v},\mathbf{y}|\mathbf{x},\boldsymbol{\alpha})=\frac{p(\mathbf{x}|\mathbf{u},\mathbf{v},\mathbf{y})p(\mathbf{u},\mathbf{v},\mathbf{y}|\boldsymbol{\alpha})}{p(\mathbf{x}|\boldsymbol{\alpha})}$
(1)

Our ultimate task is to derive the posterior distribution of
$\boldsymbol{\theta}$, the parameters determining dust attenuation as a
function of the physical properties for the full population of galaxies. In
order to do that, we must determine a method to compute the probability
density of the entire set of observations used in the Prospector analyses
$\\{\mathbf{x}_{k}\\}$ given any $\boldsymbol{\theta}$, from which we can then
apply Bayes’ Theorem to compute the posterior density
$p({\boldsymbol{\theta}}|\\{\mathbf{x}_{k}\\})$.

Through the reasonable assumption that (the SEDs of) individual galaxies are
independent of one another, we can write
$p(\\{\mathbf{x}_{k}\\}|\boldsymbol{\theta})=\prod_{k}p(\mathbf{x}_{k}|\boldsymbol{\theta})$
as long as the normalization is done properly. We begin with a statement of
the law of total probability (once again omitting the subscript $k$).

$p(\mathbf{x}|\boldsymbol{\theta})=\int
p(\mathbf{x}|\mathbf{u},\mathbf{v},\mathbf{y})p(\mathbf{u},\mathbf{v},\mathbf{y}|\boldsymbol{\theta})d\mathbf{u}d\mathbf{v}d\mathbf{y}$
(2)

In Equations 1 and 2, we have made the reasonable assumption that the
observations for a galaxy $\mathbf{x}$ do not depend directly on the interim
prior parameters $\boldsymbol{\alpha}$ or the distribution of physical
parameters $\boldsymbol{\theta}$, respectively, but rather just on the set of
true values of the parameters $\mathbf{u},\mathbf{v},\mathbf{y}$. In fact,
there may be shared systematic effects and modeling choices that could break
the independence assumptions we have made, but we do not treat those effects
here.

Equation 2 presents a high-dimensional integral that is difficult to compute
as is, especially when considering that the calculation would need to be
repeated for all galaxies in the sample. We will make approximations and
adjustments to alleviate our task. First, we will multiply the integrand in
Equation 2 by 1 in the following fashion.

$p(\mathbf{x}|\boldsymbol{\theta})=\\\ \int
p(\mathbf{x}|\mathbf{u},\mathbf{v},\mathbf{y})p(\mathbf{u},\mathbf{v},\mathbf{y}|\boldsymbol{\theta})\frac{p(\mathbf{u},\mathbf{v},\mathbf{y}|\mathbf{x},\boldsymbol{\alpha})}{p(\mathbf{u},\mathbf{v},\mathbf{y}|\mathbf{x},\boldsymbol{\alpha})}d\mathbf{u}d\mathbf{v}d\mathbf{y}$
(3)

We can use Equation 1 to simplify Equation 3. We also use the fact that
$p(\mathbf{x}|\boldsymbol{\alpha})$ is independent of the variables of
integration.

$\frac{p(\mathbf{x}|\boldsymbol{\theta})}{p(\mathbf{x}|\boldsymbol{\alpha})}=\int\frac{p(\mathbf{u},\mathbf{v},\mathbf{y}|\boldsymbol{\theta})}{p(\mathbf{u},\mathbf{v},\mathbf{y}|\boldsymbol{\alpha})}p(\mathbf{u},\mathbf{v},\mathbf{y}|\mathbf{x},\boldsymbol{\alpha})d\mathbf{u}d\mathbf{v}d\mathbf{y}$
(4)

However, as we explained while defining terms, the population distribution we
want to fit is actually $p(\mathbf{u},\mathbf{t}|\boldsymbol{\theta})$ rather
than $p(\mathbf{u},\mathbf{v},\mathbf{y}|\boldsymbol{\theta})$. As defined
earlier, $\mathbf{y}$ consists of variables such as AGN fraction that are not
involved in our population model. In Prospector, these parameters have priors
that are independent of the priors on $\mathbf{u}$ and $\mathbf{v}$.
Furthermore, as $\mathbf{y}$ is unrelated to $\boldsymbol{\theta}$, the priors
on our population model do not provide any additional constraints on it
compared to the interim prior:
$p(\mathbf{y}|\boldsymbol{\theta})=p(\mathbf{y}|\boldsymbol{\alpha})$
Therefore, we have the following.

$\frac{p(\mathbf{u},\mathbf{v},\mathbf{y}|\boldsymbol{\theta})}{p(\mathbf{u},\mathbf{v},\mathbf{y}|\boldsymbol{\alpha})}=\frac{p(\mathbf{u},\mathbf{v}|\boldsymbol{\theta})}{p(\mathbf{u},\mathbf{v}|\boldsymbol{\alpha})}\frac{p(\mathbf{y}|\boldsymbol{\theta})}{p(\mathbf{y}|\boldsymbol{\alpha})}=\frac{p(\mathbf{u},\mathbf{v}|\boldsymbol{\theta})}{p(\mathbf{u},\mathbf{v}|\boldsymbol{\alpha})}$
(5)

In addition, we can write,

$\frac{p(\mathbf{u},\mathbf{v}|\boldsymbol{\theta})}{p(\mathbf{u},\mathbf{v}|\boldsymbol{\alpha})}=\frac{p(\mathbf{u},\mathbf{t}|\boldsymbol{\theta})}{p(\mathbf{u},\mathbf{t}|\boldsymbol{\alpha})}$
(6)

since the function $\mathbf{t}=f(\mathbf{u},\mathbf{v})$ does not depend on
$\boldsymbol{\alpha}$ or $\boldsymbol{\theta}$, implying that the Jacobian of
the transform in the numerator and denominator will be the same and therefore
cancel in this ratio. From these considerations, we can rewrite Equation 4 as
follows.

$\frac{p(\mathbf{x}|\boldsymbol{\theta})}{p(\mathbf{x}|\boldsymbol{\alpha})}=\int\frac{p(\mathbf{u},\mathbf{t}|\boldsymbol{\theta})}{p(\mathbf{u},\mathbf{t}|\boldsymbol{\alpha})}p(\mathbf{u},\mathbf{v},\mathbf{y}|\mathbf{x},\boldsymbol{\alpha})d\mathbf{u}d\mathbf{v}d\mathbf{y}$
(7)

We can numerically estimate $p(\mathbf{u},\mathbf{t}|\boldsymbol{\alpha})$ by
sampling directly from Prospector with no data. We discuss this process in
Appendix B and show examples of interim prior probability distributions.
Meanwhile, from Equation 1, we know that the Prospector posterior samples are
taken from the distribution
$p(\mathbf{u},\mathbf{v},\mathbf{y}|\mathbf{x},\boldsymbol{\alpha})$.
Therefore, we can use the Monte Carlo integral approximation to estimate the
integral in Equation 7.

$\frac{p(\mathbf{x}|\boldsymbol{\theta})}{p(\mathbf{x}|\boldsymbol{\alpha})}\approx\frac{1}{M}\sum_{m=1}^{M}\frac{\left(\mathbf{u}^{(m)},\mathbf{t}^{(m)}|\boldsymbol{\theta}\right)}{\left(\mathbf{u}^{(m)},\mathbf{t}^{(m)}|\boldsymbol{\alpha}\right)}$
(8)

To fit the set of all galaxies, based on the galaxy independence assumption
stated earlier, we simply take the product of probabilities for all galaxies.

$\frac{p(\\{\mathbf{x}_{k}\\}|\boldsymbol{\theta})}{p(\\{\mathbf{x}_{k}\\}|\boldsymbol{\alpha})}=\prod_{k=1}^{N}\frac{p(\mathbf{x}|\boldsymbol{\theta})}{p(\mathbf{x}|\boldsymbol{\alpha})}$
(9)

The quantity $p(\\{\mathbf{x}_{k}\\}|\boldsymbol{\theta})$ is called the
likelihood, which we will label $\mathcal{L}$. For computational purposes, we
take the logarithm and ignore the constant
$p(\\{\mathbf{x}_{k}\\}|\boldsymbol{\alpha})$ since it is independent of
$\boldsymbol{\theta}$.

$\ln{\mathcal{L}}=\sum_{k=1}^{N}\ln\left(\sum_{m=1}^{M}\frac{p\left(\mathbf{t}_{k}^{(m)},\mathbf{u}_{k}^{(m)}|\boldsymbol{\theta}\right)}{p\left(\mathbf{t}_{k}^{(m)},\mathbf{u}_{k}^{(m)}|\boldsymbol{\alpha}\right)}\right)$
(10)

In §3.4, we make Equation 10 more specific to our model.

###  3.2 Dust Attenuation Parameterization

As mentioned in 3.1, the goal of this work is to learn how the dust
attenuation parameters vary as a function of other galaxy properties,
including stellar mass, SFR, stellar metallicity, axis ratio, and redshift.
Before getting into the details of our population model, we will describe the
parameterization of dust attenuation that we are using. Prospector uses a two-
component model (modified from Charlot & Fall 2000): diffuse dust, affecting
all sources in each galaxy, and birth cloud dust, affecting just the most
recently formed stars. The birth cloud dust attenuation curve is modeled as a
$\lambda^{-1}$ power law while the diffuse dust is modeled in the form
introduced by Noll et al. (2009).

$\tau_{\lambda,1}=\tau_{1}\left(\frac{\lambda}{5500\AA}\right)^{-1}$ (11)

$\tau_{\lambda,2}=\frac{\tau_{2}}{4.05}\left(k^{\prime}(\lambda)+D(\lambda)\right)\left(\frac{\lambda}{5500\AA}\right)^{n}$
(12)

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

$\lambda_{0},~{}\Delta\lambda,~{}E_{b}=2175\AA,~{}350\AA,~{}0.85-1.9n$ (14)

Here, $k^{\prime}(\lambda)$ is the Calzetti Law (see Calzetti et al., 2000).
The formula $E_{b}=0.85-1.9n$ is given by Kriek & Conroy (2013). The
parameters that are allowed to vary are $\tau_{1}$, $\tau_{2}$, and $n$.
However, Prospector puts a normal prior on $\tau_{1}/\tau_{2}$ with a mean
$\mu=1$ and width $\sigma=0.3$ based on empirical evidence (e.g., Calzetti et
al., 1994; Price et al., 2014). This prior places strong constraints on
$\tau_{1}$. Therefore, we focus on $\tau_{2}$ (diffuse dust normalization) and
$n$ (relating to the slope of the dust attenuation curve, with negative values
indicating steeper than the Calzetti Law and positive values indicating
shallower) as the primary dependent variables of the dust attenuation model.

In addition, we compute dust attenuation models for a one-component dust model
with the same parameterization as diffuse dust (Equation 12). To ensure the
same treatment in the population Bayesian modeling, we calculated the
effective dust attenuation curve for every Prospector posterior sample by
generating the dust-free and dust-included spectra using the Prospector
parameter values and fitting an optical depth model with Equation 12 using a
basic curve fitting package from SciPy. The one-component model is useful in
its simplicity, especially in cases where not much is known about the stellar
populations comprising a given galaxy.

While we provide different models for single- and two-component dust models,
it is useful and interesting to study the direct relationship between the two
as we can learn more about how the stellar and nebular components contribute
individually to the overall picture. In Figure 2, we provide the comparison
between diffuse and effective dust.

The left panel shows $n_{\rm eff}$ vs $n$ (diffuse). We see that the effective
slope is systematically steeper than the diffuse dust slope, with a difference
of $0.05-0.2$, which is consistent with the theoretical expectations and
observational findings of Inoue (2005),Salim et al. (2018), and Salim &
Narayanan (2020), as stars affected by both diffuse and birth cloud dust will
be heavily attenuated.

However, we do observe objects above the one-to-one line in Figure 2. As we
did not impose the requirement that $n_{\rm eff}<n_{\rm diff}$, there are
cases in which the fitted effective optical depth is considerably larger than
the diffuse optical depth and the fitted effective slope is shallower than the
diffuse slope. Nevertheless, such cases account for fewer than 1% of all
posterior samples, making them have minimal effect on the conclusions we draw
throughout the paper.

The fact that the best-fit slope for $n_{\rm eff}$ vs $n_{\rm diff}$ is not
one is a bit more complex. As discussed later in §5.2, higher optical depths
are linked to shallower slopes (Salim & Narayanan, 2020, and references
therein). In addition, as we show in Figure 3 (§4), birth cloud dust optical
depth increases super-linearly with diffuse dust optical depth. Therefore, at
lower optical depths, which are associated with steeper slopes, the difference
between birth cloud dust and diffuse dust is smaller, leading to a smaller
change between diffuse and effective slope.

The right panel shows $\tau_{\rm eff}$ vs $\tau_{2}$. The relationship between
these two is straightforward: as long as the effective slope is not radically
different than the diffuse slope, the effective dust optical depth must be
greater than the diffuse dust optical depth. As optical depth increases, the
birth cloud contribution increases, leading to a greater divergence between
effective and diffuse.

Figure 2: Direct comparison of effective and diffuse dust attenuation curves.
Left panel: slope parameter comparison. In agreement with Inoue (2005), Salim
et al. (2018), and Salim & Narayanan (2020), we find that the effective slope
is systematically steeper than the diffuse slope, which can be attributed to
the unevenly distributed effects of birth cloud dust on UV light from a (star-
forming) galaxy. Right panel: optical depth comparison. As expected, the
effective optical depth is greater than the diffuse optical depth due to the
contribution from birth cloud dust.

###  3.3 Interpolation Configuration

The functional form we have chosen for the dust attenuation curve is an N-D
linear interpolation with an intrinsic Gaussian scatter. Given that the slope
and normalization of the dust attenuation curve have been found to be
correlated (Salim & Narayanan, 2020, and references therein), we consider the
joint modeling of the slope ($n$) and optical depth ($\tau$)111Henceforth,
when generalizing for both one-component and two-component dust models, $n$
and $\tau$ without any subscripts will be used. For the one-component model
specifically, we will use $n_{\rm eff}$ and $\tau_{\rm eff}$, and for two-
component $n$ and $\tau_{2}$. While the interpolation components of the model
are considered independent, we treat the intrinsic scatter as a bivariate
Gaussian with a possible covariance $\rho$ between $n$ and $\tau$.

Below, we show the equations for our population model, where $I_{1}$ and
$I_{2}$ are linear N-D interpolations, $\mathcal{Z}_{1}$ and $\mathcal{Z}_{2}$
are samples from the standard normal distribution, and $\sigma_{1}$ and
$\sigma_{2}$ are the intrinsic widths of the relation in $n$ and $\tau_{2}$,
respectively.

$\displaystyle n$ $\displaystyle=I_{1}(\log~{}M_{*},\log{\rm
SFR},...)+\sigma_{1}\mathcal{Z}_{1}$ (15) $\displaystyle\tau$
$\displaystyle=I_{2}(\log~{}M_{*},\log{\rm
SFR},...)+\sigma_{2}\left(\rho\mathcal{Z}_{1}+\sqrt{1-\rho^{2}}\mathcal{Z}_{2}\right)$
(16)

In addition to this bivariate model, we compute univariate models in which
$\tau$ is an independent variable and $n$ is the sole dependent variable
(Equation 15).

The parameters for the interpolation $I_{1}$ and $I_{2}$ are the values of $n$
and $\tau$, respectively, at a given set of grid points in the N-D independent
parameter space. We employ unevenly spaced grids where points are closer
together in regions where more data is present. This process helps maximize
the utility of the data. We use five grid points per dimension, leading to
3125 total grid points for $n$ and $\tau$ (6250 parameters total) for the 5-D
models.

###  3.4 Modeling Details

In §3.1, we laid out the mathematical foundation for our population Bayesian
model. In §3.3, we described the model parameterization we adopted. In this
section, we connect §3.1 and §3.3 and provide details about the modeling
process itself.

To run our population model, we use the Python package pymc3 (Salvatier et
al., 2016), which employs an efficient Hamiltonian Monte Carlo algorithm
called No U-Turn Sampler (NUTS; Hoffman & Gelman, 2011) to move through the
parameter space. We typically utilize four Markov chains, each with 1,000 tune
steps (initial steps that are subsequently thrown away) and 1,000 accepted
steps. In all cases for our simulated data (Appendix C) and most cases for our
actual results (§4), these numbers are sufficient for chain convergence.

The model parameters are the values of $n$ and $\tau$ at the interpolation
grid points (6250 parameters for the 5-D models), $\log\sigma_{1}$,
$\log\sigma_{2}$, and $\rho$. The prior for $n$ and $\tau$ grid points is
simply the Prospector prior for those variables: a uniform distribution
between $[-1.0,0.4]$ and a truncated Gaussian distribution with
$(\mu=0.3,~{}\sigma=1.0)$ and bounds $[0.0,4.0]$, respectively (with slightly
different values for the one-component dust model). The priors for
$\log\sigma_{1}$ and $\log\sigma_{2}$ are Student-T distributions with
$(\nu=5,\mu=-3.0,\sigma=0.5)$ while the prior for $\rho$ is a uniform
distribution in $[-1,1]$. For simplicity, we will refer to the $n$ and $\tau$
grid parameters using the interpolation functions $I_{1}$ and $I_{2}$.

If we let $I_{1}^{*}$, $I_{2}^{*}$, $\sigma_{1}^{*}$, $\sigma_{2}^{*}$, and
$\rho^{*}$ represent the current model parameters, where ∗ denotes the
particular model, we can define the likelihood of the $m^{\rm th}$ posterior
of the $k^{\rm th}$ galaxy in the sample in the following manner. Here,
$I_{1,k,m}^{*}\equiv I_{1}^{*}(\log M_{*k,m},...)$, etc.

$\displaystyle\mathcal{Z}_{k,m}$
$\displaystyle=\frac{{\left(n_{k,m}-I_{1,k,m}^{*}\right)}^{2}}{\sigma_{1}^{2}}+\frac{\left(\tau_{2,k,m}-I_{2,k,m}^{*}\right)^{2}}{\sigma_{2}^{2}}$
$\displaystyle+\frac{\left(n_{k,m}-I_{1,k,m}^{*}\right)\left(\tau_{2,k,m}-I_{2,k,m}^{*}\right)}{\sigma_{1}\sigma_{2}}$
(17) $\displaystyle\ln\mathcal{L}_{{\rm hyper},k,m}$
$\displaystyle=-0.5\frac{\mathcal{Z}_{k,m}}{1-\rho^{2}}-\ln{\left(\sigma_{1}\sigma_{2}\sqrt{1-\rho^{2}}\right)}$
(18)

From Equation 10, we see that to compute the true likelihood, we still need to
account for the Prospector priors. As the process of computing the prior
probabilities in a consistent manner with the likelihood in Equation 18 is not
straightforward but is a minor concern compared to the details of the
population model, we have kept details of our process in Appendix B.

Based on Equation 10, we compute the contribution to the likelihood from the
$k$th galaxy’s $m$th Prospector sample

$\ln\mathcal{L}_{k,m}=\ln\mathcal{L}_{{\rm hyper},k,m}-\ln\pi_{k,m},$ (19)

where $\pi_{k,m}$ is the Prospector prior on the subset of variables used in
our population model evaluated at the given Prospector sample, indexed by
$k,m$, as described in Appendix B.

We have one more consideration for the overall likelihood calculation that is
particular to our model formulation. While interpolations provide a flexible
and essentially non-parametric form for the models, they often struggle with
smoothness, especially when the data has strong variability or noise. This
last situation applies to the Prospector data, so we add a constraint to help
remove unnaturally large variations from the model.

For the constraint, we stipulate that the differences between the dependent
variables at adjacent grid points for the interpolation follow Gaussian
distributions with mean $\mu=0$ and a fixed standard deviation $\sigma$. While
the value of $\sigma$ could in principle be inferred as an additional model
parameter, in practice this presents a challenge to the NUTS sampler. For
simplicity we have instead tried several fixed values, and we find that for
our 5-D models, $\sigma=1.0$ allows for a good balance between smoothness and
detail. Note that the dimensions of $\sigma$ are $n$ or $\tau$ divided by the
(typically logarithmic and adaptive) spacing between adjacent grid points. We
sum over the log probability densities of every difference to get the total
“pair potential,” which we denote $\gamma_{\mathrm{pair}}$. This quantity is
not galaxy-dependent and is instead a constant for a particular choice of
$\boldsymbol{\theta}$.

The total likelihood $\ln\mathcal{L}_{\rm tot}$ can now be calculated with the
following equations. Here, $M$ is once again the total number of posterior
samples per galaxy, $N$ is the total number of galaxies, and
$\ln\mathcal{L}_{k,{\rm max}}$ is the maximum likelihood posterior sample for
a given galaxy. This last quantity is used to prevent extreme numbers in the
exponentiation process.

$\displaystyle\ln\mathcal{L}_{k,{\rm marg}}$
$\displaystyle=\ln\mathcal{L}_{k,{\rm
max}}+\ln\left[\sum_{i=1}^{M}\exp{\left(\ln{\mathcal{L}_{k,i}}-\ln\mathcal{L}_{k,{\rm
max}}\right)}\right]$ (20) $\displaystyle\ln\mathcal{L}_{\rm tot}$
$\displaystyle=\gamma_{\mathrm{pair}}+\sum_{k=1}^{N}\mathcal{L}_{k,{\rm
marg}}$ (21)

###  3.5 Summary of Models

Using the posterior samples from Prospector of nearly $\sim 30,000$ 3D-HST
galaxies at $0.5<z<3.0$ (Leja et al., 2019, 2020) and the hierarchical
Bayesian framework described throughout this section, we have modeled dust
attenuation curves as a function of stellar mass, SFR, metallicity, redshift
and axis ratio. In Appendix C, we demonstrate how successfully our framework
works by applying it to a case in which the true relation is known.

Given the connection between dust attenuation slope and normalization (Salim &
Narayanan, 2020, and references therein), we model the slope ($n$) and optical
depth ($\tau$) of the attenuation curves jointly, with a bivariate Gaussian
scatter. Going forward, we refer to such models as “bivariate.”

We also calculate models with only $n$ as the dependent variable and $\tau$ as
an independent variable. Such cases employ a univariate Gaussian error
distribution rather than the more complex bivariate distribution described in
§3.3. We refer to such models as “univariate.” We find that in general the
best model for $n$ includes $\tau$ as an independent variable, as expected
from Salim & Narayanan (2020).

Our models come in two flavors: single- and two-component dust models. The
latter models follow directly from the dust treatment of Prospector. Diffuse
dust, parameterized by Equation 12, affects all stars in the galaxy, whereas
birth cloud dust, parameterized by Equation 11, affects only stars under 10
Myr old. When the star formation history of a galaxy (whether observed or
modeled) is not well defined or known, it is useful to have a one-component
effective dust model. Once again, we use the flexible parameterization of
Equation 12.

Therefore, we have a total of four 5-D models:

  1. 1.

Univariate (dependent variable $n$), two-component dust model

  2. 2.

Univariate (dependent variable $n_{\rm eff}$), one-component dust model

  3. 3.

Bivariate (dependent variables $n$ and $\tau$), two-component dust model

  4. 4.

Bivariate (dependent variables $n_{\rm eff}$ and $\tau_{\rm eff}$), one-
component dust model

All of these models are functions of log stellar mass, log SFR, log stellar
metallicity, redshift, and a fifth parameter. For the univariate cases, the
fifth parameter is the dust optical depth $\tau$, whereas for the bivariate
cases the parameter is the axis ratio $b/a$, which is a proxy for inclination.

In the case of the two-component dust models, we model only the diffuse dust
with this 5D configuration. Given the close relation between diffuse and birth
cloud dust (e.g., Calzetti et al., 1994; Price et al., 2014), which is
enforced in Prospector priors, as well as the difficulty to distinguish its
relatively small effects on unresolved galaxy spectra, we model birth cloud
dust (optical depth) as a function solely of diffuse dust (optical depth).
Though this is a simpler model with only one independent variable, it is still
a proper population model with the setup described in this section.

##  4 Results

In this section, we explore the relationships among dust optical depth,
attenuation slope, stellar mass, SFR, metallicity, redshift, and axis ratio
based on population Bayesian models. As detailed in §3.5, we have four 5-D
linear interpolation models for diffuse and effective dust as well as one 1-D
model for birth cloud dust.

In §4.1, we highlight the quasi-linear relation between diffuse dust and birth
cloud dust optical depth. In §4.2, we show that our population model produces
different but more accurate trends than a simple model that ignores
uncertainties in individual galaxy parameter values. Then, in §4.3, we compare
our derived attenuation slope vs $A_{V}$ relation to the literature. Next, we
present the relation among $\tau_{\rm eff}$, stellar mass, and SFR in §4.4.
Subsequently, in §4.5, we study the effects of inclination on dust
attenuation. We then exhibit the cosmic evolution of attenuation in §4.6.
Finally, we examine the intricacies of dust attenuation slope in §4.7.

###  4.1 Birth Cloud Dust

In Figure 3, we show our model for diffuse dust optical depth $\tau_{2}$ vs
birth cloud dust optical depth $\tau_{1}$. Note that although uncertainties on
the mean relation are shown on the plot, they are too small to be visible. We
also show the $1-\sigma$ Prospector prior region in $\tau_{2}$ vs $\tau_{1}$
space in orange.

We can see that the relationship is close to linear, but the values are not
consistent with a slope of 1 (except when $\tau_{2}\lesssim 1$). Birth cloud
dust tends to have a higher normalization (optical depth) than diffuse dust.
Moreover, as $\tau_{1}$ increases, the model curve moves farther from the
Prospector prior, showing that the result in Figure 3 is not merely a
restatement of the Prospector prior but rather one that is informed by
photometric data, and indeed qualitatively similar to results using the same
data (Price et al., 2014). We discuss the physical context and implications of
the birth cloud – diffuse dust connection in §5.1.

Figure 3: Result of modeling diffuse dust optical depth $\tau_{2}$ solely as a
function of birth cloud optical depth $\tau_{1}$. The model is nearly linear
but does not have a slope of 1 (except when $\tau_{2}\lesssim 1$). Note that
uncertainties on the mean relation are shown but too small to be visible. We
also show the Prospector $1-\sigma$ prior region in $\tau_{2}$ vs $\tau_{1}$
space in orange. As $\tau_{1}$ increases, the model moves farther from the
prior, highlighting the influence of the photometry on the result.

###  4.2 Benefits of the Population Approach

We find that the more sophisticated population modeling approach with a proper
error treatment gives a qualitatively and quantitatively different answer than
a simpler one in which parameters of individual galaxies are treated as
singular points. This corroborates the conclusion from Qin et al. (2022) that
degeneracies between $n$, $A_{V}$, and the intrinsic UV slope $\beta_{0}$ from
SED fitting strongly influence the resulting slope vs $A_{V}$ relation.

In Figure 4, We show the diffuse dust slope parameter $n$ over 2-D cross-
section of $\tau_{2}$ vs log stellar mass in three different situations. In
the left panel, we take the maximum likelihood Prospector fit for each galaxy,
bin the space into a $25\times 25$ grid, and record the median $n$ value in
each bin as long as at least 10 galaxies are within the bin. This is by far
the simplest and most direct way of exploring the relationships between
galaxies, but if the posterior volumes of individual galaxies are large (i.e.,
significant measurement errors), which is the case here, galaxies are likely
to end up in the wrong bins. Furthermore, if posteriors are correlated, the
population trends may reflect those degeneracies rather than any true trends.
In addition, the maximum likelihood solution is simply a particularly
interesting moment of the full parameter posterior; using the posterior
directly produces more reliable results.

In terms of identifiable features, we see that the relationship between $n$,
stellar mass, and $\tau_{2}$ is very pronounced, with the shallowest
attenuation curves at high $\tau_{2}$ and
$\log\left(M_{*}/M_{\odot}\right)\sim 9.5$.

In the middle panel, we show the results of a simple Bayesian model in which
we continue to ignore uncertainties in the galaxy properties, just limiting
ourselves to the maximum likelihood Prospector fit, which simplifies the
likelihood calculations. We marginalize (see Appendix D) over the other three
parameters (log SFR, log metallicity, and redshift).

We can see the strongest trend from the left panel reproduced in our model:
the shallowest slopes are found at higher $\tau_{2}$ values and stellar masses
around $\log\left(M_{*}/M_{\odot}\right)\sim 9.8$ whereas the attenuation
curves at $\log\left(M_{*}/M_{\odot}\right)\sim 11$ and low $\tau_{2}$ are
quite steep.

In the right panel, we show the mean results (not including the intrinsic
scatter) of our univariate model for diffuse dust. The marginalization is
performed in the same way as for the simple model. We find largely the same
results as in the simple model, but with a slightly more pronounced shallow-
slope peak for $\log~{}M_{*}\sim 10$ and $\tau_{2}\gtrsim 0.5$.

Before we analyze Figure 4, we must mention a couple of caveats. First, the
left panel is not an indication of the truth as it is subject to the errors in
measurements of individual galaxy parameters, though the median binning method
is robust to some non-systematic sources of error. Second, the degree of
similarity between the binned parameters and the models depends on the
particular cross-section. For the $\tau$ vs stellar mass space (Figure 4), the
population model shows a stronger trend, slightly more in accordance with the
binned parameters, but for a few other cross-sections, the simple model shows
stronger trends.

Nevertheless, what we observe from Figure 4 is that the simple model and
population model yield qualitatively and quantitatively different results.
This shows that properly accounting for covariances in the galaxy parameters
affects the result. In other words, employing a population or hierarchical
model will produce different and more accurate results.

The fact that the smaller-scale and more minor trends in the left panel are
not reproduced by the models is not surprising. Our smoothing operation (see
§3.4) allows our interpolation models to converge and avoid extreme variations
on small scales due to noisy data. However, it does also prevent the recovery
of small scale or extremely minor trends. Of course, such trends as found
directly from data (i.e., binning) are quite suspect and are often not
statistically significant.

One thing not conveyed through Figure 4 is the relative amounts of uncertainty
in the simple vs population models. We find that the simple Bayesian model is
unable to recover as much of the true trends and indeed predicts a larger
intrinsic scatter ($\log\sigma=-1.21$ compared to $\log\sigma=-1.72$ in the
population model). In other words, our population model is more successful at
attributing causation to the input parameters than the simpler model.

Figure 4: Prediction of diffuse dust attenuation slope parameter $n$ vs
stellar mass and $\tau_{2}$ using three methods. Left panel: we bin the
$\tau_{2}$ vs $\log~{}M_{*}$ parameter space into $25\times 25$ and record the
median $n$ value from Prospector within the bin as long as at least 10
galaxies are in the bin. Maximum likelihood values for each parameter are
assumed. Middle panel: we show the result of a simple Bayesian model in which
we take the maximum likelihood solutions for each parameter and ignore
measurement uncertainty. Right panel: result of our univariate diffuse dust
model. The red points in the middle and right panels show the interpolation
grid in the $\tau_{2}$ vs $\log~{}M_{*}$ space. From all three panels, we
observe that the relation among $n$, $\tau_{2}$, and $\log~{}M_{*}$ is strong.
We find a qualitatively and quantitatively different result in the population
model (right) than the simple Bayesian model (middle), highlighting that
properly accounting for uncertainties in individual galaxy parameters
substantially affects the results. While it is not shown here, the population
model is able to explain more of the scatter in the data than the simpler
model.

###  4.3 Relationship between Slope and $A_{V}$

In Figure 5, we show the relation between effective slope and $A_{V}$ for our
one-component univariate dust model. The quantity $S$ is defined as
$S=A_{1500}/A_{V}$. For convenience, we show the parameter $n$ used in our
models on the right y-axis.

We compare our results to Buat et al. (2012), Kriek & Conroy (2013), Reddy et
al. (2015), Battisti et al. (2019), Álvarez-Márquez et al. (2019), as
calculated and plotted by Salim & Narayanan (2020), in Figure 5. We also
include the results from Salim et al. (2018), as plotted by Salim & Narayanan
(2020), for the local universe, and plot the Calzetti et al. (2000) relation
as a baseline. For reference, the average extinction curves for the Milky Way
and Small Magellanic Cloud are included.

In the univariate model, $n_{\rm eff}$ is a direct function of $\tau_{\rm
eff}$. We vary $\tau_{\rm eff}$, or more accurately $A_{V}\approx
1.086\tau_{\rm eff}$, on a grid from $0$ to roughly $1.4$, and marginalize
over the other four independent variables.

We randomly sample from the model posterior and apply each instance (a single
posterior sample per galaxy) to all “galaxies” to create a curve per sample
set. We take the average curve (black line) and calculate the standard
deviation of the average curve (dark shaded region) and of all realizations of
the model posterior in general (light shaded region).

Our average curve tends to have a less pronounced (shallower) relation between
slope and optical depth than most previous results, with the notable
exceptions of Buat et al. (2012) and Reddy et al. (2015). We discuss the
implications of our results in §5.2.

Figure 5: Dust attenuation slope vs $A_{V}$ as predicted by our one-component
univariate population model in comparison to previous studies: Calzetti et al.
(2000), Buat et al. (2012), Kriek & Conroy (2013), Reddy et al. (2015), Salim
et al. (2018), Battisti et al. (2019), Álvarez-Márquez et al. (2019). The
literature values are taken from Figure 10b in Salim & Narayanan (2020).
Average extinction-law values for the Milky Way and Small Magellanic Cloud are
included. We show the best-fit curve as a solid black line with the mean
standard deviation shaded darkly and the sample standard deviation shaded
lightly. We see that our model predicts a flatter relation between slope and
$A_{V}$ compared to the literature, with the exceptions of Buat et al. (2012)
and Reddy et al. (2015). We provide a possible explanation in §5.2.

###  4.4 Dust Optical Depth vs Stellar Mass and SFR

We find that the dust effective (and diffuse) optical depth $\tau$ depends
most strongly on SFR. The dependence on the other four parameters is not as
significant but still substantial. In Figure 6, we show the relation between
$\tau_{\rm eff}$, SFR, and stellar mass.

The center panel shows the 2-D cross-section of the model marginalized over
the other three parameters. The behavior of the model is quite complex, making
it quite clear that the dependence of $\tau$ on stellar mass and SFR cannot be
separated. In other words, dust optical depths in massive galaxies are
generally higher at a fixed SFR than in low-mass galaxies.

The left and right panels show the 1-D cross-sections of the model at five
fixed values of the other variable. The left panel shows $\tau_{\rm eff}$ vs
log stellar mass. We observe that at the three lowest values of SFR, the
relation in stellar mass is nearly flat, suggesting little to no correlation
between $\tau_{\rm eff}$ and stellar mass in galaxies with low star formation
rates. On the other hand, $\tau_{\rm eff}$ clearly increases with stellar mass
for star-forming galaxies for $\log~{}M_{*}\lesssim 10.5$. For
$\log~{}M_{*}\gtrsim 10.5$, $\tau_{\rm eff}$ seems to slowly decrease with
stellar mass for most values of SFR.

Meanwhile, the right panel shows $\tau_{\rm eff}$ vs log SFR. Interestingly,
for all five stellar masses, $\tau_{\rm eff}$ stays nearly constant until
$\log~{}{\rm SFR}\sim 0$, above which it quickly increases.

Figure 6: Effective dust optical depth $\tau_{\rm eff}$ as a function of
stellar mass and SFR, marginalized over $\log~{}(Z/Z_{\odot})$, $z$, and
$b/a$. The center panel shows the full 2-D cross-section while the left and
right panels show the 1-D cross sections at five values of the other variable.
The red points represent the interpolation grid while the faint black pixels
show the Prospector mass-complete sample. We observe that $\tau_{\rm eff}$ is
nearly constant at low SFRs, while it increases with stellar mass in
vigorously star forming galaxies until $\log~{}M_{*}\sim 10.5$, above which it
decreases slightly. $\tau_{\rm eff}$ increases with SFR at all stellar masses
when $\log~{}{\rm SFR}\gtrsim 0$.

###  4.5 The Effects of Inclination on Dust Attenuation

In our bivariate model, we include the axis ratio $b/a$ as a proxy for
inclination. As mentioned in §2, the axis ratio data comes from van der Wel et
al. (2014) through GALFIT (Peng et al., 2002, 2010a) fits.

In our model, $b/a$ is roughly as important as stellar mass in determining the
dust optical depth. Figure 7 shows $\tau_{\rm eff}$ over the 2-D cross-section
of $b/a$ vs $\log~{}M_{*}$ (center) in addition to the 1-D cross-sections of
$b/a$ (right) and $\log~{}M_{*}$ (left).

In the $\tau_{\rm eff}$ vs $\log~{}M_{*}$ plot, we can see that for
$\log~{}M_{*}\lesssim 10$, the curves at all values of $b/a$ are overlapping.
As mass increases, the curves spread out, and we can see that edge-on galaxies
(represented by the $b/a=0.09$ curve) have higher $\tau_{\rm eff}$ values than
the face-on galaxies (represented by the $b/a=0.97$ curve).

Similarly, in the $\tau_{\rm eff}$ vs $b/a$ plot, we see that for the highest
two $\log~{}M_{*}$ values, $\tau_{\rm eff}$ decreases as we move toward face-
on galaxies. On the other hand, the curves for the lowest three $\log~{}M_{*}$
values are closer to constant.

Figure 7: Effective dust optical depth $\tau_{\rm eff}$ vs $\log~{}M_{*}$
(left) and axis ratio $b/a$ (right), marginalized over $\log$ SFR,
$\log~{}(Z/Z_{\odot})$, and $z$. The center panel shows the full relationship
while the two side plots show the 1-D cross-sections for five fixed values of
the other parameter. We find that edge-on galaxies tend to have larger optical
depths than face-on galaxies, but only for galaxies with $\log~{}M_{*}\gtrsim
10$. Low mass galaxies in the distant universe have still not developed a
distinct disk-like morphology, resulting in no relationship between $b/a$ and
viewing $\tau_{\rm eff}$.

To probe the question of inclination farther, we show how our modeled
$\tau_{\rm eff}$ and $n_{\rm eff}$ correlate with $b/a$ as a function of
$\log~{}M_{*}$ and $\log$ SFR in Figure 8. The left panels of the plots show
the $R^{2}$ value of the model $\tau_{\rm eff}$ (top) and $n_{\rm eff}$
(bottom) vs $b/a$ in 225 bins in the $\log$ SFR vs $\log~{}M_{*}$ parameter
space. In other words, we are determining how important inclination is in the
dust attenuation model for different values of stellar mass and SFR. In the
right panels, we show the predicted model for $\tau_{\rm eff}$ (top) and
$n_{\rm eff}$ (bottom) vs $b/a$ at the highest and lowest $R^{2}$ values in
the parameter space.

We see that inclination plays the biggest role in determining dust attenuation
curves for massive, vigorously-star-forming galaxies. For $n_{\rm eff}$, the
peak in $R^{2}$ is at slightly lower SFRs than for $\tau_{\rm eff}$, but the
overall effect on the attenuation curve is as stated. As expected from Figure
7, for massive, star-forming galaxies, $\tau_{\rm eff}$ decreases as we move
toward face-on galaxies. Given the strong relation between $\tau_{\rm eff}$
and $n_{\rm eff}$ (Figure 5), the steepening of the attenuation curve as we
move toward face-on galaxies is a natural consequence.

The only region where inclination plays a minimal role in attenuation is at
low masses and low SFRs, which is consistent with our findings from Figure 7.
In this region, we observe that $\tau_{\rm eff}$ is nearly constant while
$n_{\rm eff}$ has a quasi-sinusoidal relation with $b/a$. The net effect on
the overall dust attenuation curve is small. We discuss the physical
implications of these findings in §5.3.

Figure 8: The left panels show the correlation ($R^{2}$) between model
$\tau_{\rm eff}$ (top) and $n_{\rm eff}$ (bottom) vs axis ratio $b/a$ as a
function of stellar mass and SFR. The right panels show the modeled relation
between $\tau_{\rm eff}$ (top) or $n_{\rm eff}$ (bottom) and $b/a$ at the
locations in parameter space with the highest and lowest $R^{2}$ values, as
indicated by the black and red arrows, respectively. We find that inclination
is most important for dust attenuation (at least the optical depth parameter)
for massive galaxies, especially those with a high SFR. As expected from
Figure 7 and our physical understanding (§5.3), for such galaxies $\tau_{\rm
eff}$ decreases as we move toward face-on orientations. We also indirectly
observe the strong connection between $\tau_{\rm eff}$ and $n$ (Figure 5) in
the bottom right panel with the steepening of the attenuation curve as we move
toward face-on orientations. In parallel, we see that inclination has little
impact on dust attenuation curves in low-mass, low-SFR galaxies, further
supporting the picture that high-redshift, low-mass galaxies lack disk-like
morphologies.

###  4.6 Redshift Evolution of Dust Attenuation

In all 1-D and 2-D cross-sections of parameter space, there are no strong
relations between dust attenuation curves and redshift. However, the situation
changes in 3-D, as shown in Figure 9. On the left we plot the correlation
($R^{2}$) between the model parameters $\tau_{\rm eff}$ (top) or $n_{\rm eff}$
(bottom) vs redshift as a function of stellar mass and redshift. The setup of
these panels is identical to those in Figure 8.

The right panels show the model result for the dust attenuation parameters
$\tau_{\rm eff}$ (top) and $n_{\rm eff}$ (bottom) as a function of redshift
given the stellar mass and SFR values where the highest and lowest $R^{2}$
values are recorded, as indicated by the colored arrows.

We find that redshift is highly influential for galaxies with
$\log~{}M_{*}\lesssim 9.7$ and $\log~{}{\rm SFR}\lesssim 0.5$: low-mass, low-
SFR galaxies. It suggests that the dust in such galaxies evolves the most over
cosmic time. From the right panels, we observe that for the low-mass, low-SFR
galaxies, dust attenuation curves tend to have higher optical depths and
become steeper with increasing redshift, but only for $z\gtrsim 1$.

However, we can be certain about these results only for redshifts
$0.5<z\lesssim 2$, as our mass-complete sample does not extend to such low
stellar masses at the highest redshifts. To put it in quantitative terms, at
$z=2.1$ and $z=3.0$, the minimum mass in our samples is $\log~{}M_{*}=9.79$
and $\log~{}M_{*}=10.15$, respectively. The differential mass-complete range
may contribute to the seeming importance of redshift to the relation or at
least prevent the creation of a sensible model from lack of data. Pushing down
to lower masses at higher redshift, whether through better data or a more
sophisticated treatment of the selection function, may change the model
predictions at low stellar mass and SFR.

Figure 9: The left panels show the correlation ($R^{2}$) between model
$\tau_{\rm eff}$ (top) and $n_{\rm eff}$ (bottom) vs $z$ as a function of
stellar mass and SFR. The right panels show the modeled relation between
$\tau_{\rm eff}$ (top) or $n_{\rm eff}$ (bottom) and redshift at the locations
in parameter space with the highest and lowest $R^{2}$ values, as indicated by
the black and red arrows. We find that the dust attenuation curve undergoes
the strongest redshift evolution for low-mass, low-SFR galaxies. For such
galaxies, dust attenuation curves have higher optical depths and become
steeper with increasing redshift. However, the picture may be convoluted by
the differential mass-completeness limit in our sample: at $z=2.1$ and
$z=3.0$, the limit is $\log~{}M_{*}=9.79$ and $\log~{}M_{*}=10.15$,
respectively. Therefore, for $z\gtrsim 2$ the dust evolution must be confirmed
by studies that push down to lower masses.

###  4.7 High-level correlations between $n_{\rm eff}$, $\tau_{\rm eff}$,
Stellar Mass, SFR, and Metallicity

In Figure 10, we show $n_{\rm eff}$ as a function of $\tau_{\rm eff}$ and
$\log~{}M_{*}$. Once again, the center panel shows the full function
(marginalized over $\log$ SFR, $\log~{}(Z/Z_{\odot})$, and $z$), while the
left and right panels show 1-D cross-sections with $\log~{}M_{*}$ and
$\tau_{\rm eff}$, respectively. From the $n_{\rm eff}$ vs $\tau_{\rm eff}$
plot, we can very clearly see that as $\tau_{\rm eff}$ increases, $n_{\rm
eff}$ also increases for all stellar masses (with a slight decrease at high
$\tau_{\rm eff}$ when $\log~{}M_{*}\sim 9$), which is consistent with the
relationship between $A_{V}$ and $S$ in Figure 5.

The relationship between $n_{\rm eff}$ and $\log~{}M_{*}$ is much more complex
and is highly dependent on $\tau_{\rm eff}$, but in general, we see that $n$
increases (becomes shallower) and then decreases (becomes steeper) with
stellar mass.

Figure 10: Effective attenuation slope $n_{\rm eff}$ as a function of
$\log~{}M_{*}$ and $\tau_{\rm eff}$, marginalized over $\log$ SFR,
$\log~{}(Z/Z_{\odot})$, and $z$. The center panel shows the full relation,
while the left and right panels show the 1-D cross-sections for five values of
the other parameter. From the right panel, we once again confirm the result
that increasing $\tau$ leads to a shallower attenuation curve. From the left
panel we see that the dust attenuation curve first gets shallower and then
steeper with stellar mass. This can be explained by the relative amounts of
obscured and unobscured stars of different masses and dust clumpiness. The
full plausible scenario is presented in §5.4.

Next, in Figure 11, we show the correlation ($R^{2}$) of $n_{\rm eff}$ vs
$\tau_{\rm eff}$ as a function of stellar metallicity and SFR (left). In the
right panel we show the modeled relation between $n_{\rm eff}$ and $\tau_{\rm
eff}$ at the highest and lowest $R^{2}$ values. Interestingly, even though we
have repeatedly seen the dependence of $n$ on $\tau$ (which is still the case
for the highest $R^{2}$ location), from this plot we observe that the
dependence is not universal. Indeed, $\tau_{\rm eff}$ has a small but complex
and non-linear effect on $n_{\rm eff}$ when $\log~{}{\rm SFR}\lesssim 0.5$ and
$\log~{}(Z/Z_{\odot})\lesssim-0.4$.

Figure 11: Left: correlation ($R^{2}$) between $n_{\rm eff}$ and $\tau_{\rm
eff}$ as a function of stellar metallicity and SFR. Right: modeled relation
between $n_{\rm eff}$ and $\tau_{\rm eff}$ at the locations in parameter space
with the highest and lowest $R^{2}$ values, as indicated by the black and red
arrows. The relationship at the highest $R^{2}$ location is exactly what we
expect: the attenuation curve becomes flatter as we increase optical depth. On
the other hand, we observe that $\tau_{\rm eff}$ is uncorrelated with $n_{\rm
eff}$ when $\log~{}{\rm SFR}\lesssim-0.5$ and
$\log~{}(Z/Z_{\odot})\lesssim-1.2$. Specifically, the modeled relation in this
regime is quasi-sinusoidal with an amplitude small compared to the range of
$n$ (under $10\%$). This suggests that dust attenuation at the low-
metallicity, low-SFR regime behaves differently, although a strong statement
cannot be made with just the evidence we have available.

##  5 Discussion

In this section, we discuss the scientific implications of a few key results
from §4. In §5.1, we connect the modeled relationship between diffuse and
birth cloud dust optical depth to a simple physical picture depicted by Reddy
et al. (2015). Then, in §5.2, we discuss the origins of the $S$ vs $A_{V}$
relation as well as why our $S$ vs $A_{V}$ curve is shallower than that
predicted by the majority of the literature. Next, we consider the
relationship between inclination and dust attenuation over a wide range of
galaxies in §5.3. Finally, in §5.4, we analyze the complex relationships
between $n_{\rm eff}$ and the other galaxy properties.

###  5.1 Diffuse vs Birth Cloud Dust: A Closer Look

In Figure 3 we show the relation between diffuse and birth cloud dust optical
depth. As noted in §4, the relation is nearly linear but inconsistent with a
slope of unity: as the optical depth increases, birth cloud dust grows faster
than diffuse dust.

This can be explained by the connection between dust attenuation and SFR: the
optical depth tends to increase with SFR (e.g., Buat et al., 2009; Reddy et
al., 2010, 2015). Reddy et al. (2015) suggest a simple model to explain the
trend found in Figure 3: as the SFR increases, dust enrichment in birth clouds
is more prominent than in the diffuse ISM, meaning that on average massive
stars face more dust attenuation. We refer the reader to their §4.3 for a more
thorough discussion on this point (e.g., Figure 19).

###  5.2 Relationship between $A_{V}$ and Slope: Origin and Literature
Comparison

Theory and empirical evidence suggests that the dust attenuation slope becomes
shallower with increasing optical depth (Salim & Narayanan, 2020, and
references therein). The physical explanation, originating from Chevallard et
al. (2013), is related to the fact that redder light scatters more
isotropically than bluer light.

At low optical depth, a significant fraction of the galaxy integrated light we
observe arises from the equatorial plane of the galaxy, where blue light has a
harder time escaping than redder light, leading to a steeper dust attenuation
law. At high optical depth, we are observing the $\tau=1$ surface, which
includes more light outside the equatorial plane and is therefore less
dependent on the scattering differences. Also, massive stars in the $\tau=1$
vicinity may or may not be obscured by dust, leading to a shallower
attenuation slope.

In Figure 5, we show effective slope vs $A_{V}$ for our univariate model. As
mentioned in §4.3, our results are consistent with the expectations mentioned
above that dust attenuation curves become flatter as the optical depth
increases. However, our curve is shallower than the majority of literature
results, with the exceptions of Buat et al. (2012) and Reddy et al. (2015).

Qin et al. (2022) suggest that the $S$ vs $A_{V}$ relation is influenced by
degeneracies between parameters. As our hierarchical framework properly treats
correlated parameters, one may suspect that the difference between our curve
and the rest of literature originates in the statistical approach. However, we
find that even using Prospector maximum likelihood parametric values directly,
a very similar result (albeit with larger uncertainty) is derived. As the
difference in statistical modeling is ruled out as an explanation, we move
toward differences in SED modeling.

In that vein, a plausible explanation for the shallower curve is that, as
mentioned in the introduction, Prospector measured higher stellar masses and
lower current SFRs in comparison to previous studies (Leja et al., 2019),
lessening the influence of young stars and birth cloud dust. In other words,
the stars dominating the UV spectrum are slightly older and redder on average
and/or are surrounded by less birth cloud dust, leading to a shallower
effective dust attenuation curve.

###  5.3 The Role of Inclination in Regulating Dust Attenuation

For disk galaxies, inclination naturally plays a role in the observed dust
attenuation curve: edge-on lines of sight will be subject to significantly
larger dust optical depths than face-on sight lines. Chevallard et al. (2013)
found that in theoretical models, dust geometry and galaxy orientation played
a major role in the resulting dust attenuation curves. Doore et al. (2021) and
Zuckerman et al. (2021) have found observational evidence for the importance
of inclination in determining dust attenuation. Like in Zuckerman et al.
(2021), we use the more-easily-measurable axis ratio $b/a$ as a proxy for
inclination.

In §4.5, we show that our results are in agreement with expectations and the
previous results mentioned for massive galaxies. In Figures 7 and 8, we
observe that optical depth increases significantly ($0.8$ to $1.4$) from face-
on to edge-on galaxies (or more accurately, $b/a=1$ to $b/a=0$), for galaxies
with $\log~{}M_{*}\gtrsim 10.5$. In this regime, the effects of inclination
are roughly equal to those of stellar mass (over the whole range of masses) on
the dust attenuation curve.

At the same time, the lack of relation between $\tau_{\rm eff}$ and $b/a$ for
lower stellar masses is quite simple: in the high-redshift universe, low-mass
galaxies have not yet formed a disk-like shape (van der Wel et al., 2014),
thus obfuscating the relation between $b/a$ and viewing inclination. At higher
mass, especially for star-forming galaxies, the disk structure is more
pronounced, and thus the paradigm of edge-on vs face-on galaxies is
established and the importance of viewing angle comes into play.

Furthermore, the lack of distinct morphological features for high-redshift,
low-mass galaxies, especially with low SFRs explains the model prediction in
Figure 8 that inclination has the least impact on dust attenuation for low-
mass, low-SFR galaxies.

###  5.4 Interpreting the Effective Slope

In §4.3, we showed that the effective attenuation slope $n_{\rm eff}$ becomes
shallower as the optical depth $\tau_{\rm eff}$ increases, although the
relation we find is slightly flatter than that found by previous studies. In
this section, we discuss further the properties of the slope in our population
model.

The slope determines the wavelength dependence of dust attenuation, thus
playing an important role in the observed galaxy color. The reasons for
different slopes are complex and may include star-dust geometry considerations
and even chemical composition and size distribution of the dust particles
themselves.

Figure 10 shows $n_{\rm eff}$ vs $\tau_{\rm eff}$ and stellar mass. The
$n_{\rm eff}$ vs $\tau_{\rm eff}$ cross-section simply reflects what we
already showed in Figure 5: dust attenuation curves get shallower as the
optical depth increases. However, the $n_{\rm eff}$ vs stellar mass cross-
section is more complex: the attenuation curve first gets shallower and then
steeper with stellar mass, with the turning point becoming shallower and
occurring at higher stellar mass as $\tau$ increases. Theoretical studies will
need to be conducted to better understand this phenomenon.

In Figure 11, we show that $\tau$ has a small but complex effect on $n$ for
low-metallicity, low-SFR galaxies. One implication from this result is that
dust attenuation quite possibly behaves differently in the low-metallicity,
low-SFR regime. Such a finding is not surprising as chemical composition of
dust as well as production and destruction pathways likely depend on
metallicity. In addition, most of the dust in low SFR galaxies was likely
created in previous epochs, which may change the overall picture.

However, caution must be used before making too general a statement. The low-
SFR, low-metallicity regime is not well represented in the 3D-HST catalog, and
metallicities are difficult to measure from SED fitting codes. Long-exposure
spectroscopy focusing on stellar metallicity indicators such as the absorption
line rest-frame equivalent width between $1496$ and $1506\AA$ (Sommariva et
al., 2012) can help better pinpoint metallicities of galaxies in the sample
and thus better inform a dust-metallicity study.

It may also be useful to study the relation between dust attenuation and gas-
phase metallicity given, for example, the correlation in the local universe
between polycyclic aromatic hydrocarbons (PAHs)—dust molecules containing up
to 15-20% of all carbon in dust (Draine, 2003b; Zubko et al., 2004)—and gas-
phase metallicity (e.g., Engelbracht et al., 2005; Madden et al., 2006). A
study using strong line indicators (e.g., Kewley & Ellison, 2008) or better
yet [O III] $\lambda 4363$ (if not too faint) to measure gas-phase metallicity
would help investigate the connection between dust attenuation and gas-phase
metallicity.

##  6 Code for Users

By coupling the posterior samples of our four population models (§4) with a
code that parses them and calculates dust attenuation curves, we have created
the Python-based module DustE222https://github.com/Astropianist/DustE.

Users supply values for stellar mass, SFR, stellar metallicity, redshift, axis
ratio, and/or dust optical depth and select which model to be used. Any
variables not provided are marginalized over accordingly. Functions are
provided to calculate and plot dust attenuation curves for each galaxy.

This code can be used to add empirically motivated dust attenuation to
theoretical models. With a bit of tweaking, it can also be used to generate
narrower priors on dust attenuation parameters $n$ and $\tau$ for SED codes.
With better and more restrictive priors on attenuation, especially with
limited photometry and/or spectroscopy (compared to the 3D-HST photometric
catalog), SED fits should have an easier time managing the dust-age-
metallicity degeneracy.

##  7 Conclusion

By utilizing posterior samples of SED fits from the state-of-the-art code
Prospector of the rich 3D-HST data set, we have been able to create
sophisticated population Bayesian models of dust attenuation as a function of
stellar mass, SFR, metallicity, redshift, and axis ratio (proxy to
inclination).

We derive a qualitatively and quantitatively different relationship when
taking the full galaxy posteriors and the Prospector priors into account
(rather than simply adopting the maximum likelihood values), as shown in
Figure 4. This shows that modeling highly correlated and prior-sensitive
parameters is best done within a hierarchical framework. Our population model
is also able to explain more of the scatter in the data than the simple
Bayesian model.

We show the result of our model for diffuse dust optical depth as a function
of birth cloud optical depth in Figure 3. The relation is quasi-linear but
considerably different from a one-to-one relation. The birth cloud dust
optical depth increases faster than diffuse dust optical depth, which is in
agreement with the picture painted in Figure 19 of Reddy et al. (2015) in
conjunction with the connection between higher optical depths and higher SFRs
(e.g., Buat et al., 2009; Reddy et al., 2010, 2015).

We then highlight some of the interesting results and implications of our
bivariate (dependent variables optical depth $\tau$ and slope $n$) and
univariate (variable $n$) models. We mainly focus on the effective dust
attenuation models, but the diffuse dust generally behaves quite similarly.
Here are the important points gleaned from our studies.

  1. 1.

The effective dust attenuation slope flattens with increasing optical depth,
in agreement with the literature (e.g., Salim & Narayanan, 2020, and
references therein). However, the expected relation between slope and $A_{V}$
is flatter than in previous studies, which is possibly due to Prospector’s
novel empirical, non-parametric star formation history, which led to lower
current SFRs (Leja et al., 2019). As a result, the birth cloud dust
contribution is less pronounced, signifying shallower attenuation curves.

  2. 2.

The optical depth $\tau$ is most strongly correlated with SFR (positive).
However, for $\log~{}{\rm SFR}\lesssim 0$, $\tau$ is nearly constant, even
over a wide range of stellar masses.

  3. 3.

Edge-on galaxies tend to have higher $\tau$ than face-on galaxies, but only
for $\log~{}M_{*}\gtrsim 10$. We certainly expect such a relation, as for an
edge-on galaxy, the line-of-sight penetrates through a significantly larger
dust column density. However, low-mass galaxies lack disk-like morphologies.
In other words, the relation between $b/a$ and viewing inclination is
obfuscated.

  4. 4.

Dust attenuation curves evolve strongly with redshift at $z\gtrsim 1$ for low-
mass, low-SFR galaxies, with stronger but flatter curves at increasing
redshift. However, for $z>2$ we would need to push down to lower mass limits
in the data to make a firm statement.

  5. 5.

The relation between $n$ and $\log~{}M_{*}$ is multifaceted (the slope first
flattens and then steepens with mass) and difficult to explain.

Dust attenuation is an important aspect of integrated galaxy SEDs that must be
properly accounted for in order to accurately measure other physical
properties like SFR. However, its effects on the SED are similar to those of
metallicity and age. With often limited photometry, SED fitting codes struggle
to break these degeneracies. Meanwhile, ensemble dust attenuation studies are
difficult to apply to individual galaxies due to the large variations found in
the universe. Finally, there is no theoretical model of galaxy formation and
evolution that treats dust from first principles all the way from the
beginning to the current stage of the galaxy and predicts a reliable dust
attenuation curve.

By creating more sophisticated models of dust attenuation based on a
statistically rigorous framework that accounts for the large errors and often
strong covariances between parameters of SED fits, we are advancing the
knowledge of how dust attenuation curves vary from galaxy to galaxy based on
their physical properties. Our results, accessible through the Python module
DustE (§6), can be used directly by modelers who wish to include the effects
of dust when predicting observables without performing (full or approximate)
radiative transfer. In addition, they can be used as stronger priors for
individual SED fits in order to ensure more well-behaved posteriors, even when
the likelihoods are poorly constrained by a lack of narrow- or medium-band
multiwavelength photometry.

We have used the 3D-HST data set, with galaxies at $0.5<z<3.0$ over a wide
variety of stellar masses, SFRs, and metallicities. In this paper, we have
begun to explore the scientific implications of our population models
regarding dust attenuation throughout the universe. With five dimensions and a
flexible, essentially non-parametric functional form, we are able to register
both expected and surprising trends. Our results can be used as tests of dust
evolution models.

Further questions remain. First, a thorough theoretical analysis needs to be
performed to explain some of the intriguing results we presented. Beyond the
scope of this work, we restricted our analysis to the mass-complete sample,
but employing treatments for incompleteness and pushing down to faint, lower-
mass galaxies may lead to surprising changes in our picture of dust
attenuation. Also, applying our hierarchical framework to a local set of
galaxies, perhaps including a greater fraction of massive quiescent galaxies,
could possibly yield different and equally interesting science results. In
addition, simultaneously fitting dust emission using far-IR data (e.g.,
Nagaraj et al., 2021a) may provide further constraints on or even different
results for dust attenuation curves.

This material is based upon work supported by the National Science Foundation
Graduate Research Fellowship under Grant No. DGE1255832. Any opinion,
findings, and conclusions or recommendations expressed in this material are
those of the authors(s) and do not necessarily reflect the views of the
National Science Foundation. The Flatiron Institute is supported by the Simons
Foundation. We thank Rachel Somerville, Tjitske Starkenburg, and the CCA’s
Astronomical Data Group for helpful discussions.

##  Appendix A Library of Terms

In this paper, we have defined many symbols with varying levels of
abstraction. To help the reader keep track, we present a glossary that
includes concrete examples of abstract terms in Table 1.

Symbol | Meaning | Examples | §
---|---|---|---
$\\{\mathbf{w}_{k}\\}$ | Set of all true parameters of all individual galaxies |  $\log~{}M_{*}$, SFR, AGN Fraction | 3.1
$\\{\mathbf{t}_{k}\\}$ | Set of (true) parameters in our population model that are derived from other parameters in Prospector |  $\log~{}M_{*}$, SFR | 3.1
$\\{\mathbf{u}_{k}\\}$ | Set of (true) parameters in our population model that are fitted directly in Prospector |  $\tau_{2}$, $n$ | 3.1
$\\{\mathbf{v}_{k}\\}$ | Set of (true) parameters in Prospector not used in our population model that are used to compute $\\{\mathbf{t}_{k}\\}$ | SFH Parameters | 3.1
$\\{\mathbf{y}_{k}\\}$ | Set of (true) parameters in Prospector completely unrelated to those we care about in the population model | AGN Fraction | 3.1
$\\{\mathbf{x}_{k}\\}$ | Set of all data used in Prospector to fit individual galaxies | U-band photometry | 3.1
$\boldsymbol{\theta}$ | Set of all hyperparameters dealing with the population of galaxies |  $\sigma_{1}$ (see below), $M_{*}$ distribution parameters | 3.1
$\boldsymbol{\alpha}$ | Set of all parameters defining Prospector prior distributions (fixed) | Min, max values for $n$ uniform distribution | 3.1
$\tau_{1}$ | Birth cloud dust optical depth: normalization of Equation 11 | – | 3.2
$\tau_{2}$/$\tau_{\rm eff}$ | Diffuse/Effective dust optical depth: normalization of Equation 12 | – | 3.2
$n$/$n_{\rm eff}$ | Diffuse/Effective dust slope parameter: see Equation 12. Higher values mean flatter curves | – | 3.2
S | Alternate slope parameter: $S=A_{1600}/A_{V}$. Higher values mean steeper curves | – | 4.3
$I_{1}$, $I_{2}$ | Linear interpolation functions used in population dust attenuation model | – | 3.4
$\sigma_{1}$, $\sigma_{2}$, $\rho$ | Parameters of bivariate Gaussian distribution describing intrinsic scatter in the interpolation model | – | 3.4
$\mathcal{L}$ | Statistical likelihood | – | 3.4

Table 1: Glossary of mathematical symbols used in the paper. Intermediate
quantities used to simplify equations and rarely mentioned terms are not
included. For convenience, we also include the first section of appearance.

##  Appendix B Prospector Prior Sampling

In §3.1, we discussed that while the Prospector interim priors are defined as
$p(\mathbf{u},\mathbf{v},\mathbf{y}|\boldsymbol{\alpha})$ (see Appendix A for
definitions), we would like to compute
$p(\mathbf{u},\mathbf{t}|\boldsymbol{\alpha})$ instead considering
$\mathbf{u}$ and $\mathbf{t}$ are the variables involved in our population
model. To estimate $p(\mathbf{u},\mathbf{t}|\boldsymbol{\alpha})$, we sample
directly from the Prospector priors without any data and empirically compute
the probability distribution function (PDF)
$p(\mathbf{u},\mathbf{t}|\boldsymbol{\alpha})$.

We use kernel density estimation (KDE) techniques to calculate
$p(\mathbf{u},\mathbf{t}|\boldsymbol{\alpha})$ for the variables we are
considering in the particular population model. Parameters with uniform priors
are ignored as they would simply add constants in log space. For example, for
the bivariate model of $(n_{\rm eff},\tau_{\rm eff})$ as a function of stellar
mass, SFR, stellar metallicity, redshift, and axis ratio; we calculate the
joint prior PDF for stellar mass, SFR, metallicity, and $\tau_{\rm eff}$,
which have non-uniform priors. We use a Gaussian kernel with a bandwidth of
$0.2*n-0.1$, where $n$ is the number of variables in the prior. This
functional form was decided empirically based on optimization studies on prior
sub-samples.

We have sampled from the Prospector priors ($N=500,000$) at six evenly-spaced
redshifts from $z=0.5$ to $z=3.0$. We construct a PDF at each redshift using
KDE and use a linear interpolation to calculate the effective joint prior PDF
for any given galaxy sample given the redshift.

In Figure 12, we show the prior PDF for stellar mass (single variable) at
redshift $z=0.5$ (left panel) and the PDF for specific SFR (sSFR) and
metallicity (two variables) at redshift $z=3.0$ in the right panel. The
stellar mass PDF is nearly uniform, which is what we expect since stellar mass
is closely related to the total mass formed over a galaxy’s history, which has
a uniform prior distribution in Prospector. In the right panel, we see that
according to the Prospector priors, at $z=3$ higher sSFRs and metallicities
are more probable than lower ones.

Figure 12: Prior PDF with respect to stellar mass (as a single variable) at
redshift $z=0.5$ (left panel) and with respect to sSFR and metallicity (as a
double variable) at $z=3$ (right panel). In the left panel, we see that the
PDF (red line) traces the prior sample histogram (blue) very well.

##  Appendix C Validation of the Population Methodology

To ensure that our population model works properly, we simulated data and
tested how well the model could reproduce the truth. We selected anywhere from
one to four independent variables (from the list of stellar mass, sSFR,
metallicity, redshift, and axis ratio, though only the number of variables is
important). The independent variable values were chosen from a uniform random
distribution.

The true values of the dependent variables $n$ and $\tau$ were formulated with
N-D degree-2 polynomials of the independent variables such that the values
were always within the bounds introduced in §3.4. Scatter was added to the
polynomials based on the true values of $\log\sigma_{1}$, $\log\sigma_{2}$,
and $\rho$ as indicated by the equations in §3.3.

To create “posterior samples” for the simulated objects, we added a Gaussian
random variable with mean of zero and specified $\sigma$ to shift the true
values of both independent and dependent variables to an observed mean. Then,
we added another Gaussian random variable with mean of zero and the same
$\sigma$ to create samples around the observed mean. The ratio of $\sigma$
values for the different variables was set to mimic the Prospector data (based
on the average standard deviations of the true samples), but the absolute
value was determined on the command line in order to test how well our
algorithm works for both clean and noisy data.

We then modeled the simulated dust attenuation parameters (§3.2) using the
equations and techniques described in §3.3, using the interpolation model with
the same set of independent variables.

In essentially all cases, our population model is able to retrieve the true
parameter values within reasonable error bounds. Here, we show the results for
a dust attenuation model of degree 2 (truth) as a function of stellar mass and
sSFR using 6000 simulated galaxies. In Figure 13, we show the posterior
probability (left panels) of two parameters in the model (values of $n$ on the
grid and the fixed intrinsic Gaussian scatter for $n$) as well as their values
throughout the NUTS Markov chains (right panels). Given the uniformity in each
parameter within each chain (clear posterior probability peaks with Gaussian
tails and fraction-of-a-percent change in the parameters over the chains’
duration) and between chains (overlap between differently-styled lines of the
same color), we can see that the code has converged successfully.

Figure 13: Kernel density estimation (KDE) plots (left panels) and trace plots
(right panels) for simulation of dust attenuation as a degree-2 polynomial
(truth) in stellar mass and sSFR. The KDE plots show the posterior probability
of the given parameter as a function of that value. Different chains are shown
with different line styles. The top parameter, ”Grid n values,” is in fact a
conglomeration of several variables (the value of $n$ at different parts of
the grid), each one with a different color. The convergence for each parameter
within each chain and between chains is made clear through well-defined peaks
with Gaussian tails (especially visible for the scatter parameter) and the
overlap of same-color curves with different line styles. The trace plots show
the value of each variable over the duration of each chain. The variables vary
within a fraction of a percent over the entire duration of the chains, showing
a very stable solution.

Next, we show the actual fitted models and comparison between true and model
parameter values in Figure 14. The model recovers the truth very well, as
evidenced by the middle and bottom panels. The bottom panels show the true
(“posteriors”) vs measured $n$ and $\tau$ values, which are consistent with a
one-to-one relationship. The middle panels show the normalized residuals
(difference between modeled and true divided by the quadratic sum of intrinsic
scatter in the truth and model).

There are no visible systematics. We do slightly overestimate the intrinsic
scatter in the relationship ($\log\sigma_{1}=\log\sigma_{2}=-5.30$ truth vs
recovered values of $\log\sigma_{1}=-4.95$ and $\log\sigma_{2}=-5.09$). This
$\sim 30\%$ overestimate in intrinsic scatter arises because the simulation is
unable to perfectly distinguish between measurement errors and intrinsic
scatter when the errors (taken to be of the same order of magnitude as the
intrinsic scatter) are only $5\%$ of the Prospector measurement errors. We
find for measurement errors on the order of those in Prospector (not plotted
in this paper), as long as the intrinsic scatter is over $\sim 20\%$ of the
typical Prospector errors, we recover even the scatter with no bias.

Figure 14: Test of our population modeling capabilities. For the true dust
attenuation, we set both $n$ and $\tau_{2}$ as degree-2 polynomials of log
stellar mass and log sSFR and create 6000 “galaxies.” We simulate with our
interpolation model. The top two panels show the mean model for $n$ (left) and
$\tau_{2}$ (right) as a function of $\log~{}M_{*}$ and $\log{\rm sSFR}$. The
middle two panels show the normalized residual maps (difference between model
and truth divided by overall intrinsic scatter) in the same parameter space.
The red dots in this case represent the grid on which the polynomial model was
calculated. The lack of visible systematics suggests that the model is
successful. In this particular model, we do slightly overestimate the
intrinsic scatter (simulated with $\log\sigma_{1}=\log\sigma_{2}=-5.30$ but
recovered $\log\sigma_{1}=-4.95$ and $\log\sigma_{2}=-5.09$), but given how
small the “measurement errors” are ($\sim 5\%$ of the Prospector errors), this
is not a cause for concern. The bottom two panels show the comparison between
true $n$ (left) and $\tau_{2}$ (right) parameters and model parameters. We see
excellent agreement between the truth and the model, especially for $n$.

##  Appendix D Marginalization

At various points in this work, we have visualized our model in a lower-
dimensional space than its native 5 dimensions of independent parameters. Here
we explain the not-necessarily-straightforward choices we make to do so.

The hyperparameters of our model, $\boldsymbol{\theta}$, are mostly the values
of the dependent variable(s) at control points used for interpolation. Given a
particular choice of $\boldsymbol{\theta}$, our model predicts a distribution
of some dependent variable(s) at every point in the five-dimensional space of
independent variables. The independent variables are things like $M_{*}$, SFR,
and so on, and the dependent variables are $n$, $n_{\mathrm{eff}}$, $n$ and
$\tau$, or $n_{\mathrm{eff}}$ and $\tau_{\mathrm{eff}}$. When we visualize the
model’s predictions as a function of, say, just $M_{*}$ and SFR, we have to
make a choice about what to do with the remaining dimensions, e.g. $z$, $b/a$,
and $Z_{*}$. In analogy with visualizing a 3D cube of data in 2 dimensions,
one could imagine 2 basic options. The first is to fix the values of these
other parameters to create a 2D slice of the prediction. The second is to
integrate over the other parameters to create a projection, possibly with some
weight other than simply the “volume” of the parameter space, which may be
somewhat arbitrary, e.g. when integrating over $M_{*}$, should one integrate
over $M_{*}$ or $\log M_{*}$?

To avoid slices which would presumably be quite sensitive to the particular
fixed values of the suppressed independent parameters, we have chosen to use a
projection. In some sense, therefore, our task now is to choose a reasonable
weight function and elucidate exactly how we approximate the corresponding
integrals as sums. We also need to keep in mind that the model predicts a
gaussian distribution at every point in the 5-dimensional independent
parameter space, not just a single value.

Let’s define the model’s predicted PDF at such a point in parameter space to
be $f(\delta|{\bf X},\boldsymbol{\theta})$. The dust parameters are denoted
$\delta$ (which will be one- or two-dimensional depending on the flavor of the
model), the point in the 5-dimensional space of independent variables is ${\bf
X}$, and the parameters of the population model are $\boldsymbol{\theta}$. At
a fixed value of ${\bf X}$, there is no ambiguity about the mean value of
$\delta$:

$E_{\boldsymbol{\theta}|\\{\mathbf{x}_{k}\\}}[\delta|{\bf X}]=\int\delta\
f(\delta|{\bf
X},\boldsymbol{\theta})p(\boldsymbol{\theta}|\\{\mathbf{x}_{k}\\})d\boldsymbol{\theta}d\delta.$
(D1)

The integral over $\delta$ is necessary because $f$ is itself a PDF. The
second factor, $p(\boldsymbol{\theta}|\\{\mathbf{x}_{k}\\})$ is the posterior
distribution of the model’s hyperparameters $\boldsymbol{\theta}$. Using the
usual monte carlo approximation of this integral, we would estimate

$E_{\boldsymbol{\theta}|\\{\mathbf{x}_{k}\\}}[\delta|{\bf X}]\approx\int
M^{-1}\sum_{m=1}^{M}\delta\ f(\delta|{\bf
X},\boldsymbol{\theta}^{(m)})d\delta\approx
M^{-1}\sum_{m=1}^{M}\delta_{\mathrm{avg},\boldsymbol{\theta}^{(m)}},$ (D2)

where $m$ indexes the $M$ posterior samples of our population model, and the
second approximation holds if the width of the distribution $f$ is insensitive
to $\boldsymbol{\theta}$, which we find to be the case in practice.
$\delta_{\mathrm{avg}}$ is a straightforward function of
$\boldsymbol{\theta}$, i.e. the result of the N-dimensional interpolations, so
this approximation is particularly convenient.

If we now split ${\bf X}$ into ${\bf X}_{v}$ and ${\bf X}_{s}$, which will be
respectively visualized and suppressed, the quantity we should compute is of
the form

$E_{{\bf X}_{s},\boldsymbol{\theta}|\\{\mathbf{x}_{k}\\}}[\delta|{\bf
X}_{v}]=\int\delta\ f(\delta|{\bf
X},\boldsymbol{\theta})p(\boldsymbol{\theta}|\\{\mathbf{x}_{k}\\})W({\bf
X}_{s}|\boldsymbol{\theta},\\{\mathbf{x}_{k}\\})d\boldsymbol{\theta}d\delta
d{\bf X}_{s}.$ (D3)

A reasonable choice for the weight function $W({\bf
X}_{s}|\boldsymbol{\theta},\\{x_{i}\\})$, which in the general case may depend
on $\boldsymbol{\theta}$ and $\\{x_{i}\\}$, is the density of data points. In
principle this should have the advantage of not incorporating regions of
parameter space with large uncertainty, i.e. regions not constrained by data,
into the visualized averages. However, the true density of galaxies at any
point ${\bf X}$ is nontrivial to determine (see Foreman-Mackey et al., 2014)
given the uncertainties in the data. These uncertainties are encapsulated in
the posterior samples we have from Prospector for each galaxy, which are taken
to be drawn from the distribution
$p(\\{\mathbf{t}_{k}\\},\\{\mathbf{u}_{k}\\}|\\{\mathbf{x}_{k}\\},\boldsymbol{\alpha})$
(related to Equation 1), where again $\boldsymbol{\alpha}$ represents the
parameters controlling the interim prior used by Prospector. Note that ${\bf
X}$ is a subset of the parameters sets $\mathbf{t}_{k}$ and $\mathbf{u}_{k}$,
so $p({\bf X}_{s}|\\{\mathbf{x}_{k}\\},\boldsymbol{\alpha})$ is a marginal
distribution of
$p(\\{\mathbf{t}_{k}\\},\\{\mathbf{u}_{k}\\}|\\{\mathbf{x}_{k}\\},\boldsymbol{\alpha})$.
However, just as in the main text, we are interested in this distribution with
the effects of the interim prior $\boldsymbol{\alpha}$ removed.

To approximate this integral, we employ similar methodology to our discussion
in Section 3. In particular, we rely on sums over the Prospector samples to
perform the integral over ${\bf X}_{s}$, sums over the posterior samples of
the population model to integrate over $\boldsymbol{\theta}$, and given the
relative lack of posterior uncertainty in the parameters controlling the width
of the normal distributions $f(\delta|{\bf X},\boldsymbol{\theta})$, the
integral over $\delta$ is approximated by neglecting this width as in the
second approximation in Equation (D2). We therefore approximate this average
as

$E_{{\bf X}_{s},\boldsymbol{\theta}|\\{\mathbf{x}_{k}\\}}[\delta|{\bf
X}_{v}]\approx
M^{-1}Q^{-1}\sum_{m=1}^{M}\sum_{q=1}^{Q}\delta_{\mathrm{avg},\boldsymbol{\theta}^{(m)}}.$
(D4)

Here again $m$ indexes samples from the population model. The $Q$ points,
meanwhile, are drawn from the sample of $N$ galaxies times the (up to) 500
Prospector samples available per galaxy, with a weight equal to
$1/p(\mathbf{t}_{q},\mathbf{u}_{q}|\boldsymbol{\alpha})$, i.e. the inverse of
the Prospector priors at the given sample point. In the language of slices and
projections, this is a projection constructed by averaging many slices through
the parameter space using a fine grid of ${\bf X}_{v}$ values, where the fixed
values of the suppressed dimensions ${\bf X}_{s}$ are drawn from the
population of Prospector posterior points with the Prospector priors removed.

## References

  * Acquaviva et al. (2011) Acquaviva, V., Gawiser, E., & Guaita, L. 2011, ApJ, 737, 47
  * Álvarez-Márquez et al. (2019) Álvarez-Márquez, J., Burgarella, D., Buat, V., Ilbert, O., & Pérez-González, P. G. 2019, A&A, 630, A153
  * Arnouts et al. (2013) Arnouts, S., Le Floc’h, E., Chevallard, J., et al. 2013, A&A, 558, A67
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123
  * Battisti et al. (2016) Battisti, A. J., Calzetti, D., & Chary, R. R. 2016, ApJ, 818, 13
  * Battisti et al. (2017) —. 2017, ApJ, 840, 109
  * Battisti et al. (2019) Battisti, A. J., da Cunha, E., Grasha, K., et al. 2019, ApJ, 882, 61
  * Behroozi et al. (2013) Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57
  * Boquien et al. (2019) Boquien, M., Burgarella, D., Roehlly, Y., et al. 2019, A&A, 622, A103
  * Bouwens et al. (2012) Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2012, ApJ, 754, 83
  * Bouwens et al. (2015) —. 2015, ApJ, 803, 34
  * Bowman et al. (2020) Bowman, W. P., Zeimann, G. R., Nagaraj, G., et al. 2020, ApJ, 899, 7
  * Brammer et al. (2012) Brammer, G. B., van Dokkum, P. G., Franx, M., et al. 2012, ApJS, 200, 13
  * Buat et al. (2009) Buat, V., Takeuchi, T. T., Burgarella, D., Giovannoli, E., & Murata, K. L. 2009, A&A, 507, 693
  * Buat et al. (2012) Buat, V., Noll, S., Burgarella, D., et al. 2012, A&A, 545, A141
  * Burgarella et al. (2005) Burgarella, D., Buat, V., & Iglesias-Páramo, J. 2005, MNRAS, 360, 1413
  * Calzetti (2001) Calzetti, D. 2001, PASP, 113, 1449
  * Calzetti et al. (2000) Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682
  * Calzetti et al. (1994) Calzetti, D., Kinney, A. L., & Storchi-Bergmann, T. 1994, ApJ, 429, 582
  * Calzetti et al. (1997) Calzetti, D., Meurer, G. R., Bohlin, R. C., et al. 1997, AJ, 114, 1834
  * Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763
  * Charlot & Fall (2000) Charlot, S., & Fall, S. M. 2000, ApJ, 539, 718
  * Chevallard & Charlot (2016) Chevallard, J., & Charlot, S. 2016, MNRAS, 462, 1415
  * Chevallard et al. (2013) Chevallard, J., Charlot, S., Wandelt, B., & Wild, V. 2013, MNRAS, 432, 2061
  * Conroy (2013) Conroy, C. 2013, ARA&A, 51, 393
  * Conroy & Gunn (2010) Conroy, C., & Gunn, J. E. 2010, ApJ, 712, 833
  * Conroy et al. (2009) Conroy, C., Gunn, J. E., & White, M. 2009, ApJ, 699, 486
  * da Cunha et al. (2008) da Cunha, E., Charlot, S., & Elbaz, D. 2008, MNRAS, 388, 1595
  * Daddi et al. (2007) Daddi, E., Dickinson, M., Morrison, G., et al. 2007, ApJ, 670, 156
  * Doore et al. (2021) Doore, K., Eufrasio, R. T., Lehmer, B. D., et al. 2021, arXiv e-prints, arXiv:2109.05039
  * Draine (2003a) Draine, B. T. 2003a, ARA&A, 41, 241
  * Draine (2003b) —. 2003b, ApJ, 598, 1017
  * Driver et al. (2018) Driver, S. P., Andrews, S. K., da Cunha, E., et al. 2018, MNRAS, 475, 2891
  * Elbaz et al. (2011) Elbaz, D., Dickinson, M., Hwang, H. S., et al. 2011, A&A, 533, A119
  * Engelbracht et al. (2005) Engelbracht, C. W., Gordon, K. D., Rieke, G. H., et al. 2005, ApJ, 628, L29
  * Erb et al. (2006) Erb, D. K., Shapley, A. E., Pettini, M., et al. 2006, ApJ, 644, 813
  * Finkelstein et al. (2012) Finkelstein, S. L., Papovich, C., Salmon, B., et al. 2012, ApJ, 756, 164
  * Forbes et al. (2019) Forbes, J. C., Krumholz, M. R., & Speagle, J. S. 2019, MNRAS, 487, 3581
  * Foreman-Mackey et al. (2014) Foreman-Mackey, D., Hogg, D. W., & Morton, T. D. 2014, ApJ, 795, 64
  * Fritz et al. (2007) Fritz, J., Poggianti, B. M., Bettoni, D., et al. 2007, A&A, 470, 137
  * Galliano (2018) Galliano, F. 2018, MNRAS, 476, 1445
  * Grogin et al. (2011) Grogin, N. A., Kocevski, D. D., Faber, S. M., et al. 2011, ApJS, 197, 35
  * Han & Han (2014) Han, Y., & Han, Z. 2014, ApJS, 215, 2
  * Hoffman & Gelman (2011) Hoffman, M. D., & Gelman, A. 2011, arXiv e-prints, arXiv:1111.4246
  * Inoue (2005) Inoue, A. K. 2005, MNRAS, 359, 171
  * Johnson et al. (2021) Johnson, B. D., Leja, J., Conroy, C., & Speagle, J. S. 2021, ApJS, 254, 22
  * Johnson et al. (2007) Johnson, B. D., Schiminovich, D., Seibert, M., et al. 2007, ApJS, 173, 392
  * Jones et al. (2001–) Jones, E., Oliphant, T., Peterson, P., et al. 2001–, SciPy: Open source scientific tools for Python
  * Kelly et al. (2012) Kelly, B. C., Shetty, R., Stutz, A. M., et al. 2012, ApJ, 752, 55
  * Kewley & Ellison (2008) Kewley, L. J., & Ellison, S. L. 2008, ApJ, 681, 1183
  * Koekemoer et al. (2011) Koekemoer, A. M., Faber, S. M., Ferguson, H. C., et al. 2011, ApJS, 197, 36
  * Kriek & Conroy (2013) Kriek, M., & Conroy, C. 2013, ApJ, 775, L16
  * Leja et al. (2017) Leja, J., Johnson, B. D., Conroy, C., van Dokkum, P. G., & Byler, N. 2017, ApJ, 837, 170
  * Leja et al. (2020) Leja, J., Speagle, J. S., Johnson, B. D., et al. 2020, ApJ, 893, 111
  * Leja et al. (2019) Leja, J., Johnson, B. D., Conroy, C., et al. 2019, ApJ, 877, 140
  * Lilly et al. (2013) Lilly, S. J., Carollo, C. M., Pipino, A., Renzini, A., & Peng, Y. 2013, ApJ, 772, 119
  * Madau & Dickinson (2014) Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415
  * Madden et al. (2006) Madden, S. C., Galliano, F., Jones, A. P., & Sauvage, M. 2006, A&A, 446, 877
  * Mannucci et al. (2010) Mannucci, F., Cresci, G., Maiolino, R., Marconi, A., & Gnerucci, A. 2010, MNRAS, 408, 2115
  * Meurer et al. (1999) Meurer, G. R., Heckman, T. M., & Calzetti, D. 1999, ApJ, 521, 64
  * Momcheva et al. (2016) Momcheva, I. G., Brammer, G. B., van Dokkum, P. G., et al. 2016, ApJS, 225, 27
  * Nagaraj et al. (2021a) Nagaraj, G., Ciardullo, R., Bowman, W. P., & Gronwall, C. 2021a, ApJ, 913, 34
  * Nagaraj et al. (2021b) Nagaraj, G., Ciardullo, R., Lawson, A., et al. 2021b, ApJ, 912, 145
  * Noll et al. (2009) Noll, S., Pierini, D., Cimatti, A., et al. 2009, A&A, 499, 69
  * Oesch et al. (2013) Oesch, P. A., Bouwens, R. J., Illingworth, G. D., et al. 2013, ApJ, 773, 75
  * Peng et al. (2002) Peng, C. Y., Ho, L. C., Impey, C. D., & Rix, H.-W. 2002, AJ, 124, 266
  * Peng et al. (2010a) —. 2010a, AJ, 139, 2097
  * Peng et al. (2010b) Peng, Y.-j., Lilly, S. J., Kovač, K., et al. 2010b, ApJ, 721, 193
  * Prevot et al. (1984) Prevot, M. L., Lequeux, J., Maurice, E., Prevot, L., & Rocca-Volmerange, B. 1984, A&A, 132, 389
  * Price et al. (2014) Price, S. H., Kriek, M., Brammer, G. B., et al. 2014, ApJ, 788, 86
  * Qin et al. (2022) Qin, J., Zheng, X. Z., Fang, M., et al. 2022, arXiv e-prints, arXiv:2201.05467
  * Reddy et al. (2010) Reddy, N. A., Erb, D. K., Pettini, M., Steidel, C. C., & Shapley, A. E. 2010, ApJ, 712, 1070
  * Reddy et al. (2015) Reddy, N. A., Kriek, M., Shapley, A. E., et al. 2015, ApJ, 806, 259
  * Salim et al. (2018) Salim, S., Boquien, M., & Lee, J. C. 2018, ApJ, 859, 11
  * Salim & Narayanan (2020) Salim, S., & Narayanan, D. 2020, ARA&A, 58, 529
  * Salim et al. (2016) Salim, S., Lee, J. C., Janowiecki, S., et al. 2016, ApJS, 227, 2
  * Salmon et al. (2016) Salmon, B., Papovich, C., Long, J., et al. 2016, ApJ, 827, 20
  * Salvatier et al. (2016) Salvatier, J., Wiecki, T. V., & Fonnesbeck, C. 2016, PeerJ Computer Science, 2, e55
  * Santini et al. (2015) Santini, P., Ferguson, H. C., Fontana, A., et al. 2015, ApJ, 801, 97
  * Shapley et al. (2001) Shapley, A. E., Steidel, C. C., Adelberger, K. L., et al. 2001, ApJ, 562, 95
  * Shivaei et al. (2015) Shivaei, I., Reddy, N. A., Shapley, A. E., et al. 2015, ApJ, 815, 98
  * Skelton et al. (2014) Skelton, R. E., Whitaker, K. E., Momcheva, I. G., et al. 2014, ApJS, 214, 24
  * Somerville et al. (2015) Somerville, R. S., Popping, G., & Trager, S. C. 2015, MNRAS, 453, 4337
  * Sommariva et al. (2012) Sommariva, V., Mannucci, F., Cresci, G., et al. 2012, A&A, 539, A136
  * Stecher (1965) Stecher, T. P. 1965, ApJ, 142, 1683
  * Tal et al. (2014) Tal, T., Dekel, A., Oesch, P., et al. 2014, ApJ, 789, 164
  * Tress et al. (2018) Tress, M., Mármol-Queraltó, E., Ferreras, I., et al. 2018, MNRAS, 475, 2363
  * van der Wel et al. (2014) van der Wel, A., Franx, M., van Dokkum, P. G., et al. 2014, ApJ, 788, 28
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261
  * Walcher et al. (2011) Walcher, J., Groves, B., Budavári, T., & Dale, D. 2011, Ap&SS, 331, 1
  * Wild et al. (2011) Wild, V., Charlot, S., Brinchmann, J., et al. 2011, MNRAS, 417, 1760
  * Witt & Gordon (2000) Witt, A. N., & Gordon, K. D. 2000, ApJ, 528, 799
  * Zubko et al. (2004) Zubko, V., Dwek, E., & Arendt, R. G. 2004, ApJS, 152, 211
  * Zuckerman et al. (2021) Zuckerman, L. D., Belli, S., Leja, J., & Tacchella, S. 2021, ApJ, 922, L32