The MOSDEF Survey: Implications of the Lack of Evolution in the Dust Attenuation-Mass Relation to z~2

We investigate the relationship between dust attenuation and stellar mass ($M_*$) in star-forming galaxies over cosmic time. For this analysis, we compare measurements from the MOSFIRE Deep Evolution Field (MOSDEF) survey at $z\sim2.3$ and the Sloan Digital Sky Survey (SDSS) at $z\sim0$, augmenting the latter optical dataset with both UV Galaxy Evolution Explorer (GALEX) and mid-infrared Wide-field Infrared Survey Explorer (WISE) photometry from the GALEX-SDSS-WISE Catalog. We quantify dust attenuation using both spectroscopic measurements of H$\alpha$ and H$\beta$ emission lines, and photometric measurements of the rest-UV stellar continuum. The H$\alpha$/H$\beta$ ratio is used to determine the magnitude of attenuation at the wavelength of H$\alpha$, $A_{{\rm H}\alpha}$. Rest-UV colors and spectral-energy-distribution fitting are used to estimate $A_{1600}$, the magnitude of attenuation at a rest wavelength of 1600\AA. As in previous work, we find a lack of significant evolution in the relation between dust attenuation and $M_*$ over the redshift range $z\sim0$ to $z\sim2.3$. Folding in the latest estimates of the evolution of $M_{{\rm dust}}$, $({M_{{\rm dust}}}/{M_{{\rm gas}}})$, and gas surface density at fixed $M_*$, we find that the expected $M_{{\rm dust}}$ and dust mass surface density are both significantly higher at $z\sim2.3$ than at $z\sim0$. These differences appear at odds with the lack of evolution in dust attenuation. To explain the striking constancy in attenuation vs. $M_*$, it is essential to determine the relationship between metallicity and $({M_{{\rm dust}}}/{M_{{\rm gas}}})$, the dust mass absorption coefficient, and dust geometry, and the evolution of these relations and quantities from $z\sim0$ to $z\sim2.3$.


INTRODUCTION
Tracing the effects of dust attenuation on starlight is crucial for obtaining a complete census of star formation over cosmic time (e.g., Madau & Dickinson 2014).There are many different methods for quantifying the role played by dust in starforming galaxies over a wide range in redshift.These rely on multi-wavelength data probing the ratio of far-infrared (far-IR) to UV emission -i.e., re-radiated vs. unobscured starlight; and various measures of UV/optical reddening and/or attenuation.Measuring the dust content of star-forming galaxies over cosmic time in systems spanning a range of stellar and gas masses also places constraints on models for the formation and destruction of dust grains in the interstellar medium (ISM) (e.g., Popping et al. 2017).
One particularly striking observation regarding dust in star-forming galaxies is that the relationship between dust attenuation and stellar mass (M * ) does not significantly evolve between z ∼ 0 and z ∼ 2 (and perhaps to even higher redshift).Here dust attenuation has been quantified as the ratio of far-IR to UV star-formation rates (SFRs) or luminosities, which is also called "IRX" (e.g., Meurer et al. 1999;Heinis et al. 2014;Bouwens et al. 2016;Bourne et al. 2017); the magnitude of far-UV (i.e., 1600Å) attenuation, or A 1600 (e.g., McLure et al. 2018;Pannella et al. 2015); the fraction of star formation that is obscured, f obscured (Whitaker et al. 2017), and the nebular attenuation based on the Balmer decrement (i.e., Hα/Hβ ratio; Kashino et al. 2013;Domínguez et al. 2013;Price et al. 2014).There is less consensus regarding the form of the attenuation vs. M * relation at z > 3, with some evidence that it may evolve towards lower attenuation at fixed M * (e.g., Fudamoto et al. 2020).However, at least out to z ∼ 2, multiple results sug-gest a constant relation between dust attenuation and M * in mass-complete samples (Whitaker et al. 2017;McLure et al. 2018).
The degree of dust attenuation in a galaxy reflects multiple key features of its ISM.First, there is the total dust content, M dust , and its relationship with the gas content of the galaxy M gas .This dust content is intimately connected with the degree of metal enrichment in the galaxy, given that dust grains form from heavy elements (Rémy-Ruyer et al. 2014;De Vis et al. 2019).However, dust attenuation reflects not only the total dust content, M dust , but also its spatial distribution, which can be quantified in the most simplistic manner in terms of a characteristic radius, r dust (i.e., a dustcontinuum half-light radius).In more detail, the non-uniformity of the spatial distribution of dust must also be taken into account (e.g., Witt & Gordon 2000;Charlot & Fall 2000;Seon & Draine 2016).Finally, there is the question of the very properties of the dust grains, including their chemical composition, size distribution and morphologies, which determines the relationship between M dust and opacity.
Thus far, analyses of the (lack of) evolution in the attenuation vs. M * relation have not incorporated what is known about change in these other ISM components: M gas , metallicity (12 + log(O/H), and M dust .Yet, they must be considered in order to gain a full understanding of why a measurement of M * at z ∼ 0 to z ∼ 2 is so determinative of the degree of dust attenuation.Here, we analyze z ∼ 2.3 dust attenuation based on rest-optical spectroscopic measurements of the Hα/Hβ Balmer decrement and rest-UV continuum measures of dust reddening.We also fold in independent results on the evolution of galaxy metallicities and gas and dust content, in order to gain a complete picture of the ISM of star-forming galaxies over the past ∼ 10 billion years.In §2, we describe our samples and observations.In §3, we present results on the observed relationship between dust attenuation and stellar mass at both z ∼ 2.3 and z ∼ 0. In §4, we discuss the surprising implications of the lack of strong evolution in the attenuation vs. mass relation.Throughout, we adopt cosmological parameters of H 0 = 70 km s −1 Mpc −1 , Ω m = 0.30, and Ω Λ = 0.7, and a Chabrier (2003) IMF.

MOSDEF
The analysis presented here is based on data from the MOSFIRE Deep Evolution Field (MOS-DEF) survey, a large observing program using the Multi-Object Spectrometer for Infrared Exploration (MOSFIRE; McLean et al. 2012) on the 10 m Keck I telescope.As described in Kriek et al. (2015), with the MOSDEF survey we obtained rest-optical spectra for a sample of ∼ 1500 galaxies within three distinct redshift intervals spanning 1.4 ≤ z ≤ 3.8.These intervals are 1.37 ≤ z ≤ 1.70, 2.09 ≤ z ≤ 2.61, and 2.95 ≤ z ≤ 3.80, where the strongest rest-optical emission lines can be observed within windows of atmospheric transmission.MOSDEF targets fell in the COSMOS, GOODS-N, AEGIS, GOODS-S, and UDS fields, in regions covered by the CANDELS and 3D-HST surveys (Grogin et al. 2011;Koekemoer et al. 2011;Momcheva et al. 2016).These fields feature extensive multi-wavelength datasets spanning the electromagnetic spectrum, which can be used to infer a wide range of galaxy properties.The data used for fitting the spectral energy distributions (SEDs) of MOSDEF galaxies have been cataloged by the 3D-HST survey (Skelton et al. 2014), and include optical and near-infrared ground-based and Hubble Space Telescope photometry, as well as Spitzer/IRAC mid-infrared measurements.MOS-DEF targets are selected based on existing photometric or spectroscopic redshifts, with a magnitude limit in the rest-optical (observed H band).This limit is H AB =24, 24, and 25, respectively, for the lowest, middle, and highest redshift interval of the survey.
Here we focus on MOSDEF star-forming galaxies within the central target redshift range, i.e., 2.09 ≤ z ≤ 2.61.Our z ∼ 2 sample is very similar to the one analyzed in Sanders et al. (2021), with the additional constraint of ≥ 3σ detections of both Hα and Hβ line fluxes.Each galaxy has a robust estimate of nebular oxygen abundance (12 + log(O/H)) based on the subset of detected strong nebular emission lines drawn from [OII]λλ3726, 3729, [NeIII]λ3869, Hβ, and [OIII]λ5007 (at the very least [OII]λλ3726, 3729, Hβ, and [OIII]λ5007), as described in Sanders et al. (2021), and a measure of nebular dust attenuation (A Hα ) based on the ratio Hα/Hβ and assuming the Cardelli et al. (1989) extinction law.We also estimated stellar masses by modeling the multi-wavelength photometric spectral energy distributions (SEDs) cataloged by the 3D-HST team (Skelton et al. 2014), where the near-infrared photometry was corrected for the contribution of strong nebular emission lines (Sanders et al. 2021).For SED modeling, we used the program FAST (Kriek et al. 2009), assuming the stellar population synthesis models of Conroy et al. (2009), a Chabrier (2003) IMF, a Calzetti et al. (2000) dust law, and delayed-τ star-formation histories, where SFR(t) ∝ t exp(−t/τ ).Here, t is the time since the onset of star formation and τ is the characteristic star-formation timescale.We note that Balmer emission-line fluxes were corrected for the underlying stellar absorption implied by the best-fit stellar population model, with typical Balmer absorption corrections of ∼ 1% for Hα and ∼ 7% for Hβ.Finally, AGNs were identified and removed based on their X-ray and infrared properties, as well as those with log([NII]/Hα) > −0.3 (Coil et al. 2015;Azadi et al. 2017;Leung et al. 2019).In total, our MOSDEF sample includes 210 galaxies with a median redshift of z med = 2.28, a median stellar mass of log(M * /M ) med = 9.88, and a median dust-corrected Hα-based SFR of SFR med = 29M yr −1 .These median properties are very well-matched to those of the z ∼ 2.3 sample analyzed in Sanders et al. (2021).
As an additional measure of dust attenuation, which probes the stellar continuum, we estimated the UV slope, β, directly from broadband photometry.β is calculated by fitting a power law of the form f λ ∝ λ β to the photometric bands spanning the rest-wavelength range 1268 − 2580Å.This fit is typically determined based on 4−5 bands for galaxies in all fields except COSMOS, where the fit is typically based on 20 bands.The values of β were translated into estimates of rest-UV continuum attenuation (A 1600 ) for our sample based on a few different prescriptions.Following Reddy et al. (2018) and assuming an intrinsic stellar population from the Binary Population and Spectral Synthesis (BPASS) code (Eldridge et al. 2017), including binaries, an upper mass cut-off of 300M , Z * = 0.02, and log(age/yr) = 8.0, we find the relation A 1600 = 2.13 × β + 5.04 (1) for a Calzetti et al. (2000) dust attenuation law, and for an SMC extinction law (Gordon et al. 2003).We note that the relationship above for the Calzetti et al. (2000) law is based on assuming an intrinsically bluer UV slope of β int = −2.37 for z ∼ 2 star-forming galaxies than that found in earlier work for z ∼ 0 starbursts.Specifically, in Meurer et al. (1999), the relationship between A 1600 and β (A 1600 = 1.99 × β + 4.43) assumes an intrinsic UV slope of β int = −2.23.We also note that if BPASS models like the ones used to derive equations (1) and (2) are adopted to estimate stellar masses, we obtain results extremely consistent with those based on FAST models.The same holds if we adopt Bruzual & Charlot (2003), constant star formation, solar-metallicity models to estimate stellar masses.
Recent results from the MOSDEF survey (Shivaei et al. 2020) suggest that a Calzetti et al. (2000)type curve is appropriate at metallicities of 12 + log(O/H) ≥ 8.5, while one resembling the SMC curve seems to apply at 12 + log(O/H) < 8.5.These results echo previous evidence from Reddy et al. (2010) and Reddy et al. (2012) that an SMC curve best describes the youngest (age<100 Myr) systems among a large sample of UV-selected galaxies at z ∼ 2. Other recent work (e.g., McLure et al. 2018) presents evidence that the Calzetti et al. (2000) applies over a wide range of stellar masses (log(M * /M ) ≥ 9.75), and, correspondingly, metallicities.Accordingly, our two favored methods for translating β into A 1600 values for MOSDEF galaxies consist of (1) assuming the Calzetti et al. (2000) curve for the entire sample; (2) assuming the Calzetti et al. (2000) curve at 12 + log(O/H) ≥ 8.5 and the SMC curve at 12 + log(O/H) < 8.5.

SDSS Comparison Sample
In order to perform an evolutionary comparison, we selected a sample of local galaxies from the Sloan Digital Sky Survey (SDSS) Data release 7 (DR7; Abazajian et al. 2009).Stellar masses and emission-line measurements corrected for stellar absorption were drawn from the MPA-JHU catalog of measurements for DR7 1 .In this catalog, SDSS stellar masses were estimated by fitting a grid of Bruzual & Charlot (2003) models spanning a wide range in star-formation histories to u, g, r, i, and z emission-line-corrected photometry.We restricted the SDSS sample to galaxies at 0.04 ≤ z ≤ 0.10 to reduce aperture effects, and, following Andrews & Martini (2013), required 5σ detections for [OII]λλ3726, 3729, Hβ, Hα, and [NII]λ6584, and a 3σ detection for [OIII]λ5007.We also removed galaxies satisfying the optical emission-line AGN criterion of Kauffmann et al. (2003), yielding a z ∼ 0 comparison sample of 73,492 galaxies with log(M * /M ) med = 9.85.The SDSS sample is very well-matched to the z ∼ 2.3 MOSDEF sample in terms of median stellar mass.However, the median SFR for the SDSS sample is SFR med = 1.3M yr −1 , i.e., a factor of ∼ 20 lower, which reflects the evolution of the starforming main sequence (e.g., Förster Schreiber & Wuyts 2020, and references therein).Metallicities for SDSS galaxies were estimated using the calibrations presented in Figure 3  In order to compare measures of rest-UV attenuation, we used the GALEX-SDSS-WISE Legacy catalog (GSWLC) presented in Salim et al. (2016).As described in Salim et al. (2016), UV/optical galaxy SEDs including Galaxy Evolution Explorer (GALEX) and SDSS photometry were modeled using the CIGALE code (Boquien et al. 2019), including two-component exponentially-declining star-formation histories generated using Bruzual & Charlot (2003), and a range of dust attenuation laws with varying UV slopes and UV-bump strengths.Constraints on the preferred dust law were obtained by forcing agreement between SEDbased and mid-IR (from the Wide-field Infrared Survey Explorer; WISE) dust luminosities.One of the key outputs of the SED modeling using the GSWLC is the rest-UV attenuation, A 1600 .GALEX near-UV and far-UV coverage is available for 30,204 galaxies in our SDSS comparison sample.We have confirmed that this GALEX subsam-ple is representative of the larger SDSS-only sample, in terms of the relationship between Balmer line ratios and stellar mass (see Section 3).

RESULTS
The relationship between attenuation and stellar mass has been considered using several different measures of the effects of dust.These include UV attenuation, A 1600 (McLure et al. 2018), restoptical continuum and line attenuation, A V and A Hα (Cullen et al. 2018;Garn & Best 2010), the ratio of FIR to UV luminosity, "IRX" (Heinis et al. 2014;Bourne et al. 2017), the fraction of star formation that is obscured, f obscured (Whitaker et al. 2017), and the Balmer decrement (Kashino et al. 2013;Domínguez et al. 2013;Price et al. 2014).Here we analyze the attenuation of both rest-optical lines and rest-UV continuum as a function of stellar mass, and compare these relations for star-forming galaxies at z ∼ 0 and z ∼ 2.3.
First, we consider attenuation estimated from the ratio of Balmer emission lines, Hα/Hβ.Figure 1 (left) shows that there is a significant correlation between Hα/Hβ and M * , such that higher Hα/Hβ ratios are associated with higher M * .This correlation applies to both the z ∼ 0 SDSS and z ∼ 2.3 MOSDEF samples.We also plot the running median of Hα/Hβ in bins of M * .The z ∼ 2.3 running median relation between Hα/Hβ and M * is entirely consistent with that of the z ∼ 0 sample.The Hα/Hβ ratio can be translated into the nebular attenuation at the wavelength of Hα using a specific dust extinction curve.For star-forming galaxies in the local universe, the Milky Way curve of Cardelli et al. (1989) is typically used to interpret nebular reddening.Recently, Reddy et al. (2020) demonstrated that the Cardelli et al. (1989) curve is also appropriate for dust-correcting nebular emission lines in MOSDEF galaxies at 1.4 ≤ z ≤ 2.6.Using the Cardelli et al. (1989) curve for both z ∼ 0 and z ∼ 2.3 samples, we cast the Hα/Hβ vs. M * plot in terms of A Hα , the magnitude of nebular attenuation at the wavelength of Hα (Figure 1, right).As in the case of Hα/Hβ, the running medians of A Hα in bins of M * are also consistent between z ∼ 0 and z ∼ 2.3.We also note that the z ∼ 0 relationship presented here between A Hα and M * is entirely consistent with that of Garn & Best (2010), when the same nebular attenuation law is assumed.Finally, the stacked z ∼ 1.6 A Hα measurements of Kashino et al. (2013), in three bins of M * , are consistent with both our z ∼ 2.3 and the SDSS z ∼ 0 median relations in Figure 1, when the same stellar IMF and nebular attenuation law are assumed.
Dust attenuation is also commonly quantified in terms of its effect on the rest-UV stellar continuum.Such estimates of dust attenuation can be based on measurements of the UV slope, β, or else from SED fitting over a wider wavelength range.As described in Section 2.1, for the z ∼ 2.3 MOS-DEF sample, we use measurements of β and either the Calzetti et al. (2000) dust curve at all metallicities (and masses), or else a metallicity-dependent dust curve (Calzetti et al. 2000 at high metallicity and SMC at low metallicity) to estimate A 1600 .For the z ∼ 0 SDSS sample, A 1600 is inferred from UV/optical SED fitting and energy balance considerations.
Figure 2 shows the relationship between A 1600 and M * .For the sake of simplifying the plot while still conveying key information, we do not display individual z ∼ 2.3 datapoints under different assumptions regarding the dust curve, but rather only the corresponding running medians.Over the full range in stellar mass, the z ∼ 2.3 and z ∼ 0 medians in A 1600 are consistent within the errors when the Calzetti et al. (2000) curve is exclusively assumed (here z ∼ 2.3 points are connected by a solid curve).In this case, the z ∼ 2.3 median A 1600 value fall slightly higher than their z ∼ 0 counterparts at the same stellar mass.Our results are consistent with those of McLure et al. (2018), equation ( 17), in the highest-and lowest-mass bins, though, on average those authors find slightly higher A 1600 at fixed mass than we do here.If we use a metallicitydependent dust law to convert β to A 1600 (here z ∼ 2.3 points are connected by a dashed curve), Here attenuation is quantified as A 1600 , the magnitude of attenuation at a rest wavelength of 1600Å.As in Figure 1, the local SDSS sample is indicated as a grey histogram, and its running median A 1600 estimated in bins of M * is shown with red, connected symbols.A 1600 is determined for SDSS galaxies as described in Salim et al. (2016).The running median A 1600 for the z ∼ 2.3 MOSDEF sample in bins of M * is shown in turquoise, using three different prescriptions to translate the measured UV slope value, β, to A 1600 .The solid curve shows A 1600 based on assuming a Calzetti et al. (2000) dust attenuation curve; the dashed curve shows A 1600 assuming a Calzetti et al. (2000) curve at 12 + log(O/H) ≥ 8.5 and an SMC curve at 12 + log(O/H) < 8.5; finally, the dotted curve shows A 1600 assuming an SMC curve for the entire MOSDEF z ∼ 2.3 sample (a scenario not favored by the MOSDEF data, but included for completeness).As in Figure 1, median error bars are shown for the z ∼ 2.3 sample, while those for the SDSS sample are smaller than the symbols.
the z ∼ 2.3 A 1600 median values fall below the corresponding datapoints at z ∼ 0 and fixed stellar mass.The z ∼ 2.3 medians are consistent with those at z ∼ 0 down to a stellar mass of ∼ 10 10 M (at the ∼ 0.4σ and 2.5σ level, respectively, for the highest-and second-highest z ∼ 2.3 stellar mass bins).In the two lower-mass bins, the z ∼ 2.3 me-Figure 3. UV Attenuation (A 1600 ) vs. UV slope β for local SDSS galaxies.The SDSS sample is indicated as a grey histogram, and its running median A 1600 estimated in bins of β is shown with red, connected symbols.Also shown as bracketing cases are the expressions relating A 1600 and β from equations ( 1) and ( 2).The solid, turquoise line indicates the A 1600 -β relation assuming the Calzetti et al. (2000) dust curve, while the dotted turquoise line shows the A 1600 -β relation assuming the SMC dust curve.In addition, we show with a green, solid line the A 1600 -β relation derived for z ∼ 0 starbursts by Meurer et al. (1999).This relation is based on a slightly redder intrinsic UV slope than the one we have assumed for z ∼ 2 star-forming galaxies and is consistent with the Calzetti et al. (2000) curve for z ∼ 0 starbursts.The median relation between A 1600 and β for SDSS galaxies, based on SED fitting and energy balance, more closely resembles that for an SMC dust curve than either Calzetti et al. (2000) relation, at least in the UV regime.dians diverge downwards from the z ∼ 0 median curve at the 6.5 − 9.4σ level.For completeness, we also show the running median A 1600 vs. stellar mass at z ∼ 2.3 if the SMC dust curve is assumed for the entire sample (Reddy et al. 2018), in which case the z ∼ 2.3 sample running median A 1600 falls significantly below the z ∼ 0 trend (here z ∼ 2.3 points are connected by a dotted curve).However, we argue below that this final set of assumptions is inconsistent with other previous MOSDEF results.Based on the two more likely, bracketing, prescriptions for converting β to A 1600 at z ∼ 2.3, we also infer no significant evolution from z ∼ 0 to z ∼ 2.3 in the relationship between A 1600 and M * .
It is worth emphasizing the distinct methodologies used for inferring A 1600 at z ∼ 0 and z ∼ 2.3, each of which based on the available information at that redshift.In Figure 3, we show the relationship between A 1600 and β for the z ∼ 0 SDSS sample, where β and A 1600 are inferred from SED fitting spanning from the UV to optical, along with longer-wavelength energy balance considerations as described in Salim et al. (2016).Along with the running median A 1600 for SDSS in bins of β, we overplot the relations from equations ( 1) and (2) as two bracketing cases, along with the z ∼ 0 Calzetti relation from Meurer et al. (1999).The median SDSS relation falls in between those dictated by the Calzetti et al. (2000) and SMC curves, but significantly closer to the relationship traced by the SMC curve (see also Salim et al. 2018).Accordingly, while the SDSS sample is described by a wide range of dust attenuation curves (Salim et al. 2016), on average these curves are steeper than the Calzetti et al. (2000) curve, and similar in UV slope to the SMC curve (Salim et al. 2018).
We do not have direct A 1600 estimates for the majority of MOSDEF galaxies due to a lack of individual far-IR photometric detections, and rely on β measurements and the assumption of different dust curves to infer A 1600 .At least within the high-mass (M * ∼ 10 10 − 10 11 M ), high-SFR (30 − 250M yr −1 ) range of the z ∼ 2.3 MOSDEF sample, where Spitzer/MIPS, Herschel/PACS and SPIRE detections have been achieved, SED fitting by Shivaei et al. (2016) suggests that the Calzetti et al. (2000) (and not an SMC) law provides an accurate description of the energy balance between rest-frame UV and far-IR.Similar dust-continuum measurements are now required for lower-mass, lower-luminosity z ∼ 2.3 MOSDEF galaxies, as in Reddy et al. (2018).

DISCUSSION
We find no significant evolution in the relationship between dust attenuation and stellar mass between z ∼ 0 and z ∼ 2.3.Our results join a list of several that arrive at similar conclusions, based on different proxies for dust attenuation (L IR /L UV , or IRX; the fraction of obscured star formation, f obscured ; A V ; and A 1600 ) and different multi-wavelength datasets (e.g., Heinis et al. 2014;Pannella et al. 2015;Bourne et al. 2017;Whitaker et al. 2017;McLure et al. 2018;Cullen et al. 2018).What is new here is the estimate of dust attenuation at z ∼ 2.3 based on a large set of individual Balmer decrement measurements, as well as a direct comparison of UV attenuation, A 1600 , at z ∼ 2.3 and z ∼ 0 based on GALEX data for the local SDSS comparison sample.As we now discuss, the lack of strong evolution in the attenuation vs. stellar mass relation has striking implications, based on what else is known about the evolution of the ISM in star-forming galaxies between z ∼ 0 and z ∼ 2.3.

The Connection between Attenuation and other ISM Properties
We start with the expression for dust attenuation at a given wavelength, λ.For a given dust optical depth, τ λ , and attenuation A λ , with the latter in magnitudes, we find: exp(−τ λ ) = 10 −0.4A λ (3) which corresponds to: We then recall the expression for τ λ , as a function of the dust mass absorption coefficient, κ λ , in units of m 2 kg −1 , the dust density, ρ dust , in units of kg m −3 , and the differential path length along the line of sight, ds: In a simplified model where κ λ is spatially independent and dust is smoothly distributed through-out the galaxy disk, the integral can be reexpressed as the product of the dust mass absorption coefficient and the dust mass surface density, (M dust /(πr 2 dust )), where M dust is the dust mass and r dust is the scale of the disk over which dust is distributed: Folding in equation ( 4), we can express A λ as a function of dust properties: Within the context of the simplified model presented here, the lack of significant evolution in A λ at fixed stellar mass from z ∼ 0 to z ∼ 2.3 implies that the product κ λ × (M dust /(πr 2 dust )) remains constant at fixed stellar mass.

Direct M dust Measurements
First, we consider what is known from direct measurements of dust masses at high redshift.The Atacama Large Millimeter/submillimeter Array (ALMA) is now beginning to enable M dust estimates for galaxies in the luminous end of the Luminous Infrared Galaxy (LIRG) regime (i.e., 10 11 L ≤ L IR ≤ 10 12 L ; Aravena et al. 2020;Shivaei et al. 2021, in prep), but does not yet cover the full range of stellar masses and SFRs in our sample.However, initial ALMA results from Magnelli et al. (2020) and Donevski et al. (2020) indicate signficant evolution in the (M dust /M * ) ratio at fixed M * from the local universe out to z ∼ 2. Based on stacked measurements of galaxies covered by the ALMA Spectroscopic Survey (ASPECS) large program, Magnelli et al. (2020) find a best-fit factor of 10 evolution in (M dust /M * ) for star-forming galaxies at a fiducial stellar mass of M * = 10 10.7 M between z = 0.45 and z = 2.0.Likewise, for starforming galaxies individually detected by ALMA, at a median redshift of z med = 2.39 and median stellar mass of M * = 10 11 , Donevski et al. (2020) find an order of magnitude increase in (M dust /M * ) relative to that observed in the local galaxy sample of Andreani et al. (2018).
This observed evolution in (M dust /M * ) exists in tension with model predictions from, e.g., Popping et al. (2017), which include a roughly constant relationship between M dust and M * over the redshift range z ∼ 0 to z ∼ 2. As highlighted by Popping et al. (2019), such models also significantly (by a factor of 2 − 3) underpredict the molecular gas content of z ∼ 2 galaxies, which likely contributes to the discrepancy between their predicted and observed dust masses.The tension between the predicted and observed evolution of ISM gas and dust masses needs to be addressed.4.1.2.The (M dust /M gas ) Ratio and Gas Surface

Density
We can also consider the question of evolution in A λ by re-expressing the dust mass surface density, (M dust /(πr 2 dust )), i.e., Σ dust , as the product (M dust /M gas ) × Σ gas .Here, Σ gas is the gas surface density.Accordingly, we can rewrite equation (7) as: With this equation, we can gain indirect constraints on the evolution of Σ dust from estimates of the evolution of galaxy metallicity and gas surface density.For this analysis, we consider the evolution in galaxy properties at a fiducial stellar mass of 10 10 M , close to the median stellar mass of both the z ∼ 2.3 MOSDEF and z ∼ 0 SDSS samples.
In the first step, we can use the evolution in metallicity from z ∼ 0 to z ∼ 2.3 to infer the evolution in (M dust /M gas ).As shown in Figure 4 and previous works (e.g., Steidel et al. 2014;Sanders et al. 2015Sanders et al. , 2018Sanders et al. , 2021)), z ∼ 2.3 galaxies have lower metallicity at fixed stellar mass relative to galaxies at z ∼ 0. In addition to the local and z ∼ 2.3 samples analyzed here, we plot the bestfit mass-metallicity relations from Sanders et al. (2021), which are consistent with our measure-Figure 4. 12 + log(O/H) vs. M * , i.e., the massmetallicity relation (MZR).Symbols and running medians for z ∼ 2.3 MOSDEF galaxies and local SDSS galaxies are as in Figure 1. 12 + log(O/H) is estimated for z ∼ 0 and z ∼ 2.3 galaxies following the prescriptions in Sanders et al. (2021).Also plotted are the bestfit z ∼ 0 and z ∼ 2.3 MZR relations from Table 3 of Sanders et al. (2021).The z ∼ 2.3 best-fit relation exactly follows the running median 12 + log(O/H) in bins of M * , while the z ∼ 0 best-fit relation falls slightly below our running median at log(M * /M ) < 10, due to small differences in SDSS comparison sample selection between the two works.At M * ∼ 10 10 M , however, there is agreement between the best-fit SDSS MZR and our corresponding running median value.There is clear evolution towards lower 12 + log(O/H) at fixed M * , which should be accompanied by a lower (M dust /M gas ) according to the (M dust /M gas ) vs. 12+log(O/H) relation in De Vis et al. (2019).
ments and show that our dust attenuation sample is representative of MOSDEF star-forming galaxies.The offset in metallicity is ∆12 + log(O/H) = 0.26 ± 0.02 dex towards lower metallicity at z ∼ 2.3.This difference in metallicity can be translated into a reduction in (M dust /M gas ) from z ∼ 0 to z ∼ 2.3 via the relationship between (M dust /M gas ) and 12 + log(O/H).Most recently, De Vis et al. (2019) constructed this relationship for a large sample of local galaxies with M dust , M gas , and 12 + log(O/H) measurements.While the slope, a, of the relationship varies in detail depending on which empirical strong-line metallicity indictator is adopted, Table 4 of De Vis et al. (2019) shows that the values for a are distributed around 2. Shapley et al. (2020) reported evidence for a lack of evolution in the log(M dust /M gas ) vs. 12 + log(O/H) relationship, at least at the high-mass end (10 10.5 M M * 10 11 M ), and, therefore, we assume that the local relation can be applied to the z ∼ 2.3 MOS-DEF sample.Accordingly, we adopt a slope of a = 2.0 ± 0.2 for the (M dust /M gas ) vs. 12 + log(O/H) relationship, and find that the decrease in 12 + log(O/H) implies a decrease of 0.52 ± 0.05 dex in (M dust /M gas ) at fixed mass (i.e., a linear factor of 3.3 ± 0.4).
Next, we consider evolution in the gas surface density of star-forming galaxies out to z ∼ 2.3.The molecular gas surface density evolves dramatically at fixed stellar mass over this redshift range.This evolution can be traced either by direct CO measurements (Tacconi et al. 2013), or else measurements of the change in SFR surface density, Σ SFR , coupled with an inversion of the Kennicutt-Schmidt (K-S) Law.For example, based on the dust-corrected Hα SFRs and rest-optical half-light radii for both our z ∼ 2.3 MOSDEF and z ∼ 0 SDSS comparison samples, we find an increase of a factor of ∼ 65 in the median Σ SFR (see also Shapley et al. 2019).Assuming a linear powerlaw slope for the K-S law (Tacconi et al. 2013), we infer an increase in molecular Σ gas by the same factor of 65 from z ∼ 0 to z ∼ 2.3.While the gas content of z ∼ 2 star-forming galaxies is wellapproximated by the molecular component (Tacconi et al. 2018), in local galaxies at the median stellar mass of our SDSS sample the molecular component only comprises ∼ 20% of the total molecular plus atomic gas mass (Catinella et al. 2018, Table 3).The evolution in total Σ gas from z ∼ 0 to z ∼ 2.3 is therefore reduced by a factor of ∼ 5 relative to the inferred evolution in molecular Σ gas , since the total Σ gas for z ∼ 0 galaxies is a factor of ∼ 5 higher than the molecular component alone.Total Σ gas is thus inferred to increase by a factor of ∼ 13.The net effect of the inferred evolution in (M dust /M gas ) and Σ gas is a factor of (M dust /M gas ) × Σ gas ∼ (1/3.3)× 13, i.e., a factor of ∼ 4.This factor corresponds to the increase in Σ dust , assuming that the spatial extent of dust and molecular gas is the same.ALMA has been used to obtain spatially-resolved dust-continuum and CO maps for small samples of massive (M * > 10 11 M ) and luminous (L > 10 12 L ) galaxies at z ∼ 2 (Tadaki et al. 2017;Calistro Rivera et al. 2018;Kaasinen et al. 2020), but a clear picture has yet to emerge from these measurements regarding the relative extents of dust and molecular gas emission.A larger sample of spatially resolved measurements of both dust continuum and CO emission is required for less luminous z ∼ 2.3 main sequence galaxies in the LIRG regime in order to understand if the extent of dust continuum emission evolves in the same manner as that of the molecular gas.

Explaining the Lack of Evolution
In the extremely simplified picture presented in Section 4.1.1,in order to maintain a fixed attenuation A λ at fixed stellar mass there must be evolution in either κ λ , r dust , or both, in the sense that κ λ is smaller and r dust is effectively larger at z ∼ 2.3 than at z ∼ 0. Alternatively (or in addition), following the discussion in Section 4.1.2,if there is a stronger dependence of (M dust /M gas ) on metallicity at z ∼ 2.3 than at z ∼ 0 (De Vis et al. 2019), such that a reduction of 0.26 dex in metallicity corresponds to a more extreme decrease in (M dust /M gas ), this effect would also help to explain the lack of evolution in A λ at fixed stellar mass.While additional data at lower metallicities is required to show the actual form of the relation at z ∼ 2.3, initial results from Shapley et al. (2020) suggest that the normalization in the (M dust /M gas ) vs. 12 + log(O/H) relation does not evolve at solar metallicity.
We now consider one of the first two factors: the spatial extent of dust, r dust .Evolution in r dust comprises another possibility for explaining the constant attenuation vs. stellar mass relation in the face of significant (M dust /M * ) evolution.If πr 2 dust is a factor of ∼ 10 larger at z ∼ 2.3 at fixed M * , corresponding to an increase of a factor of ∼ 3 in r dust , then (M dust /(πr 2 dust )) would remain constant.However, an evolution towards larger r dust at higher redshift and fixed stellar mass is in conflict with both recent numerical simulations of galaxy formation including dust radiative transfer (Popping et al. 2021), as well as preliminary resolved ALMA measurements of the evolution of dust sizes.Specifically, Fujimoto et al. (2017) show for a sample of luminous (L FIR ≥ 10 12 ) and massive (log(M * /M ) med ∼ 11) galaxies drawn from the DANCING-ALMA survey, that rest-frame far-IR sizes measured with ALMA (tracing dust-continuum) decrease over the range 1.5 ≤ z ≤ 4. Gómez-Guijarro et al. (2021) find a similar evolution towards smaller rest-frame far-IR continuum sizes as redshift increases over z ∼ 2−4, based on GOODS-ALMA-2.0,a blind survey conducted at 1.1 mm covering a similar luminosity range.For a sample of four z ∼ 1.5 − 2.0 galaxies with lower far-IR luminosities, in the LIRG range, Cheng et al. (2020) finds comparable r dust values to those of local LIRGS from the KINGFISH (Kennicutt et al. 2011) and GOALS (Armus et al. 2009) surveys.However, there is no evidence to date for larger r dust at higher redshift.
It may be that simply using an effective size, "r dust ," is insufficient to capture differences in the typical spatial distributions of dust at z ∼ 0 and z ∼ 2.3.For example, if dust distributions are patchier and clumpier at z ∼ 2 than locally, the observed dust attenuation for a given M dust will be lower (Witt & Gordon 2000;Seon & Draine 2016).Spatially-resolved maps of dust-continuum emission extending up to z ∼ 2 will be crucial for addressing this question.In addition, a detailed analysis of the shapes of dust attenuation curves for galaxies of similar mass at low and high redshift can be used to determine the dust geometry indirectly in such systems (e.g., Chevallard et al. 2013).
One final possibility for explaining the constant attenuation vs. stellar mass relation consists of evolution in κ λ , the wavelength-dependent dust mass absorption coefficient.This dust crosssection per unit dust mass enscapsulates many different dust properties, including dust-grain size distribution, grain morphology, density, and chemical composition.Recently, Clark et al. (2019) empirically determined maps of κ λ in two nearby face-on spiral galaxies.Clark et al. (2019) not only found significant variation of κ λ within the individual galaxies targeted (see also Bianchi et al. 2019), but also that κ λ is inversely correlated with gas surface density.Such an anti-correlation is not predicted by standard dust models, in which denser ISM regions are conducive to the growth of larger grains, which have higher emissivity (i.e., κ λ ) per unit mass.However, if lower κ λ is generally associated with higher Σ gas , the significantly higher Σ gas values at z ∼ 2.3 described above may result in a lower κ λ for these high-redshift galaxies than their low-redshift counterparts.
In summary, the roughly constant relationship between attenuation and stellar mass from z ∼ 0 to z ∼ 2.3 poses an important puzzle, given the significant evolution in the gas and dust content of the ISM at fixed stellar mass over the same redshift range.We have highlighted multiple possiblities for explaining the lack of evolution in attenuation vs. M * .In particular, these include a steeper relationship between (M dust /M gas ) and 12 + log(O/H) at z ∼ 2.3, such that the decrease in 12 + log(O/H) translates into lower (M dust /M gas ) at z ∼ 2.3 than in the local universe; more extended dust distributions at z ∼ 2.3 for galaxies at fixed stellar mass (though such a possibility seems inconsisent with both theory and preliminary observations), or, on the other hand, clumpier dust distributions at z ∼ 2.3; and a lower dust mass absorption coefficient κ λ .Directly measuring κ λ at z ∼ 2.3 seems beyond the reach of current facilities.However, determining (M dust /M gas ) vs. 12 + log(O/H) at subsolar metallicities, and obtaining spatially-resolved maps of the dust continuum emission for such galaxies is well within the scope of ALMA.Such observations should be highly prioritized in order to solve the puzzle of the non-evolving attenution vs. stellar mass relation.

Figure 1 .
Figure 1.Attenuation vs. M * , based on the Balmer-line ratio, Hα/Hβ.In each panel, z ∼ 2.3 MOSDEF galaxies are indicated with blue points, and a median z ∼ 2.3 error bar is shown in the corner of the plot.The greyscale histogram corresponds to the distribution of local SDSS galaxies.Running median Hα/Hβ line ratios and the corresponding magnitude of attenuation at the wavelength of Hα, A Hα , are calculated in bins of stellar mass.The z ∼ 2.3 running median is indicated in turquoise, while that for SDSS is plotted in red.Median error bars are shown for the z ∼ 2.3 sample, while those for the SDSS sample are smaller than the symbols.Left: Hα/Hβ ratio vs. M * .Right: A Hα vs. M * .

Figure 2 .
Figure 2. Attenuation vs. M * , based on the rest-UV continuum.Here attenuation is quantified as A 1600 , the magnitude of attenuation at a rest wavelength of 1600Å.As in Figure1, the local SDSS sample is indicated as a grey histogram, and its running median A 1600 estimated in bins of M * is shown with red, connected symbols.A 1600 is determined for SDSS galaxies as described inSalim et al. (2016).The running median A 1600 for the z ∼ 2.3 MOSDEF sample in bins of M * is shown in turquoise, using three different prescriptions to translate the measured UV slope value, β, to A 1600 .The solid curve shows A 1600 based on assuming aCalzetti et al. (2000) dust attenuation curve; the dashed curve shows A 1600 assuming aCalzetti et al. (2000) curve at 12 + log(O/H) ≥ 8.5 and an SMC curve at 12 + log(O/H) < 8.5; finally, the dotted curve shows A 1600 assuming an SMC curve for the entire MOSDEF z ∼ 2.3 sample (a scenario not favored by the MOSDEF data, but included for completeness).As in Figure1, median error bars are shown for the z ∼ 2.3 sample, while those for the SDSS sample are smaller than the symbols.