{'fraction_non_alphanumeric': 0.03522762548035238, 'fraction_numerical': 0.02420226715620685, 'mean_word_length': 4.27599254532106, 'pattern_counts': {'":': 0, '<': 71, '<?xml version=': 0, '>': 40, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Which galaxies dominate the neutral gas content of the Universe?

We study the contribution of galaxies with different properties to the global densities of star formation rate (SFR), atomic (HI) and molecular hydrogen (H2) as a function of redshift. We use the GALFORM model of galaxy formation, which is set in the LCDM framework. This model includes a self-consistent calculation of the SFR, which depends on the H2 content of galaxies. The predicted SFR density and how much of this is contributed by galaxies with different stellar masses and infrared luminosities are in agreement with observations. The model predicts a modest evolution of the HI density at z<3, which is also in agreement with the observations. The HI density is predicted to be always dominated by galaxies with SFR<1Msun/yr. This contrasts with the H2 density, which is predicted to be dominated by galaxies with SFR>10Msun/yr. Current high-redshift galaxy surveys are limited to detect carbon monoxide in galaxies with SFR>30Msun/yr, which in our model make up, at most, 20% of the H2 in the universe. In terms of stellar mass, the predicted H2 density is dominated by massive galaxies, Mstellar>10^10Msun, while the HI density is dominated by low mass galaxies, Mstellar<10^9Msun. In the context of upcoming neutral gas surveys, we suggest that the faint nature of the galaxies dominating the HI content of the Universe will hamper the identification of optical counterparts, while for H2, we expect follow up observations of molecular emission lines of already existing galaxy catalogues to be able to uncover the H2 density of the Universe.


INTRODUCTION
Observations of local and high-redshift galaxies indicate that stars form from molecular gas, albeit with a low efficiency (e.g. Bigiel et al. 2008;Genzel et al. 2010). On scales of kilo-parsecs, this fact is observed as a close to linear correlation between the star formation rate (SFR) surface density and the surface density of molecular hydrogen (H2) Bigiel et al. 2010;Leroy et al. 2008;Schruba et al. 2011). The observational and theoretical evidence gathered point to a scenario in which gas cools down and becomes molecular before collapsing and fragmenting to form stars, with the small efficiency explained as being due to self-regulation as the results of star formation at the level of giant molecular clouds. Therefore, the presence of neutral gas by itself does not guarantee star formation (see Kennicutt & Evans 2012 for a recent review). This emphasises the importance of H2 as the main tracer of the dense regions where stars form. Regarding the atomic hydrogen (HI), studies of local Universe galaxies show a strong correlation between the HI mass and the stellar mass, indicating that the evolution of the two quantities is related (e.g. Catinella et al. 2010;Cortese et al. 2011;Huang et al. 2012). In addition, the morphology of the atomic hydrogen in galaxies is closely related to baryonic processes such as gas accretion and outflows in galax-ies (Fraternali et al. 2002;Oosterloo et al. 2007;Boomsma et al. 2008). Observations have shown that the presence of stellar driven outflows depends on the star formation rate density (e.g. Chen et al. 2010;Newman et al. 2012). Lagos et al. (2013) and Creasey et al. (2013), using hydrodynamical modelling of the growth of supernovae driven bubbles in the ISM of galaxies, show that the fundamental property setting the outflow rate is the gas surface density, as it affects both the SFR, which sets the energy input from supernovae, and the time bubbles take to escape the galaxy disk.
All of the evidence above points to the fact that a key step in understanding galaxy formation is the observation of multiple gas phases in the ISM and their relation to the presence of SF. In particular, observations of atomic and molecular hydrogen, which make up most of the mass in the interstellar medium (ISM) of galaxies, are key to understanding how SF proceeds and how gas abundances can be linked to accretion and outflow of gas. Only through tracing both molecular and atomic gas components at the same time as the SFR and stellar mass in galaxies will we be able to develop a better understanding of the processes of SF and feedback and to put strong constraints on galaxy formation simulations.
Currently, observational constraints on the HI and H2 contents of galaxies are available for large local samples, and for increasingly large samples at high-redshift. In the case of HI, ac-curate measurements of the 21 cm emission in large surveys of local galaxies have been presented by Zwaan et al. (2005) using the HI Parkes All-Sky Survey (HIPASS; Meyer et al. 2004) and by Martin et al. (2010) using the Arecibo Legacy Fast ALFA Survey (ALFALFA; Giovanelli et al. 2005). The global HI mass density at z = 0 has been estimated from these surveys to be in the range ΩHI = 3.6 − 4.2 × 10 −4 , with little evolution up to z ≈ 5 (e.g. Péroux et al. 2003;Noterdaeme et al. 2009). This lack of evolution is fundamentally different from the evolution inferred in the SFR density, which shows a strong increase from z ≈ 0 to z ≈ 2 − 3, followed by a slow decline to higher redshifts (e.g. Madau et al. 1996;Lilly et al. 1996;Hopkins & Beacom 2006). A plausible explanation for the different evolution seen in the SFR and HI densities is that HI is not the direct fuel of SF but instead this gas still needs to cool down further to form stars. This indicates that in order to better understand SF and galaxy evolution, good observations of the denser gas in the ISM, i.e. H2, are required. Molecular hydrogen is commonly traced by the 12 CO (hereafter 'CO') molecule, which is the second most abundant molecule after H2 and is easily excited. However, direct CO detections are very scarce and complete samples are limited to the local Universe. Keres et al. (2003) reported the first attempt to derive the local CO(1 − 0) luminosity function from which they inferred the H2 mass function and the local ΩH 2 = 1.1 ± 0.4 × 10 −4 h −1 , adopting a Milky-Way type CO(1 − 0)-H2 conversion factor. It has not yet been possible to estimate the cosmic H2 abundance at high redshift. Observational samples that detect CO at high-redshifts are limited to moderately and highly star-forming galaxies with SFR 20 M⊙ yr −1 (Tacconi et al. 2013;. To estimate a density of H2 from these galaxies is difficult given the selection effects and volume corrections. Recently, Berta et al. (2013) inferred the molecular gas mass function at high-redshift by using the ratio between the UV to mid-IR emission and its empirical relation to the molecular gas mass . However, these estimates are subject to the uncertain extrapolation of empirical relations to galaxies and redshifts where they have not been measured. For these reasons the construction of representative samples of H2 in galaxies is still an unfinished task that is crucial to our understanding of galaxy formation.
From the theoretical point of view, more sophisticated treatments of SF and its relation to the abundance of HI and H2 have only recently started to be explored in ab-initio cosmological galaxy formation simulations (e.g. Fu et al. 2010;Cook et al. 2010;Lagos et al. 2011b;Kuhlen et al. 2012;Duffy et al. 2012a;Popping et al. 2013a), as well as in simulations of individual galaxies (e.g. Gnedin et al. 2009;Pelupessy & Papadopoulos 2009;Kim et al. 2011;Mac Low & Glover 2012). The improved modelling of the ISM and SF in galaxies that has been applied in cosmological galaxy formation simulations has allowed a better understanding of the neutral content of galaxies; for instance, of the effect of H2 self-shielding on the column density distribution of HI (Altay et al. 2011;Duffy et al. 2012a), the lack of evolution of the global HI density (Lagos et al. 2011a;Davé et al. 2013), and the increasing molecular gas-fractions of galaxies with redshift Lagos et al. 2012;Fu et al. 2012).
The issue motivating this paper is the contribution of starforming galaxies of different properties to the atomic and molecular gas content of the Universe. Note that this question is similar to asking about the steepness of the faint-end of the gas and stellar mass functions; i.e. the steeper the mass function, the larger the contribution from low-mass galaxies to the total density of mass in the Universe. This exercise is important for two reasons. First, blind CO surveys are very challenging and expensive. For instance, the Atamaca Large Millimeter Array can easily detect galaxies at high-redshift (e.g. Vieira et al. 2013;Hodge et al. 2013), but has the downside of having a very small field of view (< 1 arcmin), which is not ideal for large surveys. A possible solution is to follow-up samples of galaxies selected by SFR or stellar mass in molecular emission, which can give insight into the relation between SF and the different gas phases in the ISM. Second, because a key science driver of all large HI programs is to understand how the HI content relates to multi-wavelength galaxy properties (see for instance the accepted proposals of the ASKAP HI All-Sky Survey 1 , WALLABY, and the Deep Investigation of Neutral Gas Origins survey 2 , DINGO, Johnston et al. 2008). DINGO will investigate this directly as it overlaps with the Galaxy And Mass Assembly multiwavelength survey (Driver et al. 2009).
In this paper we bring theoretical models a step closer to the observations and examine how star-forming galaxies of different properties contribute to the densities of SFR, atomic and molecular hydrogen, to inform future neutral gas galaxy surveys about the expected contribution of star-forming galaxies to the HI and H2 contents of the Universe. We also seek to understand how far we currently are from uncovering most of the gas in galaxies in the Universe. For this study we use three flavours of the semianalytical model GALFORM in a ΛCDM cosmology , namely those of Lagos et al. (2012), Gonzalez-Perez et al. (2013), and Lacey et al. (2014, in prep.). The three models include the improved treatment of SF implemented by Lagos et al. (2011b). This extension explicitly splits the hydrogen content of the ISM of galaxies into HI and H2. The advantage of using three different flavours of GALFORM is the ability to characterise the robustness of the trends found. The outputs of the three models shown in this paper will be publicly available from the Millennium database 3 (detailed in Appendix C). This paper is organised as follows. In §2 we present the galaxy formation model and describe the realistic SF model used. We also describe the main differences between the three models that use this SF law and the dark matter simulations used. In §3 we compare the model predictions with observations of the evolution of the SFR, HI and H2 densities, and the SF activity in normal and highly star-forming galaxies. None of the observations in §3 were used to constrain model parameters, and therefore provide independent verifications of the predictions of the models. Given the success of our model predictions, we investigate the contribution from star-forming galaxies to the densities of HI and H2 in §4 and discuss the physical drivers of the predictions. In §5 we discuss the implications of our predictions for the next generation of HI and H2 surveys and present the main conclusions in §6.

MODELLING THE TWO-PHASE COLD GAS IN GALAXIES
In this section, we briefly describe the main aspects of the GALFORM semi-analytical model of galaxy formation and evolution , focusing on the key features of the three flavours adopted in this study, which are described in Lagos et al. (2012;hereafter 'Lagos12'), Gonzalez-Perez et al. (2013; hereafter 'Gonzalez-Perez14') and Lacey et al. (2014; hereafter 'Lacey14').
The three GALFORM models above take into account the main physical processes that shape the formation and evolution of galaxies. These are: (i) the collapse and merging of dark matter (DM) halos, (ii) the shock-heating and radiative cooling of gas inside DM halos, leading to the formation of galactic disks, (iii) quiescent star formation in galaxy disks, (iv) feedback from supernovae (SNe), from active galactic nucleus (AGN) heating and from photoionization of the inter-galactic medium (IGM), (v) chemical enrichment of stars and gas, and (vi) galaxy mergers driven by dynamical friction within common DM halos which can trigger bursts of SF, and lead to the formation of spheroids (for a review of these ingredients see Baugh 2006;Benson 2010). Galaxy luminosities are computed from the predicted star formation and chemical enrichment histories using a stellar population synthesis model (see Gonzalez-Perez et al. 2013).
The three flavours of GALFORM used in this study adopt the same SF law, which is a key process affecting the evolution of the gas content of galaxies. In § 2.1 we describe this choice of SF law and how this connects to the two-phase ISM, in § 2.2 we describe the differences between the three flavours of GALFORM and in §2.3 we briefly describe the dark matter simulations used for this study and the cosmological parameters adopted.

Interstellar medium gas phases and the star formation law
The three flavours of GALFORM used in this study adopt the SF law developed in Lagos et al. (2011a, hereafter 'L11'), in which the atomic and molecular phases of the neutral hydrogen in the ISM are explicitly distinguished. L11 found that the SF law that gives the best agreement with the observations without the need for fine tunning is the empirical SF law of Blitz & Rosolowsky (2006). Given that the SF law has been well constrained in spiral and dwarf galaxies in the local Universe, L11 decided to implement this molecularbased SF law only in the quiescent SF mode (SF due to gas accretion onto the disk), keeping the original prescription of Cole et al. (2000) for starbursts (driven by galaxy mergers and global disk instabilities). We provide more details below. Quiescent Star Formation. The empirical SF law of Blitz & Rosolowsky has the form, where ΣSFR and Σgas are the surface densities of SFR and total cold gas mass (molecular and atomic), respectively, νSF is the inverse of the SF timescale for the molecular gas and f mol = Σ mol /Σgas is the molecular to total gas mass surface density ratio. The molecular and total gas contents include the contribution from helium, while HI and H2 only include hydrogen (which in total corresponds to a fraction XH = 0.74 of the overall cold gas mass). The ratio f mol depends on the internal hydrostatic pressure as ΣH 2 /ΣHI = f mol /(f mol − 1) = (Pext/P0) α . To calculate Pext, we use the approximation from Elmegreen (1989), in which the pressure depends on the surface density of gas and stars. We give the values of the parameters involved in this SF law in § 2.2.
Starbusts. In starbursts the SF timescale is proportional to the bulge dynamical timescale above a minimum floor value and involves the whole cold gas content of the galaxy, SFR = M cold /τSF (see Granato et al. 2000 andLacey et al. 2008 for details). The SF timescale is defined as where τ dyn is the bulge dynamical timescale, τmin is a minimum duration of starbursts and f dyn is a parameter. We give the values of the parameters involved in this SF law in § 2.2. Throughout this work we assume that in starbursts, the cold gas content is fully molecular, f mol = 1. Note that this is similar to assuming that the relation between the ratio ΣH 2 /ΣHI and Pext holds in starbursts given that large gas and stellar densities lead to f mol ≈ 1. Throughout the paper we refer to galaxies going through a starburst, which are driven by galaxy mergers or disk instabilities, as 'starbursts'.
Another key component of the ISM in galaxies is dust. To compute the extinction of starlight we adopt the method described in Lacey et al. (2011) and Gonzalez-Perez et al. (2012), which uses the results of the physical radiative transfer model of Ferrara et al. (1999) to calculate the dust attenuation at different wavelength. The dust model assumes a two-phase interstellar medium, with starforming clouds embedded in a diffuse medium. The total mass of dust is predicted by GALFORM self-consistently from the cold gas mass and metallicity, assuming a dust-to-gas ratio which is proportional to the gas metallicity. The radius of the diffuse dust component is assumed to be equal to the half-mass radius of the disk, in the case of quiescent SF, or the bulge, in the case of starbursts. With this model, we can estimate the attenuation in the wavelength range from the far-ultraviolet (FUV) to the near-infrared (IR), including emission lines. In GALFORM the extinction for emission lines is calculated as due to the diffuse dust component only, so it is the same as the extinction of older stars (i.e. outside their birth clouds) at the same wavelength. We define the total IR luminosity as the total luminosity emitted by interstellar dust, free from contamination by starlight, which approximates to the integral over the rest-frame wavelength range 8−1000 µm.

Differences between the Lagos12, Gonzalez-Perez14 and Lacey14 models
The Lagos12 model is a development of the model originally described in Bower et al. (2006), which was the first variant of GALFORM to include AGN feedback as the mechanism suppressing gas cooling in massive halos. The Lagos12 model assumes a universal initial mass function (IMF), the Kennicutt (1983) IMF 4 . Lagos12 extend the model of Bower et al. by including the self-consistent SF law described in § 2.1, and adopting νSF = 0.5 Gyr −1 , log(P0/kB[cm −3 K]) = 4.23, where kB is is Boltzmann's constant, and α = 0.8, which correspond to the volume of the parameters reported by Leroy et al. (2008) for local spiral and dwarf galaxies. This choice of SF law greatly reduces the parameter space of the model and also extends its predictive power by directly modelling the atomic and molecular hydrogen content of galaxies. All of the subsequent models that make use of the same SF law have also the ability to predict the HI and H2 gas contents of galaxies. Lagos12 adopt longer duration starbursts (i.e. larger f dyn ) compared to Bower et al. to  Eq. 2. The Lagos12 model was developed in the Millennium simulation, which assumed a WMAP1 cosmology (Spergel et al. 2003).
The Gonzalez-Perez14 model updated the Lagos12 model to the WMAP7 cosmology (Komatsu et al. 2011) and vary few model parameters to help recover the agreement between the model predictions and the observed evolution of both the UV and K-band luminosity functions. These changes include a slightly shorter starburst duration, i.e. τmin = 50 Myr and f dyn = 10, and weaker supernovae feedback.
The Lacey14 models is also developed on the WMAP7 cosmology, as the Gonzalez-Perez14 model, but it differs from both the Lagos12 and Gonzalez-Perez14 models in that it adopts a bimodal IMF. The IMF describing SF in disks (i.e. the quiescent mode) is the same as the universal IMF in the other two models, but a topheavy IMF is adopted for starbursts (i.e. with an IMF slope x = 1). This is inspired by Baugh et al. (2005) who used a bimodal IMF to recover the agreement between the models predictions and observations of the number counts and redshift distribution of submillimeter galaxies. We notice, however, that Baugh et al. adopted a more top-heavy IMF for starbursts with x = 0. The stellar population synthesis model for Lacey14 is also different. While both Lagos12 and Gonzalez-Perez14 use Bruzual & Charlot (2003), the Lacey14 model makes use of Maraston (2005). Another key difference between the Lacey14 model and the other two GALFORM flavours considered here, is that Lacey14 adopt a slightly larger value of the SF timescale, νSF = 1.2 Gyr −1 , still within the range allowed by the most recent observation compilation of Bigiel et al. (2011), making SF more efficient.

Halo merger trees and cosmological parameters
GALFORM uses the formation histories of DM halos as a starting point to model galaxy formation (see Cole et al. 2000). In this paper we use Monte-Carlo generated merger histories from Parkinson et al. (2008). We adopt a minimum halo mass of M halo,min = 5×10 8 h −1 M⊙ and model a representative sample of dark matter halos such that the predictions presented for the mass functions and integrated distribution functions are robust. The minimum halo mass chosen enables us to predict cold gas mass structures down to the current observational limits (i.e. MHI ≈ 10 6 M⊙, Martin et al. 2010). At higher redshifts, this minimum halo mass is scaled with redshift as (1 + z) −3 to roughly track the evolution of the break in the halo mass function, so that we simulate objects with a comparable range of space densities at each redshift. This allows us to follow a representative sample of dark matter halos. We aim to study in detail the HI and H2 contents of galaxies in the universe. To achieve this, the use of Monte-Carlo trees is compulsory as the range of halo masses needed to do such a study is much larger than the range in currently available N -body simulations, such as the Millennium (Springel et al. 2005) and the Millennium-II 5 (Boylan- Kolchin et al. 2009). In Appendix C, we discuss in detail the effect of resolution on the neutral gas content of galaxies and show that the chosen minimum halo mass for this study is sufficient to achieve our goal.
Lagos12 adopt the cosmological parameters of the Millennium N-body simulation (Springel et al. 2005

THE EVOLUTION OF SFR, ATOMIC AND MOLECULAR GAS DENSITIES
In this section we explore the predictions for the evolution of the Hα and UV luminosity functions, and the global densities of SFR, HI and H2 ( § 3.1). The aim of this is to establish how well the models describe the star-forming galaxy population at different redshifts and if the overall predicted gas content of the universe is in agreement with observational estimates. We then focus in one model, the one that best reproduces the observations, to study the contribution to the SFR density from galaxies of different stellar masses and IR luminosities in § 3.2 and § 3.3.

The Hα and UV luminosity functions
A fair comparison between the observations and the model predictions for the SFR can be made by directly contrasting the observed Hα and UV luminosity functions with the predicted ones, which correspond to two of the most popular SFR tracers. Fig. 1 shows the intrinsic Hα luminosity function from z ≈ 0 to z ≈ 1.5 for the three models and observations. The observations have been separated into those which present the Hα luminosity function including a correction for dust attenuation (filled symbols), and those which do not (open symbols). For the three models we show both the intrinsic and dust extincted Hα luminosity functions. The intrinsic luminosity functions predicted by the models should be compared to the filled symbols. The three models have predictions in reasonable agreement with the observations, within the errorbars. At z = 0 and z = 0.3, the models predict Hα luminosity functions in broad agreement with the observations. At z = 1, z = 1.5 and z = 2 we only show direct observations of the Hα luminosity function from Sobral et al. (2013), who present the largest galaxy samples to date, greatly exceeding the number of galaxies and observed sky area from previous works (i.e. 14, 011 emitters in 2.2 deg 2 ). These observations are directly comparable to the model predictions that include dust. The three models are consistent with the z ≈ 1 observations, but are below the observations at z ≈ 1.5 and z ≈ 2 at the faintend by a factor of ≈ 2 − 3. This may be due either to a low number of highly star-forming galaxies in the models or to the dust extinction modelling. The former would change the intrinsic Hα luminosity function, while the latter would affect only the extincted Hα emission. We expect both to contribute to some extent. Differences between the three models in Fig. 1 are mainly seen at the bright-end where the errorbars on the observations are the largest. In general, our dust model predicts that dust attenuation is luminosity dependent, which is a natural consequence of the different gas metallicities in our model galaxies. The faintend is expected to be only slightly affected by dust, but bright Hα galaxies can be highly attenuated by dust. However, the predicted extinctions for z 1 are typically smaller than that assumed in observations (AHα = 1 mag). The Lacey14 model predicts higher number densities of luminous Hα galaxies, with intrinsic   Ly et al. (2007) and Gunawardhana et al. (2013) are shown in the top two panels, and from Sobral et al. (2013) in the bottom three panels using symbols, as labelled. Filled symbols represent observational data that has been corrected for dust extinction (i.e. comparable to the model intrinsic luminosity function), and open symbols represent the observed luminosity function (i.e. comparable to the model extinct luminosity function). Note that we express densities and luminosities in comoving units for both the observations and the models. To transform these units to physical units the reader needs to evaluate h with the value given in § 2.3.  The rest-frame UV(1500Å) luminosity function at z = 2 (top panel), z = 3 (middle panel) and z = 6 (bottom panel) for the Lagos12, Lacey14 and Gonzalez-Perez14 models, as labelled. Symbols represent the observational estimates from Sawicki & Thompson (2006), Bouwens et al. (2007), Reddy &Steidel (2009) andMcLure et al. (2009), as labelled. The predictions shown for the models already include the effects of dust extinction and should be compared to the observed luminosity functions without any correction for extinction.The units of density here are as in Fig. 1. However, to convert the magnitudes from the adopted units to M AB (UV) the reader needs to add 5log h to the quantity in the x-axis, where h has the value showed in § 2.3.
LHα > 5 × 10 42 h −2 erg s −1 , than the other two models. This is due the predominance of starbursts at these high luminosities and the top-heavy IMF adopted by Lacey14 in bursts. Note that the three models predict that the normalisation and break of the Hα luminosity function increase with increasing redshift between z = 0 and z = 1, while at z > 1, the evolution is mainly in the number density of bright galaxies. This overall evolution agrees well with that observed.
At high redshifts, a popular way of studying SF in galaxies is to use the rest-frame FUV emission as a tracer. We compare the predicted rest-frame FUV luminosity function in the three flavours of GALFORM to observations from z = 2 to z = 6 in Fig. 2. The observed luminosity function corresponds to the rest-frame 1500Å luminosity, without any attempt to correct for dust extinc-tion. This is directly comparable to the model predictions for the rest-frame 1500Å luminosity, after dust has been included in the model (lines in Fig. 2). The predictions from the three flavours of GALFORM agree well with the observations, even up to z = 6. This is partially due to the assumed duration of starbursts. Lacey et al. (2011) show that very short starbursts, such as the ones assumed in Bower et al. (2006), result in the bright-end of the UV luminosity function being largely over predicted. Longer starbursts allow the model to have more dusty starbursts and reproduce the break in the UV LF. The duration of the bursts is of about few 100 Myr, and therefore in better agreement with the latest estimates of the duration from observational data .

The cosmic density of SFR
We are interested in connecting the evolution of the gas content of galaxies with that of the SFR. With this in mind, we show in Fig. 3 the evolution of the global density of SFR, HI and H2 for the Lagos12, Lacey14 and Gonzalez-Perez14 models. In the case of the SFR, the observational estimates have been corrected to our choice of IMF (Kennicutt 1983) (see Appendix A for the conversions). Note that we express the SFR density in comoving units as in Fig. 1. The predictions of the three models are in reasonable agreement with the observational estimates at z < 0.3 and at z > 2, within the uncertainties, but significantly lower than the observational estimates in the redshift range 0.3 < z < 2. The comparison between the model predictions and the observed density of SFR needs to be undertaken carefully as there are many systematic uncertainties which are not taken into account in the errorbars of the observations. One example is the dust absorption correction applied when ultraviolet bands or emission lines are used to trace SFR (e.g. Hopkins & Beacom 2006). Typically, empirical dust corrections calibrated in the local Universe are extrapolated to highredshifts (e.g. Calzetti et al. 2007). Another source of uncertainty is related to the extrapolation needed to correct the observed SFR density to account for faint galaxies that are not detected. These two factors drive most of the dispersion observed in Fig. 3 between the different estimates.
To help illustrate the impact of the extrapolation to fainter luminosities and the dust correction, we show in the top panel of Fig. 3 the SFR density from Sobral et al. (2013) which includes both the dust correction and the extrapolation to faint galaxies (squares), the SFR density with the dust correction but integrated only over the range covered by the measurements (circles with inner crosses) and the latter but without the dust correction (squares with inner crosses). The overall effect of the dust correction is at most of ≈ 0.4 dex and of the extrapolation can be an increase of up to ≈ 0.6 dex. If the true slope of the Hα LF is shallower and/or the dust correction is less dramatic than assumed in the observations the inferred ρSFR would easily move down closer to our predictions. Recently, Utomo et al. (2014) show evidence that the dust corrections applied in most of the observations shown in the top panel of Fig. 3 are overestimated, leading to an overestimation of the SFRs. Other systematics are hidden in the adopted IMF. Although we scale the inferred observations to our adopted IMF, this conversion is valid only in the case of a simple SF history, for instance a constant SFR, and a constant metallicity, at a certain stellar population age. This is not the case for our simulated galaxies that undergo starbursts and gas accretion which can change the available gas content at any time (Mitchell et al. 2013). Otí-Floranes & Mas-Hesse (2010) argue that the systematic uncertainty due to the IMF choice can be up to a factor of ≈ 4.  Wilkins et al. (2012) show that the ratio between the SFR and the FUV luminosity is correlated with SFR, and that systematic errors relating to adopting a universal ratio can be a factor of ≈ 1.5 or more. A good example of the uncertainties in the observationally inferred ρSFR is the fact that the three models predict a Hα luminosity function at z 1 in broad agreement with the observations (see Fig. 1), but a ρSFR that is a factor of ≈ 2 − 4 lower than the Hα inferred ρSFR. Thus, we consider the SFR density inferred from observations as a rough indicator of the SF activity as opposed to a stringent constrain to our model predictions. The fact that the three models have ρSFR at z ≈ 1 below the observations is not currently of great concern for two reasons: (i) given the difficulty of quantifying the effect of systematic biases on the observationally inferred ρSFR, and (ii) because the predicted Hα and UV luminosity functions are in reasonable agreement with observations. Overall, the three models predict ρSFR of a similar shape and normalisation. However, in the detail, we can distinguish two differences that are worth exploring: (i) the Lacey14 model predicts a lower peak of ρSFR compared to the other two models, and (ii) at z > 6, there is an offset of a factor of 2−3 between ρSFR predicted by the Lagos12 model and the other two models. Feature (i) is due to the adopted top-heavy IMF in starbursts in the Lacey14 model, which results in lower SFRs to drive the same level of FUV, IR and emission line fluxes compared to the choice of a Kennicutt IMF. In fact, we can calculate the difference between the two different IMFs in the flux predicted for the different SFR tracers (see Appendix A). We can scale the contribution from starbursts up by 1.8 to account for the different IMF; this SFR density would roughly correspond to the inferred one if a Kennicutt IMF was adopted as universal and UV emission was used as a tracer of SF. This is shown by the dotted line in the top panel of Fig. 3. This scaling accounts for most of the difference between the Lacey14 and the other two models that adopt a universal Kennicutt IMF. Feature (ii) is due to the different cosmologies adopted in the Lagos12 (WMAP1) and the other two models (WMAP7). The WMAP7 cosmology has a lower σ8, which delays the collapse of the first structures, giving rise to the first galaxies at lower redshifts compared to the WMAP1 universe (see Gonzalez-Perez et al. 2013 andGuo et al. 2013 for an analysis of the effect of these cosmologies on galaxy properties). Note that the difference in Ωm reported by WMAP1 and WMAP7 is irrelevant at these high redshifts as Ωm ≈ 1.

The cosmic densities of HI and H2
The middle and bottom panels of Fig. 3 shows that the atomic and molecular gas densities are predicted to have very different evolution in the three GALFORM models. Note that we have converted all of the observations to our adopted cosmology ( § 2.3). The predicted atomic gas density displays very little evolution at z 2.5. At higher redshifts, ρHI slowly decreases with increasing redshift. This behaviour is common in the three models. Quantitatively, ρHI decreases by only ≈ 11% between z = 0 and z = 2.5 in the three models, while between z = 2.5 and z = 6 it decreases by 68% in the Lagos12 model, 82% in the Gonzalez-Perez14 model and 80% in the Lacey14 model. At z 3, the predicted evolution of the three models agrees well with observations of ρHI that come either from direct HI detections at z ≈ 0 − 0.2 or inferences from spectral stacking, intensity mapping and damped-Lyα systems (DLAs) at z 0.2. At z > 3, the predictions for ρHI from the three models systematically deviate from the observed ρHI. One reason for this is that we only model neutral gas inside galaxies, i.e. in the ISM, and we do not attempt to estimate the neutral fraction of the gas outside galaxies (i.e. in the IGM). This argument is supported by the results of hydrodynamical simulations which explicitly follow the neutral fraction in different environments. For instance, van de Voort et al. (2012) and Davé et al. (2013) show that the neutral content of the universe in their simulations is dominated by gas outside galaxies at z 3.
The weak evolution of ρHI shown in Fig. 3 contrasts with the evolution of the molecular hydrogen global density, ρH 2 , which shows a steep increase between z = 0 and z = 3 followed by a slow decrease as redshift increases in the three models. Quantitatively, ρH 2 increases by a factor ≈ 7 between z = 0 and z = 3 in three models, followed by a decline of a factor of ≈ 3 for the La-gos12 and Lacey14 models, and of ≈ 5 for the Gonzalez-Perez14 model between z = 3 and z = 6. At z = 0, the Lagos12 and the Gonzalez-Perez14 models predict a density of H2 that is in good agreement with the inference of Keres et al. (2003), while the Lacey14 model predicts a H2 density that is only marginally consistent with the observations. In Appendix B we show the predicted H2 mass function in the three models compared to observations of Keres et al. and show that this tension between the Lacey14 model and the observationally inferred ρH 2 rises from the lower number density of galaxies at the knee of the mass function. We warn the reader, however, that the sample of Keres et al. is not a blind CO survey and also makes use of a constant conversion factor between CO and H2. Therefore it is unclear how much of this tension is real. (see Lagos et al. 2012 for an illustration of how the CO-H2 ratio can vary.) We have also included in the middle panel of Fig. 3 recent inferences of ρH 2 from Berta et al. (2013). Berta et al. use the SFR attenuation (the ratio between the emission from the UV and from the mid-IR) to convert to molecular gas content using the empirical relation of Nordon et al. (2013). The two sets of symbols correspond to integrating the inferred H2 mass function in the range derived by the observed SFRs and integrating down to H2 masses of 10 7 M⊙, which requires extrapolation of the inferred H2 mass function. These two sets of data serve as an indication of the uncertainties in ρH 2 , although no systematic effects are included in the errorbars. The Lagos12 and Gonzalez-Perez14 models are broadly consistent with the Berta et al. inference, while the Lacey14 model lies below. In § 4.2, we describe how the physical mechanisms included in the models interplay to drive the evolution of the three quantities shown in Fig. 3.
The Lagos12, Lacey14 and Gonzalez-Perez14 models predict interesting differences, that can be tested in the future. However, with current observations of galaxies, we cannot distinguish between the models. Given that the Lagos12 model has been extensively explored in terms of the gas abundance of galaxies (see for instance Geach et al. 2011 for the evolution of the molecular gas fractions, Lagos et al. 2012 for the evolution of the CO-FIR luminosity relation and luminosity functions, and Kim et al. 2013 for an analysis of the HI mass function and clustering), we use this model in the next sections to explore the contribution from star-forming galaxies to the densities of SFR, HI and H2. In § 4.3, we explore how robust these predictions are by comparing with the Lacey14 and Gonzalez-Perez14 models. However, we remark that the other two models give similar results to the Lagos12 model in the comparisons presented in § 3.2 and 3.3.

The evolution of normal star-forming galaxies
An important test of the model is the observed relation between the SFR density at fixed stellar mass, ρSFR,M, and stellar mass . Star formation rate density as a function of stellar mass for three different redshifts, as labelled, for the Lagos12 model. Solid lines show the model predictions and dashed lines show, for reference, the effect of including an error of 0.2 dex in the stellar mass estimate (see Mitchell et al. 2013 for a complete discussion on the uncertainties in stellar mass estimates). We show the observational results from Brinchmann et al. (2004) at z ≈ 0.2, Gilbank et al. (2010) at z ≈ 1, Karim et al. (2011) at z ≈ 1 and z ≈ 2, and Sawicki (2012) at z ≈ 2. Note that we express densities and masses in comoving units for both the observations and the models. To transform these units to physical units the reader needs to divide the mass by the value of h given in § 2.3, and multiply the density by the h value squared. at a given redshift. Several works have analysed this relation and concluded that it is very difficult for current semi-analytic models and simulations of galaxy formation to predict the right trends. The argument is that this relation is very sensitive to feedback mechanisms and how inflow and outflow compensate to each other (e.g. Dutton et al. 2010;Fontanot et al. 2009;Lagos et al. 2011b;Davé et al. 2011;Weinmann et al. 2012).
Observationally, ρSFR,M is constructed using galaxies with detected SFRs (i.e. above a sensitivity limit in SFR). These galaxies are refer to as 'star-forming', and lie on a tight relationship between the stellar mass and the SFR (Schiminovich et al. 2007;Noeske et al. 2007;Elbaz et al. 2007), with outliers typically being starbursts, which have disturbed dynamics (e.g. Wuyts et al. 2011). However, in the model we include all the galaxies in a given stellar mass bin, even if they have low SFRs, to calculate ρSFR,M given that passive galaxies, which lie below the sequence of star-forming galaxies, contribute very little to the global SFR (Lagos et al. 2011b). We show the ρSFR,M-stellar mass distribution function at three different redshifts in Fig. 4. Observational data have been scaled to our choice of IMF using the scalings described in Appendix A.
There are three interesting aspects of Fig. 4. First, the Lagos12 model predicts that the highest contribution to ρSFR comes from relatively massive galaxies and that this peak mass, Ms, increases slightly with increasing redshift, from Ms ≈ 2 × 10 10 M⊙ h −1 at z = 0 to Ms ≈ 4 × 10 10 M⊙ h −1 at z = 2. The stellar mass corresponding to the peak of ρSFR,M is also consistent with the recent observational inferences of the ρSFR,M distribution function by Sobral et al. (2013). In the model, this increase of a factor of 2 is driven by the higher frequency of bright starbursts at high redshifts, which are also more frequent in massive galaxies, M stellar > 10 10 M⊙ h −1 (see Lagos et al. 2011b). Secondly, the drop in the contribution to ρSFR from galaxies with stellar masses above Ms becomes steeper with increasing redshift, which is related to the low number density of galaxies more massive than Ms at high-redshift; i.e. the break in the stellar mass function becomes smaller than Ms at z 2 (Behroozi et al. 2013). Finally, the slope of the relation between ρSFR,M and M stellar below Ms is very weakly dependent on redshift, but the normalisation increases with increasing redshift. Below Ms, most of the galaxies lie on the SFRstellar mass sequence of star-forming galaxies. Thus the slope of the ρSFR,M-M stellar is set by the slope of the SFR-stellar mass relation. Lagos et al. (2011b) show that this slope is sensitive to how quickly the gas that is expelled from the galaxy by stellar feedback reincorporates into the halo and becomes available for further cooling. Lagos et al. show that the model predicts a slope close to the observed one if this reincorporation timescale is short (i.e. of the order of a few halo dynamical timescales). The observations do not have yet the volumes and sensitivity limits needed to probe these three features, but with the current observations and the comparison with the model predictions we can conclusively say that the normalisation of the ρSFR,M-M stellar relation increases with increasing redshift, and that there is a stellar mass where the relation peaks.
Overall the model predicts a ρSFR,M-M stellar relation in reasonable agreement with the observations. The differences between the predictions of the model and the observations are in the range 0.1 − 0.3 dex. When integrating over stellar mass to give ρSFR, as presented in Fig. 3, the differences add up to make the large deviations shown in Fig. 3. Note that at z = 2, the Lagos12 model predicts a lower Ms than observed. This could be due to the fact that we are not including any error source for the stellar mass and SFRs in the model (see Marchesini et al. 2009;Mitchell et al. 2013). We show in Fig. 4 the convolution of our predicted stellar masses with a Gaussian of width 0.2 dex to illustrate the effect of this uncertainty on the predictions. Once uncertainties are taken into account, the predicted Ms agrees better with observations.

The frequency of highly star-forming galaxies
A commonly used tracer for SF is the IR luminosity, which approximates to the total luminosity emitted by interstellar dust, which, in media that are optically thick to UV radiation, is expected to correlate closely with the SFR. Infrared galaxy surveys have changed our view of the contribution from IR luminous galaxies to the IR background, or the global SFR density (e.g. Le Floc'h et al. 2005). From these surveys it is now clear that the contribution from the most luminous IR galaxies to ρSFR increases dramatically with in- creasing redshift. New mid-and far-IR surveys have allowed the characterisation of the contribution to the IR background and ρSFR from galaxies with IR luminosities down to LIR ≈ 10 10 L⊙ between z = 0 and z ≈ 1.5. These observations allowed the quantification of how more frequent bright IR galaxies are at high-redshift compared to the local Universe. We compare in Fig. 5 the predicted fractional contributions to ρSFR from galaxies with different IR luminosities in the Lagos12 model to recent observational estimates from Murphy et al. (2011), which combine measurements at 24µm and 70µm, Casey et al. (2012) using Herschel-SPIRE 250µm, 350µm and 500µm bands and Magnelli et al. (2013), which use the Herschel-PACS 70µm, 100µm and 160µm bands. The use of several IR bands ensures the detection of the peak of the IR spectral energy distribution and an accurate estimate of the total IR luminosity. Our model predicts that at z 3, the three bins of IR luminosity contribute similarly to ρSFR. At lower redshifts the situation is very different. At z < 1, ρSFR is dominated by galaxies with LIR < 10 11 L⊙, while between 1 < z < 2, ρSFR is dominated by galaxies with LIR > 10 11 L⊙. These bright IR galaxies in the Lagos12 model correspond to a combination of normal star-forming galaxies and starbursts. Note that the brightest IR galaxies make a negligible contribution to ρSFR at z < 1.5, but their contribution increases up to 34% of the total at z ≈ 3. We find that the Lagos12 model predicts fractional contributions to ρSFR that are in good agreement with the observations. Although there are systematic differences in the inferred contribution from the brightest IR galaxies to ρSFR between the observational samples, they are still consistent with each other within the errors. Magnelli et al. (2013) argue that one possible driver of such systematics is the spectroscopic redshift incompleteness, which is different for each sample.

THE CONTRIBUTION FROM STAR-FORMING GALAXIES TO THE DENSITIES OF ATOMIC AND MOLECULAR GAS
The main goal of this paper is to explore the contribution from galaxies with different properties to the overall HI and H2 content of the universe. From this we gain an insight into how far observations currently are from tracing the bulk of the neutral gas in galaxies in the universe. In this section we start by describing the contribution from star-forming galaxies to ρHI and ρH 2 and then we analyse the physical drivers behind the trends found in § 4.2.
In § 4.3 we investigate how robust these trends are by comparing the three GALFORM models (described in § 2). We will refer to the HI and H2 densities of galaxies selected by their SFR as ρHI,SFR and ρH 2 ,SFR, respectively, and by stellar mass as ρHI,m and ρH 2 ,m, respectively. We remind the reader that both gas components correspond exclusively to gas in galaxies, and that any neutral gas outside galaxies is not accounted for. As in Fig. 3, we have converted all of the observations to our adopted cosmology ( § 2.3).

Stellar mass and SFR dependence of ρHI and ρH 2
We first discuss the HI content of galaxies, then their H2 content and finish with what we expect for galaxies with different IR luminosities. HI in galaxies. The top left-hand panel of Fig. 6 shows the evolution of the contribution from galaxies with different SFRs to ρHI for the Lagos12 model. Most of the galaxies that account for the HI content of the universe have modest SFRs, SFR< 1 M⊙ yr −1 , at any time. This galaxy population alone makes up ≈ 80% of ρHI at z = 0 and ≈ 94% at z = 6. The remaining HI is found mainly in galaxies with SFRs in the range 1 < SFR/M⊙ yr −1 < 10. Galaxies with high SFRs, SFR> 10 M⊙ yr −1 , only contribute ≈ 0.7% to ρHI at z = 0 and reach a maximum contribution of ≈ 2.3% at z = 3. The form of ρHI,SFR vs. redshift for galaxies with low and high SFRs is very different. The former monotonically decreases with look-back time, while the latter have ρHI,SFR slowly increasing from z = 0 to z ≈ 2 − 3, followed by a gentle decline.
The middle left-hand panel of Fig. 6 shows the evolution of the contribution from galaxies of different stellar masses to ρHI for the Lagos12 model. All galaxies, passive and star-forming, have been included in the figure. Model galaxies with M stellar < 10 9 M⊙ dominate the HI content of the universe at z 1, while at z 1 most of the HI is locked up in galaxies with M stellar < 10 8 M⊙. Our model predicts that galaxies with M stellar > 10 10 M⊙ build up ≈ 14% of the HI content of the universe at z = 0. Schiminovich et al. (2010) and Lemonias et al. (2013) estimated ρHI,m in the observed galaxies of the GALEX Arecibo SDSS Survey (GASS) survey, which is a stellar mass selected catalogue of local Universe galaxies (see top panel of Fig. 7). For our adopted IMF, galaxies in the GASS sample have M stellar > 8.9 × 10 9 M⊙. We calculate ρHI,m of model galaxies selected using the same stellar mass cut at z = 0 and find that our prediction is a factor ≈ 0.1 − 0.2 dex lower than the observational inference. The exact factor depends on the treatment of the non-detections in the observations. This difference could be partially explained as being driven by cosmic variance; our simulated volume is relatively small, and we lack the very gas-rich galaxies with HI masses comparable to the most gas-rich galaxies included in Schiminovich et al. and Lemonias et al. to infer ρHI,m of massive galaxies. Another source of discrepancy is related to the crude treatment of stripping of HI gas in galaxy groups and clusters. In our model, when galax-  L IR total H 2 log(L IR )<10 10 10 10 <log(L IR )<10 11 10 11 <log(L IR )<10 12 log(L IR )>10 12 ies become satellites, they instantaneously lose their hot gas reservoir (i.e. strangulation of hot gas), transferring it to the central halo while slowly consuming ther remaining gas in the ISM. A more accurate treatment may be needed to explain the neutral gas content of galaxies in groups and clusters (Lagos, Davis et al. in preparation).
With the aim of making a more quantitative comparison to lower stellar mass bins, we calculate ρHI,m for galaxies with M stellar > 10 8 M⊙, M stellar > 10 9 M⊙ and M stellar > 8.9 × 10 9 M⊙ by combining the observed HI mass function of Zwaan et al. (2005) with the HI-stellar mass relation found by Huang et al. (2012) for the ALFALFA sample, corrected to a Kennicutt IMF. The Lagos12 model predictions show good agreement with these observationally inferred ρHI,m in the three mass bins.
The scatter of the observationally inferred ρHI,m is rather large due to the scatter in the HI-stellar mass relation.
Compared to previous theoretical estimates, we find that our predictions for ρHI,m agree well with the predicted trends in Davé et al. (2013). Davé et al. find that most of HI in the universe, which in their case corresponds to both HI in galaxies and in the intergalactic medium, is locked up in galaxies with stellar masses 10 8 M⊙ at z > 1.5. At z < 1.5, Davé et al. find that galaxies with stellar masses in the range 10 8 − 10 9 M⊙ increase their contribution to ρHI to similar values as the lower stellar mass galaxies. These trends are similar to the ones we find for the Lagos12 model. Popping   (2013) of carbon monoxide in normal star-forming galaxies with and without the incompleteness correction are shown as filled and empty squares, respectively. Note that we express masses in physical units. For these we corrected all of the observations to our adopted cosmology ( § 2.3). Lagos et al. (2011b) and found that the HI density is dominated by galaxies with M stellar < 10 7 , M⊙, which is even lower than in our case. These predictions from Popping et al. may be in slight tension with the observationally inferred contribution from galaxies with M stellar > 10 8 , M⊙ shown in Fig. 6, which can account for ≈ 70% of the observed HI mass density at z = 0, albeit with large errorbars. H2 in galaxies. The behaviour of star-forming galaxies in the ρH 2 -z plane, shown in the top and middle right-hand panels of Fig. 6, contrasts with that in the ρHI-z plane. Galaxies with SFR< 1 M⊙ yr −1 represent ≈ 35% of the H2 at z = 0 and only ≈ 15% at z = 3. Galaxies with 1 < SFR/M⊙ yr −1 < 10 contribute ≈ 56% of the H2 at z = 0, and their contribution also decreases with increasing redshift, yielding ≈ 40% at z = 2. One of the most interesting predictions of the top right-hand panel of Fig. 6 is that galaxies with 10 < SFR/M⊙ yr −1 < 50 contain a large fraction of the H2 of the universe at high-redshift, reaching a maximum contribution of ≈ 40% at z = 2. Highly star-forming galaxies are also more important in H2 than HI, reaching a maximum contribution of ≈ 9% in the redshift range z ≈ 3.5 − 4.5. As in ρHI, the functional form of ρH 2 ,SFR of highly star-forming galaxies is very different from that of galaxies with modest SFRs. Galaxies with SFR< 10 M⊙ yr −1 give a ρH 2 ,SFR slowly increasing with redshift up to ≈ 3.5, followed by a very slow decline at higher redshifts, while galaxies with SFR> 10 M⊙ yr −1 show ρH 2 ,SFR increasing by more than a factor ≈ 30 between z = 0 and z ≈ 2.
Observational estimates of ρH 2 ,SFR are very scarce and subject to strong systematics, such as the choice of CO-H2 conversion factor, the dust-to-gas mass ratio and the characterisation of the selection function. Even though these uncertainties limit current possibilities to derive the quantities above, it is possible to use available CO surveys of highly star-forming galaxies to infer ρH 2 ,SFR to within a factor of ≈ 2 − 3. Swinbank & et al. (2013) infer ρH 2 ,SFR for galaxies with SFR 30 M⊙ yr −1 . Swinbank et al. analyse the H2 abundance of sub-millimeter galaxies, together with three other samples of galaxies with SFRs Note that the SFR lower limit in the observed samples is approximate and can vary by a factor of ≈ 2. We compare the estimates of Swinbank et al. with our predicted ρH 2 ,SFR for these highly star-forming galaxies in the bottom panel of Fig. 7. To illustrate the effect of varying the SFR limits by a factor of 2, we show in Fig. 7 ρH 2 ,SFR for galaxies with SFR 50 M⊙ yr −1 and SFR 100 M⊙ yr −1 , which are compared to the inferences of to the CO-to-H2 conversion factor. This implies that the errorbars on the observational data in Fig. 7 are likely to be lower limits to the true uncertainty. Note that the samples of Daddi et al., Tacconi et al. and Combes et al. correspond to observations of CO, while Swinbank et al. used inferred dust masses from their FIR observations to estimate H2 masses by adopting a dust-to-gas mass ratio. The model predictions are in good agreement with the observations, within the errorbars, over the whole redshift range. Note that the data here gives independent constrains on our model than those provided by Berta et al. and shown in Fig. 3. Galaxies used to calculate ρH 2 ,SFR in the observations, as in our model predictions, are a mixture of merging systems and normal star-forming galaxies (which are referred to in our model as quiescent galaxies), which has been concluded from the measured velocity fields of the CO emission lines in some of the observed galaxies.
The trends between ρH 2 ,m and stellar mass in the middle righthand panel of Fig. 6 are again very different from those obtained for ρHI. We find that ρH 2 is dominated by relatively massive galaxies, 10 10 M⊙ < M stellar < 10 11 M⊙, from z = 0 up to z ≈ 3. Note that the contribution from these galaxies to ρH 2 increases from 42% at z = 0 to 52% at z = 1.5. At z 3, galaxies with lower stellar masses, 10 8 M⊙ < M stellar < 10 10 M⊙, dominate the H2 content. It is also interesting to note that Karim et al. (2011) and Sobral et al. (2013) find that the SFR density is dominated by galaxies with stellar masses in the range 1.4×10 10 M⊙ < M stellar < 10 11 M⊙, which overlaps with the stellar mass range we find to dominate ρH 2 . Note that the contribution from galaxies with M stellar < 10 8 M⊙ is resolved, as higher resolution runs produce the same small contribution agreeing to a factor of better than 5% (see Appendix C).
In Fig. 8, we compare the relation between the molecular gas fraction and the stellar mass in the model with observations of CO emission lines (Saintonge et al. 2011;Tacconi et al. 2013) and dust masses (Santini et al. 2013). From this figure we can draw two main conclusions. First, our model predicts a strong evolution of the molecular gas fraction that agrees well with the observationally inferred one (see also Geach et al. 2011 for a former comparison), and second, the comparison with observations needs to be done carefully, as the predicted evolution strongly depends on galaxy properties, such as the gas fraction itself, stellar mass or halo mass. Observations other than at z = 0 are strongly biased towards high molecular gas fractions, and therefore observations at z = 1.1 should be compared with the evolution predicted for galaxies with molecular gas fractions that are above 1%. Another interesting aspect is the difference between central and satellite galaxies, which comes mainly from the lack of cooled gas accretion in satellites. Lagos et al. (2011b) show that most of the galaxies in the main sequence of the SFR-M stellar plane correspond to central galaxies in the model, and therefore observed galaxies, which also fall in this main sequence, should be compared to our expectations for centrals.
The dependence on IR luminosity. In the bottom panels of Fig. 6 we show the contribution to the HI and H2 densities from galaxies with different IR luminosities. This quantity is more easily inferred from observations than the SFR or stellar masses, which, as we discussed in § 3, are affected by systematic uncertainties that are not accurately constrained. The IR luminosity, on the other hand, is well estimated from photometry in few IR bands (e.g. Elbaz et al. 2010). The trends seen with IR luminosity in Fig. 6 are similar to those obtained for the SFR (top panels of Fig. 6). This is due to the good correlation between the IR luminosity and the SFR predicted by the model. However, the normalisation of the LIR-SFR relation is predicted to vary with SFR and redshift, due to the underlying evolution of gas metallicity and the gas density in the ISM of galaxies.

The physical drivers of the gas content of star-forming galaxies
There are three key physical processes that drive the trends with SFR and stellar mass in Fig. 6; (i) the regulation between gas outflows and accretion of gas onto galaxies, (ii) the size evolution of normal star-forming galaxies (i.e. which undergo quiescent SF) and its effect on the gas surface density, and (iii) the frequency of starbursts and its evolution. Below we summarise the main effects each of these physical processes have on ρSFR, ρHI and ρH 2 . Evolution in ρ cold . The evolution of the total density of cold gas ρ cold = ρHI + ρH 2 is shown in Fig. 9. Also shown are the evolution of ρSFR and ρH 2 for all galaxies, for all quiescent galaxies, and the median half-mass radius of all galaxies with stellar masses in the range 7 × 10 9 M⊙ < M stellar < 10 10 M⊙.
Quiescent SF dominates the SFR in the universe at z 2.8, and at higher redshift it contributes ≈ 30 − 50% to ρSFR. The peak of ρSFR for quiescent galaxies takes place at lower redshift than  for starbursts. The evolution of ρSFR of quiescent galaxies is very similar to the evolution of ρH 2 , reassuring our previous statement that starbursts contribute very little to ρH 2 . The latter closely follows ρ cold at z 2. The evolution of ρ cold is driven by the balance between accretion of gas and outflows from galaxies. On average, accretion and outflows in galaxies at a given time self-regulate in a way that the SFR-M stellar relation is characterised by a small dispersion and a power-law index of ≈ 0.8 − 0.9 (Lagos et al. 2011b). Small perturbations to this self-regulation phase of galaxy growth drive the slow evolution of ρ cold and of the normalisation of the SFR-M stellar relation of normal star-forming galaxies. In the regime where ρHI ≫ ρH 2 , which is the case at z 1, ρHI follows very closely the evolution of ρ cold . At z 4, the ratio between ρH 2 and ρHI evolves slowly, and therefore ρH 2 and ρHI decrease with increasing redshift at a similar rate as ρ cold decreases. This can be seen in the evolution of the former two quantities shown in Fig. 3 and of ρ cold shown in Fig. 9. In the intermediate redshift regime, 2 < z < 4, ρH 2 ρHI, and therefore ρH 2 peaks at a similar redshift as ρ cold .
At z < 3, the evolution of ρH 2 and ρSFR of quiescent galaxies is more dramatic than it is for ρ cold . This is due to the fact that the ratio ρH 2 /ρHI for quiescent galaxies is strongly evolving due to the change in the surface density of gas plus stars, which change the hydrostatic pressure in the midplane of disks. A factor that strongly affects the evolution of the surface density of gas and stars is the evolution in the sizes of galaxies.
Evolution in size. On average, galaxy sizes increase with time (see Fig. 9). What drives the growth in the sizes of galaxy disks with time is the angular momentum of the accreted gas. As the universe expands, gas that gets to galaxies comes from further out, which means that it brings higher angular momentum compared to gas that was accreted in the past. If the angular momentum loss is small, the newly cooled gas that is accreted settles down at an outer radius, increasing the half-mass radius (see Cole et al. 2000 for a quantitative description of the size evolution of galaxies in GALFORM). We compare the evolution of sizes in the model galaxies with the observed increase reported. Note that in the case of Dahlen et al. (2007), Oesch et al. (2010) and Mosleh et al. (2011Mosleh et al. ( ,2012 there is no morphological selection and the size evolution is measured for all galaxies, while in the case of Dutton et al. (2011), the reported half-light radius evolution is for late-type galaxies only. In all cases we plot the reported evolution for galaxies with stellar masses M stellar ∼ 8.9 × 10 9 M⊙. The stellar masses of the observed samples have been scaled from a Chabrier (2003) to a Kennicutt (1983) IMF. In the model we select all late-type galaxies (i.e. bulge-to-total stellar mass ratio < 0.5) with stellar masses in the range 7 × 10 9 M⊙ < M stellar < 10 10 M⊙ and calculate the median half-mass radius of their disks. We find that the growth rate of galaxies agrees well with the observations, although with small discrepancies at higher redshifts towards larger radii. The observed samples at these high-redshifts are very small and encompass all galaxies, late-and early-type, which can bias the size-stellar mass relation toward smaller sizes. In addition, only Poisson errors have been reported which makes it difficult to estimate the significance of the deviations above. At z > 5 there are only a handful of galaxies with stellar masses in the range we plot and therefore we do not show them in this figure.
The size evolution has an effect on the gas surface density, and therefore on ρH 2 /ρHI and ρSFR. Quantitatively, ρH 2 decreases by ≈ −0.9 dex from z ≈ 2.5 to z = 0. In the same redshift range, ρ cold decreases by ≈ −0.3 dex and the median size of late-type galaxies in the model increases by 0.4 dex. From combining the evolution in ρ cold and r50 in the regime where Σ cold ≫ Σ stellar , we can derive the expected change in ρH 2 to be −1.2 dex. This number is an upper limit on the magnitude of the variations in ρH 2 between ≈ 2.5 and z = 0, as the contribution from stars to the hydrostatic pressure, on average, increases with time (see L11).
The frequency of starbursts. Fig. 10 shows ρSFR and ρH 2 for all galaxies, for all starburst galaxies, and starbursts with SFR> 50 M⊙ yr −1 . Starbursts make an important contribution to the global SFR at z 3. However, their contribution to ρH 2 is less important. The reason for this is that the SF law adopted for starbursts in our model is typically more efficient at converting gas into stars than that adopted for quiescent SF (see § 2). Starbursts with SFR> 50 M⊙ yr −1 contribute about half of the SFR and H2 mass density coming from all starbursts. The contribution from these high SFR starbursts to ρH 2 ,SFR of all galaxies, quiescent and starbursts, with SFR> 50 M⊙ yr −1 is large, and this can be seen from the similar shape and normalisation of the curve for starbursts with SFR> 50 M⊙ yr −1 in Fig. 10 and all galaxies with SFR> 50 M⊙ yr −1 in Fig. 6.
In fact, the ρH 2 ,SFR from all galaxies with SFR> 50 M⊙ yr −1 and SFR> 100 M⊙ yr −1 shown in the top right-hand panel of Fig. 6, increases by ≈ 2 − 2.3 dex from z = 0 to z ≈ 3.5. This increase is slightly larger than that of ≈ 2 dex in ρSFR and ρH 2 ,SFR of starburst galaxies with SFR> 50 M⊙ yr −1 shown in Fig. 10. The small difference between the two is due to the contribution from quiescent galaxies with SFR> 50 M⊙ yr −1 , which also have ρH 2 ,SFR increasing with increasing redshift, but more slowly compared to starbursts. Note that the strong increase of ρSFR and ρH 2 of starbursts is also largely responsible for the evolution of galaxies with LIR > 10 12 L⊙ in Fig. 5. The contribution from starbursts to ρHI is negligible.

The HI and H2 gas contents of galaxies in the Lacey14 and Gonzalez-Perez14 models
One of the aims we have when analysing different GALFORM models is to get a better insight into how much variation we expect in the predicted contributions from galaxies selected according to their stellar mass or SFR to the densities of HI and H2. This also helps to illustrate how some physical processes, other than the SF law, affect the predictions described in § 4. From this, it is also possible to identify how robust the trends described in § 4 are. Fig. 11 shows the fractional contribution from galaxies in two different bins of SFR to the densities of HI, H2 and SFR in the Lagos12, Lacey14 and Gonzalez-Perez14 models. The bins were chosen to match those of Fig. 6. We regard the following trends as robust as they are present in all three models: (i) galaxies with low SFRs dominate ρHI, while for ρH 2 and ρSFR, they make a modest contribution that decreases from z = 0 to z ≈ 2.5, (ii) galaxies with high SFRs have a contribution to ρHI that increases with increasing redshift, at least up to z ≈ 2.5, but that never reach more than few percent, for ρH 2 . Their contribution is of ≈ 8 − 15% at z ∼ 0, increasing to ≈ 25 − 40% at z ≈ 1 − 2.5, and for ρSFR their contribution increases with redshift up to z ≈ 1.5 and stays relatively constant at higher redshifts. Note that the exact peak and normalisation of the contributions of each model sample to the densities of HI, H2 and SFR depend on the model. Fig. 12 shows the fractional contribution from galaxies in two different bins of stellar mass to ρHI, ρH 2 and ρSFR in the Lagos12, Lacey14 and Gonzalez-Perez14 models. The bins were chosen to match those of Fig. 6. Here, the following trends are robust: (i) ρHI has always a large contribution from low mass galaxies that increases with increasing redshift; for H2, this contribution is negligible at z 2, but becomes systematically more important at higher redshifts, and for ρSFR they contribute less than 10% at z 3.5, increasing slightly their contribution at higher redshifts; (ii) more massive galaxies (M stellar > 10 10 M⊙) are never important contributors to ρHI, but host most of the H2 in the universe at z 3; in the case of ρSFR, they make more than 30% of the total SFR at z 3.5, percentage that slowly decreases at higher redshifts. The exact normalisations of ρHI,m and ρH 2 ,m slightly vary between the three different models. The largest differences obtained for the contribution to ρHI, ρH 2 and ρSFR from different selections in SFR are seen at z 4.

IMPLICATIONS FOR THE NEXT GENERATION OF GALAXY SURVEYS
In this section we discuss the implications of the results presented in § 4 for the current and next generation of galaxy surveys that aim to detect gas in galaxies. We first focus on the atomic content of galaxies in § 5.1 and then on the molecular gas in § 5.2.

Atomic hydrogen in high-redshift galaxies
There are two ways of studying HI, (i) through its emission at 21 cm, either by individually detecting galaxies or stacking them to reach a high signal-to-noise ratio and (ii) through the absorption of the light of a background object by the presence of neutral hydrogen in the line-of-sight. We discuss future surveys in the context of HI in emission and absorption. HI in emission. Large blind surveys of the 21 cm emission will be carried out by the next generation of radio telescopes, such as the APERITIF project at the Westerbork Synthesis Radio Telescope, the Australian SKA Pathfinder 6 (ASKAP), the South-African SKA pathfinder MeerKAT 7 (MeerKAT), the Jansky Very Large Array 8 (JVLA) and in the longer term future, the Square Kilometer Array 9 (SKA). One of the key surveys planned with ASKAP is DINGO, in which galaxies up to z ≈ 0.43, with the key science goal of understanding how the HI content of galaxies co-evolves with other major galactic components, such as stars, dust and halo mass. WAL-LABY will cover much larger areas in the sky, allowing the understanding of the effect of environment on the HI content of galaxies. At higher redshifts, the Looking At the Distant Universe with the MeerKAT Array survey 10 (LADUMA) plans to detect HI emission in individual galaxies up to z ≈ 1.
Although these future surveys will provide unprecedented knowledge of the neutral gas in galaxies and its cycle, an important fraction of the HI in the universe is expected to be hosted in galaxies with small HI contents (MHI < 10 6 M⊙) due to the steep low-mass end of the HI mass function. These masses may even be too small to be detected by ASKAP and MeerKAT. To overcome this problem and estimate the HI mass abundances in galaxies that are not individually detected, it is possible to use intensity mapping or stacking of the 21 cm emission line (see Pritchard & Loeb 2011 for a review). The main difference in the estimates that can be made with intensity mapping is that what it is inferred is ρHI · bHI (e.g. Chang et al. 2010;Masui et al. 2013), where bHI is the bias in the clustering of HI selected sources with respect to the expected underlying dark matter distribution, while stacking can infer a value Table 1. The fraction of galaxies in the Lagos12 model that have HI velocity-integrated fluxes larger than 3 times the 1σ sensitivity limits, S HI,lim (column 1), of 4 different surveys discussed in § 5.1, and that have stellar masses < 10 8 M ⊙ , SFR< 0.1 M ⊙ yr −1 or IR luminosities < 10 9 L ⊙ , for different redshifts that cover the redshifts explored by the surveys (columns 2 to 6). The numbers in each column are the fraction of galaxies fulfilling the condition of a > 3σ detection at the given redshift, given the flux limit, [M stellar < M ⋆,lim ; SFR<SFR lim ; L IR < L IR,lim ], where M ⋆,lim = 10 8 M ⊙ , SFR lim = 0.1 M ⊙ yr −1 and L IR,lim = 10 9 L ⊙ . The sensitivity limits of the ASKAP surveys given in (1) were calculated by Duffy et al. (2012b) for a bandwith of 3.86 km s −1 , and for LADUMA by Sarah Blyth (priv. communication) for a bandwidth of 3.2 km s −1 . (1) (2) (3) (4) (5) 1σ sensitivity in mJy km s −1 z = 0.08 z = 0.17 z = 0.5 z = 0. for ρHI from a parent galaxy sample (e.g. Delhaize et al. 2013;Rhee et al. 2013).
Another additional problem that these future surveys may encounter is the identification of UV, optical and/or IR counterparts. We argue that although to measure the total HI density at a given redshift is an important quantity for galaxy formation theory and cosmology, to be able to measure HI in galaxy populations which are selected by other means, either optical emission lines, broadband photometry, etc., is also valuable as additional constraints that can be put on galaxy formation models. However, even if a galaxy is detected in HI emission, its stellar mass and/or SFR may be too small to be measured, making the task of learning about galaxy formation from these HI surveys more challenging.
HI in absorption. Another solution to the problem of faint HI emission is to observe neutral hydrogen in absorption. This has been extensively used to measure the number density of absorbers of different classes, such as DLAs, Lyman limit systems, etc. (see for example Péroux et al. 2003). The first blind absorption survey in HI is going to be carried out by ASKAP 11 (FLASH). This survey will allow measurements of ρHI as well as constraints on the evolution of faint HI galaxies by combining the FLASH survey with optical galaxy surveys as discussed above, in the redshift range 0.5 < z < 1. Recent HI absorption surveys have started to sample statistically the difference between absorbers around late-and early-type galaxies, as well as galaxies in different environments (Tejos et al. 2013;Tumlinson et al. 2013).
We expect that only a combination of the surveys above (direct detection, absorption studies and stacking) will be able to account for most of the HI in galaxies and determine the exact contribution to ρHI from galaxies with different properties. We summarise in Table 1 the fraction of galaxies that have low stellar masses, SFRs or IR luminosities but that have HI velocity-integrated fluxes that would be detected at a > 3σ level by the deep HI surveys we described above. These numbers are intended to show the predicted frequency of HI detections that are likely to have no optical counterpart in large optical and near-IR galaxy surveys. A possible solution to this problem is to follow up these HI sources with deep optical and near-IR observations. 11 http://www.physics.usyd.edu.au/sifa/Main/FLASH

Molecular hydrogen in high-redshift galaxies
To trace molecular gas, lower order molecule transitions are better, as they are most commonly excited predominantly through collisions (in optically thick media), and they are more easily excited than the higher order molecule transitions. The most used molecule to trace H2 is CO, which is the most abundant molecule after H2. However, the condition of optically thick media, that defines the excitation means, may break-down in very low metallicity gas, where the CO lines become optically thin (see  for a recent review). Lagos12 find that the number of galaxies with gas metallicities Zgas < 0.01 Z⊙ is < 0.05% for galaxies with IR luminosities > 10 9 L⊙ in the redshift range 0 z 6. This IR luminosity limit corresponds roughly to a SFR∼ 0.1 M⊙ yr −1 . Thus, our model indicates that CO can be used as a good tracer of H2 for the galaxies that dominate the H2 content of the universe, which are those with SFR> 1 M⊙ yr −1 or with M stellar > 10 10 M⊙. Galaxies in which CO is not a good tracer of H2 (i.e. those with low gas metallicities) represent only a small contribution to the H2 content of the universe. The downside of CO is the CO-H2 conversion factor, which is a strong function of metallicity and other physical conditions in the ISM of galaxies (Boselli et al. 2002). The positive side is that the CO-H2 conversion factor has started to be systematically explored in cosmological simulations of galaxy formation with the aim of informing how much variation with redshift is expected (Narayanan et al. 2012;Lagos et al. 2012;Popping et al. 2013b).
Blind molecular emission surveys are the best way of studying the cold gas content of galaxies in an unbiased way. However, they require significant time integrations to ensure good signal-to-noise. Alternatively, molecular emission line surveys can be guided by well-known UV, optical, IR and/or radio surveys, with well identified positions and spectroscopic or photometric redshifts. Even if only highly star-forming galaxies are used as guide for follow-up molecular emission line surveys, it is possible to place strong constraints on models of galaxy formation, as we show in § 4, and to provide new insight into the co-evolution of the different baryonic components of galaxies.
The Atacama Large Millimeter Observatory (ALMA) can currently detect CO(1-0) and CO(2-1) up to z ≈ 0.84. At z 1.3, instruments that are primarily designed to detect atomic hydrogen, such as the JVLA are able to detect the emission from the redshifted frequencies of low-J CO transitions.A good example of this is presented in Aravena et al. (2012), who performed a JVLA blind survey in a candidate cluster at z ∼ 1.5. Aravena et al. obtained two Table 2. The fraction of galaxies in the Lagos12 model that have CO emission line velocity-integrated fluxes in the rotational transitions from 1-0 to 6-5 (column 1) above a 3σ flux, with 1σ = 1 mJy km s −1 , and that have stellar masses < 10 9 M ⊙ , SFR< 10 M ⊙ yr −1 or IR luminosities < 10 10 L ⊙ , for different redshifts up to z = 2.8 (columns 2 to 7). The numbers in each column are the fraction of galaxies at the given redshift fulfilling the conditions [M stellar < M ⋆,lim ; SFR<SFR lim ; L IR < L IR,lim ], where M ⋆,lim = 10 9 M ⊙ , SFR lim = 10 M ⊙ yr −1 and L IR,lim = 10 10 L ⊙ . (1) (2) (3) (4) (5) (6) CO transition z = 0.1 z = 0.5 z = 1 z = 1.5 z = 2 z = 2.8 detections and were able to place constraints on the number density of bright CO(1 − 0) galaxies. Similarly, blind CO surveys in small areas of the sky are ongoing using the JVLA and the Plateau de Bureau Interferometer . For faint sources, the intensity mapping technique described in § 5.1 can also be applied to molecular emission. The line confusion from galaxies at different redshifts can be overcome by crosscorrelating different molecular and/or atomic emission lines. Different emission lines coming from the same galaxies would show a strong correlation, while emission lines coming from galaxies at different redshift would not be correlated. From this, the quantity that is derived is the two-point correlation function or power spectrum weighted by the total emission in the spectral lines being correlated (see Pritchard & Loeb 2011 for a review). The power of using intensity mapping for CO lines have been shown and discussed for instance by Carilli (2011) and Gong et al. (2011). Carilli argues that intensity mapping of CO lines can be performed with telescopes of small field-of-view at redshifts close or at the epoch of reionisation (z ∼ 3 − 8).
Unlike in HI, optical counterparts could be easily identified in the case of individual detection of molecular emission, since they are expected to be bright and have large stellar masses. Table 2 shows the fraction of galaxies in model samples selected by their CO velocity-integrated flux and that have stellar masses < 10 9 M⊙, SFR< 10 M⊙ yr −1 or IR luminosities < 10 10 L⊙. These fractions are encouraging as they are generally low compared to HI mass selected samples, particularly at high redshifts; this means that we expect effectively all of the molecular emission lines detected to have detectable optical/IR counterparts. From this, one can conclude that the most efficient strategy to uncover the H2 in the universe is the follow-up of already existing optical/IR/radio catalogues with submillimeter and/or radio telescopes with the aim of detecting molecular emission lines.

CONCLUSIONS
We have presented predictions for the contribution to the densities of atomic and molecular hydrogen from galaxies with different stellar masses, SFRs and IR luminosities in the context of galaxy formation in a ΛCDM framework. We use three flavours of the GALFORM semi-analytic model of galaxy formation, the Lagos12, Gonzalez-Perez14 and Lacey14 models. For quiescent star formation, the three models use the pressure-based SF law of Blitz & Rosolowsky (2006), in which the ratio between the gas surface density of H2 and HI is derived from the radial profile of the hydrostatic pressure of the disk, and calculates the SFR from the surface density of H2. The advantage of this SF law is that the atomic and molecular gas phases of the ISM of galaxies are explicitly distinguished. Other physical processes in the three models are different, such as the adopted IMF and the strength of both the SNe and the AGN feedback, as well as the cosmological parameters adopted.
We identify the following trends in the three models and regard them as robust: (i) The HI density shows almost no evolution at z 2, and slowly decreases with increasing redshift for z 2.
(ii) The H2 and the SFR densities increase strongly from z = 0 up to z ≈ 2 − 3, followed by a slow decline at higher redshifts. The models predict an increase in SFR larger than for H2 due to the contribution from starbursts, which is larger for the SFR (close to 50% in the three models at z ≈ 3) than for H2 (less than 10% at z ≈ 3). The latter is due to starbursts being on average more efficient at converting gas into stars, and therefore contributing more to the SFR density than the gas density of the universe.
(iii) The three models predict densities of SFR, HI and H2 in broad agreement with the observations at z 3. A higher redshifts important deviations are observed between the predicted ρHI and the observationally inferred values in the three models. We argue that this is not unexpected as the HI in our simulations correspond exclusively to HI in galaxies and no contribution from the HI in the IGM is included. The latter is expected to become important at z 3 (van de Voort et al. 2012).
(iv) The density of HI is always dominated by galaxies with low SFR (SFR< 1 M⊙ yr −1 ), and low stellar masses (M stellar < 10 9 M⊙), while the H2 density is dominated by galaxies with relatively large SFRs (SFR> 1 M⊙ yr −1 ) and large stellar masses (M stellar > 10 10 M⊙). The latter is also true for the global SFR density.
(v) The predicted evolution of the global SFR density observed in the universe can be largely explained as driven by the steep decline of the molecular mass towards z = 0. The combination of the evolution of the total neutral gas, atomic and molecular, and the increasing galaxy sizes with decreasing redshift explain the evolution of the H2 density. The global SFR density evolution can therefore be linked to the evolution of the neutral gas surface den-sity of the galaxies dominating the SFR in the Universe at a given time. The evolution of the neutral gas content of galaxies is set by the balance between inflows and outflows of gas in galaxies, and therefore also plays a key role in the evolution of the molecular gas, that ultimately set the SFR. Thus, one must understand that the H2 content and SFR of galaxies in the star-forming sequence of the SFR-stellar mass plane are set by the self-regulation of inflow and outflows and that small deviations to this self-regulation produce the time evolution of this sequence.
(vi) We find that a group of observations can be connected and understood in the models, including the Hα and UV luminosity functions, distribution function of the global SFR density on stellar mass, the global SFR density from highly star-forming galaxies, the inferred HI density from massive galaxies in the local Universe, and the inferred H2 density of highly star-forming galaxies at highredshifts. This broad set of comparisons between observations and model predictions that test at the same time the stellar mass and SFR and the neutral gas content of galaxies is unprecedented in semi-analytic models of galaxy formation.
We discuss our findings in the context of future surveys and suggest that optical identification of faint HI counterparts is going to be difficult due to the low stellar masses and SFRs we expect for them. These faint HI galaxies are, however, key to our understanding of galaxy formation as we expect them to dominate the HI content of the universe. We suggest that stacking techniques are promising to measure the contribution to HI from galaxies selected by their stellar mass or SF activity. Such measurements have the potential to place strong constraints on galaxy formation models. Regarding H2, the easiest way to make progress is to use current surveys with well identified optical, IR and/or radio positions and well known spectroscopic or photometric redshifts and follow them up with current facilities, with the downside that the selection function becomes more complex to model. The fact that most of the H2 is predicted by our model to be hosted by relatively massive, starforming galaxies, which are already detectable by optical and IR surveys, it is implied that follow-up surveys will be able to uncover the H2 density of the universe without the need of blind surveys. From such observations, it is implied that both direct detection and stacking would provide valuable information on the cold gas content of galaxies and place strong constraints on galaxy formation models.  Kennicutt (1983) IMF and M IMF stellar the stellar mass inferred in observations that adopted a different IMF. These are taken from the work of Bell et al. (2003) and Gilbank et al. (2010).
The case of SFRs is more complex than the stellar mass, since the conversion depends on the SFR tracer used. We here follow the scalings calculated by Gonzalez-Perez et al. (2013, in prep.). Gonzalez-Perez et al. use the PEGASE (Fioc & Rocca-Volmerange 1999), and assume a constant solar etallicity and SFR to calculate the Hα, Hβ, [OII], FUV and radio luminosities for each choice of IMF. Note, however, that the conversions are subject to systematics related to the assumptions of SF history and metallicity (e.g. Wilkins et al. 2012). We summarise in Table A2 the scalings applied to each observationally inferred SFR depending on the adopted IMF and used tracer to convert to a Kennicutt (1983) IMF. Here, FUV corresponds to a filter centered at 1500Å and width 400Å and NUV to a filter centered at 2000Å and width 400Å.

APPENDIX B: THE LOCAL UNIVERSE HI AND H2 MASS FUNCTIONS
In Fig. B1 we show the predicted HI and H2 mass functions at z = 0 for the Lagos12, Lacey14 and Gonzalez-Perez14 models. We compare with observations of the HI mass function from Zwaan et al. (2005) and Martin et al. (2010), and of the inferred H2 mass function from Keres et al. (2003). The latter observations correspond to CO(1 − 0) molecular emission line in galaxies that were previously selected from their emission in the B-band or 60µm. We then use a Milky-Way type conversion between the CO(1−0) emission and the H2 mass, and infer a H2 mass function. We showed in Lagos et al. (2012) that doing the opposite exercise of converting predicted H2 masses to CO(1−0) emission using radiative transfer, photon dominated region codes to calculate the CO-H2 conversion factor, leads to similar CO(1 − 0) luminosity function than the one obtained by assuming a constant conversion factor. This is due to the fact that most model galaxies at z = 0, which have CO(1 − 0) luminosities above the observational limit of Keres et al. (2003), have gas metallicities and SFR densities that are not too different from normal spiral galaxies. The three models give predictions that are in good agreement with the observations.

C1 The Lagos12, Lacey14 and Gonzalez-Perez14 models in the Millennium database
The predictions of the three GALFORM models used in this paper will be publicly available from the Millennium database. However, the reader should keep in mind that the runs performed in this paper are not 'standard' in the sense that we use Monte-Carlo realisations of a wider dynamical range of halo masses compared to the N -body simulations Millennium (Springel et al. 2005), used by the Lagos12 model, and MS-W7 (same as the Millennium simulation but with WMAP7 cosmology), used by Gonzalez-Perez14 and Lacey14. This means that the runs shown in this paper follow the formation of galaxies in dark matter halos that extend to much lower masses than the two N -body simulations above. The minimum halo mass of the Millennium and MS-W7 simulations is fixed on 1.72×10 10 h −1 M⊙. This has the limitation that the results presented here for the HI density cannot be fully reproduced with the galaxies in the Millennium database. This is due to the faint nature of the galaxies that dominate ρHI, which are hosted by dark matter halos that are not resolved in the Millennium and MS-W7 N -body simulations. All the results related to H2 and SFR density (including the Hα and UV luminosity functions) can be reproduced using the galaxies in the Millennium database at z < 4. This is because the galaxies that dominate these statistical measurements are relatively massive, and hosted by halos that are well resolved in the N -body simulations. At z > 4, the SFR density is not fully resolved in the resolution of the Millennium and MS-W7 N -body simulations.
We have also performed runs using the dark matter halo merger trees of the Millennium-II simulation, which has a minimum halo mass fixed at 1.38×10 8 h −1 M⊙ (Boylan-Kolchin et al. 2009). We find that the dynamical range of the Millennium-II is sufficient to resolve the HI in galaxies at z 2.5. However, at higher redshifts, the Millennium-II misses the lowest halo masses that still contribute to the HI in the universe (10 7 − 10 8 h −1 M⊙), Here we show observational inferences from Keres et al. (2003) using direct detection of the CO(1 − 0) emission line in B-band (triangles) and 60 µm (filled circles) selected samples of galaxies. The observational mass function is then calculated assuming the Milky Way H 2 -to-CO conversion factor, N H 2 /cm −2 = 2 × 10 −20 I CO /K km s −1 . Here N H 2 is the column density of H 2 and I CO is the integrated CO(1 − 0) line intensity per unit surface area.
that are covered by our Monte-Carlo merger trees thanks to the redshift scaling applied to the minimum halo mass.

C2 Resolution tests
A questions that comes immediately from the lack of resolution in the Millennium and MS-W7 N -body simulations is whether the results presented here for the HI density are converged or not. To answer this question, we performed a simple exercise: we run the same model under 6 different resolutions, with and without the scaling with redshift of the minimum halo mass. Table C1 summarises the different runs. The scaling with redshift is also changed Table C1. Simulations used to test the convergence of the predictions of the models. Column (1) gives the name of each simulation, (2) the minimum halo mass adopted, M halo,min , and (3) the scaling of the minimum halo mass with redshift. The standard simulation used throughout the manuscript is r3.
Simulation M halo,min /M ⊙ h −1 Scaling r1 10 10 (1 + z) −3 r2 5 × 10 9 (1 + z) −3 r3 5 × 10 8 (1 + z) −3 r4 2.5 × 10 8 (1 + z) −3 r5 10 10 no scaling r6 5 × 10 9 no scaling to study the effect of this scaling on the HI density to reproduce the conditions of the N -body runs (i.e. fixed minimum halo mass) and also to show the importance of adding the redshift scaling in the Monte-Carlo merger trees. Fig. C1 shows a visualization of the different resolutions of all the runs. We also show in Fig. C1 the resolution of the Millennium I and Millennium II N -body simulations. We show in Fig. C2 the global HI density and the contribution to it from galaxies with M stellar < 10 8 M⊙ for the Lagos12 model in the five Monte-Carlo tree resolutions, all normalised to our standard choice, the 'r3' Monte-Carlo merger trees. By comparing the r3 with the r4 resolution, we can conclude that our standard resolution (r3) is converged to a factor better than 2%. When comparing the resolutions r1 and r2 with our standard r3 resolution, we conclude that although the global HI density is converged to a factor better than 10%, the contribution from low mass galaxies is not. This is shown in the bottom panel of Fig. C2, where r1 and r2 predict lower contributions from the low mass galaxies to ρHI by ≈ 60% and ≈ 10%, respectively. The importance of the redshift scaling introduced in the minimum halo mass in our stan-   Figure C2. The ratio between the HI density in a given resolution and the r3 run, as a function of redshift. The top panel shows the global HI density ratio and the bottom panel shows the HI density in galaxies with M stellar < 10 8 M ⊙ in the Lagos12 model. Each line shows a different resolution test as labelled with respect to our standard run, the r3 resolution. dard run is crucial to make our predictions converge. This can be concluded from the large offset between the r3 and the runs r5 and r6, which increases with increasing redshift. Lagos et al. (2011a) argue that this is the reason why previous models predict a monotonically increasing H2/HI global density ratio with increasing redshift (e.g. Popping et al. 2013a), and miss the turnover in this ratio, which takes place at z ≈ 5 (i.e. in the regime where the fixed resolution run largely underestimates the global HI density). Also note that the runs r5 ad r6 predict a ρHI that converges to the high-resolution ρHI at z 1. The contribution from low mass galaxies (M stellar < 10 8 M⊙) to the density of H2, shown in Fig. C3 is converged to better than 10% at the resolution of r2, meaning that the contribution from small galaxies is resolved in or standard simulation. The fact that the convergence is achieved at lower resolutions compared to the case of HI is due to the flatter low-mass end of the H2 mass function compared to the steepness of the low-mass end of the HI mass function.