11institutetext:  INAF-Italian National Institute of Astrophysics, Observatory
of Trieste, via G. Tiepolo 11, 34143 Trieste, Italy 22institutetext: Max
Planck Institute for Astrophysics, Karl-Schwarzschild-Str. 1, 85748 Garching
bei München, Germany 33institutetext: European Southern Observatory, Karl-
Schwarzschild-Str. 2, 85748 Garching bei München, Germany; 44institutetext:
Aix Marseille Université, CNRS, LAM (Laboratoire d’Astrophysique de Marseille)
UMR 7326, 13388, Marseille, France

# Atomic and molecular gas from the epoch of
reionization down to redshift 2

Umberto Maio  Céline Péroux 1122 Benedetta Ciardi 334422

###### Abstract

Context. Cosmic gas makes up about 90% of baryonic matter in the Universe and
H2 is the closest molecule to star formation.

Aims. In this work we study cold neutral gas and its H2 component at different
epochs, exploiting state-of-the-art hydrodynamic simulations that include
time-dependent atomic and molecular non-equilibrium chemistry coupled to star
formation, feedback effects, different UV backgrounds presented in the recent
literature and a number of additional processes occurring during structure
formation (ColdSIM).

Methods. We predict gas evolution and contrast the mass density parameters and
gas depletion timescales, as well as their relation to cosmic expansion, in
light of the latest IR and (sub-)mm observations in the redshift range
$2\lesssim z\lesssim 7$.

Results. By performing updated non-equilibrium chemistry calculations we are
able to broadly reproduce the latest HI and H2 observations. We find neutral-
gas mass density parameters $\Omega_{\rm neutral}\simeq 10^{-3}$ and
increasing from lower to higher redshift, in agreement with available HI data.
Because of the typically low metallicities during the epoch of reionization,
time-dependent H2 formation is mainly led by the H- channel in self-shielded
gas, while H2 grain catalysis becomes important in locally enriched sites at
any redshift. UV radiation provides free electrons and facilitates H2 build-up
while heating cold metal-poor environments. Resulting H2 fractions can be as
high as $\sim 50\%$ of the cold gas mass at $z\sim 4$-8, in line with the
latest measurements from high-redshift galaxies. The H2 mass density parameter
increases with time until a plateau of $\Omega_{\rm H_{2}}\simeq 10^{-4}$ is
reached. Quantitatively, we find an agreement between the derived $\Omega_{\rm
H_{2}}$ values and the observations up to $z\sim 7$ and both HI and H2 trends
are better reproduced by our non-equilibrium H2-based star formation
modelling. The predicted gas depletion timescales decrease at lower $z$ in the
whole time interval considered, with H2 depletion times remaining below the
Hubble time and comparable to the dynamical time at all $z$. This implies that
non-equilibrium molecular cooling is efficient at driving cold-gas collapse in
a broad variety of environments and since the very early cosmic epochs. While
the evolution of chemical species is clearly affected by the details of the UV
background and gas self-shielding, the assumptions on the adopted initial mass
function, different parameterizations of H2 dust grain catalysis,
photoelectric heating and cosmic-ray heating can affect the results in a non-
trivial way. In appendix, we show detailed analyses of individual processes,
as well as simple numerical parameterizations and fits to account for them.

Conclusions. Our findings suggest that, in addition to HI, non-equilibrium H2
observations are pivotal probes for assessing cold-gas cosmic abundances and
the role of UV background radiation at different epochs.

###### Key Words.:

Cosmology: theory - structure formation; cosmic gas

††offprints: U. Maio, e-mail<EMAIL_ADDRESS>
##  1 Introduction

About 90% of the baryons in the Universe are in the form of gas and most of
this resides either in the cosmic space, constituting the so-called inter-
galactic medium (IGM), or around galaxies, where it shapes the circum-galactic
medium (CGM). Only a tiny, but extremely relevant, fraction of gas populates
neutral high-density regions, where the local medium cools and fragments.
There, a large portion of the chemical composition is represented by neutral
hydrogen, HI, and molecular hydrogen, H2, the most abundant molecule in the
Universe. Different gas phases are not necessarily separated, as structure
formation brings diffuse material into dense molecular regions, as well as
spreads initially dense parcels into the hotter CGM/IGM. How the involved
processes work in detail is still matter of debate and the interplay between
different phases has drastic consequences for the ionization state of H and
the global amount of H2.
In the whole $z<6$ epoch, the measured neutral-gas mass density normalized to
the present-day critical density (neutral-gas density parameter) is
$\Omega_{\rm neutral}\sim 10^{-3}$ and increases towards higher redshifts
Noterdaeme et al. (2009); Zafar et al. (2013); Rhee et al. (2013); Sánchez et
al. (2016); Jones et al. (2018); Hu et al. (2019); Wolz et al. (2021); Chen et
al. (2021). H2 molecules have a more complex behaviour and are notoriously
more difficult to probe. Until a couple of years ago there was little reliable
information at intermediate redshifts and no data at all about cold-gas H2
content in the distant Universe. Only recently first data from a number of
observational campaigns based on e.g. VLA,111http://www.vla.nrao.edu/
ALMA222https://www.almaobservatory.org, NOEMA333https://www.iram-
institute.org/EN/noema-project.php or UKIRT444http://wsa.roe.ac.uk/UKIRT.html
have started to set constraints on its cosmological evolution.
The exact processes driving molecular build-up in different cosmological
environments are not fully understood, yet. For this reason, H2 mass in
targeted objects is usually inferred through empirical fits, based, for
instance, on emission from CO Scoville & Solomon (1975); Dickman et al.
(1986); Maloney & Black (1988); Hughes et al. (2017) or [C II] Madden et al.
(1997); Pak et al. (1998); Hollenbach & Tielens (1999) calibrated in the local
Universe. Heavy elements and C-rich molecules have been detected in metal-
enriched collapsing sites, suggesting that CO cores exist surrounded by a [C
II] envelope where H2 is self-shielded from photodissociating radiation Röllig
et al. (2006); Wolfire et al. (2010); Genzel et al. (2012); Bolatto et al.
(2013); Zanella et al. (2018); Madden et al. (2020); Gong et al. (2020).
Notwithstanding the uncertainties related to calibration techniques,
completeness and conversion factors, the currently estimated H2 molecular-mass
density in the local Universe is $\rho_{\rm H_{2}}\simeq 6.8\times 10^{6}\,\rm
M_{\odot}/Mpc^{3}$ and accounts for about 20 per cent of the total abundance
of cold gas Fletcher et al. (2021). Despite the large scatter of one order of
magnitude or more among different studies, the H2 mass density parameter
$\Omega_{\rm H_{2}}$ is observed to increase at $z\simeq 0$-3, assembling its
mass with practically no observed environmental dependence Darvish et al.
(2018); Tadaki et al. (2019); Garratt et al. (2021), and to mimick the
cosmological star formation rate (SFR) density at higher $z$ Decarli et al.
(2020); Riechers et al. (2020a). These observations have stressed the
longstanding need to explore the origin of the large supply of H2 gas required
to sustain star formation at different $z$ Rodighiero et al. (2019); Tacconi
et al. (2020); Hunt et al. (2020) and questioned the ability to explain the
large H2 fractions (above 50%) detected in galaxies around $z\sim 4$-6
Dessauges-Zavadsky et al. (2020); Boogaard et al. (2021), in a primordial,
presumably metal-poor, gas.
From a theoretical point of view, cold gas plays a key role for cooling and
fragmentation. It hosts the conversion from atomic to molecular phase, through
which it fuels collapse and star formation (Saslaw & Zipoy 1967, Peebles &
Dicke 1968). Despite the well-determined HI observational trend, a full
understanding of its physical behaviour over cosmological history still needs
to be consolidated. H-derived molecules are crucial for gas chemistry
evolution (see e.g. Bromm & Yoshida 2011; Girichidis et al. 2020, for
reviews), because they supply low-temperature coolants to cosmic gas and are
the main drivers of gas depletion timescales. A full assessment of H2
evolution in different epochs is elusive, because of both the lack of
observational data at our disposal until recently and the large computing
power needed to perform detailed chemistry calculations within three-
dimensional simulations. A correct evaluation of cold-gas fractions is
sensitive to a number of cooling and heating processes due to either small-
scale physics or large-scale events Maio (2021). Numerical analyses typically
predict lower-than-observed neutral/HI mass densities and have difficulties in
reproducing their increasing trend in redshift at $z\simeq 2$-6. Accurate
chemistry calculations following H2 formation and the consequent run-away
cooling regime require a full time-dependent ”non-equilibrium” approach that
is computationally expensive and hence not widespread within numerical
implementations of cosmological simulations. It has so far been employed to
study mainly primordial epochs Abel et al. (1997); Yoshida et al. (2003); Maio
et al. (2010); Wise et al. (2012b); Muratov et al. (2013); Skinner & Wise
(2020), while non-equilibrium evaluations for $\Omega_{\rm neutral}$,
$\Omega_{\rm H_{2}}$ and depletion times during and after the epoch of
reionization are still missing.
The main channels for H2 production in pristine media, composed mainly by $\rm
H$, He and electrons, are based on catalysis of intermediate species
interacting with H, with possible photon emission accompanying recombination
Galli & Palla (1998), i.e. H- catalysis ($\rm H+e^{-}\rightarrow H^{-}+\gamma$
and $\rm H^{-}+H\rightarrow H_{2}+e^{-}$) and H${}_{2}^{+}$ catalysis ($\rm
H+H^{+}\rightarrow H_{2}^{+}+\gamma$ and $\rm H_{2}^{+}+H\rightarrow
H_{2}+H^{+}$). While these processes are dominant at intermediate densities
below $10^{4}\,\rm K$ or in recombining gas where free charges are available,
three-body interactions ($\rm 3H\rightarrow H_{2}+H$) could further contribute
in dense neutral media.
For metal-enriched environments the picture is even more complex, since
metallicity ($Z$) evolution and dust grain catalysis might influence H2
abundances ($\rm 2H\rightarrow H_{2}$ on dust grains), also depending on the
efficiency of photoelectric heating (e.g. Bakes & Tielens 1994), whose
implications on H2 production at different cosmological epochs are currently
uncertain. Cosmic-ray protons losing their energy in dense molecular gas
Goldsmith & Langer (1978) or in thinner photodissociating regions Shaw et al.
(2009) could heat the medium and affect early H2 formation in a non-trivial
way. Their possible role for pristine primordial gas has been pointed out with
basic semi-analytic calculations by Jasche et al. (2007). However, also
because of the poor information on cosmic-ray heating yield and rate at that
time, their impact has not been simulated within three-dimensional numerical
simulations of early structure formation and the implications on (non-
equilibrium) HI and H2 masses are still unknown.
Feedback mechanisms can remove local gas mass and inject entropy into the
surrounding medium (e.g. Nagamine et al. 2004). As UV photons are capable to
dissociate molecules and ionize neutral atoms in non-shielded regions, the
build-up of a UV background determined by forming structures in the first Gyr,
at redshift $z\gtrsim 6$, affects the gas thermal state, as well as HI and H2
density parameters. Gas self-shielding (e.g. Draine & Bertoldi 1996; Rahmati
et al. 2013; Hartwig et al. 2015; Wolcott-Green et al. 2017; Ploeckinger &
Schaye 2020; Luo et al. 2020; Walter et al. 2020), limits the ability of UV
radiation to penetrate dense material and preserves HI or H2 fractions in star
forming regions. The relevance and the relative role of the two main (HI and
H2) shielding processes for chemical abundances in cosmic environments need to
be further explored, though.
In addition, stellar evolution modulates the whole baryon cycle depending on
stellar initial mass function (IMF), metal yields and lifetimes. The shape of
the IMF has implications for the resulting stellar population, its chemical
patterns and the galaxy luminosity function Leitherer et al. (1999); Bruzual &
Charlot (2003); Vazdekis et al. (2010); Tescari et al. (2014); Mancini et al.
(2015, 2016); Ma et al. (2015, 2017a); Valentini et al. (2019b).
Modifications in the supernova (SN) efficiency and/or in the mechanical
feedback can affect the whole energetics of the picture at all times and
induce variations in the amount of gas mass expelled from the collapsing
regions into the IGM Maio et al. (2011a); Campisi et al. (2011), with possible
consequences on atomic and molecular species (in particular $\Omega_{\rm
neutral}$ and $\Omega_{\rm H_{2}}$).
More or less powerful sources alter the ionization state of the IGM and
suppress gas cooling depending on the strength of radiative feedback and HI or
H2 shielding parameterisations Draine & Bertoldi (1996); O’Shea & Norman
(2008); Petkova & Springel (2011); Petkova & Maio (2012); Wise et al. (2012a);
Rahmati et al. (2013, 2015); Sternberg et al. (2014); Maio et al. (2016);
Wolcott-Green et al. (2017); Jones & Bate (2018).
How much all these processes affect the overall molecular chemistry in
realistic environments and at different redshifts is matter of debate Genel et
al. (2014); Furlong et al. (2015); Pallottini et al. (2017); Davé et al.
(2017, 2019); Donnari et al. (2019). In this respect, dedicated numerical
efforts to explore the role of the various physical phenomena that condition
the manifold regimes of cosmic environments must be performed to interpret
available data. Throughout this work we will shed light on the evolution of
cosmic atomic and molecular gas by interpreting new, recent, state-of-the-art
observations and by modelling HI and H2 species in the frame of high-
resolution N-body hydrodynamical simulations that include full time-dependent
non-equilibrium chemistry, star formation feedback and stellar evolution. A
further goal of our study is to show how neutral and molecular gas densities
and corresponding depletion times are modulated by different UV radiation
fields, HI or H2 gas self-shielding, stellar population parameters and the
local processes affecting cosmic chemistry in different epochs.
We adopt a flat $\rm\Lambda$CDM cosmological model, consistently with cosmic
microwave background data Komatsu et al. (2011), with present-day expansion
parameter normalised to 100 $\rm km/s/Mpc$ of $h=0.7$. Baryon, matter and
cosmological-constant density parameters are assumed to be
$\Omega_{0,\rm{b}}=0.045$, $\Omega_{0,\rm{m}}=0.274$ and
$\Omega_{0,\rm{\Lambda}}=0.726$, respectively. The spectral parameters are
$\sigma_{8}=0.8$ for the $z=0$ mass variance within 8 $\rm Mpc/{\it h}$ radius
and $n=0.968$ for the slope of the primordial power spectrum, consistent with
the latest Planck results Planck Collaboration (2020). The present-day
cosmological critical density is $\rho_{0,\rm{crit}}\simeq 277.4\,h^{2}\,\rm
M_{\odot}/kpc^{3}$.
The paper is organized as follows. Details on observational data, numerical
implementations and data analysis are given in Sect. 2; results are presented
in Sect. 3 and discussed in Sect. 4; conclusions are summarised in Sect. 5.

##  2 Method

In the following we present the most recent observations of cold atomic and
molecular gas, the theoretical models implemented (comprising several physical
processes newly included in the simulations), as well as a comprehensive list
of all our runs.

###  2.1 Observations

To quantify the amount of HI and H2 masses at various epochs, different
authors use either mass densities, $\rho_{\rm HI}$ and $\rho_{\rm H_{2}}$ (in
$\rm M_{\odot}/Mpc^{3}$), or the corresponding dimensionless $\Omega_{\rm HI}$
and $\Omega_{\rm H_{2}}$ mass density parameters. The evolution of the entire
amount of cold (below $10^{4}\,\rm K$) neutral gas density, $\rho_{\rm
neutral}$, can be obtained, as usually done in the literature, correcting
H-based measurements by helium content, i.e. scaling HI masses by a factor of
1.3 (e.g. Crighton et al. 2015). In the following we will discuss neutral and
molecular gas in terms of $\Omega_{\rm neutral}=\rho_{\rm
neutral}/\rho_{0,\rm{crit}}$ and $\Omega_{\rm H_{2}}=\rho_{\rm
H_{2}}/\rho_{0,\rm{crit}}$. In practice, these parameters correspond to the
mass fraction of the considered species times $\Omega_{0,\rm{b}}$.
We consider HI data at $z\lesssim 5$, as collected by Péroux & Howk (2020)
(PH20) in their Tab. 1, which rely on low-$z$ 21 cm emission ($z\lesssim 0.4$)
and high-$z$ quasar absorption data. At low redshift, 21 cm emission analysis
and spectral stacking are performed on data from 2dFGRS Delhaize et al.
(2013), WSRT Rhee et al. (2013); Hu et al. (2019), GMRT Rhee et al. (2016,
2018); Bera et al. (2019), Arecibo/ALFA Hoppmann et al. (2015) and ALFALFA
Jones et al. (2018) observational programs. The present-day estimated value of
the HI mass density parameter is $\simeq(4.02\pm 0.26)\times 10^{-4}$ Hu et
al. (2019). High-redshift measurements are obtained using data from high-
resolution spectroscopic surveys of tens of objects performed at VLT Zafar et
al. (2013), Calar Alto/CALIFA Sánchez et al. (2016), HST Rao et al. (2017),
and lower-resolution spectra of hundreds or thousands of objects recorded by
SDSS Noterdaeme et al. (2009, 2012); Crighton et al. (2015). The HI density
evolution is quite smooth, increasing by a factor of $\sim 3$ from the present
up to $z\simeq 5.3$ (see Fig. 3). Mass density parameters for cold neutral gas
vary within $\sim 0.5$ and $\sim 1.4\times 10^{-3}$ and can be fitted by eq.
(13) in Péroux & Howk (2020), i.e.: $\Omega_{\rm neutral}=(4.6\pm 0.2)\times
10^{-4}\,(1+z)^{0.57\pm 0.04}$.
Up-to-date H2 mass densities at different $z$ (for any given observational
programme) are found in Riechers et al. (2020a, b) (VLASPEC, COLDz), Decarli
et al. (2020) (ASPECS), Lenkić et al. (2020) (PHIBBS-2), Garratt et al. (2021)
(UKIDSS-UDS), and Hamanowicz et al. (2021) (ALMACAL-CO). Some of these are
also present in the compilation by Péroux & Howk (2020). Upper limits at
$z\simeq 0$-2 are provided by Klitsch et al. (2019) (ALMACAL-absorption,
ALMACAL-abs hereafter) and impose $\rho_{\rm H_{2}}<10^{8.39}\,\rm
M_{\odot}/Mpc^{3}$, corresponding to $\Omega_{\rm H_{2}}<1.81\times 10^{-3}$.
This value is much higher than early Herschel PEP lower limits of molecular-
mass density parameters $\gtrsim 10^{-5}$ inferred through indirect methods
for gas at $z\lesssim 3$ Berta et al. (2013). The currently established
present-day value of $\rho_{\rm H_{2}}\simeq 6.8\times 10^{6}\,\rm
M_{\odot}/Mpc^{3}$ is obtained by Fletcher et al. (2021) (xCOLDGASS) and
corresponds to $\Omega_{\rm H_{2}}\simeq 5\times 10^{-5}$ with a factor-of-2
relative error. ALMA and VLA observational programmes have determined an
increasing $\Omega_{\rm H_{2}}$ trend from such $z\simeq 0$ value to values
$\sim 10^{-4}$ at $z\simeq 2$-3 Decarli et al. (2020); Riechers et al.
(2020a). Infrared UKIDSS-UDS constraints Garratt et al. (2021) have confirmed
existing data, finding molecular-mass densities $\rho_{\rm H_{2}}\simeq
2\times 10^{7}\rm M_{\odot}/Mpc^{3}$ (i.e. $\Omega_{\rm H_{2}}\simeq 1.5\times
10^{-4}$) at $z\simeq 2$. H2 fractions up to 50-60% in galaxies at redshift
$z\sim 4$-6 have been lately reported by e.g. Dessauges-Zavadsky et al.
(2020), although H2 densities within 1$\sigma$ confidence levels decrease down
to $\Omega_{\rm H_{2}}\sim 10^{-5}$ towards $z\simeq 7$ Riechers et al.
(2020b). The derived $\Omega_{\rm H_{2}}$ values at low $z$ span more than two
orders of magnitude, from $10^{-6}$ to few $10^{-4}$, but they are more robust
at $z\gtrsim 2$. In practice, the overall ensemble of recent $\Omega_{\rm
H_{2}}$ observational findings show a variety of results comprised between
$10^{-3}$ and $10^{-6}$ in the $z\simeq 0$-7 range and the shape of the
resulting $\Omega_{\rm H_{2}}$ evolution is completely different from the
smooth $\Omega_{\rm neutral}$ one. A list with the detailed values derived by
different observational programs is presented in Appendix A (Tab. 2).
From $\rho_{\rm neutral}$ and $\rho_{\rm H_{2}}$, it is possible to derive the
corresponding depletion times, $t_{\rm depl,neutral}=\rho_{\rm
neutral}/\dot{\rho}_{\star}$ and $t_{\rm depl,H_{2}}=\rho_{\rm
H_{2}}/\dot{\rho}_{\star}$, where $\dot{\rho}_{\star}$ is the cosmic star
formation rate density.

###  2.2 Theoretical models

To perform a solid interpretation of observational data, we run and analyse a
number of numerical simulations, based on N-body and hydrodynamical
calculations with an extended version of the parallel code P-Gadget3 Springel
(2005) modified to address in detail cold-gas atomic and molecular content
over different cosmological epochs (ColdSIM). Besides gravity and smoothed-
particle hydrodynamics, the code contains a self-consistent time-dependent
treatment of atomic and molecular chemistry for e-, H, H+, H-, He, He+, He++,
H2, H${}^{+}_{2}$, D, D+, HD, HeH+, which follows, in non-equilibrium
conditions, a network of reactions for species ionization, recombination and
dissociation Abel et al. (1997); Yoshida et al. (2003); Maio et al. (2007).
Here the network has been extended to account for the main processes taking
place from the epoch of reionization down to lower redshift and these are
coupled to star formation, stellar evolution and heavy-element enrichment of
He, C, N, O, Ne, Mg, Si, S, Ca, Fe, etc. from SN II, AGB and SN Ia phases
Tornatore et al. (2007, 2010); Maio et al. (2010, 2011a). In the following
sections we describe in more detail the chemical processes included in the
original version of the code, as well as the additional processes included for
a better modelling of cold-gas evolution.

####  2.2.1 Basic chemistry evolution

Non-equilibrium chemistry is followed through a set of first-order
differential equations accounting for creation and destruction processes, as
well as corresponding heating or cooling, as described in e.g. Maio et al.
(2010) and Maio et al. (2016). Number density evolution for each species and
for each gas particle is computed by solving the differential equations
through a backward difference scheme with an integration time of 0.1 the
actual timestep (to grant numerical convergence; Anninos et al. 1997).
Stellar evolution is followed for a range of stellar masses and initial
metallicities. Stars with masses above $8\,\rm M_{\odot}$ explode as SNe II
and inject a typical kinetic energy of $10^{51}\,\rm erg$ in the surrounding
medium, while lower-mass stars evolve through AGB or SNe Ia phase with
consequent mass loss Tornatore et al. (2007). Star formation is modelled
stochastically (e.g. Katz 1992; Springel & Hernquist 2003) in particles with
density above a threshold of $10\,\rm cm^{-3}$ and temperature below
$10^{4}\,\rm K$ after the ”loitering” regime and the catastrophic run-away
cooling phase Maio et al. (2009). Wind feedback takes place with a reference
wind velocity of 350 km/s, while wind particles are decoupled from
hydrodynamics and treated in a collisionless fashion (i.e. they do not cool;
Maio et al. 2011a) for a delay time after the initial kick of 0.025-0.1 the
Hubble time, $t_{\rm H}$. Spreading of heavy elements during different stellar
phases is traced for He, C, N, O, Ne, Mg, Si, S, Ca, Fe, etc. Metal-dependent
yields are taken from Woosley & Weaver (1995), van den Hoek & Groenewegen
(1997), Woosley et al. (2002) and Thielemann et al. (2003), while stellar
lifetimes follow Padovani & Matteucci (1993) (see also Sect. 2.1 in Maio et
al. 2016, for more details). The reference initial mass function (IMF) has a
Salpeter (1955) shape over the stellar-mass range
${\rm[0.1,100]\leavevmode\nobreak\ M_{\odot}}$.
The chemical network is coupled to the insurgence of UV radiation according to
available $z$-dependent tabulated photoheating and photoionization rates.
These are added to the local gas heating caused by ongoing exothermic chemical
reactions and physical processes. Works in the recent literature have
suggested different epochs for the establishment of a UV background and
different temporal evolution of the adopted UV rates, which we test by means
of dedicated runs. The main physical differences among the UV models are
related to the treatment of the possible sources of cosmic re-ionization and
the availability of observational data to constrain model parameters at
various $z$. We adopt the standard literature prescriptions by Haardt & Madau
(1996) (HM), as well as the most up-to-date modelling including latest
high-$z$ data by Puchwein et al. (2019) (P19) and Faucher-Giguère (2020)
(FG20). While the HM case is a simple model based on early QSO data only, the
FG20 model include also the contributions from early galaxies with an onset of
HI and HeII photoheating at $z\simeq 8$ and 4, respectively (similarly to the
Haardt & Madau 2012, where HI and HeII photoheating kicks in at $z\simeq 9$
and 5) and seems to reproduce better the Ly$\alpha$ forest mean flux decrement
in the low-redshift IGM Christiansen et al. (2020). HM rates have a fairly
smooth evolution until $z\simeq 6$, although updated analyses by P19,
including contributions from both QSOs and galaxies, suggest rates that have a
tail extending up to $z\simeq 15$.
Although this set-up is already enough to follow the time-dependent HI and H2
evolution with a precision that is not commonly obtained in standard
implementations of cosmic structure formation models, to evaluate even more
accurately the properties of cold gas in different environments, additional
processes that might influence chemical abundances in complex ways should be
included, such as: HI and H2 self-shielding from UV radiation, H2 grain
catalysis, photoelectric as well as cosmic-ray heating. In the following we
will discuss separately all these newly implemented processes.

####  2.2.2 Gas HI and H2 self-shielding

While UV radiation can easily penetrate thin media, gas self-shielding effects
should be taken into account when estimating HI and H2 fractions in dense gas.
For H2 self-shielding we adopt the formulae by Draine & Bertoldi (1996), who
fit their results either by means of a single power law (their eq. 36) or by a
slightly more accurate expression (their eq. 37). We also use the more recent
results of the work by Wolcott-Green et al. (2017), which, similarly to
Sternberg et al. (2014), is a reformulation of the ones by Draine & Bertoldi
(1996). Most runs discussed here include H2 shielding according to eq. 36 of
Draine & Bertoldi (1996) if not otherwise stated. In the case of Wolcott-Green
et al. (2017), the optimal scale to estimate the local H2 column density,
$N_{\rm H_{2}}$, and related shielding have been investigated in depth by Luo
et al. (2020). They showed that $N_{\rm H_{2}}$ can be quantified as $N_{\rm
H_{2}}=n_{\rm H_{2}}\lambda_{\rm J}/4$, where $n_{\rm H_{2}}$ is the the
three-dimensional H2 number density and $\lambda_{\rm J}$ the Jeans length.
This has been proven to be an excellent approximation. In Appendix B we show
that different H2 shielding implementations have a negligible effect on our
results.
For HI self-shielding we rely on the redshift and density dependent
corrections proposed by Rahmati et al. (2013). We note that, while the authors
provide corrections at $z\leq 5$, both observational data and theoretical
investigations suggest the presence of cold shielded gas also at higher
redshifts. Despite the numerical uncertainties discussed in the original
paper, the redshift evolution of self-shielding corrections as a function of
density presents only little (non-monotonic) variations and deviations from
the $z$-averaged mean/median values for any given density are modest.
Therefore, to capture the essential behaviour of HI, we fit the median values
of those shielding curves. As the redshift dependence is weak and in the
absence of a better prescription, we extend such relation to $z>5$. The
resulting density-dependent fitting expression is given in Appendix B (eq. 1).
This nicely shows that gas with H number densities $n_{\rm H}\lesssim
10^{-3}\,\rm cm^{-3}$ is essentially un-shielded (i.e. the correction factor
is close to unity), while at larger number densities the effective
photoionization rate decreases by several orders of magnitude.
Thus, while we find that the implementation of gas self-shielding in the
simulations affects the resulting $\Omega_{\rm neutral}$ and $\Omega_{\rm
H_{2}}$ (this will be discussed in detail in Sect. 3), variations in its exact
formulation have a minor impact.

####  2.2.3 H2 dust grain catalysis

We have coupled H2 gas-grain catalysis to the non-equilibrium chemical network
and followed it for different gas temperatures and metallicities. Condensation
of metals into dust grains can enhance H2 formation via $\rm H$ abstraction,
although the quantum modelling of H2 formation on dust grains bears
uncertainties of a factor of 2 Bron et al. (2014). Chemical rates for H2 gas-
grain catalysis are expected to be of the order of $\sim 10^{-17}\,\rm
cm^{3}/s$ Duley (1996) and generally scale with the sticking coefficient (i.e.
the fraction of particles that adsorb to the grain surface), $S_{\rm H}$, and
the square root of gas temperature, $\sqrt{T}$. Cazaux & Tielens (2004)
suggested a rate of $3\times 10^{-17}\,\rm cm^{3}/s$ and a sticking
coefficient following a second-order polynomial depending on the adopted grain
temperature, $T_{\rm gr}$ Hollenbach & McKee (1979). Based on Tielens &
Hollenbach (1985), Omukai et al. (2005) propose an H2 rate of $6\times
10^{-17}\,\rm cm^{3}/s$ for gas at 300 K (corresponding to about $3.5\times
10^{-17}\,\rm cm^{3}/s$ at 100 K) and the same sticking coefficient as
Hollenbach & McKee (1979). Similar approaches have been used later on e.g. by
Tomassetti et al. (2015) and more recently by Sillero et al. (2021). Updated
quantum-mechanical calculations by Thi et al. (2020) find a value of
$3.74\times 10^{-17}\,\rm cm^{3}/s$ for a gas at 80 K assuming a constant
$S_{\rm H}$ of unity. Broadly speaking, all these results are in line with
observational constraints of 3-$4.5\times 10^{-17}\,\rm cm^{3}/s$ Jura (1975a,
b); Gry et al. (2002), although the value by Thi et al. (2020) does not
incorporate a temperature-dependent $S_{\rm H}$. Gas-grain energy losses scale
proportionally to $\sqrt{T}\,(T-T_{\rm gr})$, hence contributing to gas
heating or cooling depending on both $T$ and $T_{\rm gr}$. On the base of
Milky Way calibrations, the normalization is expected to be $\sim
10^{-33}\,\rm erg/s/K^{3/2}cm^{3}$, with a roughly 1-dex error at low $Z$.
Such large deviations introduce a remarkable degree of freedom in the
parameter study, suggesting for example that variations of the dust-to-gas
ratio within a factor of 10 would be still consistent with the constraints.
For this reason, the exact shape and $Z$ scaling of dust-induced H2 formation
and cooling/heating rates have been long debated in the literature Hollenbach
& McKee (1989); Omukai et al. (2005); Draine & Li (2007); Goldsmith et al.
(2011); Krumholz & Gnedin (2011); Glover & Clark (2012).
Here, we adopt an H2 formation rate of $3.5\times 10^{-17}\,\rm cm^{3}/s$ at
100 K, scaled by $\sqrt{T}$ and $S_{\rm H}$ Hollenbach & McKee (1979); Cazaux
& Tielens (2004) (see Appendix C) to update the H2 abundances employed in the
cooling/heating calculations Omukai et al. (2005). Since observational
determinations are based mainly on data from the Milky Way, we scale the rates
linearly with $Z$ at metallicities different from solar. As a fiducial value
we assume $T_{\rm gr}=75\,\rm K$ at all $z$. Nevertheless, we additionally
check alternative hypotheses (see Appendix C), i.e. a lower and a higher value
of 40 K and 120 K, respectively, both consistent with observational
constraints and theoretical estimates of grain properties555 Local
determinations of $T_{\rm gr}$ are related to gas surrounding nearby early-
type stars (60-120 K) and photodissociation regions (65 K for $\theta^{1}$ Ori
C), instead theoretical estimates refer to H astration modeling on a perfect
surface (75 K) or completely imperfect surface (44 K). Updated observational
best-fit values at different redshifts are $\sim 41$ K at $z\simeq 4.5$ and
38-43 K at $z\simeq 5.5$ Khusanova et al. (2020); Faisst et al. (2020). . To
study the effects of an evolving $T_{\rm gr}$, we also consider either grain
temperatures varying as the CMB temperature, i.e. $T_{\rm gr}(z)=T_{\rm
0,CMB}(1+z)$ with $T_{\rm 0,CMB}\simeq 2.7\,\rm K$, or varying in accordance
to the energy balance between CMB radiation and dust grain power-law emission
Draine & Lee (1984), i.e. $T_{\rm gr}^{4+\beta}(z)=T_{\rm
0,gr}^{4+\beta}+T_{\rm 0,CMB}^{4+\beta}[(1+z)^{4+\beta}-1],$ with $\beta\simeq
2$ emissivity index666 Common $\beta$ values are between roughly 1 and
slightly greater than 2, but the resulting $T_{\rm gr}(z)$ is not very
sensitive to particular choices. chosen in accordance to the latest ALMA
determinations da Cunha et al. (2021) and $T_{\rm 0,gr}=18\,\rm K$ present
grain temperature da Cunha et al. (2013). We explore degeneracies and
metallicity effects at different cosmological epochs implementing also the
possibility to adopt fixed values for $Z$ (e.g. 0.01-$1\,\rm Z_{\odot}$).
Different $T_{\rm gr}$, either constant or $z$-dependent, cause local
fluctuations of comparable amplitude on H2 chemistry, but have little impact
on the global $\Omega_{\rm H_{2}}$ behaviour. Metallicity is relevant for
values close to $Z_{\odot}$, which boost both local H2 fractions (by one dex
since $z\simeq 8$) and density parameters (by more than two dex at $z\gtrsim
10$ and by a factor of a few at $z\simeq 4$).
In practice, H2 grain catalysis has a relevant local impact in metal enriched
sites, while, generally speaking, its global role decreases at higher $z$ as a
result of the decreasing metallicity. Finally, it does not affect the global
HI chemical history and $\Omega_{\rm neutral}$ values do not vary
significantly from the reference case.

####  2.2.4 Photoelectric heating

Gas-grain interactions are usually accompanied by photoelectric heating. This
is one of the main heating processes of the diffuse interstellar medium (ISM)
and is caused by the absorption of UV photons by dust grains which then re-
emit electrons and heat up the surrounding gas. According to standard
expressions commonly used in ISM investigations Weingartner & Draine (2001),
the photoelectric heating rate density can be written as a function of gas
temperature, Habing parameter $G_{0}$ quantifying the local UV field (Habing
1968) and electron fraction Draine & Sutin (1987); Bakes & Tielens (1994).
With the exception of Milky Way calibrations, the details of such process in
different cosmological contexts are still unknown. We thus either consider the
fiducial constant case of $G_{0}=1.7$, or, because of the close relation
between star formation and UV photon production, we scale it by the local SFR
Narayanan & Krumholz (2017), i.e. $G_{0}\rightarrow G_{0}\times\rm
SFR/SFR_{\odot}$, where $\rm SFR_{\odot}\sim 2\,M_{\odot}/yr$ is the value of
the Galactic SFR Chomiuk & Povich (2011); Licquia & Newman (2015).777 A
precise quantification of the local SFR depends on the IMF and different
assumptions can lead to errors up to $40\%$. This is smaller than typical
uncertainties of some chemistry rates (known within up to a factor of 2),
hence induced variations should not be dramatic.  We have verified that the HI
and H2 densities obtained with both prescriptions are very similar, with
differences visible only on local scales (see details in Appendix D). We note
that, while heating the surrounding gas, the photoelectric effect provides
additional free electrons, useful for the H- channel. These charges react with
H via radiative attachment and increase H2 formation via associative
detachment.
The H2 abundances reached in the presence of photoelectric heating are higher
by a factor of a few mostly at early times, while the overall effect on
$\Omega_{\rm neutral}$ is modest.

####  2.2.5 Cosmic-ray heating

Another source that likely plays a role for H2 evolution is the cosmic-ray
heating. Cosmic rays are mostly constituted by high-energy protons usually
formed in explosive SN events that can impact on molecular-mass build-up
during and after the epoch of reionization. It is known that heating by cosmic
rays is related to ionization of atoms and consequent production of free
electrons. Each cosmic-ray proton ionizes the medium generating slow electrons
and these lead to second or later-generation ionization in 2/3 of the cases
(causing an increase in the yield of liberated energy by 5/3) Spitzer & Scott
(1969); Opal et al. (1971). About $50\%$ of the energy released goes into gas
heating. Because of their relatively large abundances, the species mainly
involved in the ionization process and accounting for most of the gas heating
are H, He and H2 Glassgold et al. (2012), while others are almost
irrelevant888 The probability that cosmic-ray protons interact with e.g. dust
or heavy molecules is expected to be orders of magnitude smaller than
interactions with H2 molecules. As a reference, the dust cross section per H
nucleus in the diffuse ISM is $\sim 10^{-21}\,\rm cm^{-2}$, while in the case
of H2 it is $\sim 10^{-17}\,\rm cm^{2}$. Even if interactions with dust should
occur, only $1\%$ of them would produce photoelectrons Dwek & Smith (1996). .
Gas heating due to cosmic rays is usually parameterized as a function of the
heating yield and ionization rate. A wide range of values exists in the
literature for the cosmic-ray heating yield, which is estimated to be $\sim
7$-20 eV (Stahler & Palla 2004; Goldsmith & Langer 1978, respectively) when
considering electron energies higher than 1 keV, but reach about 44 eV if
including also lower electron energies Glassgold et al. (2012). As the
normalisation of the transferred energy is weakly dependent on the local
environment, we always assume 20 eV as typical value. The ionization rate is
expected to be $\zeta\sim 10^{-17}\,\rm s^{-1}$ per H atom, although IR
measurements of $\zeta\gtrsim 10^{-16}{\rm s}^{-1}$ in dense molecular sites
Indriolo et al. (2007); Indriolo & McCall (2012); Hollenbach et al. (2012)
demonstrated that it can vary with density and reach larger values. New
analysis by Thi et al. (2020) has suggested a typical ionization rate of
$1.7\times 10^{-17}\,\rm s^{-1}$ per H2 molecule. Recent theoretical results
by Padovani et al. (2018) have shown that $\zeta$ could be described by a
decreasing function of the density, peaking at a value of $\sim 10^{-16}\,\rm
s^{-1}$ at gas column densities of $10^{19}\,\rm cm^{-2}$. This estimate is
more than a factor of ten larger than the one by Thi et al. (2020), but still
consistent with SN energies of $10^{51}\,\rm erg$ and with uncertainties
related to the different derivation technique.
In our implementation we couple cosmic-ray heating with non-equilibrium
chemistry to handle both a constant $\zeta=10^{-17}\,\rm s^{-1}$ and a
density-dependent $\zeta$, like the one expected by Padovani et al. (2018). In
this latter case, we fit tabulated results of their model $\mathcal{L}$ in a
friendly form (given in Appendix E). To constrain extreme behaviours, cosmic-
ray heating is applied either to all gas particles in the simulated volume or
to star forming particles only. Because of the lack of data for ionization
rates in regimes different from the local one, $\zeta$ is generally scaled by
the local SFR, as it is linked to SN shocks in star forming regions.
While we refer the reader to Appendix E for more details on the impact of
different prescriptions, here we just mention that the inclusion of our
reference cosmic-ray heating (i.e. a density-dependent $\zeta$ applied to star
forming particles) affects only mildly the density parameters at all
redshifts.

####  2.2.6 Star formation model

Finally, commonly used star formation implementations that convert gas into
stars stochastically estimate the SFR relying on the ratio between the gas
density above a given threshold, $\rho_{\rm th}$, and a typical star formation
time that scales with gas density. In particular, the star formation timescale
can be written as $t_{\rm sf}=t_{\rm sf,0}(\rho/\rho_{\rm th})^{1/2}$, where
the maximum star formation time $t_{\rm sf,0}$ is an empirical parameter that
varies substantially in different observations. At each integration time,
collisionless star particles are spawned stochastically with a probability
given by $p=(m/m_{\rm star})\\{1-\exp[-(1-\beta)x\Delta t/t_{\rm sf}]\\}$,
with $m$ gas particle, $m_{\rm star}$ stellar particle, $\beta$ fraction of
short-lived stars that die as SNae, $x$ cold fraction and $\Delta t$ numerical
timestep Springel & Hernquist (2003). The resulting star formation rate
associated to each gas parcel is then $\dot{M}_{\rm sf}=(1-\beta)xm/t_{\rm
sf}$. Here we adopt a reference value of $t_{\rm sf,0}=2\leavevmode\nobreak\
\rm Gyr$ (in line with Leroy et al. 2013; Tacconi et al. 2013, 2020). The
value of $m_{\rm star}$ is estimated by assuming that each gas particle forms
4 stellar particles (4 ‘generations’), $\beta$ is computed according to the
adopted IMF and $x$ follows Springel & Hernquist (2003), resulting usually
close to unity in high-density regions. In practice, all this means that at
the threshold value the actual star formation time is $t_{\rm sf}=t_{\rm
sf,0}$ and the star formation probability $p$ is small, being $\Delta t/t_{\rm
sf}\ll 1$ typically. At larger densities, $t_{\rm sf}$ becomes shorter, while
the star formation probability $p$ increases. Since we are able to follow
detailed non-equilibrium H2 evolution, we additionally implement an H2-based
approach in which the SFR is driven by the non-equilibrium H2 density in each
gas parcel above the star formation threshold of $10\leavevmode\nobreak\ \rm
cm^{-3}$. In general, the H2 star formation efficiency can be thought as
implicitly coded in the star formation timescale, thus efficiency values equal
to (lower than) unity correspond to timescale values equal to (higher than)
the reference timescale. We run corresponding simulations and discuss the
implications for HI and H2 evolution in the next sections.

###  2.3 Numerical simulations

All the processes mentioned above and newly implemented in our theoretical
framework have an impact on the amount of free electrons or protons available
for H2 formation, in addition to increasing gas heating or cooling and
associated species destruction or formation. What their net result is in
different environments and at different epochs needs to be investigated.
To this goal, we run a number of numerical simulations in a cubic volume of
length 10 co-moving Mpc/h (ColdSIM), as listed in Tab. 1. Gas and dark-matter
fields are each sampled with 5123 particles at an initial redshift of 100
(different initial redshifts have a modest effect on the results; e.g. Maio et
al. 2011b). This set-up corresponds to a resolution of $7.89\times 10^{4}\,\rm
M_{\odot}$ and $4.72\times 10^{5}\,\rm M_{\odot}$, for gas and dark-matter
particles, respectively. In our reference simulation (HM-HISSmed in Tab. 1)
star formation happens in gas above 10 cm-3 with a Salpeter (1955) IMF. It has
a SN efficiency of 0.1, wind velocity of 350 km/s, wind delay time of
$0.025\,t_{\rm H}$, HM UV background, density-dependent HI self-shielding
according to eq. 1 and H2 self-shielding by Draine & Bertoldi (1996) (their
eq. 36). Additional simulations have been run to explore and quantify the role
of different assumptions (for UV background, star formation and feedback
effects related to IMF, SNe or winds), as well as gas physics and chemistry
processes (HI or H2 shielding, H2 formation channels, H2 dust grain processes,
cosmic-ray heating), as detailed in Tab. 1. Implications from the different
HM, P19 and FG20 backgrounds (highlighted with different colors in the table
and in the next figures) have been addressed with dedicated HI-shielded or un-
shielded simulations. We note that most HI-shielded runs include HI self-
shielding according to eq. 1 (and are denoted by the HISSmed suffix after the
UV background name, i.e. HM-HISSmed, P19-HISSmed and FG20-HISSmed), while HI
self-shielding according to Rahmati et al. (2013)’s tables is tested only in
selected runs (denoted via the HISS suffix after the UV background name, i.e.
HM-HISS, P19-HISS and FG20-HISS). The gas physics and chemistry modelling
described above has been investigated by always using the HM background. For
example, this choice has been done when adopting (non-equilibrium) H2-based
star formation (HM-HISSmed-H2). Although our fiducial star formation threshold
is $10\,\rm cm^{-3}$, comparisons have been made to investigate the impact on
the H2 mass density parameter of a standard low-density threshold $\sim
0.1\,\rm cm^{-3}$ with an HM background (HM-std and HM-HISSmed-std), finding
that in this latter case H2 fractions are not resolved properly. Similarly, we
have adopted the HM scenario to quantify the contribution of each H2 formation
channel by including in the reference run the three-body processes together
with the H- and H${}_{2}^{+}$ channels (HM-HISSmed-3b). H2 self-shielding by
Draine & Bertoldi (1996) (eq. 37) or Wolcott-Green et al. (2017) are tested
via HM-HISSmed-3b-DB2 and HM-HISSmed-3b-WG runs. To isolate the implications
of individual channels, two further tests have been implemented, the first
including the H${}_{2}^{+}$ channel only (HM-HISSmed-H${}_{2}^{+}$) and the
second including none of the H2 channels (HM-HISSmed-none). A similar analysis
has been done for H2 grain catalysis at different grain temperatures and
metallicities (see Tab. 1), with or without HI shielding (HM-HISSmed-
cat-75-Zevol, HM-cat-75-Zevol), or photoelectric heating (HM-HISSmed-
cat-75-pe). Cosmic-ray heating with constant or density-dependent Padovani et
al. (2018) ionization rate in shielded gas has been explored by means of HM-
HISSmed-cr and HM-HISSmed-crP18 simulations, respectively. To have a more
realistic analysis, we additionally consider the possibility that cosmic-ray
heating is injected only in star forming parcels (HM-HISSmed-crP18sfr).
Effects linked to variations of the SN efficiency and IMF have been tested by
simply employing un-shielded UV backgrounds to maximize the impact of stellar
feedback. Other features related to powerful winds (at 700 km/s) or winds re-
coupled after a longer time (wind delay time of $0.1\,t_{\rm H}$, 4 times the
reference $0.025\,t_{\rm H}$) are explored, too. For each output time, cosmic
structures are identified by means of a friends-of-friends algorithm with a
linking length of $20$ per cent the mean inter-particle separation.
Substructures are identified by the Subfind algorithm (Dolag et al. 2009, and
references therein) and are post-processed to trace: masses, positions, radii,
velocities, star formation rates, temperatures, chemical abundances (of e-, H,
H+, H-, He, He+, He++, H2, H${}^{+}_{2}$, D, D+, HD, HeH+, C, N, O, Ne, Mg,
Si, S, Ca, Fe, etc.), substructures and all the relevant physical properties
of each object. Our choices for initial conditions, box size and resolution
are determined by the necessary trade-off between the required accuracy of the
physical descriptions implemented in the code and the numerical feasibility of
the runs.

Table 1:  Inputs to simulation models.

Properties | Simulation run | UV | HI shielding |  H2 shielding | IMF | $\varepsilon_{\rm SN}$ | Extra modules
---|---|---|---|---|---|---|---
(a) | (b) | (c) | (d) | (e) | (f) | (g) | (h)
| HM | HM | – | DB96 | S | 0.1 | –
| HM-HISS | HM | R13 | DB96 | S | 0.1 | –
|  HM-HISSmed (ref) | HM | 1 | DB96 | S | 0.1 | –
| P19 | P19 | – | DB96 | S | 0.1 | –
UV rate | P19-HISS | P19 | R13 | DB96 | S | 0.1 | –
modeling | P19-HISSmed | P19 | 1 | DB96 | S | 0.1 | –
| FG20 | FG20 | – | DB96 | S | 0.1 | –
| FG20-HISS | FG20 | R13 | DB96 | S | 0.1 | –
| FG20-HISSmed | FG20 | 1 | DB96 | S | 0.1 | –
| HM-std | HM | – | DB96 | S | 0.1 |  threshold $0.1\,\rm cm^{-3}$
star formation | HM-HISSmed-std | HM | 1 | DB96 | S | 0.1 |  threshold $0.1\,\rm cm^{-3}$
| HM-HISSmed-H2 | HM | 1 | DB96 | S | 0.1 |  non-equilibrium H2-based
| HM-HISSmed-3b | HM | 1 | DB96 | S | 0.1 | 3-body processes
| HM-HISSmed-3b-DB2 | HM | 1 | DB96 (eq. 37) | S | 0.1 | 3-body processes
| HM-HISSmed-3b-WG | HM | 1 | WG17 | S | 0.1 | 3-body processes
|  HM-HISSmed-H${}_{2}^{+}$ | HM | 1 | DB96 | S | 0.1 |  H${}_{2}^{+}$ channel only
| HM-HISSmed-none | HM | 1 | DB96 | S | 0.1 |  no H2 channels
| HM-HISSmed-cat-75-Zevol | HM | 1 | DB96 | S | 0.1 | catalysis at 75 K
| HM-HISSmed-cat-75-pe | HM | 1 | DB96 | S | 0.1 | catalysis at 75 K and
|  |  |  |  |  |  | photoelectric heating
gas physics | HM-HISS-cat-40-Zevol | HM | R13 | DB96 | S | 0.1 | catalysis at 40 K
& chemistry | HM-HISS-cat-75-Zevol | HM | R13 | DB96 | S | 0.1 | catalysis at 75 K
| HM-HISS-cat-120-Zevol | HM | R13 | DB96 | S | 0.1 | catalysis at 120 K
| HM-HISS-cat-Tcmb-Zevol | HM | R13 | DB96 | S | 0.1 |  catalysis at $T_{\rm CMB}(z)$
| HM-HISS-cat-Tbeta-Zevol | HM | R13 | DB96 | S | 0.1 |  catalysis at $T_{\beta}(z)$
| HM-HISS-cat-75-Zsun | HM | R13 | DB96 | S | 0.1 |  catalysis at 75 K and $Z=\rm Z_{\odot}$
| HM-HISS-cat-75-0.01Zsun | HM | R13 | DB96 | S | 0.1 |  catalysis at 75 K and $Z=10^{-2}\rm Z_{\odot}$
| HM-HISS-cat-75-Zt | HM | R13 | DB96 | S | 0.1 |  catalysis at 75 K and $Z=Z(t)$
| HM-cat-75-Zevol | HM | – | DB96 | S | 0.1 | catalysis at 75 K
| HM-HISSmed-cr | HM | 1 | DB96 | S | 0.1 |  constant $\zeta$ , all gas
| HM-HISSmed-cr-P18 | HM | 1 | DB96 | S | 0.1 |  $\zeta$ by P18,∗ all gas
| HM-HISSmed-cr-P18sfr | HM | 1 | DB96 | S | 0.1 |  $\zeta$ by P18,∗ SFR-scaled
| HM-1 | HM | – | DB96 | S | 1.0 | –
| HM-0.5 | HM | – | DB96 | S | 0.5 | –
| HM-0.01 | HM | – | DB96 | S | 0.01 | –
| HM-w700 | HM | – | DB96 | S | 0.1 | wind velocity of 700 km/s
feedback | HM-wf | HM | – | DB96 | S | 0.1 |  winds decoupled for $0.1\,t_{\rm H}$
& IMF | HM-Chabrier | HM | – | DB96 | Ch | 0.1 | –
| HM-1-Chabrier | HM | – | DB96 | Ch | 1.0 | –
| P19-1 | P19 | – | DB96 | S | 1.0 | –
| P19-Chabrier | P19 | – | DB96 | Ch | 0.1 | –
| P19-1-Chabrier | P19 | – | DB96 | Ch | 1.0 | –

Blocks in the first column (a) collect, from top to bottom, runs with
different implementation properties. Following columns refer to: simulation
name (b), cosmological UV rates (c), HI shielding (d), H2 shielding (e), IMF
(f), SN energy fraction (g), and description of extra modules adopted besides
basics ones (h). UV backgrounds are by Haardt & Madau (1996, HM), Puchwein et
al. (2019, P19) or Faucher-Giguère (2020, FG20); HI shielding by Rahmati et
al. (2013, R13) or eq. 1 in Appendix B; H2 shielding by Draine & Bertoldi
(1996, DB96, their eq. 36 if not otherwise specified) or Wolcott-Green et al.
(2017, WG17); IMF by Salpeter (1955, S) or Chabrier (2003, Ch).

$\ast$ Simple fit in Appendix E to the $\mathcal{L}$ model by Padovani et al.
(2018).

##  3 Results

In the following we present the main results obtained by the analysis of the
simulations and their behaviour with respect to the observations described
before.

###  3.1 Gas thermal state

P19 FG20 HM

Figure 1: Phase diagrams at redshift $z=15.0$ (upper row) and $z=8.1$ (lower
row) for the simulations in Tab. 1 run with UV rates by Puchwein et al. (2019)
(P19, left), Faucher-Giguère (2020) (FG20, center) and Haardt & Madau (1996)
(HM, right). The simulation data for temperature, $T$, and gas overdensity
with respect to the gas mean, $\delta$, have been gridded on a base-10
logarithmic scale. The color coding refers to the total probability density
function (PDF).

0.01, 350 km/s 1, 350 km/s 0.1, 700 km/s

Figure 2: Phase diagrams at redshift $z=6.1$ in scenarios with HM background
and, from left to right a SN efficiency (wind velocity) of: 0.01 (350 km/s), 1
(350 km/s), 0.1 (700 km/s). These are the runs HM-0.01, HM-1 and HM-w700,
respectively.

We check the impact of the different UV backgrounds by looking at the effects
on the thermal state of cosmic gas. In Fig. 1 we show the density-temperature
diagram (phase space) for indicative runs, in terms of gas overdensity,
$\delta=\rho/\overline{\rho}$, where $\rho$ is the mass density of each gas
parcel and $\overline{\rho}$ the mean gas density. More specifically, phase
space diagrams for three representative runs with P19, FG20 and HM UV rates,
are displayed and the color code reflects the corresponding total probability
distribution function. In these models UV photoionization is switched on at
$z=$ 15.1 (P19), 8.3 (FG20) and 6 (HM). At early times ($z=15$), trends in
simulations run with FG20 and HM are similar, with most of the gas being cold,
below $10^{4}\,\rm K$. This results from the natural evolution of gas
chemistry in absence of UV radiation. Low-density gas is shock-heated and the
molecular fraction stays close to the initial-condition value (a few
$10^{-6}$). At higher densities, instead, non-equilibrium chemistry induces a
rapid formation of H2 molecules and consequent runaway cooling. In the P19
case, IGM/CGM-like material (around and below the mean density) has just been
warmed up by the newly injected UV photoheating and experiences an inversion
of its equation of state. High-density regions are little affected, though, as
a result of gas clumpiness in such regimes. Obviously, how photoheating and
photodissociation proceed depend on the chosen scenario, as highlighted by the
lower panels at $z=8.1$, when the gas thermal evolution has reached different
stages: the traditional HM background has not heated cold, thin, cosmic gas,
yet; in the FG20 and P19 panels, instead, photoheating is clearly effective
since $z=8.3$ for FG20 and $z=15.1$ for P19. We note that the runaway
molecular cooling branch is only mildly affected by UV radiation, as a
consequence of the larger densities that shield H2 molecules. At $z\lesssim 6$
all three models present phase diagrams similar to the latter case. We note
that, irrespectively from the chosen UV background rates, high-density regions
behave in a similar way, dictated by local thermal processes rather than
background radiation. Indeed, the star-forming branch exhibits the equation of
state of SN-heated gas based on the multiphase ISM model of Springel &
Hernquist (2003). This behaviour is caused by the SN energy released in the
surrounding medium and is responsible for the increase in gas temperatures
visible on the right of each the panel.

Collapsing star forming regions are more sensitive to feedback effects and the
implications of different assumptions for feedback parameters might be
relevant for the thermodynamical evolution of cosmic gas, independently from
the adopted UV prescriptions. As indicative examples, in Fig. 2 we show phase
diagrams at redshift $z=6.1$ for the runs with SN efficiencies of 0.01-1 and
wind velocity of 350 or 700 km/s in the HM case. A low SN efficiency of 0.01
(left panel, HM-0.01) is incapable of removing dense material from the star
forming regions (at $\rm Log_{10}\delta\gtrsim 1$), where gas retains
temperatures of $10^{4}$-$10^{5}\,\rm K$ without cooling further (red/blue
branch). A larger efficiency of 1 (central panel, HM-1) removes more gas, as
clear from the probability distribution of the star forming branch (red region
at $\rm Log_{10}\delta\sim 1$-2.5). Such material circulates through the
IGM/CGM again and becomes available for further cooling and fragmentation. The
0.1 case with the same winds (not plotted) features an intermediate behaviour.
The results displayed both in the left and in the central panel have a wind
velocity of 350 km/s. Stronger winds with velocities of e.g. 700 km/s (right
panel with SN efficiency of 0.1, HM-w700) impact gas removal from star forming
sites. In the case shown, dense material is ejected to outer regions more
efficiently and populates abundantly the thin IGM at $\rm
Log_{10}\delta\lesssim-1$ and $T\gtrsim 10^{5}\,\rm K$. Besides the role of
wind velocity, other wind parameters, such as the wind delay time before re-
coupling to hydrodynamics (e.g. $0.1\,t_{\rm H}$ adopted in HM-wf, instead of
the reference $0.025\,t_{\rm H}$ adopted in Fig. 2) are expected to have less
vital effects. The picture at different redshifts is qualitatively similar, as
it is when assuming a Chabrier IMF (HM-Chabrier and HM-1-Chabrier), instead of
a standard Salpeter IMF. Alternative UV rates for these models (e.g. P19,
P19-1, P19-Chabrier, P19-1-Chabrier) show analogous behaviours, except for the
mentioned (temporal) differences. These preliminary considerations highlight
the main impact of different assumptions on the UV background and feedback
effects, specifically on the thin cosmic medium. The implications for cold
dense media (privileged sites hosting neutral and molecular gas) are analyzed
in the next Sect. 3.2.

###  3.2 Density parameters

In the following we focus on the evolution of mass density parameters for
neutral gas, $\Omega_{\rm neutral}$, and molecular H2 gas, $\Omega_{\rm
H_{2}}$, as extracted from simulation outputs at different times, and compare
them to state-of-the-art observational data (Sect. 2.1) in Fig. 3.

Figure 3: Upper panel: redshift evolution of $\Omega_{\rm neutral}$ compared
to inferred values from HI observations with $1\sigma$ error bars (grey
points) and corresponding fit (from darker to lighter shades) at 1, 2 and
3$\sigma$ confidence levels (see Péroux & Howk 2020). The simulations
considered are: HM-HISSmed (thick black solid lines), P19-HISSmed (thick blue
dashed lines) and FG20-HISSmed (thick orange dot-dot-dot-dashed lines). They
all include HI self-shielding. As a comparison, corresponding un-shielded runs
with the same UV backgrounds (HM, P19 and FG20) are displayed by thin lines of
the same color and style. Lower panel: corresponding evolution of $\Omega_{\rm
H_{2}}$ compared to the observational constraints listed in Tab. 2 (PH20 1, 2
and 3$\sigma$ confidence levels, shades; ALMACAL-CO, light grey; VLASPEC,
magenta; PHIBBS-2, purple; COLDz, yellow; ASPECS, green; UKIDSS-UDS, cyan;
ALMACAL-abs and Hershel PEP limits, dark grey).

Theoretical trends refer to runs with the reference implementation including
time-dependent non-equilibrium abundance computations, a high-density
threshold of $10\,\rm cm^{-3}$ for stochastic star formation, stellar
evolution, metal spreading, and H2 self-shielding. We distinguish between
simulations accounting for HI self-shielding, (HM-HISSmed, P19-HISSmed and
FG20-HISSmed denoted by thick lines) and analogous un-shielded results (HM,
P19 and FG20 denoted by thin lines).

In the top panel, we focus on $\Omega_{\rm neutral}$. We can see that
theoretical results increase smoothly with redshift at $z\simeq 2$-5, in broad
agreement with observational data. UV radiation impacts our results after the
on-set of the UV background and leaves the signal from earlier epochs
unaffected (as visible for all the UV models e.g. in Appendix G). In both HI
shielded and un-shielded HM scenarios, the cold mass surviving at $z\lesssim
6$, after photoionization is turned on, is only 5% of the original one,
corresponding to $\Omega_{\rm neutral}\sim 10^{-3}$. As HI shielding preserves
HI mass, the resulting $\Omega_{\rm neutral}$ at $z\simeq 2$-3 is larger by a
factor of up to 2 with respect to the un-shielded case. Although both curves
are broadly consistent with the available data at $z\lesssim 5$, the un-
shielded run seems to predict a steeper evolution of cold gas. The shielded HM
scenario is slightly above the 1$\sigma$ error bars of observational data
points. In the P19 and FG20 cases the UV model is turned on earlier, resulting
in a larger effect of the associated feedback and flatter curves (compared to
those obtained with an HM background). In these cases, the predictions are
slightly below 1$\sigma$ observational determinations for the un-shielded
runs, while the corresponding shielded cases are roughly consistent with
observational lower limits. The addition of further physical processes, as we
will see in the following, does not strongly affect $\Omega_{\rm neutral}$,
whose values remain robustly close to the ones displayed in the figure and
within the error bars of available data. We note that the displayed fitting
expression by Péroux & Howk (2020), i.e. $\Omega_{\rm
neutral}\propto(1+z)^{0.57}$, is similar to other independent formulations
finding $\Omega_{\rm neutral}\propto(1+z)^{0.60}$ Crighton et al. (2015) or
$\Omega_{\rm neutral}\propto(1+z)^{0.64}$ Rao et al. (2017). Thus,
uncertainties on the numerical fit do not affect our conclusions.

In the bottom panel of Fig. 3 we explore the outcome for $\Omega_{\rm H_{2}}$,
compared to the observational constraints (listed in Tab. 2) by Hamanowicz et
al. (2021) (ALMACAL-CO), Riechers et al. (2020a) (VLASPEC), Lenkić et al.
(2020) (PHIBBS-2) Riechers et al. (2020b) (COLDz), Decarli et al. (2020)
(ASPECS), and Garratt et al. (2021) (UKIDSS-UDS). For ALMACAL-CO maximum and
minimum values of preliminary estimates at $z<2$ are shown. Darker to lighter
shades correspond to 1, 2 and 3$\sigma$ confidence levels by Péroux & Howk
(2020). ALMACAL-abs upper limit in the $z\simeq 0$-2 range is taken by Klitsch
et al. (2019) while Herschel PEP most probable lower limit at z$\simeq 2-$3 is
by Berta et al. (2013).
H2 mass forms in cold high-density regions with rates depending on gas
temperature and ionization state. Therefore, as opposed to $\Omega_{\rm
neutral}$, predicted values for $\Omega_{\rm H_{2}}$ are very sensitive to the
underlying chemical and physical modelling. We neatly see that HI shielding
effects in the HM-HISSmed, P19-HISSmed and FG20-HISSmed runs are responsible
for largely increasing the amount of H2 with respect to the corresponding un-
shielded HM, P19 and FG20 results. Indeed, the $\Omega_{\rm H_{2}}$ parameter
is enhanced by $\sim$1 dex, while the HI mass increases only by a factor up to
2. This reveals the highly non-linear nature of the involved chemical
processes, as to a change of a factor of a few in HI correspond much larger
variations in H2. Clearly, the un-shielded HM, P19 and FG20 models produce
excessive H2 dissociation and are not consistent with the data (the related
$\Omega_{\rm H_{2}}$ values fall down to $\lesssim 2\times 10^{-6}$ at
$z\lesssim 3$, far from observational limits by almost to 2 orders of
magnitude). Despite the inclusion of HI shielding, though, the P19 and FG20
scenarios feature an $\Omega_{\rm H_{2}}$ evolution that flattens at $z\simeq
4$-7, hardly reaching $\Omega_{\rm H_{2}}\sim 10^{-4}$. This is a result of
the earlier onset of the UV background and more powerful rates associated to
these two models at such epochs. Conversely, the HI-shielded HM model produces
values closer to the observed ones by employing only the essential H2
formation channels, which are able to enhance the H2 mass density parameter at
$z\simeq 4$-6 in the presence of a milder UV background. Indeed, the evolution
of the H2 density parameter is affected by the evolution of the
photoionization and photoheating rates adopted in the simulations. The P19
rates evolve smoothly from high redshift until $z\simeq 6$, when they get
enhanced by more than 3 orders of magnitude, while the FG20 rates at $z\simeq
8$ go from zero to levels comparable to those of the P19 model. Although newly
available free charges induce H2 formation, on longer timescales the overall
effect is H2 disruption, in shielded as well as un-shielded cases. The decline
observed in $\Omega_{\rm H_{2}}$ at lower redshift is linked to the stronger
effect of both UV radiation and star formation, which destroy molecules
through photoheating and stellar feedback. The latter is increasingly
effective at later times (see e.g. lower-$z$ trends for un-shielded models),
while the former depends on the details of the UV background at different $z$.
For an earlier and/or stronger UV background (e.g. P19 and FG20) H2
suppression is more significant, albeit modulated by gas self-shielding.
By looking at $\Omega_{\rm H_{2}}$ data in light of our findings, one might
think that the large dispersion in the observations at low redshifts (e.g.
Klitsch et al. 2019; Decarli et al. 2019, 2020; Riechers et al. 2020b; Garratt
et al. 2021) could be related to targeted objects found in conditions that are
heterogeneous in terms of gas shielding, heating, ionization or
photodissociating radiation (hence, the need for large statistical samples).
In practice, the right combination of UV model and HI self-shielding is needed
to match both HI and H2 data. A UV background turning on at $z\gtrsim 8$ seems
to induce an excessive feedback and predict too little cold gas. This suggests
that either UV rates at those times should be lower than those adopted here or
that high-$z$ gas self-shielding should be stronger.

Typically, stochastic conversion of gas into stars is responsible for removing
very dense gas parcels from collapsing material. Nevertheless, $\Omega_{\rm
neutral}$ and $\Omega_{\rm H_{2}}$ parameters might be affected if gas-to-
stars conversion involves particles that are still at relatively low density
(below $\sim$1 particle $\rm cm^{-3}$). In this case the abundance of HI and
H2 cannot grow significantly during their cooling process. To make sure that
gas evolution through the cooling branch is followed properly and the amount
of HI and H2 formed is realistic, we have always used a density threshold of
10 cm-3, as time-dependent non-equilibrium runaway cooling usually takes place
in gas with densities $>0.1\,\rm cm^{-3}$ and temperatures $\sim
10^{2}-10^{4}\,\rm K$. Besides the density threshold, though, gas-to-star
conversion is often assumed to depend on gas density, $\rho$, and a typical
star formation timescale, $t_{\rm sf}$ (see Sect. 2.2.6). This means that the
particle SFR is $\propto\rho/t_{\rm sf}$. Since H2 molecules are the main
drivers of star formation (not the gas tout court) and since we are able to
follow their non-equilibrium chemical abundance precisely, in the following we
employ the local H2 density, $\rho_{\rm H_{2}}$, as derived by our non-
equilibrium calculations to estimate the particle SFR $\propto\rho_{\rm
H_{2}}/t_{\rm sf}$ and related feedback effects (non-equilibrium H2-based star
formation). In Fig. 4 we compare the resulting $\Omega_{\rm neutral}$ and
$\Omega_{\rm H_{2}}$ parameters for such scenario (HM-HISSmed-H2 run) to the
reference results (HM-HISSmed run). Although the evolution of $\Omega_{\rm
neutral}$ and $\Omega_{\rm H_{2}}$ obtained with the H2-based approach is
similar to the reference one, it seems to be a better fit to observations of
both gas components. With respect to the reference trends, $\Omega_{\rm
H_{2}}$ values are slightly higher as consequence of the less effective star
formation feedback. Consequently, HI is more easily converted into H2 and
resulting $\Omega_{\rm neutral}$ values are smaller.

Figure 4: Upper panel: $\Omega_{\rm neutral}$ evolution for the non-
equilibrium H2-based star formation model (red dashed line, HM-HISSmed-H2)
with an HM background and compared to the reference run (solid black line, HM-
HISSmed). Lower panel: Corresponding $\Omega_{\rm H_{2}}$ evolution.
Observational determinations are listed in Tab. 2 (see also text and Fig. 3).

We have additionally investigated the role of individual H2 formation
channels, finding that H- is a fundamental driver at all times, while three-
body interactions are less important on large scales. Nevertheless, they are
useful to reach large H2 fractions of at least $50\%$ (i.e. $\sim 67\%$ of the
available H mass) even during the epoch of reionization and explain the
$50-60\%$ values at $z\sim 4$-6 lately reported by e.g. Dessauges-Zavadsky et
al. (2020) and Tacconi et al. (2018, 2020) (Appendix F). Including gas-grain
processes coupled to the non-equilibrium molecular network and metal spreading
leads to a more complete treatment of cold-gas chemistry, without
substantially affecting our conclusions. Also dust grains can enhance H2
fractions up to $\sim 50\%$ by $z\simeq 7$, in line with observations (for a
more detailed discussion see Appendix C and D), but are effective at
relatively large metallicities. Heating due to cosmic rays formed in star
forming regions is rather localized and, no matter the details of the
ionization rates, its impact remains at a ten-per-cent level (Appendix E).
Feedback due to different SN efficiencies or wind parameterizations impact
$\Omega_{\rm neutral}$ and $\Omega_{\rm H_{2}}$, inhibiting more or less
efficiently the amount of cold neutral and molecular gas available. Material
in the distant IGM/CGM suffers only mild effects from these events (Appendix H
and I).
To summarize, we find that to model accurately cold gas and to reproduce
observed atomic and molecular abundances it is crucial to follow non-
equilibrium chemistry of HI and H2 self-shielded material down to densities of
at least $\sim 10\,\rm cm^{-3}$ and to resolve the gas runaway cooling. In
this way, our results are in better agreement with observations than previous
simulation predictions at high $z$.

###  3.3 Gas depletion times

Gas depletion times are useful tools to investigate the baryon cycle and its
relation to cosmic expansion (often parameterized by the Hubble time, $t_{\rm
H}$). As hot gas is unlikely to collapse, it is natural to focus on cold
neutral gas and on molecular H2-rich gas.
Fig. 5 shows cold-neutral-gas depletion time for different UV models, compared
to observational data and the Hubble time. For all the UV models, depletion
times decrease with decreasing redshift, suggesting that typical cosmic
structures can accumulate cold material that will eventually collapse and form
stars. This qualitative behaviour is basically universal and little dependent
on model details. Quantitatively, results from simulations including HI
shielding are more consistent with the data, with values becoming lower than
$t_{\rm H}$ at $z\lesssim 4$. Un-shielded models with P19 and FG20 background
radiation are only marginally off observational constraints.

Figure 5: Redshift evolution of the neutral-gas depletion times in simulations
with HM, P19 and FG20 (solid, dashed and dot-dot-dot-dashed lines,
respectively) UV backgrounds compared to depletion times with $1\sigma$ errors
(grey points) derived from $\Omega_{\rm neutral}$ observational data and 1, 2
and 3$\sigma$ confidence levels (from darker to lighter shades) (see also Fig.
3). Results are for HI shielded (thick) and HI un-shielded (thin) scenarios.
The Hubble time (grey dotted line) is shown, as well.  Figure 6: Redshift
evolution of the H2 depletion time in simulations with HM, P19 and FG20
(solid, dashed and dot-dot-dot-dashed lines, respectively) UV backgrounds. H2
depletion times within $1\sigma$ range (slanted shapes), 1, 2 and 3$\sigma$
confidence levels (from darker to lighter shades), as well as maximum and
minimum values (grey arrows) are derived from $\Omega_{\rm H_{2}}$
observational data and limits in Tab. 2 and Fig. 3. The Hubble time (thick
grey dotted line) and the dynamical time, $t_{\rm dyn}$ (thin grey dotted
line) are shown, as well. Results refer to both HI shielded (thick) and HI un-
shielded (thin) scenarios.

Fig. 6 shows the corresponding H2 depletion times compared to observational
data. In this case, besides the Hubble time, we also plot the ‘dynamical’
time, $t_{\rm dyn}$, defined as a fraction of 10 per cent of the Hubble time,
$t_{\rm dyn}=0.1\,t_{\rm H}$. This expression can be explained with simple
dynamical arguments (see Appendix J) and is often quoted in observational
studies Lilly et al. (2013); Tacconi et al. (2020) as indicative of the
evolution of a marginally stable rotating disk in which molecular cooling, gas
fragmentation and star formation should take place. In Fig. 6, H2 global
estimates agree well with observational results from alternative galaxy-based
approaches, that suggest a best-fitting trend of $t_{\rm
depl,H_{2}}(z)=(1.6\pm 0.5)(1+z)^{-(0.98\pm 0.12)}\,\rm Gyr$ within 2$\sigma$
error bars (very close to $t_{\rm dyn}$) and data dispersion of almost 1 dex
(see Tacconi et al. 2020, who used existing literature data and ALMA archive
molecular-mass detections up to $z\simeq 5.3$). In line with studies of cosmic
cold gas showing no explicit dependence on galaxy mass, size or environment
for galaxy integrated depletion times Darvish et al. (2018); Tadaki et al.
(2019) and notwithstanding the operational costs to determine molecular
masses, H2 depletion times can track cosmic gas evolution if the star
formation history is well known Tacconi et al. (2020). Despite the relatively
low scatter at $z\simeq 2$-4 (between 0.1 and 1 Gyr), lower-$z$ determinations
have a larger spread. Using available observational data, at $z\simeq 0$-2 we
find $t_{\rm depl,H_{2}}$ comprised between $10\,\rm Gyr$ (ALMACAL-abs upper
limits) and $0.01\,\rm Gyr$ (ASPECS and Herschel PEP lower limits). The latest
UKIDSS data at $z\simeq 2$ suggest values around 0.1-0.3 Gyr. Interestingly,
expected H2 depletion times are always found below the Hubble time and, often,
below the dynamical time. This is a more stringent constraint than the one
derived from neutral-gas depletion times and it supports the idea that cosmic
structures form much beyond redshift $z\simeq 6$. Thus, future observational
efforts in the (sub-)mm regimes leading to larger samples of high-$z$ star
forming galaxies, similar to those recently detected at $z\simeq 6$ Gruppioni
et al. (2020); Talia et al. (2020); Zavala et al. (2021), are required.

##  4 Discussion

In the previous sections we have investigated the role of different physical
and chemical processes that could affect atomic and molecular cold-gas
evolution. We have tested three different UV backgrounds (HM, P19, FG20)
coupled to non-equilibrium chemistry and found that UV radiation from forming
structures has a clear impact on thermodynamics, composition and time
evolution of cosmic gas. This is generally true for all the UV scenarios,
although the implications of each individual model vary remarkably and appear
to be extreme for large UV onset redshifts ($z\gtrsim 8$).

Inclusion of HI and H2 self-shielding is relevant to address the role of UV
photoionization at high (above 1 particle cm-3) and intermediate densities (in
the range between $\sim 10^{-2}$ and 1 particle cm-3). HI shielding has an
important indirect impact during and after reionization, since it
significantly preserves the amount of HI, indispensable for H2 formation, from
UV photoionization. Thus, self-shielding ensures that enough amounts of HI and
H2 are obtained to consistently match the observations. We remind that H2
self-shielding is effective mainly on dense gas and has a small impact on
IGM/CGM-like material (e.g. Petkova & Maio 2012), thus HI self-shielding
results to play a major role in cosmic gas chemistry. The effect of varying H2
shielding prescriptions, instead, is modest (Appendix B). As an example, the
recently updated model by Wolcott-Green et al. (2017) results in excellent
agreement with previous works by Draine & Bertoldi (1996) as long as the
length scale used to derive the local H2 column density is approximated by 1/4
the Jeans length. The effects of alternative choices are small in thin media,
but become important in dense sites (e.g. when estimating molecular cooling of
collapsing gas or critical fluxes of impinging radiation during the formation
of massive black-hole seeds; Luo et al. 2020; Maio et al. 2019). Neglecting HI
shielding corrections at high redshift ($z>5$) leads to a large degree of H
ionization and H2 photodissociation, causing under-estimates of HI and H2
cosmological densities and leading to a trend of the atomic and molecular
density parameters in tension with data even at $z<5$. Furthermore, gas self-
shielding seems crucial to justify the existence of molecular gas at high
redshift, as observed ubiquitously in outflows of massive galaxies at
$z\gtrsim 4$ with loading factors of the order of unity Spilker et al. (2020).

Gas-grain processes have been implemented adopting different values for a
constant grain temperature ($T_{\rm gr}=40$, 75 and 120 K) and a $z$-dependent
one, based either on the CMB scaling or on the dust emissivity index, $\beta$
da Cunha et al. (2021). Since $T_{\rm gr}$ evolution is not very sensitive to
$\beta$, we find that details of its implementation do not have a large impact
on the bulk of structure formation and the overall H2 mass formed. They can be
relevant locally, though, in metal-rich sites (in agreement with expectations
from high-$z$ starburst galaxies; Overzier et al. 2011) and boost molecular
fractions above $50\%$ (Appendix C).

Our results show that the origin of H2 in different epochs is typically linked
to free electrons made available during structure formation that boost the H-
channel in most cosmic environments. Such a relevant role of the H- channel is
in line with conclusions from analytic investigations that stressed the need
to reduce the density of H- to suppress star formation in early mini-haloes
Cen (2017), as well as with recent semi-analytical ISM studies that
highlighted the important role played by H- Thi et al. (2020). This species is
effective in electron-rich, albeit not extremely dense ($\sim 0.1$-1 $\rm
cm^{-3}$) gas, as is the case for most of the cosmic medium during or after
reionization. Free electrons are captured by H atoms via radiative attachment
and increase the amount of H- available for H2 formation. The H${}_{2}^{+}$
channel can provide an additional contribution (of up to $\sim 1/5$), but it
is typically smaller than that from the H- channel, due to the lower
collisional rates (Appendix F). Reactions involving H- and H${}_{2}^{+}$
catalysts to form molecular hydrogen ($\rm H^{-}+H_{2}^{+}\rightarrow
H_{2}+H$) are less likely than H${}_{2}^{+}$ or H- impacts with H and the
resulting H2 enhancement is not substantial, given also their relatively small
abundances. Three-body processes can provide some H2 mass content. Their
effects are subdominant for $\rm\Omega_{\rm neutral}$ and $\Omega_{\rm H_{2}}$
estimates, but help reach high H2 fractions of roughly $50\%$ at $z\sim 4$-6
in metal-poor sites (i.e. values similar to those reached by H2 grain
catalysis in metal-rich sites).

These processes explain well the origin of the large supply of molecular mass
within galaxies at $z\gtrsim 2$ Garratt et al. (2021). In fact, large
observationally inferred H2 fractions at early times Dessauges-Zavadsky et al.
(2020) can be explained either via standard H2 formation channels in metal-
poor/pristine gas or via H2 grain catalysis in metal-rich sites. We warn the
reader that, from a purely statistical point of view, large H2 values,
although detected, might not necessarily be typical for the whole galaxy
population at those epochs. Indeed, there is a clear observational bias
because of which galaxies with large H2 fractions are easier to observe via
their strong CO emission. Previous results from sub-grid models for isolated
galaxies at $z\simeq 2$ focusing on H2 formation on dust grains only
Tomassetti et al. (2015) confirm the leading role of H2 grain catalysis in
enriched environments. The possibility to form H2 molecules via the two
mentioned chemistry paths is also consistent with the lack of observed
environmental (metallicity) dependence up to $z\sim 3.5$ in the assembly of H2
mass (e.g. Darvish et al. 2018; Tadaki et al. 2019). This hides the
fundamental physical difference between the nearby Universe, where H2
formation is led by grain catalysis in polluted media (in agreement with
results of local molecular clouds, Wakelam et al. 2017), and the distant one,
where H2 production is mainly led by the H- channel in shielded
low-$Z$/pristine gas. We warn that, in principle, H2 production could be
enhanced by grain processes at any $z$, as long as interested sites have local
metallicities of $Z\gtrsim 0.1\,\rm Z_{\odot}$. This conclusion is consistent
with detailed ISM analyses by e.g. Gong et al. (2020) (who consider $Z$ within
0.5-2 $Z_{\odot}$) and Hu et al. (2021) (who explore $Z$ within 0.1-3
$Z_{\odot}$).

We have improved upon previous semi-analytical studies by Jasche et al. (2007)
about the possible heating due to cosmic rays at early epochs (Appendix E).
Our constant cosmic-ray ionisation rate per H atoms, $\zeta=10^{-17}\,\rm
s^{-1}$, is consistent with traditional values, with the latest determinations
by Thi et al. (2020) and with the rate assumed in Narayanan & Krumholz (2017),
who correct by a factor of 0.1 the suggested $\zeta\simeq 10^{-16}\,\rm
s^{-1}$ Indriolo & McCall (2012). We additionally consider the calculations by
Padovani et al. (2018) for a density-dependent $\zeta$ and adopt it to assess
the implications on star forming gas. We find that typical values for the
cosmic-ray ionization rate induce heating in addition to the one from
structure formation shocks and cause a deficit of $\Omega_{\rm H_{2}}$ at
early times. In the realistic scenario of cosmic-ray heating linked to star
forming gas parcels, the effects observed at high $z$ are not very prominent
and the resulting $\Omega_{\rm H_{2}}$ and $\Omega_{\rm neutral}$ (or
$\Omega_{\rm HI}$) are close to those without cosmic-ray heating. At lower
$z$, the formation of molecular gas is slightly suppressed and $\Omega_{\rm
H_{2}}$ is a bit smaller than in the reference scenario without cosmic-ray
heating. Overall, the effects are of the order of ten per cent. These results
are consistent with the small (up to a factor of 2) star formation suppression
found in latest studies of isolated disk galaxies Semenov et al. (2021) and
hint at the limited role of cosmic-ray heating in high-$z$ structures.

The additional processes mentioned above (metallicity evolution, metal-
dependent H2 grain catalysis, photoelectric effect, cosmic-ray heating, etc.)
may have an impact on the local environment, depending on the chemical and
thermal state of hosting gas. H2 grain catalysis and photoelectric heating can
boost $\Omega_{\rm H_{2}}$ up to a factor of 3 at $z\simeq 6$, still within
observational constraints Riechers et al. (2020b), while large metallicities
enhance H2 grain catalysis up to $\Omega_{\rm H_{2}}\simeq 2$-$3\times
10^{-4}$ for $Z\sim\,Z_{\odot}$. We note that we have scaled H2 grain reaction
and cooling/heating rates linearly with $Z$, but dust-to-gas ratios could be
steeper in low-$Z$ environments (like in blue compact dwarves, which can be
thought as nearby probes of the physics of distant primeval sources). In such
cases the suggested scaling could be $\propto Z^{1.5}$ Galametz et al. (2011)
or even $\propto Z^{2}$ Herrera-Camus et al. (2012). This would lower the
contribution of gas-grain cooling/heating at low $Z$, imply metallicities
closer to solar to get H2 grain catalysis effective and indicate a much lesser
role of heavy elements at early times. Non-equilibrium results about the
cosmic-ray heating rely on the ionization rate and yield derived from
realistic H, He and H2 gas mixtures. Although improved with respect to
previous works, those quantities are still obtained by neglecting heavy
elements, thus analyses of the energy transfer in more general conditions are
still required. In our study we varied it by up to one order of magnitude and
found small effects for HI and H2 mass densities in star forming regions at
$z<6$. Furthermore, we tested the implications of different IMFs and wind
assumptions on mass density parameters, without finding significant changes
(Appendix I).

The actual role of some physical processes, like photoelectric effect, cosmic
rays or dust cycle in galaxies during the epoch of reionisation, is still
unclear and more data need to be collected Burgarella et al. (2020). Despite
these uncertainties, gas self-shielding, that preserves HI and H2 abundances,
is likely the principal reason for cool molecular gas to extend further out
than dust in observations of star forming galaxies. Indeed, we checked that,
without self-shielding, lower H2 values are obtained even with the inclusion
of H2 grain catalysis. This is in line with recent interpretations of ALMA
observations Calistro Rivera et al. (2018) and should hold a fortiori at high
redshift, when metals and dust grains are expected to be under-abundant. In
this respect observational analyses by Calistro Rivera et al. (2018) support
our conclusions.

Earlier cosmological simulations have highlighted the difficulty in
reproducing the behaviour of chemical density parameters. In particular, the
observed increasing trend of $\rm\Omega_{\rm HI}$ and $\rm\Omega_{\rm
neutral}$ has been challenging to model. Numerical efforts by Nagamine et al.
(2004) using similar star formation and feedback models to Springel &
Hernquist (2003) in a dozen Gadget simulations highlighted the relevant role
of wind feedback and obtained $\rm\Omega_{\rm HI}$ results consistent with
observations at $z=2$-5 with a ”weak” wind speed of 242 km/s. This value is in
line with with our reference wind speed of 350 km/s. The consistency in the
general trend is encouraging and corroborates the results presented in this
manuscript. Nagamine et al. (2004) could not study $\Omega_{\rm H_{2}}$
evolution because of the lack of observational data at the time. Davé et al.
(2013) proposed a treatment for HI and H2 gas content based on H ionization
balance to account for HI and on ISM pressure equilibrium to account for H2999
They compute H2 fractions in star forming particles only. Gas that is not star
forming is assumed to have zero molecular content.. These approximations were
required because their star formation threshold was $0.1\,\rm cm^{-3}$ (too
low to resolve runaway molecular cooling in dense cold gas) and HI/H2 self-
shielding effects were not included. Gas properties were corrected for HI
shielding in post-processing, but they found an HI mass density parameter
decreasing with increasing redshift, with $\rm\Omega_{\rm HI}<10^{-3}$ at
$z<5$ and $\rm\Omega_{\rm HI}\simeq$ a few $10^{-4}$ at $z\simeq 5$. These
predictions led the authors to question the ability of equilibrium-based
models to capture the fundamental processes shaping the HI and H2 behaviour,
as remarked also in Davé et al. (2020). Similarly, low HI mass densities and
decreasing $\rm\Omega_{\rm HI}$ redshift evolution have been obtained by other
independent groups based on semi-analytical models run on top of N-body
simulations (e.g. Spinelli et al. 2020). The deviations of H2 abundances from
equilibrium in cosmic environments highlight the need to follow a time-
dependent non-equilibrium approach to obtain reliable molecular fractions (as
stressed lately by Hu et al. 2017, 2021, and as done throughout this work).
Other popular hydrodynamic simulations (such as Owls or Eagle; Schaye et al.
2015, and references therein) adopt UV rates similar to the ones used here
(HM), but employ low-density thresholds for star formation and do not
partition hydrogen into its ionised, neutral and molecular components (as also
done in the NIHAO project; Obreja et al. 2019). So, post-processing
assumptions on how to distribute these species are needed (Blitz & Rosolowsky
2006 or Gnedin & Kravtsov 2011). Post-processing modelling of hydrogen phases
for the Eagle simulations Rahmati et al. (2015) leads to a roughly flat HI
curve at $z\sim 2-$6, while Illustris-TNG results in 100 Mpc/$h$ side box
feature a non-monotonic behaviour attributed to lack of resolution
Villaescusa-Navarro et al. (2018). Being usually based on local calibrations,
this approach is effective at studying the nearby Universe or enriched
systems, but might present issues at higher redshift or in low-metallicity
environments, such as the whole epoch of reionization and the first few
cosmological Gyrs. In practice, since they are not included in hydro and
chemistry run-time calculations, H2 molecules cannot play any role in gas
cooling. Probably for this reason, semi-analytical computations by Lagos et
al. (2018, and references therein) performed on the outputs of numerical
simulations to estimate HI and H2 fractions with post-processing scaling
relations, find HI densities which are not fully consistent with observations.
Their HI (neutral) mass density parameter decreases from a few $10^{-4}$ (at
$z\simeq 1$) down to $10^{-5}$ (at $z\simeq 3.5$) and also their values never
reach $\Omega_{\rm HI}\sim 10^{-3}$ Rhee et al. (2018). Despite that, Lagos et
al. (2018) can roughly reproduce $\Omega_{\rm H_{2}}$ values at $z\lesssim 3$.
The various Illustris simulations (e.g. Springel et al. 2018, and references
therein) address HI and H2 properties via the empirical model of Gnedin &
Draine (2014, 2016), but do not implement non-equilibrium atomic and molecular
chemistry. When post-processed with semi-analitical models Somerville et al.
(2015) results predict H2 masses that are too low at $z>1$ and cosmic
densities ($\Omega_{\rm H_{2}}$) in tension with ASPECS data, independently
from the adopted H2 partition recipe, as their simulated gas (for which the
star formation density threshold is $\sim 0.1\,\rm cm^{-3}$) is not dense
enough to become molecular Popping et al. (2019). The smaller $\Omega_{\rm
H_{2}}$ values found by Popping et al. (2019) might also be linked to the
chosen UV background by Faucher-Giguère et al. (2009), that (similarly to
Faucher-Giguère 2020, FG20) assumes an earlier onset and hence causes a longer
photoheating injection with respect to the HM model. Although this latter is
adopted in Lagos et al. (2018), their latest semi-analytical estimates Lagos
et al. (2020) still find $\rho_{\rm H_{2}}$ values systematically below the
observational constraints. Other recent studies (such as Gjergo et al. 2018,
2020, Vogelsberger et al. 2019, Graziani et al. 2020, Granato et al. 2021)
have performed dust mass modelling from gas mass metallicities within three-
dimensional numerical simulations, however their implementations are not
coupled to non-equilibrium H2 chemistry. Finally, recent H2 modelling by
Schäbe et al. (2020) that includes H2 formation via dust grain catalysis only
in a 12 Mpc/$h$ side box has suggested the existence of undetected molecular
amounts at high redshift, since their predicted H2 mass densities at $z\simeq
4-$6 were higher than current observational determinations (as in Fig. 9).

We note that, as in any such work, numerical accuracy is important, but the
different star formation schemes and feedback implementations adopted often
have a larger impact. Feedback efficiency and star formation timescales
modulate the whole evolution of H2 mass. Changing the shape of the IMF from
Salpeter to Chabrier alters slightly star formation related quantities (in
agreement with e.g. Cassarà et al. 2013), while varying wind velocities Hassan
& Gronke (2021) or SN efficiency has an impact of up to a dex when these
parameters are varied by a factor of 5-10, independently from other model
assumptions Steinwandel et al. (2020). For HI the impact is less dramatic and
results vary by less than half dex. Chemistry rates and yields are still
affected by factor-of-two uncertainties (for examples compare H2 grain
catalysis rates in Hollenbach & McKee 1979; Tielens & Hollenbach 1985; Cazaux
& Tielens 2004; Tomassetti et al. 2015) and this influences resulting
abundance ratios Maio & Tescari (2015); Ma et al. (2015, 2017b); Valentini et
al. (2019a).

From an observational point of view, HI mass densities are relatively well
established up to $z\lesssim 6$ with errors below a factor of 2. Moreover,
exploiting [C II] emission in high-$z$ enriched environments (as lately done
by Heintz et al. 2021) would possibly open the door to observing $\Omega_{\rm
neutral}$ at $z>6$, where our results predict different scenarios depending on
the chosen UV background. On the other hand, estimates of the reservoir of H2
masses in star forming environments from [C II] or CO observational data are
challenging and usually adopted conversion factors bear uncertainties and
biases Madden et al. (2020); Gong et al. (2020); Breysse et al. (2021). The [C
II] 158 micron line is a workhorse for (sub-)mm observations and is resolved
on kpc-scale by ALMA Rybak et al. (2019). Since extended [C II] halos are
found more and more commonly in high-$z$ galaxies (early star forming sources
of $10^{11}$-$10^{12}$ M⊙ fit with classical ‘normal’ star forming galaxies;
Ginolfi et al. 2019; Hodge & da Cunha 2020) and since [C II] emissions tracing
a combination of molecular, ionised and atomic gas are still very little
explored, future assessments of H2 via [C II] could improve significantly the
current picture for $\Omega_{\rm H_{2}}$.

Throughout this work we have assumed a $\rm\Lambda$CDM cosmology. We briefly
point out that in some particular cases (such as for warm dark matter, Maio &
Viel 2015) the background cosmological model influences the baryonic-structure
evolution, mostly at higher redshift. This plays a role during the beginning
of the onset of star formation, but in the epochs of interest here baryon
evolution tends to dominate and alleviate any discrepancy due either to
alternative cosmologies Maio et al. (2006) and matter non-Gaussianities Maio &
Iannuzzi (2011); Maio (2011); Maio et al. (2012); Maio & Khochfar (2012) or to
a different dark-matter nature Maio & Viel (2015). Finally, the baryon census
($\Omega_{0,\rm{b}}$) is crucial for assessing atomic and molecular gas
evolution. Alternative estimates of $\Omega_{0,\rm{b}}$ have been proposed
(using e.g. FRBs, Macquart et al. 2020) and, although uncertainties still
persist, they are consistent with the cosmic microwave background and the Big
Bang nucleosynthesis.

##  5 Conclusions

We have quantified the evolution of cold atomic and molecular cosmic gas, from
the epoch of reionization down to $z\sim 2$. We have performed cosmological
N-body hydrodynamic chemistry simulations to model the origin and evolution of
cold-gas density parameters and timescales. Our implementation follows time-
dependent non-equilibrium atomic and molecular abundance calculations and
considers the impact of different UV backgrounds, as well as a number of
physical processes that could affect the build-up of H2 and neutral-gas mass.
Star formation in dense HI and H2 self-shielded regions, stellar feedback and
production of heavy-elements from stars with different masses and
metallicities are also accounted for. We constrain the effects of the major
processes involved in cold-gas evolution and interpret state-of-the-art
observations in the (sub-)mm and IR ranges. This is among the first studies
addressing cosmic HI and H2 evolution at these epochs by three-dimensional
hydrodynamic simulations and including detailed non-equilibrium chemistry
networks. In summary, we find that:

  * –

the evolution of cold atomic and molecular gas obtained with basic time-
dependent non-equilibrium chemical reactions in dense gas is broadly
consistent with the latest HI and H2 observational data once HI and H2 self-
shielding from UV radiation is taken into account within the chemical network;

  * –

the resulting mass of neutral gas slowly increases with redshift up to
$\Omega_{\rm neutral}\sim 2\times 10^{-3}$ at $z\simeq 6$, while H2 features a
plateau at $\Omega_{\rm H_{2}}\sim 10^{-4}$ at $z\gtrsim 2$;

  * –

details of the UV background model adopted (HM, P19 and FG20) have a modest
impact on $\Omega_{\rm neutral}$ at $z\lesssim 6$, while they are relevant for
$\Omega_{\rm H_{2}}$: the onset epoch of the UV background in particular seems
to be crucial, suggesting that early ($z\gtrsim 8$) onsets could be in tension
with H2 data and requiring either lower UV rates or larger self-shielding at
those times;

  * –

resolving dense shielded gas and adopting a high star formation density
threshold is necessary to properly estimate the amount of cold gas and to
obtain results in line with HI and H2 observational data;

  * –

our non-equilibrium H2-based star formation model leads to results that are in
slightly better agreement with available HI and H2 determinations;

  * –

H2 evolution is globally driven by the H- channel, with a relevant
contribution from H2 grain catalysis at high metallicities and on local
scales, while other formation channels are negligible: this explains the
observationally inferred lack of environmental dependence in H2 mass build-up
at $z\lesssim 3.5$;

  * –

cosmic-ray heating in star forming regions, quantified in accordance to
different models for the cosmic-ray ionization rate, affects the H2 abundance
at ten-per-cent level, globally;

  * –

wind feedback impacts results mostly at lower redshift, while SN explosions
affect early star forming regimes, with effects from IMF variations remaining
usually small at all times;

  * –

recently detected molecular fractions of $\sim 50\%$ at $z\simeq 4$-6 can be
explained by essential non-equilibrium H2 formation channels in metal-poor
media or H2 grain catalysis in metal-rich gas: in this latter case, a precise
determination of grain temperature has little relevance and is generally
subdominant to metallicity effects;

  * –

despite the insurgence of UV radiation, H2 depletion times reach values lower
than the Hubble time, even comparable to the dynamical time, within the first
half Gyr, meaning that cosmic gas is able to collapse already at such early
epochs.

The ability to form significant amounts of H2 and to reach short depletion
times during the epoch of reionization means that future discoveries of
molecular-rich star forming galaxies at early times will be possible in the
next years. New data from upcoming international facilities, such as JWST
Gardner & JWST Science Working Group (2009), ELT Puech et al. (2010), SKA
Koopmans et al. (2015), Roman/WFIRST Whalen et al. (2013) or ATHENA Nandra et
al. (2013) will be decisive to shed light on the still unanswered questions
about the origin of chemical species in the cosmic gas and the whole baryon
cycle.

###### Acknowledgements.

We are grateful to the anonymous referee for his/her careful reading and
suggestions that improved the original manuscript. We are thankful to F.
Haardt, S. Borgani, L. Tornatore and R. Szakacs for their constructive
comments. U.M. is indebted to L. Böss and U. Steinwandel for insightful
conversations and to M. Padovani for kindly sharing the full table of his
ionisation rate calculations used throughout this work. We thank A. Hamanowicz
for sharing with us preliminary ALMACAL-CO results. We also thank the Max
Planck Computation and Data Facility (MPCDF) of the Max Planck Society, where
simulation runs and analyses have been performed. We acknowledge the NASA
Astrophysics Data System and the JSTOR archive for their bibliographic tools.

## References

  * Abel et al. (1997) Abel, T., Anninos, P., Zhang, Y., & Norman, M. L. 1997, New A, 2, 181
  * Abel et al. (2002) Abel, T., Bryan, G. L., & Norman, M. L. 2002, Science, 295, 93
  * Anninos et al. (1997) Anninos, P., Zhang, Y., Abel, T., & Norman, M. L. 1997, New A, 2, 209
  * Bakes & Tielens (1994) Bakes, E. L. O. & Tielens, A. G. G. M. 1994, ApJ, 427, 822
  * Bera et al. (2019) Bera, A., Kanekar, N., Chengalur, J. N., & Bagla, J. S. 2019, ApJ, 882, L7
  * Berta et al. (2013) Berta, S., Lutz, D., Nordon, R., et al. 2013, A&A, 555, L8
  * Binney & Tremaine (2008) Binney, J. & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton University Press)
  * Blitz & Rosolowsky (2006) Blitz, L. & Rosolowsky, E. 2006, ApJ, 650, 933
  * Bolatto et al. (2013) Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013, ARA&A, 51, 207
  * Boogaard et al. (2021) Boogaard, L. A., Bouwens, R. J., Riechers, D., et al. 2021, arXiv e-prints, arXiv:2105.12489
  * Breysse et al. (2021) Breysse, P. C., Yang, S., Somerville, R. S., et al. 2021, arXiv e-prints, arXiv:2106.14904
  * Bromm & Yoshida (2011) Bromm, V. & Yoshida, N. 2011, ARA&A, 49, 373
  * Bron et al. (2014) Bron, E., Le Bourlot, J., & Le Petit, F. 2014, A&A, 569, A100
  * Bruzual & Charlot (2003) Bruzual, G. & Charlot, S. 2003, MNRAS, 344, 1000
  * Burgarella et al. (2020) Burgarella, D., Nanni, A., Hirashita, H., et al. 2020, A&A, 637, A32
  * Calistro Rivera et al. (2018) Calistro Rivera, G., Hodge, J. A., Smail, I., et al. 2018, ApJ, 863, 56
  * Campisi et al. (2011) Campisi, M. A., Maio, U., Salvaterra, R., & Ciardi, B. 2011, MNRAS, 416, 2760
  * Cassarà et al. (2013) Cassarà, L. P., Piovan, L., Weiss, A., Salaris, M., & Chiosi, C. 2013, MNRAS, 436, 2824
  * Cazaux & Tielens (2004) Cazaux, S. & Tielens, A. G. G. M. 2004, ApJ, 604, 222
  * Cen (2017) Cen, R. 2017, MNRAS, 465, L69
  * Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763
  * Chen et al. (2021) Chen, Q., Meyer, M., Popping, A., et al. 2021, arXiv e-prints, arXiv:2104.07973
  * Chomiuk & Povich (2011) Chomiuk, L. & Povich, M. S. 2011, AJ, 142, 197
  * Christiansen et al. (2020) Christiansen, J. F., Davé, R., Sorini, D., & Anglés-Alcázar, D. 2020, MNRAS, 499, 2617
  * Crighton et al. (2015) Crighton, N. H. M., Murphy, M. T., Prochaska, J. X., et al. 2015, MNRAS, 452, 217
  * da Cunha et al. (2013) da Cunha, E., Groves, B., Walter, F., et al. 2013, ApJ, 766, 13
  * da Cunha et al. (2021) da Cunha, E., Hodge, J. A., Casey, C. M., et al. 2021, arXiv e-prints, arXiv:2106.08566
  * Darvish et al. (2018) Darvish, B., Scoville, N. Z., Martin, C., et al. 2018, ApJ, 860, 111
  * Davé et al. (2019) Davé, R., Anglés-Alcázar, D., Narayanan, D., et al. 2019, MNRAS, 486, 2827
  * Davé et al. (2020) Davé, R., Crain, R. A., Stevens, A. R. H., et al. 2020, MNRAS, 497, 146
  * Davé et al. (2013) Davé, R., Katz, N., Oppenheimer, B. D., Kollmeier, J. A., & Weinberg, D. H. 2013, MNRAS, 434, 2645
  * Davé et al. (2017) Davé, R., Rafieferantsoa, M. H., & Thompson, R. J. 2017, MNRAS, 471, 1671
  * Decarli et al. (2020) Decarli, R., Aravena, M., Boogaard, L., et al. 2020, ApJ, 902, 110
  * Decarli et al. (2019) Decarli, R., Walter, F., Gónzalez-López, J., et al. 2019, ApJ, 882, 138
  * Delhaize et al. (2013) Delhaize, J., Meyer, M. J., Staveley-Smith, L., & Boyle, B. J. 2013, MNRAS, 433, 1398
  * Dessauges-Zavadsky et al. (2020) Dessauges-Zavadsky, M., Ginolfi, M., Pozzi, F., et al. 2020, A&A, 643, A5
  * Dickman et al. (1986) Dickman, R. L., Snell, R. L., & Schloerb, F. P. 1986, ApJ, 309, 326
  * Dolag et al. (2009) Dolag, K., Borgani, S., Murante, G., & Springel, V. 2009, MNRAS, 399, 497
  * Donnari et al. (2019) Donnari, M., Pillepich, A., Nelson, D., et al. 2019, MNRAS, 485, 4817
  * Draine & Bertoldi (1996) Draine, B. T. & Bertoldi, F. 1996, ApJ, 468, 269
  * Draine & Lee (1984) Draine, B. T. & Lee, H. M. 1984, ApJ, 285, 89
  * Draine & Li (2007) Draine, B. T. & Li, A. 2007, ApJ, 657, 810
  * Draine & Sutin (1987) Draine, B. T. & Sutin, B. 1987, ApJ, 320, 803
  * Duley (1996) Duley, W. W. 1996, MNRAS, 279, 591
  * Dwek & Smith (1996) Dwek, E. & Smith, R. K. 1996, ApJ, 459, 686
  * Faisst et al. (2020) Faisst, A. L., Fudamoto, Y., Oesch, P. A., et al. 2020, MNRAS, 498, 4192
  * Faucher-Giguère (2020) Faucher-Giguère, C.-A. 2020, MNRAS, 493, 1614
  * Faucher-Giguère et al. (2009) Faucher-Giguère, C.-A., Lidz, A., Zaldarriaga, M., & Hernquist, L. 2009, ApJ, 703, 1416
  * Fletcher et al. (2021) Fletcher, T. J., Saintonge, A., Soares, P. S., & Pontzen, A. 2021, MNRAS, 501, 411
  * Flower & Harris (2007) Flower, D. R. & Harris, G. J. 2007, MNRAS, 377, 705
  * Forrey (2013) Forrey, R. C. 2013, ApJ, 773, L25
  * Furlong et al. (2015) Furlong, M., Bower, R. G., Theuns, T., et al. 2015, MNRAS, 450, 4486
  * Galametz et al. (2011) Galametz, M., Madden, S. C., Galliano, F., et al. 2011, A&A, 532, A56
  * Galli & Palla (1998) Galli, D. & Palla, F. 1998, A&A, 335, 403
  * Gardner & JWST Science Working Group (2009) Gardner, J. P. & JWST Science Working Group. 2009, in Bulletin of the American Astronomical Society, Vol. 41, American Astronomical Society Meeting Abstracts #213, 254
  * Garratt et al. (2021) Garratt, T. K., Coppin, K. E. K., Geach, J. E., et al. 2021, arXiv e-prints, arXiv:2103.08613
  * Genel et al. (2014) Genel, S., Vogelsberger, M., Springel, V., et al. 2014, MNRAS, 445, 175
  * Genzel et al. (2012) Genzel, R., Tacconi, L. J., Combes, F., et al. 2012, ApJ, 746, 69
  * Ginolfi et al. (2019) Ginolfi, M., Schneider, R., Valiante, R., et al. 2019, MNRAS, 483, 1256
  * Girichidis et al. (2020) Girichidis, P., Offner, S. S. R., Kritsuk, A. G., et al. 2020, Space Sci. Rev., 216, 68
  * Gjergo et al. (2018) Gjergo, E., Granato, G. L., Murante, G., et al. 2018, MNRAS, 479, 2588
  * Gjergo et al. (2020) Gjergo, E., Palla, M., Matteucci, F., et al. 2020, MNRAS, 493, 2782
  * Glassgold et al. (2012) Glassgold, A. E., Galli, D., & Padovani, M. 2012, ApJ, 756, 157
  * Glover & Clark (2012) Glover, S. C. O. & Clark, P. C. 2012, MNRAS, 421, 9
  * Glover & Savin (2009) Glover, S. C. O. & Savin, D. W. 2009, MNRAS, 393, 911
  * Gnedin & Draine (2014) Gnedin, N. Y. & Draine, B. T. 2014, ApJ, 795, 37
  * Gnedin & Draine (2016) Gnedin, N. Y. & Draine, B. T. 2016, ApJ, 830, 54
  * Gnedin & Kravtsov (2011) Gnedin, N. Y. & Kravtsov, A. V. 2011, ApJ, 728, 88
  * Goldsmith & Langer (1978) Goldsmith, P. F. & Langer, W. D. 1978, ApJ, 222, 881
  * Goldsmith et al. (2011) Goldsmith, P. F., Liseau, R., Bell, T. A., et al. 2011, ApJ, 737, 96
  * Gong et al. (2020) Gong, M., Ostriker, E. C., Kim, C.-G., & Kim, J.-G. 2020, ApJ, 903, 142
  * Granato et al. (2021) Granato, G. L., Ragone-Figueroa, C., Taverna, A., et al. 2021, MNRAS, 503, 511
  * Graziani et al. (2020) Graziani, L., Schneider, R., Ginolfi, M., et al. 2020, MNRAS, 494, 1071
  * Gruppioni et al. (2020) Gruppioni, C., Béthermin, M., Loiacono, F., et al. 2020, A&A, 643, A8
  * Gry et al. (2002) Gry, C., Boulanger, F., Nehmé, C., et al. 2002, A&A, 391, 675
  * Haardt & Madau (1996) Haardt, F. & Madau, P. 1996, ApJ, 461, 20
  * Haardt & Madau (2012) Haardt, F. & Madau, P. 2012, ApJ, 746, 125
  * Habing (1968) Habing, H. J. 1968, Bull. Astron. Inst. Netherlands, 19, 421
  * Hamanowicz et al. (2021) Hamanowicz, A., Péroux, C., Zwaan, M., & et al. 2021, submitted
  * Hartwig et al. (2015) Hartwig, T., Glover, S. C. O., Klessen, R. S., Latif, M. A., & Volonteri, M. 2015, MNRAS, 452, 1233
  * Hassan & Gronke (2021) Hassan, S. & Gronke, M. 2021, ApJ, 908, 219
  * Heintz et al. (2021) Heintz, K. E., Watson, D., Oesch, P., Narayanan, D., & Madden, S. C. 2021, arXiv e-prints, arXiv:2108.13442
  * Herrera-Camus et al. (2012) Herrera-Camus, R., Fisher, D. B., Bolatto, A. D., et al. 2012, ApJ, 752, 112
  * Hodge & da Cunha (2020) Hodge, J. A. & da Cunha, E. 2020, Royal Society Open Science, 7, 200556
  * Hollenbach et al. (2012) Hollenbach, D., Kaufman, M. J., Neufeld, D., Wolfire, M., & Goicoechea, J. R. 2012, ApJ, 754, 105
  * Hollenbach & McKee (1979) Hollenbach, D. & McKee, C. F. 1979, ApJS, 41, 555
  * Hollenbach & McKee (1989) Hollenbach, D. & McKee, C. F. 1989, ApJ, 342, 306
  * Hollenbach & Tielens (1999) Hollenbach, D. J. & Tielens, A. G. G. M. 1999, Reviews of Modern Physics, 71, 173
  * Hoppmann et al. (2015) Hoppmann, L., Staveley-Smith, L., Freudling, W., et al. 2015, MNRAS, 452, 3726
  * Hu et al. (2017) Hu, C.-Y., Naab, T., Glover, S. C. O., Walch, S., & Clark, P. C. 2017, MNRAS, 471, 2151
  * Hu et al. (2021) Hu, C.-Y., Sternberg, A., & van Dishoeck, E. F. 2021, arXiv e-prints, arXiv:2103.03889
  * Hu et al. (2019) Hu, W., Hoppmann, L., Staveley-Smith, L., et al. 2019, MNRAS, 489, 1619
  * Hughes et al. (2017) Hughes, T. M., Ibar, E., Villanueva, V., et al. 2017, MNRAS, 468, L103
  * Hunt et al. (2020) Hunt, L. K., Tortora, C., Ginolfi, M., & Schneider, R. 2020, A&A, 643, A180
  * Indriolo et al. (2007) Indriolo, N., Geballe, T. R., Oka, T., & McCall, B. J. 2007, ApJ, 671, 1736
  * Indriolo & McCall (2012) Indriolo, N. & McCall, B. J. 2012, ApJ, 745, 91
  * Jasche et al. (2007) Jasche, J., Ciardi, B., & Enßlin, T. A. 2007, MNRAS, 380, 417
  * Jones et al. (2018) Jones, M. G., Haynes, M. P., Giovanelli, R., & Moorman, C. 2018, MNRAS, 477, 2
  * Jones & Bate (2018) Jones, M. O. & Bate, M. R. 2018, MNRAS, 480, 2562
  * Jura (1975a) Jura, M. 1975a, ApJ, 197, 581
  * Jura (1975b) Jura, M. 1975b, ApJ, 197, 575
  * Katz (1992) Katz, N. 1992, ApJ, 391, 502
  * Khusanova et al. (2020) Khusanova, Y., Béthermin, M., Le Fèvre, O., et al. 2020, arXiv e-prints, arXiv:2007.08384
  * Klitsch et al. (2019) Klitsch, A., Zwaan, M. A., Péroux, C., et al. 2019, MNRAS, 482, L65
  * Komatsu et al. (2011) Komatsu, E., Smith, K. M., Dunkley, J., et al. 2011, ApJS, 192, 18
  * Koopmans et al. (2015) Koopmans, L., Pritchard, J., Mellema, G., & and 41 co-authors. 2015, Advancing Astrophysics with the Square Kilometre Array (AASKA14), 1
  * Krumholz & Gnedin (2011) Krumholz, M. R. & Gnedin, N. Y. 2011, ApJ, 729, 36
  * Lagos et al. (2020) Lagos, C. d. P., da Cunha, E., Robotham, A. S. G., et al. 2020, MNRAS, 499, 1948
  * Lagos et al. (2018) Lagos, C. d. P., Tobar, R. J., Robotham, A. S. G., et al. 2018, MNRAS, 481, 3573
  * Leitherer et al. (1999) Leitherer, C., Schaerer, D., Goldader, J. D., et al. 1999, ApJS, 123, 3
  * Lenkić et al. (2020) Lenkić, L., Bolatto, A. D., Förster Schreiber, N. M., et al. 2020, AJ, 159, 190
  * Leroy et al. (2013) Leroy, A. K., Walter, F., Sandstrom, K., et al. 2013, AJ, 146, 19
  * Licquia & Newman (2015) Licquia, T. C. & Newman, J. A. 2015, ApJ, 806, 96
  * Lilly et al. (2013) Lilly, S. J., Carollo, C. M., Pipino, A., Renzini, A., & Peng, Y. 2013, ApJ, 772, 119
  * Luo et al. (2020) Luo, Y., Shlosman, I., Nagamine, K., & Fang, T. 2020, MNRAS, 492, 4917
  * Ma et al. (2015) Ma, Q., Maio, U., Ciardi, B., & Salvaterra, R. 2015, MNRAS, 449, 3006
  * Ma et al. (2017a) Ma, Q., Maio, U., Ciardi, B., & Salvaterra, R. 2017a, MNRAS, 466, 1140
  * Ma et al. (2017b) Ma, Q., Maio, U., Ciardi, B., & Salvaterra, R. 2017b, MNRAS, 472, 3532
  * Macquart et al. (2020) Macquart, J. P., Prochaska, J. X., McQuinn, M., et al. 2020, Nature, 581, 391
  * Madden et al. (2020) Madden, S. C., Cormier, D., Hony, S., et al. 2020, A&A, 643, A141
  * Madden et al. (1997) Madden, S. C., Poglitsch, A., Geis, N., Stacey, G. J., & Townes, C. H. 1997, ApJ, 483, 200
  * Maio (2011) Maio, U. 2011, Classical and Quantum Gravity, 28, 225015
  * Maio (2021) Maio, U. 2021, arXiv e-prints, arXiv:2111.01804
  * Maio et al. (2019) Maio, U., Borgani, S., Ciardi, B., & Petkova, M. 2019, PASA, 36, e020
  * Maio et al. (2010) Maio, U., Ciardi, B., Dolag, K., Tornatore, L., & Khochfar, S. 2010, MNRAS, 407, 1003
  * Maio et al. (2009) Maio, U., Ciardi, B., Yoshida, N., Dolag, K., & Tornatore, L. 2009, A&A, 503, 25
  * Maio et al. (2007) Maio, U., Dolag, K., Ciardi, B., & Tornatore, L. 2007, MNRAS, 379, 963
  * Maio et al. (2006) Maio, U., Dolag, K., Meneghetti, M., et al. 2006, MNRAS, 373, 869
  * Maio & Iannuzzi (2011) Maio, U. & Iannuzzi, F. 2011, MNRAS, 415, 3021
  * Maio & Khochfar (2012) Maio, U. & Khochfar, S. 2012, MNRAS, 421, 1113
  * Maio et al. (2011a) Maio, U., Khochfar, S., Johnson, J. L., & Ciardi, B. 2011a, MNRAS, 414, 1145
  * Maio et al. (2011b) Maio, U., Koopmans, L. V. E., & Ciardi, B. 2011b, MNRAS, 412, L40
  * Maio et al. (2016) Maio, U., Petkova, M., De Lucia, G., & Borgani, S. 2016, MNRAS, 460, 3733
  * Maio et al. (2012) Maio, U., Salvaterra, R., Moscardini, L., & Ciardi, B. 2012, MNRAS, 426, 2078
  * Maio & Tescari (2015) Maio, U. & Tescari, E. 2015, MNRAS, 453, 3798
  * Maio & Viel (2015) Maio, U. & Viel, M. 2015, MNRAS, 446, 2760
  * Maloney & Black (1988) Maloney, P. & Black, J. H. 1988, ApJ, 325, 389
  * Mancini et al. (2016) Mancini, M., Schneider, R., Graziani, L., et al. 2016, MNRAS, 462, 3130
  * Mancini et al. (2015) Mancini, M., Schneider, R., Graziani, L., et al. 2015, MNRAS, 451, L70
  * Muratov et al. (2013) Muratov, A. L., Gnedin, O. Y., Gnedin, N. Y., & Zemp, M. 2013, ApJ, 772, 106
  * Nagamine et al. (2004) Nagamine, K., Springel, V., & Hernquist, L. 2004, MNRAS, 348, 421
  * Nandra et al. (2013) Nandra, K., Barret, D., Barcons, X., et al. 2013, ArXiv:1306.2307
  * Narayanan & Krumholz (2017) Narayanan, D. & Krumholz, M. R. 2017, MNRAS, 467, 50
  * Noterdaeme et al. (2012) Noterdaeme, P., Petitjean, P., Carithers, W. C., et al. 2012, A&A, 547, L1
  * Noterdaeme et al. (2009) Noterdaeme, P., Petitjean, P., Ledoux, C., & Srianand, R. 2009, A&A, 505, 1087
  * Obreja et al. (2019) Obreja, A., Macciò, A. V., Moster, B., et al. 2019, MNRAS, 490, 1518
  * Omukai et al. (2005) Omukai, K., Tsuribe, T., Schneider, R., & Ferrara, A. 2005, ApJ, 626, 627
  * Opal et al. (1971) Opal, C. B., Peterson, W. K., & Beatty, E. C. 1971, J. Chem. Phys., 55, 4100
  * Orel (1987) Orel, A. E. 1987, J. Chem. Phys., 87, 314
  * O’Shea & Norman (2008) O’Shea, B. W. & Norman, M. L. 2008, ApJ, 673, 14
  * Overzier et al. (2011) Overzier, R. A., Heckman, T. M., Wang, J., et al. 2011, ApJ, 726, L7
  * Padovani et al. (2018) Padovani, M., Galli, D., Ivlev, A. V., Caselli, P., & Ferrara, A. 2018, A&A, 619, A144
  * Padovani & Matteucci (1993) Padovani, P. & Matteucci, F. 1993, ApJ, 416, 26
  * Pak et al. (1998) Pak, S., Jaffe, D. T., van Dishoeck, E. F., Johansson, L. E. B., & Booth, R. S. 1998, ApJ, 498, 735
  * Palla et al. (1983) Palla, F., Salpeter, E. E., & Stahler, S. W. 1983, ApJ, 271, 632
  * Pallottini et al. (2017) Pallottini, A., Ferrara, A., Bovino, S., et al. 2017, MNRAS, 471, 4128
  * Peebles & Dicke (1968) Peebles, P. J. E. & Dicke, R. H. 1968, ApJ, 154, 891
  * Péroux & Howk (2020) Péroux, C. & Howk, J. C. 2020, ARA&A, 58, 363
  * Petkova & Maio (2012) Petkova, M. & Maio, U. 2012, MNRAS, 422, 3067
  * Petkova & Springel (2011) Petkova, M. & Springel, V. 2011, MNRAS, 415, 3731
  * Planck Collaboration (2020) Planck Collaboration. 2020, A&A, 641, A6
  * Ploeckinger & Schaye (2020) Ploeckinger, S. & Schaye, J. 2020, MNRAS, 497, 4857
  * Popping et al. (2019) Popping, G., Pillepich, A., Somerville, R. S., et al. 2019, ApJ, 882, 137
  * Puchwein et al. (2019) Puchwein, E., Haardt, F., Haehnelt, M. G., & Madau, P. 2019, MNRAS, 485, 47
  * Puech et al. (2010) Puech, M., Rosati, P., Toft, S., et al. 2010, MNRAS, 402, 903
  * Rahmati et al. (2013) Rahmati, A., Pawlik, A. H., Raičević, M., & Schaye, J. 2013, MNRAS, 430, 2427
  * Rahmati et al. (2015) Rahmati, A., Schaye, J., Bower, R. G., et al. 2015, MNRAS, 452, 2034
  * Rao et al. (2017) Rao, S. M., Turnshek, D. A., Sardane, G. M., & Monier, E. M. 2017, MNRAS, 471, 3428
  * Rhee et al. (2018) Rhee, J., Lah, P., Briggs, F. H., et al. 2018, MNRAS, 473, 1879
  * Rhee et al. (2016) Rhee, J., Lah, P., Chengalur, J. N., Briggs, F. H., & Colless, M. 2016, MNRAS, 460, 2675
  * Rhee et al. (2013) Rhee, J., Zwaan, M. A., Briggs, F. H., et al. 2013, MNRAS, 435, 2693
  * Riechers et al. (2020a) Riechers, D. A., Boogaard, L. A., Decarli, R., et al. 2020a, ApJ, 896, L21
  * Riechers et al. (2020b) Riechers, D. A., Hodge, J. A., Pavesi, R., et al. 2020b, ApJ, 895, 81
  * Rodighiero et al. (2019) Rodighiero, G., Enia, A., Delvecchio, I., et al. 2019, ApJ, 877, L38
  * Röllig et al. (2006) Röllig, M., Ossenkopf, V., Jeyakumar, S., Stutzki, J., & Sternberg, A. 2006, A&A, 451, 917
  * Rybak et al. (2019) Rybak, M., Calistro Rivera, G., Hodge, J. A., et al. 2019, ApJ, 876, 112
  * Salpeter (1955) Salpeter, E. E. 1955, ApJ, 121, 161
  * Sánchez et al. (2016) Sánchez, S. F., García-Benito, R., Zibetti, S., et al. 2016, A&A, 594, A36
  * Saslaw & Zipoy (1967) Saslaw, W. C. & Zipoy, D. 1967, Nature, 216, 976
  * Schäbe et al. (2020) Schäbe, A., Romano-Díaz, E., Porciani, C., Ludlow, A. D., & Tomassetti, M. 2020, MNRAS, 497, 5008
  * Schaye et al. (2015) Schaye, J., Crain, R. A., Bower, R. G., et al. 2015, MNRAS, 446, 521
  * Scoville & Solomon (1975) Scoville, N. Z. & Solomon, P. M. 1975, ApJ, 199, L105
  * Semenov et al. (2021) Semenov, V. A., Kravtsov, A. V., & Caprioli, D. 2021, ApJ, 910, 126
  * Shaw et al. (2009) Shaw, G., Ferland, G. J., Henney, W. J., et al. 2009, ApJ, 701, 677
  * Sillero et al. (2021) Sillero, E., Tissera, P. B., Lambas, D. G., et al. 2021, arXiv e-prints, arXiv:2104.04357
  * Skinner & Wise (2020) Skinner, D. & Wise, J. H. 2020, MNRAS, 492, 4386
  * Somerville et al. (2015) Somerville, R. S., Popping, G., & Trager, S. C. 2015, MNRAS, 453, 4337
  * Spilker et al. (2020) Spilker, J. S., Phadke, K. A., Aravena, M., et al. 2020, ApJ, 905, 85
  * Spinelli et al. (2020) Spinelli, M., Zoldan, A., De Lucia, G., Xie, L., & Viel, M. 2020, MNRAS, 493, 5434
  * Spitzer & Scott (1969) Spitzer, Lyman, J. & Scott, E. H. 1969, ApJ, 158, 161
  * Springel (2005) Springel, V. 2005, MNRAS, 364, 1105
  * Springel & Hernquist (2003) Springel, V. & Hernquist, L. 2003, MNRAS, 339, 289
  * Springel et al. (2018) Springel, V., Pakmor, R., Pillepich, A., et al. 2018, MNRAS, 475, 676
  * Stahler & Palla (2004) Stahler, S. W. & Palla, F. 2004, The Formation of Stars (Weinheim: Wiley-VCH)
  * Steinwandel et al. (2020) Steinwandel, U. P., Moster, B. P., Naab, T., Hu, C.-Y., & Walch, S. 2020, MNRAS, 495, 1035
  * Sternberg et al. (2014) Sternberg, A., Le Petit, F., Roueff, E., & Le Bourlot, J. 2014, ApJ, 790, 10
  * Tacconi et al. (2018) Tacconi, L. J., Genzel, R., Saintonge, A., et al. 2018, ApJ, 853, 179
  * Tacconi et al. (2020) Tacconi, L. J., Genzel, R., & Sternberg, A. 2020, ARA&A, 58, 157
  * Tacconi et al. (2013) Tacconi, L. J., Neri, R., Genzel, R., et al. 2013, ApJ, 768, 74
  * Tadaki et al. (2019) Tadaki, K.-i., Kodama, T., Hayashi, M., et al. 2019, PASJ, 71, 40
  * Talia et al. (2020) Talia, M., Cimatti, A., Giulietti, M., et al. 2020, arXiv e-prints, arXiv:2011.03051
  * Tescari et al. (2014) Tescari, E., Katsianis, A., Wyithe, J. S. B., et al. 2014, MNRAS, 438, 3490
  * Thi et al. (2020) Thi, W. F., Hocuk, S., Kamp, I., et al. 2020, A&A, 634, A42
  * Thielemann et al. (2003) Thielemann, F.-K., Argast, D., Brachwitz, F., et al. 2003, Nuclear Physics A, 718, 139
  * Tielens & Hollenbach (1985) Tielens, A. G. G. M. & Hollenbach, D. 1985, ApJ, 291, 722
  * Tomassetti et al. (2015) Tomassetti, M., Porciani, C., Romano-Díaz, E., & Ludlow, A. D. 2015, MNRAS, 446, 3330
  * Tornatore et al. (2007) Tornatore, L., Borgani, S., Dolag, K., & Matteucci, F. 2007, MNRAS, 382, 1050
  * Tornatore et al. (2010) Tornatore, L., Borgani, S., Viel, M., & Springel, V. 2010, MNRAS, 402, 1911
  * Valentini et al. (2019a) Valentini, M., Borgani, S., Bressan, A., et al. 2019a, MNRAS, 485, 1384
  * Valentini et al. (2019b) Valentini, M., Chiappini, C., Bossini, D., et al. 2019b, A&A, 627, A173
  * van den Hoek & Groenewegen (1997) van den Hoek, L. B. & Groenewegen, M. A. T. 1997, A&AS, 123, 305
  * Vazdekis et al. (2010) Vazdekis, A., Sánchez-Blázquez, P., Falcón-Barroso, J., et al. 2010, MNRAS, 404, 1639
  * Villaescusa-Navarro et al. (2018) Villaescusa-Navarro, F., Genel, S., Castorina, E., et al. 2018, ApJ, 866, 135
  * Vogelsberger et al. (2019) Vogelsberger, M., McKinnon, R., O’Neil, S., et al. 2019, MNRAS, 487, 4870
  * Wakelam et al. (2017) Wakelam, V., Bron, E., Cazaux, S., et al. 2017, Molecular Astrophysics, 9, 1
  * Walter et al. (2020) Walter, F., Carilli, C., Neeleman, M., et al. 2020, ApJ, 902, 111
  * Weingartner & Draine (2001) Weingartner, J. C. & Draine, B. T. 2001, ApJS, 134, 263
  * Whalen et al. (2013) Whalen, D. J., Fryer, C. L., Holz, D. E., et al. 2013, ApJ, 762, L6
  * Wise et al. (2012a) Wise, J. H., Abel, T., Turk, M. J., Norman, M. L., & Smith, B. D. 2012a, MNRAS, 427, 311
  * Wise et al. (2012b) Wise, J. H., Turk, M. J., Norman, M. L., & Abel, T. 2012b, ApJ, 745, 50
  * Wolcott-Green et al. (2017) Wolcott-Green, J., Haiman, Z., & Bryan, G. L. 2017, MNRAS, 469, 3329
  * Wolfire et al. (2010) Wolfire, M. G., Hollenbach, D., & McKee, C. F. 2010, ApJ, 716, 1191
  * Wolz et al. (2021) Wolz, L., Pourtsidou, A., Masui, K. W., et al. 2021, arXiv e-prints, arXiv:2102.04946
  * Woosley et al. (2002) Woosley, S. E., Heger, A., & Weaver, T. A. 2002, Reviews of Modern Physics, 74, 1015
  * Woosley & Weaver (1995) Woosley, S. E. & Weaver, T. A. 1995, ApJS, 101, 181
  * Yoshida et al. (2003) Yoshida, N., Abel, T., Hernquist, L., & Sugiyama, N. 2003, ApJ, 592, 645
  * Zafar et al. (2013) Zafar, T., Popping, A., & Péroux, C. 2013, A&A, 556, A140
  * Zanella et al. (2018) Zanella, A., Daddi, E., Magdis, G., et al. 2018, MNRAS, 481, 1976
  * Zavala et al. (2021) Zavala, J. A., Casey, C. M., Manning, S. M., et al. 2021, arXiv e-prints, arXiv:2101.04734

##  Appendix A Data samples

Data samples for H2 mass density values, $\rho_{\rm H_{2}}$ (and corresponding
$\Omega_{\rm H_{2}}$), are summarised in Tab. 2, where redshift range, ${\rm
H_{2}}$ measurements and publication references are quoted. Many of them are
present in Péroux & Howk (2020).

Table 2:  Observational determinations of molecular (H2) gas considered in
this work.

$z_{-}$ | $z_{+}$ |  Log $\rho_{{\rm H_{2}},-}$ |  Log $\rho_{{\rm H_{2}},+}$ |  Log $\Omega_{{\rm H_{2}},-}$ |  Log $\Omega_{{\rm H_{2}},+}$ | Reference | Survey |
---|---|---|---|---|---|---|---|---
0.0 | 0.5 | 5.903 | 6.940 | -5.230 | -4.194 |  Hamanowicz et al. (2021)⋆ | ALMACAL-CO |
0.5 | 1.0 | 6.991 | 7.510 | -4.142 | -3.623 | ” | ” |
1.0 | 1.5 | 6.832 | 7.420 | -4.301 | -3.714 | ” | ” |
0.0 | 0.5 | 6.362 | 7.228 | -4.772 | -3.906 | ” | ” |
0.0 | 0.5 | 5.08 | 6.57 | -6.05 | -4.56 | Garratt et al. (2021) | UKIDSS-UDS |
0.5 | 1.0 | 7.10 | 7.33 | -4.03 | -3.80 | ” | ” |
1.0 | 1.5 | 7.28 | 7.44 | -3.85 | -3.69 | ” | ” |
1.5 | 2.0 | 7.47 | 7.60 | -3.66 | -3.53 | ” | ” |
2.0 | 2.5 | 7.21 | 7.40 | -3.92 | -3.73 | ” | ” |
- | 0 | 6.581 | 6.992 | -4.553 | -4.142 | Fletcher et al. (2021) | xCOLDGASS |
0.003 | 0.369 | 5.17 | 6.36 | -5.96 | -4.77 | Decarli et al. (2020) | ASPECS |
0.271 | 0.631 | 6.78 | 7.33 | -4.35 | -3.80 | ” | ” |
0.695 | 1.174 | 7.44 | 7.87 | -3.69 | -3.26 | ” | ” |
1.006 | 1.738 | 7.61 | 7.87 | -3.52 | -3.26 | ” | ” |
2.008 | 3.107 | 7.27 | 7.65 | -3.86 | -3.48 | ” | ” |
3.011 | 4.475 | 7.20 | 7.52 | -3.93 | -3.61 | ” | ” |
2.0 | 3.0 | 6.96 | 7.80 | -4.17 | -3.33 | Riechers et al. (2020a) | VLASPEC |
1.95 | 2.85 | 7.041 | 7.748 | -4.093 | -3.386 | Riechers et al. (2020b) | COLDz |
4.90 | 6.70 | 6.146 | 7.041 | -4.988 | -4.093 | ” | ” |
0.4799 | 1.4799 | 7.47 | 7.72 | -3.66 | -3.41 | Lenkić et al. (2020) | PHIBBS-2 |
1.0056 | 2.0056 | 7.56 | 7.78 | -3.57 | -3.35 | ” | ” |
2.0084 | 3.0084 | 7.42 | 7.67 | -3.71 | -3.46 | ” | ” |
3.1404 | 4.1404 | 6.91 | 7.34 | -4.22 | -3.79 | ” | ” |
4.2451 | 5.2451 | 6.58 | 7.22 | -4.55 | -3.91 | ” | ” |
0.003 | 0.369 | - | ¡8.26 | - | ¡-2.87 | Klitsch et al. (2019) | ALMACAL-abs |
0.2713 | 0.6306 | - | ¡8.32 | - | ¡-2.81 | ” | ” |
0.6950 | 1.1744 | - | ¡8.39 | - | ¡-2.74 | ” | ” |
1.006 | 1.738 | - | ¡8.21 | - | ¡-2.92 | ” | ” |
2.0 | 3.0 | 6.00 | - | -5.14 | - | Berta et al. (2013) | Herschel PEP |

From left to right, columns refer to lower and upper redshift range limits, H2
mass density $\rho_{\rm H_{2}}\rm[M_{\odot}/Mpc^{3}]$ with corresponding
$\Omega_{\rm H_{2}}$ range limits (1$\sigma$), publication reference from
which H2 data are taken or converted and observational survey.

⋆ Submitted. Private communication.

##  Appendix B HI and H2 self-shielding

We write the HI self-shielding, $f_{\rm shield,HI}$, as a function of H number
density, $n_{\rm H}$, as:

$f_{\rm shield,HI}=0.985\left[1+\left(\frac{n_{\rm
H}}{n_{0}}\right)^{1.9}\right]^{-2}+0.015\left(1+\frac{n_{\rm
H}}{n_{0}}\right)^{-0.72},$ (1)

where $n_{0}=0.006\,\rm cm^{-3}$ is a parameter. Asymptotically $f_{\rm
shield,HI}$ converges to 1 for $n_{\rm H}\ll n_{0}$ and to 0 for $n_{\rm H}\gg
n_{0}$. Besides eq. 1, we have checked that the alternative use of the tables
by Rahmati et al. (2013) does not change substantially our findings, although
truncation at $z=5$ causes unstable behaviours in abundance calculations
around that particular redshift.
For H2 shielding we refer to Draine & Bertoldi (1996) (their eq. 36 in most
cases and their eq. 37 in selected runs). Alternative formulations for H2
shielding, such as the revisitations by Sternberg et al. (2014) or Wolcott-
Green et al. (2017) of Draine & Bertoldi (1996)’s expressions, have a similar
impact on the $\Omega_{\rm H_{2}}$ mass density parameter. A detailed
comparison of different H2 shielding models is also done in Luo et al. (2020).
Maximum deviations from the reference run, reached locally at different
epochs, can be estimated by looking at the maximum H2 mass fractions recorded
at different $z$. These are more sensitive than $\Omega_{\rm H_{2}}$ and are
plotted in Fig. 7 for simulations including H-, H${}_{2}^{+}$ and three-body
formation channels, with H2 self-shielding parameterizations by Draine &
Bertoldi (1996) and the more recent Wolcott-Green et al. (2017) one. Although
we observe up to an order of magnitude differences at selected redshifts,
overall the maximum H2 mass fraction reached is not substantially affected by
the prescription used to model the H2 self-shielding.

Figure 7: Redshift evolution of the maximum H2 abundance reached in
simulations with different H2 self-shielding formulations: eq. 36 (violet dot-
dot-dot-dashed line; HM-HISSmed-3b in Tab. 1) and eq. 37 (green dotted line;
HM-HISSmed-3b-DB2) in Draine & Bertoldi (1996) and the formulation by Wolcott-
Green et al. (2017) (magenta dashed line; HM-HISSmed-3b-WG).

##  Appendix C H2 grain catalysis

The rate of H2 formation due to dust grain catalysis at gas temperature $T$
and grain temperature $T_{\rm gr}$ is:

$k_{\rm H_{2},gr}=3.5\times 10^{-17}\,S_{\rm
H}\sqrt{\frac{T}{100\leavevmode\nobreak\ {\rm K}}}\,{\rm cm^{3}/s},$ (2)

where the sticking coefficient $S_{\rm H}$ is given by Hollenbach & McKee
(1979); Cazaux & Tielens (2004):

$\frac{1}{S_{\rm H}}=1+0.4\sqrt{\frac{T+T_{\rm gr}}{100\leavevmode\nobreak\
{\rm K}}}+0.2\left(\frac{T}{100\leavevmode\nobreak\ {\rm
K}}\right)+0.08\left(\frac{T}{100\leavevmode\nobreak\ {\rm K}}\right)^{2}.$
(3)

In the absence of precise determinations at different epochs and $Z$, $k_{\rm
H_{2},gr}$ is linearly scaled by $Z$, via $k_{\rm H_{2},gr}\rightarrow k_{\rm
H_{2},gr}\times\,Z/\rm Z_{\odot}$.
The normalization of eq. 2 is similar to the one in Thi et al. (2020) and
Sternberg et al. (2014), while Omukai et al. (2005) divide eq. 2 by $1+{\rm
exp}{[750(1/75-1/T_{\rm gr})]}$, which suppresses catalysis by a factor $\geq
2$ at grain temperatures $T_{\rm gr}\geq 75\,\rm K$. The energy density
transfer rate, $\Lambda_{\rm gr}$, is computed accordingly. In the canonical
form, calibrated in the near Universe,

$\Lambda_{\rm gr}\simeq 3\times 10^{-33}\sqrt{T}\left(T-T_{\rm
gr}\right)\,n_{\rm H}^{2}\,\rm erg/s/cm^{3}.$ (4)

Also in this case there are alternative formulations, as e.g. in Omukai et al.
2005, where the expression is multiplied by the dimensionless factor
$1-0.8{\rm exp}{(\rm-75\,K/T)}$. Similarly to the H2 formation rate, the