Reproducing the UVJ Color Distribution of Star-forming Galaxies at 0.5

We analyze the distribution of rest-frame U-V and V-J colors for star-forming galaxies at 0.5<z<2.5. Using stellar population synthesis, stochastic star formation histories, and a simple prescription for the dust attenuation that accounts for the shape and inclination of galaxies, we construct a model for the distribution of galaxy colors. With only two free parameters, this model is able to reproduce the observed galaxy colors as a function of redshift and stellar mass remarkably well. Our analysis suggests that the wide range of dust attenuation values measured for star-forming galaxies at a given redshift and stellar mass is almost entirely due to the effect of inclination; if all galaxies were observed edge-on, they would show very similar dust attenuation. This result has important implications for the interpretation of dust attenuation measurements, the treatment of UV and IR luminosity, and the comparison between numerical simulations and observations.


INTRODUCTION
Much of our knowledge about the physical properties of distant galaxies derives from fitting models to multiband photometric observations (Conroy 2013, and references therein). By combining high-quality panchromatic photometry with modern fitting techniques it is now possible to study individual galaxies in great detail, and constrain realistic physical models with a large number of free parameters (e.g., Johnson et al. 2021).
A complementary approach is to consider a small number of filters or colors and study their distribution for a galaxy population. This is typically done using color-color diagrams; among these, the most popular one is the U V J diagram, i.e. rest-frame U − V versus V − J (Wuyts et al. 2007;Williams et al. 2009). The combination of these two colors is able to break the degeneracy between age and dust reddening, and is therefore particularly effective at separating quiescent and star-forming galaxies. However, the U V J diagram has been used for more than merely grouping galaxies into two categories. For example, U V J colors have been used to infer the star formation rate and dust attenuation for star-forming systems (Fang et al. 2018), and the stellar ages for quiescent systems (Belli et al. 2019). This is possible because galaxies follow scaling relations, so that their properties are tightly correlated and span only a small region of the available parameter space. In other words, the U V J diagram reveals more about the physical properties than what can be gathered by the analysis of the U − V and the V − J colors alone (Leja et al. 2019).
The regularity observed in the distribution of UVJ galaxy colors, however, is still poorly understood. For example, the most advanced simulations of galaxy formation still struggle to reproduce the observed color distribution of star-forming galaxies (e.g., Donnari et al. 2019;Akins et al. 2021). The goal of the present work is to develop a simple, empirical model that is able to reproduce the observed distribution of star-forming galaxies on the U V J diagram as a function of redshift and stellar mass. We use the 3D-HST catalog (Brammer et al. 2012;Skelton et al. 2014;Momcheva et al. 2016) of galaxies detected in Hubble Space Telescope (HST) imaging from the five fields of the CANDELS survey (Grogin et al. 2011;Koekemoer et al. 2011). We select galaxies with redshift (obtained by a joint fit to grism spectroscopy and photometry) in the range 0.5 < z < 2.5 and with stellar mass (derived from a fit to the photometry) in the range 9.5 < log(M * /M ) < 11. This selection yields a sample of 19553 galaxies that is more than 90% complete (Tal et al. 2014). We then remove galaxies with bad photometry (via the use phot flag provided in the catalog), bad fit (χ 2 > 5), or a bad measurement of the axis ratio q in the van der Wel et al. (2014) catalog. This leaves us with a final sample of 15073 galaxies.
For these galaxies we retrieve the rest-frame U − V and V − J colors from the 3D-HST catalog, which are derived from the observed photometry using the EAZY code (Brammer et al. 2008). We show the distribution of the sample on the U V J diagram in Figure 1, split into three stellar mass bins and four redshift bins, and colorcoded by axis ratio. The black line shows the definition of the quiescent box: (Muzzin et al. 2013;Belli et al. 2019). We use these criteria to select the star-forming population, which lies below the black line.

MODEL
In order to understand the distribution of galaxies on the U V J diagram we develop a simple Monte Carlo simulation (described below) which includes stochastic star formation histories (SFHs), stellar population synthesis, and a geometric model for the effect of dust attenuation.

Star Formation Histories
We adopt the stochastic model of Caplar & Tacchella (2019) to generate SFHs that are realistic but require minimal physical assumptions. In this model, galaxies are assumed to fluctuate about the star-forming main sequence, which describes the typical star formation rate as a function of stellar mass (Noeske et al. 2007). Since we are interested only in galaxy colors and not in the absolute fluxes, the normalization of the main sequence does not matter and we assume that the typical star formation rate for all galaxies, at all redshifts, is 1 M yr −1 .
This approach neglects the fact that the normalization of the main sequence declines with cosmic time. As we explain below, this assumption should not have a strong effect on the distribution of synthetic colors. Other work supports this model, for example Iyer et al. (2020) find that at fixed stellar mass, galaxies in the IllustrisTNG simulation have roughly a constant star formation rate (except at high stellar masses where the effect of quenching, which we do not consider in this work, is important).
In the Caplar & Tacchella (2019) model, a stochastic SFH is fully characterized by the power spectrum density (i.e., the power as a function of frequency) of its fluctuations around the main sequence. The power spectrum density is assumed to be constant on timescales longer than a break timescale τ x , below which the power spectrum follows a power law with exponent α = −2, corresponding to a damped random walk. Tacchella et al. (2020) show that this power spectrum density can be obtained by solving the gas regulator model (Lilly et al. 2013) when the gas accretion is a white noise process. Moreover, a broken power law with a break timescale of approximately 1 Gyr is consistent with the results of numerical simulations (Iyer et al. 2020). For this reason we adopt a fixed value of τ x = 1 Gyr.
We begin the Monte Carlo simulation by generating 100 stochastic SFHs with a time sampling of 1 Myr and a duration equal to the current age of the universe. This library of SFHs constitutes the core of our model.

Stellar Population Synthesis
To produce synthetic galaxy colors we rely on the Flexible Stellar Population Synthesis (FSPS) library (Conroy et al. 2009;Conroy & Gunn 2010). We assume a Kroupa (2001) initial mass function, zero dust attenuation, and we include nebular emission as described in Byler et al. (2017). Then we calculate the synthetic photometry in the rest-frame Johnson U , Johnson V , and 2MASS J filters (the same ones used for the observed galaxies) at each time step of our stochastic SFHs. The blue line in the first panel of Figure 1 shows an example of synthetic evolutionary track obtained in this way. The track features a high degree of variability, reflecting the fact that the SFH is, by construction, a stochastic process. The model galaxy reaches the bluest colors when experiencing a peak in star formation rate, while in times of relative quiescence its colors become redder.
Observations of galaxies, however, are unable to constrain the track of individual systems across cosmic time, but can only measure the color distribution for a galaxy population. We therefore build a synthetic population of galaxies by drawing 10 different random values of stellar metallicity (which has a smaller effect on the colors compared to the SFH), uniformly distributed over −0.5 < log Z/Z < 0.2, for each of the 100 stochastic SFHs. This gives a sample of 1000 synthetic galaxy tracks, each of which includes hundreds of points on the U V J diagram. The resulting two-dimensional distribution is shown in red in Figure 1. In each bin we only consider ages between zero and the age of the universe 9.5 < log(M*/M0) < 10 0.5 < z < 1 1 < z < 1.5 1.5 < z < 2 2 < z < 2. at that redshift. The model distributions are quite similar at all redshifts despite the substantial difference in the maximum stellar age considered. This is not surprising, since the photometry of star-forming galaxies is mainly driven by the younger, more luminous stars, and is not very sensitive to the past history. It is likely, then, that the simplifying assumption of a non-evolving main sequence has a negligible effect on our conclusions.

Dust Attenuation
The synthetic distributions in Figure 1 clearly do not provide a good fit to the data, and this is because the effect of dust reddening has so far been neglected. We illustrate the expected color shift caused by an attenuation of one magnitude in the V band (A V = 1) in the bottom left panel of Figure 1. The exact direction of this dust arrow depends on the dust attenuation law; we assume the widely used Calzetti et al. (2000) law. This attenuation law appears to produce a color shift along the direction connecting the synthetic colors to the observed ones. For this reason, and to keep the model as simple as possible, we assume the Calzetti law for all the synthetic galaxies, even though in principle some variation between galaxies is expected (Salim & Narayanan 2020).
Dust attenuation is clearly fundamental in shifting galaxies along the diagonal direction of the U V J diagram. However, Figure 1 also shows a strong trend of the observed axis ratio along the diagonal direction. In particular, among star-forming galaxies those that are found at the reddest edge of the UVJ distribution are mostly edge-on, a fact first noticed by Patel et al. (2012).This suggests that the dust attenuation in each galaxy strongly depends on the geometry of the system, consistent with previous work (e.g., Giovanelli et al. 1995).
Several studies have sought to determine the effect of inclination in a variety of ways, for example using the thickness of the Tully-Fisher relation (e.g., Verheijen 2001), galaxy SED fitting (e.g., Kauffmann et al. 2003), or by examining how the color of individual galaxies change as a function of inclination (e.g., Maller et al. 2009). Instead, in this work we choose to capture this correlation with a simple geometric model in which each galaxy is represented by an oblate spheroid with halfradii a, b, c with a = b > c, and where D = c/a is the galaxy thickness. We define the inclination angle θ so that a fully face-on galaxy has θ = 0 • and a fully edge-on galaxy has θ = 90 • . We then make the assumption that the dust is uniformly distributed in the galaxy, so that the optical depth as a function of inclination can be calculated from the geometry of the model: τ (θ) ∝ a (sin 2 θ + D −2 cos 2 θ) −1/2 . By assuming that the attenuation A V is proportional to the optical depth along the line-of-sight (and therefore neglecting the emission from stars that are mixed with the dust), we obtain the following relation: where A Vmax is the attenuation when the galaxy is viewed edge-on.
To illustrate the effect of galaxy geometry on the dust attenuation, in Figure 2 we show the results of our model for a thin disk (D = 0.2, top row) and for a thick disk (D = 0.8, bottom row), assuming A Vmax = 2. The first column illustrates the relation between A V and inclination given in Equation 1. As expected, the inclination effect is strong for thin galaxies and less important for thicker galaxies; in the limit of a spherical galaxy (D = 1) there would be no difference between face-on and edge-on attenuation. The central panels show the expected distribution of dust attenuation for an ensemble of identical galaxies observed from random lines of sight. Thin galaxies feature a peak around the minimum value of attenuation (which, in our model, is given by A V = D · A Vmax ), while thick galaxies experience a substantially smaller dynamic range in attenuation. Finally, the panels on the right show the effect of dust attenuation on the U V J colors. We conclude that, in this simple model, the galaxy shape (i.e. the thickness parameter D) and the dust content (parameterized with A Vmax ) play a fundamental role in determining the distribution of a galaxy population on the U V J diagram.
Our geometric model of dust attenuation is substantially more flexible than the often-used slab model (where τ ∝ 1/ cos θ), while at the same time being conceptually simpler and requiring less assumptions than some of the more advanced models (e.g., Chevallard et al. 2013). The key feature of our approach is the presence of a free parameter describing the galaxy thickness; see Padilla & Strauss (2008) for a similar approach.

Monte Carlo Modeling of Galaxy Populations
Central to our analysis is the application of the geometric model to the synthetic U V J colors: this requires the parameters θ, D, and A Vmax for each synthetic track. Since galaxies are randomly oriented in space, we draw the inclination θ from a random distribution uniform in cos θ. For the galaxy thickness D we use the results of van der Wel et al. (2014), who derived the intrinsic shape of star-forming galaxies from the distribution of axis ratios measured for the 3D-HST sam-

Rest-frame V-J color
Rest-frame U-V color ple. They provide the mean and standard deviation for the distribution of D, assumed to be Gaussian, in each redshift and stellar mass bin. We use these distributions to draw a random value of D for each synthetic track. Finally, we assume that the A Vmax parameter of all the galaxies in a given bin follows a Gaussian distribution with free mean and dispersion values, which represent the only two free parameters of our Monte Carlo model.

Model Fitting
For each bin in mass and redshift, we divide the U V J plane into n cells and calculate the observed probability distribution p i (where i = 1, 2, ..., n is the cell index) representing the fraction of observed galaxies that lie in each cell. Then we generate a grid of Monte Carlo models by varying the two free parameters: A Vmax = 0, 0.3, 0.6, 0.9, 1, 1.5, 2, 2.5 and σ max = 0.02, 0.06, 0.1, 0.15, 0.2, 0.3, 0.4, 0.5, 0.6, controlling respectively the mean and dispersion of the maximum (i.e., edge-on) dust attenuation. At each point of the model grid we calculate the model probability distribution q i (A Vmax , σ max ) using the same U V J cells as for the observations. We also include a small Gaussian scatter of width 0.1 magnitude that represents the typical observational uncertainty on the colors. We then compare the data and the model using the K-L divergence (Kullback & Leibler 1951): The K-L divergence is infinity for cells where no model galaxies, but at least one observed galaxy is present. In principle this eliminates models that are clearly inconsistent with the data, but we need some flexibility to account for systematic uncertainties (such as photometric redshift errors). We therefore assign, in each model, a small probability to all U V J cells that represents the rare possibility of substantial outliers. We can then calculate the total K-L divergence between the data and a specific model by integrating over all the U V J cells: For a large number of galaxies N , where L is the likelihood of obtaining the data if the model is correct (see, e.g., Shlens 2014). Thus, the model that minimizes the K-L divergence is also the model that maximizes the likelihood, and we call this the "best fit".  We show the results of the fit in Figure 3. In each redshift and stellar mass bin the observed galaxies are shown in blue and the best-fit model in red. Considering that our simple model has only two free parameters, it is able to reproduce the data remarkably well across the mass and redshift range. The best-fit values of the model parameters are listed in each panel, and are also plotted in the top panel of Figure 4 as a function of redshift, with the vertical bars representing the scatter σ max . We find that the scatter is low compared to A Vmax , particularly for systems with intermediate and high stellar masses. This means that, at fixed mass and redshift, all starforming galaxies have roughly the same intrinsic amount of dust, and the wide observed range in A V and color is mostly due to inclination, confirming the result of Patel et al. (2012). We also find that the dust content of galaxies grows with cosmic time and stellar mass, in agreement with previous studies (e.g., Fang et al. 2018).

Dust Attenuation across Mass and Redshift
The bottom two panels of Figure 4 show the redshift trend for the galaxy triaxiality and thickness measured by van der Wel et al. (2014). Most galaxies are relatively thin disks with D ∼ 0.25 and a small scatter. However, low-mass systems at high redshift are also highly triaxial, which suggests that our oblate spheroid model (which has zero triaxiality) is not appropriate in this regime. Interestingly, Figure 3 shows that the model does not fit this population of low-mass, high-redshift galaxies well, with the best-fit model extending on the diagonal rather than the horizontal direction of the U V J diagram. However, the fact that the best-fit dust attenuation is negligible for this population may also indicate a problem with the assumed SFH, such as the strength or timescale of its stochastic fluctuations. Additionally, our model appears insufficient to reproduce the observations at the opposite end of the parameter space, i.e. for high-mass, low-redshift galaxies. In this case the model produces too many galaxies with low A V . This is likely connected to the simplistic assumptions regarding the spatial distribution of stars in our model. If one were to consider the contribution from stars that are uniformly mixed with the dust, then the model would need a much higher A Vmax to reproduce the dustier systems, and as a consequence the face-on galaxies would also appear more attenuated.

Dust attenuation versus dust mass
The main result of this work is that the observed variation in UVJ colors of star-forming systems at a given stellar mass can be attributed almost entirely to variation in viewing angle rather than variation in intrinsic dust mass. This brings several implications for galaxy studies. Most importantly, we conclude that the U V J diagram is not a good tool for identifying "dusty" galaxies. Colors that appear reddened might only indicate higher inclinations (or thicker galaxies, as illustrated in Figure 2), and not intrinsically higher dust content. On the other hand, if color is a useful indicator of inclination, our model could be used to select galaxies of a given inclination solely based on the photometry, instead of more difficult morphological measurements.

The interpretation of IR and UV emission
Another important consequence of our results is related to the use of L IR /L UV (i.e., the ratio of IR and UV luminosity) to estimate the proportion of star formation that is obscured. This method is based on the assumption that L IR represents the obscured star formation while L UV traces the unobscured star formation. However, in our model the dust attenuation in the UV is strongly dependent on inclination, while the IR luminosity is not affected by inclination since dust attenuation in the IR is negligible. This means that the ratio L IR /L UV at a given stellar mass is mostly a measurement of inclination, and not of the intrinsic fraction of obscured star formation. For example, adopting Galaxy properties are derived by fits to the UV-to-IR photometry (Leja et al. 2020). Only galaxies with at least two Herschel detections are shown. D = 0.25, A Vmax = 2.5, and the Calzetti law, our model predicts a difference in the edge-on versus face-on attenuation of 4 magnitudes in the UV, meaning that the L IR /L UV ratio can vary by a factor of ∼ 40 solely due to inclination effects. We test this prediction of our model using a small sample of galaxies with available far-IR data. We select galaxies from the COSMOS2015 catalog (Laigle et al. 2016) with 0.5 < z < 0.8 and 10.5 < log(M * /M ) < 11 that are detected in Herschel. Figure 5 shows this sample on the U V J diagram, color-coded by L IR /M * (left) and L IR /L UV (right). Moving along the diagonal direction towards redder, more dust-attenuated systems, the L IR /L UV ratio varies by almost two orders of magnitude, while the total infrared luminosity (normalized by stellar mass) is virtually constant. This finding is consistent with our main result: at a given stellar mass, a galaxy's position along the dust vector in U V J space is determined primarily by galaxy orientation and does not correlate with the true amount of dust in the galaxy (see also Arnouts et al. 2013 for a similar analysis).
The fact that inclination has dramatically different effects on L IR and L UV should be taken into account in many other contexts; for example, when deriving star formation rates using the "UV+IR" method. Similarly, when imposing energy balance in photometric fits that include both UV and IR observations, if the inclination angle can be estimated from imaging then it is possible, in principle, to calculate inclination-dependent corrections (see Doore et al. 2021 for a first attempt in this direction).

Understanding the U V J colors of galaxies
Advanced cosmological simulations are still unable to reproduce the full extent of the U V J colors of starforming galaxies, particularly at the red end (Donnari et al. 2019;Akins et al. 2021). Our results can greatly help in understanding this discrepancy, and suggest that both the line-of-sight used to calculate the colors and the galaxy shapes are crucial for a meaningful comparison. Moreover, the fact that we are able to reproduce most of the observations using stochastic SFHs suggests that galaxy colors have limited constraining power on the history of star-forming galaxies, in agreement with the results of recent studies (Chaves-Montero & Hearin 2020; Hahn et al. 2021).
Our work can be considered the first step towards a fully quantitative understanding of the distribution of galaxies on the U V J diagram. Future developments should include a more realistic treatment of SFH, stellar metallicity, dust extinction law, and galaxy geometry, and should be extended to the quiescent population. Such a comprehensive empirical model of how galaxy colors evolve as a function of mass and redshift will provide new constraints on the formation, growth, and quenching of galaxies.