# The Evolution of the Baryons Associated with Galaxies Averaged over Cosmic Time and Space

Fabian Walter, Chris Carilli, Marcel Neeleman, Roberto Decarli, Gergö Popping, Rachel S Somerville, Manuel Aravena, Frank Bertoldi, Leindert Boogaard, Pierre Cox, et al.

► To cite this version:

Fabian Walter, Chris Carilli, Marcel Neeleman, Roberto Decarli, Gergö Popping, et al. The Evolution of the Baryons Associated with Galaxies Averaged over Cosmic Time and Space. The Astrophysical Journal, 2020, 902, 10.3847/1538-4357/abb82e. hal-02996651

HAL Id: hal-02996651

https://hal.science/hal-02996651

Submitted on 9 Nov 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# The Evolution of the Baryons Associated with Galaxies Averaged over Cosmic Time and Space

Fabian Walter¹,² ID, Chris Carilli² ID, Marcel Neeleman¹ ID, Roberto Decarli³ ID, Gergö Popping⁴ ID, Rachel S. Somerville⁵,⁶,
Manuel Aravena⁷ ID, Frank Bertoldi⁸ ID, Leindert Boogaard⁹ ID, Pierre Cox¹⁰ ID, Elisabete da Cunha¹¹, Benjamin Magnelli⁸ ID,
Danail Obreschkow¹¹, Dominik Riechers¹² ID, Hans-Walter Rix¹ ID, Ian Smail¹³ ID, Axel Weiss¹⁴ ID, Roberto J. Assef⁷ ID,
Franz Bauer¹⁵,¹⁶,¹⁷ ID, Rychard Bouwens⁹ ID, Thierry Contini¹⁸ ID, Paulo C. Cortes¹⁹,²⁰ ID, Emanuele Daddi²¹ ID,
Tanio Diaz-Santos⁷,²²,²³ ID, Jorge González-López⁷ ID, Joseph Hennawi²⁴ ID, Jacqueline A. Hodge⁹ ID, Hanae Inami²⁵ ID,
Rob Ivison⁴ ID, Pascal Oesch²⁶,²⁷ ID, Mark Sargent²⁸ ID, Paul van der Werf²⁹, Jeff Wagg³⁰, and L. Y. Aaron Yung⁵

¹Max Planck Institute for Astronomy, Königstuhl 17, D-69117 Heidelberg, Germany; walter@mpia.de

²National Radio Astronomy Observatory, Pete V. Domenici Array Science Center, P.O. Box O, Socorro, NM 87801, USA

³INAFOsservatorio di Astrofisica e Scienza dello Spazio, via Gobetti 93/3, I-40129, Bologna, Italy

⁴European Southern Observatory, Karl-Schwarzschild-Strasse 2, D-85748, Garching, Germany

⁵Rutgers University, 136 Frelinghuysen Road, Piscataway, NJ 08854-8019, USA

⁶Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA

⁷Núcleo de Astronomia, Facultad de Ingeniería y Ciencias, Universidad Diego Portales, Av. Ejército 441, Santiago, Chile

⁸Argelander-Institut für Astronomie, Universität Bonn, Auf dem Hügel 71, D-53121 Bonn, Germany

⁹Leiden Observatory, Leiden University, P.O. Box 9513, NL-2300 RA Leiden, The Netherlands

¹⁰Institut d'Astrophysique de Paris, Sorbonne Université, CNRS, UMR 7095, 98 bis Blvd. Arago, F-75014 Paris, France

¹¹International Centre for Radio Astronomy Research, The University of Western Australia, 35 Stirling Highway, Crawley WA 6009, Australia

¹²Cornell University, 220 Space Sciences Building, Ithaca, NY 14853, USA

¹³Centre for Extragalactic Astronomy, Durham University, Department of Physics, South Road, Durham DH1 3LE, UK

¹⁴Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany

¹⁵Instituto de Astrofísica, Facultad de Física, Pontificia Universidad Católica de Chile Av. Vicuña Mackenna 4860, 782-0436 Macul, Santiago, Chile

¹⁶Millennium Institute of Astrophysics (MAS), Nuncio Monseñor Sótero Sanz 100, Providencia, Santiago, Chile

¹⁷Space Science Institute, 4750 Walnut Street, Suite 205, Boulder, CO 80301, USA

¹⁸Institut de Recherche en Astrophysique et Planétologie (IRAP), Université de Toulouse, CNRS, UPS, F-31400 Toulouse, France

¹⁹Joint ALMA Office, Alonso de Cordova 3107, Vitacura, Santiago, Chile

²⁰National Radio Astronomy Observatory, Charlottesville, VA 22903, USA

²¹Laboratoire AIM, CEA/DSM-CNRS-Université Paris Diderot, Irfu/Service d'Astrophysique, CEA Saclay, Orme des Merisiers, F-91191 Gif-sur-Yvette cedex, France

²²Chinese Academy of Sciences South America Center for Astronomy (CASSACA), National Astronomical Observatories, CAS, Beijing 100101, People's Republic of China

²³Institute of Astrophysics, Foundation for Research and Technology-Hellas (FORTH), Heraklion, GR-70013, Greece

²⁴Department of Physics, Broida Hall, University of California, Santa Barbara, CA 93106-9530, USA

²⁵Hiroshima Astrophysical Science Center, Hiroshima University, 1-3-1 Kagamiyama, Higashi-Hiroshima, Hiroshima, 739-8526, Japan

²⁶Department of Astronomy, University of Geneva, Ch. des Maillettes 51, 1290 Versoix, Switzerland

²⁷International Associate, Cosmic Dawn Center (DAWN) at the Niels Bohr Institute, University of Copenhagen and DTU-Space, Technical University of Denmark, Copenhagen, Denmark

²⁸Astronomy Centre, Department of Physics and Astronomy, University of Sussex, Brighton, BN1 9QH, UK

²⁹Leiden Observatory, Leiden University, P.O. Box 9513, NL-2300 RA Leiden, The Netherlands

³⁰SKA Organization, Lower Withington Macclesfield, Cheshire SK11 9DL, UK

Received 2019 December 20; revised 2020 August 13; accepted 2020 August 13; published 2020 October 19

## Abstract

We combine the recent determination of the evolution of the cosmic density of molecular gas (H₂) using deep, volumetric surveys, with previous estimates of the cosmic density of stellar mass, star formation rate and atomic gas (HⅠ), to constrain the evolution of baryons associated with galaxies averaged over cosmic time and space. The cosmic HⅠ and H₂ densities are roughly equal at z ~ 1.5. The H₂ density then decreases by a factor $\delta_{+2}^{\pm 3}$ to today's value, whereas the HⅠ density stays approximately constant. The stellar mass density is increasing continuously with time and surpasses that of the total gas density (HⅠ and H₂) at redshift z ~ 1.5. The growth in stellar mass cannot be accounted for by the decrease in cosmic H₂ density, necessitating significant accretion of additional gas onto galaxies. With the new H₂ constraints, we postulate and put observational constraints on a two-step gas accretion process: (i) a net infall of ionized gas from the intergalactic/circumgalactic medium to refuel the extended HⅠ reservoirs, and (ii) a net inflow of HⅠ and subsequent conversion to H₂ in the galaxy centers. Both the infall and inflow rate densities have decreased by almost an order of magnitude since z ~ 2. Assuming that the current trends continue, the cosmic molecular gas density will further decrease by about a factor of two over the next 5 Gyr, the stellar mass will increase by approximately 10%, and cosmic star formation activity will decline steadily toward zero, as the gas infall and accretion shut down.

Unified Astronomy Thesaurus concepts: Galaxy evolution (594); High-redshift galaxies (734); Interstellar medium (847); Molecular gas (1073)

<figure><img src="image_3.png" /></figure>
# 1. Introduction

The principal goal in galaxy evolution studies is to understand how the cosmic structure and galaxies that we see today emerged from the initial conditions imprinted on the cosmic microwave background (CMB). In the hierarchical structure formation paradigm, galaxies grow both through the smooth accretion of dark matter and baryons, and through distinct mergers of dark matter halos (and their associated baryons). The accretion of gas eventually leads to the formation of stars in galaxies in the centers of the individual dark matter halos (e.g., White & Rees 1978; Blumenthal et al. 1984; White & Frenk 1991). The winds, UV photons, and supernovae from the ensuing star formation, along with possible episodic accretion onto the supermassive black hole at the center (active galactic nuclei), provide effective “feedback” to the surrounding gas. This may—at least temporarily—suppress the formation of further stars, or may even expel the cold gas from the centers of the potential wells (e.g., Dekel & Silk 1986; Silk & Rees 1998; Croton et al. 2006; Somerville et al. 2008). Together, this leads to a baryon cycle through different gas phases and galactocentric radii (e.g., Tumlinson et al. 2017). Of particular interest in this baryon cycle is the question: how much gas was present both within and around galaxies to explain the formation of stars in galaxies through cosmic times?

Over the past decades, deep sky surveys of star formation and stars in the optical and (near-)infrared bands have put tight constraints on the buildup of the stellar mass in galaxies from early cosmic times to the present (e.g., review by Madau & Dickinson 2014). In parallel, the atomic hydrogen content has been derived through H I emission in the local universe (e.g., Zwaan et al. 2005), and quasar absorption spectroscopy at high redshift (e.g., Prochaska & Wolfe 2009). The molecular gas content of galaxies, the immediate fuel for star formation, has now also been constrained as a function of redshift through measurements of the molecular transitions of carbon monoxide, CO, as well as the far-infrared dust continuum (e.g., reviews by Carilli & Walter 2013; Tacconi et al. 2020; Hodge & da Cunha 2020; Peroux & Howk 2020). These include recent measurements from the ALMA Spectroscopic Survey in the Hubble Ultra-deep Field (ASPECS; Decarli et al. 2019, 2020; Magnelli et al. 2020). Together, the available data have now reached the point that we can account for the total cold gas content (H I and H₂) that is associated with galaxies as a function of cosmic time.

In this paper we discuss how these new molecular gas constraints impact our view of the cosmic baryon cycle of galaxies, and in particular how they affect our view of gas accretion to sustain the observed star formation rate density in the centers of galaxies. Throughout this paper we only consider densities that are averaged over cosmic space and wide time bins to characterize the cosmic baryon cycle (Section 3). We argue that such an approach is justified as molecular gas, star formation, and stellar mass are found to be approximately spatial in galaxies, and because the averaging times are significantly longer than the physical processes under consideration (Section 2). We thus stress that many conclusions of this paper, including the accretion and inflow rates, will not be applicable to individual galaxies, but only to volume-averaged galaxy samples (Section 4). Given the available observational constraints, we here focus on redshifts below $z \sim 4$ (when the universe was older than 1.5 Gyr).

We adopt a “cosmic concordance cosmology” with the following parameters: a reduced Hubble constant $h = H_0/(100 \text{ km s}^{-1} \text{ Mpc}^{-1}) = 0.7$, a matter density parameter

$\Omega_m = 0.31$ (which is the sum of the dark matter density parameter $\Omega_c = 0.259$ and the baryon density parameter $\Omega_b = 0.048$), and a dark energy density parameter $\Omega_\Lambda = (1 - \Omega_m) = 0.69$, similar to Planck constraints (Planck Collaboration et al. 2016), and those used in the review on cosmic star formation rates and associated stellar mass buildup by Madau & Dickinson (2014). All the volume-averaged, cosmological densities quoted in this paper are in comoving units.

# 2. A Simple Schematic

Figure 1 shows a schematic of the different baryonic components that are present within the dark matter halo of a galaxy. The central region of the galaxy contains the majority of the stars, molecular gas, and star formation at any given time (Sections 3.2.1, 3.2.3). In this region, stars form out of giant molecular clouds with a typical timescale of order $10^7$ yr (e.g., Kawamura et al. 2009; Meidt et al. 2015; Schinnerer et al. 2019) and molecular gas is expected to form out of atomic gas on a similar timescale (depending on metallicity, e.g., Fukui et al. 2009; Glover & Mac Low 2011; Clark et al. 2012; Walch et al. 2015). These periods are significantly shorter than the Gyr-averaged timescales discussed in this study (Section 4).

Throughout this paper the term “disk” is used to define this region (with a typical$^{31}$ radius $r_{\text{stars}} < 10 \text{ kpc}$). Note that the term “disk” should not be taken literally: for example, low-mass galaxies may not form well-defined disks, and many massive disk galaxies will transition to elliptical galaxies through mergers over time. We thus consider the “disk” nomenclature to define the main stellar components of galaxies, which, for main-sequence star-forming galaxies at high redshift, can be considered disklike in many cases (e.g., Förster Schreiber et al. 2009; Wuyts et al. 2011; Law et al. 2012; Salmi et al. 2012).

This nominal “disk” region is surrounded by a reservoir of atomic gas (H I) with radii $r_{\text{HI}} < 50 \text{ kpc}$ (Section 3.2.2), as demonstrated by observations in the local universe (e.g., Walter et al. 2008; Leroy et al. 2009), high-redshift observations (Krogager et al. 2017; Neeleman et al. 2017, 2019) and simulations (e.g., Bird et al. 2014; Rahmati & Schaye 2014). Outside the atomic gas region is the circumgalactic medium (CGM), defined to be located within the virial radius ($r_{\text{vir}} \sim 50–300 \text{ kpc}$), meaning gravitationally bound to the dark matter halo, and decoupled from the expansion of the universe (e.g., Tumlinson et al. 2017). The CGM consists of predominantly ionized gas at a range of temperatures ($T \sim 10^4–10^6 \text{ K}$). The timescale to accrete material from the cool $T \sim 10^4 \text{ K}$ CGM is comparable to the dynamical time ($\sim 10^8$ yr), orders of magnitudes shorter than the cooling time of the hot $T \sim 10^6 \text{ K}$ CGM ($\gg 1$ Gyr, Section 4.5.2). The medium outside this gravitationally collapsed/bound structure (i.e., beyond $r_{\text{vir}}$) is referred to as the intergalactic medium (IGM).

The above defined regions are not static, and gas can be exchanged between these regions. The most important gas flows are also included in the schematic shown in Figure 1, i.e., outflows as well as gas accretion. As detailed below (Section 4.3), the accretion process can be described as: (i) the net infall of ionized material from the CGM and/or IGM onto the extended H I reservoir, and (ii) the net inflow of H I

<footnote>$^{31}$ The physical scales quoted here in kiloparsecs are only given as examples for typical $M_*$ star-forming galaxies, and will scale as a function of the actual mass of a given dark matter halo. For a dependence of $r_{\text{stars}}$ on $r_{\text{vir}}$, see, e.g., Somerville et al. (2018).</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Schematic of the different baryonic components that are present within the dark matter halo of a galaxy (defined as $r < r_{vir}$). The central “disk” region ($r < r_{stars}$), contains the vast majority of stars and molecular gas, and stars form here at a rate $\psi_{stars}$. This region is surrounded by a reservoir of atomic gas (H I), with $r < r_{HI}$. The predominantly ionized material (H II) beyond this radius, but within $r_{vir}$, constitutes the circumgalactic medium (CGM). Beyond $r_{vir}$ is the realm of the intergalactic medium (IGM). Blue arrows indicate the (net) infall of ionized gas to the H I reservoir ($\tilde{\psi}_{HII \to HI}$) as well as the (net) inflow of atomic gas to the molecular gas (H$_2$) reservoir ($\tilde{\psi}_{HI \to H2}$). The red arrow indicates the material entrained in outflows that can reach the CGM and possibly the IGM (here assumed to be proportional to $\psi_{stars}$).</figcaption></figure>

from the H I reservoir (within $r_{HI}$), with the subsequent conversion to H$_2$, onto the central region of the galaxy (within $r_{stars}$). We also note that our schematic does not include the accretion of mass through galaxy mergers. Their contribution to the mass buildup in galaxies is significantly smaller than that from accretion (e.g., van de Voort et al. 2011).

We emphasize that the demarcation of IGM versus CGM versus “disk” is not a simple geometric one, with material necessarily transitioning from one region to the other over time. For instance, the H I and warm/hot halo gas may mix substantially through streams, Galactic fountains and outflows, as well as filaments. Likewise, many galaxies reside in groups or clusters, where the dark matter halos may overlap, and defining whether gas is in the IGM versus CGM may be ambiguous. However, for the purpose of the analysis presented in this paper, where we focus on the evolution of the baryonic components of the “disk” structure, the proposed simple schematic in Figure 1 should suffice as a representative guide.

### 3. Mass Components

To put the different baryonic mass components in galaxies in context, we compile current literature estimates of their “cosmic mean density” as a function of redshift. The total number of baryons is conserved over time, and therefore, by definition, the density of baryons does not change with time when considering comoving volumes.32

#### 3.1. The $z \sim 0$ Census

For the low-redshift universe ($z \lesssim 0.3$), an almost complete census of the baryons is available (e.g., Shull et al. 2012; Tumlinson et al. 2017; Nicastro et al. 2018). The latest studies place the large majority ($\sim$82%) of the cosmic baryons in the IGM (e.g., Shull et al. 2012). These baryons are highly ionized (temperatures between $10^5$ and $10^7$ K), and detected via O VI and O VII absorption features, and as the Ly$\alpha$ forest. The distribution is thought to be highly filamentary, with the majority of the IGM residing in the “cosmic web.” Recent work on fast radio bursts has shown promise to detect this hard-to-trace component (e.g., McQuinn 2014; Shull & Danforth 2018; Macquart et al. 2020).

The remaining 18% of the baryons at $z \approx 0$ then belong to the “collapsed phase” (Shull et al. 2012, see also their Figure 10), gravitationally bound to galaxies, groups, and clusters that we will discuss in the following. The hot intercluster medium (ICM; $\ge 10^7$ K), seen in X-rays, comprises 4% of the cosmic total.33 The stars in all types and masses of galaxies comprise 7% of the total baryon density. The cold gas (H I and H$_2$) comprises a little more than one percent at $z = 0$, $\sim$85% of which is in H I. The CGM (also called “hot halos”), comprises about 5% of the cosmic total, although again, the exact demarcation of the CGM remains somewhat subjective (see also Shull et al. 2012; Tumlinson et al. 2017, and Section 2).

There are other mass components in galaxies, but they only marginally contribute to the total mass budget, as briefly summarized in the following. As their combined contribution is of the order a few percent, we do not consider them further in our analysis.

**Warm ionized medium:** the warm ionized medium (WIM) is visible in H$\alpha$ and X-rays, and makes up less than 1% of the total baryon mass in galaxy disks (e.g., Anderson & Bregman 2010; Putman et al. 2012; Werk et al. 2014).

**Black holes:** the majority of galaxies are thought to host central supermassive black holes (SMBH). Various studies put this ratio at $\sim$0.1% of the total stellar mass in galaxies (Kormendy & Ho 2013). The remnants of massive stars are by definition included in the stellar initial mass function (IMF) determinations (e.g., IMF review by Bastian et al. 2010). Some black holes may be ejected entirely from galaxies via interactions with other black holes, but this net mass effect is minor (e.g., Loeb 2007).

**Dust:** although dust plays a central role in the formation of stars, dust only makes up about 1% of the total ISM mass (e.g., Sandstrom et al. 2013). The cosmic evolution of the dust content in the universe has recently been discussed in Driver et al. (2018) and Magnelli et al. (2020).

#### 3.2. Redshift Evolution

In the following we discuss the key baryonic mass components in galaxies, and their evolution with cosmic time.

##### 3.2.1. Star Formation and Stars

The evolution of the cosmic star formation rate density ($\psi_{stars}$; Figure 2, top left) has been constrained through various multiwavelength studies of large samples of individual galaxies over the last decades (as summarized in the review by

<footnote n="32">32 Strictly speaking, the baryon density decreases with time due to fusion, as some of the mass is converted to energy, e.g., in the case of the fusion of two hydrogen atoms to form helium, 0.7% of the mass is lost to radiation. During a complete CNO cycle, approximately the same amount of energy is being released. As only a small fraction of all baryons, those within the centers of stars, take part in the fusion process, we ignore this mass loss here.</footnote>

<footnote n="33">33 The ICM is not labeled in the schematic shown in Figure 1 as it only applies for cluster environments.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Redshift evolution of different baryonic components in galaxies compiled from the literature. The measurements of the cosmic star formation rate density (top left) and stellar mass density (top right) are from the compilation in Madau & Dickinson (2014; their Tables 1 and 2). The solid line is the best-fit functional form to the data (Section 3.3) with the parameters given in Table 1, and the shaded region marks the 1σ region (16th to 84th percentile) from a Monte Carlo Markov Chain analysis (see Section 3.3 for details). The orange dashed line in the cosmic stellar mass density panel is the integration of the best-fit function form to the star formation rate density. The discrepancy between this curve and the measurements is described by the return fraction (see the text and Madau & Dickinson 2014). Observational constraints on $\rho_{HI}$ (bottom left) are from a compilation given in Neeleman et al. (2016) updated with some recent constraints at low redshift (Section 3.2.2). Gray points indicate measurements at <6σ and black points are measurements at >6σ (see Appendix B). Constraints on $\rho_{H2}$ (bottom right) are from ASPECS (Decarli et al. 2019, 2020) and other CO surveys (black points; see Appendix B). Gray points indicate measurements obtained through dust continuum observations (Appendix B; including the ASPECS measurements by Magnelli et al. 2020).</figcaption></figure>

Madau & Dickinson 2014). Early studies were based on rest-frame UV observations (e.g., Lilly et al. 1996; Madau et al. 1996; Bouwens et al. 2012a, 2012b; Cucciati et al. 2012, see the above review for a complete list of references), and are complemented through observations at longer wavelengths (e.g., Magnelli et al. 2011, 2013; Gruppioni et al. 2013; Sobral et al. 2013; Bouwens et al. 2016, 2020; Novak et al. 2017; Dudzevičiūtė et al. 2020; Khusanova et al. 2020). These estimates indicate that the peak of cosmic star formation occurred at $z \sim 2$, with a subsequent decline by a factor of $\sim 8$

to the present day. Integrating $\psi_{stars}$ gives the stellar mass formed at a given cosmic time, and this integral is shown as a dashed orange line in Figure 2 (top right).

This integral can be compared to the independently measured stellar mass density $\rho_{stars}$ (shown as a red line in Figure 2, top right). This stellar mass density has been determined by numerous studies (e.g., Pérez-González et al. 2008; Marchesini et al. 2009; Caputi et al. 2011; Ilbert et al. 2013; Muzzin et al. 2013), as compiled and homogenized in the review by Madau & Dickinson (2014). Both the stellar mass
and the star formation rates depend on the choice of the IMF, and SED-fitting method (e.g., Bastian et al. 2010; Kennicutt & Evans 2012; Leja et al. 2020).<fnref n="34" />

This temporal integral of the star formation rate density lies above the measured stellar mass density $\rho_{\text{stars}}$ by a factor $1.4 \pm 0.1$. This is due to the fact that not all stellar mass that is formed will stay locked in stars; some fraction will be returned to the ISM, CGM, or IGM (depending on the mass of the galaxy). The cosmic-averaged star formation rate density $\psi_{\text{stars}}(z)$ is thus the first time derivative of $\rho_*(z)$, modulo the return fraction<fnref n="35" /> of stars R to the interstellar medium through stellar winds and/or supernova explosions (e.g., Madau & Dickinson 2014), i.e.,

$$ \dot{\rho}_{\text{stars}}(z) = (1 - R) \psi_{\text{stars}}(z). $$

The fact that the integral of $\psi_{\text{stars}}$, after accounting for the return fraction, is in reasonable agreement with $\rho_{\text{stars}}$ is remarkable, as highlighted in Madau & Dickinson (2014), if one considers the number of assumptions that go into each measurement.<fnref n="36" /> These mass estimates do not include stars that are found outside galaxy disks, e.g., in stellar streams around galaxies, and the intracluster environment. This stellar mass component, however, only constitutes a small fraction of the stellar mass present in the galaxy disks (e.g., Behroozi et al. 2013), and we therefore do not consider this component further. For completeness it should be noted that some of the stellar mass growth can occur through mergers of galaxies, but this gain (“ex-situ”) through merging of the existing stellar masses is small compared to the actual star formation process (“in situ”), at least for galaxies around $L^*$ (e.g., Behroozi et al. 2019).

### 3.2.2. Atomic Gas

The evolution of the cosmic density of atomic gas associated with galaxies ($\rho_{\text{HI}}(z)$; Figure 2, bottom left) has been constrained with several different approaches depending on redshift range. At $z \approx 0$, large surveys aimed at measuring the H I 21 cm emission from local galaxies can constrain the H I mass function (e.g., Zwaan et al. 2005; Braun 2012; Jones et al. 2018) whose integral provides an estimate of $\rho_{\text{HI}}$. At higher redshifts ($0.3 \lesssim z \lesssim 1$), where the H I 21 cm emission becomes increasingly faint and therefore single sources are below the detection threshold of current radio-wavelength facilities, stacking of H I 21 cm emission from a large sample of galaxies provides an alternative approach to measuring $\rho_{\text{HI}}(z)$ (e.g., Lah et al. 2007; Delhaize et al. 2013; Rhee et al. 2013; Kanekar et al. 2016; Bera et al. 2019). In addition, the cross-correlation between 21 cm intensity maps and the large-scale structure (so-called 21 cm intensity mapping) provide an independent measurement of $\rho_{\text{HI}}$ at these redshifts (e.g., Masui et al. 2013; Switzer et al. 2013).

At $z \gtrsim 1.6$, H I can be observed using ground-based optical telescopes through its Ly$\alpha$ transition. Quasar absorption spectroscopy of the strongest Ly$\alpha$ absorbers, the so-called damped Ly$\alpha$ systems (DLAs; Wolfe et al. 2005) has yielded

estimates of $\rho_{\text{HI}}$ up to $z \sim 5$ (e.g., Crighton et al. 2015). The $\rho_{\text{HI}}$ estimate obtained from DLA surveys is simply the total H I column density detected in DLAs divided by the path length of the survey. Here the main uncertainties come from relatively poorly understood systematics between varying methods of measuring DLAs and a potential bias against dusty, high H I column density systems (Ellison et al. 2001; Jorgenson et al. 2006; Krogager et al. 2019). Most numerical simulations predict DLAs to probe gas near galaxies (Rahmati & Schaye 2014), which is supported by observations of the cross-correlation function between DLAs and the Ly$\alpha$ forest (Pérez-Ràfols et al. 2018).

These measurements do not include any contributions from systems below the DLA column density threshold, because these systems contain less than 20% of the total cosmic atomic gas density (Péroux et al. 2003; O'Meara et al. 2007; Noterdaeme et al. 2012; Berg et al. 2019), and their connection with galaxies is less certain. However, we do account for the contribution of helium, which corresponds to a correction factor of $\mu = 1.3$.

The emerging picture is that the cosmic density of neutral atomic gas remains approximately constant with redshift, with a decline by a factor of $\sim 2$ from $z \sim 3$ to $z = 0$. We remind the reader that the H I is coming from a more extended reservoir compared to the stellar mass and star formation measurements of galaxies (see discussion in Section 2).

### 3.2.3. Molecular Gas

A number of approaches have been followed in the past to constrain the evolution of the cosmic molecular gas density ($\rho_{\text{H2}}$; Figure 2, bottom right). Here we focus on methods that are not merely based on stellar mass and star formation rate determinations with subsequent application of scaling relations. In particular, we include the recent results from ASPECS, which perform deep frequency scans to detect redshifted CO lines without any preselection. This approach has been successfully applied in a number of studies (Decarli et al. 2014, 2016, 2019, 2020; Walter et al. 2014, 2016; Pavesi et al. 2018; Klitsch et al. 2019; Riechers et al. 2019; Lenkić et al. 2020). Molecular gas constraints derived from dust emission (frequently using scaling relations based on stellar mass or star formation rates) and other approaches show a consistent evolution (e.g., Berta et al. 2013; Scoville et al. 2017; Driver et al. 2018; Liu et al. 2019; Dudzevičiūtė et al. 2020; Magnelli et al. 2020).

In order to convert CO to H$_2$ measurements, the detected CO emission has to be corrected for excitation and a CO-to-H$_2$ conversion factor has to be applied (that also accounts for helium). The CO-to-H$_2$ conversion factor is the main systematic uncertainty in the analysis. For the ASPECS measurement, the majority of the molecular gas mass density comes from individually detected galaxies (Decarli et al. 2019). Their metallicities (consistent with solar) and stellar masses ($M_{\text{stars}} \gtrsim 10^{10} M_{\odot}$) justify the choice of a Galactic conversion factor to determine the molecular gas mass (Boogaard et al. 2019). The uncertainties in molecular gas excitation, as derived for the ASPECS galaxies in Boogaard et al. (2020), have been anchored based on CO(1–0) observations out to $z \sim 3$ (VLASPECS, Riechers et al. 2020), and were folded into the ASPECS measurements (Decarli et al. 2020). Converting dust measurements to molecular gas masses requires the choice of a

<footnote n="34">³⁴ Madau & Dickinson (2014) assume a Salpeter IMF and a lower fixed threshold in luminosity of 0.03 $L_*$.</footnote>

<footnote n="35">³⁵ The return fraction is $R = 0.27$ for a Salpeter IMF and $R = 0.41$ for a Chabrier IMF that is more weighted toward massive stars (Madau & Dickinson 2014).</footnote>

<footnote n="36">³⁶ But see Hopkins et al. (2018), who argues that this overall agreement does not necessarily imply that the IMF has to be universal.</footnote>
dust temperature, emissivity, and a dust-to-gas ratio (see the above references).

Stacking and intensity mapping techniques (Inami et al. 2020; Uzgil et al. 2019) indicate that the majority of all CO emission in the UDF is captured by the current observations, i.e., the faint-end slope of the CO luminosity functions is such that extrapolating to lower masses would not significantly (less than 50%) increase the total emission (see also Decarli et al. 2020). These high-redshift measurements are anchored at $z = 0$ through detailed studies of the molecular gas content in the local universe (Keres et al. 2003; Boselli et al. 2014; Saintonge et al. 2017; Fletcher et al. 2020).

The emerging picture based on the abovementioned molecular gas and dust studies is that the cosmic density of molecular gas decreased by a factor of $6^{+3}_{-2}$ from the peak of cosmic star formation ($z \sim 2$) to today (see also recent reviews by Hodge & da Cunha 2020; Peroux & Howk 2020; Tacconi et al. 2020). There is evidence that the molecular gas density increased from $z \sim 6$ to $z \sim 2$ (Riechers et al. 2019; Decarli et al. 2019, 2020), but the associated uncertainties are significant for $z > 3$.

### 3.3. Fitting Functions

In order to capture the global trends in the cosmic density measurements discussed in the previous paragraphs, we have fitted the observational data with functional forms (data given in Appendix B). In particular, for $\rho_{HI}$, $\rho_{stars}$, and $\psi_{stars}$, we adopt a smooth double power law, similar to that defined in Madau & Dickinson (2014):

$$ \rho_{x}(z) = \frac{A(1 + z)^{B}}{1 + [(1 + z)/C]^{D}}. \qquad (1) $$

In order to capture the apparent flattening of the evolution of $\rho_{HI}$ at both low and high redshift, we adopted a hyperbolic tangent function (as in Prochaska & Neeleman 2018):

$$ \rho_{HI}(z) = A \tanh(1 + z - B) + C. \qquad (2) $$

These functional forms are not physically motivated and are simply meant to capture the general trends of the data points. To estimate the best-fit parameters and associated uncertainties, we fit the data using a Monte Carlo Markov Chain approach utilizing the `emcee` package (Foreman-Mackey et al. 2013). For all cosmic densities, we marginalize over a nuisance parameter to account for intrinsic scatter within the data points not accounted for by the uncertainties of the individual points. To take into account systematic uncertainties within the varying data sets, we symmetrically (in log scale) increase the formal uncertainties derived from the fitting procedure such that >68% of all measurements are contained within the $1\sigma$ boundaries (16th to 84th percentile). The best fits are shown as solid lines in Figures 2 and 3, whereas the $1\sigma$ boundaries of the fitting functions are shown as colored regions. The fitting parameters are summarized in Table 1. We note that a fit to $\rho_{HI}$ based on just the ASPECS data gives almost identical parameters as those shown in Table 1.

### 3.4. Cosmic Averages

In the analysis that follows, we will consider the above volume-averaged measurements (Section 3.2) to derive volume-averaged properties (such as depletion times, gas accretion rates). The fundamental assumption is that, statistically speaking, the galaxies are similar to the picture discussed

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Census of baryons inside and outside galaxies using the fitting functions shown in Figure 2 and presented in Section 3.3. The colors have the same definitions as those in Figure 2. The orange line shows the sum of the $H_1$ and $H_2$ components, whereas the black line shows the sum of all of the baryons (stars, $H_1$ and $H_2$) associated with galaxies. The dotted line is the total cosmic baryon content and the dashed line is the total dark matter content for the given $\Lambda$-CDM universe. The same curves are plotted as a function of (linear) time in Figure 7.</figcaption></figure>

in Section 2 and Figure 1. One can express the quantities discussed here as a function of the well-characterized stellar mass function (SMF) $\Phi_{\star}(z, M)$ (e.g., Davidzon et al. 2017). Then the cosmic stellar mass density can be written as

$$ \rho_{\star}(z) = \int \Phi_{\star}(z, M_{\star}) dM_{\star}, $$

where $M_{\star}$ is the stellar mass. The gas ($H_2$ or $H_1$) density can then be expressed as

$$ \rho_{\text{gas}}(z) = \int \Phi_{\star}(z, M_{\star}) \times f_{\text{gas}}(z, M_{\star}) dM_{\star}, $$

where $f_{\text{gas}}$ is the gas-to-stellar mass fraction ($f_{H_2}$ or $f_{HI}$).

By definition, these functions are volume averages that marginalize over dependencies of baryonic components on other parameters (such as, e.g., environment, metallicity, and feedback processes).

## 4. Discussion

We now discuss the density evolution of the various mass components in the universe and implications for gas accretion rates. As stressed before, our measurements are volume- and time-averaged. The timescales of the individual mass conversion processes ($\lesssim 0.1$ Gyr, Section 2) are smaller than the cosmic timescales over which we are averaging ($\Delta z = 1$
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Left: the ratio of cosmic molecular to atomic gas density as a function of redshift. The ratio peaks at $z \sim 1.5$, close to the peak of the star formation rate density. Middle: the ratio of molecular gas to stellar mass density as a function of redshift. Right: cosmic gas depletion timescale, defined as the density in molecular gas divided by the cosmic star formation rate density. The gray dashed curve is the Hubble time vs. redshift. In all panels the thick solid line is derived from the functional form to the data (Section 3.3) with the parameters given in Table 1. The shaded region marks the 1σ region (16th to 84th percentile) of all the curves from a Monte Carlo Markov Chain analysis. Thin lines show several random realizations of this analysis.</figcaption></figure>

<table><caption>**Table 1**

Fitting Functions to the Observed Cosmic Density Measurements Shown in Figure 2</caption><thead><tr><th></th><th>Fitting function</th><th>A</th><th>B</th><th>C</th><th>D</th></tr></thead><tbody><tr><td>ρ<sub>H2</sub>(z)[M<sub>☉</sub> Mpc<sup>-3</sup>]</td><td>Equation (1)</td><td>(1.00 ± 0.14) × 10<sup>7</sup></td><td>3.0 ± 0.6</td><td>2.3 ± 0.3</td><td>5.1 ± 0.5</td></tr><tr><td>ρ<sub>stars</sub>(z)[M<sub>☉</sub> Mpc<sup>-3</sup>]</td><td>Equation (1)</td><td>(1.3<sup>+1.0</sup><sub>-0.6</sub>) × 10<sup>10</sup></td><td>-4.1 ± 0.4</td><td>2.5 ± 0.4</td><td>-3.8 ± 0.3</td></tr><tr><td>ψ<sub>stars</sub>(z)[M<sub>☉</sub> yr<sup>-1</sup> Mpc<sup>-3</sup>]</td><td>Equation (1)</td><td>0.0158 ± 0.0010</td><td>2.88 ± 0.16</td><td>2.75 ± 0.11</td><td>5.88 ± 0.15</td></tr><tr><td>ρ<sub>HI</sub>(z)[M<sub>☉</sub> Mpc<sup>-3</sup>]</td><td>Equation (2)</td><td>(4.5 ± 0.5) × 10<sup>7</sup></td><td>2.8 ± 0.4</td><td>(1.01 ± 0.07) × 10<sup>8</sup></td><td>...</td></tr></tbody></table>

corresponds to a time period of ~0.6 Gyr at $z = 3.5$, ~2.5 Gyr at $z = 1.5$, and ~5.5 Gyr at $z = 0.5$). Therefore, our conclusions will not be applicable to all individual galaxies.

## 4.1. The Evolution of the Cosmic Baryon Density

Figure 3 summarizes the evolution of the baryon content in stars, H I, and H$_2$ associated with galaxies together with the cosmic dark matter and the total baryon density. As discussed in Section 3, the large discrepancy between the total baryon density (dotted curve in Figure 3) and the baryon density inside galaxies $\rho_{bar,gal}$ (black curve in Figure 3) indicates that most baryons are not inside galaxies, but are in the predominantly ionized IGM (and CGM). The stellar mass density is increasing continuously with time, and surpasses that of the total gas density (H I and H$_2$) at redshift $z \sim 1.5$.

In Figure 4 (left) we plot the ratio of molecular to atomic gas density as a function of redshift. This ratio peaks at $z \sim 1.5$, close to the peak of the star formation rate density. Figure 4 (middle) shows the ratio of molecular gas to stellar mass as a function of redshift. At redshifts $z \lesssim 2$ the stellar mass density starts to dominate over the molecular gas density.

The last panel in Figure 4 (right) shows the molecular gas depletion time, i.e., how long will it take to deplete the molecular gas reservoir at the current rate of star formation. The depletion time ($\rho_{H2}/\psi_{stars}$) is approximately constant above redshifts $z \gtrsim 2$, and then increases slightly from $\tau_{depl} \sim (4 \pm 2) \times 10^8$ yr at $z \sim 2$ to $\tau_{depl} = (7 \pm 3) \times 10^8$ yr at $z = 0$, and is shorter than the Hubble time at all redshifts. This immediately implies that the molecular gas reservoir needs to be continuously replenished (i.e., through accretion). Both

the ratio of molecular gas to stellar mass and the depletion times for the molecular gas phase are similar to what is found in scaling-relation studies of individual galaxies (e.g., Daddi et al. 2010; Genzel et al. 2010; Bothwell et al. 2013; Tacconi et al. 2018; Aravena et al. 2019, 2020).

## 4.2. The Need for Accretion

The need for gas accretion onto galaxies from the cosmic web to sustain the observed star formation activity has been noted numerous times before (e.g., Bouché et al. 2010; Bauermeister et al. 2010; Davé et al. 2012; Behroozi et al. 2013; Béthermin et al. 2013; Conselice et al. 2013; Lilly et al. 2013; Tacconi et al. 2013; Peng & Maiolino 2014; Rathaus & Sternberg 2016; Scoville et al. 2017; Tacconi et al. 2018). Prior to the availability of direct measurements of the H$_2$ density it was occasionally argued that, given the approximate constancy of the H I density through cosmic time, the net gas accretion rate density needed to be approximately equal to the star formation rate density. Now that the molecular density is directly observed, this topic can be revisited (see also the recent reviews by Hodge & da Cunha 2020; Peroux & Howk 2020; Tacconi et al. 2020).

We first ask how much stellar mass could in principle be formed by looking at the decrease in the molecular gas density since the peak of the cosmic molecular gas density at $z \sim 1.5$. If we assume that the net loss in H$_2$ since that time is fully due to the formation of stellar mass, we can derive the maximum stellar mass growth due to this conversion. This is shown in Figure 5 as the blue curve. For completeness, we also show the loss in H I (orange line: sum of H I and H$_2$ loss) that eventually may also end up as stellar mass over the same cosmic time via a transition
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Cumulative gain of the stellar mass density (red line) compared to the cumulative loss of the gas mass density (H₂: blue line, total gas: orange line), starting at a redshift of z = 1.5 (T<sub>Univ</sub> ~ 4 Gyr), i.e., approximately the peak of the molecular gas density. The lower (upper) abscissa shows cosmic age (redshift). Even assuming that 100% of the gas will end up in stars, the gas observations cannot account for the observed stellar mass buildup. The remaining mass to build up the stellar mass must be accreted onto the galaxy.</figcaption></figure>

through the molecular gas phase (Section 4.3). We compare this to the total observed gain in stellar mass over the same cosmic time (shown as the red curve in Figure 5, based on the red curve show in Figure 2). Even assuming that all the molecular gas ends up in stars, the observed decline in H I and H₂ is only able to account for ~25% of the total stellar mass formed during this time. Also note that the above stellar mass measurement ignores the return of stellar mass to the ISM, CGM, and IGM. If this additional stellar mass is accounted for, the observed H I and H₂ can only account for ~20% of the total stellar mass formed. The difference in mass is thus the minimum amount of material that needs to be accreted by the galaxies from the IGM/CGM since the universe was 4 Gyr old.

### 4.3. H II Infall and H I Inflow Rates

Most of the stars are thought to form out of H₂ and not atomic hydrogen (e.g., Schruba et al. 2011), at least at the redshifts considered in this paper. However, the presence of H I is a prerequisite to form H₂. In nearby galaxies it is found that H I is significantly more extended than the stellar component, which also harbors most of the star formation and the H₂ (Walter et al. 2008; Leroy et al. 2009). At high redshift, the situation is likely very similar, as indicated by the fact that the impact parameter for the DLAs found in quasar spectra are ≤50 kpc (Section 3.2.2), whereas the stellar components are typically ≤10 kpc in size (e.g., Fujimoto et al. 2017; Elbaz et al. 2018; Jiménez-Andrade et al. 2019). The fact that DLAs show little to no Lyman-Werner absorption from molecules also points toward the fact that the H I is more extended than the H₂, i.e., that the DLAs contain negligible molecules (Noterdaeme et al. 2008; Jorgenson et al. 2014; Muzahid et al. 2015).

We consider the accretion of material to the central star-forming “disk” as a two-step process. The first is the net infall of ionized gas (H II) onto the extended H I reservoir, $\tilde{\psi}_{HII \to HI}$, e.g., through cold mode accretion (Section 4.5). In a second step the gas further cools and settles in the central region where it forms H₂, which we refer to as net inflow, $\tilde{\psi}_{HI \to H2}$. We stress that we can only consider net rates: it is also possible that H₂ (or H I) is dissociated/photoionized to form H II through feedback processes. Our data do not allow us to differentiate between inflows and outflows, and we here define the net flow rates in such a direction that they are likely positive, i.e., $\tilde{\psi}_{HII \to HI} > 0$, $\tilde{\psi}_{HI \to H2} > 0$. We note that strictly speaking we refer to net flow rate densities (averaged over cosmic volume) throughout this work. For simplicity, however, we refer to these as rates throughout.

As detailed in Appendix A, the rate at which the observed H₂ density $\rho_{H2}$ is used up for star formation, lost due to feedback (both stellar or AGN) to the CGM, and is being replenished by H I can be written as

$$ \dot{\rho}_{H2}(z) = - \frac{\tilde{\psi}_{\text{stars}}(z)}{\text{star formation rate}} - \frac{\xi \tilde{\psi}_{\text{stars}}(z)}{\text{H2 loss due to feedback}} + \frac{\tilde{\psi}_{\text{HI} \to \text{H2}}(z)}{\text{H2 gain from HI reservoir}} \quad (3) $$

where $\psi_{\text{stars}}$ is the star formation rate density, and $\tilde{\psi}_{HI \to H2}$ is the net conversion rate of H I to H₂; $\dot{\rho}_{H2}(z)$ is the time derivative of $\rho_{H2}(z)$. The (unknown) mass loading factor $\xi$ accounts for mass loss due to outflow driven by active star formation and AGN activity that is a function of the environment and mass distribution(s) within a galaxy. We simplistically assume that this outflow/mass loading is linearly correlated with the star formation rate density, $\psi_{\text{stars}}$ (e.g., Spilker et al. 2018; Schroetter et al. 2019) with a universal proportionality factor $\xi$.

The material that is required to replenish the H I reservoir ($\dot{\rho}_{HI}(z)$ being the time derivative of $\rho_{HI}(z)$) can be expressed as

$$ \dot{\rho}_{HI}(z) = - \frac{\tilde{\psi}_{HI \to H2}(z)}{\text{loss to H2}} + \frac{\tilde{\psi}_{HII \to HI}(z)}{\text{net HI gain from HII reservoir}} \quad (4) $$

where $\psi_{HII \to HI}$ is the net infall of gas from the ionized gas phase. As described in Appendix A this expression for $\dot{\rho}_{HI}(z)$ (unlike the one for $\dot{\rho}_{H2}(z)$) does not include a mass loading term, as it is included in the net flow term $\tilde{\psi}_{HII \to HI}$.

We can solve Equations (3) and (4) for the net inflow rate $\tilde{\psi}_{HI \to H2}$ and the net infall rate $\tilde{\psi}_{HII \to HI}$ as a function of observables $\rho_{HI}(z)$, $\rho_{H2}(z)$ and $\psi_{\text{stars}}$:

$$ \tilde{\psi}_{HI \to H2} = \dot{\rho}_{H2}(z) + (1 + \xi) \psi_{\text{stars}}(z) \quad (5) $$

and

$$ \tilde{\psi}_{HII \to HI} = \dot{\rho}_{HI}(z) + \dot{\rho}_{H2}(z) + (1 + \xi) \psi_{\text{stars}}(z). \quad (6) $$

In Figure 6 we plot these net flows rates (Equations (5) and (6)), along with the star formation rate density $\psi_{\text{stars}}$. We also show the time derivatives of $\rho_{HI}$ and $\rho_{H2}$, derived as the proper time derivatives of the measured relations with redshift, as parameterized in Equations (1) and (2). The differences between the net flow rates and the star formation rate density are due to the building up, or depletion, of gas in the neutral atomic and molecular phase, as dictated by the time derivative curves.

At high redshift ($z > 3$), both the net H II infall rate ($\tilde{\psi}_{HII \to HI}$) and H I inflow rate ($\tilde{\psi}_{HI \to H2}$) are larger than the star
<figure><img src="image_9.png" /><figcaption>**Figure 6.** The H II net infall rate ($\tilde{\psi}_{HII \to HI}$, Equation (6), orange curve) and H I inflow rate ($\tilde{\psi}_{HI \to H2}$, Equation (5), black curve) are plotted together with the cosmic star formation rate density ($\psi_{stars}$, red curve), assuming a mass loading factor of $\xi = 0$. When including feedback/mass loading (i.e., $\xi > 0$), the inflow and accretion rate would have to increase correspondingly, to account for the extra loss of gas. We also show the time derivatives of the H I and H$_2$ densities ($\dot{\rho}_{HI}(z)$ and $\dot{\rho}_{H2}(z)$), as derived from the temporal gradients of the measured density curves in Figure 2, as parameterized in Equations (1) and (2). The curves of $\tilde{\psi}_{HII \to H2}$ and $\tilde{\psi}_{HI \to HI}$ are a linear combination of the measured quantities: $\psi_{stars}$, $\dot{\rho}_{HI}(z)$, and $\dot{\rho}_{H2}(z)$, as per Equations (5) and (6). Below $z \approx 1.5$ the inflow rate $\tilde{\psi}_{HI \to H2}$ drops below $\psi_{stars}$, as the cosmic H$_2$ reservoir is used up to form stars (negative $\dot{\rho}_{H2}(z)$).</figcaption></figure>

formation rate density, which is reflected in the build up of molecular gas over time, with the HI being a pass-through phase (close to zero derivative). At $z \gtrsim 1.5$ the net inflow rate $\tilde{\psi}_{HI \to H2}$ is higher than $\psi_{stars}$. At these redshifts, the H$_2$ cosmic density is still increasing with time. Therefore on top of the flow of H$_2$ into stars, additional accretion is needed to build up $\rho_{H2}$, while HI is slowly being depleted. Conversely, at $z \lesssim 1.5$ the net inflow rate $\tilde{\psi}_{HI \to H2}$ is lower than $\psi_{stars}$. This is because the H$_2$ reservoir is decreasing with time in this redshift range, and therefore less H$_2$ needs to be replenished.

## 4.4. The Cosmic Future

Under the assumption of continuity, and that the physical process currently in play continue to dominate, we can use our empirical fitting functions (Section 3.3) to forecast the evolution of the baryon content associated with galaxies over the next few Gyr. This is shown in Figure 7 where we plot the same information as in Figure 3 but as a function of (linear) cosmic time. Assuming that our fits can be extrapolated to the future, the molecular mass density will decrease by about a factor of two over the next 5 Gyr, the HI mass density will remain approximately constant, and the stellar mass density will increase by about 10%. The star formation rate density will follow the decrease of H$_2$. Consequently, the total cold gas content in galaxies will be dominated by diffuse atomic gas even more than today. In this scenario, the ionized gas in the ICM/CGM will stay in this state and will not enter the main body of the galaxies. The inflow and infall rates (Equations (5) and (6)) will decrease correspondingly. Figure 7 shows that the universe has entered

“Cosmic Twilight,” during which the star formation activity in galaxies inexorably declines, as the gas inflow and accretion shuts down (see also Salcido et al. 2018). Over this same time period, the majority of stars with masses greater than the Sun will have exceeded their main-sequence lifetimes, leaving increasingly cooler, low-mass stars to illuminate the universe.

## 4.5. Theory Connection

Thus far, we have taken a strictly phenomenological approach to the trends observed in the data. We now discuss if cosmological simulations provide a sufficient amount of (dark and baryonic) matter to be accreted onto galaxy halos, to account for the observed net flows ($\tilde{\psi}_{HII \to HI}$ and $\tilde{\psi}_{HI \to H2}$). We also consider the potential role of preventive feedback mechanisms (such as virial shocks, AGN feedback, and cosmic expansion).

### 4.5.1. Accretion onto Dark Matter Halos

We estimate the amount of baryonic matter that is accreted onto galaxy halos using the results from cosmological simulations. More specifically, we estimate the matter (dark and baryonic combined) accretion rate onto halos $\dot{M}_{matter}(M_{vir}, z)$ as a function of halo virial mass and redshift using the fitting function presented in Rodríguez-Puebla et al. (2016, their Equation (11) adopting the dynamically averaged scenario). The authors obtained this fitting function by measuring the growth of halos in the Bolshoi-Planck and MultiDark-Planck $\Lambda$CDM cosmological simulations (Klypin et al. 2016). The cosmic (dark + baryonic) matter accretion rate $\psi_{matter}(z)$ is then obtained by taking the integral (over the virial masses considered) of the product between the matter accretion rate $\dot{M}_{matter}(M_{vir}, z)$ and the number density of halos with that mass $\Phi_{vir}(M_{vir}, z)$, such that

$$ \psi_{\text{matter}}(z) = \int_{M_{\text{vir,min}}}^{M_{\text{vir,max}}} \dot{M}_{\text{matter}}(M_{\text{vir}}, z) \times \Phi_{\text{vir}}(M_{\text{vir}}, z) dM_{\text{vir}}, \quad (7) $$

where the number density of halos as a function of virial mass and redshift is from Equation (23) in Rodríguez-Puebla et al. (2016). These accretion curves are shown in Figure 8 as dashed lines. The different lines show the accretion rates assuming different dark matter halo mass ranges, where the lowest mass considered here ($M_{vir} = 10^{10} M_{\odot}$) corresponds to the mass resolution in the simulations considered (corresponding to a stellar mass of a few times $10^7 M_{\odot}$). The resulting accretion rates are similar to matter accretion rates estimated in earlier works (e.g., Dekel et al. 2009).

As we are not primarily interested in the accretion of dark matter, but of the baryonic matter, we multiply the total matter accretion rate with the constant baryonic matter fraction to obtain the baryonic accretion rate onto halos $\dot{M}_{baryon}(M_{vir}, z)$. This assumes a perfect mixing between dark and baryonic matter in the IGM. The resulting baryonic accretion rates are shown as solid lines in Figure 8.

It is interesting to note that these accretion curves have a shape similar to our derived net infall/inflow rates ($\tilde{\psi}_{HI \to H2}$, $\tilde{\psi}_{HII \to HI}$): the accretion rates rise from high redshift to about $z \sim 2$ (depending on the virial masses considered). This increase in accretion to its peak value is dominated by gravitationally driven growth of the halo mass function. The subsequent decline toward $z = 0$ is due to the fact that the
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Same information as Figure 3, but with the following changes: (a) the lower abscissa shows cosmic time on a linear scale (redshift on the upper abscissa), (b) we extrapolate our fitting functions to the future (the present day is indicated by a vertical line, $z=0$), (c) we add units on the ordinate axis in $g \text{ cm}^{-3}$. As in all other plots we start plotting our functions at $z=4$. Under the assumption that our extrapolations are valid, the molecular gas density will decline by about a factor two over the next 5 Gyr, the stellar mass will increase by approximately 10%, and the inflow and accretion rates will decline correspondingly.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Comparison of the observed net accretion rates (orange/black curves) and predictions from theory (blue curves). The observed net infall and net inflow rates onto the central disk are the same as in Figure 6, but are shown here on a logarithmic ordinate axis. The predictions from theory, based on the Bolshoi-Planck and MultiDark-Planck ΛCDM cosmological simulations, of the accretion rate of the total (dark and baryonic) matter onto the dark matter halo are shown as dashed blue curves for different virial mass ranges. The solid blue curves show the accretion rates for baryonic matter only (see discussion in Section 4.5). The predicted baryonic accretion rates onto the galaxy halos are larger than the observationally required net infall rates onto the central disk, indicating that most of the accreted baryons do not end up in the centers of galaxies.</figcaption></figure>

universe expands and to the gradual decrease in the availability of accretable (dark) matter.<fnref n="37" />

### 4.5.2. Accretion Onto Central Disks

So far we have only considered the accretion of matter on a dark matter halo. We now compare these rates to the actual accretion to the central disk, and add our H II net infall and H I inflow rates to Figure 8 (same curves as in Figure 6, but on a logarithmic scale). A comparison to the total baryonic matter that is being accreted onto galaxies (solid blue lines) immediately implies that the material that is needed for the infall/inflow rates can be easily accounted for: of the total matter that is being accreted onto the dark matter halos of galaxies, only about 10%–30% is needed to explain the infall/ inflow rates that are inferred by the observations (Section 4.3). Consequently, the majority of the accreted baryons will not make it to the central galaxy “disks.”

An extensive body of literature has addressed the question of how the material that is accreted onto dark matter halos ends up in the centers of galaxies. In the standard picture, baryons from the IGM accrete onto dark matter halos, converting their gravitational energy into kinetic energy, which is subsequently shock-heated to the virial temperature of the halo. In addition to the formation of this hot halo, a large body of work suggests that dense filaments permeate the halos, leading to the formation of cold streams that feed the cold gas reservoir, and thus star formation, in the centers of galaxies. This process,

<footnote n="37">³⁷ As pointed out by Salcido et al. (2018), a universe without an accelerated (Λ-dominated) expansion does not significantly change the accretion rates, i.e., the *accelerated* expansion of the universe is not the reason for the observed decline in the accretion rates.</footnote>
referred to as “cold mode accretion,” occurs on timescales of the order of a free-fall or the dynamical crossing time of a spherical halo (~10⁸ yr, depending on mass), where the separation between the “cold” and “hot” phase is around 10⁵.⁵ K. The fraction of gas that is accreted in this cold phase depends on both the halo mass and the redshift, but it is the cold mode that appears to be the dominant accretion mechanism throughout all redshifts in most simulations for all but the most massive halos (Kereš et al. 2005; Dekel & Birnboim 2006; Dekel et al. 2009; van de Voort et al. 2011; Pan et al. 2019; Nelson et al. 2013). Once the gas is in a cool phase, it cools quickly to lower temperatures that are typical of the atomic/molecular interstellar medium, on timescales <10⁷ yr (e.g., Cornuault et al. 2018). Our observations do not allow us to distinguish between the different accretion mechanisms (“cold” versus “hot”).

We note that the decline (from the peak to $z=0$) in the baryonic accretion rate onto the halo (factor of a few) is smaller than that of our observationally derived net infall/inflow rates onto the disks (decline by almost an order of magnitude). This implies that additional mechanisms are suppressing the accretion of material, and these mechanisms become more dominant at $z < 2$. In a simple picture, the baryonic material that does not make it to the galaxy centers is heated by a number of processes, e.g., shocks, photoionization, through, e.g., stellar and AGN feedback. This hot material has very long cooling times. Assuming a typical temperature ~10⁶ K and a density ~10⁻⁵ cm⁻³, with substantial variation (e.g., Shull et al. 2012, 2017; Nicastro et al. 2018), the nominal bremsstrahlung cooling time is about 8.5 × 10¹¹ yr, or more than 60 times the Hubble time (Rosati et al. 2002).

We note that the cosmic density of AGN and star formation has decreased by about an order of magnitude since its peak, and will continue to do so in the future. Hence, feedback must play a less important role at late cosmic times. To explain the continued decline in the cosmic star formation rate at late times, we conjecture that only the densest gas in IGM filaments has been able to cool and stream into galaxy potential wells, and that these dense regions have been effectively “tapped out” over the eons. In this picture, most of the gas in the IGM that was predestined to fall into galaxies has done so already, and what is left will diffuse away with cosmic expansion.<fnref n="38" />

## 5. Concluding Remarks

We have used measurements of the cosmic molecular gas density to put new constraints on the baryon cycle and the gas accretion process for gas that is gravitationally bound to galaxies. We find that the cosmic H₂ density is less than or equal to the cosmic H I density over all times, briefly approaching equality at $z \sim 1.5$. Below a redshift of $z \sim 1.5$, the stellar mass density begins to dominate over all gas components (H₂ and H I), and completely dominates the baryon content within the main body of galaxies by $z=0$. The average cosmic gas depletion time, defined as the molecular gas density divided by the star formation rate density, is approximately constant above redshifts $z \gtrsim 2$, and then increases slightly from $\tau_{depl} \sim (4 \pm 2) \times 10^8$ yr at $z \sim 2$ to $\tau_{depl} = (7 \pm 3) \times 10^8$ yr at $z=0$. Significant accretion of gas onto galaxies is needed to

form the bulk of the stellar mass at $z < 1.5$: assuming that the maximum molecular gas density (seen at $z \sim 1.5$) will be fully transformed to stellar mass can only account for at most a quarter of the stellar mass seen at $z=0$.

The new H₂ constraints can be used to break up the gas accretion process onto galaxies in two steps. (i) First is the net inflow of atomic gas, and conversion to molecular gas, from the extended reservoirs to the centers of the dark matter halos (Equation (5)). (ii) Second is the net infall of mostly diffuse (ionized) gas to refuel the H I reservoirs (Equation (6)). We find that both flow processes decrease sharply at redshifts $z \lesssim 1.5$, following, to first order, the star formation rate density.

Zooming out, we can describe the gas cycle in galaxies as follows: an extended reservoir of atomic gas (H I) is formed by a (net) infall of gas from the IGM/CGM at a rate of $\bar{\psi}_{HII \to HII}$. This extended H I component is not immediately associated with the star formation process. Further (net) inflow from the H I reservoir at a rate of $\bar{\psi}_{HII \to H2}$ then leads to a molecular gas phase in the centers of the dark matter potentials. As the extended H I density remains approximately constant, these two net rates are similar. Stars are then formed out of the molecular gas phase, and the resulting star formation rate surface density in a galaxy is expected to be correlated with the molecular gas surface density (e.g., Bigiel et al. 2008; Leroy et al. 2013). The functional shape of the star formation rate density $\psi_{stars}$ is thus mostly driven by the availability of molecular gas, which in turn is defined by (net) infall rates of gas from larger distances. A comparison to numerical simulations shows that there is ample material that is being accreted onto dark matter halos. The decrease in gas accretion since $z \sim 1.5$ is then a result of the decreased growth of dark matter halos (partly due to the expansion of the universe), combined with the effects of feedback from stars and accreting black holes.

Lastly, by extrapolating our empirical fitting functions for the evolution of the stellar mass, H I, and H₂, we find that the molecular gas density will decrease by about a factor of two in the next 5 Gyr, the H I mass density will remain approximately constant, and the stellar mass density will increase by about 10%. The inflow and accretion rates will decrease correspondingly, and the cosmic star formation rate density will continue its steady descent to the infinitesimal.

We thank the referee for a very constructive report that helped to improve the paper. We thank Annalisa Pillepich and Andrea Ferrara for useful discussions. F.W. and M.N. acknowledge support from the ERC Advanced Grant 740246 (Cosmic_Gas). B.M. acknowledges support from the Collaborative Research Centre 956, sub-project A1, funded by the Deutsche Forschungsgemeinschaft (DFG)—project ID 184018867. T.D.-S. acknowledges support from the CAS-SACA and CONICYT fund CAS-CONICYT Call 2018. R.J.A. was supported by FONDECYT grant 1191124. D.R. acknowledges support from the National Science Foundation under grant numbers AST-1614213 and AST-1910107 and from the Alexander von Humboldt Foundation through a Humboldt Research Fellowship for Experienced Researchers. I.R.S. acknowledges support from STFC (ST/T000244/1). H.I. acknowledges support from JSPS KAKENHI grant No. JP19K23462. J.H. acknowledges support of the VIDI research program with project number 639.042.611, which is (partly) financed by the Netherlands Organisation for Scientific Research (NWO). D.O. is a recipient of an Australian Research

<footnote n="38">³⁸ This situation is similar to the conclusion of the pioneering work by Toomre & Toomre (1972) on the predestiny of galaxy mergers, in which they conclude: 'Hence one must presume that the partners in most cases were already bound to each other prior to their latest encounter.'</footnote>
Council Future Fellowship (FT190100083) funded by the Australian Government. This paper makes use of the following ALMA data: ADS/JAO.ALMA# 2017.1.00118.S, ADS/JAO.ALMA# 2015.1.01115.S. 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.

Facility: ALMA.

## Appendix A
### Background for Equations (3) and (4)

Here we provide some more background for the derivations of Equations (3) and (4) in the main body of the text. We consider the four main baryonic phases that are introduced in the schematic (Figure 1). These are measured by their respective space densities: $\rho_{\text{stars}}$ of the gas in disks and bulges, $\rho_{\text{HII}}$ of the ionized gas in the CGM and IGM, $\rho_{\text{HI}}$ of the atomic gas within galaxy disks and their environment, and $\rho_{\text{H2}}$ of the molecular gas in galaxy disks. We neglect all other minor mass components discussed in Section 3.1. We can then express the phase evolution of the universe in terms of 12 flow rates $\psi_{x \to y}(z)$ (with $x, y$ = [stars, HII, HI, H2], and $x \neq y$) that describe the following four phases $\rho_x$. These flow rates are a function of redshift $z$, but for simplicity we omit the (z) notation in the following.

$$ \dot{\rho}_{\text{stars}} = \frac{\psi_{\text{H2} \to \text{stars}} - \psi_{\text{stars} \to \text{H2}}}{\text{H2} \to \text{stars}} - \frac{\psi_{\text{stars} \to \text{HII}} + \psi_{\text{HII} \to \text{stars}}}{\text{HII} \to \text{stars}} \times \frac{-\psi_{\text{stars} \to \text{HI}} + \psi_{\text{HI} \to \text{stars}}}{\text{HI} \to \text{stars}}, \quad (\text{A1}) $$

$$ \dot{\rho}_{\text{HII}} = \frac{\psi_{\text{H2} \to \text{HII}} - \psi_{\text{HII} \to \text{H2}}}{\text{H2} \to \text{HII}} - \frac{\psi_{\text{HII} \to \text{HI}} + \psi_{\text{HI} \to \text{HII}}}{\text{HI} \to \text{HII}} \times \frac{-\psi_{\text{HII} \to \text{stars}} + \psi_{\text{stars} \to \text{HII}}}{\text{stars} \to \text{HII}}, \quad (\text{A2}) $$

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Diagram of the flows between the different main baryonic phases (stars, ionized gas H II, atomic gas H I, molecular gas H2) in the universe. The phase evolution can be expressed with 12 flow rates $\psi_{x \to y}$ that describe the 4 phases $\rho_x$ as indicated in the left panel (Equations (A1)–(A4)). Darker arrows indicate the main flows in the diagram. The right panel shows how the flows are simplified, by setting some flow terms to zero (see text) and introducing *net* flows that are marked with a tilde.</figcaption></figure>

$$ \dot{\rho}_{\text{HI}} = \frac{\psi_{\text{H2} \to \text{HI}} - \psi_{\text{HI} \to \text{H2}}}{\text{H2} \to \text{HI}} + \frac{\psi_{\text{HII} \to \text{HI}} - \psi_{\text{HI} \to \text{HII}}}{\text{HII} \to \text{HI}} \times \frac{+\psi_{\text{stars} \to \text{HI}} - \psi_{\text{HI} \to \text{stars}}}{\text{stars} \to \text{HI}}, \quad (\text{A3}) $$

$$ \dot{\rho}_{\text{H2}} = \frac{\psi_{\text{HII} \to \text{H2}} - \psi_{\text{H2} \to \text{HII}}}{\text{HII} \to \text{H2}} + \frac{\psi_{\text{stars} \to \text{H2}} - \psi_{\text{H2} \to \text{stars}}}{\text{stars} \to \text{H2}} \times \frac{+\psi_{\text{HI} \to \text{H2}} - \psi_{\text{H2} \to \text{HI}}}{\text{HI} \to \text{H2}}. \quad (\text{A4}) $$

These flow rates are also visualized in Figure 9, and the “+” and “-” signs in the equations above denote “gains” and “losses,” as indicated by the arrows in that figure. We have measurements for three of these quantities, i.e., $\dot{\rho}_{\text{stars}}$, $\dot{\rho}_{\text{HI}}$, and $\dot{\rho}_{\text{H2}}$, but many more unknowns, leaving the problem unsolvable. We can however simplify the above equations given our assumptions discussed in Section 2 and the corresponding schematic (Figure 1).

$\dot{\rho}_{\text{stars}}$: in Equation (A1) the first two terms correspond to the mass flows between the stars and the H2, the second two terms to the flows between the stars and the H II, and the last two terms the flows between the stars and the H I. We here assume that stars do not directly form out of H I or H II, and thus set these terms ($\psi_{\text{HI} \to \text{stars}}$, $\psi_{\text{HII} \to \text{stars}}$) to zero. Likewise, we assume that stars do not produce atomic hydrogen nor molecular gas, and we thus set these two terms ($\psi_{\text{stars} \to \text{HI}}$, $\psi_{\text{stars} \to \text{H2}}$) to zero as well. Equation (A1) therefore simplifies to:

$$ \dot{\rho}_{\text{stars}} = \psi_{\text{H2} \to \text{stars}} - \psi_{\text{stars} \to \text{HII}}. \quad (\text{A5}) $$

$\dot{\rho}_{\text{HII}}$: in Equation (A2) the first two terms correspond to the mass flows between the H2 and the ionized gas (H II), the second two terms to the flows between the atomic and ionized gas, and the last two terms the flows between the stars and the ionized gas. We assume that ionized gas cannot directly form molecular gas, and that ionized gas cannot directly form stars and set $\psi_{\text{HII} \to \text{stars}}$ and $\psi_{\text{HII} \to \text{H2}}$ to zero. The other flows are in principle plausible, and Equation (A2) thus becomes

$$ \dot{\rho}_{\text{HII}} = \psi_{\text{H2} \to \text{HII}} - \psi_{\text{HII} \to \text{HI}} + \psi_{\text{HI} \to \text{HII}} + \psi_{\text{stars} \to \text{HII}}, \quad (\text{A6}) $$

$$ \dot{\rho}_{\text{HII}} = \psi_{\text{H2} \to \text{HII}} - \tilde{\psi}_{\text{HII} \to \text{HI}} + \psi_{\text{stars} \to \text{HII}}, \quad (\text{A7}) $$

where in a second step we introduced $\tilde{\psi}_{\text{HII} \to \text{HI}}$ as being the *net* flow between the ionized gas phase and the atomic gas. As we

<figure><img src="image_13.png" /></figure>
assume that $|\psi_{HII \to HI}| > |\psi_{HI \to HII}|$ we assign a minus sign to the net flow $\tilde{\psi}_{HII \to HI}$ in Equation (A7).

$\dot{\rho}_{HI}$: in Equation (A3) the first two terms correspond to the mass flows between the $H_2$ and the atomic gas, the second two terms to the flows between the atomic and ionized gas, and the last two terms to the flows between the stars and the atomic gas. As before, we assume that no stars can form out of H I, and that stars cannot form observable (i.e., long-lived) H I, i.e., both $\psi_{HI \to stars}$ and $\psi_{stars \to HI}$ are set to zero. Equation (A3) then yields

$$ \dot{\rho}_{HI} = \psi_{H2 \to HI} - \psi_{HI \to H2} + \psi_{HII \to HI} - \psi_{HI \to HII}, \quad (A8) $$

$$ \dot{\rho}_{HI} = -\tilde{\psi}_{HI \to H2} + \tilde{\psi}_{HII \to HI}, \quad (A9) $$

where in a second step we again introduce an additional net flow between the atomic gas and the $H_2$ ($\tilde{\psi}_{HI \to H2}$). Its sign is determined as we assume $|\psi_{H2 \to HI}| < |\psi_{HI \to H2}|$. This equation is identical to Equation (4) in the main body of this manuscript.

$\dot{\rho}_{H2}$: in Equation (A4) the first two terms correspond to the mass flows between the $H_2$ and the ionized gas, the second two terms to the flows between the $H_2$ and the stars, and the last two terms to the flows between the $H_2$ and the atomic gas. As before, we set $\psi_{HII \to H2}$ and $\psi_{stars \to H2}$ to zero. This yields

$$ \dot{\rho}_{H2} = -\psi_{H2 \to HII} - \psi_{H2 \to stars} + \psi_{HI \to H2} - \psi_{H2 \to HI}, \quad (A10) $$

$$ \dot{\rho}_{H2} = -\psi_{H2 \to HII} - \psi_{H2 \to stars} + \tilde{\psi}_{HI \to H2}, \quad (A11) $$

where we again use the previously defined net flow between the atomic gas and the $H_2$, $\tilde{\psi}_{HI \to H2}$, in the second step.

We can now further simplify Equation (A11) by setting $\psi_{H2 \to stars}$ to the star formation rate density, i.e., $\psi_{H2 \to stars} = \psi_{stars}$, i.e., the fundamental process that forms stars out of molecular gas. We can also set the feedback rate (molecular gas that will be ionized), $\psi_{H2 \to HII}$, to be proportional to the star formation rate density, i.e., $\psi_{H2 \to HII} = \xi \psi_{stars}$. This yields Equation (3) in the main text.

As a side note, assuming that the return rate of the stars to the ionized medium ($\psi_{stars \to HII}$) is proportional to the star formation rate density $\psi_{stars}$ with a proportionality factor $\tilde{R}$, i.e., $\psi_{stars \to HII} = \tilde{R} \psi_{stars}$, we can rewrite Equation (A5) that expresses the change in the stellar mass density as follows

$$ \dot{\rho}_{stars} = \psi_{stars} - \tilde{R}\psi_{stars}, \quad (A12) $$

$$ = (1 - \tilde{R}) \psi_{stars}. \quad (A13) $$

Note that the factor $\tilde{R}$ would be equal to the classical return factor $R$ (see Section 3.2.1) if we assume that all loss of stellar mass (stellar winds, SN explosions) would end up in the ionized phase of the interstellar medium.

## Appendix B
### Observational Data for the Cosmic Density Measurements

In this appendix, we give the observational measurements for both the cosmic H I mass density (Table 2) and the cosmic $H_2$ mass density (Table 3) from the literature. The observational data used to fit the cosmic stellar mass density and the cosmic star formation rate density are taken from the compilation in Madau & Dickinson (2014).

<table><caption>**Table 2**
Measurements of the Cosmic H I Mass Density</caption><thead><tr><th>Redshift</th><th>ρ<sub>HI</sub><br>(10<sup>8</sup>M<sub>c</sub> Mpc<sup>-3</sup>)</th><th>Method</th><th>Reference</th></tr></thead><tbody><tr><td>0.0</td><td>0.60 ± 0.10</td><td>21 cm</td><td>Zwaan et al. (2005)</td></tr><tr><td>0.0</td><td>0.51 ± 0.09</td><td>21 cm</td><td>Jones et al. (2018)</td></tr><tr><td>0.0</td><td>1.03 ± 0.17</td><td>21 cm</td><td>Braun (2012)</td></tr><tr><td>0.028</td><td>0.65<sup>+0.13</sup><sub>-0.06</sub></td><td>21 cm stacked</td><td>Delhaize et al. (2013)</td></tr><tr><td>0.096</td><td>0.73<sup>+0.13</sup><sub>-0.10</sub></td><td>21 cm stacked</td><td>Delhaize et al. (2013)</td></tr><tr><td>0.06</td><td>0.37 ± 0.11</td><td>21 cm stacked</td><td>Hoppmann et al. (2015)</td></tr><tr><td>0.1</td><td>0.53 ± 0.08</td><td>21 cm stacked</td><td>Rhee et al. (2013)</td></tr><tr><td>0.2</td><td>0.55 ± 0.14</td><td>21 cm stacked</td><td>Rhee et al. (2013)</td></tr><tr><td>0.24</td><td>1.11 ± 0.49</td><td>21 cm stacked</td><td>Lah et al. (2007)</td></tr><tr><td>0.2–0.4</td><td>0.62 ± 0.09</td><td>21 cm stacked</td><td>Bera et al. (2019)</td></tr><tr><td>1.265</td><td>&lt;0.342</td><td>21 cm stacked</td><td>Kanekar et al. (2016)</td></tr><tr><td>0.2–1.8</td><td>0.99<sup>+0.55</sup><sub>-0.36</sub></td><td>21 cm cross-<br>correlation</td><td>Masui et al. (2013)</td></tr><tr><td>0.11–0.61</td><td>1.11 ± 0.37</td><td>Mg II-selection</td><td>Rao et al. (2017)</td></tr><tr><td>0.61–0.89</td><td>0.96 ± 0.23</td><td>Mg II-selection</td><td>Rao et al. (2017)</td></tr><tr><td>0.89–1.65</td><td>1.08 ± 0.42</td><td>Mg II-selection</td><td>Rao et al. (2017)</td></tr><tr><td>0.01–0.48</td><td>0.40<sup>+0.31</sup><sub>-0.19</sub></td><td>DLA</td><td>Shull et al. (2017)</td></tr><tr><td>0.01–1.6</td><td>0.33<sup>+0.26</sup><sub>-0.16</sub></td><td>DLA</td><td>Neeleman et al. (2016)</td></tr><tr><td>1.55–2.0</td><td>1.01<sup>+0.34</sup><sub>-0.29</sub></td><td>DLA</td><td>Péroux et al. (2003)</td></tr><tr><td>2.0–2.7</td><td>1.45<sup>+0.34</sup><sub>-0.30</sub></td><td>DLA</td><td>Péroux et al. (2003)</td></tr><tr><td>2.7–3.5</td><td>1.33<sup>+0.35</sup><sub>-0.30</sub></td><td>DLA</td><td>Péroux et al. (2003)</td></tr><tr><td>3.5–4.85</td><td>0.82<sup>+0.30</sup><sub>-0.27</sub></td><td>DLA</td><td>Péroux et al. (2003)</td></tr><tr><td>1.55–2.73</td><td>2.13<sup>+1.28</sup><sub>-0.86</sub></td><td>DLA</td><td>Sánchez-Ramírez et al.<br>(2016)</td></tr><tr><td>2.73–3.21</td><td>0.79<sup>+0.63</sup><sub>-0.38</sub></td><td>DLA</td><td>Sánchez-Ramírez et al.<br>(2016)</td></tr><tr><td>3.21–4.5</td><td>1.78<sup>+0.49</sup><sub>-0.42</sub></td><td>DLA</td><td>Sánchez-Ramírez et al.<br>(2016)</td></tr><tr><td>2.0–2.3</td><td>1.42 ± 0.07</td><td>DLA</td><td>Noterdaeme et al. (2012)</td></tr><tr><td>2.3–2.6</td><td>1.25 ± 0.06</td><td>DLA</td><td>Noterdaeme et al. (2012)</td></tr><tr><td>2.6–2.9</td><td>1.50 ± 0.07</td><td>DLA</td><td>Noterdaeme et al. (2012)</td></tr><tr><td>2.9–3.2</td><td>1.58 ± 0.11</td><td>DLA</td><td>Noterdaeme et al. (2012)</td></tr><tr><td>3.2–3.5</td><td>1.83 ± 0.19</td><td>DLA</td><td>Noterdaeme et al. (2012)</td></tr><tr><td>2.2–2.4</td><td>0.74 ± 0.13</td><td>DLA</td><td>Prochaska &<br>Wolfe (2009)</td></tr><tr><td>2.4–2.7</td><td>1.00 ± 0.11</td><td>DLA</td><td>Prochaska &<br>Wolfe (2009)</td></tr><tr><td>2.7–3.0</td><td>1.00 ± 0.10</td><td>DLA</td><td>Prochaska &<br>Wolfe (2009)</td></tr><tr><td>3.0–3.5</td><td>1.40 ± 0.12</td><td>DLA</td><td>Prochaska &<br>Wolfe (2009)</td></tr><tr><td>3.5–4.0</td><td>1.62 ± 0.28</td><td>DLA</td><td>Prochaska &<br>Wolfe (2009)</td></tr><tr><td>4.0–5.5</td><td>1.58<sup>+0.34</sup><sub>-0.30</sub></td><td>DLA</td><td>Prochaska &<br>Wolfe (2009)</td></tr><tr><td>2.55–3.4</td><td>1.48<sup>+0.66</sup><sub>-0.50</sub></td><td>DLA</td><td>Guimarães et al. (2009)</td></tr><tr><td>3.4–3.83</td><td>1.31<sup>+0.51</sup><sub>-0.45</sub></td><td>DLA</td><td>Guimarães et al. (2009)</td></tr><tr><td>3.83–5.03</td><td>1.16<sup>+0.52</sup><sub>-0.42</sub></td><td>DLA</td><td>Guimarães et al. (2009)</td></tr><tr><td>3.56–4.45</td><td>1.87<sup>+0.42</sup><sub>-0.40</sub></td><td>DLA</td><td>Crighton et al. (2015)</td></tr><tr><td>4.45–5.31</td><td>1.56<sup>+0.31</sup><sub>-0.27</sub></td><td>DLA</td><td>Crighton et al. (2015)</td></tr></tbody></table>

**Note.** Values have been converted, where necessary, to the adopted cosmology and have been corrected for the contribution of helium. For the DLA measurements, no correction for lower H I column density systems has been applied. The table has been subdivided into categories based on the observational method used to estimate $\rho_{HI}$.
<table><caption>**Table 3**
Measurements of the Cosmic H₂ Mass Density</caption><thead><tr><th>Redshift</th><th>ρ<sub>H2</sub><br>(10<sup>8</sup>M<sub>☉</sub> Mpc<sup>-3</sup>)</th><th>Method</th><th>Reference</th></tr></thead><tbody><tr><td>0.01–0.05</td><td>0.104<sup>+0.009</sup><sub>-0.009</sub></td><td>CO/<br>xCOLD GASS</td><td>Fletcher et al. (2020)</td></tr><tr><td>0.00–0.37<sup>a</sup></td><td>0.009<sup>+0.019</sup><sub>-0.007</sub></td><td>CO/ASPECS</td><td>Decarli et al. (2019)</td></tr><tr><td>0.27–0.63</td><td>0.11<sup>+0.10</sup><sub>-0.05</sub></td><td>CO/ASPECS</td><td>Decarli et al. (2020)</td></tr><tr><td>0.69–1.17</td><td>0.46<sup>+0.27</sup><sub>-0.18</sub></td><td>CO/ASPECS</td><td>Decarli et al. (2020)</td></tr><tr><td>1.01–1.74</td><td>0.55<sup>+0.20</sup><sub>-0.15</sub></td><td>CO/ASPECS</td><td>Decarli et al. (2019)</td></tr><tr><td>2.01–3.11</td><td>0.29<sup>+0.15</sup><sub>-0.11</sub></td><td>CO/ASPECS</td><td>Decarli et al. (2019)</td></tr><tr><td>3.01–4.47</td><td>0.24<sup>+0.09</sup><sub>-0.07</sub></td><td>CO/ASPECS</td><td>Decarli et al. (2019)</td></tr><tr><td>0.48–1.48</td><td>0.41<sup>+0.11</sup><sub>-0.12</sub></td><td>CO/PHIBBS2</td><td>Lenkić et al. (2020)</td></tr><tr><td>1.01–2.01</td><td>0.48<sup>+0.12</sup><sub>-0.11</sub></td><td>CO/PHIBBS2</td><td>Lenkić et al. (2020)</td></tr><tr><td>2.01–3.01</td><td>0.37<sup>+0.10</sup><sub>-0.11</sub></td><td>CO/PHIBBS2</td><td>Lenkić et al. (2020)</td></tr><tr><td>3.14–4.14</td><td>0.15<sup>+0.07</sup><sub>-0.07</sub></td><td>CO/PHIBBS2</td><td>Lenkić et al. (2020)</td></tr><tr><td>4.25–5.25<sup>b</sup></td><td>0.10<sup>+0.07</sup><sub>-0.06</sub></td><td>CO/PHIBBS2</td><td>Lenkić et al. (2020)</td></tr><tr><td>1.95–2.85</td><td>0.27<sup>+0.16</sup><sub>-0.11</sub></td><td>CO/COLDz</td><td>Riechers et al.<br>(2019)</td></tr><tr><td>4.9–6.7<sup>b</sup></td><td>0.047<sup>+0.034</sup><sub>-0.023</sub></td><td>CO/COLDz</td><td>Riechers et al.<br>(2019)</td></tr><tr><td>1.95–2.85</td><td>0.28<sup>+0.18</sup><sub>-0.12</sub></td><td>CO/VLASPECS</td><td>Riechers et al.<br>(2020)</td></tr><tr><td>0.3–0.6</td><td>0.20<sup>+0.16</sup><sub>-0.14</sub></td><td>Dust continuum</td><td>Magnelli et al.<br>(2020)</td></tr><tr><td>0.6–1.0</td><td>0.30 ± 0.18</td><td>Dust continuum</td><td>Magnelli et al.<br>(2020)</td></tr><tr><td>1.0–1.6</td><td>0.51 ± 0.11</td><td>Dust continuum</td><td>Magnelli et al.<br>(2020)</td></tr><tr><td>1.6–2.3</td><td>0.33 ± 0.07</td><td>Dust continuum</td><td>Magnelli et al.<br>(2020)</td></tr><tr><td>2.3–3.2</td><td>0.42 ± 0.09</td><td>Dust continuum</td><td>Magnelli et al.<br>(2020)</td></tr><tr><td>0.2–0.5</td><td>0.43<sup>+0.06</sup><sub>-0.05</sub></td><td>Dust continuum</td><td>Scoville et al. (2017)</td></tr><tr><td>0.5–0.8</td><td>0.47<sup>+0.08</sup><sub>-0.07</sub></td><td>Dust continuum</td><td>Scoville et al. (2017)</td></tr><tr><td>0.8–1.1</td><td>0.59<sup>+0.10</sup><sub>-0.09</sub></td><td>Dust continuum</td><td>Scoville et al. (2017)</td></tr><tr><td>1.1–1.5</td><td>0.84<sup>+0.14</sup><sub>-0.13</sub></td><td>Dust continuum</td><td>Scoville et al. (2017)</td></tr><tr><td>1.5–2.0</td><td>0.90<sup>+0.16</sup><sub>-0.13</sub></td><td>Dust continuum</td><td>Scoville et al. (2017)</td></tr><tr><td>2.0–2.5</td><td>0.74<sup>+0.12</sup><sub>-0.11</sub></td><td>Dust continuum</td><td>Scoville et al. (2017)</td></tr><tr><td>2.5–3.0</td><td>0.63<sup>+0.11</sup><sub>-0.09</sub></td><td>Dust continuum</td><td>Scoville et al. (2017)</td></tr><tr><td>3.0–4.0</td><td>0.36<sup>+0.06</sup><sub>-0.05</sub></td><td>Dust continuum</td><td>Scoville et al. (2017)</td></tr><tr><td>0.2–0.6</td><td>0.23<sup>+0.06</sup><sub>-0.05</sub></td><td>Dust continuum</td><td>Berta et al. (2013)</td></tr><tr><td>0.7–1.0</td><td>0.60<sup>+0.18</sup><sub>-0.17</sub></td><td>Dust continuum</td><td>Berta et al. (2013)</td></tr><tr><td>1.0–2.0</td><td>0.48<sup>+0.09</sup><sub>-0.07</sub></td><td>Dust continuum</td><td>Berta et al. (2013)</td></tr></tbody></table>

**Notes.** The table has been subdivided into categories based on the observational method used to estimate ρ<sub>H2</sub>, either through measurement of a CO line (from a blind survey or from targeted observations of individual galaxies), or through the measurement of the dust continuum.

<sup>a</sup> The lowest-redshift bin from the ASPECS survey is not used in the fit due to the small cosmic volume probed (Decarli et al. 2020).

<sup>b</sup> As we only consider redshifts z < 4 in this study, these measurements are not used in our fit.

## ORCID iDs

Fabian Walter https://orcid.org/0000-0003-4793-7880

Chris Carilli https://orcid.org/0000-0001-6647-3861

Marcel Neeleman https://orcid.org/0000-0002-9838-8191

Roberto Decarli https://orcid.org/0000-0002-2662-8803

Gergö Popping https://orcid.org/0000-0003-1151-4659

Manuel Aravena https://orcid.org/0000-0002-6290-3198

Frank Bertoldi https://orcid.org/0000-0002-1707-1775

Leindert Boogaard https://orcid.org/0000-0002-3952-8588

Pierre Cox https://orcid.org/0000-0003-2027-8221

Benjamin Magnelli https://orcid.org/0000-0002-6777-6490

Dominik Riechers https://orcid.org/0000-0001-9585-1462

Hans-Walter Rix https://orcid.org/0000-0003-4996-9069

Ian Smail https://orcid.org/0000-0003-3037-257X

Axel Weiss https://orcid.org/0000-0003-4678-3939

Roberto J. Assef https://orcid.org/0000-0002-9508-3667

Franz Bauer https://orcid.org/0000-0002-8686-8737

Rychard Bouwens https://orcid.org/0000-0002-4989-2471

Thierry Contini https://orcid.org/0000-0003-0275-938X

Paulo C. Cortes https://orcid.org/0000-0002-3583-780X

Emanuele Daddi https://orcid.org/0000-0002-3331-9590

Tanio Diaz-Santos https://orcid.org/0000-0003-0699-6083

Jorge González-López https://orcid.org/0000-0003-3926-1411

Joseph Hennawi https://orcid.org/0000-0002-7054-4332

Jacqueline A. Hodge https://orcid.org/0000-0001-6586-8845

Hanae Inami https://orcid.org/0000-0003-4268-0393

Rob Ivison https://orcid.org/0000-0001-5118-1313

Pascal Oesch https://orcid.org/0000-0001-5851-6649

Mark Sargent https://orcid.org/0000-0003-1033-9684

## References

Anderson, M. E., & Bregman, J. N. 2010, ApJ, 714, 320

Aravena, M., Boogaard, L., Gónzalez-López, J., et al. 2020, ApJ, 901, 79

Aravena, M., Decarli, R., Gónzalez-López, J., et al. 2019, ApJ, 882, 136

Bastian, N., Covey, K. R., & Meyer, M. R. 2010, ARA&A, 48, 339

Bauermeister, A., Blitz, L., & Ma, C.-P. 2010, ApJ, 717, 323

Behroozi, P., Wechsler, R. H., Hearin, A. P., & Conroy, C. 2019, MNRAS, 488, 3143

Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57

Bera, A., Kanekar, N., Chengalur, J. N., & Bagla, J. S. 2019, ApJL, 882, L7

Berg, T. A. M., Ellison, S. L., Sánchez-Ramírez, R., et al. 2019, MNRAS, 488, 4356

Berta, S., Lutz, D., Nordon, R., et al. 2013, A&A, 555, L8

Béthermin, M., Wang, L., Doré, O., et al. 2013, A&A, 557, A66

Bigiel, F., Leroy, A., Walter, F., et al. 2008, AJ, 136, 2846

Bird, S., Vogelsberger, M., Haehnelt, M., et al. 2014, MNRAS, 445, 2313

Blumenthal, G. R., Faber, S. M., Primack, J. R., & Rees, M. J. 1984, Natur, 311, 517

Boogaard, L., van der Werf, P., Weiss, A., et al. 2020, ApJ, 902, 109

Boogaard, L. A., Decarli, R., González-López, J., et al. 2019, ApJ, 882, 140

Boselli, A., Cortese, L., Boquien, M., et al. 2014, A&A, 564, A66

Bothwell, M. S., Smail, I., Chapman, S. C., et al. 2013, MNRAS, 429, 3047

Bouché, N., Dekel, A., Genzel, R., et al. 2010, ApJ, 718, 1001

Bouwens, R., González-López, J., Aravena, M., et al. 2020, ApJ, 902, 112

Bouwens, R. J., Aravena, M., Decarli, R., et al. 2016, ApJ, 833, 72

Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2012a, ApJ, 754, 83

Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2012b, ApJL, 752, L5

Braun, R. 2012, ApJ, 749, 87

Caputi, K. I., Cirasuolo, M., Dunlop, J. S., et al. 2011, MNRAS, 413, 162

Carilli, C. L., & Walter, F. 2013, ARA&A, 51, 105

Clark, P. C., Glover, S. C. O., Klessen, R. S., & Bonnell, I. A. 2012, MNRAS, 424, 2599

Conselice, C. J., Mortlock, A., Bluck, A. F. L., Grützbauch, R., & Duncan, K. 2013, MNRAS, 430, 1051

Cornuault, N., Lehnert, M. D., Boulanger, F., & Guillard, P. 2018, A&A, 610, A75

Crighton, N. H. M., Murphy, M. T., Prochaska, J. X., et al. 2015, MNRAS, 452, 217

Croton, D. J., Springel, V., White, S. D. M., et al. 2006, MNRAS, 365, 11

Cucciati, O., Tresse, L., Ilbert, O., et al. 2012, A&A, 539, A31

Daddi, E., Bournaud, F., Walter, F., et al. 2010, ApJ, 713, 686

Davé, R., Finlator, K., & Oppenheimer, B. D. 2012, MNRAS, 421, 98

Davidzon, I., Ilbert, O., Laigle, C., et al. 2017, A&A, 605, A70

Decarli, R., Aravena, M., Boogaard, L., et al. 2020, ApJ, 902, 110

Decarli, R., Walter, F., Aravena, M., et al. 2016, ApJ, 833, 69
Decarli, R., Walter, F., Carilli, C., et al. 2014, ApJ, 782, 78

Decarli, R., Walter, F., Gómez-López, J., et al. 2019, ApJ, 882, 138

Dekel, A., & Birnboim, Y. 2006, MNRAS, 368, 2

Dekel, A., Birnboim, Y., Engel, G., et al. 2009, Natur, 457, 451

Dekel, A., & Silk, J. 1986, ApJ, 303, 39

Delhaize, J., Meyer, M. J., Staveley-Smith, L., & Boyle, B. J. 2013, MNRAS, 433, 1398

Driver, S. P., Andrews, S. K., da Cunha, E., et al. 2018, MNRAS, 475, 2891

Dudzevičiūtė, U., Smail, I., Swinbank, A. M., et al. 2020, MNRAS, 494, 3828

Elbaz, D., Leiton, R., Nagar, N., et al. 2018, A&A, 616, A110

Ellison, S. L., Yan, L., Hook, I. M., et al. 2001, A&A, 379, 393

Fletcher, T. J., Saintonge, A., Soares, P. S., & Pontzen, A. 2020, arXiv:2002.04959

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, M. J. 2013, PASP, 125, 306

Förster Schreiber, N. M., Genzel, R., Bouché, N., et al. 2009, ApJ, 706, 1364

Fujimoto, S., Ouchi, M., Shibuya, T., & Nagai, H. 2017, ApJ, 850, 83

Fukui, Y., Kawamura, A., Wong, T., et al. 2009, ApJ, 705, 144

Genzel, R., Tacconi, L. J., Gracia-Carpio, J., et al. 2010, MNRAS, 407, 2091

Glover, S. C. O., & Mac Low, M. M. 2011, MNRAS, 412, 337

Gruppioni, C., Pozzi, F., Rodighiero, G., et al. 2013, MNRAS, 432, 23

Guimarães, R., Petitjean, P., de Carvalho, R. R., et al. 2009, A&A, 508, 133

Hodge, J. A., & da Cunha, E. 2020, arXiv:2004.00934

Hopkins, P. F., Wetzel, A., Kereš, D., et al. 2018, MNRAS, 480, 800

Hoppmann, L., Staveley-Smith, L., Freudling, W., et al. 2015, MNRAS, 452, 3726

Ilbert, O., McCracken, H. J., Le Fèvre, O., et al. 2013, A&A, 556, A55

Inami, H., Decarli, R., Walter, F., et al. 2020, ApJ, 902, 113

Jiménez-Andrade, E. F., Magnelli, B., Karim, A., et al. 2019, A&A, 625, A114

Jones, M. G., Haynes, M. P., Giovanelli, R., & Moorman, C. 2018, MNRAS, 477, 2

Jorgenson, R. A., Murphy, M. T., Thompson, R., & Carswell, R. F. 2014, MNRAS, 443, 2783

Jorgenson, R. A., Wolfe, A. M., Prochaska, J. X., et al. 2006, ApJ, 646, 730

Kanekar, N., Sethi, S., & Dwarakanath, K. S. 2016, ApJ, 818, L28

Kawamura, A., Mizuno, Y., Minamidani, T., et al. 2009, ApJS, 184, 1

Kennicutt, R. C., & Evans, N. J. 2012, ARA&A, 50, 531

Kereš, D., Katz, N., Weinberg, D. H., & Davé, R. 2005, MNRAS, 363, 2

Keres, D., Yun, M. S., & Young, J. S. 2003, ApJ, 582, 659

Khusanova, Y., Le Fèvre, O., Cassata, P., et al. 2020, A&A, 634, A97

Glitsch, A., Péroux, C., Zwaan, M. A., et al. 2019, MNRAS, 490, 1220

Klypin, A., Yepes, G., Gottlöber, S., Prada, F., & Heß, S. 2016, MNRAS, 457, 4340

Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511

Krogager, J.-K., Fynbo, J. P. U., Møller, P., et al. 2019, MNRAS, 486, 4377

Krogager, J. K., Møller, P., Fynbo, J. P. U., & Noterdaeme, P. 2017, MNRAS, 469, 2959

Lah, P., Chengalur, J. N., Briggs, F. H., et al. 2007, MNRAS, 376, 1357

Law, D. R., Steidel, C. C., Shapley, A. E., et al. 2012, ApJ, 759, 29

Leja, J., Speagle, J. S., Johnson, B. D., et al. 2020, ApJ, 893, 111

Lenkić, L., Bolatto, A. D., Förster Schreiber, N. M., et al. 2020, AJ, 159, 190

Leroy, A. K., Walter, F., Bigiel, F., et al. 2009, AJ, 137, 4670

Leroy, A. K., Walter, F., Sandstrom, A., et al. 2013, AJ, 146, 19

Lilly, S. J., Carollo, C. M., Pipino, A., Renzini, A., & Peng, Y. 2013, ApJ, 772, 119

Lilly, S. J., Le Fèvre, O., Hammer, F., & Crampton, D. 1996, ApJL, 460, L1

Liu, D., Schinnerer, E., Groves, B., et al. 2019, ApJ, 887, 235

Loeb, A. 2007, PhRvL, 99, 041103

Macquart, J. P., Prochaska, J. X., McQuinn, M., et al. 2020, Natur, 581, 391

Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415

Madau, P., Ferguson, H. C., Dickinson, M. E., et al. 1996, MNRAS, 283, 1388

Magnelli, B., Boogaard, L., Decarli, R., et al. 2020, ApJ, 892, 66

Magnelli, B., Elbaz, D., Chary, R. R., et al. 2011, A&A, 528, A35

Magnelli, B., Popesso, P., Berta, S., et al. 2013, A&A, 553, A132

Marchesini, D., van Dokkum, P. G., Förster Schreiber, N. M., et al. 2009, ApJ, 701, 1765

Masui, K. W., Switzer, E. R., Banavar, N., et al. 2013, ApJL, 763, L20

McQuinn, M. 2014, ApJL, 780, L33

Meidt, S. E., Hughes, A., Dobbs, C. L., et al. 2015, ApJ, 806, 72

Muzahid, S., Srianand, R., & Charlton, J. 2015, MNRAS, 448, 2840

Muzzin, A., Marchesini, D., Stefanon, M., et al. 2013, ApJ, 777, 18

Neeleman, M., Kanekar, N., Prochaska, J. X., et al. 2017, Sci, 355, 1285

Neeleman, M., Kanekar, N., Prochaska, J. X., Rafelski, M. A., & Carilli, C. L. 2019, ApJL, 870, L19

Neeleman, M., Prochaska, J. X., Ribaudo, J., et al. 2016, ApJ, 818, 113

Nelson, D., Vogelsberger, M., Genel, S., et al. 2013, MNRAS, 429, 3353

Nicalastro, F., Kaastra, J., Krongold, Y., et al. 2018, Natur, 558, 406

Noterdaeme, P., Ledoux, C., Petitjean, P., & Srianand, R. 2008, A&A, 481, 327

Noterdaeme, P., Petitjean, P., Carithers, W. C., et al. 2012, A&A, 547, L1

Novak, M., Smolčić, V., Delhaize, J., et al. 2017, A&A, 602, A5

O'Meara, J. M., Prochaska, J. X., Burles, S., et al. 2007, ApJ, 656, 666

Pan, Z., Peng, Y., Zheng, X., Wang, J., & Kong, X. 2019, ApJ, 885, 14

Pavesi, R., Sharon, C. E., Riechers, D. A., et al. 2018, ApJ, 864, 49

Peng, Y.-j., & Maiolino, R. 2014, MNRAS, 443, 3643

Pérez-González, P. G., Rieke, G. H., Villar, V., et al. 2008, ApJ, 675, 234

Pérez-Rafols, I., Font-Ribera, A., Miralda-Escudé, J., et al. 2018, MNRAS, 473, 3019

Peroux, C., & Howk, C. 2020, ARA&A, 58, 021820

Péroux, C., McMahon, R. G., Sorrie-Lombardi, L. J., & Irwin, M. J. 2003, MNRAS, 346, 1103

Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A13

Prochaska, J. X., & Neeleman, M. 2018, MNRAS, 474, 318

Prochaska, J. X., & Wolfe, A. M. 2009, ApJ, 696, 1543

Putman, M. E., Peek, J. E. G., & Joung, M. R. 2012, ARA&A, 50, 491

Rahmati, A., & Schaye, J. 2014, MNRAS, 438, 529

Rao, S. M., Turnshek, D. A., Sardane, G. M., & Monier, E. M. 2017, MNRAS, 471, 3428

Rathaus, B., & Sternberg, A. 2016, MNRAS, 458, 3168

Rhee, J., Zwaan, M. A., Briggs, F. H., et al. 2013, MNRAS, 435, 2693

Riechers, D. A., Boogaard, L. A., Decarli, R., et al. 2020, ApJ, 896, 21

Riechers, D. A., Pavesi, R., Sharon, C. E., et al. 2019, ApJ, 872, 7

Rodríguez-Puebla, A., Behroozi, P., Primack, J., et al. 2016, MNRAS, 462, 893

Rosati, P., Borgani, S., & Norman, C. 2002, ARA&A, 40, 539

Saintonge, A., Catinella, B., Tacconi, L. J., et al. 2017, ApJS, 233, 22

Salcido, J., Bower, R. G., Barnes, L. A., et al. 2018, MNRAS, 477, 3744

Salmi, F., Daddi, E., Elbaz, D., et al. 2012, ApJL, 754, L14

Sánchez-Ramírez, R., Ellison, S. L., Prochaska, J. X., et al. 2016, MNRAS, 456, 4488

Sandstrom, K. M., Leroy, A. K., Walter, F., et al. 2013, ApJ, 777, 5

Schinnerer, E., Hughes, A., Leroy, A., et al. 2019, ApJ, 887, 49

Schroetter, I., Bouché, N. F., Zabl, J., et al. 2019, MNRAS, 490, 4368

Schruba, A., Leroy, A. K., Walter, F., et al. 2011, AJ, 142, 37

Scoville, N., Lee, N., Vanden Bout, P., et al. 2017, ApJ, 837, 150

Shull, J. M., & Danforth, C. W. 2018, ApJL, 852, L11

Shull, J. M., Danforth, C. W., Tilton, E. M., Moloney, J., & Stevans, M. L. 2017, ApJ, 849, 106

Shull, J. M., Smith, B. D., & Danforth, C. W. 2012, ApJ, 759, 23

Silk, J., & Rees, M. J. 1998, A&A, 331, L1

Sobral, D., Smail, I., Best, P. N., et al. 2013, MNRAS, 428, 1128

Somerville, R. S., Behroozi, P., Pandya, V., et al. 2018, MNRAS, 473, 2714

Somerville, R. S., Hopkins, P. F., Cox, T. J., Robertson, B. E., & Hernquist, L. 2008, MNRAS, 391, 481

Spilker, J. S., Aravena, M., Béthermin, M., et al. 2018, Sci, 361, 1016

Switzer, E. R., Masui, K. W., Bandura, K., et al. 2013, MNRAS, 434, L46

Tacconi, L. J., Genzel, R., Saintonge, A., et al. 2018, ApJ, 853, 179

Tacconi, L. J., Genzel, R., & Sternberg, A. 2020, arXiv:2006.06245

Tacconi, L. J., Neri, R., Genzel, R., et al. 2013, ApJ, 768, 74

Toomre, A., & Toomre, J. 1972, ApJ, 178, 623

Tumlinson, J., Peeples, M. S., & Werk, J. K. 2017, ARA&A, 55, 389

Uzgil, B. D., Carilli, C., & Lidz, A. 2019, ApJ, 887, 37

van de Voort, F., Schaye, J., Booth, C. M., & Dalla Vecchia, C. 2011, MNRAS, 415, 2782

Walch, S., Girichidis, P., Naab, T., et al. 2015, MNRAS, 454, 238

Walter, F., Brinks, E., de Blok, W. J. G., et al. 2008, AJ, 136, 2563

Walter, F., Decarli, R., Aravena, M., et al. 2016, ApJ, 833, 67

Walter, F., Decarli, R., Sargent, M., et al. 2014, ApJ, 782, 79

Werk, J. K., Prochaska, J. X., Tumlinson, J., et al. 2014, ApJ, 792, 8

White, S. D. M., & Frenk, C. S. 1991, ApJ, 379, 52

White, S. D. M., & Rees, M. J. 1978, MNRAS, 183, 341

Wolfe, A. M., Gawiser, E., & Prochaska, J. X. 2005, ARA&A, 43, 861

Wuys, S., Förster Schreiber, N. M., van der Wel, A., et al. 2011, ApJ, 742, 96

Zwaan, M. A., Meyer, M. J., Staveley-Smith, L., & Webster, R. L. 2005, MNRAS, 359, L30