# Empirical Dust Attenuation Model Leads to More Realistic UVJ Diagram for
TNG100 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 Gravitation and the Cosmos, The Pennsylvania State University, University
Park, PA 16802, USA  Institute for Computational & Data Sciences, The
Pennsylvania State University, University Park, PA 16802, USA  Dan 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 attenuation varies substantially from galaxy to galaxy and as of yet
cannot be reproduced from first principles in theoretical models. In Nagaraj
et al. (2022), we developed the first Bayesian population model of dust
attenuation as a function of stellar population properties and projected
galaxy shape, built on spectral energy distribution (SED) fits of nearly
30,000 galaxies in the 3D-HST grism survey with broadband photometric coverage
from the rest-frame UV to IR. In this paper, we apply the model to galaxies
from the large-volume cosmological simulation TNG100. We produce a UVJ diagram
and compare it with one obtained in previous work by applying approximate
radiative transfer to the simulated galaxies. We find that the UVJ diagram
based on our empirical model is in better agreement with observations than the
previous effort, especially in the number density of dusty star forming
galaxies. We also construct the intrinsic dust-free UVJ diagram for TNG and
3D-HST galaxies at $z\sim 1$, finding qualitative agreement but residual
differences at the $10-20\%$ level. These differences can be largely
attributed to the finding that TNG galaxies have, on average, $29\%$ younger
stellar populations and $0.28$ dex higher metallicities than observed
galaxies.

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., 2018), SciPy (Virtanen et al.,
2020), FSPS (Conroy et al., 2009; Conroy & Gunn, 2010), Prospector (Leja et
al., 2017; Johnson et al., 2021)

##  1 Introduction

Dust plays a major role in shaping the ultraviolet (UV) through infrared (IR)
portion of galaxy spectra, absorbing and scattering UV through near-IR (NIR)
light and re-radiating it in the mid- and far-IR (MIR and FIR, e.g., Draine,
2003). As dust preferentially absorbs and scatters bluer light, it causes not
only extinction but also reddening.

Due to the small angular size covered by cosmologically distant galaxies, we
typically are able to observe only integrated photometry and spectroscopy,
rather than resolved lines of sight. In a galaxy, the complex configuration of
stars and dust causes differential levels of obscuration, compounded by
scattering both into and out of different lines of sight. By re-imagining the
scenario as a dust screen (i.e., a sheet) in front of the stars (e.g.,
Calzetti et al., 1994, 2000; Chevallard et al., 2013), we can create an
effective dust extinction law for the galaxy, designated as dust attenuation.

Dust attenuation has been shown to vary widely from galaxy to galaxy (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). For a given galaxy, dust
attenuation is difficult to measure given the similar reddening effects it has
on the galaxy spectrum as increasing metallicity or average stellar age (e.g.,
Conroy, 2013; Santini et al., 2015). In fact, dust attenuation is a major
source of uncertainty in high redshift star formation studies (e.g., Bouwens
et al., 2012; Finkelstein et al., 2012; Oesch et al., 2013).

While our understanding of dust grains in the local universe is substantial
(Galliano et al., 2018, and references therein), the reduced level of
knowledge at high redshift coupled with the endless intricate possibilities of
star-dust geometry severely limits our theoretical foundations for dust
attenuation in the early universe. Therefore, theoretical efforts in
understanding galaxy evolution, from analytic models (e.g., Forbes et al.,
2014, 2019) to hydrodynamical simulations (e.g., Genel et al., 2014; Pillepich
et al., 2018a; Davé et al., 2019), often exclude dust physics. For such
models, an empirically based dust attenuation prescription written as a
function of physical properties of galaxies such as stellar mass, star
formation rate (SFR), metallicity, and redshift would be extremely useful.

Many astronomers have explored the connection between dust attenuation curves
and galaxy properties, 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). Nevertheless, a complete understanding of dust
attenuation is still far from realization (Salim & Narayanan, 2020, and
references therein).

Nagaraj et al. (2022), hereafter referred to as Paper I, describes the first
Bayesian population model for dust attenuation (as a function of stellar
population properties and projected galaxy shape) developed from state-of-the-
art spectral energy distribution (SED) fits (Prospector, Leja et al., 2017,
2019b, 2020) of a mass-complete set of $\sim 29,000$ $0.5<z<3.0$ galaxies from
the 3D-HST grism survey (Brammer et al., 2012; Skelton et al., 2014; Momcheva
et al., 2016). By combining the grism redshifts with copious multiwavelength
photometry from CANDELS (Grogin et al., 2011; Koekemoer et al., 2011), the
3D-HST survey allows for more accurate inferred galaxy properties. The multi-
layered aspect of the model properly accounts (Foreman-Mackey et al., 2014a)
for the large uncertainties and covariances among galaxy parameters (e.g.,
Conroy, 2013; Santini et al., 2015), making the model more robust against
complex correlations between inferred galaxy properties. We demonstrate this
improved accuracy in a test in Paper I.

In this paper, we highlight how to use the (unresolved) population model and
demonstrate that it is more effective at reproducing the observed galaxy
colors than a first-principles approach (Nelson et al., 2018; Donnari et al.,
2019) by applying it to subhalos (galaxies) from the IllustrisTNG simulations
(Pillepich et al., 2018a) at redshift $z=1$ and comparing the $U-V$ and $V-J$
colors we derive to those of Donnari et al. (2019) and to an observed set of
colors from 3D-HST. Having an accurate working model for dust attenuation
allows us to identify specific areas of improvement for theoretical models of
galaxy evolution by comparing to observations since cosmological simulations
are not yet at a level to reproduce dust attenuation from first principles.

In §2, we present a visualization of the model to show how it can be directly
applied to theoretical models to produce attenuated spectra. In §3, we
introduce the UVJ diagram and discuss how we calculated colors for TNG100
galaxies. Then, in §4 we present the UVJ colors of TNG100 and compare them to
both previous results from Donnari et al. (2019) and the measured 3D-HST
colors. We discuss our results in §5 and conclude in §6.

To ensure consistency with Nelson et al. (2018) and Donnari et al. (2019), we
use the cosmological parameters from Planck Collaboration et al. (2016),
namely $\Omega_{m,0}=0.3089$, $\Omega_{\Lambda,0}=0.6911$,
$\Omega_{b,0}=0.0486$, and $H_{0}=67.74$ $\text{km s}^{-1}\text{ Mpc}^{-1}$;
the Chabrier (2003) IMF; and the AB magnitude system (Oke, 1974).

##  2 Using the Population Model

Our population models of dust
attenuation111https://github.com/Astropianist/DustE, as described in detail in
Paper I, describe both the two-component dust attenuation curve of Prospector
(Leja et al., 2017; Johnson et al., 2021), inspired by Charlot & Fall (2000),
and an effective attenuation curve (single component). In the two-component
model, the birth cloud dust attenuation ($\tau_{\lambda,1}$), which affects
only stars less than 10 million years old, is modeled as
$\tau_{\lambda,1}\propto\lambda^{-1}$. The diffuse dust attenuation
($\tau_{\lambda,2}$) is formulated according to Noll et al. (2009) and Kriek &
Conroy (2013). This is a flexible extension of the attenuation curve from
Calzetti et al. (2000), or “Calzetti Law” (represented as
$k^{\prime}(\lambda)$ below) with parameters for slope and 2175 Å bump
strength. Prospector fixes the bump strength $E_{b}$ based on results from
Kriek & Conroy (2013).

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

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

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

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

The parameters fitted in the two-component model are $n$ (related to slope),
$\tau_{2}$, and $\tau_{1}$. In the effective dust attenuation model, the curve
is parameterized by Equation 2, with parameters $n_{\rm eff}$ and $\tau_{\rm
eff}$. In this work, we use the two-component dust attenuation model, in which
$\tau_{1}$ is a linear interpolation function of $\tau_{2}$ while $n$ and
$\tau_{2}$ are fitted as a linear interpolation function of stellar mass, SFR,
stellar metallicity, redshift, and axis ratio (a proxy for inclination). Both
are models described in Paper I.

The package we have created to access the population dust models is quite
simple to use, and we provide a number of demos in the footnote link given
earlier. A description of the code is also presented in Paper I.

As a demonstration of our models, Figure 1 shows the evolution of the diffuse
dust attenuation curve over the stellar mass – SFR parameter space. We observe
that the attenuation generally increases (larger optical depth and typically
steeper slope) with stellar mass and SFR, but the exact relation between the
curves and the parameter space are quite complex. Discussions on the
scientific implications of our models can be found in Paper I.

Also, in Paper I, using mock tests, we show that the population model is
substantially more accurate than the common procedure of binning in, or even
fitting a Bayesian model to, point estimates of each galaxy’s properties. By
averaging the likelihood over the full posterior of each object, we can avoid
our results being driven by well-known correlations between parameters in
individual fits (e.g., $\tau_{2}$ and $n$) and thus produce more reliable
inferred relationships. With 5-D dust attenuation models, there is a plethora
of interesting trends that are still to be gleaned, making the package we
provide useful for learning more about dust attenuation in addition to the
purpose it serves in providing an interface between theory and observations as
demonstrated in this work.

Figure 1: Variation in the diffuse dust attenuation curve over the stellar
mass – SFR parameter space. The attenuation curves correspond to the model’s
prediction at the location in parameter space at the center of each small plot
(as highlighted by the dotted black lines). The $1-\sigma$ errors on the mean
attenuation curves are shaded but are too small to be visible. The faint
points in the background show the mean Prospector posterior sample for stellar
mass and SFR and are colored by sSFR, as are the attenuation curves
themselves. For easier interpretation, in each attenuation curve plot, we show
the curve for the other eight locations as gray dashed lines. We can see that
attenuation generally increases (larger optical depth and typically steeper
attenuation slope) with stellar mass and SFR, but the detailed behavior is
rather complex. See Paper I for extensive discussions of the science results
of our models.

##  3 Methods

###  3.1 Introduction to the UVJ Diagram

It is important to be able to easily distinguish between galaxies that are red
due to quiescence and those that are red due to dust. One simple observational
technique is to place the galaxies on a rest-frame $U-V$ vs. $V-J$ color-color
plot (Wuyts et al., 2007). While a single color ($U-V$ or $V-J$) produces a
degeneracy between intrinsically red objects and systems that have been
reddened by dust, Wuyts et al. (2007) showed that the two galaxy classes
separate in this two-color space, as dusty star forming galaxies tend to have
redder $V-J$ colors than their quiescent counterparts. Many studies have used
rest-frame $UVJ$ colors to distinguish quiescent and star-forming galaxies
(e.g., Williams et al., 2009; Patel et al., 2012; Whitaker et al., 2012;
Muzzin et al., 2013; Fumagalli et al., 2014; Bowman et al., 2019).

More fundamentally, quiescence can be defined via specific SFR (sSFR), though
observational inferences of sSFR are more subject to systematic effects than
rest-frame colors. Of course, in simulations, the true sSFR is known. However,
there is no precise definition for quiescence, and sSFR limits range from
$10^{-10}$ yr-1 (e.g., Whitaker et al., 2017; Wu et al., 2018) to $10^{-11}$
yr-1 (e.g., Cecchi et al., 2019). Nevertheless, Leja et al. (2019a) have shown
that the standard UVJ cuts correspond to $\log~{}{({\rm
sSFR/yr}^{-1})}\sim-10.0$.

###  3.2 Measuring Magnitudes for TNG100 Subhalos

To exhibit our population dust attenuation models’ validity, we apply them to
simulated galaxies from the IllustrisTNG project222https://www.tng-
project.org/ (Pillepich et al., 2018a), namely TNG100 (Marinacci et al., 2018;
Naiman et al., 2018; Nelson et al., 2018; Pillepich et al., 2018b; Springel et
al., 2018), and compare the derived UVJ colors to those of Donnari et al.
(2019) and the observations of 3D-HST galaxies from Prospector (Leja et al.,
2019b, 2020). The results are shown in §4.

Donnari et al. (2019) apply the Nelson et al. (2018) resolved dust model C to
TNG simulated galaxies and measure the magnitudes $U$, $V$, and $J$ of the
resulting spectra. To ensure a suitable comparison to their work, we employed
a similar process to that outlined in Nelson et al. (2018) and Donnari et al.
(2019) to measure colors of TNG galaxies, with the major exception being the
treatment of dust. We took all TNG100 $z=1$ (snapshot 50) subhalos with
stellar masses of at least $10^{9}\hbox{$M_{\odot}$}$ and summed up the
magnitudes in the U and V (Bessell, 1990) and J (Skrutskie et al., 2006)
filters for the stellar particles within $2r_{1/2}$ in each subhalo, where
$r_{1/2}$ is the half-mass radius.

To model spectra, we used the FSPS stellar population synthesis code (Conroy
et al., 2009; Conroy & Gunn, 2010; Foreman-Mackey et al., 2014b) with a
Chabrier IMF 2003, the Padova isochrones, and the MILES spectral library. To
speed up the computation process, we calculated U, V, and J magnitudes using
FSPS for each of the metallicity and age steps in the Padova database. We
treated all of these grid points as simple stellar populations (single-age,
single-metallicity populations). We included nebular emission (Byler et al.,
2017), with the gas-phase metallicity equal to the stellar metallicity for
self consistency and with ionization parameter fixed at $\log~{}U=-2.0$.

Again following the steps outlined in Nelson et al. (2018) and Donnari et al.
(2019), for each stellar particle within $2r_{1/2}$ (twice the stellar half-
mass radius) of a given subhalo, we used a bicubic spline over log metallicity
and log age to approximate the magnitudes and then added up the fluxes in
proportion to the mass of each particle. Metallicities and ages outside the
grid bounds yield the magnitudes at the nearest boundary point.

The aforementioned process measures dust-free magnitudes for the stellar
populations in each subhalo. To include dust, we looked up the attenuation
parameters of our population model given the values of stellar mass, SFR,
metallicity, redshift, and axis ratio $b/a$ (proxy for inclination) for each
subhalo. We used the average SFR over the last 100 million years of the star
formation history (SFH) and the mass-weighted metallicity for the calculation
to emulate observations of each quantity.

Axis ratio depends on viewing angle, which is arbitrary for the simulated
galaxies, so we used the following procedure to select axis ratios for TNG
galaxies. The connection between axis ratio and inclination becomes tenuous
for galaxies with triaxial morphologies (van der Wel et al., 2014). Galaxies
with lower SFRs and/or stellar masses tend to fall in this category.
Therefore, the values chosen for axis ratio in this regime do not matter.
Thus, for galaxies with $\log({\rm
SFR}/[\hbox{$M_{\odot}~{}\text{yr}^{-1}$}])<1$, we simply drew values for
$b/a$ from a uniform distribution between $0$ and $1$.

For galaxies with $\log({\rm SFR}/[\hbox{$M_{\odot}~{}\text{yr}^{-1}$}])\geq
1$, we employed the simple triaxial model from Zuckerman et al. (2021). The
sources are treated as ellipsoids with axes parameters $a^{\prime}$,
$b^{\prime}=a^{\prime}$, and $c^{\prime}$ (thickness). We let the ratio of
thickness $c^{\prime}$ to the semi-major axis $a^{\prime}$ be $D\equiv
c^{\prime}/a^{\prime}=0.1$. Then, for a given inclination $\phi$, where $0$
corresponds to face-on, the relationship between axis ratio $b/a$ and
inclination $\phi$ is given by the following formula.

$\frac{b}{a}=\frac{D}{\sqrt{\sin^{2}{\phi}+D^{2}\cos^{2}{\phi}}}$ (5)

We drew values of $\cos{\phi}$ from a uniform distribution between $[0,1]$
(consistent with a random viewing angle) and calculated axis ratio $b/a$ for
each inclination.

We corrected the measured U, V, and J magnitudes of stellar particles for dust
attenuation and then added the fluxes of particles within $2r_{1/2}$ to get
the overall colors for each subhalo.

###  3.3 Measuring Colors for 3D-HST Galaxies

To measure rest-frame colors for observed galaxies at cosmological distances,
we need a predicted spectrum conditioned on existing data. This means that
even observationally-inferred rest-frame UVJ-colors are model-dependent.

In our case, we are using spectra generated through FSPS during Prospector
fits of 3D-HST galaxies. Prospector uses nested sampling methods (Speagle,
2020) to efficiently probe the physical parameter space in order to create a
spectrum that best fits the observed photometry.

This process is complicated by degeneracies between parameters such as age,
metallicity, and dust attenuation (e.g., Conroy, 2013; Santini et al., 2015).
As changes in these parameters have similar effects on the generated spectrum,
we can have different parameter combinations with the same rest-frame colors.

In this work, we also compare dust-free TNG colors to the inferred rest-frame
UVJ colors when dust is removed to directly study properties of TNG galaxies
without the added layer of our dust model. For the measurement of dust-free
colors, the degeneracies between dust, age, and metallicities become
particularly important as we set the dust attenuation to zero.

In order to reduce the effects of degeneracies and generate more reliable
dust-free colors, we use the dust model from Paper I hierarchically to re-
weight the Prospector posterior samples of individual galaxies. This process
is called hierarchical shrinkage: we use what we have learned about dust
attenuation for the entire sample to help restrict the parameter space for
individual galaxies.

We now derive the weights used in the process described above. We consider a
single galaxy. The physical parameters that influence dust attenuation (i.e.,
stellar mass, SFR, stellar metallicity, redshift, and axis ratio, in our
model) are denoted as $\mathbf{w}$. Meanwhile, the parameters of our dust
model for all galaxies in the sample are called $\boldsymbol{\theta}$.

Next, we label the (fixed) parameters defining the Prospector priors as
$\boldsymbol{\alpha}$. While most priors are flat over the parameters of
interest, the exact limits of those priors belong to $\boldsymbol{\alpha}$. In
addition, there are informative priors for some variables, including diffuse
dust optical depth $\tau_{2}$. For $\tau_{2}$, the mean, standard deviation,
and limits of the truncated normal distribution serving as the prior are part
of $\boldsymbol{\alpha}$. We denote the set of all photometric observations
for the galaxy $\mathbf{X}$. In all equations below $p(\ldots)$ represents the
probability distribution of the quantity in the parentheses.

We can then write the following formula for the expected value of a function
$f$ for the galaxy. In our case, $f$ is the value of the dust attenuation
parameters $n$ and $\tau$ as a function of $\mathbf{w}$ given population
parameters $\boldsymbol{\theta}$. We show that $f$ depends on the particular
population parameters $\boldsymbol{\theta}$ using a subscript:
$f_{\boldsymbol{\theta}}$.

$\langle f_{\boldsymbol{\theta}}\rangle=\int
f_{\boldsymbol{\theta}}(\mathbf{w})p(\mathbf{w}|\mathbf{X},\boldsymbol{\theta})d\mathbf{w}$
(6)

In this equation, we have assumed that the galaxy observations $\mathbf{X}$ do
not depend directly on the population model parameters $\boldsymbol{\theta}$.
Next, from Bayes’ Theorem, we can write the following equation.

$p(\mathbf{w}|\mathbf{X},\boldsymbol{\Omega})=\frac{p(\mathbf{X}|\mathbf{w})p(\mathbf{w}|\boldsymbol{\Omega})}{p(\mathbf{X}|\boldsymbol{\Omega})}$
(7)

In Equation 7, we have made the reasonable assumption that $\mathbf{X}$ does
not depend directly on the interim priors $\boldsymbol{\alpha}$, nor, again,
on the population parameters $\boldsymbol{\theta}$, but only on the galaxy
properties $\mathbf{w}$. Let $\boldsymbol{\Omega}$ represent either quantity,
so this assumption is equivalent to
$p(\mathbf{X}|\boldsymbol{\Omega},\mathbf{w})=p(\mathbf{X}|\mathbf{w})$. Using
Equation 7, we can substitute for
$p(\mathbf{w}|\mathbf{X},\boldsymbol{\theta})$ and multiply the integrand in
Equation 6 by 1 in the following manner.

$\langle f_{\boldsymbol{\theta}}\rangle=\int
f_{\boldsymbol{\theta}}(\mathbf{w})\frac{p(\mathbf{X}|\mathbf{w})p(\mathbf{w}|\boldsymbol{\theta})}{p(\mathbf{X}|\boldsymbol{\theta})}\left[\frac{p(\mathbf{w}|\mathbf{X},\boldsymbol{\alpha})p(\mathbf{X}|\boldsymbol{\alpha})}{p(\mathbf{w}|\boldsymbol{\alpha})p(\mathbf{X}|\mathbf{w})}\right]d\mathbf{w}$
(8)

Therefore, we end up with the simplified result below.

$\langle f_{\boldsymbol{\theta}}\rangle=\int
f_{\boldsymbol{\theta}}(\mathbf{w})\frac{p(\mathbf{w}|\boldsymbol{\theta})p(\mathbf{X}|\boldsymbol{\alpha})}{p(\mathbf{w}|\boldsymbol{\alpha})p(\mathbf{X}|\boldsymbol{\theta})}p(\mathbf{w}|\mathbf{X},\boldsymbol{\alpha})d\mathbf{w}$
(9)

For a given set of parameters $\boldsymbol{\theta}$, the quantity
$\frac{p(\mathbf{X}|\boldsymbol{\alpha})}{p(\mathbf{X}|\boldsymbol{\theta})}$
is constant over $\mathbf{w}$ and therefore can be moved outside the integral.
We hereafter refer to this quantity as $1/Z_{\boldsymbol{\theta}}$. For
convenience, we label
$\frac{p(\mathbf{w}|\boldsymbol{\theta})}{p(\mathbf{w}|\boldsymbol{\alpha})}$
as $g_{\boldsymbol{\theta}}(\mathbf{w})$. We thus rewrite Equation 9 as the
following.

$\langle f_{\boldsymbol{\theta}}\rangle=\frac{1}{Z_{\boldsymbol{\theta}}}\int
f_{\boldsymbol{\theta}}(\mathbf{w})g_{\boldsymbol{\theta}}(\mathbf{w})p(\mathbf{w}|\mathbf{X},\boldsymbol{\alpha})d\mathbf{w}$
(10)

If we let the function $f$ be unity, the integral above must yield a value of
$1$. Therefore, we find that the normalization constant
$Z_{\boldsymbol{\theta}}$ is given by the following formula.

$Z_{\boldsymbol{\theta}}=\int
g_{\boldsymbol{\theta}}(\mathbf{w})p(\mathbf{w}|\mathbf{X},\boldsymbol{\alpha})d\mathbf{w}$
(11)

Since we have at our disposal posterior samples from Prospector fits of 3D-HST
galaxies, which are sampled from the distribution $\mathbf{w}_{n}\sim
p(\mathbf{w}|\mathbf{X},\boldsymbol{\alpha})$, we can use the Monte Carlo
integration technique to approximate the integrals above. If we have $N$
posterior samples for a given galaxy, we get the following result.

$\langle
f_{\boldsymbol{\theta}}\rangle\approx\frac{1/N\sum_{n=1}^{N}f_{\boldsymbol{\theta}}(\mathbf{w}_{n})g_{\boldsymbol{\theta}}(\mathbf{w}_{n})}{1/N\sum_{n=1}^{N}g_{\boldsymbol{\theta}}(\mathbf{w}_{n})}$
(12)

Equation 12 gives a straightforward method to measure the weight for any given
Prospector posterior sample and set of population model parameters.

Ultimately, we want to get the expected value of $f$ not just conditioned on
any one value of the population parameters $\boldsymbol{\theta}$, but rather
the full posterior distribution of $\boldsymbol{\theta}$ conditioned
simultaneously on the photometry from all $\sim 30,000$ 3D-HST galaxies which
we obtained in Paper I. In that case, we have another integral over
$d\boldsymbol{\theta}$ and can do another Monte Carlo integration using the
posterior samples from the population modeling efforts to find the effective
weights.

By applying the population-model-posterior-averaged weights to the Prospector
posterior samples, we measure values of $n$ and $\tau$ for each galaxy that
are less impacted by parameter correlations in the SED fitting process. Thus,
we are able to compute more reliable dust-attenuated and especially dust-free
UVJ colors.

In Section 4, we compare data to simulations in both the dust-free and with-
dust spaces (reverse and forward modeling, respectively). In order to
calculate dust-free 3D-HST colors, we use the equations above to reweight the
Prospector posterior samples for each galaxy by our hierarchical population
model, take the weighted average of the samples, set the dust attenuation
parameters to zero, and then run FSPS to get the resulting spectrum. For
comparisons with the simulated galaxies, we limit the redshift range to
$0.8\leq z\leq 1.2$ to better mirror the $z=1$ selection in TNG. We also
impose the mass-completeness limits in the 3D-HST sample on the TNG galaxies
(at $z=1$).

##  4 TNG100 vs 3D-HST UVJ Diagrams

We first compare dust-free colors of our TNG100 subhalos and Prospector dust-
free rest-frame colors conditioned on 3D-HST photometry. There is scatter
between the TNG dust-free colors measured by Donnari et al. (2019) and ours on
the order of $0.05$ magnitudes with negligible bias. These do not affect our
conclusions.

In Figure 2, we show the distribution of dust-free TNG100 UVJ colors in red
(whose calculation was described in §3.2). We juxtapose the distribution of
colors from the dust-free Prospector 3D-HST galaxies in blue (discussed in
§3.3).

The boundaries between quiescent and star forming are taken from Whitaker et
al. (2011), with modifications at the upper right corner inspired by Belli et
al. (2019), while the boundary between dusty and unobscured star forming
galaxies (DSFGs and USFGs) is from Fumagalli et al. (2014).

In the top and right panels, we show the 1-D projected distributions of $V-J$
and $U-V$, respectively, for both sources. The Prospector 3D-HST galaxies tend
to have slightly bluer colors, especially $V-J$, which we discuss briefly in
§5.

Figure 2: Dust-free UVJ diagram for TNG100 vs 3D-HST. The distribution of
dust-free colors from TNG100 at $z=1$—calculated for each galaxy by summing
the fluxes from every star particle within $2r_{1/2}$ to mimic the setup in
Donnari et al. (2019)–is shown in red. The distribution of dust-free colors of
3D-HST galaxies (whose calculation is described in §3.3) is shown in blue. The
68% (darker) and 95% (lighter) contours are shaded. See §3.2 for more details
on the procedure to calculate colors. The boundaries between quiescent and
star forming are taken from Whitaker et al. (2011), with modifications at the
upper right corner inspired by Belli et al. (2019), while the boundary between
dusty and unobscured is from Fumagalli et al. (2014). The observed galaxies
tend to have slightly bluer colors than the simulated galaxies, especially
$V-J$.

Next, we compare the colors with dust included to Prospector rest-frame colors
in Figure 3. Once again, the distribution of values from this work are shown
in red and that of observed galaxies in orange. In addition, we show the
distribution of colors for Donnari et al. (2019) using the theoretical
resolved dust model (C) from Nelson et al. (2018) in orange.

The differences between the TNG color distributions measured via two different
dust models is quite clear. In particular, we have a much more extended
quiescent strip at high $V-J$, a larger range of unobscured star-forming
galaxies, as well as a more populated DSFG region. Indeed, the Donnari et al.
(2019) distribution includes almost no DSFGs, while we can see they make up a
prominent part of the 3D-HST galaxies.

On the other hand, the Donnari et al. (2019) results have a wider distribution
of unobscured star forming galaxies (USFGs) than in this work when
$0.5\lesssim V-J\lesssim 1.3$, which is in better agreement with the 3D-HST
colors.

For both TNG efforts, the locus of USFGs is slightly shifted to redder colors
than in 3D-HST. However, this likely reflects the same difference in the dust-
free colors (Figure 2)—rather than limitations in the dust models. We discuss
this notion further in §5. Overall, our results are closer to the 3D-HST rest-
frame colors than are those from Donnari et al. (2019).

Figure 3: UVJ diagram for TNG100 vs 3D-HST (including dust). The distribution
of values from this work are shown in red, 3D-HST galaxies in blue, and those
from Donnari et al. (2019) in orange, with 68% (darker) and 95% (lighter)
contours shaded. Our dust model yields a longer quiescent strip of galaxies as
well as many more dusty star forming galaxies (DSFGs) than the resolved Nelson
et al. (2018) model used by Donnari et al. (2019), thereby giving it greater
resemblance to the observed UVJ distribution.

##  5 Discussion

We now examine the origins of the differences between the distribution of
dust-free colors for TNG100 $z=1$ and 3D-HST $0.8\leq z\leq 1.2$ galaxies. We
observe in Figure 2 that the observed galaxies have a slightly bluer locus of
points, with a greater effect in $V-J$ than $U-V$. In addition, they occupy a
larger part of parameter space.

In Figure 4, we show the distributions of average stellar population age
(left) and stellar metallicity (right) for the TNG simulated galaxies at $z=1$
(blue) and 3D-HST galaxies at $0.8\leq z\leq 1.2$ (red). We immediately notice
that the observed galaxies occupy wider distributions in age and metallicity
as well as systematically older and more metal-poor populations.

The larger spread in properties for the observed galaxies is presumably in
part due to the wider redshift range $0.8\leq z\leq 1.2$, a range chosen to
include a large enough sample ($9,000$) of Prospector galaxies, compared with
a precise $z=0.9973$ for TNG100. However, we find this to be a nearly
negligible effect. When we change the range from $0.8\leq z\leq 1.2$ to
$0.9\leq z\leq 1.1$, the means and standard deviations of the Prospector $U-V$
and $V-J$ distributions are shifted by less than $1\%$. This suggests that the
wider Prospector age and metallicity distributions, which are similarly
unchanged for the two redshift ranges, are the primary cause of the spread.

Another important consideration is the fact that the dust-free 3D-HST galaxy
colors are not “directly” observed but rather calculated with an SED fitting
code. Therefore, they are subject to modeling uncertainties and not just
galaxy variations. We find that this uncertainty is on the order of $0.1$
magnitudes. By employing our dust attenuation model to re-sample the
Prospector posterior distribution (§3.3), we get a $\sim 15\%$ reduction in
the standard deviation of the uncertainties over all galaxies: in other words,
we reduce the influence of outliers with very large uncertainties.

The systematically bluer 3D-HST colors, primarily $V-J$, can be explained by
the older, more metal-poor stellar populations. Lower metallicities yield
bluer colors from UV through NIR while older populations lead to redder
colors. Both of these effects are functions of wavelength, meaning the balance
between the opposing effects is different for $U-V$ and $V-J$.

Figure 4: Distributions of average stellar population age (left) and stellar
metallicity (right) for TNG simulated galaxies at $z=1$ (blue) vs 3D-HST
galaxies at $0.8\leq z\leq 1.2$. The wider distributions of observed galaxies
for age and metallicity may contribute to their larger distribution of UVJ
colors. Also, we observe that on average, TNG galaxies have younger stars and
higher metallicities. The combination of these two differences can help
explain the bluer 3D-HST colors. Higher metallicities signify redder colors,
while younger populations lead to bluer colors (opposite effect). The balance
of these effects as a function of wavelength can result in the disparity we
observe. We offer a more detailed explanation in the text.

In Figure 3, we compare the dust-included UVJ colors from this work, Donnari
et al. (2019), and 3D-HST rest-frame colors. As shown in §4, the TNG UVJ
diagram we calculate in this work is more similar to the 3D-HST observed
diagram than is that derived by Donnari et al. (2019). This is especially true
in the higher number densities of dusty star forming galaxies (DSFGs), which
are almost completely absent in the Donnari et al. (2019) diagram. The lack of
DSFGs has also been noted in the MUFASA simulation (Davé et al., 2017) and is
potentially related to resolution limitations.

Nevertheless, important differences remain between our TNG UVJ diagram and the
3D-HST rest-frame colors. For example, the distribution we derive is more
compact, with a smaller portion of color-color space being occupied by
unobscured star forming galaxies. In addition, it is slightly redder than the
observed distribution.

One of the reasons for this phenomenon is the intrinsic dust-free color
differences. The 3D-HST distribution is considerably wider than TNG100.
Naturally, that intrinsic thickness will be imprinted on the dust-included
colors. Similarly, the redder dust-free TNG colors help lead to redder dust-
included colors.

In addition to the dust-free considerations, we can attribute some of the
differences to contrasts between the population dust model and the underlying
data on which it was built. When we replace the individually fit dust
parameters (used to reproduce the rest-frame spectra) with our dust
attenuation model (not portrayed in the paper), the resulting distribution is
more compact and very slightly redder than the observed distribution,
suggesting that there are still some limitations in the dust model that affect
the UVJ distribution.

Perhaps the most important point is that in the 5-D linear interpolation
model, we are able to include only five grid points per dimension ($3,125$
model parameters total) to achieve convergence. Therefore, small-scale
variations in the dust attenuation function will be smoothed over. In
addition, we assume a single intrinsic scatter in the relationship for
simplicity; in reality, the scatter is likely a function of parameters.
Finally, there may be additional variables needed to fully explain the scatter
in attenuation curves. Nevertheless, we have demonstrated in both Paper I and
this work the prowess of the population dust attenuation model in its high
level of statistical accuracy and ability to provide a more reliable interface
between simulations and observations.

Analyzing how the differences between the UVJ colors derived in this work and
Donnari et al. (2019) vary as a function of physical properties like stellar
mass will help us understand the underlying differences in the dust models
used. In Figure 5, we show the difference in $U-V$ (blue) and $V-J$ (red) as a
function of stellar mass (top left), sSFR (top right), metallicity (bottom
left), and average age of stellar populations (bottom right). To remove traces
of differences in forward-modeling the TNG SEDs, we subtract the differences
in derived dust-free colors from the result. For example, for $U-V$, the
quantity plotted is $(U-V)_{\rm Don}-(U-V)_{\rm
Nag}-\left[(U-V)_{\rm{Don,~{}dust-free}}-(U-V)_{\rm{Nag,~{}dust-
free}}\right]$.

From Figure 5, we observe that for both $U-V$ and $V-J$, the color residuals
are systematically negative for high stellar mass, low sSFR, high metallicity,
and/or older galaxies. This equates to a prediction of bluer galaxies in that
part of parameter space from the resolved dust model of Nelson et al. (2018).
Comparing a resolved and unresolved dust model is difficult as the effects of
geometry are complex, but we can make a simple statement: if the resolved
model were replaced by an effective dust screen, the resulting attenuation
curve would either have a smaller optical depth and/or be flatter than our
empirical attenuation model for high-mass, high-metallicity, low-sSFR, and/or
old galaxies.

Figure 5: Galaxy-by-galaxy differences between dust-attenuated UVJ colors for
TNG100 galaxies derived in this work and Donnari et al. (2019) as a function
of stellar mass (top left), sSFR (top right), metallicity (bottom left), and
average age of stellar populations (bottom right). Very minor differences
between the dust-free UVJ colors are removed to create a fair baseline for the
comparisons. We see that the UVJ differences are systematically negative for
high stellar mass, low sSFR, high metallicity, and/or older galaxies. In other
words, Donnari et al. (2019) predicts bluer galaxies in that part of parameter
space. If the effective attenuation curve obtained by Nelson et al. (2018)
based on the simulation particle data were replaced by an effective dust
screen, the predicted attenuation curve would either have a lower optical
depth or be flatter than our empirical curve for such galaxies.

##  6 Conclusion

We have created a population Bayesian model of dust attenuation (Paper I). In
this paper, we detail its usage in theoretical scenarios and demonstrate its
success by applying it to TNG100 galaxies at $z=1$ and comparing the resulting
UVJ colors to both previous efforts to model the effects of dust in these
galaxies (Nelson et al., 2018; Donnari et al., 2019) and rest-frame colors of
3D-HST galaxies.

Overall, our population Bayesian dust attenuation model results in a more
observation-like UVJ diagram for TNG100 galaxies with stellar masses $M\geq
10^{9}\hbox{$M_{\odot}$}$ than previous efforts.

Furthermore, by using the population model to reweight posterior samples of
individual galaxy SED fits (hierarchical shrinkage), we are able to calculate
more reliable dust-free rest-frame colors for observed galaxies and thus
uncover more fundamental differences between TNG galaxies and the observed
universe. In particular, we find that on average, TNG galaxies have a 29%
younger and 28% narrower (stellar population) age distribution and a 45%
narrower and $0.28$ dex more metal-rich stellar metallicity distribution,
leading to a slightly redder dust-free distribution for TNG100 galaxies than
Prospector 3D-HST galaxies.

Dust plays a major role in determining the spectra of galaxies, but, as of
yet, its effects cannot be reproduced from first principles in galaxy
formation and evolution simulations. Compounding the issue, its effects on
integrated galaxy spectra are hard to distinguish from those of metallicity
and age, and thus confound comparisons with observations. By using a
statistically rigorous empirical approach to create our model, we are able to
provide a simple but accurate interface between theoretical (dust-free) models
and observed SEDs. Through the efforts described in this paper, we are able to
better understand the limitations of such models (e.g., issues in the
treatment of chemical evolution) without most of the complications from dust.

We would like to thank Sandro Tacchella for computing star formation
histories, metallicity histories, and stellar masses of TNG100 subhalos for
us. We would also like to thank Dylan Nelson for providing us with UVJ colors
used in Donnari et al. (2019) with and without dust for TNG100 subhalos. We
also thank Julianne Dalcanton for insightful conversations. 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.

## References

  * Arnouts et al. (2013) Arnouts, S., Le Floc’h, E., Chevallard, J., et al. 2013, A&A, 558, A67
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123
  * Belli et al. (2019) Belli, S., Newman, A. B., & Ellis, R. S. 2019, ApJ, 874, 17
  * Bessell (1990) Bessell, M. S. 1990, PASP, 102, 1181
  * Bouwens et al. (2012) Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2012, ApJ, 754, 83
  * Bowman et al. (2019) Bowman, W. P., Zeimann, G. R., Ciardullo, R., et al. 2019, ApJ, 875, 152
  * Brammer et al. (2012) Brammer, G. B., van Dokkum, P. G., Franx, M., et al. 2012, ApJS, 200, 13
  * 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
  * Byler et al. (2017) Byler, N., Dalcanton, J. J., Conroy, C., & Johnson, B. D. 2017, ApJ, 840, 44
  * 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
  * Cecchi et al. (2019) Cecchi, R., Bolzonella, M., Cimatti, A., & Girelli, G. 2019, ApJ, 880, L14
  * Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763
  * Charlot & Fall (2000) Charlot, S., & Fall, S. M. 2000, ApJ, 539, 718
  * 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
  * Davé et al. (2019) Davé, R., Anglés-Alcázar, D., Narayanan, D., et al. 2019, MNRAS, 486, 2827
  * Davé et al. (2017) Davé, R., Rafieferantsoa, M. H., & Thompson, R. J. 2017, MNRAS, 471, 1671
  * Donnari et al. (2019) Donnari, M., Pillepich, A., Nelson, D., et al. 2019, MNRAS, 485, 4817
  * Draine (2003) Draine, B. T. 2003, ARA&A, 41, 241
  * Finkelstein et al. (2012) Finkelstein, S. L., Papovich, C., Salmon, B., et al. 2012, ApJ, 756, 164
  * Forbes et al. (2014) Forbes, J. C., Krumholz, M. R., Burkert, A., & Dekel, A. 2014, MNRAS, 443, 168
  * Forbes et al. (2019) Forbes, J. C., Krumholz, M. R., & Speagle, J. S. 2019, MNRAS, 487, 3581
  * Foreman-Mackey et al. (2014a) Foreman-Mackey, D., Hogg, D. W., & Morton, T. D. 2014a, ApJ, 795, 64
  * Foreman-Mackey et al. (2014b) Foreman-Mackey, D., Sick, J., & Johnson, B. 2014b, Python-Fsps: Python Bindings To Fsps (V0.1.1), v0.1.1, Zenodo
  * Fumagalli et al. (2014) Fumagalli, M., Labbé, I., Patel, S. G., et al. 2014, ApJ, 796, 35
  * Galliano et al. (2018) Galliano, F., Galametz, M., & Jones, A. P. 2018, ARA&A, 56, 673
  * Genel et al. (2014) Genel, S., Vogelsberger, M., Springel, V., et al. 2014, MNRAS, 445, 175
  * Grogin et al. (2011) Grogin, N. A., Kocevski, D. D., Faber, S. M., et al. 2011, ApJS, 197, 35
  * Johnson et al. (2021) Johnson, B. D., Leja, J., Conroy, C., & Speagle, J. S. 2021, ApJS, 254, 22
  * 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. (2019a) Leja, J., Tacchella, S., & Conroy, C. 2019a, ApJ, 880, L9
  * Leja et al. (2019b) Leja, J., Johnson, B. D., Conroy, C., et al. 2019b, ApJ, 877, 140
  * Marinacci et al. (2018) Marinacci, F., Vogelsberger, M., Pakmor, R., et al. 2018, MNRAS, 480, 5113
  * Momcheva et al. (2016) Momcheva, I. G., Brammer, G. B., van Dokkum, P. G., et al. 2016, ApJS, 225, 27
  * Muzzin et al. (2013) Muzzin, A., Marchesini, D., Stefanon, M., et al. 2013, ApJ, 777, 18
  * Nagaraj et al. (2022) Nagaraj, G., Forbes, J. C., Leja, J., Foreman-Mackey, D., & Hayward, C. C. 2022, arXiv e-prints, arXiv:2202.05102
  * Naiman et al. (2018) Naiman, J. P., Pillepich, A., Springel, V., et al. 2018, MNRAS, 477, 1206
  * Nelson et al. (2018) Nelson, D., Pillepich, A., Springel, V., et al. 2018, MNRAS, 475, 624
  * 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
  * Oke (1974) Oke, J. B. 1974, ApJS, 27, 21
  * Patel et al. (2012) Patel, S. G., Holden, B. P., Kelson, D. D., et al. 2012, ApJ, 748, L27
  * Pillepich et al. (2018a) Pillepich, A., Springel, V., Nelson, D., et al. 2018a, MNRAS, 473, 4077
  * Pillepich et al. (2018b) Pillepich, A., Nelson, D., Hernquist, L., et al. 2018b, MNRAS, 475, 648
  * Planck Collaboration et al. (2016) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A13
  * 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
  * Skelton et al. (2014) Skelton, R. E., Whitaker, K. E., Momcheva, I. G., et al. 2014, ApJS, 214, 24
  * Skrutskie et al. (2006) Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163
  * Speagle (2020) Speagle, J. S. 2020, MNRAS, 493, 3132
  * Springel et al. (2018) Springel, V., Pakmor, R., Pillepich, A., et al. 2018, MNRAS, 475, 676
  * 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., Chang, Y.-Y., Bell, E. F., et al. 2014, ApJ, 792, L6
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261
  * Whitaker et al. (2012) Whitaker, K. E., Kriek, M., van Dokkum, P. G., et al. 2012, ApJ, 745, 179
  * Whitaker et al. (2011) Whitaker, K. E., Labbé, I., van Dokkum, P. G., et al. 2011, ApJ, 735, 86
  * Whitaker et al. (2017) Whitaker, K. E., Bezanson, R., van Dokkum, P. G., et al. 2017, ApJ, 838, 19
  * Wild et al. (2011) Wild, V., Charlot, S., Brinchmann, J., et al. 2011, MNRAS, 417, 1760
  * Williams et al. (2009) Williams, R. J., Quadri, R. F., Franx, M., van Dokkum, P., & Labbé, I. 2009, ApJ, 691, 1879
  * Wu et al. (2018) Wu, P.-F., van der Wel, A., Bezanson, R., et al. 2018, ApJ, 868, 37
  * Wuyts et al. (2007) Wuyts, S., Labbé, I., Franx, M., et al. 2007, ApJ, 655, 51
  * Zuckerman et al. (2021) Zuckerman, L. D., Belli, S., Leja, J., & Tacchella, S. 2021, ApJ, 922, L32