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

TL;DR

This study uses ALMA dust continuum observations to analyze how interstellar medium mass, star formation, and gas accretion evolve in galaxies from redshift 0.3 to 4.5, revealing that early galaxy activity is driven by increased gas mass and efficiency.

Contribution

It provides new empirical relations between ISM mass, star formation, and redshift, highlighting the role of gas accretion in galaxy evolution at high redshift.

Findings

ISM masses increase with redshift and sSFR

Star formation efficiency rises with redshift and above the main sequence

Gas accretion rates are >100 solar masses per year at z > 2.5

Abstract

ALMA observations of the long wavelength dust continuum are used to estimate the interstellar medium (ISM) masses in a sample of 708 galaxies at z = 0.3 to 4.5 in the COSMOS field. The galaxy sample has known far-infrared luminosities and, hence, star formation rates (SFRs), and stellar masses (M$_{\rm *}$) from the optical-infrared spectrum fitting. The galaxies sample SFRs from the main sequence (MS) to 50 times above the MS. The derived ISM masses are used to determine the dependence of gas mass on redshift, M$_{\rm *}$, and specific SFR (sSFR) relative to the MS. The ISM masses increase approximately 0.63 power of the rate of increase in SFRs with redshift and the 0.32 power of the sSFR/sSFR$_MS$. The SF efficiencies also increase as the 0.36 power of the SFR redshift evolutionary and the 0.7 power of the elevation above the MS; thus the increased activities at early epochs are driven by both increased ISM masses and SF efficiency. Using the derived ISM mass function we estimate the accretion rates of gas required to maintain continuity of the MS evolution ($>100$ \msun yr$^{-1}$ at z $>$ 2.5). Simple power-law dependences are similarly derived for the gas accretion rates. We argue that the overall evolution of galaxies is driven by the rates of gas accretion. The cosmic evolution of total ISM mass is estimated and linked to the evolution of SF and AGN activity at early epochs.

Tables

Table 1. Table 1 ALMA Datasets used for Continuum Measurements

Project Code	PI	Band		$\#$ pointings	HPBW$>$	${\sigma }_{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

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

	Eq. #
${𝐌}_{\mathrm{𝐈𝐒𝐌}}=7.07\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟗}}{𝐌}_{\odot }\times {\left(\mathrm{𝟏}+𝐳\right)}^{1.84}\times {\left(\mathrm{𝐬𝐒𝐅𝐑}/{\mathrm{𝐬𝐒𝐅𝐑}}_{\mathrm{𝐌𝐒}}\right)}^{0.32}\times {\left(\frac{{𝐌}_{\ast }}{{\mathrm{𝟏𝟎}}^{\mathrm{𝟏𝟎}}{𝐌}_{\odot }}\right)}^{0.30}$	6
$\mathrm{𝐒𝐅𝐑}=0.31{𝐌}_{\odot }{\mathrm{𝐲𝐫}}^{-\mathrm{𝟏}}\times \left(\frac{{𝐌}_{\mathrm{𝐈𝐒𝐌}}}{{\mathrm{𝟏𝟎}}^{\mathrm{𝟗}}{𝐌}_{\odot }}\right)\times {\left(\mathrm{𝟏}+𝐳\right)}^{1.04}\times {\left(\mathrm{𝐬𝐒𝐅𝐑}/{\mathrm{𝐬𝐒𝐅𝐑}}_{\mathrm{𝐌𝐒}}\right)}^{0.70}\times {\left(\frac{{𝐌}_{\ast }}{{\mathrm{𝟏𝟎}}^{\mathrm{𝟏𝟎}}{𝐌}_{\odot }}\right)}^{0.01}$	7
${\dot{𝐌}}_{𝒂𝒄𝒄}=2.27{𝐌}_{\odot }{\mathrm{𝐲𝐫}}^{-\mathrm{𝟏}}\times {\left(\mathrm{𝟏}+𝐳\right)}^{3.60}\times \left({\left(\frac{{𝐌}_{\ast }}{{\mathrm{𝟏𝟎}}^{\mathrm{𝟏𝟎}}{𝐌}_{\odot }}\right)}^{0.56}-0.56\times {\left(\frac{{𝐌}_{\ast }}{{\mathrm{𝟏𝟎}}^{\mathrm{𝟏𝟎}}{𝐌}_{\odot }}\right)}^{0.74}\right)$	13
${𝝉}_{\mathrm{𝐝𝐞𝐩}}=3.23\mathrm{𝐆𝐲𝐫}\times {\left(\mathrm{𝟏}+𝐳\right)}^{-1.04}\times {\left(\mathrm{𝐬𝐒𝐅𝐑}/{\mathrm{𝐬𝐒𝐅𝐑}}_{\mathrm{𝐌𝐒}}\right)}^{-0.70}\times {𝐌}_{\mathrm{𝟏𝟎}}^{-0.01}$	8
${𝐌}_{\mathrm{𝐈𝐒𝐌}}/{𝐌}_{\ast }=0.71\times {\left(\mathrm{𝟏}+𝐳\right)}^{1.84}\times {\left(\mathrm{𝐬𝐒𝐅𝐑}/{\mathrm{𝐬𝐒𝐅𝐑}}_{\mathrm{𝐌𝐒}}\right)}^{0.32}\times {𝐌}_{\mathrm{𝟏𝟎}}^{-0.70}$	9
${𝐟}_{\mathrm{𝐠𝐚𝐬}}\equiv \frac{{𝐌}_{\mathrm{𝐈𝐒𝐌}}}{{𝐌}_{\ast }+{𝐌}_{\mathrm{𝐈𝐒𝐌}}}={\left\{\mathrm{𝟏}+1.41\times {\left(\mathrm{𝟏}+𝐳\right)}^{-1.84}\times {\left(\mathrm{𝐬𝐒𝐅𝐑}/{\mathrm{𝐬𝐒𝐅𝐑}}_{\mathrm{𝐌𝐒}}\right)}^{-0.32}\times {𝐌}_{\mathrm{𝟏𝟎}}^{0.70}\right\}}^{-\mathrm{𝟏}}$	10

Table 3. Table 3 Relations Normalized to z = 2 and 𝐌 ∗ = 𝟓 × 𝟏𝟎 𝟏𝟎 subscript 𝐌 5 superscript 10 10 \rm M_{*}=5\times 10^{10} 𝐌 ⊙ subscript 𝐌 direct-product \rm M_{\odot}

	Eq. #
$\mathrm{𝐍𝐨𝐫𝐦𝐚𝐥𝐢𝐳𝐞𝐝}\mathrm{𝐭𝐨}𝐳=\mathrm{𝟐}\mathrm{𝐚𝐧𝐝}{𝐌}_{\ast }=\mathrm{𝟓}\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟏𝟎}}$${𝐌}_{\odot }$
${𝐌}_{\mathrm{𝐈𝐒𝐌}}=8.65\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟏𝟎}}{𝐌}_{\odot }\times \left[{\left(\mathrm{𝟏}+𝐳\right)}_{\mathrm{𝟐}}^{1.84}\times {\left(\mathrm{𝐬𝐒𝐅𝐑}/{\mathrm{𝐬𝐒𝐅𝐑}}_{\mathrm{𝐌𝐒}}\right)}^{0.32}\times {𝐌}_{\ast \mathrm{𝟓}}^{0.30}\right]$	6
$\mathrm{𝐒𝐅𝐑}=9.9\times \left(\frac{{𝐌}_{\mathrm{𝐈𝐒𝐌}}}{{\mathrm{𝟏𝟎}}^{\mathrm{𝟏𝟎}}{𝐌}_{\odot }}\right){𝐌}_{\odot }{\mathrm{𝐲𝐫}}^{-\mathrm{𝟏}}\times \left[{\left(\mathrm{𝟏}+𝐳\right)}_{\mathrm{𝟐}}^{1.04}\times {\left(\mathrm{𝐬𝐒𝐅𝐑}/{\mathrm{𝐬𝐒𝐅𝐑}}_{\mathrm{𝐌𝐒}}\right)}^{0.70}\times {𝐌}_{\ast \mathrm{𝟓}}^{0.01}\right]$	7
$=\mathrm{𝟖𝟓}{𝐌}_{\odot }{\mathrm{𝐲𝐫}}^{-\mathrm{𝟏}}\times \left[\left(\frac{{𝐌}_{\mathrm{𝐈𝐒𝐌}}}{8.65\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟏𝟎}}{𝐌}_{\odot }}\right)\times {\left(\mathrm{𝟏}+𝐳\right)}_{\mathrm{𝟐}}^{1.04}\times {\left(\mathrm{𝐬𝐒𝐅𝐑}/{\mathrm{𝐬𝐒𝐅𝐑}}_{\mathrm{𝐌𝐒}}\right)}^{0.70}\times {𝐌}_{\ast \mathrm{𝟓}}^{0.01}\right]$	7
${\dot{𝐌}}_{𝒂𝒄𝒄}=\mathrm{𝟕𝟑}{𝐌}_{\odot }{\mathrm{𝐲𝐫}}^{-\mathrm{𝟏}}\times \left[2.3\times {\left(\mathrm{𝟏}+𝐳\right)}_{\mathrm{𝟐}}^{3.60}\times \left({𝐌}_{\ast \mathrm{𝟓}}^{0.56}-0.56\times {𝐌}_{\ast \mathrm{𝟓}}^{0.74}\right)\right]$	13
${𝝉}_{\mathrm{𝐝𝐞𝐩}}\equiv \frac{{𝐌}_{\mathrm{𝐈𝐒𝐌}}}{\mathrm{𝐒𝐅𝐑}}=1.01\mathrm{𝐆𝐲𝐫}\times \left[{\left(\mathrm{𝟏}+𝐳\right)}_{\mathrm{𝟐}}^{-1.04}\times {\left(\mathrm{𝐬𝐒𝐅𝐑}/{\mathrm{𝐬𝐒𝐅𝐑}}_{\mathrm{𝐌𝐒}}\right)}^{-0.70}\times {𝐌}_{\ast \mathrm{𝟓}}^{-0.01}\right]$	8
$\mathrm{𝐠𝐚𝐬}/\mathrm{𝐬𝐭𝐞𝐥𝐥𝐚𝐫}\equiv \frac{{𝐌}_{\mathrm{𝐈𝐒𝐌}}}{{𝐌}_{\ast }}=1.74\times \left[{\left(\mathrm{𝟏}+𝐳\right)}_{\mathrm{𝟐}}^{1.84}\times {\left(\mathrm{𝐬𝐒𝐅𝐑}/{\mathrm{𝐬𝐒𝐅𝐑}}_{\mathrm{𝐌𝐒}}\right)}^{0.32}\times {𝐌}_{\ast \mathrm{𝟓}}^{-0.70}\right]$	9
${𝐟}_{\mathrm{𝐠𝐚𝐬}}\equiv \frac{{𝐌}_{\mathrm{𝐈𝐒𝐌}}}{{𝐌}_{\ast }+{𝐌}_{\mathrm{𝐈𝐒𝐌}}}=0.63\times \left[1.58/\left(\mathrm{𝟏}+0.58\times {\left(\mathrm{𝟏}+𝐳\right)}_{\mathrm{𝟐}}^{-1.84}\times {\left(\mathrm{𝐬𝐒𝐅𝐑}/{\mathrm{𝐬𝐒𝐅𝐑}}_{\mathrm{𝐌𝐒}}\right)}^{-0.32}\times {𝐌}_{\ast \mathrm{𝟓}}^{0.70}\right)\right]$	10

Table 4. Table 4 Sources with 240 and 345 GHz Measurements

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

Equations

$$\mathrm{SF}{\mathrm{R}}_{\mathrm{MS}}$$

$$\displaystyle\rm\mbox{$\zeta_{\rm MS}$}(z)\equiv~{}SFR_{\rm MS}(z)/SFR_{\rm MS}(z=0)$$

$$\mathrm{SF}{\mathrm{R}}_{\mathrm{MS}}$$

$$\displaystyle\rm\mbox{$\zeta_{\rm MS}$}(M_{\rm*})$$

$$S_{\nu }$$

$$M_D$$

$$M_{ISM}$$

$${\alpha }_{\nu }$$

$${\mathrm{M}}_{\mathrm{ISM}}\mathrm{and}\mathrm{SFR}$$

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

$$\mathrm{SFR}$$

$${\tau }_{\mathrm{dep}}$$

$${\mathrm{M}}_{\mathrm{ISM}}/{\mathrm{M}}_{\ast }$$

$${\mathrm{f}}_{\mathrm{gas}}$$

$$\frac{\mathrm{d}{\mathrm{M}}_{\mathrm{ISM}}}{dt}={\dot{\mathrm{M}}}_{\mathrm{acc}}-(1-{\mathrm{f}}_{\mathrm{mass}\mathrm{return}})\times \mathrm{SFR},$$

$$\frac{\mathrm{dM}}{dt}$$

$${\dot{\mathrm{M}}}_{acc}$$

$$M_{\mathrm{ISM}}$$

Full text

1. Galaxy Evolution 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 (AGN) 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 AGN.

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 $\lesssim$ 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 $\sim$ 5% of the population but have 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 (SB), 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_{\rm*}$ 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; Schinnerer et al., 2016; Berta 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 two accuracy in the derived ISM masses – provided one restricts the galaxy samples to high stellar mass galaxies (>10^{10}$$\rm 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 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 2 (Section 4.3, see also (Scoville et al., 2016). We note also that in contrast to the CO, which typically constitutes just $\sim$0.1% of the gas mass, the dust is a 1% mass tracer. CO is also susceptible to depletion by photo-dissociation 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 non-contiguous but their fields of view (FOVs), totaling 102.9 arcmin2 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:

How does the gas content depend on the stellar mass of the galaxies? 2. 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. 3.

In the starburst populations, is the prodigious SF activity driven by increased gas supply or increased efficiency for converting the existing gas into stars? 4. 4.

How rapidly is the ISM being depleted? The depletion timescale is characterized by the ratio $M_{\rm 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. 5.

If the ISM depletion times are as short as they appear to be (several $\times 10^{8}$ yrs), yet the SF galaxy population persists for a much longer span of cosmic time ($\gtrsim 5\times 10^{9}$ yrs), 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. 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

We first provide a brief background to our investigation and develop the logical framework for tracking the galaxy evolution in Section 2. We then describe the datasets used for this investigation and the galaxy sample selection (Section 3). A brief background on the use of the RJ dust emission to probe ISM masses is given in Section 4; a more complete description is provided in Appendix A and Scoville et al. (2016).

In Section 6, analytic fits are obtained for the dependence of the ISM masses and SFRs on galaxy properties and redshift. These expressions are then used to elucidate the evolution of the ISM contents and the efficiencies for star formation as a function of redshift, stellar mass and sSFR relative to the MS. Section 7 uses these results to determine the gas depletion timescales and gas mass fractions. Section 8 provides estimates of the gas accretion rates required to maintain the SF in the galaxy population and Section 9 provides a brief discussion of the differences between the starburst and the MS galaxies. The derived equations are summarized in Tables 2 and 3.

The overall cosmic evolution of the ISM contents and the gas mass fractions integrated over the galaxy population from $M_{*}=10^{10}-10^{12}$ $\rm M_{\odot}$  is presented in Section 10. A check on the self-consistency of these relations is provided in Section 11 and a comparison with previous results based on CO and dust measurements is in given in Section 12. Section 13 provides a summary of the main results and their implications.

2. Evolution of the star formation rates and stellar mases

The bulk of the SF galaxy population can be located on a ‘Main Sequence’ (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).111See Eales et al. (2016) for some reservations regarding the MS.

2.1. Evolution of the MS with redshift

The MS locus (SFR(z,M${}_{\rm*})$) 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

[TABLE]

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

[TABLE]

For fitting 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 which have a $(1+z)^{2.9}$ dependence.

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${}_{\rm*}>5\times 10^{10}$ $\rm 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,

[TABLE]

and normalize to unity at log M∗ = 10.5 $\rm M_{\odot}$,

[TABLE]

Figure 1-Right shows the adopted shape function $\zeta_{\rm MS}$(M∗). The normalization constant is a SFR = 3.23 $\rm M_{\odot}$ $yr^{-1}$ at z = 0 and logM${}_{\rm*}=10.5$ $\rm 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${}_{\rm*})\sim$ 10.5 $\rm M_{\odot}$  which 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 Main Sequence 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 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% 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{}_{\rm*}>2\times 10^{11}$$\rm M_{\odot}  at z $>1$) and in dense environments at lower redshift (Peng et al., 2010; Darvish et al., 2016). At z $>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_{\rm*}/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 $\times 10^{10}$ $\rm 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:

the gas accretion-dominated and stellar mass buildup phase at z $>2$ corresponding to cosmic age less than 3.3 Gyr (see Section 8); 2. 2.

the transition phase where gas accretion approximately balances SF consumption and the evolution becomes diagonal and, 3. 3.

the epoch of ISM exhaustion at z $\lesssim$ 0.1 (age 12.5 Gyr) where the evolution will be vertically downward in the SFR versus M∗ plane.

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 3 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

The SFR evolution for the MS (relative to the MS at z = 0) is shown in Figure 1. Figure 2 shows the evolution for three fiducial stellar masses. In our analysis to fit the evolution of the ISM, SFR and accretion rates (Section 6), we will fit directly for 1+z dependence. The derived evolutionary trends with redshift can then be referenced to the well-known redshift evolution SFRs.

2.5. Starbursts

At each epoch there exists a much smaller population ($\sim 5$% by number at z = 2) which have SFRs 2 to 100 times that of the MS at the same stellar mass. Do these starburst galaxies quickly exhaust their supply of star forming gas, thus evolving rapidly back to the MS, or are their ISM masses systematically larger so that their depletion times differ little from the MS galaxies? These SB galaxies must be either a short-duration, but common evolutionary phase for the galaxies, or of long-duration, but a phase not undergone by the majority of the galaxy population. Despite their small numbers, their significance in the overall cosmic evolution of SF is greater than 5%, since they have 2 to 100 times higher SFRs.

3. Datasets

3.1. ALMA Band 6 and 7 Continuum Data

Within the COSMOS survey field, there now exist extensive observations from ALMA for the dust continuum of high redshift galaxies. Here, we make use of all those data which 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 datasets is 1011, covering a total area 0.0286 deg2 or 102.9 arcmin2 within the Half-Power Beam Width (HPBW).

3.2. IR Source Catalog

Our source finding used a positional prior: the Herschel-based catalog of far-infrared sources in the COSMOS field (Lee et al., 2013, 2015, 13597 objects from[). COSMOS was observed at 100 $\mu$m and 160 $\mu$m by Herschel PACS (Poglitsch et al., 2010) as part of the PACS Evolutionary Probe program (PEP; Lutz et al., 2011)), and down to the confusion limit at 250 $\mu$m, 350 $\mu$m, and 500 $\mu$m by Herschel SPIRE (Griffin et al., 2010) as part of the Herschel Multi-tiered Extragalactic Survey (HerMES; Oliver et al., 2012)).

In order to measure accurate flux densities of sources in the confusion-dominated SPIRE mosaics, it is necessary to extract fluxes using prior-based methods, as described in Lee et al. (2013). In short, we begin with a prior catalog that contains all COSMOS sources detected in the Spitzer 24$\mu$m and VLA 1.4 GHz catalogs (Le Floc’h et al., 2009; Schinnerer et al., 2010), which have excellent astrometry. Herschel fluxes are then measured using these positions as priors. The PACS 100 and 160 $\mu$m prior-based fluxes were provided as part of the PEP survey (Lutz et al., 2011), while the SPIRE 250, 350, and 500 $\mu$m fluxes were measured using the XID code of Roseboom et al. (2010, 2012), which uses a linear inversion technique of cross-identification to fit the flux density of all known sources simultaneously (Lee et al., 2013). From this overall catalog of infrared sources, we select reliable far-infrared bright sources by requiring at least 3$\sigma$ detections in at least 2 of the 5 Herschel bands. This should greatly limits the number of false positive sources in the catalog.

An in-depth analysis of the selection function for this particular catalog is provided in Lee et al. (2013), but in short, the primary selection function is set by the 24 $\mu$m and VLA priors catalog. As with many infrared-based catalogs, there is a bias toward bright, star-forming galaxies, but the requirement of detections in multiple far-infrared bands leads to a flatter dust temperature selection function than typically seen in single-band selections.

Since the selection function is biased to IR bright and massive galaxies, the sample is not representative of the high redshift SF galaxy population. However, in the analytic fitting below we obtain analytic dependencies for the ISM masses and SFRs on the sSFR, the stellar mass and redshift. These analytic fits can then be used to analyze the more representative populations. This approach is used in Section 10 to estimate the cosmic evolution of ISM.

3.3. Redshifts, Stellar Masses and Star Formation Rates

Spectroscopic redshifts were used for 5066 of these sources; the remainder had photometric redshifts from Ilbert et al. (2009) and Laigle et al. (2016) (This catalog does not include objects for which the photometric redshift fitting indicated a possible AGN.) For the final sample of sources falling within the ALMA pointings, 38% had spectroscopic redshifts.

The primary motivation for using the Herschel IR catalog for positional priors is the fact that once one has far-infrared detections of a galaxy, the SFRs can be estimated more reliably (including the dominant contributions of dust-obscured SF) rather than relying on optical/UV continuum estimation, which often have corrections by factors $\gg 5$ for dust obscuration. This said, the SFRs derived from the far-infrared are still probably individually uncertain by a factor 2, given uncertainties in the stellar IMF and the assumed timescale over which the young stars remain dust-embedded.

The conversion from IR (8-1000$\mu$m) luminosity makes use of SFR_{IR}=8.6\times 10^{-11}L_{IR}/\mbox{L_{\odot}}~{} using a Chabrier stellar IMF from 0.1 to 100 $\rm M_{\odot}$(Chabrier, 2003). The scale constant is equivalent to assuming that 100% of the stellar luminosity is absorbed by dust for the first $\sim$100 Myr and 0% for later ages. For a shorter dust enshrouded timescale of 10 Myr the scaling constant is $\sim$1.5 times larger (Scoville, 2013). In 706 of the 708 sources in our measured sample, the IR SFR was greater than the optical/UV SFR. The final SFRs are the sum of the opt/UV and the IR SFRs.

The stellar masses of the galaxies are taken from the photometric redshift catalogs (Ilbert et al., 2009; Laigle et al., 2016); these are also uncertain by at least factors of 2 due to uncertainties in the spectral energy distribution (SED) modeling and extinction corrections. Their uncertainties are probably less than those for the optically derived SFRs, since the stellar mass in galaxies is typically more extended than the SF activity, and therefore is likely to be less extincted.

The other galaxy property we wish to correlate with the derived ISM masses is the elevation of the galaxy above or below the SF MS. This enhanced SFR is quantified by sSFR/sSFR${}_{\rm MS}(z,M_{\rm*})$ with the MS definition taken from the combination of Speagle et al. (2014) and Lee et al. (2015) (see Section 2.1 and 2.2 and Figure 1).

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

The galaxy sample analyzed here includes 708 galaxies, yielding 575 to be used in the fitting after redshift and steller mass selection. This large number of objects can then overcome uncertainties in the z, SFR and M∗ for each individual object. This sample also has calibrations that are uniform across the full sample without the need for zero-point corrections. All the objects are in the COSMOS survey field and thus have the same photometry and scheme for evaluating the stellar masses and the redshifts; they also all have similar depth Herschel and Spitzer infrared observations.

For the ISM measurements, we use exclusively continuum observations from ALMA – these are consistently calibrated and with resolution ($\sim 1$′′) such that source confusion is not an issue. Lastly, our analysis involving the RJ dust continuum avoids the issue of variable excitation which causes uncertainty when using different CO transitions across galaxy samples. The excitation and brightness per unit mass for the different CO transitions is likely to vary by factors of 2 to 3 from one galaxy to another and within individual galaxies (see Carilli & Walter, 2013).

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 ($>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 $\sim 2$% of the detections at the $2\sigma$ limit could be spurious. Since there are $\sim 240$ galaxies detected at 2-3$\sigma$, we can expect $\sim 4$ of the detections could be false. 222We have also used detection thresholds of 3 and 4$\sigma$ and did not find significant changes in the fitting coefficients in the analysis below, that is less than 10% change; the 2$\sigma$ limit was therefore used since it yields measurements and uncertainty estimates for the complete sample of sources.

Some of the sources are likely to be somewhat extended relative to the synthesized beams (typically $\sim 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 SNR was $>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 $\sim 1$′′, and the interferometry should have good flux recovery out to sizes $\sim 4$ times this; since the galaxy sizes are typically $\leq 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 2 when the source is near the HPBW radius for the 12m telescopes.

A total of 708 of the Herschel sources were found within the ALMA FOVs, 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 $>2$ and $>7\sigma$ detections, respectively. Thus at 2$\sigma$, 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.

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.

4. ISM Masses Estimated from RJ Dust Emission

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 AGN; hence the far infrared luminosity provides a measure of that portion of the SF and nuclear activities which is dust obscured. For the standard Galactic ISM at solar metallicity, a modest column of gas and dust, $N_{H}\gtrsim 2\times 10^{21}$ cm*-2*, 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 submm wavelengths ($\kappa_{\nu}\propto\lambda^{\sim-1.8}$ at $\lambda\gtrsim 200\mu$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 of 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

[TABLE]

where $T_{\rm D}$ is the temperature of the emitting dust grains, $\kappa_{D}(\nu)$ is the dust opacity per unit mass of dust, $M_{\rm 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 observed specific luminosity $L_{\nu}$ on the RJ tail:

[TABLE]

Here $<T_{D}>_{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-3$ submm galaxies (SMGs) to empirically calibrate the RJ luminosity-to-mass ratio

[TABLE]

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_{CO}=6.5~{}\mbox{\rm M_{\odot}}/\rm{K~{}km~{}s^{-1}pc^{2}} or $X_{CO}=3\times 10^{20}$ N(H2) cm*-2* (K km s${}^{-1})^{-1}$. These molecular gas masses include a factor 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 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_{*}>2\times 10^{10}$ $\rm M_{\odot}$(and most are above $5\times 10^{10}$ $\rm M_{\odot}$).

Probing lower mass galaxies, which presumably would have significantly sub-solar 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~{}cm^{-2}/\rm A_{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 to approximately solar metallicity HI and H2 regions in the galaxies, such as the inner disks. For this reason, we will henceforth refer to our mass estimates as MISM. 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 $\rm T_{D}$

It is important to emphasize that the dust temperature entering Equation 3 is a mass-weighted dust temperature $\rm<T_{D}>_{M}$, not the luminosity-weighted temperature $\rm<T_{D}>_{L}$. The latter is what would be derived by fitting the overall SED of the emergent IR radiation to a black body or modified black body curve. In such spectral fitting, the derived TD is largely determined by the shape and location of the SED peak, and hence the dust emitting the bulk of the emergent luminosity.

Often, the far infrared SEDs are analyzed by fitting either modified black body curves or libraries of dust SEDs to the observed SEDs (e.g. Draine et al., 2007; da Cunha & Charlot, 2011; Magdis et al., 2012a, b). In essentially all instances the intrinsic SEDs used for fitting are taken to be optically thin. They thus do not include the attenuation expected on the short wavelength side of the far infrared peak, instead attributing the drop at short wavelengths solely to a lack of high temperature grains. The $\rm T_{D}$ determined in these cases is thus not even a luminosity-weighted $\rm T_{D}$ of all the dust, but just the dust above $\tau_{\nu}\sim 1$.

Even if the source were optically thin in the mid infrared, it should be clear that the temperature determined by SED fitting will be the luminosity-weighted temperature of grains undergoing strong radiative heating, i.e. those closest to the sources of luminosity – not a linear sampling of the mass of ISM in each cloud, the bulk of which is at lower temperature. To illustrate the difference between mass and luminosity weighted T’s, we refer to resolved studies of local GMCs: Orion (Lombardi et al., 2014), Auriga (Harvey et al., 2013) and W3 (Rivera-Ingraham et al., 2013).

4.3.1 Local GMCs and Galaxies

For the two Orion GMCs, the dominant infrared luminosity sources are the KL nebula/Orion A and Orion B. Lombardi et al. (2014) present a comprehensive analysis of Herschel and Planck data. The peak dust temperatures are $\rm T_{D}\sim 40-50$ K at IRAS and Herschel angular resolutions. Figures 2 and 7 in Lombardi et al. (2014) clearly show that these areas of peak dust temperature are confined to the immediate vicinity of the two OB star formation regions ($\sim 2$ arcmin in size) while the bulk of the GMCs, extending over 6 and 4∘ at $\lambda>200\mu$m, is at $\rm T_{D}=14-20$ K. A similar dichotomy between the luminosity- and mass-weighted TD is seen clearly in the Herschel data for the W3 GMC complexes (Rivera-Ingraham et al., 2013).

Harvey et al. (2013) provide a quantitative analysis of the relative masses of low and high TD dust based on their Herschel observations of the Auriga-California GMC . In this GMC, the major luminosity source is LkH$\alpha$101, which is the dominant source at 70$\mu$m, has $\rm T_{D}$ = 28K, and extends over $\sim 10$′. On the other hand the $\lambda=250\mu$m (SPIRE) image shows the GMC extending $\sim 6$∘ with a median dust temperature $\rm T_{D}=14.5$ K (Harvey et al., 2013). Only $\sim 2$% of the GMC mass is contained in the higher temperature region associated with LkH$\alpha$101.

Global SED fits for 11 Kingfish Survey galaxies using two temperature components were done by Galametz et al. (2011). The warm components were 51 - 59 K and the cold components had fitted T = 17 - 24 K and the fraction of the total far infrared luminosity in the warm component was 21 - 81%, with a mean value $\sim 50$% (see also Skibba et al., 2011; Dale et al., 2012).

Unfortunately, the extragalactic observations do not have the spatial resolution to resolve hot and cold dust regions in GMCs in other galaxies, and thus determine the relative masses of hot and cold ISM dust. One is thus forced to learn from the lessons of the above cited studies of Galactic GMCs or apply physical understanding.

With regard to the latter, the reason for the localization of the hot dust emitting most of the far infrared luminosity but containing little of the mass, is due to the strong wavelength dependence of the dust opacity. The short wavelength primary photons from the energy sources are absorbed in the first column of A${}_{V}\simeq$ 1 mag or $2\times 10^{21}$ H cm*-2*. On the other hand, the luminosity of the re-radiated secondary or tertiary photons will require much larger columns of gas and dust to be absorbed, due to their longer wavelength (since the dust absorption coefficient is much lower at the longer wavelength). Thus, it becomes obvious that the mass of the cold dust component will be 10-100 times larger than that of the warm dust. This must be the case unless the primary luminosity sources are uniformly distributed on the scale of A${}_{V}\sim 1$ within the clouds – which is clearly not the case in Galactic GMCs. One can’t rule out the possibility that the luminosity sources in the high z galaxies are more uniformly distributed on scales of A${}_{V}\sim 1$ mag; however, it would seem implausible that all the SF in high z galaxies occurs in such small columns of gas.

4.3.2 Adopted constant $\rm T_{D}=25$K

Here we have adopted a constant value for the mass-weighted temperature of the ISM dust to be used in our estimations of ISM masses. This value is slightly higher than the mean of 18 K determined by Planck observation of the cold dust in the Galaxy (Planck Collaboration, 2011b, a); we use the higher value in recognition of the fact that the mean radiation fields are likely to be somewhat higher in the high z galaxies. Our adoption of a constant value is different from the approach used by others to estimate ISM dust and gas masses (da Cunha & Charlot, 2011; Magdis et al., 2012a, b; Genzel et al., 2015; Schinnerer et al., 2016; Berta et al., 2016).

For the reasons outlined just above, we believe that using a $\rm T_{D}$ derived from the SED fits yields a luminosity-weighted TD, not the physically correct mass-weighted $\rm T_{D}$. That is not to say that there will not be some variations in the mass-weighted $\rm T_{D}$; however, with extragalactic observations which do not resolve the star forming regions, it is impossible to determine the mass-weighted TD. It is expected that variations in the mass-weighted $\rm T_{D}$ will be much less than the variations in the luminosity-weighted $\rm T_{D}$. The latter is dependent on the concentration of energy sources, which may vary considerably, whereas an increase in the overall mass-weighted $\rm T_{D}$ requires both a higher concentration of primary energy sources, and that they be uniformly distributed on scales comparable with the mean free path of the primary dust-heating photons.

Lastly, it is important to realize that the linear dependence on $\rm T_{D}$ for the flux to mass conversion and the expected small dynamic range (at most 15 - 40 K) in the mass-weighted $\rm T_{D}$ implies that the conversion factor will not have large errors with the adopted constant $\rm T_{D}=25$ K.

5. Observational Results

The measured ALMA fluxes were converted to estimates of the ISM masses using Equation A1. The distributions of flux and derived MISM are shown in Figure 5. The estimated MISM 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 MISM; 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 strong dependence on z. Yet, this can not 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 MISM). 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.

6. Dependence of MISM and SFR on z, M∗, and sSFR/sSFR${}_{MS}(\rm M_{\rm*},z)$

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 MISM on the most likely parameters: z, M∗, and MS-Ratio ($\rm=sSFR/sSFR_{\rm MS}$); then we use this power law expression for MISM in order to elucidate the dependences of the SFRs/MISM on z, $M_{\rm*}$, 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_{\rm*}$, 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${}_{\rm*}>3\times 10^{10}$ $\rm M_{\odot}$, thus reducing the areas where 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:

[TABLE]

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. MISM

The result of the MCMC fitting for the dependence of the MISM on redshift, MS-Ratio and stellar mass is:

[TABLE]

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., non-degenerate) 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 (filled dots) are shown in the Left panel. Their fractional differences between the fit and the observations are shown in Middle and Right panels. Figure 7-Middle and Right show a scatter roughly equivalent to the observed values; however, given the large dynamic range (a factor $\sim 100$) in the MISM and the fitting parameters, this scatter is within the combined uncertainties of those parameters.

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

The ISM masses clearly increase toward higher z, depending on $(1+z)^{1.84}=((1+z)^{2.9})^{0.63}$, that is, not evolving as rapidly as the SFRs which vary as $(1+z)^{2.9}$. 2. 2.

Above the MS, the ISM content increases, but not as rapidly as the SFRs (0.32 versus unity power laws). 3. 3.

The ISM contents increase as $M_{\rm*}^{0.30}$, indicating that the gas mass fractions must decrease in the higher $M_{\rm*}$ galaxies.

The first conclusion clearly implies 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 MISM, using MISM taken from Equation 6 rather than going back to the observed MISM values. Effectively, we are then fitting for the star formation efficiencies ($\rm SFR/M_{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.

The MCMC solution is:

[TABLE]

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 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 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 that:

The SFR per unit ISM mass must increase approximately as $\rm(1+z)^{1}$ as compared to the $\rm(1+z)^{2.9}$ dependence of the SFRs, i.e. the efficiency for converting gas to stars (the SFR per unit mass of gas) is clearly increasing at high-z. [The increase in ISM contents is separately represented by the linear term $\rm M_{\rm ISM}$.] 2. 2.

The SFR efficiency similarly increases as the 0.7 power of sSFR/sSFRMS, 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. 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.

7. Implications of the ISM and SFR Relations

Using the empirically based relations for MISM 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_{\rm dep}=\rm M_{\rm ISM}$ / SFR) and the gas mass fraction ($f_{gas}$= MISM/(M∗ + MISM):

[TABLE]

[TABLE]

and

[TABLE]

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

$\tau_{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${}_{\rm*}=5\times 10^{10}$ $\rm M_{\odot}$(instead of z = 0 used earlier).

Depletion times at z $>1$ are $\sim(2-10)\times 10^{8}$ yrs, much shorter than for low-z galaxies (e.g. $\sim 2\times 10^{9}$ yrs for the Milky Way). The depletion times have no dependence on $M_{\rm*}$.

The gas mass fractions shown in Figure 10-Right are very high ($>50$%) at z $>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 the MS and M∗).

For MS galaxies Saintonge et al. (2013) and Genzel et al. (2015) obtain slightly lower depletion times $\sim 450$ Myr and gas mass fractions $\sim 40$% at z $\sim 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.

8. Accretion Rates

Using the MS Continuity Principle (Section 2.3) and the ISM contents obtained from Equation 6 with $\rm sSFR/sSFR_{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:

[TABLE]

assuming that major merging events are rare. $\rm f_{mass~{}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

[TABLE]

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 $\sim 100$ $\rm M_{\odot}$yr*-1* at z $\sim$ 2.5.

The power-law fit to the accretion rates is

[TABLE]

where M_{10}=M_{*}/10^{10}\mbox{\rm M_{\odot}}. The combination of two separate mass terms is required to match the curvatures shown in Figure 11; they don’t 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${}_{\rm*}=5\times 10^{10}$ $\rm M_{\odot}$  is 73 $\rm M_{\odot}$ yr*-1* 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 $\sim 2$ of $33$ $\rm M_{\odot}$yr*-1* with a (1+z)2.5 dependence (Dekel et al., 2013), that is a somewhat slower redshift evolution.

Two important points to emphasize are: 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}$ yrs; 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 $\rm M_{\rm ISM}$ and its change with time in order to estimate the first term on the right 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.

9. SB versus MS Galaxy Evolution

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 $\sqrt{\rm 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 pre-existing ISM mass back down to that of a MS galaxy. The same $\sqrt{\rm 2}$ factor of increase in the ISM mass will arise from the merging of the pre-existing ISM masses of two galaxies of half the observed mass. This follows from the dependence of MISM on stellar mass, varying only as M${}_{\rm*}^{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${}_{\rm*}=10^{10}$ $\rm M_{\odot}$and initial SFR a factor 6 above the MS. For reference, the evolution of MS galaxies with the same stellar mass are 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 a MS galaxy with the same $\rm M_{\rm*}$. This calculation shows that ultimately the SB ends up with approximately a factor 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.

10. Cosmic Evolution of ISM and Stellar Mass

Using the mass functions (MF) of SF and passive galaxies (Ilbert et al., 2013), we estimate the total cosmic mass density of 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_{\rm*}=10^{10}$ to $10^{12}$ $\rm 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 stellar mass equal to $10^{9}$ $\rm M_{\odot}$, the stellar (and presumably the ISM masses) are increased by 10 to 20% (Ilbert et al., 2013).

The evolution of the ISM mass density shown in Figure 14-Left is similar in magnitude to the theoretical predictions based on semi-analytic models by Obreschkow et al. (2009); Lagos et al. (2011); 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 $\simeq$ 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_{\rm*}=10^{10}$ to $10^{12}$ $\rm 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 $\sim 80$% when averaged over the galaxy population. Thus, the evolution of ISM contents which peak at z $\simeq$ 2, is likely responsible for the peak in SF and AGN activity at that epoch (see Figure 15 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 $\sim 0.3$ is uncertain since it relies on extrapolation of Equation 6 to low M∗ and low z, where there exist relatively few galaxies in our sample.)

11. A Self-Consistency  Test

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 a functions of z, sSFR and M∗ by Equations 6, 7, and 13. 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 2nd 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${}_{\rm*}=10^{10}$ $\rm M_{\odot}$, we find a similar consistency for initial M${}_{\rm*}=5\times 10^{10}$ $\rm 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.]

12. Comparison with Previous Work

There are now over 200 detections of CO line emission at high redshift (z $>$ 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, and 2) many of those studies have adopted a Galactic $\alpha_{\rm CO(1-0)}=2\times 10^{20}\rm cm^{-2}(K~{}\hbox{km~{}s}^{-1})^{-1}$, based largely on Galactic gamma ray observations, see review by Bolatto et al. (2013). Here we have used $\alpha_{\rm CO(1-0)}=3\times 10^{20}\rm cm^{-2}(\rm K~{}\hbox{km~{}s}^{-1})^{-1}$ in our calibration of the dust emission to gas masses (see Appendix A). This $\alpha_{\rm CO(1-0)}$ is derived from correlation of the CO line luminosities and virial masses for resolved Galactic Giant Molecular Clouds (GMCs). We believe the former value is not as reliable – it entails an assumption that the cosmic rays which produce the $\sim 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 Appendix in Scoville et al. (2016)). Below, we note where the works have used a different $\alpha_{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 $\sim 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 6 CO (2-1) detections at z $\sim 1.5$ (Daddi et al., 2010) and 6 at z $\sim 1$ (Magdis et al., 2012a). They find a single best-fit depletion time of 0.7$\pm 0.2$ Gyr at z =1 to 2.2. which would correspond to 1.1$\pm 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 Gyr at z = 2 (see Table 3).

These estimates can be compared with the z $\sim 0$ estimate of 1.24$\pm 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(sSFR/sSFR_{MS})^{-0.5}\times M_{*}^{\sim 0}$, compared with $(1+z)^{-1.04}\times(sSFR/sSFR_{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${}_{\rm gas}=\rm M_{\rm gas}/(M_{\rm gas}+M_{\rm*})$) found by Tacconi et al. (2013) were 0.49 and 0.47 at z = 1.2 and 2.2, respectively. At z $\sim 2$, we find f${}_{\rm gas}\simeq 0.6$ (Figure 16) and 0.5 at z =1 for MS galaxies. Tacconi et al. (2013) used a lower $\alpha_{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${}_{H2}/M_{*}\propto(1+z)^{3}$ whereas we find M$\rm{}_{H2}/M_{*}\propto\rm(1+z)^{1.84}\times sSFR/sSFR_{\rm 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 $\rm M_{*}$ were very different between the two samples.

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 H2 content. Figure 6 in Popping et al. (2015) indicates M$\rm{}_{H2}/M_{*}\simeq 0.3$ and 1 at z = 1 and 3, averaged over $M_{*}=10^{10}-10^{12}$ $\rm M_{\odot}$; Figure 16 indicates corresponding values of 0.42 and 3.5 after converting from M$\rm{}_{ISM}/(M_{*}+M_{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).

Magdis et al. (2012a) also examined the variations above the MS, obtaining M${}_{\rm gas}/\rm M_{\rm*}=2.05\pm 0.32\times(\rm sSFR/sSFR_{\rm MS})^{0.87}$. This can be compared with our relation in Table 3, M${}_{\rm gas}/\rm M_{\rm*}=1.69\pm 0.1\times(\rm sSFR/sSFR_{\rm 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_{dep}=0.68$ Gyr. Presumably this value would be larger if the non-detected sources were included. Their median gas-to-stellar mass ratio was $1.65\pm 0.17$ for the 45 objects. [Their ALMA dataset 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 ($\sim 2$ min per galaxy with ALMA) and with uniform classification of the individual galaxy properties (redshift, stellar mass and SFR relative to the MS), have allowed us to thoroughly explore the variations of gas content with these properties.

13. Summary

We have analyzed the ISM gas contents of a sample of 708 galaxies at z $>0.3$, having stellar masses determined from optical/NIR SED fitting and SFRs well-constrained 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 MISM 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${}_{\rm*}=5\times 10^{10}$ $\rm M_{\odot}$.

We find:

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, ($SFR_{MS}\propto(1+z)^{2.9}$). 2. 2.

Similarly, the ISM contents increase as the 0.36 power of the sSFR above the MS. 3. 3.

The efficiency for forming stars per unit gas mass increases as the $\sim$0.32 power of redshift evolution of the SFRs and the 0.7 power of the elevation above the MS. . 4. 4.

Combining # 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. 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. 6.

The ISM contents increase as $\rm M_{\rm*}^{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. 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$\rm M_{\odot}$yr*-1* above z =2. The specific accretion rates (normalized to the stellar mass of the galaxy) decrease for higher mass galaxies. 8. 8.

An analytic fit for the accretion rates shows that the accretion increases even more rapidly at high redshift than the MS SFRs. Thus, it is the evolution of this accretion which 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. 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

The variations of ISM masses, accretion and their relation to star formation have been explored with the most extensive sample yet of high redshift galaxies. Although the deduced estimates are consistent with existing studies using the CO lines (Tacconi et al., 2010; Daddi et al., 2010; Genzel et al., 2010; Riechers et al., 2011; Ivison et al., 2011; Magdis et al., 2012b; Saintonge et al., 2013; Carilli & Walter, 2013; Tacconi et al., 2013; Bolatto et al., 2015; Genzel et al., 2015), the sample of galaxies used here has the virtue that it maps the parameter space of z, $M_{\rm*}$, and sSFR quite effectively out to z = 3 using high quality and uniform ancillary data from the COSMOS survey field. We thus can simultaneously constrain the functional dependencies on redshift, sSFR relative to the MS and stellar mass. Our technique also does not suffer from the uncertainties introduced by variable excitation in the higher-J CO lines.

There are obvious extensions which need to be done in this field: extending with larger samples at z $<1$ and z $>3$. The former is straightforward since the lower z galaxies can be observed on the RJ tail in ALMA band 7 (for which the ALMA sensitivity is excellent and the fluxes are high); the latter will be more time-consuming since the observations must shift to Band 6 or even lower frequencies to stay on the Rayleigh-Jeans tail and the fluxes will be lower. In addition, the number of z $>$ 3 galaxies with measured far-infrared luminosities is much less.

Obtaining high quality CO (1-0) data on a subset of the galaxies analyzed here is also a high priority in order to firm up the calibration of dust flux to ISM masses at high redshift, and to determine the range of stellar masses above which it is valid. This will be both time-consuming and fraught with difficulties in the analysis since the CO line has its own calibration issues. But the spectral line data will also provide very helpful dynamical mass estimates – these can certainly be a useful reality check.

A major uncertainty for both the dust and the CO line studies is, of course, their dependence on metallicity (Z). Both are probably more robust than is generally assumed in the community. The CO line is heavily saturated (even in Galactic GMCs which have typical $\rm\tau_{CO(1-0)}\gtrsim 10$). To underscore this point, we note that the 13CO emission line is typically $\sim 1/5$ of the CO line flux in Galactic GMCs despite the much lower 13C/C abundance ($\sim$1/60 to 1/90). Thus, the line luminosity must scale as the $\sim$1/3 power of the CO abundance, and hence the metallicity.

With respect to the dust emission as a probe of ISM, it is reassuring that the dust-to-gas abundance ratio in low redshift galaxies appears fairly constant at $\sim$1% by mass from solar down to 1/5 solar metallicity (see Draine & Li, 2007) and (Berta et al., 2016, Figure 16) (although why this is the case is not understood).

Our finding that the ISM-to-stellar mass ratio and the accretion rates are both generally higher for lower mass galaxies has implications for the gas-phase metallicites of galaxies. Assuming that the metallicity of freshly accreting gas is significantly lower than that of the internal gas or the stars formed out of the gas, one would expect the gas phase metallicity to increase in higher stellar mass galaxies. This is, of course, known to be true; and it is a major motivation for our restriction to galaxies with relatively high M∗. It is also clear that a so-called ‘closed box’ model for the evolution of metal content has little physical justification in light of the extremely large accretion rates derived here.

We thank Zara Scoville for proof reading 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. RJI acknowledges support from the European Research Council (ERC) in the form of Advanced Grant, 321302, COSMICISM. ST acknowledges support from the European Research Council (ERC) 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.

Appendix A Long Wavelength Dust Continuum as an ISM Mass Tracer

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 3 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 $\sim$ 2 submm 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}cm^{-2}(K\hbox{km~{}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 $\sim 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$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:

[TABLE]

In Equation A1, $\Gamma$ is a correction for departures from strict $\nu^{2}$ of the RJ continuum, and \alpha_{850}=6.7\pm 1.7\times 10^{19}\rm erg~{}sec^{-1}Hz^{-1}{\mbox{\rm M_{\odot}}}^{-1} is the derived calibration constant between 850$\mu$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, 2011b, a). [Berta et al. (2016) and Bianchi (2013) provide extensive discussions of possible variations in $\beta$.]

In the present work, the conversion of the fluxes measured with ALMA to ISM masses is done with Equation A1. Using Equation A1, the predicted fluxes for a fiducial ISM mass of 10^{10}$$\rm M_{\odot} in the ALMA bands are shown in Figure 18 as a function of redshift.

In the current work we restricted the observed galaxies to be relatively massive ($M_{*}>10^{10}$ $\rm 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 $\sim 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 photo-dissociation as the metallicity drops.

Appendix B Dust Emission Spectral Index and Dust Temperature

The 19 objects which have both Band 6 and 7 measurements can be used to check for consistency with the spectral index of the submm dust emissivity ($\kappa_{\nu}$ = 1.8) and adopted dust temperature ($T_{D}=25$ K) used for translating fluxes to masses. The SNR-weighted mean value of the Band 7 / Band 6 flux ratio is $2.52\pm 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.

Appendix C Fitting Co-variances

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.