Evolution of Interstellar Medium, Star Formation, and Accretion at High Redshift

Galaxy evolution in the early universe is dominated by three processes: the conversion of interstellar gas into stars, the accretion of intergalactic gas to replenish the interstellar gas reservoir, and the merging of galaxies. The latter redistributes the gas within the galaxies, promotes starburst activity, fuels active galactic nuclei (AGNs), and transforms the stellar morphology from disk-like (rotation dominated) to ellipsoidal systems. In all of these processes, the interstellar gases play a determining role—including in the outcome of galaxy mergers, since the gas is dissipative and becomes centrally concentrated, fueling starbursts and AGNs.

At present, the gas supply and its evolution at high redshifts are only loosely constrained (see the reviews Solomon & Vanden Bout 2005; Carilli & Walter 2013). Only ≲200 galaxies at z > 1 have been observed in the CO lines (and most not in the CO (1–0) line, which has been calibrated as a mass tracer). To properly chart and understand the evolution, large samples are required, probing multiple characteristics: (1) the variation with redshift or cosmic time, (2) the dependence on galaxy mass, and (3) the differences between the galaxies with "normal" SF activity on the main sequence (MS) and the starbursts (SB). The latter constitute only ∼5% of the population but have star formation rates (SFRs) elevated to 2–100 times higher levels than the MS. The contribution of the SB galaxies to the total SF at z < 2 is 8%–14% (Sargent et al. 2012). Although the high-z galaxies above the MS will here be referred to as starbursts, nearly all high-z SF galaxies would be classified as SB galaxies if they were at low redshift.

Properly constraining the evolution of the gas contents would require time-consuming CO line observations spanning and populating the full ranges of redshift, ${M}_{* }$, and MS versus SB populations. In addition, translating the observed excited-state CO emission fluxes into reliable estimates of the gas contents remains a problem (Carilli & Walter 2013).

As an alternative, we here develop a formulation for the high-redshift galaxy evolution, using very extensive observations of the Rayleigh–Jeans (RJ) dust continuum emission and applying the calibration of this technique developed in Scoville et al. (2016; a brief summary is provided here in Section 4 and Appendix A; see also Magnelli et al. 2014; Santini et al. 2014; Genzel et al. 2015; Berta et al. 2016; Schinnerer et al. 2016). The calibration used here is based on observations of the RJ dust emission and CO (1–0) emission in low-redshift galaxies; the technique provides roughly a factor of 2 accuracy in the derived interstellar medium (ISM) masses—provided one restricts the galaxy samples to high stellar mass galaxies ($\gt {10}^{10}\,{M}_{\odot }$), which are expected to have near-solar gas-phase metallicity (e.g., Erb et al. 2006). Such galaxies are expected to have gas-to-dust abundance ratios and dust properties similar to those for galaxies in our calibration sample at low redshift (Draine et al. 2007).

The long wavelength dust emission is nearly always optically thin and is only mildly susceptible to "excitation" variations, since the RJ emission depends linearly on the mass-weighted dust temperature, and there is no dependence on the volume density of the gas and dust. In fact, the mass-weighted ISM dust temperatures are likely to vary less than a factor of 2 (Section 4.3; see also Scoville et al. 2016). We note also that, in contrast to the CO, which typically constitutes just ∼0.1% of the gas mass, the dust is a 1% mass tracer. CO is also susceptible to depletion by photodissociation in strong UV radiation fields, whereas the dust abundance is less sensitive to the ambient radiation field.

The galaxy sample used here has ALMA continuum observations in Band 6 (240 GHz) and Band 7 (345 GHz); they are all within the COSMOS survey field and thus have excellent ancillary data. The ALMA pointings are noncontiguous, but their fields of view (FOVs), totaling 102.9 arcmin², include 708 galaxies measured at far-infrared wavelengths by Herschel. All of the Herschel sources within the ALMA FOVs are detected by ALMA, so this is a complete sampling of the IR-bright galaxies at z = 0.3 to 4.5. Given the positional prior from the COSMOS Herschel catalog, all sources are detected in the far-infrared; the sample therefore has reliable estimates of the dust-embedded SFR activity. The dusty SF activity is in virtually all cases dominant (5–10 times) over the unobscured SF probed in the optical/UV.

The major questions we address are as follows:

• 1.  
  How does the gas content depend on the stellar mass of the galaxies?
• 2.  
  How do these gas contents evolve with cosmic time, down to the present, where they are typically less than 10% of the galactic mass?
• 3.  
  In the starburst populations, is the prodigious SF activity driven by increased gas supply or increased efficiency of converting the existing gas into stars?
• 4.  
  How rapidly is the ISM being depleted? The depletion timescale is characterized by the ratio ${M}_{\mathrm{ISM}}$/SFR, but to date this has not been measured in broad samples of galaxies due to the paucity of quantitative ISM measurements at high redshift.
• 5.  
  If the ISM depletion times are as short as they appear to be (several ×10⁸ years), yet the SF galaxy population persists for a much longer span of cosmic time (≳5 × 10⁹ years), then there must be large rates of accretion of new gas from the halo or intergalactic medium (IGM) to maintain the activity. At present there are virtually no observational constraints on these accretion rates, only theoretical predictions (Dekel et al. 2013), so a major unknown is, what are these accretion rates?
• 6.  
  How does the efficiency of star formation from a given mass of ISM (the so-called gas-depletion timescale) vary with cosmic epoch, the stellar mass of the galaxy, and whether the galaxy is on the MS or undergoing a starburst?

1.1. Outline of This Work

The bulk of the SF galaxy population can be located on an MS locus in the plane of SFR versus ${M}_{* }$ (Noeske et al. 2007b; Peng et al. 2010; Rodighiero et al. 2011; Whitaker et al. 2012; Lee et al. 2015).¹⁹

2.1. Evolution of the MS with Redshift

The MS locus (SFR(z, M_*)) evolves to higher SFRs at earlier epochs as shown in Figure 1, left (Speagle et al. 2014, their favored model #49). This MS definition was based on an extensive reanalysis of all previous work; similar MS definitions have been obtained by Béthermin et al. (2012), Whitaker et al. (2014), Schreiber et al. (2015), Lee et al. (2015), and Tomczak et al. (2016). Model #49 is expressed analytically by

$$\begin{eqnarray*}&&{\mathrm{SFR}}_{\mathrm{MS}}=\,{10}^{(0.59\mathrm{log}({M}_{* })-5.77)}\times {(1+z)}^{(0.22\mathrm{log}({M}_{* })+0.59)}\,.\end{eqnarray*}$$

We use this function, evaluated at log(${M}_{* })=10.5\,{M}_{\odot }$, to describe the redshift evolution of the MS:

$$\begin{eqnarray}&&{\zeta }_{\mathrm{MS}}(z)\equiv \,{\mathrm{SFR}}_{\mathrm{MS}}(z)/{\mathrm{SFR}}_{\mathrm{MS}}(z=0)=\,{(1+z)}^{2.9}\,.\end{eqnarray} \tag{ 1 }$$

To fit the redshift dependence of ISM masses in Section 6.1, we fit for an exponent of (1 + z) so that we can directly compare the evolution of the ISM masses with that of the SFRs that have a ${(1+z)}^{2.9}$ dependence.

Zoom In Zoom Out Reset image size

2.2. Shape of the MS with M_*

Although the early descriptions of the MS used single power laws as a function of stellar mass, the more recent work (Whitaker et al. 2014; Lee et al. 2015; Tomczak et al. 2016) indicated a break in the slope of the MS at ${M}_{* }\gt 5\times {10}^{10}\,{M}_{\odot }$, having lower sSFR at higher ${M}_{* }$. Here, we use the shape of the MS from Lee et al. (2015) at z = 1.2 for the stellar mass dependence of the MS,

$$\begin{eqnarray*}&&\,{\mathrm{SFR}}_{\mathrm{MS}}=\,{10}^{(1.72-\mathrm{log}(1+{({10}^{(\mathrm{log}{M}_{* }-10.31)})}^{-1.07}))}\,,\end{eqnarray*}$$

and normalize to unity at log ${M}_{* }=10.5\,{M}_{\odot }$:

$$\begin{eqnarray}&&{\zeta }_{\mathrm{MS}}({M}_{* })\equiv \,{\mathrm{SFR}}_{\mathrm{MS}}({M}_{* })/{\mathrm{SFR}}_{\mathrm{MS}}(\mathrm{log}{M}_{* }=10.5)\,.\end{eqnarray} \tag{ 2 }$$

Figure 1 right shows the adopted shape function ${\zeta }_{\mathrm{MS}}$(${M}_{* }$). The normalization constant is SFR = 3.23 $\,{M}_{\odot }$ yr⁻¹ at z = 0 and log M${}_{* }=10.5$ $\,{M}_{\odot }$.

We have investigated the use of alternative MS formulations, and the conclusions derived here do not depend qualitatively on the adopted MS formulation, although the numerical values of the power-law exponents can change by $\sim \pm 0.1$ with adoption of one of the other MS definitions. We used the combination of Speagle et al. (2014) and Lee et al. (2015), since the former covers the complete redshift range covered here (but has only power-law dependence on mass), whereas the latter has the break in the MS at log(${M}_{* })\sim 10.5\,{M}_{\odot }$ that is seen in the latest determinations of the MS.

2.3. Continuity of the MS Evolution

In our analysis, we make use of a principle we refer to as the continuity of MS evolution: simply stated, the temporal evolution of the SF galaxy population may be followed by Lagrangian integration of the MS galaxy evolution. This follows from the fact that approximately 95% of the SF galaxies at each epoch lie on the MS with SFRs dispersed only a factor of 2 above or below the MS (e.g., Rodighiero et al. 2011). (A similar approach has been used by Noeske et al. (2007a), Renzini (2009), and Leitner (2012) and references therein.)

This continuity assumption ignores the galaxy buildup arising from major mergers of similar-mass galaxies since they can depopulate the MS population in the mass range of interest. In fact, the major mergers may be responsible for some of the 5% of galaxies in the SB population above the MS (see Section 9). On the other hand, minor mergers may be considered simply as one element of the average accretion process considered in Section 8.

We are also neglecting the SF quenching processes in galaxies. This occurs mainly in the highest-mass galaxies (${M}_{* }\gt 2\times {10}^{11}\,{M}_{\odot }$  at $z\gt 1$) and in dense environments at lower redshift (Peng et al. 2010; Darvish et al. 2016). At $z\gt 1.2$ the quenched red galaxies are a minor population (see Figure 14), and the quenching processes are of lesser importance.

The paths of evolution in the SFR versus ${M}_{* }$ plane can be easily derived since the MS loci give ${{dM}}_{* }/{dt}$ = SFR (${M}_{* }$). One simply follows each galaxy in a Lagrangian fashion as it builds up its mass.

Using this continuity principle to evolve each individual galaxy over time, the evolution for MS galaxies across the SFR–${M}_{* }$ plane is as shown in Figure 2. Here we have assumed that 30% of the SFR is eventually put back into the ISM by stellar mass loss. This is appropriate for the mass loss from a stellar population with a Chabrier IMF (see Leitner & Kravtsov 2011). In this figure, the curved horizontal lines are the MS at fiducial epochs or redshifts, while the downward curves are the evolutionary tracks for fiducial ${M}_{* }$ from 1 to 10 × 10¹⁰ $\,{M}_{\odot }$. At higher redshift, the evolution is largely toward increasing ${M}_{* }$, whereas at lower redshift the evolution is in both SFR and ${M}_{* }$. In future epochs, the evolution is likely to be still more vertical as the galaxies exhaust their gas supplies. Thus there are three phases in the evolution:

• 1.  
  the gas accretion dominated and stellar mass buildup phase at $z\gt 2$ corresponding to a cosmic age less than 3.3 Gyr (see Section 8),
• 2.  
  the transition phase where gas accretion approximately balances SF consumption and the evolution becomes diagonal, and
• 3.  
  the epoch of ISM exhaustion at z ≲ 0.1 (age 12.5 Gyr) where the evolution will be vertically downward in the SFR versus ${M}_{* }$ plane.

Zoom In Zoom Out Reset image size

These evolutionary phases are all obvious (and not a new development here), but in Section 8 we make use of the continuity principle to derive the accretion rates and hence substantiate the three phases as separated by their accretion rates relative to their SFRs. When these phases begin and end is a function of the galaxy stellar mass—the transitions in the relative accretion rates take place much earlier for the more massive galaxies.

2.4. MS Galaxy Evolution Paths

2.5. Starbursts

3.1. ALMA Band 6 and 7 Continuum Data

Within the COSMOS survey field, extensive observations from ALMA now exist for the dust continuum of high-redshift galaxies. Here, we make use of all those data that are publicly available, in addition to our own still-proprietary observations. The 18 projects are listed in Table 1 together with summary listings of the observed bands, the number of pointings, and the average frequency of observation, the synthesized beam sizes, and the typical rms noise. The total number of pointings in these data sets is 1011, covering a total area 0.0286 deg² or 102.9 arcmin² within the half-power beam width (HPBW).

Table 1.  ALMA Data Sets Used for Continuum Measurements

Project Code	PI	Band	$\langle {\nu }_{\mathrm{obs}}\rangle$	$\#$ pointings	$\langle \mathrm{HPBW}\rangle$	${\sigma }_{\mathrm{rms}}$
			(GHz)		(arcsec)	(mJy)
2011.0.00097.S	Scoville	7	341	105	0.64	0.300
2011.0.00964.S	Riechers	7	296	1	0.69	0.091
2012.1.00076.S	Scott	6	244	20	1.27	0.062
2012.1.00323.S	Popping	6	226	2	1.01	0.093
2012.1.00523.S	Capak	6	291	3	0.74	0.041
"		7	298	5	0.77	0.039
2012.1.00978.S	Karim	7	338	6	0.37	0.084
2013.1.00034.S	Scoville	6	245	60	0.66	0.063
"		7	343	98	0.51	0.141
2013.1.00118.S	Aravena	6	240	128	1.31	0.124
2013.1.00151.S	Schinnerer	6	240	86	1.71	0.075
2013.1.00208.S	Lilly	7	343	15	1.20	0.129
2013.1.00276.S	Martin	6	290	1	1.49	0.029
2013.1.00668.S	Weiss	6	260	1	1.50	0.029
2013.1.00815.S	Willott	6	263	1	0.97	0.017
2013.1.00884.S	Alexander	7	343	53	1.29	0.188
2013.1.01258.S	Riechers	7	296	7	1.05	0.062
2013.1.01292.S	Leiton	7	344	45	1.05	0.253
2015.1.00137.S	Scoville	6	240	68	1.09	0.081
"		7	343	300	0.57	0.158
2015.1.00695.S	Freundlich	6	261	6	0.05	0.030

Download table as:  ASCIITypeset image

Table 2.  Summary of Relations for ISM, SFR, and Accretion

	Equation #
${M}_{\mathrm{ISM}}=7.07\,\times {10}^{9}\,\,{M}_{\odot }\times {(1+z)}^{1.84}\times {(\mathrm{sSFR}/{\mathrm{sSFR}}_{\mathrm{MS}})}^{0.32}\,\times$ ${\left(\tfrac{{M}_{* }}{{10}^{10}{M}_{\odot }}\right)}^{0.30}$	(6)
$\mathrm{SFR}=0.31\,\,{M}_{\odot }\,{\mathrm{yr}}^{-1}\times \left(\tfrac{{M}_{\mathrm{ISM}}}{{10}^{9}\,{M}_{\odot }}\right)\times {(1+z)}^{1.04}\,\times$ ${(\mathrm{sSFR}/{\mathrm{sSFR}}_{\mathrm{MS}})}^{0.70}\times {\left(\tfrac{{M}_{* }}{{10}^{10}{M}_{\odot }}\right)}^{0.01}$	(7)
${\dot{M}}_{\mathrm{acc}}=2.27\,\,{M}_{\odot }\,{\mathrm{yr}}^{-1}\times {(1+z)}^{3.60}\,\times$ $\left({\left(\tfrac{{M}_{* }}{{10}^{10}{M}_{\odot }}\right)}^{0.56}-0.56\times {\left(\tfrac{{M}_{* }}{{10}^{10}{M}_{\odot }}\right)}^{0.74}\right)$	(13)
${\tau }_{\mathrm{dep}}=3.23\,\mathrm{Gyr}\,\times {(1+z)}^{-1.04}\times {(\mathrm{sSFR}/{\mathrm{sSFR}}_{\mathrm{MS}})}^{-0.70}\times {M}_{10}^{-0.01}\,$	(8)
${M}_{\mathrm{ISM}}/{M}_{* }=0.71\,\times {(1+z)}^{1.84}\times {(\mathrm{sSFR}/{\mathrm{sSFR}}_{\mathrm{MS}})}^{0.32}\times {M}_{10}^{-0.70}\,$	(9)
${f}_{\mathrm{gas}}\equiv \tfrac{{M}_{\mathrm{ISM}}}{{M}_{* }+{M}_{\mathrm{ISM}}}=\{1\,+\,1.41\,\times$ ${(1+z)}^{-1.84}\times {(\mathrm{sSFR}/{\mathrm{sSFR}}_{\mathrm{MS}})}^{-0.32}\times {M}_{10}^{0.70}\}{}^{-1}$	(10)

Download table as:  ASCIITypeset image

3.2. IR Source Catalog

3.3. Redshifts, Stellar Masses, and SFRs

3.4. Complete Sample of ALMA-detected IR-bright Galaxies

3.5. Source Measurements

At each IR source position falling within the ALMA primary beam HPBW (typically 20'' in Band 7), we searched for a significant emission source ($\gt 2\sigma$) within 2'' radius of the IR source position. This radius is the expected maximum size for these galaxies. The adopted detection limit implies that ∼2% of the detections at the $2\sigma$ limit could be spurious. Since there are ∼240 galaxies detected at 2–3σ, we can expect that about four of the detections could be false.²⁰

Some of the sources are likely to be somewhat extended relative to the synthesized beams (typically ∼1''); we therefore measure both the peak and integrated fluxes. The latter were corrected for the fraction of the synthesized beam falling outside the aperture. The adopted final flux for each source was the maximum of these as long as the signal-to-noise ratio (S/N) was $\gt 2$.

The noise for both the integrated and peak flux measures was estimated by placing 50 randomly positioned apertures of similar size in other areas of the FOV and measuring the dispersion of those measurements. The synthesized beams for most of these observations were ∼1'', and the interferometry should have good flux recovery out to sizes ∼4 times this; since the galaxy sizes are typically $\leqslant 2$ to 3'', we expect the flux recovery to be nearly complete, that is, there should be relatively little resolved-out emission. All measured fluxes were corrected for primary beam attenuation. The maximum correction is a factor of 2 when the source is near the HPBW radius for the 12 m telescopes.

A total of 708 of the Herschel sources were found within the ALMA FOVs, and the positions of this sample were used as priors for the ALMA flux measurements. For the sample of 708 objects, the measurements yielded 708 and 182 objects with $\gt 2$ and $\gt 7\sigma$ detections, respectively. Thus at 2σ, all sources were detected. No correction for Malmquist bias was applied since there were detections for the complete sample of sources falling within the survey area.

We then restricted the sample used for fitting to z and M_* where the sampling is good (see Section 6). A total of 772 flux measurements were made—some of the Herschel sources had multiple ALMA observations. In some cases the duplicate pointings were in both Bands 6 and 7 (see Appendix B).

In summary, all of the Herschel sources within the ALMA pointings were detected. The final sample of detected objects with their sSFRs and redshifts is shown in Figure 3; the histogram distributions of ${M}_{* }$ and SFR are in Figure 4. Figure 5 shows the distribution of measured fluxes and derived ISM masses.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In approximately 10% of the images, there is more than one detection. However, the redshifts of these secondary sources and the distributions of their offsets from the primary source indicate that most of the secondaries are not physically associated with the primary sources.

Here we briefly summarize the physical basis underlying the estimation of ISM masses from the Rayleigh–Jeans dust emission. Appendix A provides more detail, including the full equation used for mass calculations; a complete discussion is given in Scoville et al. (2016).

4.1. Physical Basis

The far-infrared dust continuum emission of galaxies arises from interstellar dust, heated by absorption of short-wavelength radiation (UV, optical, and near-infrared). The short-wavelength sources of this luminosity are recently formed stars or AGNs; hence the far-infrared luminosity provides a measure of that portion of the SF and nuclear activities that is obscured by dust. For the standard Galactic ISM at solar metallicity, a modest column of gas and dust, ${N}_{H}\gtrsim 2\times {10}^{21}$ cm⁻², will produce an extinction ${A}_{V}=1$ mag. Since the dust grains are small ($\lesssim 1$ micron), the dust opacity is highest in the UV and optical but decreases strongly at far-infrared and submillimeter wavelengths (${\kappa }_{\nu }\propto {\lambda }^{\sim -1.8}$ at $\lambda \gtrsim 200\,\mu {\rm{m}}$; Planck Collaboration 2011b).

At long wavelengths on the Rayleigh–Jeans tail, the dust emission is almost always optically thin, and the emission flux per unit mass of dust is only linearly dependent on the dust temperature. Thus the flux observed on the RJ tail provides a linear estimate of the dust mass and hence the ISM mass, provided the dust emissivity per unit mass and the dust-to-gas abundance ratio can be constrained. Fortunately, both of these prerequisites are well established from observations of nearby galaxies (e.g., Draine et al. 2007; Galametz et al. 2011).

Theoretical understanding of the dust emissivity has also significantly improved in the last two decades (Draine & Li 2007). On the optically thin RJ tail of the IR emission, the observed flux density is given by

$$\begin{eqnarray*}&&{S}_{\nu }\propto {\kappa }_{D}(\nu ){T}_{D}{\nu }^{2}\displaystyle \frac{{M}_{D}}{{d}_{L}^{2}},\end{eqnarray*}$$

where ${T}_{D}$ is the temperature of the emitting dust grains, ${\kappa }_{D}(\nu )$ is the dust opacity per unit mass of dust, ${M}_{D}$ is the total mass of dust, and d_L is the source luminosity distance. Thus, the mass of dust and ISM can be estimated from the observed specific luminosity ${L}_{\nu }$ on the RJ tail:

$$\begin{eqnarray}&&{M}_{D}\propto \displaystyle \frac{{L}_{\nu }}{\langle {T}_{D}{\rangle }_{M}{\kappa }_{D}(\nu )}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{M}_{\mathrm{ISM}}\propto \displaystyle \frac{{L}_{\nu }}{{f}_{D}\langle {T}_{D}{\rangle }_{M}{\kappa }_{D}(\nu )}.\end{eqnarray} \tag{ 4 }$$

Here, $\langle {T}_{D}{\rangle }_{M}$ is the mean mass-weighted dust temperature, and f_D is the dust-to-ISM mass ratio (typically $\sim 1/100$ for solar-metallicity ISM).

4.2. Empirical Calibration

In Scoville et al. (2016), we used a sample of 30 local star-forming spiral galaxies, 12 ultraluminous IR galaxies (ULIRGs), and 30 z $\sim 2\mbox{--}3$ submillimeter galaxies (SMGs) to empirically calibrate the RJ luminosity-to-mass ratio:

$$\begin{eqnarray}{\alpha }_{\nu } & \equiv & \langle {L}_{{\nu }_{850\mu {\rm{m}}}}/{M}_{\mathrm{mol}}\rangle \\ & = & 6.7\pm 1.7\times {10}^{19}\mathrm{erg}\,{{\rm{s}}}^{-1}{\mathrm{Hz}}^{-1}{{M}_{\odot }}^{-1}\,.\end{eqnarray} \tag{ 5 }$$

The gas masses were in all cases estimated from global measurements of the CO (1–0) emission in the galaxies assuming a single Galactic CO conversion factor ${\alpha }_{\mathrm{CO}}=6.5\,\,{M}_{\odot }/{\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1}{\mathrm{pc}}^{2}$ or ${X}_{\mathrm{CO}}=3\times {10}^{20}$ N(H₂) cm⁻² (K km ${{\rm{s}}}^{-1}{)}^{-1}$. These molecular gas masses include a factor of 1.36 to account for the associated mass of heavy elements (mostly He at 8% by number).

It is noteworthy that across the three diverse samples of galaxies the empirical calibration varied less than a factor of 2 (see Figure 17). The calibration therefore includes only massive galaxies (Milky Way or greater), which are likely to have near-solar metallicity. The samples of high-redshift galaxies in the present paper are therefore restricted to stellar mass ${M}_{* }\ \gt 2\times {10}^{10}\,{M}_{\odot }$ (and most are above 5 × 10¹⁰ $\,{M}_{\odot }$).

Probing lower-mass galaxies, which presumably would have significantly subsolar metallicity, will require careful calibration as a function of metallicity or mass. We note that in Draine et al. (2007) there is little evidence of variation in the dust-to-gas abundance ratio for the first factor of 4–5 down from solar metallicity. However, at lower metallicities, the dust-to-gas abundance does clearly decrease (see Figure 17 in Draine et al. 2007 and Figure 16 in Berta et al. 2016).

Lastly, we point out that since the standard Galactic ISM dust-to-gas abundance ratio (quantified by $2\times {10}^{21}\,{\rm{H}}\,{\mathrm{cm}}^{-2}/{{\rm{A}}}_{{\rm{V}}}$) is generally the same in both molecular and atomic hydrogen regions, one might expect that a similar calibration ${\alpha }_{\nu }$ pertains to the atomic and molecular phases of galaxies, provided one restricts the sample to approximately solar-metallicity HI and H₂ regions in the galaxies, such as the inner disks. For this reason, we will henceforth refer to our mass estimates as M_ISM. The ISM mass estimates found below for the high-z galaxies are very large (5–100 times that of the Galaxy); one could expect that the dominant phase is in fact molecular within the high-z galaxies.

4.3. Mass-weighted ${T}_{D}$

The measured ALMA fluxes were converted to estimates of the ISM masses using Equation (14). The distributions of flux and derived M_ISM are shown in Figure 5. The estimated ${M}_{\mathrm{ISM}}$ are shown along with the SFRs in Figure 6 for the 708 detected sources. Here one clearly sees a strong increase in SFRs for galaxies with large ${M}_{\mathrm{ISM}};$ however, there is also a very large dispersion, implying that the SFR per unit ISM mass depends on other variables. For example, as clearly seen in the color coding (redshift) of the points, there must be a strong dependence on z. Yet, this cannot be the entire explanation since the high-redshift galaxies with a large range of ISM masses exhibit similarly high SFRs (that is, a lot of the galaxies in each redshift range are distributed horizontally over a range of ${M}_{\mathrm{ISM}}$). For a given SFR, there are big spreads in sSFR relative to the MS and in stellar mass. The other variables linking SF and gas contents are the stellar masses and the elevation above the MS. We fit for all of these dependencies simultaneously in the following sections.

Zoom In Zoom Out Reset image size

In order to unravel the dependencies of ISM content and SFR on intrinsic and extrinsic galaxy parameters, we have fit a power-law dependence of ${M}_{\mathrm{ISM}}$ on the most likely parameters: z, ${M}_{* }$, and MS ratio ($=\mathrm{sSFR}/{\mathrm{sSFR}}_{\mathrm{MS}}$); then we use this power-law expression for ${M}_{\mathrm{ISM}}$ in order to elucidate the dependencies of the SFRs/M${}_{\mathrm{ISM}}$ on z, ${M}_{* }$, and the MS ratio. In this second stage of fitting, the terms can be viewed as SFR efficiencies per unit mass of ISM gas, as they depend separately on z, ${M}_{* }$, and the MS ratio.

To avoid areas of z and stellar mass where the sampling is low, we restrict our sample further to z = 0.3 to 3 and ${M}_{* }\gt 3\times {10}^{10}\,{M}_{\odot }$, thus reducing the areas where a Malmquist bias could be significant. The upper redshift limit is to avoid galaxies for which the flux measurement would be off the RJ tail; the lower mass limit is to avoid galaxies with expected low metallicity. The final sample used for fitting the dependencies has 575 distinct galaxies.

In all of the empirical fittings, we adopt a power law in each of the independent variables,

$$\begin{eqnarray*}&&{M}_{\mathrm{ISM}}\,\mathrm{and}\,\mathrm{SFR}=\mathrm{constant}\,P{1}^{\alpha }\times P{2}^{\beta }\times P{3}^{\gamma },\end{eqnarray*}$$

and solve for the minimum chi-square fit as a function of the independent variables (P). A Monte Carlo Markov chain (MCMC) routine (MLINMIX_ERR; in IDL) was used for the fitting. This is a Bayesian method for linear regression that takes into account measurement errors in all variables, as well as their intrinsic scatter. The Markov chain Monte Carlo algorithm (with Gibbs sampling) is used to randomly sample the posterior distribution (Kelly 2007).

6.1.  ${M}_{\mathrm{ISM}}$

The result of the MCMC fitting for the dependence of the ${M}_{\mathrm{ISM}}$ on redshift, MS ratio, and stellar mass is

$$\begin{eqnarray}{M}_{\mathrm{ISM}} & = & 7.07\pm 0.88\,\times {10}^{9}\,\,{{\rm{M}}}_{\odot }\times {(1+{\rm{z}})}^{1.84\pm 0.14}\\ & & \times {(\mathrm{sSFR}/{\mathrm{sSFR}}_{\mathrm{MS}})}^{0.32\pm 0.06}\times {\left(\displaystyle \frac{{M}_{* }}{{10}^{10}{M}_{\odot }}\right)}^{0.30\pm 0.04}.\end{eqnarray} \tag{ 6 }$$

The uncertainties were derived from the MCMC fitting; they do not account for possible calibration uncertainties such as the conversion from L_IR to SFR. In Appendix C, we provide plots showing the covariance of the fitted parameters with their distributions (see Figure 20 left). The posterior distributions of the different parameters in the MCMC fitting are single-valued (i.e., nondegenerate) and smooth. This also indicates that there are no degenerate parameters and hence that the choice of variables is appropriate.

In Figure 7 the observed ISM masses (solid dots) are shown in the left panel. Their fractional differences between the fit and the observations are shown in the middle and right panels. The Figure 7 middle and right panels show a scatter roughly equivalent to the observed values; however, given the large dynamic range (a factor of ∼100) in the ${M}_{\mathrm{ISM}}$ and the fitting parameters, this scatter is within the combined uncertainties of those parameters.

Zoom In Zoom Out Reset image size

Thus Equation (6) quantifies three major results regarding the ISM contents and their variation for high-redshift galaxies relative to low-redshift galaxies:

• 1.  
  The ISM masses clearly increase toward higher z, depending on ${(1+z)}^{1.84}={({(1+z)}^{2.9})}^{0.63}$; that is, they are not evolving as rapidly as the SFRs that vary as ${(1+z)}^{2.9}$.
• 2.  
  Above the MS, the ISM content increases, but not as rapidly as the SFRs (0.32 versus unity power laws).
• 3.  
  The ISM contents increase as ${M}_{* }^{0.30}$, indicating that the gas mass fractions must decrease in the higher-${M}_{* }$ galaxies.

The first conclusion clearly implies that the SF efficiency per unit gas mass must increase at high redshift (as discussed below). The second conclusion indicates that the galaxies above the MS have higher gas contents, but not in proportion to their elevated SFRs, and the third conclusion indicates that higher stellar mass galaxies are relatively gas-poor. Thus, galaxies with higher stellar mass likely use up their fuel at earlier epochs and have lower specific accretion rates (see Section 8) than the low-mass galaxies. This is a new aspect of the "downsizing" in the cosmic evolution of galaxies.

6.2. SFR

In fitting for the SFR dependencies, we wish to clearly distinguish between the obvious intuition that when there is more ISM there will be both more SF and a higher efficiency for converting the gas to stars. Thus, we impose a linear dependence of the SFR on ${M}_{\mathrm{ISM}}$, using ${M}_{\mathrm{ISM}}$ taken from Equation (6) rather than going back to the observed ${M}_{\mathrm{ISM}}$ values. Effectively, we are then fitting for the star formation efficiencies ($\mathrm{SFR}/{M}_{\mathrm{ISM}}$) for star formation per unit gas mass as a function of z, MS ratio, and ${M}_{* }$. The use of ISM masses from Equation (6) is necessary in order to isolate the efficiency variation with redshift, MS ratio, and ${M}_{* }$ from the variation of the ISM masses with the same three parameters. The result is shown in Figure 8 with color codings by redshift and MS ratio. The plots indicate that the dispersion in the fit is distributed around unity, and there appear to be no systematic offsets with respect to redshift or MS ratio.

Zoom In Zoom Out Reset image size

The MCMC solution is as follows:

$$\begin{eqnarray}\mathrm{SFR} & = & 0.31\pm 0.01\,\,{M}_{\odot }\,{\mathrm{yr}}^{-1}\times \left(\displaystyle \frac{{M}_{\mathrm{ISM}}}{{10}^{9}\,{M}_{\odot }}\right)\\ & & \times {(1+z)}^{1.05\pm 0.05}\\ & & \times {(\mathrm{sSFR}/{\mathrm{sSFR}}_{\mathrm{MS}})}^{0.70\pm 0.02}\\ & & \times {\left(\displaystyle \frac{{M}_{* }}{{10}^{10}{M}_{\odot }}\right)}^{0.01\pm 0.01}.\end{eqnarray} \tag{ 7 }$$

The uncertainties in each of the fit parameters were derived from the MCMC fitting (see Appendix C). The covariances for the SFR fit are shown in Figure 20 right. The posterior distributions of the different parameters in the MCMC fitting are well-behaved.

Figure 8 shows the fitted ISM masses (Equation (6)) and measured SFRs. In the left panel, the color coding is according to the redshift, and in the right panel according to their sSFR relative to the MS. Here one can see that the dominant dependence is on the ISM masses, but there is clearly a dependence on redshift and sSFR.

Equation (7) was written in a form such that the three power-law terms can be interpreted as efficiencies (SFR per unit gas mass). The resulting solution (Equation (7)) indicates the following:

• 1.  
  The SFR per unit ISM mass must increase approximately as ${(1+z)}^{1}$ as compared to the ${(1+z)}^{2.9}$ dependence of the SFRs; that is, the efficiency of converting gas to stars (the SFR per unit mass of gas) is clearly increasing at high z. (The increase in ISM content is separately represented by the linear term ${M}_{\mathrm{ISM}}$.)
• 2.  
  The SFR efficiency similarly increases as the 0.7 power of sSFR/sSFR_MS, implying that the enhanced SFRs above the MS are partially due to higher SF efficiencies (SFR per unit gas mass) in those galaxies—equivalently, the SB galaxies have a shorter gas-depletion timescale.
• 3.  
  There is no significant dependence of the SFR efficiency on stellar mass (${M}_{* }^{0.01}$).

The lack of dependence of the SF efficiency on galaxy mass is reasonable. If at high redshift the SF gas is in self-gravitating GMCs (as at low z), the internal structure of the GMCs then influences the physics of the SF, and the gas does not know that it is in a more or a less massive galaxy.

Figure 9 shows the relative evolutionary dependencies of the SFRs, the ISM masses, and the SF efficiency per unit mass of ISM gas normalized to unity at z = 0. The fundamental conclusion is that the elevated rates of SF activity at both high redshift and above the MS are due to both increased gas contents and increased efficiencies for converting the gas to stars.

Zoom In Zoom Out Reset image size

Using the empirically based relations for ${M}_{\mathrm{ISM}}$ and the SFR per unit ISM mass (Equations (6) and (7)), we now explore the consequences of these relations as obtained by simple algebraic manipulation.

7.1. Gas-depletion Timescales and Gas Mass Fraction

Combining Equations (6) and (7), one can derive the gas-depletion time (${\tau }_{\mathrm{dep}}={M}_{\mathrm{ISM}}$/SFR) and the gas mass fraction (f_gas = ${M}_{\mathrm{ISM}}$/(${M}_{* }$ + ${M}_{\mathrm{ISM}}$):

$$\begin{eqnarray}{\tau }_{\mathrm{dep}} & \equiv & {M}_{\mathrm{ISM}}/\mathrm{SFR}\\ & = & 3.23\pm 0.10\,\mathrm{Gyr}\,\times {(1+z)}^{-1.05\pm 0.05}\\ & & \times {(\mathrm{sSFR}/{\mathrm{sSFR}}_{\mathrm{MS}})}^{-0.70\pm 0.02}\times {M}_{10}^{-0.01\pm 0.01}\,,\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}{M}_{\mathrm{ISM}}/{M}_{* } & = & 0.71\pm 0.09\,\times {(1+z)}^{1.84\pm 0.14}\\ & & \times {(\mathrm{sSFR}/{\mathrm{sSFR}}_{\mathrm{MS}})}^{0.32\pm 0.06}\times {M}_{10}^{-0.70\pm 0.04}\,\end{eqnarray} \tag{ 9 }$$

and

$$\begin{eqnarray}{f}_{\mathrm{gas}} & \equiv & \displaystyle \frac{{M}_{\mathrm{ISM}}}{{M}_{* }+{M}_{\mathrm{ISM}}}\\ & = & \{1\,+\,1.41\pm 0.18\times {(1+z)}^{-1.84\pm 0.14}\\ & & {\times {(\mathrm{sSFR}/{\mathrm{sSFR}}_{\mathrm{MS}})}^{-0.32\pm 0.06}\times {M}_{10}^{0.70\pm 0.04}\}}^{-1},\end{eqnarray} \tag{ 10 }$$

where ${M}_{10}={M}_{* }/{10}^{10}\,{M}_{\odot }$. Note that Equation (8) is obtained simply by canceling out the leading ${M}_{\mathrm{ISM}}$ term in Equation (7) (since the latter already has a linear term in ${M}_{\mathrm{ISM}}$)—not by dividing the full Equation (6) by (7).

The ${\tau }_{\mathrm{dep}}$ and f_gas are shown in Figure 10. (In keeping with convention, these depletion times do not include a correction for the mass return to the ISM during stellar evolution.) In Table 3, we have normalized these equations to z = 2 and ${M}_{* }=5\times {10}^{10}\,{M}_{\odot }$ (instead of z = 0 used earlier).

Zoom In Zoom Out Reset image size

Table 3.  Relations Normalized to z = 2 and ${M}_{* }=5\times {10}^{10}\,{M}_{\odot }$

	Equations #
$\mathrm{Normalized}\,\mathrm{to}\,z\,=\,2\,\mathrm{and}\,{M}_{* }=5\times {10}^{10}\,{M}_{\odot }$
${M}_{\mathrm{ISM}}=8.65\times {10}^{10}\,\,{M}_{\odot }\,\times [{(1+z{)}_{2}^{1.84}\times (\mathrm{sSFR}/{\mathrm{sSFR}}_{\mathrm{MS}})}^{0.32}\times {M}_{* \,5}^{0.30}]$	(6)
$\mathrm{SFR}=9.9\times \left(\tfrac{{M}_{\mathrm{ISM}}}{{10}^{10}\,{M}_{\odot }}\right)\,\,{M}_{\odot }\,{\mathrm{yr}}^{-1}\,\times$ $[{(1+z{)}_{2}^{1.04}\times (\mathrm{sSFR}/{\mathrm{sSFR}}_{\mathrm{MS}})}^{0.70}\,\times {M}_{* \,5}^{0.01}]$	(7)
$\,\,\,\,\,\,\,\,=\,85\,\,{M}_{\odot }\,{\mathrm{yr}}^{-1}\,\,\times \left[\left(\tfrac{{M}_{\mathrm{ISM}}}{8.65\times {10}^{10}\,{M}_{\odot }}\right)\,\times \right.$ ${(1+z{)}_{2}^{1.04}\times (\mathrm{sSFR}/{\mathrm{sSFR}}_{\mathrm{MS}})}^{0.70}\times {M}_{* \,5}^{0.01}]$	(7)
${\dot{M}}_{\mathrm{acc}}=73\,\,{M}_{\odot }\,{\mathrm{yr}}^{-1}\times [2.3\times (1+z{)}_{2}^{3.60}\times ({M}_{* \,5}^{0.56}-0.56\times {M}_{* \,5}^{0.74})]$	(13)
${\tau }_{\mathrm{dep}}\equiv \tfrac{{M}_{\mathrm{ISM}}}{\mathrm{SFR}}=1.01\,\mathrm{Gyr}\,\times [{(1+{\rm{z}}{)}_{2}^{-1.04}\times (\mathrm{sSFR}/{\mathrm{sSFR}}_{\mathrm{MS}})}^{-0.70}\times {{\rm{M}}}_{* \,5}^{-0.01}]$	(8)
$\mathrm{gas}/\mathrm{stellar}\equiv \tfrac{{M}_{\mathrm{ISM}}}{{M}_{* }}=1.74\,\times [{(1+z{)}_{2}^{1.84}\times (\mathrm{sSFR}/{\mathrm{sSFR}}_{\mathrm{MS}})}^{0.32}\times {M}_{* \,5}^{-0.70}]$	(9)
${f}_{\mathrm{gas}}\equiv \tfrac{{M}_{\mathrm{ISM}}}{{M}_{* }+{M}_{\mathrm{ISM}}}=0.63\,\times$ $[1.58/(1\,+\,0.58\times {(1+z{)}_{2}^{-1.84}\times (\mathrm{sSFR}/{\mathrm{sSFR}}_{\mathrm{MS}})}^{-0.32}\times {M}_{* \,5}^{0.70})]$	(10)

Note. The equations are written in a form such that the quantity in $[\,]$ in each equation is equal to unity at z = 2 and ${M}_{* }=5\times {10}^{10}\,{M}_{\odot }$; ${(1+z)}_{2}$ is $(1+z)$ normalized to its value at z = 2 where (1 +z ) = 3; ${M}_{* 5}$ is the stellar mass normalized to 5 × 10¹⁰$\,{M}_{\odot }$. As noted in Section 7, the fourth relation is obtained by canceling out the ${M}_{\mathrm{ISM}}$ term in the second equation, not by division of the first equation by the second. See the original equations in the text for the uncertainties in the coefficients.

Download table as:  ASCIITypeset image

Table 4.  Sources with 240 and 345 GHz Measurements

R.A.	Decl.	z	$\langle {\nu }_{240}\rangle$	${{\rm{S}}}_{\nu }$	$\langle {\nu }_{345}\rangle$	S ${{\rm{S}}}_{\nu }$	Ratio
			GHz	mJy	GHz	mJy
149.64888	2.59813	2.94	240.0	1.12	343.5	3.43	3.07 ± 0.6
149.65558	2.71626	2.42	240.0	3.95	343.5	9.58	2.43 ± 0.3
149.66782	2.08743	0.97	240.0	1.35	341.9	5.95	4.42 ± 1.2
149.72572	2.27947	3.03	240.0	4.24	343.5	14.98	3.53 ± 0.2
149.75085	1.85219	1.85	240.0	3.50	343.5	12.48	3.57 ± 0.3
149.75095	1.85469	3.21	240.0	3.12	323.3	7.42	2.38 ± 0.2
149.87177	2.21219	2.38	240.0	2.23	343.5	20.08	9.02 ± 1.7
149.92032	2.02038	2.38	240.0	4.50	343.5	12.88	2.86 ± 0.2
150.03714	2.66954	2.93	240.0	1.37	343.5	5.73	4.19 ± 0.4
150.08232	2.53454	2.62	240.0	3.20	343.5	8.71	2.72 ± 0.3
150.09865	2.36537	2.58	240.0	0.62	343.5	3.70	6.01 ± 1.9
150.10616	2.05351	1.23	240.0	2.60	343.5	6.32	2.43 ± 0.3
150.14705	2.73147	2.21	240.0	1.94	343.5	6.16	3.17 ± 0.6
150.15022	2.47518	1.94	240.0	1.37	343.5	1.84	1.34 ± 0.2
150.17995	2.08864	1.97	240.0	0.99	343.5	1.50	1.51 ± 0.3
150.31131	2.58844	3.01	240.0	2.20	344.8	6.87	3.12 ± 0.5
150.31342	2.71619	1.80	245.0	0.54	344.8	1.57	2.91 ± 0.7
150.34572	2.33485	2.92	240.0	2.41	344.8	3.93	1.63 ± 0.2
150.34946	1.93700	2.32	240.0	2.58	343.5	7.58	2.93 ± 0.2

Note. The expected flux ratio for 25 K dust temperature and opacity spectral index of 1.8 is 2.52 ± 0.01 across the range of redshifts sampled above.

Download table as:  ASCIITypeset image

Depletion times at $z\gt 1$ are $\sim (2\mbox{--}10)\times {10}^{8}\,\mathrm{years}$, much shorter than for low-z galaxies (e.g., $\sim 2\times {10}^{9}\,\mathrm{years}$ for the Milky Way). The depletion times have no dependence on ${M}_{* }$.

The gas mass fractions shown in Figure 10 right are very high ($\gt 50$%) at $z\gt 2$ for SF galaxies, systematically decreasing at later epochs. They also show a strong dependence on stellar mass ($\propto {M}_{* }^{-0.7}$, Equation (9)) and whether a galaxy is on or above the MS (see Figure 10 right). It is therefore clear that there is no single gas mass fraction varying solely as a function of redshift (that is, independent of the sSFR relative to the MS and M_*).

For MS galaxies, Saintonge et al. (2013) and Genzel et al. (2015) obtain slightly lower depletion times of ∼450 Myr and gas mass fractions of ∼40% at z ∼ 2.8. Above the MS, they find increased ISM masses and SF efficiencies, in agreement with the results here. Their results are consistent with the estimates shown in Figure 10. But it is not really appropriate to quote single value estimates for the gas mass fraction at each redshift, neglecting the sSFR and M_* dependencies.

Using the MS continuity principle (Section 2.3) and the ISM contents obtained from Equation (6) with $\mathrm{sSFR}/{\mathrm{sSFR}}_{\mathrm{MS}}=1$, we now derive the net accretion rates of MS galaxies required to maintain the MS evolutionary tracks. Along each evolutionary track (curved lines in Figure 2), the rate balance must be given by

$$\begin{eqnarray}&&\displaystyle \frac{{{dM}}_{\mathrm{ISM}}}{{dt}}={\dot{M}}_{\mathrm{acc}}-(1-{f}_{\mathrm{mass}\mathrm{return}})\times \mathrm{SFR},\end{eqnarray} \tag{ 11 }$$

assuming that major merging events are rare. Here, ${f}_{\mathrm{mass}\mathrm{return}}$ is the fraction of stellar mass returned to the ISM through stellar mass loss, taken to be 0.3 here (Leitner & Kravtsov 2011). Since these paths are following the galaxies in a Lagrangian fashion, the time derivatives of a mass component M must be taken along the evolutionary track and

$$\begin{eqnarray}\displaystyle \frac{{dM}}{{dt}} & = & \displaystyle \frac{{dM}}{{dz}}\displaystyle \frac{{dz}}{{dt}}+\displaystyle \frac{{dM}}{{{dM}}_{* }}\displaystyle \frac{{{dM}}_{* }}{{dt}}\\ \displaystyle \frac{{dM}}{{dt}} & = & \displaystyle \frac{{dM}}{{dz}}\displaystyle \frac{{dz}}{{dt}}+\mathrm{SFR}\displaystyle \frac{{dM}}{{{dM}}_{* }}.\end{eqnarray} \tag{ 12 }$$

Figure 11 shows the accretion rates using Equation (11). The data used for this figure were generated from Equations (6) and (7) and thus are consistent with those earlier equations. These rates are ∼100 $\,{M}_{\odot }$ yr⁻¹ at z ∼ 2.5.

Zoom In Zoom Out Reset image size

The power-law fit to the accretion rates is

$$\begin{eqnarray}{\dot{M}}_{\mathrm{acc}} & = & 2.27\pm 0.24\,\,{M}_{\odot }\,{\mathrm{yr}}^{-1}\times {(1+z)}^{3.60\pm 0.26}\\ & & \times ({M}_{10}^{0.56\pm 0.04}-0.56\pm 0.04\times {M}_{10}^{0.74\pm 0.05})\end{eqnarray} \tag{ 13 }$$

where ${M}_{10}={M}_{* }/{10}^{10}\,{M}_{\odot }$. The combination of two separate mass terms is required to match the curvatures shown in Figure 11; they do not have an obvious physical justification. The first term dominates at low mass and the second at higher masses.

The value of the accretion rate at z = 2 and ${M}_{* }=5\times {10}^{10}\,{M}_{\odot }$ is 73 $\,{M}_{\odot }$ yr⁻¹ with a ${(1+z)}^{\sim 3.59}$ dependence (see Table 3). Theoretical estimates for the halo accretion rate of baryons yield lower values at z ∼2 of 33 $\,{M}_{\odot }$ yr⁻¹ with a (1 + z)${}^{2.5}$ dependence (Dekel et al. 2013), which is a somewhat slower redshift evolution.

Two important points to emphasize are that (1) these accretion rates should be viewed as net rates (that is, the accretion from the halos minus any outflow rate from SF or AGN feedback) and (2) these rates refer only to the MS galaxies where the evolutionary continuity is a valid assumption.

The derived accretion rates are required in order to maintain the SF in the early universe galaxies. Even though the existing gas contents are enormous compared to present-day galaxies, the observed SFRs will deplete this gas within $\sim 5\times {10}^{8}\,\mathrm{years};$ this is short compared to the MS evolutionary timescales. The large accretion rates, comparable to the SFRs, suggest that the higher SF efficiencies deduced for the high-z galaxies and for those above the MS may be dynamically driven by the infalling gas and mergers. These processes will shock compress the galaxy disk gas and possibly be the cause of the higher SF efficiencies.

It is worth noting that although one might think that the accretion rates could have been readily obtained simply from the evolution of the MS SFRs, this is not the case. One needs the mapping of ${M}_{\mathrm{ISM}}$ and its change with time in order to estimate the first term on the right-hand side of Equation (11). Figure 12 shows the relative evolution of each of the major rate functions over cosmic time and stellar mass. Comparing the proximity of the curves as a function of redshift, one can see modest differential change in the accretion rates and SFRs as a function of redshift and stellar mass.

Zoom In Zoom Out Reset image size

One might ask what are reasonable accretion rates to adopt for the SB galaxies above the MS, since they are not necessarily obeying the continuity assumption? Here there appear two reasonable possibilities: either the accretion rate is similar to the MS galaxy at the same stellar mass, or, if the elevation above the MS was a consequence of galactic merging, one might assume a rate equal to twice that of an individual galaxy with half the stellar mass. The latter assumption would imply a $\sqrt{2}$ higher accretion rate, thus being consistent with the higher ISM masses of the galaxies above the MS. In this case, the higher SFRs will be maintained longer than the simple depletion time it takes to reduce the preexisting ISM mass back down to that of an MS galaxy. The same $\sqrt{2}$ factor of increase in the ISM mass will arise from the merging of the preexisting ISM masses of two galaxies of half the observed mass. This follows from the dependence of ${M}_{\mathrm{ISM}}$ on stellar mass, varying only as ${M}_{* }^{0.30}$, rather than linearly. Thus, for two reasons, the notion that the SB-region galaxies are the result of major merging becomes quite attractive.

In Figure 13, we show the evolution of the SB galaxies with initial ${M}_{* }={10}^{10}\,{M}_{\odot }$ and initial SFR a factor of 6 above the MS. For reference, the evolution of MS galaxies with the same stellar mass is also shown. For the SB galaxy, the initial ISM mass is ${6}^{0.32}=1.77$ higher than the MS galaxy (see Equation (6)). The accretion rate is taken as equal to that of an MS galaxy with the same ${M}_{* }$. This calculation shows that ultimately the SB ends up with approximately a factor of 2 greater stellar mass—due to the larger initial ISM mass and the fact that the SB evolves more rapidly to higher stellar mass and thus accretes at a greater rate at high z.

Zoom In Zoom Out Reset image size

Using the mass functions (MF) of SF and passive galaxies (Ilbert et al. 2013), we estimate the total cosmic mass density of the ISM as a function of redshift using Equation (6). (This is the equivalent of the Lilly–Madau plot for the SFR density as a function of redshift.) We do this for the redshift range z = 0 to 4 and ${M}_{* }={10}^{10}$ to 10¹² $\,{M}_{\odot }$, a modest extrapolation of the ranges covered in the data presented here. Figure 14 left shows the derived cosmic mass densities of stars (SF and passive galaxies) and ISM as a function of redshift. We applied Equation (6) only to the SF galaxies and did not include any contribution from the passive galaxy population; to include the SB population, we multiplied the ISM mass of the normal SF population by a factor of 1.1. If the galaxy distribution were integrated down to a stellar mass equal to 10⁹ $\,{M}_{\odot }$, the stellar masses (and presumably the ISM masses) are increased by 10 to 20% (Ilbert et al. 2013).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The evolution of the ISM mass density shown in Figure 14 left is similar in magnitude to the theoretical predictions based on semianalytic models by Obreschkow et al. (2009), Lagos et al. (2011), and Sargent et al. (2013) (see Figure 12 in Carilli & Walter 2013). However, all of their estimations exhibit a more constant density at z > 1. The empirically based, prescriptive predictions of Popping et al. (2015) exhibit closer agreement with the evolution found by us; they predict a peak in the ISM gas at z ≃ 1.8 and a falloff at higher and lower redshift similar to that seen in Figure 14 left. Berta et al. (2013) estimated the evolution of ISM from the evolution of the SFRD by adopting a gas-depletion time from Tacconi et al. (2013). In Figure 14 right, their estimates are compared with ours and with the values from the simulation of Sargent et al. (2013).

The gas mass fractions computed for galaxies with ${M}_{* }={10}^{10}$ to 10¹² $\,{M}_{\odot }$ are shown in Figure 16. The ISM is dominant over the stellar mass down to $z\simeq 1.5$. At z = 3 to 4 the gas mass fractions get up to ∼80% when averaged over the galaxy population. Thus, the evolution of ISM contents that peak at z ≃ 2 is likely responsible for the peak in SF and AGN activity at that epoch (see Figure 15 in Madau & Dickinson 2014, and references therein). At z = 4 down to 2, the buildup in the ISM density is almost identical to that of the cosmic SFRD shown in Figure 14. (The ISM density point at z ∼ 0.3 is uncertain since it relies on extrapolation of Equation (6) to low ${M}_{* }$ and low z, where relatively few galaxies exist in our sample.)

A comparison of the MS evolutionary tracks shown in Figure 13 with those shown in Figure 2 provides a critical confirmation of the self-consistency of our derived relations for the ISM masses (Equation (6)), the SFR law (Equation (7)), and the accretion rates (Equation (13)). In Figure 2, the MS evolution tracks were derived based solely on integration of the MS stellar mass and SFR relations as a function of redshift. Thus, those tracks make no use of the relations derived here.

On the other hand, the MS evolutionary tracks in Figure 13 were obtained by integration in time of the galaxy evolution with the ISM mass, SFR, and accretion given as functions of z, sSFR, and ${M}_{* }$ by Equations (6), (7), and (13), respectively. Once the initial galaxy mass and z were specified, the initial ISM mass was taken from Equation (6). Then the SFR was calculated from Equation (7), based on the ISM mass, and the accretion rate was calculated at each time step from Equation (13). The evolution was integrated in time using a second-order Runge–Kutta method. The initial and final cosmic time for the integration was given by the redshifts for each pair of MS curves. Thus, the evolutionary tracks shown in Figure 13 rely entirely on the relations derived in this paper, specifically the dependence of ISM mass on MS location (determining the initial M_* for the integration) and, subsequently, the SFR as a function of ISM mass and the accretion rate as a function of stellar mass and redshift.

Finding that the MS evolutionary tracks in Figure 13 match those in Figure 2, one concludes that the relations derived here are entirely consistent with the a priori known MS dependencies and the evolution of the MS. Although the tracks in Figure 13 are shown in the figure only for initial ${M}_{* }={10}^{10}\,{M}_{\odot }$, we find a similar consistency for initial ${M}_{* }=5\times {10}^{10}\,{M}_{\odot }$, indicating that the dependence on ${M}_{* }$ in the equations is also self-consistent.

In summary, the evolutionary tracks shown in Figure 13 were derived using the empirical fits to the ISM mass, the SFR as a function of ISM mass, and the accretion rate function; the fact that the tracks for the MS galaxies follow closely those appearing in Figure 2 provides assurance the derived fits are consistent with the known MS evolution. (This does not provide a check on the validity of the relations above the MS, since we do not know a priori the end points of the SB evolutionary tracks.)

There are now over 200 detections of CO line emission at high redshift ($z\gt 2$), and we here compare our results with those studies. A couple caveats or cautions are required when comparing these results. (1) Most of the high-z CO detections are of CO (2–1) and (3–2); the inference of a molecular gas mass therefore requires an assumption for the scaling for the luminosity in these high-J CO lines relative to the mass-calibrated CO (1–0) line. (2) Many of those studies have adopted a Galactic ${\alpha }_{\mathrm{CO}(1--0)}=2\,$ $\times \,{10}^{20}\,{\mathrm{cm}}^{-2}$ (K $\,\mathrm{km}\,{{\rm{s}}}^{-1}$)⁻¹, based largely on Galactic gamma-ray observations; see the review by Bolatto et al. (2013). Here we have used ${\alpha }_{\mathrm{CO}(1-0)}=3\,$ $\times \,{10}^{20}\,{\mathrm{cm}}^{-2}$ (K $\,\mathrm{km}\,{{\rm{s}}}^{-1}$)⁻¹ in our calibration of the dust emission to gas masses (see Appendix A). This ${\alpha }_{\mathrm{CO}(1-0)}$ is derived from correlation of the CO line luminosities and virial masses for resolved Galactic giant molecular clouds. We believe the former value is not as reliable—it entails an assumption that the cosmic rays that produce the ∼2 MeV gamma rays by interaction with the gas fully penetrate the GMCs and, even more suspect, that their density is constant with Galactic radius (see Appendix A and the appendix in Scoville et al. 2016). Below, we note where the works have used a different ${\alpha }_{\mathrm{CO}}$ than that used in our calibration of the dust-based gas masses.

The most extensive study of high-z CO (3–2) emission is that of Tacconi et al. (2013), who detected 38 galaxies at z ∼ 1.2 and 14 at $z\sim 2.2$. (A somewhat more extensive sample including low-redshift galaxies is included in Genzel et al. (2015).) In their analysis, they make use of six CO (2–1) detections at z ∼ 1.5 (Daddi et al. 2010) and six at z ∼ 1 (Magdis et al. 2012a). They find a single best-fit depletion time of 0.7 ± 0.2 Gyr at z = 1 to 2.2, which would correspond to 1.1 ± 0.3 Gyr for the CO conversion factor used here. Allowing for uncertainties in the assumptions used here and in the CO study, this is consistent with the depletion time estimated here: $\simeq 1\,\mathrm{Gyr}$ at z = 2 (see Table 3).

These estimates can be compared with the z ∼ 0 estimate of 1.24 ± 0.06 Gyr (Saintonge et al. 2011), including HI. Genzel et al. (2015) find that the depletion times vary as ${(1+z)}^{-0.3}\ \times {(\mathrm{sSFR}/{\mathrm{sSFR}}_{\mathrm{MS}})}^{-0.5}\times {M}_{* }^{\sim 0}$, compared with ${(1+z)}^{-1.04}\ \times {(\mathrm{sSFR}/{\mathrm{sSFR}}_{\mathrm{MS}})}^{-0.70}\times {M}_{* }^{-0.01}$ from our work. Thus, we are finding a considerably steeper dependence on redshift, but similar dependencies on the elevation above the MS and the stellar mass (Table 3). They also see the gas contents varying linearly with the evolution of sSFR of the MS, whereas we find a 0.63 power-law dependence.

The molecular gas fractions (${f}_{\mathrm{gas}}={M}_{\mathrm{gas}}/({M}_{\mathrm{gas}}+{M}_{* })$) found by Tacconi et al. (2013) were 0.49 and 0.47 at z = 1.2 and 2.2, respectively. At z ∼ 2, we find ${f}_{\mathrm{gas}}\simeq 0.6$ (Figure 16) and 0.5 at z = 1 for MS galaxies. Tacconi et al. (2013) used a lower ${\alpha }_{\mathrm{CO}}$, so their gas fractions are necessarily lower than our estimates; however, we see stronger evolution with redshift. Both the CO and dust-based estimates show a decreasing gas mass fraction at higher ${M}_{* }$ (seen also by Magdis et al. 2012a). Genzel et al. (2015) find ${M}_{{\rm{H}}2}/{M}_{* }\propto {(1+z)}^{3}$, whereas we find ${M}_{{\rm{H}}2}/{M}_{* }\propto {(1+z)}^{1.84}\times \mathrm{sSFR}/{\mathrm{sSFR}}_{\mathrm{MS}}{)}^{0.32}\times {M}_{* }^{-0.70}$ (Table 3); we thus have a more gradual evolution with redshift in the gas mass fractions. Some of the difference could be understood if the sampling in sSFR and ${M}_{* }$ were very different between the two samples.

Zoom In Zoom Out Reset image size

Popping et al. (2015) predict a decreasing gas-to-stellar mass ratio for higher stellar mass galaxies similar to that in Equation (6). Their work is based on tracking the halo gas contents and then using empirical relations from low-z galaxies to model the galaxy properties such as size, SFR, and H₂ content. Figure 6 in Popping et al. (2015) indicates ${M}_{{\rm{H}}2}/{M}_{* }\simeq 0.3$ and 1 at z = 1 and 3, averaged over ${M}_{* }={10}^{10}\mbox{--}{10}^{12}\,{M}_{\odot }$; Figure 16 indicates corresponding values of 0.42 and 3.5 after converting from ${M}_{\mathrm{ISM}}/({M}_{* }+{M}_{\mathrm{ISM}})$. Thus, at z = 1 there is reasonable agreement, but at higher z we are finding larger gas mass fractions than Popping et al. (2015) did.

Magdis et al. (2012a) also examined the variations above the MS, obtaining ${M}_{\mathrm{gas}}/{M}_{* }\,=$ $\,2.05\pm 0.32\times {(\mathrm{sSFR}/{\mathrm{sSFR}}_{\mathrm{MS}})}^{0.87}$. This can be compared with our relation in Table 3, ${M}_{\mathrm{gas}}/{M}_{* }\,=$ $\,1.69\pm 0.1\times {(\mathrm{sSFR}/{\mathrm{sSFR}}_{\mathrm{MS}})}^{0.32}$ at z = 2. Their work also uses the dust continuum to estimate masses. However, they also fit the far-infrared SED to obtain a variable T_D, which we have argued against on physical grounds (see Appendix A). Schinnerer et al. (2016) observed a sample of 86 galaxies at 240 GHz with ALMA and detected 45 at z = 2.8 to 3.6. For the detected objects they obtain a median ${\tau }_{\mathrm{dep}}=0.68\,\mathrm{Gyr}$. Presumably this value would be larger if the nondetected sources were included. Their median gas-to-stellar mass ratio was 1.65 ± 0.17 for the 45 objects. (Their ALMA data set, which is now public, was included in the work presented here.)

In summary, there appears to be reasonable agreement between the results derived from CO line studies and those derived here, based on the RJ dust emission. This is indeed quite reassuring given the uncertainties in the higher-J CO line and dust continuum calibrations.

The very large galaxy sample presented here, based on the relatively short dust continuum measurements (∼2 minutes per galaxy with ALMA) and with uniform classification of the individual galaxy properties (redshift, stellar mass, and SFR relative to the MS), has allowed us to thoroughly explore the variations of gas content with these properties.

We have analyzed the ISM gas contents of a sample of 708 galaxies at $z\gt 0.3$, having stellar masses determined from optical/NIR SED fitting and well-constrained SFRs since all are detected in the far-infrared with Herschel. We quantify the evolution of the ISM contents as a function of redshift, ${M}_{* }$, and sSFR (relative to the MS) by fitting simple power-law expressions to the observed ${M}_{\mathrm{ISM}}$ for a sample of 575 galaxies at z = 0.3 to 3. The fit for the ISM contents is then combined with the observed SFRs to capture the changing efficiencies of SFR per unit gas mass with redshift, stellar mass, and sSFR relative to the MS. The redshift evolutionary dependent term in the power laws was taken to be that of the SFRs on the MS; hence, the power-law fits readily show the relative evolution of the ISM and the efficiencies for SF, compared to the well-known evolution of the MS SFRs. The derived equations are collected in Table 3 with normalization to z = 2 and ${M}_{* }=5\times {10}^{10}\,{M}_{\odot }$.

We find the following:

• 1.  
  The ISM contents of SF galaxies, both MS and SB, increase to high redshift less rapidly than the SFRs (0.63 power of the SFR evolution function, ${\mathrm{SFR}}_{\mathrm{MS}}\propto {(1+z)}^{2.9}$).
• 2.  
  Similarly, the ISM contents increase as the 0.36 power of the sSFR above the MS.
• 3.  
  The efficiency for forming stars per unit gas mass increases as the ∼0.32 power of the redshift evolution of the SFRs and the 0.7 power of the elevation above the MS.
• 4.  
  Combining points 1 to 3, it is clear that the increases in SF at high redshift and above the MS are due to both increased ISM masses and increased efficiency for converting gas to stars.
• 5.  
  The enhanced SF activity of galaxies above the MS, due to both their larger gas masses and higher SF efficiencies, suggests very plausibly that these starburst galaxies are the result of galaxy merging.
• 6.  
  The ISM contents increase as ${M}_{* }^{0.3}$, implying that the gas mass fractions decrease at higher stellar masses. The SF efficiency (SFR per unit gas mass) is virtually independent of ${M}_{* }$.
• 7.  
  We then estimate the accretion rates under the reasonable assumption of continuity from one epoch to the next in the MS galaxy populations. The derived net accretion rates (required to maintain the SF activity) are extreme, exceeding 50 $\,{M}_{\odot }$ yr⁻¹ above z = 2. The specific accretion rates (normalized to the stellar mass of the galaxy) decrease for higher mass galaxies.
• 8.  
  An analytic fit for the accretion rates shows that the accretion increases even more rapidly at high redshift than for the MS SFRs. Thus, it is the evolution of this accretion that drives the galaxy evolution in the early universe. In fact, the higher SF efficiencies at high redshift and above the MS may be linked to dynamical compression of the ISM by the infalling gas and minor mergers.
• 9.  
  To illustrate the power of these results, we can now chart the evolutionary paths of galaxies (both MS and SB galaxies) over cosmic time as they grow in mass and their SF dies out (due to decreased accretion; see Figure 13).

13.1. Comments and Implications

We thank Zara Scoville for proofreading the manuscript, and we thank the referee for a number of constructive suggestions. In addition, several useful references on the Galactic GMCs were provided by John Carpenter and John Bally, and good suggestions were provided by Fabian Walter. This paper makes use of the following ALMA data: ADS/JAO.ALMA 2011.0.00097.S, 2012.1.00076.S, 2012.1.00523.S, 2013.1.00034.S, 2013.1.00111.S, 2015.1.00137.S, 2013.1.00118.S, and 2013.1.00151.S. We plan to release the images used here via the IPAC/IRSA COSMOS archive at http://irsa.ipac.caltech.edu/data/COSMOS/ in 2017. ALMA is a partnership of ESO (representing its member states), NSF (USA), and NINS (Japan), together with NRC (Canada), NSC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. R.J.I. acknowledges support from the European Research Council (ERC) in the form of Advanced Grant 321302, COSMICISM. S.T. acknowledges support from the European Research Council in the form of Consolidator Grant 648179, ConTExt. B.D. acknowledges financial support from NASA through the Astrophysics Data Analysis Program (ADAP), grant number NNX12AE20G.

Here, we very briefly summarize the more thorough discussion in the appendix in Scoville et al. (2016), which establishes the foundation for using the long wavelength dust continuum as a tracer of ISM mass in high-redshift galaxies.

The empirical calibration of the technique (Scoville et al. 2016) is based on three different low- and high-redshift galaxy samples: (1) a sample of 30 local star-forming galaxies, (2) 12 low-z ultraluminous infrared galaxies (ULIRGs), and (3) 30 z ∼ 2 submillimeter galaxies (SMGs). These three samples with 72 galaxies are restricted to only those galaxies having good estimates of the total, source-integrated, long wavelength continuum and complete mapping of CO (1–0). (We avoid using higher-J CO lines since only the 1–0 transition has been well calibrated using large samples of Galactic GMCs with viral mass estimates; the higher CO lines have variable flux ratios with respect to the 1–0 line, so they are unlikely to be as reliable in mass estimations.) In calibrating the CO (1–0) masses, we have adopted ${\alpha }_{{CO}(1-0)}=3\,\times \,{10}^{20}\,{\mathrm{cm}}^{-2}{({\rm{K}}\mathrm{km}{{\rm{s}}}^{-1})}^{-1}$, which is derived from correlation of the CO line luminosities and virial masses for resolved Galactic GMCs. We believe this is more correct than the value obtained from Galactic gamma-ray surveys ($\alpha =2\times {10}^{20}$; see Bolatto et al. 2013) since the latter requires the questionable assumptions (1) that the cosmic rays which produce the ∼2 MeV gamma rays by interaction with the gas fully penetrate the GMCs and (2) that the cosmic-ray density is constant with Galactic radius. Obviously if one adopts the latter value, our derived scaling for the dust-based ISM masses must be reduced by a factor of 2/3. (See the appendix in Scoville et al. (2016) for a more extensive discussion.)

Figure 17 shows the ratio of specific luminosity at rest-frame $\lambda =850\,\mu {\rm{m}}$ to that of the CO (1–0) line, and one clearly sees a quite similar ratio of RJ dust continuum to CO luminosity. Using a standard Galactic CO (1–0) conversion factor, we then obtain the relation by which we convert the RJ dust continuum to ISM masses:

$$\begin{eqnarray}{M}_{\mathrm{ISM}} & = & 1.78\,{S}_{{\nu }_{\mathrm{obs}}}[\mathrm{mJy}]\,{(1+{\rm{z}})}^{-4.8}\,{\left(\displaystyle \frac{{\nu }_{850\mu {\rm{m}}}}{{\nu }_{\mathrm{obs}}}\right)}^{3.8}{({{\rm{d}}}_{L}[\mathrm{Gpc}])}^{2}\\ & & \times \,\left\{\displaystyle \frac{6.7\,\times \,{10}^{19}}{{\alpha }_{850}}\right\}\,\displaystyle \frac{{{\rm{\Gamma }}}_{0}}{{{\rm{\Gamma }}}_{\mathrm{RJ}}}\,\,{10}^{10}\,{M}_{\odot }.\end{eqnarray} \tag{ 14 }$$

Zoom In Zoom Out Reset image size

In Equation (14), Γ is a correction for departures from strict ${\nu }^{2}$ of the RJ continuum, and ${\alpha }_{850}=6.7\pm 1.7\ \times {10}^{19}\,\mathrm{erg}\,{{\rm{s}}}^{-1}{\mathrm{Hz}}^{-1}{{M}_{\odot }}^{-1}$ is the derived calibration constant between the 850 μm luminosity and ISM mass. We have adopted a dust opacity spectral index $\beta =1.8$, based on the determinations of the Planck observations in the Galaxy (Planck Collaboration 2011a, 2011b). (Berta et al. (2016) and Bianchi (2013) provide extensive discussions of possible variations in β.)

In the present work, the conversion of the fluxes measured with ALMA to ISM masses is done with Equation (14). Using Equation (14), the predicted fluxes for a fiducial ISM mass of 10¹⁰$\,{M}_{\odot }$ in the ALMA bands are shown in Figure 18 as a function of redshift.

Zoom In Zoom Out Reset image size

In the current work, we restricted the observed galaxies to be relatively massive (${M}_{* }\gt {10}^{10}$ $\,{M}_{\odot }$) since they should have close-to-solar metallicity and presumably not low dust-to-gas abundance ratios. We note that for the first factor of ∼5 down from solar metallicity in the galaxies analyzed by Draine et al. (2007), there is virtually no variation in the dust-to-gas abundance when one considers only objects with complete mapping in both CO (1–0) and the dust continuum (see Scoville et al. 2016). In fact, the dust abundance is likely to be more robust than CO, which can suffer depletion due to UV photodissociation as the metallicity drops.

The 19 objects that have both Band 6 and 7 measurements can be used to check for consistency with the spectral index of the submillimeter dust emissivity (${\kappa }_{\nu }=1.8$) and adopted dust temperature (T_D = 25 K) used for translating fluxes to masses. The S/N-weighted mean value of the Band 7/Band 6 flux ratio is 2.52 ± 0.03, consistent with the mean expected ratio of 2.52 (once one accounts for the departures from a strict Rayleigh–Jeans approximation). These measurements are shown in Figure 19. We have not explored constraining the dust temperature variations based on the range of the two-band ratios.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The covariance distributions obtained from the MCMC fitting for Equations (6) and (7) are shown. This is a Bayesian method for linear regression that takes into account measurement errors in all variables, as well as intrinsic scatter. A Markov chain Monte Carlo algorithm (Gibbs sampling) is used to randomly sample the posterior distribution.