The MOSFIRE Deep Evolution Field Survey: Implications of the Lack of Evolution in the Dust Attenuation–Mass Relation to z ∼ 2^*

The analysis presented here is based on data from the MOSFIRE Deep Evolution Field (MOSDEF) 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 multiwavelength data sets 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-IR ground-based and Hubble Space Telescope photometry, as well as Spitzer/IRAC mid-IR measurements. MOSDEF 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.5, 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+\mathrm{log}({\rm{O}}/{\rm{H}})$) based on the subset of detected strong nebular emission lines drawn from [O ii]λ λ3726, 3729, [Ne iii]λ3869, Hβ, and [O iii]λ5007 (at the very least [O ii]λ λ3726, 3729, Hβ, and [O iii]λ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 multiwavelength photometric SEDs cataloged by the 3D-HST team (Skelton et al. 2014), where the near-IR photometry was corrected for the contribution of strong nebular emission lines (Sanders et al. 2021). For SED modeling, we used the FAST (Kriek et al. 2009) program, 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 $\mathrm{SFR}(t)\propto t\exp (-t/\tau )$. 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, active galactic nuclei (AGNs) were identified and removed based on their X-ray and IR properties, as well as those with $\mathrm{log}([\mathrm{Ne}\,{\rm\small{II}}]/{\rm{H}}\alpha )\gt -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 $\mathrm{log}{({M}_{* }/{M}_{\odot })}_{\mathrm{med}}=9.88$, and a median dust-corrected Hα-based SFR of SFR_med = 29 M_⊙ yr⁻¹. 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₁₆₀₀) 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 300 M_⊙, Z_* = 0.02, and $\mathrm{log}(\mathrm{age}/\mathrm{yr})=8.0$, we find the relation

$$\begin{eqnarray}&&{A}_{1600}=2.13\times \beta +5.04\end{eqnarray} \tag{ 1 }$$

for a Calzetti et al. (2000) dust-attenuation law, and

$$\begin{eqnarray}&&{A}_{1600}=1.07\times \beta +2.52\end{eqnarray} \tag{ 2 }$$

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₁₆₀₀ and β (A₁₆₀₀ = 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 the constant-star-formation, solar-metallicity models of Bruzual & Charlot (2003) 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+\mathrm{log}({\rm{O}}/{\rm{H}})\geqslant 8.5$, while one resembling the SMC curve seems to apply at $12+\mathrm{log}({\rm{O}}/{\rm{H}})\lt 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) has presented evidence that the Calzetti et al. (2000) curve applies over a wide range of stellar masses ($\mathrm{log}({M}_{* }/{M}_{\odot })\geqslant 9.75$), and, correspondingly, metallicities. Accordingly, our two favored methods for translating β into A₁₆₀₀ 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+\mathrm{log}({\rm{O}}/{\rm{H}})\geqslant 8.5$ and the SMC curve at $12+\mathrm{log}({\rm{O}}/{\rm{H}})\lt 8.5$.

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. ¹² 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 [O ii]λ λ3726, 3729, Hβ, Hα, and [Ne ii]λ6584, and a 3σ detection for [O iii]λ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 $\mathrm{log}{({M}_{* }/{M}_{\odot })}_{\mathrm{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.3 M_⊙ yr⁻¹, i.e., a factor of ∼20 lower, which reflects the evolution of the star-forming 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 of Sanders et al. (2021), in order to provide the fairest comparison with respect to the z ∼ 2.3 MOSDEF sample. Specifically, SDSS emission-line fluxes were corrected for the contribution of diffuse ionized gas (DIG) following Sanders et al. (2017), and then the combination of [O iii]λ5007/Hβ and [O iii]λ5007/[O ii]λ λ3726, 3729 was fit simultaneously with the Sanders et al. (2021) z ∼ 0 DIG-corrected metallicity calibrations to yield $12+\mathrm{log}({\rm{O}}/{\rm{H}})$.

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 SED-based 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₁₆₀₀. GALEX near-UV and far-UV coverage is available for 30,204 galaxies in our SDSS comparison sample. We have confirmed that this GALEX subsample is representative of the larger SDSS-only sample, in terms of the relationship between Balmer line ratios and stellar mass (see Section 3).

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

$$\begin{eqnarray}&&\exp (-{\tau }_{\lambda })={10}^{-0.4{A}_{\lambda }},\end{eqnarray} \tag{ 3 }$$

which corresponds to

$$\begin{eqnarray}&&{A}_{\lambda }=\displaystyle \frac{{\tau }_{\lambda }}{0.4\mathrm{ln}(10)}=1.086\times {\tau }_{\lambda }.\end{eqnarray} \tag{ 4 }$$

We then recall the expression for τ_λ , as a function of the dust mass absorption coefficient, κ_λ , in units of m² kg⁻¹, the dust density, ρ_dust, in units of kg m⁻³, and the differential path length along the line of sight, ds:

$$\begin{eqnarray}&&{\tau }_{\lambda }=\int {\kappa }_{\lambda }\times {\rho }_{\mathrm{dust}}{ds}.\end{eqnarray} \tag{ 5 }$$

In a simplified model where κ_λ is spatially independent and dust is smoothly distributed throughout the galaxy disk, the integral can be reexpressed as the product of the dust mass absorption coefficient and the dust mass surface density, $({M}_{\mathrm{dust}}/(\pi {r}_{\mathrm{dust}}^{2}))$, where M_dust is the dust mass and r_dust is the scale of the disk over which dust is distributed:

$$\begin{eqnarray}&&{\tau }_{\lambda }\simeq {\kappa }_{\lambda }\times \left(\displaystyle \frac{{M}_{\mathrm{dust}}}{\pi {r}_{\mathrm{dust}}^{2}}\right).\end{eqnarray} \tag{ 6 }$$

Folding in Equation (4), we can express A_λ as a function of dust properties:

$$\begin{eqnarray}&&{A}_{\lambda }\simeq 1.086\times {\kappa }_{\lambda }\times \left(\displaystyle \frac{{M}_{\mathrm{dust}}}{\pi {r}_{\mathrm{dust}}^{2}}\right).\end{eqnarray} \tag{ 7 }$$

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 ${\kappa }_{\lambda }\times ({M}_{\mathrm{dust}}/(\pi {r}_{\mathrm{dust}}^{2}))$ remains constant at fixed stellar mass.

4.1.1. Direct M_dust Measurements

4.1.2. The (M_dust/M_gas) Ratio and Gas Surface Density

We can also consider the question of evolution in A_λ by reexpressing the dust mass surface density, $({M}_{\mathrm{dust}}/(\pi {r}_{\mathrm{dust}}^{2}))$, 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

$$\begin{eqnarray}&&{A}_{\lambda }\simeq 1.086\times {\kappa }_{\lambda }\times \left(\displaystyle \frac{{M}_{\mathrm{dust}}}{{M}_{\mathrm{gas}}}\right)\times {{\rm{\Sigma }}}_{\mathrm{gas}}.\end{eqnarray} \tag{ 8 }$$

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¹⁰ 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. 2015, 2018, 2021), z ∼ 2.3 galaxies have a lower metallicity at a fixed stellar mass relative to galaxies at z ∼ 0. In addition to the local and z ∼ 2.3 samples analyzed here, we plot the best-fit mass–metallicity relations from Sanders et al. (2021), which are consistent with our measurements and show that our dust-attenuation sample is representative of MOSDEF star-forming galaxies. The offset in metallicity is ${\rm{\Delta }}12+\mathrm{log}({\rm{O}}/{\rm{H}})=0.26\pm 0.02$ dex toward 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+\mathrm{log}({\rm{O}}/{\rm{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+\mathrm{log}({\rm{O}}/{\rm{H}})$ measurements. While the slope, a, of the relationship

$$\begin{eqnarray}&&\mathrm{log}\left(\displaystyle \frac{{M}_{\mathrm{dust}}}{{M}_{\mathrm{gas}}}\right)=a\times [12+\mathrm{log}({\rm{O}}/{\rm{H}})]+b\end{eqnarray} \tag{ 9 }$$

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 $\mathrm{log}({M}_{\mathrm{dust}}/{M}_{\mathrm{gas}})$ versus $12+\mathrm{log}({\rm{O}}/{\rm{H}})$ relationship, at least at the high-mass end (10^10.5 M_⊙ ≲ M_* ≲ 10¹¹ M_⊙), and, therefore, we assume that the local relation can be applied to the z ∼ 2.3 MOSDEF sample. Accordingly, we adopt a slope of a = 2.0 ± 0.2 for the (M_dust/M_gas) versus $12+\mathrm{log}({\rm{O}}/{\rm{H}})$ relationship, and find that the decrease in $12+\mathrm{log}({\rm{O}}/{\rm{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).

Zoom In Zoom Out Reset image size

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 power-law 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 well approximated 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¹¹ M_⊙) and luminous (L > 10¹² 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.

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) versus $12+\mathrm{log}({\rm{O}}/{\rm{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 versus stellar mass relation in the face of significant (M_dust/M_*) evolution. If $\pi {r}_{\mathrm{dust}}^{2}$ 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}_{\mathrm{dust}}/(\pi {r}_{\mathrm{dust}}^{2}))$ would remain constant. However, an evolution toward 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. 2022), 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¹²) and massive ($\mathrm{log}{({M}_{* }/{M}_{\odot })}_{\mathrm{med}}\sim 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. (2022) find a similar evolution toward 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 versus stellar mass relation consists of evolution in κ_λ , the wavelength-dependent dust mass absorption coefficient. This dust cross-section 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 anticorrelation 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 versus M_*. In particular, these include a steeper relationship between (M_dust/M_gas) and $12+\mathrm{log}({\rm{O}}/{\rm{H}})$ at z ∼ 2.3, such that the decrease in $12+\mathrm{log}({\rm{O}}/{\rm{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) versus $12+\mathrm{log}({\rm{O}}/{\rm{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 nonevolving attenution versus stellar mass relation.