Empirical Dust Attenuation Model Leads to More Realistic UVJ Diagram for TNG100 Galaxies

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 ~ 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.


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. 1994Calzetti et al. , 2000Chevallard 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. 2014Forbes et al. , 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. 2017Leja et al. , 2019bLeja et al. , 2020) of a masscomplete set of ∼ 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.

USING THE POPULATION MODEL
Our population models of dust attenuation 1 , as described in detail in Paper I, describe both the twocomponent 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 (τ λ,1 ), which affects only stars less than 10 million years old, is modeled as τ λ,1 ∝ λ −1 . The diffuse dust attenuation (τ λ,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 (λ) below) with parameters for slope and 2175Å bump strength. Prospector fixes the bump strength E b based on results from Kriek & Conroy (2013).
The parameters fitted in the two-component model are n (related to slope), τ 2 , and τ 1 . In the effective dust attenuation model, the curve is parameterized by Equation 2, with parameters n eff and τ eff . In this work, we use the two-component dust attenuation model, in which τ 1 is a linear interpolation function of τ 2 while n and τ 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., τ 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.

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 U V J colors to distinguish quiescent and starforming 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).

Measuring Magnitudes for TNG100 Subhalos
To exhibit our population dust attenuation models' validity, we apply them to simulated galaxies from the IllustrisTNG project 2 (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(Leja et al. , 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 M 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 , 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. 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 − σ 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.
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(SFR/[M yr −1 ]) < 1, we simply drew values for b/a from a uniform distribution between 0 and 1.
For galaxies with log(SFR/[M yr −1 ]) ≥ 1, we employed the simple triaxial model from Zuckerman et al. (2021). The sources are treated as ellipsoids with axes parameters a , b = a , and c (thickness). We let the ratio of thickness c to the semi-major axis a be D ≡ c /a = 0.1. Then, for a given inclination φ, where 0 corresponds to face-on, the relationship between axis ratio b/a and inclination φ is given by the following formula.
We drew values of cos φ 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.

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 w. Meanwhile, the parameters of our dust model for all galaxies in the sample are called θ.
Next, we label the (fixed) parameters defining the Prospector priors as α. While most priors are flat over the parameters of interest, the exact limits of those priors belong to α. In addition, there are informative priors for some variables, including diffuse dust optical depth τ 2 . For τ 2 , the mean, standard deviation, and limits of the truncated normal distribution serving as the prior are part of α. We denote the set of all photometric observations for the galaxy X. In all equations below p(. . .) 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 τ as a function of w given population parameters θ.
We show that f depends on the particular population parameters θ using a subscript: f θ .
In this equation, we have assumed that the galaxy observations X do not depend directly on the population model parameters θ. Next, from Bayes' Theorem, we can write the following equation.
p(w|X, Ω) = p(X|w)p(w|Ω) p(X|Ω) In Equation 7, we have made the reasonable assumption that X does not depend directly on the interim priors α, nor, again, on the population parameters θ, but only on the galaxy properties w. Let Ω represent either quantity, so this assumption is equivalent to p(X|Ω, w) = p(X|w). Using Equation 7, we can substitute for p(w|X, θ) and multiply the integrand in Equation 6 by 1 in the following manner.
(8) Therefore, we end up with the simplified result below.
For a given set of parameters θ, the quantity p(X|α) p(X|θ) is constant over w and therefore can be moved outside the integral. We hereafter refer to this quantity as 1/Z θ . For convenience, we label p(w|θ) p(w|α) as g θ (w). We thus rewrite Equation 9 as the following.
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 θ is given by the following formula.
Since we have at our disposal posterior samples from Prospector fits of 3D-HST galaxies, which are sampled from the distribution w n ∼ p(w|X, α), 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.
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 θ, but rather the full posterior distribution of θ conditioned simultaneously on the photometry from all ∼ 30, 000 3D-HST galaxies which we obtained in Paper I. In that case, we have another integral over dθ 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 τ 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 ≤ z ≤ 1.2 to better mirror the z = 1 selection in TNG. We also impose the masscompleteness limits in the 3D-HST sample on the TNG galaxies (at z = 1).

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.
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 starforming 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 V − J 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 re- Figure 2. Dust-free UVJ diagram for TNG100 vs 3D-HST. The distribution of dust-free colors from TNG100 at z = 1calculated 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.
sults are closer to the 3D-HST rest-frame colors than are those from Donnari et al. (2019).

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 ≤ z ≤ 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 ≤ z ≤ 1.2 (red). We immediately notice that the observed galaxies occupy wider distri-butions 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 ≤ z ≤ 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 ≤ z ≤ 1.2 to 0.9 ≤ z ≤ 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. . 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.
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 ∼ 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.
In Figure 3, we compare the dust-included UVJ colors from this work, Donnari et al. (2019), and 3D-HST restframe 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 be- 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 ≤ z ≤ 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.
ing 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, smallscale 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 forwardmodeling the TNG SEDs, we subtract the differences in derived dust-free colors from the result. For example, for U − V , the quantity plotted is ( 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.

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 ≥ 10 9 M 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 metalrich 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.