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

We combine the recent determination of the evolution of the cosmic density of molecular gas (H_2) using deep, volumetric surveys, with previous estimates of the cosmic density of stellar mass, star formation rate and atomic gas (HI), to constrain the evolution of baryons associated with galaxies averaged over cosmic time and space. The cosmic HI and H_2 densities are roughly equal at z~1.5. The H_2 density then decreases by a factor 6^{+3}_{-2} to today's value, whereas the HI density stays approximately constant. The stellar mass density is increasing continuously with time and surpasses that of the total gas density (HI and H_2) at redshift z~1.5. The growth in stellar mass cannot be accounted for by the decrease in cosmic H_2 density, necessitating significant accretion of additional gas onto galaxies. With the new H_2 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 HI reservoirs, and (ii) a net inflow of HI and subsequent conversion to H_2 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.


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 mayat 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 build-up 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;Péroux & Howk 2020;Hodge & da Cunha 2020). These include recent measurements from the ALMA Spectroscopic Survey in the Hubble Ultradeep Field (ASPECS; Decarli et al. 2019Decarli et al. , 2020Magnelli 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 2 ) 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 (Sec. 3). We argue that such an approach is justified as molecular gas, star formation, and stellar mass are found to be approximately co-spatial in galaxies, and because the averaging times are significantly longer than the physical processes under consideration (Sec. 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 (Sec. 4). Given the available observational constraints, we here focus on redshifts below z ∼ 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 kms −1 Mpc −1 ) = 0.7, a matter density parameter Ω m = 0.31 (which is the sum of the dark matter density parameter Ω c = 0.259 and the baryon density parameter Ω b = 0.048), and a dark energy density parameter Ω Λ = (1 − Ω 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 build-up by Madau & Dickinson (2014). All the volume-averaged, cosmological densities quoted in this paper are in co-moving units. Fig. 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 (Secs. 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 (Sec. 4).

A SIMPLE SCHEMATIC
Throughout this paper the term 'disk' is used to define this region (with a typical 1 radius r stars < 10 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 disk-like in many cases (e.g., Förster Wuyts et al. 2011;Salmi et al. 2012;Law et al. 2012).
This nominal 'disk' region is surrounded by a reservoir of atomic gas (H I) with radii r HI <50 kpc (Sec. 3.2.2), as demonstrated by observations in the local universe 1 The physical scales quoted here in kpc 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 rstars on r vir see, e.g., Somerville et al. (2018).
(e.g., Walter et al. 2008;Leroy et al. 2009), high redshift observations (Krogager et al. 2017;Neeleman et al. 2017Neeleman et al. , 2019 as well as 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 vir ∼ 50-300 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∼10 4 -10 6 K). The timescale to accrete material from the cool T ∼ 10 4 K CGM is comparable to the dynamical time (∼ 10 8 yr), orders of magnitudes shorter than the cooling time of the hot T ∼ 10 6 K CGM ( 1 Gyr, Sec. 4.5.2). The medium outside this gravitationally collapsed/bound structure (i.e., beyond r 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 Fig. 1, i.e. outflows as well as gas accretion. As detailed below (Sec. 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 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 build-up 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 vs. 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 Fig. 1 should suffice as a representative guide.

MASS COMPONENTS
To put the different baryonic mass components in galaxies in context, we here 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  Figure 1. Schematic of the different baryonic components that are present within the dark matter halo of a galaxy (defined as r < rvir). The central 'disk' region (r < rstars), contains the vast majority of stars and molecular gas, and stars form here at a rate ψstars. This region is surrounded by a reservoir of atomic gas (H I), with r < rHI. The predominantly ionized material (H II) beyond this radius, but within rvir, constitutes the circumgalactic medium (CGM). Beyond rvir is the realm of the intergalactic medium (IGM).
Blue arrows indicate the (net) infall of ionized gas to the H I reservoir ( ψHII→HI) as well as the (net) inflow of atomic gas to the molecular gas (H2) reservoir ( ψHI→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 ψstars).
does not change with time when considering co-moving volumes 2 .
3.1. The z ∼ 0 census For the low redshift Universe (z 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 (∼ 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-α forest. The distribution is thought to be highly filamentary, with the majority of the IGM residing in the 'cosmic web'. Recent work 2 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. 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 ≈ 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; ≥ 10 7 K), seen in X-rays, comprises 4% of the cosmic total 3 . 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, ∼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 Sec. 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α 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 ∼ 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).

Redshift evolution
In the following we discuss the key baryonic mass components in galaxies, and their evolution with cosmic time.  Table 1, and the shaded region marks the 1σ region (16 th to 84 th percentile) from a Monte Carlo Markov Chain analysis (see Sec. 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 text and Madau & Dickinson 2014

Star formation and Stars
The evolution of the cosmic star formation rate density (ψ stars ; Fig. 2, top left) has been constrained through various multi-wavelength studies of large samples of individual galaxies over the last decades (as summarized in the review by Madau & Dickinson 2014). Early studies were based on rest-frame UV observations (e.g., Madau et al. 1996;Lilly et al. 1996;Cucciati et al. 2012;Bouwens et al. 2012a,b, see above review for a complete list of references), and are complemented through observations at longer wavelengths (e.g., Mag-nelli et al. 2011Mag-nelli et al. , 2013Gruppioni et al. 2013;Sobral et al. 2013;Bouwens et al. 2016Bouwens et al. , 2020Novak 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 ∼ 2, with a subsequent decline by a factor of ∼ 8 to the present day. Integrating ψ stars gives the stellar mass formed at a given cosmic time, and this integral is shown as a dashed orange line in Fig. 2

(top right).
This integral can be compared to the independently measured stellar mass density ρ stars (shown as a red line in Fig. 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) This temporal integral of the star-formation rate density lies above the measured stellar mass density ρ stars by a factor 1.4 ± 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 ψ stars (z) is thus the first time derivative of ρ (z), modulo the return fraction 5 of stars R to the interstellar medium through stellar winds and/or supernova explosions (e.g., Madau & Dickinson 2014), i.e., The fact that the integral of ψ stars , after accounting for the return fraction, is in reasonable agreement with ρ stars is remarkable, as highlighted in Madau & Dickinson (2014), if one considers the number of assumptions that go into each measurement 6 . 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).

Atomic gas
The evolution of the cosmic density of atomic gas associated with galaxies (ρ HI (z); Fig. 2, bottom left) has 4 Madau & Dickinson (2014) assume a Salpeter IMF and a lower fixed threshold in luminosity of 0.03 L 5 The return fraction is R = 0.27 for a Salpeter IMF and R = 0.41 for a Chabrier IMF that is more weighted towards massive stars (Madau & Dickinson 2014).
6 But see Hopkins et al. (2018) who argues that this overall agreement does not necessarily imply that the IMF has to be universal. been constrained with several different approaches depending on redshift range. At z ≈ 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 ρ HI . At higher redshifts (0.3 z 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 measure ρ 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) provides an independent mea-surement of ρ HI at these redshifts (e.g., Masui et al. 2013;Switzer et al. 2013).
At z 1.6, H I can be observed using ground-based optical telescopes through its Lyα transition. Quasar absorption spectroscopy of the strongest Lyα absorbers, the so-called damped Lyα systems (DLAs; Wolfe et al. 2005) has yielded estimates of ρ HI up to z ∼ 5 (e.g., Crighton et al. 2015). The ρ 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α 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 µ = 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 ∼ 2 from z ∼ 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 Sec. 2).

Molecular gas
A number of approaches have been followed in the past to constrain the evolution of the cosmic molecular gas density (ρ H2 ; Fig. 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 here include the recent results from ASPECS, that perform deep frequency scans to detect redshifted CO lines without any pre-selection. This approach has been successfully applied in a number of studies Decarli et al. 2014Decarli et al. , 2016Decarli et al. , 2019Decarli et al. , 2020Riechers et al. 2019;Pavesi et al. 2018;Klitsch 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;Magnelli et al. 2020;Dudzevičiūtė 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 . Their metallicities (consistent with solar) and stellar masses (M stars 10 10 M ) 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 ∼ 3 (VLASPECS, Riechers et al. 2020), and were folded into the AS-PECS measurements . Converting dust measurements to molecular gas masses requires the choice of a dust temperature, emissivity, and a dust-togas ratio (see above references).
Stacking and intensity mapping techniques (Inami et al. 2020;Uzgil et al. 2020) 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 highredshift 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 above-mentioned 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 ∼ 2) to today (see also recent reviews by Tacconi et al. 2020;Péroux & Howk 2020;Hodge & da Cunha 2020). There is evidence that the molecular gas density increased from z ∼ 6 to z ∼ 2 (Riechers et al. 2019;Decarli et al. 2019Decarli et al. , 2020, but the associated uncertainties are significant for z > 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 ρ H2 , ρ stars , and ψ stars , we adopt a smooth double The ratio of cosmic molecular-to-atomic gas density as a function of redshift. The ratio peaks at z ∼ 1.5, close to the peak of the star formation rate density. Middle: The ratio of the 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 grey 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  (2014): In order to capture the apparent flattening of the evolution of ρ HI at both low and high redshift, we adopted a hyperbolic tangent function (as in Prochaska & Neeleman 2018) 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 logscale) increase the formal uncertainties derived from the fitting procedure such that >68% of all measurements are contained within the 1σ boundaries (16 th to 84 th percentile). The best fits are shown as solid lines in Figs. 2 and 3, whereas the 1σ boundaries of the fitting functions are shown as colored regions. The fitting parameters are summarized in Table 1. We note that a fit to ρ H2 based on just the ASPECS data gives almost identical parameters as those shown in Table 1.

Cosmic Averages
In the analysis that follows, we will consider the above volume-averaged measurements (Sec. 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 in Sec. 2 and Fig. 1. One can express the quantities discussed here as a function of the well-characterized stellar mass function (SMF) Φ (z, M ) (e.g. Davidzon et al. 2017). Then the cosmic stellar mass density can be written as: where M is the stellar mass. The gas (H 2 or H I) density can then be expressed as: where f gas is the gas-to-stellar mass fraction (f H2 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, feedback processes).

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 ( 0.1 Gyr, Sec. 2) are smaller than the cosmic timescales over which we are averaging (∆z = 1 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 Fig. 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 Sec. 3, the large discrepancy between the total baryon density (dotted curve in Fig. 3) and the baryon density inside galaxies ρ bar,gal (black curve in Fig. 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 ∼ 1.5.
In Fig. 4 (left) we plot the ratio of molecular to atomic gas density as a function of redshift. This ratio peaks at z ∼ 1.5, close to the peak of the star formation rate density. Fig. 4 (middle) shows the ratio of molecular gas-to-stellar mass as a function of redshift. At redshifts z 2 the stellar mass density starts to dominate over the molecular gas density.
The last panel in Fig. 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 (ρ H2 /ψ stars ) is approximately constant above redshifts z 2, and then increases slightly from τ depl ∼ (4 ± 2) × 10 8 yr at z ∼ 2 to τ depl = (7 ± 3) × 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. 2019Aravena et al. , 2020.

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.  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 Péroux & Howk 2020; Tacconi et al. 2020; Hodge & da Cunha 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 ∼ 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 Fig. 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 through the molecular gas phase (Sec. 4.3). We compare this to the total observed gain in stellar mass over the same cosmic time (shown as the red curve in Fig. 5, based on the red curve show in Fig. 2). Even assuming that all the molecular gas ends up in stars, the observed decline in H I and H 2 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 Cosmic age (Gyr) Figure 5. Cumulative gain of the stellar mass density (red line) compared to the cumulative loss of the gas mass density (H2: blue line, total gas: orange line), starting at a redshift of z = 1.5 (TUniv ∼ 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 build-up. The remaining mass to build up the stellar mass must be accreted onto the galaxy.
mass is accounted for, the observed H I and H 2 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.

H II infall and H I inflow rates
Most of the stars are thought to form out of H 2 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 2 . 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 2 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 (Sec. 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 towards the fact that the H I is more extended than the H 2 , i.e. that the DLAs contain negligible molecules (Noterdaeme et al. 2008;Jorgenson et al. 2014;Muzahid et al. 2015).
We here 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, ψ HII→HI , e.g. through cold-mode accretion (Sec. 4.5). In a second step the gas further cools and settles in the central region where it forms H 2 , which we refer to as net inflow, ψ HI→H2 . We stress that we can only consider net rates: it is also possible that H 2 (or H I) is dissociated / photo-ionized 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. ψ HII→HI > 0, ψ HI→H2 > 0. We note that strictly speaking we refer to net flow rate densities (averaged over cosmic volume) throughout this work. For simplicity, we however refer to these as rates throughout.
As detailed in Appendix A, the rate at which the observed H 2 density ρ 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: H2 loss due to feedback + ψ HI→H2 (z)

H2 gain from HI reservoir
( 3) where ψ stars is the star formation rate density, and ψ HI→H2 is the net conversion rate of H I to H 2 ;ρ H2 (z) is the time derivative of ρ H2 (z). The (unknown) mass loading factor ξ 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 here simplistically assume that this ouflow/mass loading is linearly correlated with the star formation rate density, ψ stars (e.g., Spilker et al. 2018;Schroetter et al. 2019) with a universal proportionality factor ξ. The material that is required to replenish the H I reservoir (ρ HI (z) being the time derivative of ρ HI (z)) can be expressed as: where ψ HII→HI is the net infall of gas from the ionized gas phase. As described in Appendix A this expression forρ HI (z) (unlike the one forρ H2 (z)) does not include a mass loading term, as it is included in the net flow term ψ HII→HI . . The H II net infall rate ( ψHII→HI, Eq. 6, orange curve) and H I inflow rate ( ψHI→H2, Eq. 5, black curve) are plotted together with the cosmic star-formation rate density (ψstars, red curve), assuming a mass loading factor of ξ = 0. When including feedback / mass loading (i.e. ξ > 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 H2 densities (ρHI(z) anḋ ρ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 ψHI→H2 and ψHII→HI are a linear combination of the measured quantities: ψstars,ρHI(z), anḋ ρH2(z), as per equations 5 and 6. Below z ≈ 1.5 the inflow rate ψHI→H2 drops below ψstars, as the cosmic H2 reservoir is used up to form stars (negativeρH2(z)).
We can solve equations 3 and 4 for the net inflow rate ψ HI→H2 and the net infall rate ψ HII→HI as a function of observables ρ HI (z), ρ H2 (z) and ψ stars : and ψ HII→HI =ρ HI (z) +ρ H2 (z) + (1 + ξ) ψ stars (z). (6) In Fig. 6 we plot these net flows rates (equations 5 and 6), along with the star formation rate density ψ stars . We also show the time derivatives of ρ HI and ρ 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 ( ψ HII→HI ) and H I inflow rate ( ψ HI→H2 ) are larger than the star 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 1.5 the net inflow rate ψ HI→H2 is higher than ψ 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 ρ H2 , while HI is slowly being depleted. Conversely, at z 1.5 the net inflow rate ψ HI→H2 is lower than ψ 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.

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 (Sec. 3.3) to forecast the evolution of the baryon content associated with galaxies over the next few Gyr. This is shown in Fig. 7 where we plot the same information as in Fig. 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 H I 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 (Eqs. 5 and 6) will decrease correspondingly. Fig. 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.

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 ( ψ HII→HI and ψ HI→H2 ). We also consider the potential role of preventive feedback mechanisms (such as virial shocks, AGN feedback, and cosmic expansion). Cosmic age (Gyr) Figure 7. We here plot the same information as in Fig. 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 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.

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Ṁ matter (M vir , z) as a function of halo virial mass and redshift using the fitting function presented in Rodríguez-Puebla et al.
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 Fig. 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 ) corresponds to the mass resolution in the simulations considered (corresponding to a stellar mass of a few times 10 7 M ). 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Ṁ 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 Fig. 8. It is interesting to note that these accretion curves show a similar shape as our derived net infall/inflow rates ( ψ HI→H2 , ψ HII→HI ): the accretion rates rise from high redshift to about z ∼ 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 towards z = 0 is due to the fact that the Universe expands and to the gradual decrease in the availability of accretable (dark) matter 7 .

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 Fig. 8 (same curves as in Fig. 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 (Sec. 4.3). Consequently, the majority of the accreted baryons will not make it to the central galaxy 'disks'.
An extensive 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, referred to as 'cold mode accretion', occurs on timescales of order a free-fall or the dynamical crossing time of a spherical halo (∼10 8 yr, depending on mass), where the separation between the 'cold' and 'hot' phase is around 10 5.5 K. The fraction of the 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;Pan et al. 2019;van de Voort et al. 2011;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 7 yr (e.g., Cornuault et al. 2018). Our observations do not allow to distinguish between the different accretion mechanisms ('cold' vs '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, photo-ionization, through, e.g. stellar and AGN feedback. This hot material has very long cooling times. Assuming a typical temperature ∼ 10 6 K and a density ∼ 10 −5 cm −3 , with substantial variation (e.g., Shull et al. 2012Shull et al. , 2017Nicastro et al. 2018), the nominal bremsstrahlung cooling time is about 8.5 × 10 11 yr, or more than sixty 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 'tappedout' 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 8 .

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 2 density is less than or equal to the cosmic H I density over all times, briefly approaching equality at z ∼ 1.5. Below a redshift of z ∼ 1.5, the stellar mass density begins to dominate over all gas components (H 2 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 2, and then increases slightly from τ depl ∼ (4 ± 2) × 10 8 yr at z ∼ 2 to τ depl = (7 ± 3) × 10 8 yr at z = 0. Significant ac-8 The situation is similar to the conclusion of the pioneering work by Toomre & Toomre (1972) on the pre-destiny 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.' cretion 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 ∼ 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 2 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 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 ψ HII→HI . 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 ψ HI→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 starformation 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 ψ 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 ∼ 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 2 , 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. FW and MN acknowledge support from the ERC Ad-vanced Grant 740246 (Cosmic Gas). BM acknowledges support from the Collaborative Research Centre 956, sub-project A1, funded by the Deutsche Forschungsgemeinschaft (DFG) -project ID 184018867. TD-S acknowledges support from the CASSACA and CON-ICYT fund CAS-CONICYT Call 2018. RJA was supported by FONDECYT grant 1191124. DR 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. IRS acknowledges support from STFC (ST/T000244/1). HI acknowledges support from JSPS KAKENHI Grant Number JP19K23462. JH acknowledges support of the VIDI research programme with project number 639.042.611, which is (partly) fi-nanced by the Netherlands Organisation for Scientific Research (NWO). DO is a recipient of an Australian Research 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 Likewise, we assume that stars do not produce atomic hydrogen nor molecular gas, and we thus set these two terms (ψ stars→HI , ψ stars→H2 ) to zero as well. Equation. A1 therefore simplifies to: ρ stars = ψ H2→stars − ψ stars→HII (A5) ρ HII : In Eq. A2 the first two terms correspond to the mass flows between the H 2 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 here assume that ionized gas cannot directly form molecular gas, and that ionized gas can not directly form stars and set ψ HII→stars and ψ HII→H2 to zero. The other flows are in principle plausible, and equation A2 thus becomes:ρ HII = ψ H2→HII − ψ HII→HI + ψ HI→HII + ψ stars→HII (A6) where in a second step we introduced ψ HII→HI as being the net flow between the ionized gas phase and the atomic gas.
As we assume that |ψ HII→HI | > |ψ HI→HII | we assign a minus sign to the net flow ψ HII→HI in Eq. A7. ρ HI : In Eq. 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 ψ HI→stars and ψ stars→HI are set to zero. Eq. A3 then yields: where in a second step we again introduce an additional net flow between the atomic gas and the H 2 ( ψ HI→H2 ). Its sign is determined as we assume |ψ H2→HI | < |ψ HI→H2 |. This equation is identical to Eq. 4 in the main body of this manuscript. ρ H2 : In Eq. 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 ψ HII→H2 and ψ stars→H2 to zero. This yields: ρ H2 = −ψ H2→HII − ψ H2→stars + ψ HI→H2 − ψ H2→HI (A10) ρ H2 = −ψ H2→HII − ψ H2→stars + ψ HI→H2 (A11) where we again use the previously defined net flow between the atomic gas and the H 2 , ψ HI→H2 , in the second step.
We can now further simplify equation A11 by setting ψ H2→stars to the star formation rate density, i.e. ψ H2→stars = ψ 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), ψ H2→HII , to be proportional to the star formation rate density, i.e. ψ H2→HII = ξψ stars . This yields Eq. 3 in the main text.
As a side note, assuming that the return rate of the stars to the ionized medium (ψ stars→HII ) is proportional to the star formation rate density ψ stars with a proportionality factorR, i.e. ψ stars→HII =Rψ stars , we can rewrite Eq. A5 that expresses the change in the stellar mass density as follows: ρ stars = ψ stars −Rψ stars (A12) Note that the factorR would be equal to the classical return factor R (see Sec. 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. In this appendix, we give the observational measurements for both the cosmic H I mass density and the cosmic H 2 mass density 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).