A Model Connecting Galaxy Masses, Star Formation Rates, and Dust Temperatures across Cosmic Time

For each individual galaxy, we assume a spherically symmetric system in which the dust mass density, ρ_dust, has a power-law distribution, γ. The frequency-dependent optical depth due to dust is defined

$$\begin{eqnarray}{\tau }_{\nu } & = & {\displaystyle \int }_{0}^{{R}_{\mathrm{dust}}}{\rho }_{\mathrm{dust}}(r){\kappa }_{\nu }{dr}\\ & = & {\displaystyle \int }_{0}^{{R}_{\mathrm{dust}}}{\rho }_{0}{r}^{-\gamma }{\kappa }_{\nu }{dr},\end{eqnarray} \tag{ 4 }$$

where ρ₀ is the central mass density of a given galaxy, and ${R}_{\mathrm{dust}}$ is the radial extent of dust. The total dust mass in a galaxy also depends on ρ₀ as

$$\begin{eqnarray}{M}_{\mathrm{dust}} & = & {\displaystyle \int }_{0}^{{R}_{\mathrm{dust}}}{\rho }_{0}{r}^{-\gamma }4\pi {r}^{2}{dr}\\ & = & 4\pi {\rho }_{0}\displaystyle \frac{{R}_{\mathrm{dust}}^{3-\gamma }}{3-\gamma }.\end{eqnarray} \tag{ 5 }$$

By equating Equations (4) and (5), we may express τ_dust,ν in terms of dust mass:

$$\begin{eqnarray}&&{\tau }_{\nu }=\displaystyle \frac{{M}_{\mathrm{dust}}({M}_{\star },z)}{4\pi {R}_{\mathrm{dust}}^{2}}\displaystyle \frac{3-\gamma }{1-\gamma }{\kappa }_{\nu },\end{eqnarray} \tag{ 6 }$$

where we have drawn attention to the dependency of ${M}_{\mathrm{dust}}$ on the stellar mass of the galaxy, ${M}_{\star }$, and on the redshift, z. We let γ = 0; we justify this choice in Section 4. To evaluate Equation (6), then, we need to model r, κ_ν, and ${M}_{\mathrm{dust}}$, all of which depend on z.

2.1.1. Galactic Radius

To acquire a rough approximation of galaxy disk sizes, we adopt the basic picture of Fall & Efstathiou (1980) and others (e.g., Mo et al. 1998; Sommerville & Primack 1999), in which the collapsing gas of a forming galaxy acquires the same specific angular momentum as the dark matter halo, and this angular momentum is conserved as the gas cools. The specific angular momentum is often expressed using the dimensionless spin parameter, $\lambda =J| E{| }^{1/2}{G}^{-1}{M}^{-5/2}$, where J is the angular momentum, E is the total energy of the halo, G is Newton's gravitational constant, and M is the mass (e.g., Peebles 1969; Mo et al. 1998; Somerville & Primack 1999). For a halo with a singular isothermal density profile ρ ∝ r⁻², the disk exponential scale radius is given by ${r}_{d}=\lambda {R}_{\mathrm{halo}}/\sqrt{2}$, where R_halo is the virial radius of the dark matter halo.

Somerville et al. (2018) explored λ using empirical constraints from z ∼ 0.1–3 galaxies in the GAMA survey (Driver et al. 2011; Liske et al. 2015) and CANDELS survey (Grogin et al. 2011; Koekemoer et al. 2011). Using relationships from their halo abundance matching model, they mapped galaxy stellar mass to halo mass and inferred from these results a median value of the spin parameter, λ = 0.036, corresponding to a ratio between galaxy half-mass radius and halo size of ${R}_{\mathrm{gal}}$/R_halo = 0.018. They found that λ is roughly independent of stellar mass and exhibits weak dependence on redshift. We adopt this value for λ and make the simplifying assumption that the radial extent of the dust, ${R}_{\mathrm{dust}}$, is equal to ${R}_{\mathrm{gal}}$,

$$\begin{eqnarray}&&{R}_{\mathrm{dust}}={R}_{\mathrm{gal}}=0.018\,{R}_{\mathrm{halo}}.\end{eqnarray} \tag{ 7 }$$

In Section 4, we discuss some of the limitations of assuming a single value for the spin parameter.

To determine R_halo, we adopt the R_halo–M_halo relation of Loeb & Furlanetto (2013):

$$\begin{eqnarray}{R}_{\mathrm{halo}} & = & 0.784{\left[\displaystyle \frac{{{\rm{\Omega }}}_{m}}{{{\rm{\Omega }}}_{m}(z)}\displaystyle \frac{{{\rm{\Delta }}}_{c}}{18{\pi }^{2}}\right]}^{2}\left({\displaystyle \frac{{M}_{\mathrm{halo}}}{{10}^{8}{M}_{\odot }}}^{1/3}\right)\\ & & \times \left(\displaystyle \frac{10}{1+z}\right)\,{h}^{-2/3}\,\mathrm{kpc},\end{eqnarray} \tag{ 8 }$$

where

$$\begin{eqnarray}{{\rm{\Delta }}}_{c} & = & 18{\pi }^{2}+82d-39{d}^{2}\\ d & = & {{\rm{\Omega }}}_{m}(z)-1\\ {{\rm{\Omega }}}_{m}(z) & = & \displaystyle \frac{{{\rm{\Omega }}}_{m}{(1+z)}^{3}}{{{\rm{\Omega }}}_{m}{(1+z)}^{3}+{{\rm{\Omega }}}_{{\rm{\Lambda }}}}.\end{eqnarray} \tag{ 9 }$$

For ${M}_{\mathrm{halo}}$ in Equation (8), we employ the stellar-mass–halo-mass (SMHM) relations of Behroozi et al. (2013a), discussed in detail in Section 2.2. In Figure 2, we show how the ${R}_{\mathrm{gal}}$–${M}_{\star }$ relation evolves with redshift.

Zoom In Zoom Out Reset image size

2.1.2. Opacity

For most extragalactic environments, especially for distant galaxies, we have poor knowledge of the composition and size distribution of dust grains, encapsulated in κ_ν in Equation (6). For values of ν > 10¹² Hz, we adopt the Galactic extinction laws of Mathis (1990) and Li & Draine (2001). We calculate an interpolated function for κ_ν based on a combination of these two models, since the Mathis (1990) model extends to lower frequencies, down to ∼10¹² Hz, while the Draine & Li (2001) model extends up to ∼10¹⁸ Hz. For frequencies below ν ≤ 10¹² Hz, we adopt the Beckwith et al. (1990) power-law treatment for opacity,

$$\begin{eqnarray}&&{\kappa }_{\nu }=0.1{\left(\displaystyle \frac{\nu }{1000\mathrm{GHz}}\right)}^{\beta }\,{\mathrm{cm}}^{2}\,{{\rm{g}}}^{-1},\end{eqnarray} \tag{ 10 }$$

where β = 2. The Beckwith et al. model, which originally assumed β = 1, was calibrated to match the emmisivity properties of dust around Galactic protoplanetary disks. The frequency-dependence of κ_ν is uncertain and depends on the size and composition of dust grains. Using a different slope or normalization for κ_ν (e.g., James et al. 2002; da Cunha et al. 2008; Dunne et al. 2011; Clark et al. 2016) would affect the calculation of the optical depth (Equation (6)) in our model and thus the dust temperature. We describe how changes in the normalization or slope for κ_ν affect our results in Section 4.3.

2.1.3. Dust Mass

Next, we determine ${M}_{\mathrm{dust}}$, which relates to the gas mass of a galaxy, ${M}_{\mathrm{gas}}$, via a dust-to-gas ratio, DGR:

$$\begin{eqnarray}&&{M}_{\mathrm{dust}}\equiv {M}_{\mathrm{gas}}\times \mathrm{DGR}.\end{eqnarray} \tag{ 11 }$$

A number of studies have explored the correlation between galactic gas mass and stellar mass or SFR (e.g., Sargent et al. 2014; Zahid et al. 2014). Yet, empirical relations for the total gas mass as a function of stellar mass or SFR rarely extend to redshifts higher than z ≈ 2 (e.g., Zahid et al. 2014), beyond which measurements of the H i gas mass are unreliable, and the gas fraction of galaxies is expected to be increasingly dominated by molecular gas.

After investigating different prescriptions for the gas mass, we decided to follow Zahid et al. (2014), who fit a stellar-mass–metallicity (MZ) relation for star-forming galaxies at z ≲ 1.6. They assume that ${M}_{\star }$/${M}_{\mathrm{gas}}$ ≈ (${M}_{\star }$/M₀)^γ, where M₀ is a metallicity-dependent, characteristic mass, above which ${ \mathcal Z }$ approaches a saturation limit, and γ is a power-law index. By combining this expression with their fit for the MZ relation, Zahid et al. (2014) determine

$$\begin{eqnarray}&&{M}_{\mathrm{gas}}({M}_{\star },z)=3.87\times {10}^{9}{(1+z)}^{1.35}{\left(\displaystyle \frac{{M}_{\star }}{{10}^{10}{M}_{\odot }}\right)}^{0.49}.\end{eqnarray} \tag{ 12 }$$

To determine the stellar mass in the above equation, we use the SMHM models of Behroozi et al. (2013a), who constrain average galaxy stellar masses and SFRs as a function of halo mass (see Section 2.2). We extrapolate the M_gas–${M}_{\star }$ relation of Zahid et al. (2014) to redshifts z > 1.6, as shown in Figure 2. We discuss potential consequences of this choice and investigate an alternative prescription for ${M}_{\mathrm{gas}}$ in Section 4.

Equation (11) also depends on the dust-to-gas ratio, DGR. Since dust is composed of heavy elements and traces the metal abundance of galaxies, the problem of determining the DGR can be reduced to one of determining the metallicity, ${ \mathcal Z }$. Several authors have conducted observational investigations of the MZ relation in nearby galaxies (e.g., Lequeux et al. 1979; Lee et al. 2006; Berg et al. 2012; Zahid et al. 2012, 2014) and in distant galaxies out to z ≲ 3 (Savaglio et al. 2005; Erb et al. 2006; Maiolino et al. 2008; Yabe et al. 2012; Zahid et al. 2013; Hunt et al. 2016). We adopt the prescription of Hunt et al. (2016), who compiled observations of ∼1000 galaxies up to z ∼ 3.7, with metallicities spanning 2 orders of magnitude, SFRs spanning 6 orders of magnitude, and stellar masses spanning 5 orders of magnitude. Using a principal component analysis, Hunt et al. (2016) found

$$\begin{eqnarray}&&{ \mathcal Z }=-0.14\mathrm{log}(\mathrm{SFR})+0.37\mathrm{log}({M}_{\star })+4.82,\end{eqnarray} \tag{ 13 }$$

where, by convention, ${ \mathcal Z }\equiv 12+\mathrm{log}({\rm{O}}/{\rm{H}})$ is defined in terms of the gas-phase oxygen abundance (O/H). We extrapolate Equation (13) for galaxies at z > 3.7, as shown in Figure 2.

We now relate ${ \mathcal Z }$ to the DGR, using an empirical formula determined for local galaxies. Rémy-Ruyer et al. (2014) evaluated the gas-to-dust ratio as a function of metallicity for nearby galaxies spanning the range between $1/50\,{ \mathcal Z }$ and $2\,{ \mathcal Z }$. The authors provide alternative functional forms for the gas-to-dust ratio versus metallicity relationship, depending on the CO-to-H₂ conversion factor they employed to estimate the total amount of molecular mass in a galaxy. We adopt the function they derive assuming a metallicity-dependent conversion factor (as opposed to the standard Milky Way CO-to-H₂ conversion factor). In terms of the DGR,

$$\begin{eqnarray}\mathrm{log}\left(\displaystyle \frac{\mathrm{DGR}}{{\mathrm{DGR}}_{\odot }}\right)=\left\{\begin{array}{lc}\mathrm{log}\left(\tfrac{{ \mathcal Z }}{{Z}_{\odot }}\right) & \mathrm{if}\ { \mathcal Z }\gt 0.26\,{Z}_{\odot }\\ 3.15\mathrm{log}\left(\tfrac{{ \mathcal Z }}{{Z}_{\odot }}\right)+1.25 & \mathrm{if}\ { \mathcal Z }\leqslant 0.26\,{Z}_{\odot },\end{array}\right.\end{eqnarray} \tag{ 14 }$$

where $\mathrm{log}({\mathrm{DGR}}_{\odot })=-2.21$ (Zubko et al. 2004).

In the previous subsections, we have modeled r, κ_ν, ${M}_{\mathrm{gas}}$, and the DGR. Equation (8) depends on ${M}_{\mathrm{halo}}$, Equation (12) depends on ${M}_{\star }$, and Equation (13) depends on ${M}_{\star }$ and the SFR. In the next section, we describe the self-consistent models of Behroozi et al. (2013a) that we use to determine the relation between ${M}_{\mathrm{halo}}$, ${M}_{\star }$, and the SFR at arbitrary redshifts.

Behroozi et al. (2013a) use empirical forward-modeling to constrain the evolution of the stellar-mass–halo-mass relationship (SMHM; $\mathrm{SM}({M}_{h},z)$). At fixed redshift, the adopted model for $\mathrm{SM}({M}_{h},z)$ has six parameters, which control the characteristic stellar mass, halo mass, faint-end slope, massive-end cutoff, transition region shape, and scatter of the SMHM relationship. For each parameter, there are three variables that control its redshift scaling at low (z = 0), mid (z = 1–2), and high (z > 3) redshift, with constant (i.e., no) scaling beyond z = 8.5 to prevent unphysical early galaxy formation. Additional nuisance parameters include systematic uncertainties in observed galaxy stellar masses and SFRs. Any choice of model in this parameter space gives a mapping from simulated dark matter halo catalogs (Behroozi et al. 2013b) to mock galaxy catalogs. Comparing these mock catalogs with observed galaxy number counts and SFRs from z = 0 to z = 8 results in a likelihood for a given model choice, so these constraints, combined with an MCMC algorithm, result in a posterior distribution for the allowed SMHM relationships. Average SFRs and histories for galaxies are inferred from averaged halo assembly histories (including mergers) combined with the best-fitting model for SM(M_h, z).

To determine the specific luminosity of all the stars in a galaxy, L_ν, we use version 3.0 of the Flexible Stellar Population Synthesis code (FSPS; Conroy et al. 2009; Conroy & Gunn 2010) to model the SEDs from 91 Å to 1000 μm. For each halo and each redshift, we supply the star formation history (SFH) to FSPS in tabular form, including the the time-dependent metallicity of newly born stars. Within FSPS, the SFH is linearly interpolated between the supplied time points, and the appropriate single-stellar population (SSP) weights are calculated. The spectra of the SSPs are then summed with these weights applied to produce a model galaxy spectrum corresponding to the last time-point of the supplied SFH. This model spectrum is also projected onto filter transmission curves to produce broadband rest-frame photometry in several standard filter sets. The total surviving stellar mass (including remnants) and bolometric luminosity are also calculated. No dust attenuation or IGM attenuation is applied, and we do not include nebular emission from H ii regions or dust emission.

The base SSP spectra are generated assuming a fully sampled Salpeter IMF from 0.08 to 120 ${M}_{\odot }$. We use the "Padova2007" isochrones (Bertelli et al. 1994; Girardi et al. 2000; Marigo et al. 2008) for stars less than 70 ${M}_{\odot }$. These are combined with Geneva isochrones for higher-mass stars based on the high mass-loss rate evolutionary tracks (Schaller et al. 1992; Meynet & Maeder 2000) and the post-AGB evolutionary tracks of Vassiliadis & Wood (1994). For the stellar spectra we use the BaSeL3.1 theoretical stellar library of Westera et al. (2002), augmented with the higher-resolution empirical MILES library (Sánchez-Blázquez et al. 2006) in the optical. The spectra of the OB stars are from Smith et al. (2002), and the spectra of the post-AGB stars are from Rauch (2003). The treatment of TP-AGB spectra and isochrones is described in Villaume et al. (2015). All FSPS variables that affect the isochrones and stellar spectra are set at their default values.

Figure 3 presents galactic dust mass as a function of stellar mass from z = 0 to z = 9.5. We find that for a fixed galactic stellar mass, ${M}_{\mathrm{dust}}$ increases with increasing redshift, with higher-mass galaxies displaying slightly larger increases in ${M}_{\mathrm{dust}}$ than their lower-mass counterparts. In Table 2, we summarize the predicted values of ${M}_{\mathrm{dust}}$, as well as other parameters derived in this paper, for a 10¹² ${M}_{\odot }$ galaxy—to show that our predictions compare favorably with quantities observed in the Milky Way.

Zoom In Zoom Out Reset image size

Table 2.  Parameters Derived in This Paper for Galaxies Having ${M}_{\mathrm{halo}}$ = 10¹² ${M}_{\odot }$

	z = 0	z = 1	z = 2	z = 3	z = 4	z = 6	z = 9.5
${M}_{\star }$ (1010 ${M}_{\odot }$)	2.7	2.8	2.2	2.1	2.3	1.7	0.79
SFR (${M}_{\odot }$ yr−1)	0.76	15	32	52	82	88	59
${R}_{\mathrm{gal}}$ (kpc)	4.7	2.8	2.0	1.5	1.2	0.85	0.57
${ \mathcal Z }$	8.7	8.5	8.4	8.4	8.4	8.3	8.2
DGR	1/160	1/241	1/292	1/315	1/326	1/367	1/463
AV (mag)	0.79	1.6	2.6	4.5	7.2	11	14
${M}_{\mathrm{gas}}$ (1010 ${M}_{\odot }$)	0.6	1.6	2.5	3.7	5.1	7.0	8.3
${M}_{\mathrm{dust}}$ (107 ${M}_{\odot }$)	3.9	6.7	8.5	12	16	19	18
${T}_{\mathrm{dust}}$ (K)	34	51	57	58	60	58	46

Note. From top to bottom: stellar mass, SFR, half-mass radius, metallicity, dust-to-gas ratio, visual extinction due to dust, gas mass, dust mass, and dust temperature.

Download table as:  ASCIITypeset image

As the second panel in Figure 3 shows, the dust-to-stellar mass ratio ${M}_{\mathrm{dust}}$/${M}_{\star }$ first rises, and then at a characteristic value of ${M}_{\star }$, the ratio decreases with increasing stellar mass. This turnover results from the dependency of ${M}_{\mathrm{dust}}$ on metallicity (via the DGR), which changes its functional form at ${ \mathcal Z }=0.26$ (Equation (14)). The ratio ${M}_{\mathrm{dust}}$/${M}_{\star }$ arises from a competition between the gas metallicity decreasing with redshift and the SFR increasing. For a fixed stellar mass, ${M}_{\mathrm{dust}}$/${M}_{\star }$ steadily increases with redshift. In Figure 4, we compare our results with the observations of Béthermin et al. (2015), who measure the gas and dust content of massive (∼6 × 10¹⁰ ${M}_{\odot }$) main-sequence galaxies and find little systematic variation of the dust-to-stellar mass ratio as a function of redshift up to z = 4. Although our results slightly underpredict the observations and suggest a small increase in ${M}_{\mathrm{dust}}$/${M}_{\star }$ over the same redshift range, they are compatible with the Béthermin et al. observations at the 1σ level. There are various ways to explain the lower dust-to-stellar mass ratios: low metallicities and dust-to-gas ratios, or differences in the SFRs and gas masses. It happens that the gas masses we derive for ∼6 × 10¹⁰ ${M}_{\odot }$ galaxies are in good agreement with the gas masses measured by Béthermin et al. (2015). However, if we used a different prescription for the metallicity, this would alter our results for ${M}_{\mathrm{dust}}$. For instance, the fundamental metallicity relation of Mannucci et al. (2010) yields higher metallicities than the Hunt et al. (2016) prescription we use here. Moreover, the SFRs derived by Béthermin et al. for their galaxy sample are slightly higher than the SFRs we use from the Behroozi et al. (2013) model, since the latter considers average SFRs, including contributions from both star-forming and quiescent galaxies. For a fixed stellar mass, higher SFR results correspond to lower metallicities. In Figure 4 we plot the redshift evolution ${M}_{\mathrm{dust}}$/${M}_{\star }$, where ${M}_{\mathrm{dust}}$ is derived using the metallicity prescription of Mannucci et al. (2010) and where the SFRs are corrected to include star-forming galaxies only. This correction is accomplished by dividing the SFRs by (1 − f_q), where f_q, the fraction of quiescent galaxies at a given redshift, is given by ${f}_{q}={[{({M}_{\star }/{10}^{10.2+0.5z}{M}_{\odot })}^{-1.3}+1]}^{-1}$ (Behroozi et al. 2013). The resulting curve for the evolution of ${M}_{\mathrm{dust}}$/${M}_{\star }$ now slightly overpredicts the Béthermin et al. observations at z < 3 and does not vary monotonically with redshift.

Zoom In Zoom Out Reset image size

In Figure 5, we again show ${M}_{\mathrm{dust}}$ as a function of ${M}_{\star }$, this time overplotting observations from the literature. In the first panel, we overplot observations from Rémy-Ruyer et al. (2015) and the Herschel Reference Survey (HRS; Ciesla et al. 2014; Boselli et al. 2015), and we find good agreement between our model and observed dust masses at z = 0, for stellar masses ranging from about 10⁷ to 10¹¹ ${M}_{\odot }$. At z = 1 and z = 2, the dust masses predicted by our model are in good agreement with the observations by Santini et al. (2014). We underpredict the observations by da Cunha et al. (2015) at these z = 1 and 2 by roughly 0.5 to 1 dex. This is not too surprising, however, since the da Cunha et al. (2015) observations are of submillimeter galaxies (SMGs), which have higher than typical dust infrared luminosities (>10¹² L_⊙), which are driven by high SFRs in excess of 100 ${M}_{\odot }$ yr⁻¹. Sub-millimeter selection basically selects for SFR, and the da Cunha et al. (2015) SMGs have SFRs about a factor of 3 higher than main-sequence galaxies of the same stellar mass. Galaxies such as those in the da Cunha et al. sample, which includes some starbursts, may have quickly enriched the ISM with metals and dust on timescales shorter than that for main-sequence galaxies and are not accounted for by the models adopted in this study. On the low-mass end of galaxies (${M}_{\star }$ < 10¹⁰), it will be interesting to see how well our model reproduces the ${M}_{\mathrm{dust}}$–${M}_{\star }$ trend at redshifts z ≥ 1. However, as we show in Section 3.4, detections of the dust emission in large samples of low-mass, high-redshift galaxies would be a challenging prospect for observing programs with present-day telescopes.

Zoom In Zoom Out Reset image size

We find that the relation between galaxy dust mass and stellar mass can be parameterized as a broken power law,

$$\begin{eqnarray}\mathrm{log}\left(\displaystyle \frac{{M}_{\mathrm{dust}}}{{M}_{\mathrm{dust},0}}\right)=\left\{\begin{array}{lc}{\alpha }_{1}\mathrm{log}\left(\displaystyle \frac{{M}_{\star }}{{M}_{\star ,0}}\right) & \mathrm{if}\ {M}_{\star }\leqslant {M}_{\star ,0}\\ {\alpha }_{2}\mathrm{log}\left(\displaystyle \frac{{M}_{\star }}{{M}_{\star ,0}}\right) & \mathrm{if}\ {M}_{\star }\gt {M}_{\star ,0},\end{array}\right.\end{eqnarray} \tag{ 15 }$$

where M_dust,0 is the zero-point of the dust mass, and α_1,2 are the slopes below and above M_⋆,0, the stellar mass at a given redshift where the break in the power law occurs. The break point in the ${M}_{\mathrm{dust}}$–${M}_{\star }$ relation results from the prescription we use for the dust-to-gas ratio, Equation ((14); Rémy-Ruyer et al. 2014), also a broken power law, which depends on ${M}_{\star }$ and the SFR via the metallicity. Both M_dust,0 and M_⋆,0 are functions of redshift.

The amount of dust in a galaxy is determined by the amount of available metals and gas, both of which are linked to star formation activity. Recent observational studies have shown that the dust-to-gas ratio of nearby galaxies may be characterized as a function of the gas-phase metallicity (Rémy-Ruyer et al. 2014). Theoretical work by Popping et al. (2017), who use SAMs to follow the production of interstellar dust, reproduces this observed trend and suggests that it is driven by the accretion of metals onto dust grains and the density of cold gas. Like Popping et al. (2017), we find that the normalization of the dust–mass–stellar–mass relation, M_dust,0, increases from z = 0 to z = 9.5.

We perform least-squares fits to the predicted curves in Figures 3 and 5 to determine the parameters in Equation (15). The results for z = 0 to z = 9.5 are summarized as follows:

$$\begin{eqnarray}{\alpha }_{1} & = & 1.20\pm 0.02\\ {\alpha }_{2} & = & 0.75\pm 0.02\\ \mathrm{log}{M}_{\mathrm{dust},0} & = & (6.0\pm 0.1)+(1.8\pm 0.1)\mathrm{log}(1+z)\\ \mathrm{log}{M}_{\star ,0} & = & (8.4\pm 0.1)+(1.0\pm 0.2)\mathrm{log}(1+z).\end{eqnarray} \tag{ 16 }$$

In Figure 6 we present contour maps of the optical depth due to dust in terms of visual extinction, A_V = 1.086τ_V (e.g., Draine 2011), calculated in the rest frames of the galaxies. The left panel of Figure 6 displays the redshift evolution of A_V in terms of ${M}_{\mathrm{halo}}$, and the right panel displays the evolution of A_V in terms of ${M}_{\star }$. The maps demonstrate that for constant values of ${M}_{\mathrm{halo}}$ or ${M}_{\star }$, A_V increases with increasing z. For instance, while a 10¹² ${M}_{\odot }$ halo mass galaxy has an extinction due to dust of A_V ≈ 0.8 mag at z = 0, a similar galaxy at z = 4 has an extinction of A_V ≈ 7 mag. This basic trend holds for the optical depth at other wavelengths, with the galaxies having overall higher optical depths at shorter wavelengths and lower optical depths at longer wavelengths.

Zoom In Zoom Out Reset image size

It is possible that we slightly overpredict A_V at high redshifts. As we discuss in further detail in Section 4, our approximation of spherical symmetry could lead to overestimates of the optical depth, particularly for massive galaxies, in which the bulk of interstellar dust is typically observed to reside in the disk. The rise in A_V at high redshift is mostly driven by ${R}_{\mathrm{dust}}$, since τ_ν ∝ 1/${R}_{\mathrm{dust}}$² (Equation (6)), and since ${R}_{\mathrm{dust}}$ ∝ (1 + z)⁻¹ (Equations (7) and (8); Figure 2). Thus, with our assumption of spherical geometry, and given that we have defined the optical depth as the integrated value through to the galactic center (Equation (4)), the values we have derived for A_V are most likely upper limits.

We use our results for ${M}_{\mathrm{dust}}$ and τ_ν in the previous sections to determine ${T}_{\mathrm{dust}}$ from Equation (2). In Figures 7 and 8, we present predictions for the evolution of galaxy dust temperature as a function of stellar mass. We show results for the two galaxy geometries we consider here, the "point source" and "homogeneous" models illustrated in Figure 1. We also show results for two different surface area covering fractions of dust, f_⋆ = 1 and f_⋆ = 0.25. Both models show that the ${T}_{\mathrm{dust}}$–${M}_{\star }$ relation is not monotonic, but rather peaks at characteristic values of ${M}_{\star }$. In both figures, we overplot data points of measured galaxy dust temperatures from the literature.

Zoom In Zoom Out Reset image size

The point source model (Figures 7) shows that for most galaxy stellar masses (${M}_{\star }$ ≳ 10^6.5 ${M}_{\odot }$), ${T}_{\mathrm{dust}}$ tends to increase over time, and is a poor representation of the observations. Whereas observations suggest that dust temperatures in higher-mass galaxies tend to get cooler with time, the point source model predicts the opposite. By contrast, the model in which stars and the ISM are homogeneously distributed (Figures 8) suggests that ${T}_{\mathrm{dust}}$ tends to decrease with time. At all redshifts, our homogeneous, f_⋆ = 0.25 model is in better agreement with the observations by Rémy-Ruyer et al. (2015) and da Cunha et al. (2015) than the f_⋆ = 1 model. We take the former to be our fiducial model; Table 2 lists the values of ${T}_{\mathrm{dust}}$ for a 10¹² ${M}_{\odot }$ halo mass galaxy in this model. That the f_⋆ = 0.25 model is in better agreement with the observations than the f_⋆ = 1 model is not surprising. This implies that galactic dust does not have a covering fraction of 100% and that there are regions in any given galaxy where starlight escapes without being absorbed by dust. The actual covering fraction is almost certain to vary widely from galaxy to galaxy.

Zoom In Zoom Out Reset image size

Focusing now on Figure 8, starting at about z = 6, the dust temperature tends to cool down for galaxies of fixed stellar masses, as the universe evolves. At each redshift, there are two peaks in the ${T}_{\mathrm{dust}}$–${M}_{\star }$ relation. As the redshift decreases, the first peak shifts toward higher and higher values of ${M}_{\star }$. For instance, between z = 6 and z = 3, the stellar mass at which ${T}_{\mathrm{dust}}$ peaks shifts from ${M}_{\star }$ = 10^6.6 to 10^7.7 ${M}_{\odot }$. The peak near the high-mass end is more stable and does not display a similar systematic shift over time.

Our model predicts a population of high-redshift (z ≳ 2), low-mass galaxies (${M}_{\star }$ ≈ 10⁶ to 10^8.5 ${M}_{\odot }$) with fairly hot dust. From z = 6 to z = 4, these galaxies have dust that is around the same temperature as—if not hotter than—the dust in their more massive counterparts at the same redshifts. Unfortunately, there are no observational constraints on the SFH of such galaxies, so with the current state of knowledge, indirect methods would have to be used to infer the robustness of the peak temperatures our model predicts.

In Figure 9 we plot dust temperatures as a function of redshift for fixed stellar masses. We find that dust temperatures of galaxies of all stellar masses evolve markedly with redshift. Starting the present era to z = 6, galaxies having stellar masses from 10⁸ to 10¹⁰ ${M}_{\odot }$ have dust temperatures that increase monotonically from about 25–35 to ∼55–60. Higher-mass galaxies first display increases in ${T}_{\mathrm{dust}}$, followed by decreases in the dust temperature as they evolve to higher redshifts.

Zoom In Zoom Out Reset image size

Viero et al. (2013) stacked Herschel images of a stellar-mass-selected sample of galaxies, and Béthermin et al. (2015) performed a similar stacking analysis of Spitzer, Herschel, LABOCA, and AzTEC data for galaxies taken from the COSMOS field. These authors found that the dust temperature tends to increase with redshift, up to z = 6, for galaxies of all stellar masses. While we find a similar trend for galaxies having stellar masses ${M}_{\star }$ < 10¹⁰${M}_{\odot }$, our results are at odds with the observations for massive galaxies. That the massive galaxies in our model do not show a steady increase in ${T}_{\mathrm{dust}}$ with redshift is primarily a reflection of our geometric model for the ISM. As discussed in Section 3.2, our approximation of spherical symmetry may result in overestimates of τ_ν, especially for high-redshift, massive galaxies, where most interstellar dust is observed in the disk. Overpredicting τ_ν for high-mass galaxies naturally leads to the non-monotonic evolution of ${T}_{\mathrm{dust}}$  observed in Figure 9. Nevertheless, the values we derive for ${T}_{\mathrm{dust}}$ at any fixed redshift are in general agreement with the observational results of Viero et al. (2013) and Béthermin et al. (2015), and with the theoretical models of Cowley et al. (2017). At any given redshift, we predict temperatures of a factor of roughly 1.5 higher than these authors. This systematic offset may partly be due to the choice of modified blackbody models fitted by these authors and partly a result of how they selected sources. Indeed, other authors who performed stacking analyses but fitted different models or had different selection criteria, including Pascale et al. (2009), Amblard et al. (2010), and Elbaz et al. (2010), reported higher dust temperatures in alignment with our results. We note that these latter two studies, which are based on galaxy samples selected by Herschel, may be biased in temperature. Due to the varying sensitivity of the PACS instrument with wavelength, hotter galaxies are easier to detect.

In light of recent ALMA observations that have detected large amounts of dust in galaxies during the epoch of reionization (e.g., Watson et al. 2015; Laporte et al. 2017), we are motivated to make predictions of the flux density due to dust in galaxies. We assume that the dust in a galaxy, of total mass ${M}_{\mathrm{dust}}$, will rise to an average temperature ${T}_{\mathrm{dust}}$, due to heating by starlight, and the dust will re-emit most of the light in the infrared. The galaxy will have a flux density of

$$\begin{eqnarray}&&{S}_{\nu }=\displaystyle \frac{{\kappa }_{\nu }{B}_{\nu }({T}_{\mathrm{dust}}){M}_{\mathrm{dust}}(1+z)}{{d}_{L}^{2}},\end{eqnarray} \tag{ 17 }$$

where B_ν(${T}_{\mathrm{dust}}$) is the Planck function, and d_L is the luminosity distance to the galaxy. The emission is assumed to be optically thin.

Figure 10 presents the galaxy flux density at an observed wavelength of 1.1 mm, as a function of ${M}_{\star }$ for z = 0 to z = 9.5. For fixed redshifts, S_ν increases with galaxy stellar mass. For fixed stellar masses, for sources at redshifts z ≳ 1, S_ν decreases with time. This trend is due to the combination of two effects: (1) the negative K-correction (e.g., Hughes et al. 1998; Barger et al. 1999; Blain et al. 2002; Lagache et al. 2005); and (2) our fiducial model predicts that ${T}_{\mathrm{dust}}$ increases with z. At redshifts z > 1, the far-infrared radiation emitted in distant galaxies is redshifted to submillimeter wavelengths, and the resulting negative K-correction counteracts the dimming of galaxies caused by their cosmological distances. If more distant galaxies have hotter dust, then the observed flux from these galaxies originated from radiation emitted closer to the peak of the blackbody radiation curve than nearby galaxies. For instance, photons emitted from dust in a galaxy at z = 9.5 had rest wavelengths of λ₀ = 116 μm, while photons emitted from a source at z = 2 had λ₀ = 550 μm. From Figure 8, one can see that the dust temperature of a 10⁸ ${M}_{\odot }$ stellar mass galaxy at z = 9.5 is ∼74 K, corresponding to a peak in the blackbody radiation curve at 39 μm. A 10⁸ ${M}_{\odot }$ galaxy at z = 2 has ${T}_{\mathrm{dust}}$ ≈ 45 K, corresponding to 64 μm. Thus, the observed flux at 1.1 mm from the more distant galaxy at z = 9.5 originated from dust whose emission was closer to the peak of the blackbody radiation curve, compared to the galaxy at z = 2.

Zoom In Zoom Out Reset image size

In Figure 11, we show results in terms of S_ν as a function of z, for fixed stellar masses: ${M}_{\star }$ = 10⁸, 10⁹, and 10¹⁰ ${M}_{\odot }$. For a given stellar mass, S_ν first decreases with time, then it rises sharply from z = 1 to z = 0. We overplot the observed ∼1 mm continuum fluxes of galaxies recently observed with ALMA. Watson et al. (2015) observed the lensed galaxy A1689-zD1 between 1.2 and 1.4 mm and detected a flux of S_ν =0.61 ± 0.12 mJy. Located at z ≈ 7.5 and magnified by a factor of 9.3, A1689-zD1 has a stellar mass, dust mass, and SFR of ${M}_{\star }$ ≈ 2 × 10⁹ ${M}_{\odot }$, ${M}_{\mathrm{dust}}$ ≈ 4 × 10⁷ ${M}_{\odot }$, and SFR ≈ 9 ${M}_{\odot }$ yr⁻¹. Laporte et al. (2017) observed the lensed galaxy A2744_YD4 at 0.84 mm. At a redshift of about 8.4 and magnified by a factor of ∼1.8, A2744_YD4 has a stellar mass, dust mass, and SFR of ${M}_{\star }$ ≈ 2 × 10⁹ ${M}_{\odot }$, ${M}_{\mathrm{dust}}$ ≈ 6 × 10⁶ ${M}_{\odot }$, and SFR ≈ 20 ${M}_{\odot }$ yr⁻¹. Dunlop et al. (2017) conducted the first, deep ALMA image of the Hubble Ultra Deep Field, detecting 16 sources at 1.3 mm. The sources have high stellar masses, with 13 out of 16 having ${M}_{\star }$ > 10¹⁰ ${M}_{\odot }$. Fifteen of the sources are located at redshifts 0.7 ≤ z ≤ 3, and one source is located at z = 5. The observed fluxes, and other properties of all the sources plotted in Figures 10 and 11, are summarized in Table 3.

Zoom In Zoom Out Reset image size

Table 3.  Observations Plotted in Figure 11

Galaxy Name	Redshift	Sν	${M}_{\star }$	SFRIR	References
		(mJy)	(109 ${M}_{\odot }$)	(${M}_{\odot }$ yr−1)
A1689-zD1	7.5 ± 0.2	0.61 ± 0.12	${1.7}_{-0.5}^{+0.7}$	${9}_{-2}^{+4}$	1
A2744_YD4	${8.38}_{-0.11}^{+0.13}$	0.099 ± 0.023	${1.97}_{-0.66}^{+1.45}$	${20.4}_{-9.5}^{+17.6}$	2
UDF1	3.00	0.924 ± 0.076	${50}_{-10}^{+13}$	326 ± 83	3
UDF2	2.79	0.996 ± 0.087	${126}_{-37}^{+52}$	247 ± 76	3
UDF3	2.54	0.863 ± 0.084	${20}_{-6}^{+8}$	195 ± 69	3
UDF4	2.43	0.303 ± 0.046	${32}_{-9}^{+13}$	94 ± 4	3
UDF5	1.76	0.311 ± 0.049	${25}_{-7}^{+10}$	102 ± 7	3
UDF6	1.41	0.239 ± 0.049	${32}_{-7}^{+8}$	87 ± 11	3
UDF7	2.59	0.231 ± 0.048	${40}_{-8}^{+10}$	56 ± 22	3
UDF8	1.55	0.208 ± 0.046	${159}_{-46}^{+65}$	149 ± 90	3
UDF9	0.67	0.198 ± 0.039	${10}_{-2}^{+3}$	23 ± 25	3
UDF10	2.09	0.184 ± 0.046	${16}_{-5}^{+7}$	45 ± 22	3
UDF11	2.00	0.186 ± 0.046	${6}_{-13}^{+16}$	162 ± 94	3
UDF12	5.00	0.154 ± 0.040	${4}_{-1}^{+2}$	37 ± 14	3
UDF13	2.50	0.174 ± 0.045	${63}_{-13}^{+16}$	68 ± 18	3
UDF14	0.77	0.160 ± 0.044	${5}_{-1}^{+1}$	44 ± 17	3
UDF15	1.72	0.166 ± 0.046	${8}_{-2}^{+3}$	38 ± 27	3
UDF16	1.31	0.155 ± 0.044	${79}_{-16}^{+21}$	40 ± 18	3

References. (1) Watson et al. (2015), (2) Laporte et al. (2017), (3) Dunlop et al. (2017).

Download table as:  ASCIITypeset image

To date, observations of the dust emission in sources at z ≳ 1 have been restricted to galaxies having stellar masses ≳10⁹ ${M}_{\odot }$. To achieve the sensitivities necessary to detect the dust emission of single, high-redshift L^⋆ galaxies at these masses requires total observing times of ∼2–3 hr with ALMA (e.g., Laporte et al. 2017). Since typical lower-mass galaxies (${M}_{\star }$ < 10⁹ ${M}_{\odot }$) are intrinsically fainter, the much longer observing times needed to detect their dust emission at z > 1 may be prohibitive. Yet, serendipitous occurrences, such as gravitational lensing, could possibly aid in the detection and characterization of the low-mass, star-forming population of galaxies in the early universe. If and when low-mass galaxies begin to be detected in large numbers, it may be easier (e.g., less time-consuming) to first detect galaxies at higher redshifts, since according to our model, their observed millimeter flux is expected to exceed that of lower redshift galaxies by ∼1 to 2 orders of magnitude.

In Section 2 we model galaxies as spherically symmetric with either most of the dust concentrated around the nucleus or else homogeneously mixed with gas and stars. In reality, dust is often consolidated in the disk of large galaxies, so the assumption of spherical symmetry may result in overestimates of the optical depth for these systems.

In Equation (4), we assume that galactic dust has a simple power-law density distribution, with γ = 0. The assumption that the dust profile is constant with radius may seem unfounded, given observations that dust content varies with galactic radius (e.g., Boissier et al. 2005; Muñoz-Mateos et al. 2009). Yet appropriate values for γ and ${R}_{\mathrm{dust}}$, of which we are equally ignorant, are certain to vary significantly from galaxy to galaxy. Since there are infinite combinations of γ and ${R}_{\mathrm{dust}}$ that produce identical values of τ_ν—that is, since γ and ${R}_{\mathrm{dust}}$ are essentially degenerate—we decide to absorb our ignorance about both quantities into our definition of ${R}_{\mathrm{dust}}$ in Equation (4), where we let ${R}_{\mathrm{dust}}$ = ${R}_{\mathrm{gal}}$, (recalling that ${R}_{\mathrm{gal}}$ is the half-mass galactic radius). For example, for two galaxies with identical values of κ_ν and ${M}_{\mathrm{dust}}$, τ_ν for one galaxy with ${R}_{\mathrm{dust}}$ = ${R}_{\mathrm{gal}}$ and γ = 0 is equivalent to τ_ν for the second galaxy with γ = 1/3 and ${R}_{\mathrm{dust}}$ = 2 ${R}_{\mathrm{gal}}$.

Another source of uncertainty in our model is the choice of spin parameter, λ, which is expected to link galaxy disk size and halo size (see Equation (7)), under the assumption that the collapsing baryonic matter of a galaxy inherits the same specific angular momentum as the halo (Fall & Efstathiou 1980). In Section 2.1.1, we assume that the half-mass–radius and dark matter halo radius, for all galaxy masses at all redshifts, are linked by a single value, ${R}_{\mathrm{gal}}$/R_halo = 0.018, corresponding to λ = 0.036 (Somerville et al. 2018). Some simulations suggest that there is significant scatter—about two orders of magnitude—about the mean value of λ (Teklu et al. 2015; Zavala et al. 2016), and that this scatter depends in part on galaxy morphology (e.g., Teklu et al. 2015). The results of Somerville et al. (2018), who demonstrate that λ is roughly independent mass and weakly evolves with redshift, are in general agreement with other recent studies of the relationship between galaxy and halo size (e.g., Kawamata et al. 2015; Shibuya et al. 2015; Huang et al. 2017). Yet if our adopted value of λ leads to underestimates or overestimates of ${R}_{\mathrm{gal}}$ for galaxies at certain masses and epochs, these inaccuracies will naturally propagate into our estimates of ${R}_{\mathrm{dust}}$ and τ_ν.

We use a prescription for ${M}_{\mathrm{gas}}$(${M}_{\star }$, z) determined by Zahid et al. (2014), who determined a relation between metallicity and stellar-to-gas mass ratio for galaxies at z ≲ 1.6 and with ${M}_{\star }$ ≳ 10⁹ ${M}_{\odot }$. We extrapolate this relation for higher redshifts and lower stellar masses. If the interstellar environments of low-mass galaxies manifest in significantly different relationships between ${ \mathcal Z }$, ${M}_{\star }$, and ${M}_{\mathrm{gas}}$, and if ISM conditions evolve with redshift, then the Zahid et al. relation may break down in unexpected ways in the low-${M}_{\star }$, high-z regimes. Further observational tests are required to confirm or rule out the universal metallicity relation upon which Equation (12) is based. Nevertheless, it is promising that the Zahid et al. relation is consistent with that of Andrews and Martini (2013), who measure the MZ relation down to ${M}_{\star }$ ≈ 10^7.5 ${M}_{\odot }$.

The predicted dust masses are sensitive to the functional form of ${M}_{\mathrm{gas}}$, (Equation (11)). To give a sense of how the evolution of ${M}_{\mathrm{gas}}$ affects ${M}_{\mathrm{dust}}$, we present additional calculations for ${M}_{\mathrm{dust}}$, using an alternative prescription for ${M}_{\mathrm{gas}}$. Sargent et al. (2014) compiled a sample of 131 massive (${M}_{\star }$ > 10¹⁰ ${M}_{\odot }$), star-forming galaxies at redshifts z ≲ 3. They derived a Schmidt-Kennicutt relation:

$$\begin{eqnarray}\mathrm{log}\left(\displaystyle \frac{{M}_{\mathrm{mol}}}{{M}_{\odot }}\right) & = & (9.22\pm 0.02)+(0.81\pm 0.03)\\ & & \times \mathrm{log}\left(\displaystyle \frac{\mathrm{SFR}}{{M}_{\odot }\,{\mathrm{yr}}^{-1}}\right),\end{eqnarray} \tag{ 18 }$$

where M_mol is the galactic molecular mass. Equation (18) is appealing as a comparison to the Zahid et al. (2014) formulation for ${M}_{\mathrm{gas}}$, because it is directly applicable for galaxies at higher redshifts. We do not attempt to estimate the total gas mass from Equation (18) but calculate the dust mass using ${M}_{\mathrm{gas}}$ = M_mol × DGR. Using the Sargent et al. (2014) formulation for the molecular mass of galaxies, in Figures 12 and 13 we display plots of ${M}_{\mathrm{dust}}$ as a function of ${M}_{\star }$, analogous to Figures 3 and 5. Figure 12 shows that while the Zahid et al. (2014) formulation for ${M}_{\mathrm{gas}}$ results in higher predictions for ${M}_{\mathrm{dust}}$ for most stellar masses, for galaxies with ${M}_{\star }$ ≳ 10^9.5 ${M}_{\odot }$, the Zahid et al. and Sargent et al. formulations come into better agreement. Not too surprisingly, using Equation (18) leads to a model for ${M}_{\mathrm{dust}}$ that underpredicts observed values in high-mass galaxies (Figure 13). This is because we did not account for the total galactic gas mass here, only M_mol. However, ${M}_{\mathrm{dust}}$ calculated using Equation (18) is in good agreement with low-mass galaxies with ${M}_{\star }$ ≲ 10⁸ ${M}_{\odot }$. Moreover, for redshifts z ≳ 1, ${M}_{\mathrm{dust}}$ calculated using Equation (18) comes into better agreement with the observations of high-mass galaxies, since the gas fraction in galaxies is expected to be increasingly dominated by molecular gas as redshift increases.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The total gas mass fraction, defined as f_g,tot = ${M}_{\mathrm{gas}}$/(${M}_{\mathrm{gas}}$ + ${M}_{\star }$), is the subject of a number of studies. For stellar masses in the range 10¹⁰–4 × 10¹¹ ${M}_{\odot }$, the Zahid et al. relation predicts f_g,tot ≈ 0.28, 0.40, 0.48, and 0.55 for redshifts z = 1, 2, 3, and 4, respectively, though the authors caution the extrapolation of their relation beyond z = 1.6 (private communication). These values are quite consistent with the molecular gas fractions, f_g,mol = M_mol/(M_mol + ${M}_{\star }$), derived in many studies from CO and dust observations of galaxies having a similar range of stellar masses. Daddi et al. (2010) measured f_g,mol ≈ 0.6 for six galaxies with ${M}_{\star }$ =0.33–1.1 × 10¹¹ ${M}_{\odot }$ at z = 1.5. Tacconi et al. (2010) studied 19 galaxies with ${M}_{\star }$ = 0.3–3.4 × 10¹¹ ${M}_{\odot }$. They measured f_g,mol = 0.2–0.5 at z ≈ 1.1 and f_g,mol = 0.3–0.8 at z ≈ 2.3. Later on, Tacconi et al. (2013) found similar results with a larger sample of 52 galaxies, measuring average molecular gas fractions of 0.33 and 0.47 at z ∼ 1.2 and 2.2. Magdis et al. (2012) measured f_g,mol ≈ 0.36 for a ${M}_{\star }$ =2 × 10¹¹ ${M}_{\odot }$ galaxy at z = 3.21. Saintonge et al. (2013) measured 0.45 for ${M}_{\star }$ ∼ 10¹⁰ ${M}_{\odot }$ galaxies at z = 2.8. More recently, Béthermin et al. (2015) used observations of dust emission in massive (∼6 × 10¹⁰ ${M}_{\odot }$) galaxies to measure f_g,mol = 0.16–0.35 at z < 1, 0.27–0.41 at 1 < z < 2, ∼0.5 at 2 < z < 3, and ∼0.6 at 3 < z < 4. And in an extensive study of 145 galaxies from the COSMOS survey, Scoville et al. (2016) measured molecular gas fractions of 0.16–0.67 at z ≈ 1.5, 0.24–0.75 at z ≈ 2.2, and 0.23–0.85 at z ≈ 4.4.

While it is generally agreed that high-redshift galaxies are gas-dominated, the details of the redshift evolution of f_g,mol are uncertain. By re-expressing the gas fraction (total or molecular) as $1/[1+{({t}_{\mathrm{dep}}\mathrm{sSFR})}^{-1}]$, one can see its dependence on the gas depletion time t_dep and the specific star formation rate (sSFR), both of which are expected to be redshift-dependent quantities. For instance, some studies suggest that the typical sSFR of main-sequence galaxies reaches a plateau by z ∼ 2 (e.g., González et al. 2010; Rodighiero et al. 2010; Weinmann et al. 2011), which would result in lower gas fractions than if the sSFR steadily increased beyond z = 2, as suggested by other studies (e.g., Bouwens et al. 2012; Stark et al. 2013). Thus, for example, if in this study we underestimated the amount of gas in galaxies, this would lead to underestimating ${M}_{\mathrm{dust}}$. For a fixed stellar mass, more dust means that the quantity of UV photons per unit mass of dust would be lower, resulting in a decreased dust temperature. It turns out that ${T}_{\mathrm{dust}}$ in our model is fairly robust to changes in the gas and dust mass. For instance, if ${M}_{\mathrm{dust}}$ were higher by a factor of two for all galaxy stellar masses at all redshifts, this would decrease the values of ${T}_{\mathrm{dust}}$ reported here by a factor of only ∼0.9.

In Section 2.1.2, we described how we use Galactic laws to model the opacity κ_ν. In particular, we used the Beckwith et al. (1990) relation for long wavelengths and discussed how changes in the normalization or slope of κ_ν could affect the resulting dust temperatures. We performed a series of tests to quantify these changes and found that changing β by ±1 affects the average ${T}_{\mathrm{dust}}$ by a factor of only about 1.3. On the other hand, varying the normalization by a factor of 2 affects $\langle {T}_{\mathrm{dust}}\rangle$ by a factor of only 1.1, on average.

While our results are robust to the opacity model, ${T}_{\mathrm{dust}}$ is more sensitive to the prescription for galactic metallicity and the dust-to-gas ratio. For the relationship between metallicity and DGR, we adopt the prescription of Rémy-Ruyer et al. (2014) determined from observations of galaxies at z = 0. For galaxies with ${ \mathcal Z }\gt 0.26\,{{ \mathcal Z }}_{\odot }$, Equation (14) states $\mathrm{DGR}/{\mathrm{DGR}}_{\odot }={ \mathcal Z }/{{ \mathcal Z }}_{\odot }$. This assumption is likely to break down for high-redshift, young galaxies, where the dust production sites have not yet reached equilibrium, especially if dust production by AGB stars is important (e.g., Dwek & Cherchneff 2011). Thus, such high-redshift galaxies would have higher DGRs than predicted by Equation (14), which would translate into higher dust masses.