Radiative cooling rates, ion fractions, molecule abundances and line emissivities including self-shielding and both local and metagalactic radiation fields

We use the spectral synthesis code Cloudy to tabulate the properties of gas for an extensive range in redshift (z=0 to 9), temperature (log T [K] = 1 to 9.5), metallicity (log Z/$\mathrm{Z}_{\odot}$ = -4 to +0.5, Z=0), and density (log $n_{\mathrm{H}}$ [$\mathrm{cm}^{-3}$] = -8 to +6). This therefore includes gas with properties characteristic of the interstellar, circumgalactic and intergalactic media. The gas is exposed to a redshift-dependent UV/X-ray background, while for the self-shielded lower-temperature gas (i.e. ISM gas) an interstellar radiation field and cosmic rays are added. The radiation field is attenuated by a density- and temperature-dependent column of gas and dust. Motivated by the observed star formation law, this gas column density also determines the intensity of the interstellar radiation field and the cosmic ray density. The ionization balance, molecule fractions, cooling rates, line emissivities, and equilibrium temperatures are calculated self-consistently. We include dust, cosmic rays, and the interstellar radiation field step-by-step to study their relative impact. These publicly available tables are ideal for hydrodynamical simulations. They can be used stand alone or coupled to a non-equilibrium network for a subset of elements. The release includes a C routine to read in and interpolate the tables, as well as an easy to use python graphical user interface to explore the tables.


INTRODUCTION
Radiative processes are a critical ingredient for all models that include baryons. Radiative losses are not only crucial for the formation of stars and therefore for the ignition of the cosmic baryon cycle, they are also heavily used in observational astronomy to classify the properties of gas within galaxies (interstellar medium, ISM), around galaxies (circumgalactic medium, CGM) and in between galaxies (intergalactic medium, IGM). The photons that are emitted at specific energies are important tracers for the composition, physical properties (e.g. density, temperature), and the movement of the gas (e.g. turbulence, inflow, outflows).
Determining the radiative cooling rate of a parcel of gas relies on a series of intertwined chemical reactions of species that can also interact with the ambient radiation field. A full implementation of these processes in a simulation would E-mail: ploeckinger@lorentz.leidenuniv.nl therefore require solving a large chemical network together with a coupled radiative transfer code for each resolution element at each timestep. While this can currently be done for a limited set of species in small-scale simulations of galaxy formation and evolution (e.g. isolated galaxies, or patches of the ISM), it is computationally too expensive to do in a large cosmological volume. A popular compromise to include radiative processes in simulations is to tabulate the cooling (and heating) rates. Therefore, the way these tabulated rates are calculated can be arbitrarily complicated without affecting the runtime of the simulation.
Traditionally, radiative cooling functions were calculated under the assumption of an optically thin gas in collisional ionization equilibrium (CIE, e.g. Cox & Tucker 1969;Dalgarno & McCray 1972;Raymond et al. 1976;Shull & van Steenberg 1982;Gaetz & Salpeter 1983;Boehringer & Hensler 1989;Sutherland & Dopita 1993). The ionization rates in CIE (and therefore the recombination rates and the radiative losses, Λ) are proportional to the number of colli-sions between particles (electrons and atoms). As hydrogen is the most abundant element in the Universe, the number density of particles scales with the hydrogen number density n H and the normalised cooling rate Λ/n 2 H [erg s −1 cm 3 ] is only a function of temperature (case I in Table 1). In this case, the contribution of individual elements scales linearly with their element abundance and the cooling rate for each chemical element can either be provided in terms of tables (e.g. Boehringer & Hensler 1989) or fitting functions (e.g. Dalgarno & McCray 1972).
When including photo-ionisation, the gas is over-ionised compared to the CIE case. This leads to different cooling rates caused by the different ion fractions of each individual element (e.g. Efstathiou 1992;Wiersma et al. 2009). As the photo-ionization rate scales with the number of particles to ionise (again approximated by n H ), while collisional processes scale with n 2 H , the normalised cooling rate Λ/n 2 H [erg s −1 cm 3 ] that combines both types of processes is a function of both density and temperature (case II in Table 1).
As the radiation field is not constant throughout the Universe, additional dependences are introduced. In galaxyscale or cosmological simulations, a redshift-dependent uniform UV background (UVB) is often used as radiation field, which adds a redshift dependence to the tabulated cooling rates (case IIa in Table 1, e.g. Wiersma et al. 2009, hereafter: W09). W09 showed that at the peak of the cooling function, between temperatures of 10 4 and 10 6 K, the cooling rate is reduced by up to an order of magnitude compared to the CIE case when photoionization from the UVB is included. Gnedin & Hollon (2012) tabulate the properties of photo-ionized gas for general radiation fields from stars or active galactic nuclei (AGN). They parametrize their incident spectrum with 4 independent values, one for the intensity and three that define its shape (case IIb in Table 1), resulting in an equal number of additional dimensions (7 in total, including density, temperature, and metallicity) for their tabulated cooling rates of optically thin gas in ionization equilibrium.
At temperatures between 10 4 and ≈ 5 × 10 6 K the cooling timescale can become shorter than the recombination timescale, leading to a "recombination lag". This rapidly cooling gas is therefore over-ionized compared to ionization equilibrium (e.g. Kafatos 1973). Gnat & Sternberg (2007) provide tables of radiative cooling rates including this nonequilibrium effect. Their calculations start with initially hot (T > 5 × 10 6 K) gas that cools either at constant density or constant pressure in the absence of a radiation field. The resulting cooling functions are smoother with a less pronounced peak due to hydrogen recombination compared to the equilibrium rates. As metal line cooling is important, their tabulated cooling rates also depend on the gas metallicity (case N-I, Table 1).
Similar to the approach from Gnat & Sternberg (2007), Oppenheimer & Schaye (2013) use their non-equilibrium reaction network to provide tabulated cooling and photoheating rates by following hot gas as it cools either at constant pressure or constant density. In contrast to Gnat & Sternberg (2007), Oppenheimer & Schaye (2013) include a redshift-dependent background radiation field and they tabulate their cooling rates for each ion of each considered element. As a hydrodynamic code typically does not trace the Table 1. In collisional ionization equilibrium (CIE), the cooling rate Λ scales with the hydrogen number density, n 2 H , and the term Λ/n 2 H therefore only depends on the temperature T . Including extra processes (rows) induces additional dependences on n H , redshift z, hydrogen column density N H , gas metallicity Z, as well as the intensity, I , and shape of the radiation field (RF). individual ion abundances, they also provide tabulated ion fractions of 11 elements, depending on the hydrogen number density, temperature, and metallicity of the gas, as well as the redshift for the background radiation field (case N-IIa in Table 1). Oppenheimer & Schaye (2013) conclude that non-equilibrium effects are reduced when the background radiation field is included. All above approaches focus on ionized gas with temperatures above 10 4 K, as is typical for the CGM and IGM. At these temperatures, molecules cannot form, dust grains are quickly destroyed by thermal sputtering (e.g. Tielens et al. 1994) and gas is optically thin to ionizing radiation. In the ISM, the incident radiation is however attenuated by gas and dust and can self-shield from photo-ionizing or photodissociating radiation. As the ion fractions depend on the photo-ionization rate, the cooling rate varies with the hydrogen column density N H (case III in Table 1). This can lead to structures where the illuminated side of a gas parcel is ionized, but deeper into the cloud, the gas can remain neutral and cold.
In neutral gas, the formation of molecules and the existence of dust grains further complicates the picture. Dust grains and molecules shield different parts of the intrinsic radiation field and the reaction rates of various molecules have a non-trivial dependence on the shape of the spectrum. In addition, an important formation channel of H 2 depends on dust grains as catalysts, which results in a non-trivial relation between the gas metallicity and the H 2 abundance.
The ISM gas resolution elements in hydrodynamic simulations can have a large variety of thermal and chemical histories: from gas that cools through thermal instabilities and is accreted from the galaxy halo to dense gas that is heated by nearby stellar feedback events. Tabulating nonequilibrium cooling rates and other chemical properties in a similar way (i.e. cooling from high temperatures) is therefore not a good approach. Chemical networks that solve a large number of differential equations to follow the thermo-chemical evolution of the gas particles are required to include non-equilibrium effects in the ISM.
grackle (Smith et al. 2017) is an open-source chemistry and cooling library that allows users to include hydrogen and helium non-equilibrium chemistry directly in the simulation. The primordial network includes the self-shielding approximation from Rahmati et al. (2013) (hereafter, R13). The cooling and heating rates from metals are calculated with cloudy under the assumption of ionization equilibrium, and tabulated for the UV background from Haardt & Madau (2012). In the 2017 grackle release, the metal cooling and heating rates assume that the gas is optically thin to the UV background and solar abundances. Combining the primordial network with the tabulated metal rates leads to inconsistencies in their electron fractions if self-shielding is included. The free electron fraction from hydrogen and helium used in the primordial network that includes selfshielding can be much lower than the free electron fraction from hydrogen and helium used in the optically thin cloudy tables. The metal cooling rates are therefore based on gas that is too highly ionized, leading to too high cooling rates for self-shielded gas. This has been discussed in Hu et al. (2017) and Emerick et al. (2019) produced metal cooling tables for grackle that include self-shielding of the UV background by hydrogen and helium, while shielding on metals and dust grains is neglected. Furthermore, as dust grains are not included in the cloudy calculations, grackle's shielding column lacks H 2 . The rates are calculated using a shielding length based on the local Jeans length (limited to a physical size of 100 pc) and for a single metallicity with solar relative abundances. Glover et al. (2010) present an ISM model for CO formation in turbulent molecular clouds within the ISM, including, in addition to hydrogen and helium, both carbon and oxygen chemistry. An extended chemical network with 272 chemical reactions that follow the evolution of 157 species (among these: 20 different molecules) is described in Richings et al. (2014a) for optically thin gas, and Richings et al. (2014b) for shielded gas. The effects of self-shielding compared to the assumption of optically thin gas are demonstrated in simulations of isolated galaxies by Richings & Schaye (2016). These large chemical networks that also calculate metals in non-equilibrium are computationally expensive, as the number of reactions increases rapidly with every included species.
The aim of this work is to provide tables of gas properties in ionization equilibrium that include radiation by a UV background and, optionally, from an interstellar radiation field (ISRF) for both unshielded (or optically thin to the radiation field) and self-shielded gas. The large number of dimensions would make the tables expensive to produce and would also lead to large memory requirements during any simulation run using the tables. More importantly, neither the shielding gas column density nor the local radiation field is generally known on the fly in large-scale simulations.
We solve this issue by assuming that the gas is selfgravitating and hence that its coherence scale is the local Jeans length (Schaye 2001a,b). Note that Rahmati et al. (2013) showed that this assumption reproduces the shielding lengths in their cosmological, radiative transfer simulations. By setting the shielding column density to half of the local Jeans column density we remove the N H dimension. In addition, we assume that the incident radiation field, which could be described by an arbitrary number of parameters, comes in two flavours: (1) the redshift-dependent model of the UVB by Faucher-Giguère (2020) (hereafter: FG20), modified to make the treatment of attenuation before H i and He ii reionization more self-consistent, and (2) optionally, a local interstellar radiation field (ISRF) with a constant spectral shape and a normalization that depends on the local Jeans column density of the gas as suggested by the observed Kennicutt-Schmidt (Kennicutt 1998) star formation law (case Vc,d in Table 1). For comparison, we also include the unshielded versions of the same tables (case Va,c in Table 1). We use the spectral synthesis code cloudy v17.01, last described in Ferland et al. (2017), to transmit the spectrum through the gas (for self-shielded gas), calculate the equilibrium ionization states of all atomic ions, and follow the formation and destruction of molecules to obtain the resulting cooling and heating rates for all included processes.
While we do not take non-equilibrium effects into account, the tables can easily be coupled to non-equilibrium networks that calculate the abundances of H and He, an approach that captures non-equilibrium effects on metals caused by the non-equilibrium free electron density (see Appendix A2).
The data products include a set of tables in hdf5 format of cooling and heating rates (also per element), ion and molecule fractions, and line emissivities, a routine that reads in and interpolates the cooling rates (written in C), as well as a python-based graphical user interface (gui) to visualise the table content. In equations that refer directly to individual hdf5 datasets, a line with the dataset name is added (e.g. dataset: ShieldingColumnDensityRef in Eq. 6). Figures that can be reproduced directly with the provided gui are labelled accordingly in the figure captions. The python and C routines and additional information on how to download the tables as well as examples can be found on dedicated web pages (http://radcool.strw.leidenuniv.nl/ or https://www.sylviaploeckinger.com/radcool).
The methods are described in detail in Sec. 2, including how the shielding column (Sec. 2.1), the radiation field (Sec. 2.2), and the dust content (Sec. 2.4) vary within the four-dimensional grid in redshift, temperature, metallicity, and density. We present the data products in Sec. 4, where the individual tables of the full set are presented and their information content is listed and explained. In Sec. 5 we highlight a few key results, such as the dependence of the thermal equilibrium temperatures and the most important phase transitions (ionized -neutral, neutral -molecular) on redshift, metallicity, and the radiation field. We summarise the findings in Sec. 7 and add sections on how to use (Appendix A) and reproduce (Appendix C) the provided tables.
The grid spacing in temperature, density, redshift, and metallicity as well as most of the tabulated properties are in log. Throughout the paper log refers to log 10 and we use 10 −50 as a floor value for properties that would be zero. All bins are generally equally spaced in log but the metallicity dimension has an additional bin: primordial abundances (Z = 0, referred to as log Z/Z = −50). Table 2. All gas properties are tabulated on an equally spaced grid in the dimensions: redshift z, gas temperature log T , metallicity log Z/Z , and gas density log n H . The grid points for each dimension range from their minimum (column 3) to their maximum (column 4) value, with the spacing ∆ between the grid points listed in column 5. For metallicity an additional field is included for primordial abundances (log Z/Z = −50, i.e. Z = 0). The resulting number of bins per dimension in column 6 leads to a total of 3,089,636 grid points for each

METHOD
For the calculation of cooling and heating rates as well as ionization stages and molecular fractions, we use version 17.01 of the photoionization code cloudy (Ferland et al. 1998(Ferland et al. , 2017. cloudy is a powerful and widely used open source code that can calculate the ionization, chemical, and thermal equilibrium state of gas exposed to an external radiation field. Our basic setup is a large grid of cloudy runs with dimensions z (redshift), T (temperature [K]), Z (metallicity [Z ]), and n H (total hydrogen number density [ cm −3 ]). The aim of this work is to provide tables that cover the full range of gas properties typically occurring in galaxy-scale or cosmological simulations (see Table 2). Throughout the paper we use the solar metallicity (Z = 0.0134) and solar abundance ratios from Asplund et al. (2009) (see Table 3).
For every grid point in the tables that include selfshielding, the chosen incident radiation field is set to propagate through a slab of plane-parallel gas with the given properties (T, Z/Z , n H ) until its shielding column density N sh is reached. cloudy bins the gas column into individual zones of adaptive widths and in each zone the chemical properties and ionization states are calculated and the resulting output spectrum is passed as input into the next zone. Within each gas column, the total hydrogen number density as well as the gas temperature are set to remain constant, while the ion fractions, electron densities, and molecule fractions depend on the depth into the slab of gas. For unshielded gas, the properties of the first zone are tabulated (1-zone model), while for self-shielded gas the same properties are tabulated for the last zone, where N H = N sh .
For ionized gas in ionization equilibrium, the contribution of individual elements to the total cooling and heating rates can be calculated by running the same 1-zone cloudy models once with the full abundance set and once excluding individual elements. The difference in the rates can then be attributed to the excluded element (this is done e.g. in W09). Here, this approach is not suitable, as a combination of elements is involved in the formation of molecules 1 and the shielding contributions from one species can affect the abundance of another species. We therefore run cloudy for each grid point only once and output the values for each co- oling and heating channel, as listed in Tables 9 and 10. For a full list of gas properties tabulated for each grid point, see Sec. 4.

Column density
We have already established that the gas properties of selfshielded gas depend on the depth (i.e. gas column density) into the irradiated slab of gas. While algorithms to calculate the optical depths for each resolution element in a simulation have been developed (e.g Wünsch et al. 2018 for the grid code FLASH, Fryxell et al. 2000), the shielding column or optical depth can be expensive to calculate and is not generally known on the fly. The shielding length is often approximated based on the gradient in the gas density ρ referred to as Sobolov approximation (used e.g. in Richings & Schaye 2016;Hopkins et al. 2018), or is assumed to be close to the local Jeans length for self-gravitating gas (e.g. in Schaye 2001b) where γ is the ratio of specific heats, X is the hydrogen mass fraction, k B is the Boltzmann constant, m H is the proton mass, µm H is the mean particle mass, G is the gravitational constant and n H and T are the total hydrogen density and temperature. The resulting Jeans column density is The cloudy calculation stops when a specified total hydrogen column density is reached. Using the values for γ and µ from the current cloudy zone leads to discontinuities in the shielding column densities at phase transitions as µ changes steeply (ionized -neutral: µ ≈ 0.5 → 1, atomic -molecular: µ ≈ 1 → 2). In addition, the last zone would preferentially be at the phase transition, as the shielding length suddenly decreases (N J ∝ µ −0.5 ). We therefore use the values of primordial neutral gas (from Planck Collaboration et al. 2016) for the shielding column density: X = 0.7563, γ = 5/3, and µ = 1.2328, which results in The Jeans column density is defined in Schaye (2001a,b), and R13 showed that this approximation for the shielding column reproduces the results of full radiative transfer calculations for H i in cosmological simulations that however did not include a cold ISM. We use the Jeans length approximation here, as it has the additional advantage that N J depends only on gas temperature and density, so no extra dimension ∇ρ is required in the tables.
We cover a large range of gas densities and temperatures (see Table 2) and the Jeans column density cannot be directly applied everywhere. High temperature gas (T 10 5 K) is mostly ionized and therefore optically thin, but as N J ∝ T 0.5 , its theoretical shielding length would continue to increase for a constant density. This leads to practical problems in cloudy, as the 1D shielding column starts to radiate thermal emission for high temperatures which results in an increasing radiation field with shielding column density. In addition, the shielding length can become very large for very low densities (n H < 10 −6 cm −3 ), leading to unrealistic coherence length scales corresponding to dynamical and sound crossing times that exceed the Hubble time. We therefore limit the length scale to l max = 100 kpc, the total hydrogen column densities to N max = 10 24 cm −2 and use an asymptotic function for a smooth transition between low ( 10 3 K, shielded) and high ( 10 5 K, optically thin) temperature gas.
The reference column density, N ref , is defined as dataset : ShieldingColumnRef where k = 5/ln(10) = 2.17 regulates the steepness of the transition between the asymptotes N ref reached at T √ T min T max and N min at T √ T min T max . We will use T min = 10 3 K and T max = 10 5 K. N ref is based on the length and column density limited Jeans column density N J : N min = l max n H,min = 3.08 × 10 15 cm −2 is the column density used for the minimum density in the tables n H,min = 10 −8 cm −2 . Eqs. 6 and 7 are illustrated in Fig. 1. N ref (dashed line) follows the Jeans length N J (dotted lines) for low temperatures but is limited by l > l max = 100 kpc (grey area) at low densities and by N max at high densities. Adding the asymptotic transition to N min for T √ T min T max leads to the final N ref (solid lines). A plot of log N ref as a function of density and temperature is presented in Fig. 2.
For unshielded tables, the cloudy calculation is stopped for each grid point after the first zone, where the maximum size of the first zone 2 is set to ∆r = 10 20 cm. The code stops adding zones when the total hydrogen column density exceeds the shielding column density log N sh = log(0.5 × N ref ) for shielded gas log(∆r × n H ) for optically thin gas (8)

dataset : ShieldingColumn
N sh is therefore the column density from the edge of the self-gravitating gas clump/disc to its centre/midplane, while N ref describes the total column column density through the gas clump. The shielding column density to the last zone of the cloudy column is tabulated as dataset ShieldingColumn. Depending on the thickness of the last zone, a small deviation from 0.5N ref is possible. For unshielded runs, Shiel-dingColumn is the column density of the first zone, while the values for N ref in dataset ShieldingColumnRef are the same for both unshielded and self-shielded runs. Figure 1. Temperature dependence of the discussed hydrogen column densities for selected gas densities n H (different colors). The Jeans column density N J (Eq. 5, dotted lines) increases with temperature as N J ∝ T 0.5 for a constant gas density. N ref (Eq. 7, dashed lines) follows the Jeans column density for low temperatures but is limited by a maximum column density N max = 10 24 cm −2 and a maximum length scale of l max = 100 kpc (the grey area indicates where l J > 100 kpc). The final reference column density N ref (Eq. 6, solid lines) includes the asymptotic transition in temperature from self-shielded gas to optically thin gas with a minimum column density N min . The values for N ref , which corresponds to twice the shielding column density, over the full temperaturedensity parameter space are shown in Fig. 2.

Radiation field
The incident spectrum in all tables includes a redshift dependent UVB based on FG20 (Sec. 2.2.1), as well as the cosmic microwave background (CMB) as a blackbody spectrum with a redshift dependent temperature of T CMB = T 0 (1 + z) with T 0 = 2.725 K. In some tables a local radiation field is added that models the contribution from nearby stars (Sec. 2.2.2) in the ISM.

UV background radiation field
The spectrum for the UVB from distant galaxies and active galactic nuclei (AGN) is based on FG20, which consists of three main contributions. First, the radiation field from stars (for z ≤ 8) is from the BPASS (Eldridge et al. 2017;Stanway & Eldridge 2018) spectrum for a stellar population with a constant star formation rate at an age of 300 Myr and a metallicity of 0.1Z . This spectrum is dust attenuated with a constant E(B − V) = 0.129. The redshift-dependent normalization is chosen to match the observed total UV emissivity from GALEX for z < 2 (Chiang et al. 2019) and the restframe UV galaxy luminosity function from Bouwens et al. (in prep.) for z > 2 at 8.3 eV (1500Å). Second, the radiation from AGN is modelled by combining spectral templates of obscured (75 per cent) and unobscured AGN (25 per cent) and the spectrum is normalized to the AGN ionizing emissivity from Shen et al. (2020) at 912Å. Third, recombination emission from the photons absorbed by the ISM of starforming galaxies as well as the recombination emission lines related to the "sawtooth" absorption by He ii Lyman series lines.
In addition to redshift-dependent spectra, FG20 also provide effective photo-ionization Γ x,eff and photo-heating rates q x,eff (with x = { H i, He i, He ii }) that are constructed to produce a prescribed volume-averaged ionized fraction (H ii and He iii) for each redshift, assuming photoionization equilibrium. These effective rates are lower than the rates calculated directly from the spectra before reionization in order to match the electron scattering optical depth of τ e = 0.054, the best fit value to the Planck 2018 CMB data (Planck Collaboration et al. 2018). The best model in FG20 uses reionization redshifts of z rei,HI = 7.8 and z rei,HeII = 3.5, while at z < z rei,HI − ∆z rei,HI = 7.2 and z < z rei,HeII − ∆z rei,HeII = 3, reionization is completed and the respective effective rates match the rates from the spectra.
For the UVB used in this work, we modify the spectra from FG20 for z > 3 so that they yield rates that match the effective photo-ionization and photo-heating rates before H i and He ii reionization. The details of this method as well as the resulting rates are shown in Appendix B. In short, we attenuate the UVB by H i gas for z > z rei,HI − ∆z rei,HI (before H i reionization is complete) and by He ii gas for z > z rei,HeII − ∆z rei,HeII (before He ii reionization is complete). The gas column density in both cases is chosen to match the respective effective photoionization rates from FG20.
The resulting spectrum (including the CMB) is shown in Fig. 3 for redshifts between 0 and 9. The CMB contribution dominates for energies 10 −3 Ryd while photons with higher energies originate from the UVB. For z = 9 the FG20 UVB is greatly reduced at energies −2 < log E [Ryd] < 0 as the contribution from star forming galaxies is only inclu- Figure 3. The incident background radiation field present in all tables for different redshifts consists of the CMB (lowest energy bump) and the modified FG20 UV / X-ray background (see Sec. 2.2.1 and Appendix B for details). The vertical dotted lines indicate the H i (left) and He ii (right) ionization energies. Before reionization is complete, the UVB is reduced for z > 3 (z > 7.2) at energies above the He ii (H i) ionization energy. See Fig. B1 for a comparison to the original FG20 UVB.
ded for z ≤ 8. The vertical dotted lines indicate the H i and He ii ionization energies to illustrate the onsets of the above mentioned attenuation.

Interstellar radiation field
For gas in the ISM, the diffuse ISRF from nearby stars can be stronger than the UVB. As discussed in the introduction, this type of radiation field can be described by its intensity as well as an undefined number of parameters that model its spectral shape.
Here, we fix the shape of the ISRF to that described in Black (1987) for the Milky Way Galaxy and normalize it based on the reference column density N ref , which only depends on the gas density and temperature (Eq. 6). This is motivated by the observed star formation law and it avoids an extension of the tables to additional dimensions describing the ISRF and allows their use in simulations where the strength of the ISRF is not directly traced.
In the case of a universal initial mass function (IMF), the rate at which OB stars are formed is proportional to the local SFR and the normalization of the Black (1987) radiation field is therefore assumed to scale with the SFR surface density Σ SFR as where R is a renormalization factor. The solar neighbourhood values J 0 = 4.4 × 10 5 photons cm −2 s −1 sr −1 eV −1 (or 4πνJ ν,0 = 1.4 × 10 −3 erg cm −2 s −1 ) defined at a wavelength of 1000Å (Black 1987) and Σ SFR,0 = 10 −3 M yr −1 kpc −2 (Bonatto & Bica 2011) are used to normalise the ISRF (see . Total incident radiation field at z = 0, including the CMB, the UVB and the local ISRF for different reference column densities, N ref . The shape of the ISRF is from Black (1987) and its intensity scales with N 1.4 ref as suggested by the Kennicutt-Schmidt law (see Eq. 13). The intensity of the total radiation field is defined at 1000Å (vertical dashed line), < 18, the total radiation field converges to the FG20 plus CMB radiation field (black dotted line, see Fig. 3 for other redshifts) for most energies. Table 4 for an overview of solar neighbourhood quantities used in this work).
We assume that the SFR surface density Σ SFR follows the observed Kennicutt-Schmidt relation (Kennicutt 1998): with A = 1.515 × 10 −4 M yr −1 kpc −2 and n = 1.4. The amplitude A has been decreased to account for a Chabrier (2003) IMF. These are the same values as used in the EAGLE simulations (Schaye et al. 2015) for which a similar scaling of the ISRF was used in Lagos et al. (2015) to compute molecular fractions. The gas surface density for a SFR of Σ SFR,0 is therefore As the gas surface density is expressed as a total hydrogen column density in the tables, Σ g,0 is re-written as using X = 0.7563 as for the column density in Sec. 2.1. The final scaling of the ISRF based on N ref is dataset : RadField Figure 5. Total radiation field intensity defined at 1000Å relative to the MW value of J 0 (see Table 4) at z = 0 (top panel) and z = 3 (bottom panel) for table UVB dust1 CR1 G1. The radiation field is dominated by the ISRF at high densities and temperatures 10 4 K (with log J/J 0 ∝ 1.4 × log Σ SFR ) and smoothly transitions to the redshift-dependent UVB for higher temperatures and towards low densities. For reference, the dashed line indicates where N ref = N H,0 (MW value, see Table 4). Figures made with provided gui.
For an unscaled Black (1987) radiation field (log R = 0), the molecular hydrogen fraction f H2 rises above 10 per cent only for neutral hydrogen column densities of N H i + 2N H2 10 21 cm −2 for the thermal equilibrium temperature, which is inconsistent with observations. In addition, low-metallicity gas does not significantly cool below 10 4 K at any redshift. Renormalising the ISRF and the CR rate by a factor of 1/10 (log R = −1) alleviates these tensions and the thermal equilibrium properties are in much better agreement with observations. This mismatch between the ISRF and the observed H 2 fractions has been found in previous work, where the H 2 formation rate on dust grains has to be boosted for a radiation field of 2.3 J 0 to reproduce the observations (e.g. Gnedin et al. 2009, Gnedin & Kravtsov 2011, Gnedin & Draine 2014. They relate this boost factor to the clumping factor C ρ and calibrate their model to this free parameter. This would account for unresolved higher density gas with higher molecule fractions. They find that clumping factors between C ρ ≈ 10 and 30 are necessary to reproduce the observed H i -H 2 in the MW. As cloudy solves a chemical network including several hundred molecules and atomic species, artificially increasing the formation rate of H 2 by a factor that is a free parameter is neither practical nor self-consistent. In this case, decreasing the radiation field is a more reasonable approach. Note also that the assumed scaling of the ISRF with the star formation surface density may overestimate the radiation field. A higher Σ SFR is linked to stronger local shielding of the stellar sources and may lead to a larger star formation scale height, both of which would reduce the local ISRF. Fig. 4 shows the fiducial input spectrum (log R = −1) at −4 and log E [Ryd] −2 and the intensity of the total radiation field in this energy range therefore becomes independent of the density and temperature. This transition depends on the intensity of the UVB and therefore on the redshift (see Fig. 5).

Before reionization
We apply the ISRF scaling J ∝ N 1.4 ref at all column densities, so also at column densities lower than expected for the ISM. After reionization this is not a problem, because the UVB dominates over the ISRF at low densities. However, before reionization the ionizing radiation of distant galaxies is absorbed close to the source, while local radiation fields, such as the ionizing radiation from nearby stars (Sec. 2.2.2) remain unaffected. As the local radiation field is not attenuated, its contribution can even become dominant at low column / volume densities (log n H [ cm −3 ] −5), which would not be part of the ISM and should therefore not be exposed to the ISRF. In order to avoid artificial gas heating, the ISRF normalization is drastically reduced compared to the reference value from log n H [ cm −3 ] = −2 to −6 with a similar asymptotic function as in Eq. 6: or : (10) and c J = e −8 regulate the steepness and position of the asymptotic transition to log(J/J 0 ) min = −20. Without this extra scaling the radiation field would Figure 6. The H 2 destruction rates by cosmic rays for H 2 dissociation ("CRPHOT,H2=>H,H", black lines) and ionization ("CRPHOT,H2=>e-,H2+", red lines). In cloudy v17.01 (solid lines), the ionization rate converges to the input secondary ionization rate of 4.6 × 10 −16 s −1 (red, dotted line) but the dissociation rate is a factor of ≈ 6 higher. In this work we rescale the rate for the reaction "CRPHOT,H2=>H,H"(dashed line) so that it converges to 5 × 10 −17 s −1 (black, dotted line) for high densities, the value expected from the UMIST database (see text for details).
not converge to the (modified) FG20 background radiation for very low densities. Note, that this additional scaling is only necessary for z > 7.5 as the photo-ionization rate of the UVB dominates that of the ISRF for lower redshifts at these densities.

Cosmic rays
Cosmic rays (CRs) are charged particles accelerated to high energies by diffuse shock acceleration in supernova remnants (Bell 2004). Their production rate is therefore expected to be proportional to the rate of SNe and furthermore, for a constant IMF, to the SFR. This is supported by observations of starburst galaxies, where the CR rate is several orders of magnitude higher than locally in the MW Galaxy (e.g. Suchkov et al. 1993;Veritas Collaboration et al. 2009).
Analogously to the scaling of the ISRF intensity (Sec. 2.2.2), the H i CR ionization rate ζ CR is set to where ζ CR,0 = 2 × 10 −16 s −1 is the mean CR ionization rate of neutral hydrogen in local, diffuse clouds (Indriolo et al. 2007) and log R is the same renormalization constant as in Eq. 13 (we use log R = −1 for the fiducial model).
We use the large cloudy model for the H 2 molecule (described in Shaw et al. 2005) which includes several thousand levels for the H 2 molecule. The H 2 ionization rate is therefore not constant but depends for example on the gas density (see the red solid line in Fig. 6). For high densities the H 2 destruction rate by cosmic rays due to H 2 ioni- zation (labelled "CRPHOT,H2=>e-,H2+" in cloudy) converges to the input value of 4.6 × 10 −16 s −1 (red dotted horizontal line), following the ratio between H i and H 2 ionization rates from Glassgold & Langer (1974). According to the UMIST database 3 (McElroy et al. 2013), the ratio between the rates for H 2 destruction by ionization and dissociation (labelled "CRPHOT,H2=>H,H" in cloudy) is 0.108, and the dissociation rate is therefore expected to converge towards 5×10 −17 s −1 (black, dotted horizontal line). In cloudy 17.01, the dissociation rate is a factor of 6 higher than the ionization rate (black solid line). We therefore rescale this rate to converge towards the expected value of 5 × 10 −17 s −1 (black dashed line) for high densities. This is more consistent with other chemical network codes, such as chimes (Richings et al. 2014a,b) that use the UMIST values and leads to higher H 2 fractions, in better agreement with observations (see Sec. 5.1).
The high H 2 dissociation rate in cloudy is caused by the assumed mean kinetic energy of the secondary electrons of 20 eV. For this energy the cross section for H 2 dissociation is at a maximum (Dalgarno et al. 1999). Future versions of cloudy will use a mean kinetic energy of 36 eV, more representative for the broad range of energies of secondary electrons, and this rescaling will not be necessary anymore (see Shaw et al. 2020 for details).

Metallicity and dust
One of the table dimensions is the gas metallicity Z in units of solar metallicity Z , which ranges from log Z/Z = −4 to 0.5 in steps of 0.5 dex (see Table 2). The abundances of all elements other than hydrogen and helium are scales with metallicity relative to their solar values (Table 3). For primordial abundances, an additional metallicity bin is added (log Z/Z = −50), using the primordial helium abundance (n He /n H ) prim = 0.08246 (Planck Collaboration et al. 2016). To obtain a smooth transition between zero and solar me- tallicities, the helium abundance is linearly interpolated in between with (n He /n H ) = 0.0851 (Asplund et al. 2009). The tables contain various datasets with information about the used abundance ratios (see Sec. 4.1.2). The primordial abundance set in cloudy does not only include hydrogen and helium, but also traces of lithium and beryllium. Their minor cooling contribution explains why the metal cooling rate can be non-zero for primordial abundances.

Dust
We include the "Orion" grain set from cloudy. For solar metallicity, this grain set has a default dust-to-gas mass ratio of (D/G) 0 = 5.6 × 10 −3 and the relation between the extinction A V and the total hydrogen column density N H is (A V /N H ) = 5.54 × 10 −22 mag cm 2 . The grain distribution with sizes between 300 and 2500Å consists of both graphites and silicates and matches the extinction observed in Orion. In addition to dust grains, also polycyclic aromatic hydrocarbons (PAHs, with a power-law size distribution from Abel et al. 2008) are included as they contribute to photoelectric heating, collisional processes, and stochastic heating. As PAHs are destroyed in ionised gas by ionising photons and get depleted onto larger grains in molecular gas, their abundance is set to scale with the atomic neutral hydrogen fraction n H i /n H in cloudy.
A constant dust-to-metal ratio for ISM gas with N ref > N H,0 is assumed. As for this solar neighbourhood dust-tometal ratio, some elements (e.g. Mg, Si, Ca, Fe) are already almost completely depleted on dust grains, higher dust-tometal ratios would not be consistent. For lower reference column densities, a similar scaling as for the ISRF and the CR rate is used. The dust-to-gas mass ratio D/G is then For N ref = N min grain physics is completely disabled, as even extremely low grain abundances can cause cloudy to become unstable for high temperatures (i.e. T > 10 6 K).
Here, dust grains would be destroyed by thermal sputtering on very short timescales (e.g. Tielens et al. 1994), but this process is not included in cloudy. Fig. 7 summarises the dust-to-gas ratio dependence on density and temperature. The option no qheat is a standard cloudy command to disable quantum heating. This was necessary for code stability.

Metal depletion
The abundances of the Sun, which define the solar abundances, differ from those of the local ISM. For some elements a large fraction of atoms are depleted onto dust grains and their abundance ratios measured in the gas phase are correspondingly reduced. If the gas phase abundance were set to solar and dust grains were added to the mix without depleting the gas phase, the total (dust + gas) abundances would become super-solar.
For each element i, the number fraction of atoms that are depleted on dust grains is f dust,i , while f gas,i = 1 − f dust,i is the number fraction of atoms of element i in the gas phase. The reference values for solar neighbourhood conditions log f dust,i are taken from dataset : Depletion and Fig. 8 illustrates this scaling.

HD cooling
In the cloudy version used (17.01), deuterium chemistry is disabled by default, contrary to the information in the user manual 5 . Tests with a coarser grid where the deuterium chemistry was specifically added through the command line input, revealed that cloudy aborts for an increased number of grid points due to non-convergence issues in the chemistry solver. We therefore do not use the HD cooling from cloudy but add it analytically, following the cloudy prescription as closely as possible. As in cloudy v17, the cooling rate per HD molecule, W HD [erg s −1 ], is taken from the fitting function provided by Flower et al. (2000). W HD depends on the gas density and temperature and was calculated by Flower et al. (2000) for temperatures between 30 and 3000 K and densities n H ≥ 1 cm −3 . The cooling rate from HD per unit volume is where n HD [ cm −3 ] is the number density of HD molecules. As done in earlier cloudy versions, the HD abundance is assumed to follow the abundance of H 2 with n HD /n H2 = D/H, with the primordial deuterium abundance of D/H = 1.65 × 10 −5 (Pettini & Bowen 2001). The cooling rate Λ cool,HD is extrapolated to lower densities to avoid a sharp transition at n H = 1 cm −3 . As both n HD and W HD are lower for lower cloudyastrophysics.groups.io/g/Main/message/3934 Table 5. Overview of the models and the included processes. For the column "dust", an entry "yes" means that dust grains with a dust-to-gas ratio according to Eq. 18 are included and that individual elements are depleted onto dust grains as described by Eq. 19. Cosmic rays with a rate following Eq. 16 are included for runs labelled "yes" in the column "CR". In addition to the radiation from the CMB and the modified UVB from FG20 (all tables), an interstellar radiation field (Eq. 13) is added to the runs with "yes" in the column "ISRF". The shielding column density of Eq. 8 is used for self-shielded gas ("yes" in column "shielding"). The fiducial model "UVB dust1 CR1 G1 shield1" is highlighted in italic. The symbol indicates that the corresponding rate/intensity is uncalibrated (log R = 0 in Eqs. 13 and 16). densities, Λ cool,HD is negligible for these densities. For temperatures > 10 4 K, Λ cool,HD is set to zero, as W HD increases super-linearly with temperature, which can lead to nonnegligible HD cooling rates, even if the HD abundance is very low. At low redshift the HD cooling rate is only important for very low metallicities and very high densities (it contributes between 10 and 50 per cent of the total cooling rate for log Z/Z < −3 and log n H [ cm −3 ] > 3). For redshifts z > 8 HD cooling dominates for temperatures below the CMB temperature (T CMB (z = 9) = 27.3 K). This is an artefact from the analytic fitting function that does not include a redshift dependence. It only affects the lowest temperature bins at the highest redshifts in the table, but the HD cooling rate should be used with care in this region of parameter space.

MODELS
The fiducial model that we recommend for use in simulations is "UVB dust1 CR1 G1 shield1" as it includes all processes discussed in Sec. 2 (a redshift-dependent metagalactic radiation field, Sec. 2.2, as well as a density-and temperaturedependent shielding column density, Sec. 2.1, interstellar radiation field, Sec. 2.2, cosmic ray rate, Sec. 2.3, and dust abundance, Sec. 2.4). In this model, both the CR rate as well the ISRF intensity scale with the pressure as expected for a self-gravitating disk following the Kennicutt-Schmidt star formation law, with a normalization that is reduced by 1 dex relative to the Black (1987) value for the Milky Way (log R = −1 in Eqs. 13 and 16) in order to match the observed transition from atomic to molecular hydrogen (as will be discussed in Sec. 5.1). The unshielded counterpart of this model is "UVB dust1 CR1 G1 shield0" where all other cloudy inputs, beside the shielding column density, are identical to model "UVB dust1 CR1 G1 shield1".
For comparison of the fiducial mo-del to models without CRs and ISRF, we include "UVB dust1 CR0 G0 shield0" and "UVB dust1 CR0 G0 shield1" the unshielded ("shield0") and self-shielded ("shield1") models, where the radiation field only depends on redshift and is the sum of the contributions from the CMB and the modified FG20 UVB (see Fig. 3). The unshielded table "UVB dust1 CR0 G0 shield0" is conceptually comparable to the approach in W09, and can serve as an update to the W09 tables. The main differences are that we use a more recent cloudy version (here: v17.01, W09: v07.02), a different UV background (here: FG20, W09: Haardt & Madau 2001), and include dust grains for log T [K] 5. By construction, all other tables only differ from "UVB dust1 CR0 G0 shield0" for log T [ K] 5. At higher temperatures (and also very low densities) all tables converge to this one, as all additionally included processes are relevant only for neutral and molecular gas. Table 5 presents an overview of the models for which hdf5 files are provided. In addition to the above mentioned tables whose results will be further discussed in Sec. 5, additional tables are listed separately. We do not discuss them in detail here, as we do not recommend them for use in simulations, but they can be used to explore the impact of individual processes (e.g. a model with CRs but without ISRF: "UVB dust1 CR1 G0 shield1" or a model without dust: "UVB dust0 CR0 G0 shield0"). The model with the full ISRF intensity and CR rate (log R = 0 in Eqs. 13 and 16, "UVB dust1 CR2 G2 shield1") is included as it motivates the renormalization of these two quantities.
For self-shielded models ("shield1") without CRs the individual cloudy runs can easily become unstable and crash as CRs are the only source of ionization in highly molecular regions. This is known behaviour, it is mentioned clearly in the cloudy user manual that: "the chemistry network will probably collapse if the gas becomes molecular but cosmic rays are not present". Grid points for which cloudy crashes (7,223 out of 3,089,636 or 0.2 per cent for model "UVB dust1 CR0 G0 shield1" compared to 315 or 0.01 per cent for model "UVB dust1 CR1 G0 shield1") can be identified with the hdf5 dataset GridFails (see Sec. 4.1) while for all other datasets the values for the crashed runs are interpolated in 2D (density, temperature) between their neighbouring grid points with successful runs.

DATA PRODUCTS
The data products released with this publication are hdf5 files for the different models listed in Table 5. The main datafile of each model is available as "<modelname>.hdf5" (e.g. "UVB dust1 CR1 G1 shield1.hdf5") and contains both input datasets (e.g. CR rate, radiation field intensity) as well as cloudy results (e.g. cooling and heating rates, ion fractions, H 2 and CO abundances). The content of these files is described in Sec. 4.1 and listed in Table 6.

Main hdf5 files: rates and fractions
The results for each model listed in Table 5 are stored in an hdf5 file that contains both the data as well as documentation. An overview of all datasets is given in Table 6 and additional information can be found in the attributes Dimension, Info and Unit of each dataset. We discuss a selection of entries that require more explanation below.

General structure
The hdf5 group /TableBins/ contains the values for each dimension already mentioned in Table 2 and the number of grid points in each dimension is stored in the hdf5 group /NumberOfBins/. A cloudy run was started for each redshift, temperature, metallicity, and density, which results in N z × N T × N Z × N n H = 46 × 86 × 11 × 71 = 3, 089, 636 individual cloudy runs per table.
Many hydrodynamic solvers do not use the gas temperature, T, but the internal energy per unit mass, U, as their main hydro variable. Hence, to use cooling rates tabulated as a function of temperature, the internal energy first has to be converted into the gas temperature. The relation between U and T is non-linear at phase transitions (mainly H 2 -H i and H i -H ii), where both the mean particle mass µ and the ratio of specific heats γ change rapidly. To improve the usability of the code, we provide all full data hypercuboids (Sec. 4.1.3) both as a function of temperature (hdf5 group: /Tdep/) and internal energy (hdf5 group: /Udep/). The dimension of the full grid in /Udep/ is N z × N U × N Z × N n H (with the internal energy replacing the temperature dimension) and the internal energy bins can be found in dataset: InternalEnergyBins. The hdf5 group /ThermEq/ does not include a temperature or internal energy dimension (therefore the dimensions are: N z × N Z × N n H ), as it contains the gas properties at their thermal equilibrium temperature (dataset: ThermEq/Temperature) defined at the last zone of the cloudy calculation. Note that by definition the net cooling rate (cooling -heating) is zero for the thermal equilibrium temperature.

Miscellaneous information
We select the 11 elements that contribute most to the thermal state of the gas (i.e. the dominant cooling contributions) and tabulate their properties individually and in more detail (listed in Table 3 and in dataset ElementNames). The remaining elements are combined into one entry where this is relevant (e.g. the cooling channel "OtherA", see Table 9). Several datasets contain information about these elements (i.e. their names, masses and number of ions) and their abundances.
As mentioned in Sec. 2.4 (Eq. 17), the helium abundance varies with metallicity to yield a smooth transition between the solar and primordial values. The values used for n He /n H at each metallicity are stored in dataset Abundan-ceHe. The extra scaling of helium changes the helium mass Y which results in a slightly different metal mass fraction Z = Z/(X + Y + Z) compared to the values from the dataset  Table 9 (Table 10) from model UVB dust1 CR1 G1 shield1 (see Table 5) for a gas density of n H = 100 cm −3 , solar metallicity and at redshift 0. The components are split into primordial (left panel), atomic metal (middle panel) and the remaining metal processes (right panel). The black dotted line in each panel shows the total rate of the contribution of H and He (left panel) or from metals (middle and right panel) and the solid grey line in each panel indicates the total rate. The individual cooling (heating) processes are explained in Table 9 (Table 10). Figure made with provided gui. Table 6. Overview of all datasets stored in each hdf5 table. A dataset is labelled with the prefix (in) if it contains cloudy input, and with (out) if it contains results of the cloudy calculations. Note that cloudy assumes chemical and ionization equilibrium. Column 2 shows the size and dimensions of each dataset: z (redshift), T (gas temperature), Z (gas metallicity), n H (gas density), U (gas internal energy), E (individually tabulated elements), C (cooling channels), and H (heating channels). In general, all resulting properties refer to those in the last zone in the cloudy calculation. For hydrogen and CO, average fractions over the full column density are also stored. The suffix "Col" refers to the column density fraction and the suffix "Vol" refers to the volume density fraction at the last zone.  Table 9) (out) Heating z, , Z, n H , H Heating rate log Λ heat /n 2 H for different heating channels (see Table 10) (out) AVextend z, , Z, n H Extended-source extinction, total visual extinction in mag at 5500Å (out) AVpoint z, , Z, n H Point-source extinction, total visual extinction in mag at 5500Å (out) COFractionCol z, , Z, n H CO column density fraction log N CO /N H (out) COFractionVol z, , Z, n H CO volume density fraction log n CO /n H (out) ColumnDensitiesC z, , Z, n H , 3 Selected carbon column densities: log N C i , log N C ii , log N CO (out) ColumnDensitiesH z, , Z, n H , 3 Selected hydrogen column densities: log N H i , log N H ii , log N H2 (out) ElectronFractions z, , Z, n H , E + 3 Free electron fraction log n e /n H from each individual element + prim, metal, total (out) GammaHeat z, , Z, n H Ratio of specific heats γ, varies between 7/5 and 5/3 As HydrogenFractionsVol but also: log 2n H2 + /n H , log 3n H3 + /n H , log n H − /n H (out) IonFractions/00hydrogen z, , Z, n H , 2 Fractions of hydrogen atoms in each atomic ionization state (out) IonFractions/01helium z, , Z, n H , 3 As IonFractionsVol/00hydrogen, but for helium (out) IonFractions/02carbon z, , Z, n H , 7 As IonFractionsVol/00hydrogen, but for carbon (out) IonFractions/03nitrogen z, , Z, n H , 8 As IonFractionsVol/00hydrogen, but for nitrogen (out) IonFractions/04oxygen z, , Z, n H , 9 As IonFractionsVol/00hydrogen, but for oxygen (out) IonFractions/05neon z, , Z, n H , 11 As IonFractionsVol/00hydrogen, but for neon (out) IonFractions/06magnesium z, , Z, n H , 13 As IonFractionsVol/00hydrogen, but for magnesium (out) IonFractions/07silicon z, , Z, n H , 15 As IonFractionsVol/00hydrogen, but for silicon (out) IonFractions/08sulphur z, , Z, n H , 17 As IonFractionsVol/00hydrogen, but for sulphur (out) IonFractions/09calcium z, , Z, n H , 21 As IonFractionsVol/00hydrogen, but for calcium (out) IonFractions/10iron z, , Z, n H , 27 As IonFractionsVol/00hydrogen, but for iron (out) MeanParticleMass z, , Z, n H Mean particle mass µ in units of atomic mass unit u (out) Tdep/U from T z, T, Z, n H Internal energy log U[erg g −1 ] for group /Tdep/ (out) ThermEq/U from T z, Z, n H Thermal equilibrium internal energy log U[erg g −1 ] for group /Thermeq/ (out) Udep/T from U z, U, Z, n H Temperature log T [ K] for group /Udep/ (out) ThermEq/Temperature z, Z, n H Thermal equilibrium temperature MetallicityBins, which were used to scale the metal abundances independently. The difference is very small (< 0.03 dex), but for reference the exact metallicity values are stored in TotalExactMetallicity. This array can be used instead of MetallicityBins, but has the practical disadvantage that it is not exactly uniformly spaced in log Z. The dataset ElementMasses m i is added to allow the conversion between element abundances (n i /n H ) and mass fractions, The last entry of ElementMasses is the average atomic mass of the remaining elements for solar abundances from Asplund et al. (2009).
All metallicities are tabulated relative to the solar metallicity Z and can be converted to absolute metallicities by multiplying them with the entry from the dataset SolarMetallicity. Similarly, the radiation field dataset RadField is normalised to the local value J 0 , which is stored in the dataset JMW.
The datasets IdentifierCooling and IdentifierHeating contain the labels listed in columns 2 of Tables 9 and 10. They correspond to the final dimension of the datasets Cooling and Heating (see Sec. 4.1.3).

Full data hypercuboid
There are three big groups in each hdf5 file that contain the same datasets but are tabulated in slightly different ways. In group /Tdep/ each grid point corresponds directly to a cloudy run with redshift, temperature, metallicity and density as inputs that remain constant throughout each 1D cloudy simulation. As already mentioned in Sec. 4.1.1, some hydrodynamic codes use the internal energy as the thermal variable instead of the gas temperature. For this case, using the datasets from group /Udep/ saves having to convert from T to U, as the datasets now use the internal energy instead of the temperature dimension.
To create the datasets in group /Udep/, we first calculated the internal energy with dataset : Tdep/U_from_T with the Boltzmann constant k B = 1.3806 × 10 −16 erg K −1 and the atomic mass unit m u = 1.6605 × 10 −24 g. Both the ratio of specific heats γ and the mean particle mass µ vary within the table range. For µ we use the values from the dataset MeanParticleMass, which is a direct cloudy output. We approximate γ by assuming it is dominated by the contributions from hydrogen and helium as γ = 5 (n H i + n H ii + n He + n e ) + 7n H2 3 (n H i + n H ii + n He + n e ) + 5n H2 (22) dataset : GammaHeat assuring a smooth transition between γ = 5/3 for monoatomic gas and γ = 7/5 for molecular hydrogen (as done in the krome package, Grassi et al. 2014). In a second step, the properties of all datasets are interpolated onto fine, uniformly spaced bins in internal energy (dataset: Table-Bins/InternalEnergyBins).
The large phase-space covered by the tables is not uniformly populated by the gas resolution elements in a simulation. Typically gas will spend a disproportionally large fraction of its time close to thermal equilibrium. For all datasets, the group /ThermEq/ contains the properties at the thermal equilibrium temperature at the last cloudy zone (dataset: /ThermEq/Temperature). Occasionally, there are multiple equilibrium temperatures for a given density (at constant z and Z). This occurs for example if an individual cooling contribution (e.g. molecular hydrogen) peaks around a given temperature while the heating rate stays roughly constant, which results in two thermal equilibrium tracks (at temperatures below and above the cooling peak). In this case (usually only for a narrow range in densities) we select the track that changes the equilibrium temperature the least, if the gas is moving to higher densities. The regions with multiple thermal equilibrium temperatures can be easily identified in a temperature -density plot of the net cooling rate (e.g with the provided gui). The minimum temperature in all tables is 10 K. If the cooling rate exceeds the heating rate at 10 K, we set T eq = 10 K.
Input datasets: Some of the datasets in Table 6 have the prefix (in) . These are properties that are not calculated by cloudy but are used as inputs. As they vary over different dimensions, we output their values during runtime as an independent check to verify that the scaling is implemented correctly.
Cooling and heating: cloudy follows a large number of different cooling and heating channels. We combined them in the 20 cooling and 22 heating groups listed in Tables 9 and 10. In addition, the last two entries in this dimension are the total rates from processes including hydrogen and helium (labelled as: TotalPrim) and from processes including metals (TotalMetal). This allows one to calculate the rates faster for solar relative abundances by adding only the last two entries (for more details on how to use the tables, see Appendix A).
The heating and cooling labels from cloudy are listed in column 4 of Tables 9 and 10 and we refer to the cloudy documentation for more details on the individual processes. An example of the cooling and heating rates can be found in Fig. 9.
Vol and Col: A cloudy calculation that includes selfshielding splits the gas shielding column into a varying number of zones with adaptive sizes. Each output property (labelled with (out) in Table 6) is defined at each zone throughout the gas column. In this work, every tabulated property refers to its value calculated at the last zone of the cloudy simulation. For self-shielded gas the last zone is the zone at the shielding column density and for unshielded gas the last zone is the only zone in the calculation.
If a phase transition occurs just before the last zone, the species fraction in the last zone can be very different from the species fraction integrated over the full shielding column. For the hydrogen and CO fractions we therefore store both values: n x /n H (suffix: "Vol") at the last zone as well as N x /N H (suffix: "Col") for the full shielding column.
In simulations, the shielding length l sh = N sh /n H can be compared to the length scale of the resolution element l sim to Figure 10. An example of the dataset ElectronFractions for table UVB dust1 CR1 G1 shield1 for solar metallicity gas at redshift 0. The colour coding shows the free electron fraction n e /n H in the last zone of the shielding column split up into contributions from individual elements, as well as the total contributions from hydrogen and helium (including electrons from hydrogen molecules, TotalPrim), from metals (TotalMetal) and the total electron fraction (including electrons from all molecules, Total). decide whether n x /n H (for l sim > l sh ) or N x /N H (for l sim < l sh ) is a better approximation. n x /n H will only differ significantly from N x /N H in small parts of the full table parameter space. For example for UVB dust1 CR1 G1 shield1 n H2 /n H (n H i /n H ) differs by more than a factor of 2 from N H2 /N H (N H i /N H ) in less than 10 (1) per cent of the 3,089,636 entries.
Free electron fractions: This dataset contains the number of free electrons per hydrogen atom n e,X /n H for each individually traced element X. The number of electrons per gas phase atom of element X is calculated from the fractions of atoms in ion stage i with n e,X,gas (e.g. for He: n e,He,gas /n He,gas = n He ii /n He + 2n He iii /n He ). Accounting for the total abundance of element X and depletion onto dust (for tables that include grains), the number of electrons per hydrogen atom from element X is n e,X n H = n e,X,gas n X,gas n X,gas where n X,gas is the gas phase number density and n X the total (gas + dust) number density of element X. The contributions from ionized molecules are not included here, but for metals these are negligible. The last three entries in the element dimension are the total primordial electron fractions (including species n H2 + , n H3 + and n H − ), the total fractions of free electrons from metals n e,met /n H without molecules, and n e,tot /n H from all atoms and molecules included in cloudy. Fig. 10 illustrates the individual entries for solar metallicity gas at redshift 0 (using table UVB dust1 CR1 G1 shield1). Splitting the electron fractions into individual elements allows users to combine these tables with a chemical network that follows the non-equilibrium evolution of a subset of elements (e.g. H and He). The free electrons from metals in self-shielded gas can be added to the network individually (see Appendix A2). Table 3 is included in every hdf5 file. The size of the last dimension in each ion fraction dataset (see Table 6) is the number of ions, given in the additional dataset NumberOfIons. For each element the ion fractions are defined as the number of atoms in each ionization state relative to the total number of atoms of this element in the gas phase (e.g. n H i /n H ). The sum of all ion fractions for an element does not add up to one if a nonnegligible fraction of atoms are in molecules.

Line emissivity hdf5 files:
For tables UVB dust1 CR1 G1 shield0 and UVB dust1 CR1 G1 shield1 we provide additional information on the emissivity of selected lines (see Table 7 and dataset: IdentifierLines) in a separate hdf5 file ("<modelname> lines.hdf5"). The general file structure is the same as in the hdf5 files that contain the rates and species fractions (Sec. 4.1) as the emissivities are calculated for the same grid in redshift, metallicity, density and temperature.
For each line the hdf5 file contains the line emissivity in units of erg cm −3 s −1 that is leaving the illuminated side of the cloud (using the cloudy keyword "emergent"). The luminosity of a gas particle (or gas cell) in a simulation can be calculated with following e.g. Bertone et al. (2010), where m gas is the mass of the resolution element, ρ gas is the gas density and is either the emissivity at the last zone ( Vol , dataset: Emissi-vitiesVol) or the average emissivity of the shielding column ( Col , dataset: EmissivitiesCol"). The emergent emissivities of the last zone ( Vol ) are calculated by cloudy with the command "save lines emissivity" and the column averaged emissivity is calculated by dividing the emergent line intensity of the full shielding column ("save line list absolute") by the length of the shielding column. Analogous to the "Vol" and "Col" fractions from the main hdf5 file (discussed in Sec. 4.1.3), comparing the shielding length to the length scale of the resolution element can help decide which emissivity to choose.

EXAMPLE RESULTS: THERMAL EQUILIBRIUM
The full tables can be explored with the help of the provided gui. Here, we focus on the properties of gas close to the thermal equilibrium temperature at the last zone of the gas column.
An example for a density-temperature slice through the dataset hypercuboid for solar metallicity and z = 0 is shown in Fig. 11. The plotted properties are the cooling rate (Λ cool , left panel), heating rate (Λ heat , middle panel), as well as the ratio between the two (Λ cool /Λ heat , right panel). The thermal equilibrium temperature, T eq , is defined as the temperature Figure 12. Thermal equilibrium temperatures for two different radiation fields, only CMB and UVB (UVB dust1 CR0 G0, left panel) and including CRs and ISRF (UVB dust1 CR1 G1, right panel) at redshift 0. The colours in each panel refer to different gas metallicities and the line styles indicate if the thermal equilibrium temperature is for unshielded (dashed lines, " shield0") or selfshielded (solid lines, " shield1") gas. The shaded area indicates ISM pressures of 3 < log P/k B [ K cm −3 ] < 4 as observed in the solar neighbourhood by Jenkins & Tripp (2011). where Λ cool /Λ heat = 1 6 (dashed, white line). This temperature is stored in dataset /ThermEq/Temperature and all properties in the group /ThermEq/ correspond to the relevant (for each z and metallicity bin) thermal equilibrium temperatures.
T eq typically serves as a lower temperature limit for ISM gas. Gas can have higher temperatures due to hydrodynamic processes (e.g. shocks) but temperatures below T eq are rare as the heating rate typically dominates adiabatic losses at ISM densities. For densities around and below the cosmic mean 7 T eq is not a good proxy for the actual (quasi-)equilibrium temperature because here adiabatic cooling due to the Hubble expansion typically dominates over radiative cooling. More generally, when the radiative time scales exceed the Hubble time, the gas will not be able to achieve thermal equilibrium. Fig. 12 displays an overview of the thermal equilibrium temperatures at z = 0 for a selection of metallicities. Each panel shows two tables whose names start with the labels indicated in each panel. The dashed lines represent T eq from the unshielded version (" shield0"), while the solid lines are for T eq including self-shielding (" shield1"). For tables without an ISRF (left panel) even unshielded gas can cool down to temperatures below 10 4 K as soon as it becomes dense enough (n H 1 cm −3 ). Including self-shielding moves the density of the transition from the warm (T ≈ 10 4 K) to the cold (T 10 4 K) gas phase to lower densities, depending on the gas metallicity.
Most of the ISM gas in the solar neighbourhood has thermal pressures of 3 < log P/k B [ K cm −3 ] < 4 Table 7. Line identifiers and wavelengths for the lines in the line emission files (see Sec. 4.2). The wavelengths are in units ofÅ (A), cm (c), and µm (m). "Blnd" refers to a blended line with multiple components and is listed directly after its main component(s). The line with the cloudy identifier "Cool" at 1215Å is the contribution of collisional excitation to Lyα cooling.

Line
Wavelength  (Jenkins & Tripp 2011). The shaded areas in Fig. 12 indicate this range of pressures. For solar metallicity, the minimum equilibrium temperature of the cold phase is ≈ 30 K in table UVB dust1 CR1 G1 while the warm phase (unshielded gas) has an equilibrium temperature of ≈ 10 4 K. For tables including the reduced ISRF (right panel of Fig. 12), a multi-phase medium can be established between the unshielded and the self-shielded gas for high metallicity (Z ≥ 0.1 Z ) gas over a large range of densities. The redshift dependence of T Eq is shown in Fig. 13 for solar (left set of panels) and primordial abundances (right set of panels). As in Fig. 12, the left panel does not include the ISRF and T Eq depends on redshift over the full density range. In the right panels of Fig. 13, which include the redshiftindependent ISRF, the thermal equilibrium temperatures of different redshifts converge at high densities.
T Eq increases with density for the lowest densities in the tables (see Figs. 12 and 13). This relation is caused by the combination of a photo heating rate Λ heat /n 2 H that is independent of density for highly (photo-)ionized gas and Compton cooling rate Λ cool,Compton /n 2 H that increases with decreasing density. The redshift dependence at these densi-ties is dominated by the steep dependence of Compton cooling with redshift. In practice the thermal equilibrium temperature from radiative processes is irrelevant for the lowest densities as the cooling time is longer than the Hubble time and adiabatic cooling (Hubble expansion) dominates over radiative cooling. At higher densities (log n H [ cm −3 ] −5 for z = 0) hydrogen and helium recombination cooling dominates between log T[ K] ≈ 4 − 6 and T Eq decreases with density.
The thermal equilibrium curves in Figs. 12 and 13 are reminiscent of the classical two phase structure of the neutral ISM, as described in Wolfire et al. (1995) (hereafter: W95). As it can be tempting to directly compare the thermal equilibrium curves from this work to those of W95, we discuss here the key differences in both the general approach as well as the interpretation of the results.
W95 describe the ISM densities where warm neutral ISM gas can be in pressure equilibrium with cold neutral gas. Based on the idea of a dense, cold gas cloud embedded in more tenuous, warm (10 4 K) gas, the pressure equilibrium is determined at the edge of the cold gas cloud. They assume that the ISRF is shielded by a constant column density (N w = 10 19 cm −2 in their fiducial model), where N w is the typical column density of the warm medium. In this work, the shielding from the ISRF corresponds to the self-shielding of the cold gas cloud (i.e. from its edge to its centre). The shielding column N sh therefore varies with the assumed size of the gas cloud, which is a function of gas density and temperature (or pressure, as N sh ∝ (n H T) 0.5 , see Sec. 2.1).
In addition to the varying shielding column, also the radiation field intensity and the CR rate change along the thermal equilibrium curves. Higher gas pressures lead to larger column densities which are observed to result in higher SFR densities and subsequently a stronger ISRF (see Sec. 2.2.2). While Wolfire et al. (2003) vary the ISRF to account for radial variations within the MW disk, the ISRF is constant with density for each of their equilibrium functions. Fig. 14 illustrates the variations in the shielding column (top panels) and ISRF intensity (bottom panels). The coloured lines refer to the values at the thermal equilibrium temperatures for selected metallicities at z = 0 and the grey horizontal lines show the fiducial values of the W95 model, with log N sh [ cm −2 ] = 19 and the radiation field from Draine (1978) with J D78 = 6.6 × 10 5 photons cm −2 s −1 sr −1 eV −1 and therefore J D78 /J 0 = 1.5.
The column densities for the tables that include selfshielding span several orders of magnitudes for typical ISM densities (n H 10 −2 cm −3 ). The radiation field spans an even larger range as J ∝ Σ SFR ∝ N sh 1.4 . Higher pressure environments are therefore self-consistently exposed to more intense radiation fields. W95 and Wolfire et al. (2003) focus on the MW ISM and therefore assume solar metallicity. Recently, Bialy & Sternberg (2019) extended their model to lower metallicity gas and included molecular hydrogen heating and cooling processes. They found that a multi-phase medium requires higher pressure for lower metallicities if the ISRF is the same as that assumed in solar metallicity galaxies. For extremely low metallicities (Z 10 −5 Z ) and for a total CR ionization rate per hydrogen atom of 10 −16 s −1 they found that the ISM would be single-phase, as the equilibrium temperature changes smoothly from ≈ 10 4 K to ≈ 600 K. Despite the above discussed differences, we find qualitatively similar behaviour for the thermal equilibrium temperatures of the recommended table UVB dust1 CR1 G1 shield1. The  colour scale in Fig. 15 shows the minimum pressure for unshielded (" shield0") and shielded (" shield1") gas to be in pressure equilibrium (where the former corresponds to the warm phase: log T warm [ K] > 3.5 and the latter to the cold phase: log T cold [ K] < 3). An overall trend with metallicity (i.e. lower metallicity -higher minimum pressure) is clearly visible for every redshift.
Note that for very low-metallicity gas (log Z/Z −2) the lack of metal-line cooling keeps the temperature at 800 K, even for very high density gas (n H = 10 6 cm −3 , see Fig. 12). This gas is mainly heated by CRs and by the vibrational and rotational energy of molecular hydrogen, as it absorbs UV photons.

Phase transitions
Ionized -neutral hydrogen: For each radiation field, the transition from mostly ionized to mostly neutral hydrogen marks where self-shielding becomes important. Assuming that most of the gas follows the thermal equilibrium temperature, the tables directly provide typical hydrogen species fractions for each gas density. Fig. 16 shows f H i for solar metallicity and redshift z = 0 for tables UVB dust1 CR0 G0 (UVB only, left panel) and UVB dust1 CR1 G1 (UVB, CRs and ISRF, right panel). For each radiation field, f H i is presented for both unshielded (blue) and self-shielded (red) gas. For self-shielded gas, the neutral atomic hydrogen fraction is displayed both as a column density fraction f H i = N H i /N H (solid lines) and as a volume density fraction at the end of the shielding column f H i = n H i /n H (dashed lines).
For the UVB only tables (UVB dust1 CR0 G0), we compare f H i from this work to the widely used fitting function from R13, for which n H i /n H is a function of gas density, temperature and photoionization rate Γ phot . The dotted line in the left panel of Fig. 16 indicates (n H i /n H ) R13 for the same (equilibrium) temperatures. The characteristic density above which f H i for self-shielded gas deviates  Rahmati et al. (2013) for the same T eq as the self-shielded gas. from unshielded gas is the same in the shielded table and the R13 fitting function. While the transition is very steep in the self-shielded table presented in this work, the fitting function from R13 has a more gradual increase. The difference here may be that f H i from the tables assumes that all gas follows the thermal equilibrium temperature and a fixed n H − N H relation, while f H i from R13 is a fit to cosmological radiative transfer simulations. Note, that in the tables f H i can decrease again towards higher densities as hydrogen becomes molecular. Molecular hydrogen is not included in R13.
Atomic -molecular hydrogen: The transition between atomic and molecular hydrogen is summarised in Fig. 17 for two self-shielded tables (ending with "shield1"). The left panel does not include an interstellar radiation field, while the right panel includes both CRs and the ISRF. In each panel, the H 2 fraction is shown for the thermal equilibrium temperature. The x-axis shows the column density of atomic plus molecular hydrogen through the gas cloud where N sh as well as the H i and H 2 fractions are taken from the tables along the thermal equilibrium temperature (hdf5 group /ThermEq/). For solar metallicity, the results can be compared to observations of H 2 fractions in the Milky Way in sight-lines towards stars (e.g. with Copernicus, Savage et al. 1977) or AGN (e.g. the FUSE halo survey, Gillmon et al. 2006). Inside the Milky Way Galaxy, the transition from atomic to molecular hydrogen (from log f H2 ≡ log 2N H2 /(N H i + 2N H ) −3 to log f H2 −1) has been measured to occur around log(N H i +2N H2 )[ cm −3 ] = 20.70 (Savage et al. 1977), while for high-latitude sight-lines, hydrogen becomes molecular at sli- Figure 17. Molecular hydrogen fraction f H2 , defined as 2N H2 /(N H i + 2N H2 ) (column density fraction, solid lines) or 2n H2 /(n H i + 2n H2 ) (volume density fraction at the end of the shielding column, dashed lines) assuming thermal equilibrium temperatures for different radiation fields (different panels) and different metallicities (different line colours) for the self-shielded tables. The grey band brackets the column densities for the observed H i-H 2 transition in the solar neighbourhood, from log(N H i + 2N H2 )[ cm −2 ] = 20.38 (FUSE halo survey, Gillmon et al. 2006 Gillmon et al. 2006;Copernicus, Savage et al. 1977).
Using these tables in hydrodynamical simulations, we find that the gas temperature of individual resolution elements scatters up from the thermal equilibrium temperature but rarely goes below T eq . Higher temperatures for constant density lead to an increase in N ref and typically to a decrease in f H2 . In addition, the tables assume a constant density for each column density, while a sightline in a simulation can go through a variety of gas volume densities. We find that the atomic to molecular transition typically does not scatter to lower column densities in galaxy-scale simulations (Ploeckinger et al. in prep). Therefore, table UVB dust1 CR1 G1 shield1 matches the column density where the H i-H 2 transition is observed to occur.

Ion fractions
Each table contains the ion fractions of 11 selected elements (see Table 6). Examples for table UVB dust1 CR1 G1 shield1 at z = 0 and for solar metallicity gas in thermal equilibrium are shown for carbon and oxygen (Fig. 18). Ion fractions of all 11 tabulated elements can be found at the project webpages (see Sec. 1 for links).
As hydrogen is the most abundant element and it efficiently absorbs photons with energies higher than 13.6 eV, elements with minimum ionization energies at which the hydrogen photoionization cross section is large (e.g. helium: E i = 24.59 eV, nitrogen: E i = 14.53 eV, oxygen: E i = 13.62 eV, and neon: E i = 21.56 eV), are pre-dominantly neutral in selfshielded gas (see right panel of Fig. 18 for oxygen). Other elements, such as carbon (E i = 11.26 eV) and magnesium (E i = 7.65 eV) are singly ionised over a large range of volume densities (see left panel of Fig. 18 for carbon).

H2 and CO abundance
Each tables contains the abundances of CO and H 2 , both in the last cloudy zone as well as integrated over the full shielding column. Fig. 19 shows these abundances for solar metallicity gas at z = 0 at its thermal equilibrium temperature for self-shielded tables with and without ISRF and CR rates: UVB dust1 CR0 G0 shield1 and UVB dust1 CR1 G1 shield1. While the H2 abundance decreases when the ISRF and CRs are added, at high densities (log n H [ cm −3 ] > 4) the CO abundance is largest when including an ISRF.
A closer look at the molecule abundances is presented in Fig. 20. Here, the dependence of their abundances on the depth into the gas cloud (column density N H ) is visualised for the two tables from Fig. 19. For this, an individual grid point (at log n H [ cm −3 ] = 4, see figure caption) is rerun with the additional output 8 for the volume densities of 95 molecule species in each zone. Note that all these molecules are part of the chemistry for all grid points, but this additional output is not stored for the full table. The left panel of Fig. 20 (no ISRF, no CR) explains why gas phase CO is less abundant at the centre of the cloud (here at N H ≈ 10 22 cm −2 ): the majority of CO molecules are condensed into dust grains ("COgrn"). This is also the case for the OH and H 2 O molecules ("OHgrn"and "H2Ogrn"). For the stronger radiation field (right panel) molecule fractions peak at much higher column densities and therefore close to the centre of the gas cloud.
The resulting fraction N CO /N C for the full shielding column is comparable for both radiation fields at a density of log n H [ cm −3 ] = 4 (see Fig. 19 and triangles in Fig. 20) but without an ISRF, most of the CO is depleted onto dust grains, which explains the reduced gas phase CO fraction for high densities log n H [ cm −3 ] 4 for the table without the ISRF (UVB dust1 CR0 G0 shield1).

EXAMPLE RESULTS: EMISSION LINES
Out of the 183 emission lines listed in Table 7 we show a few examples of soft X-ray emission lines in Fig. 21 and farinfrared/(sub-)mm emission lines in Fig. 22. In both figures the normalized emergent emissivities, /n 2 H , for the last cloudy zone of the shielding column ( Vol ) as well as the emissivity calculated by dividing the line intensity of the full shielding column by the length of the shielding column ( Col ) from the fiducial model "UVB dust1 CR1 G1 shield1" are displayed.    Fig 12), which explains the features (most noticeable in the C i fraction) around this density. Figure 19. CO and H 2 abundances for gas in thermal equilibrium at redshift z = 0 and solar metallicity. Column density fraction of H2 (2N H2 /N H , blue lines) and CO (N CO /N C , red lines) for two tables that include self-shielding (no ISRF and no CRs: dashed lines, including ISRF and CR: solid lines). The fraction of CO decreases with density for log n H [ cm −3 ] 3 for the weaker radiation field (dashed lines) as CO is increasingly depleted onto dust grains (see Fig. 20).
This has been traced down to differences between the output produced by the cloudy commands "save lines emissivity" and "save line list absolute" and even persists for one-zone models (e.g. "UVB dust1 CR1 G1 shield0"). The unexpected behaviour of cloudy has been posted to the cloudy group (see https://cloudyastrophysics.groups. io/g/Main/message/4301 for details). This is unlikely affecting any simulation results as this issue only appears at very low densities where the emissivities are very small.
The results for selected CO lines as well as the C ii emission line at 157 µm are presented in Fig. 22 for densities of log n H [ cm −3 ] = 2 (left panel) and log n H [ cm −3 ] = 4 (right panel). For higher densities the CO fraction increases and therefore the normalized emissivity of C ii decreases. At log n H [ cm −3 ] = 4 (right panel) CO is collisionally excited to higher energy levels, leading to brighter emission lines of Figure 20. Abundances of selected molecules (and their depletion onto dust grains, indicated by the suffix "grn") along the shielding column at an individual grid point of two tables (left panel: UVB dust1 CR0 G0 shield1, right panel: UVB dust1 CR1 G1 shield1) for the respective equilibrium temperatures at the end of the shielding columns (log T [K] = 1.00, left and 1.41, right) of z = 0, solar metallicity gas with a volume density of log n H [ cm −3 ] = 4. The H 2 abundance n H2 /n H , has been reduced by a factor of 10 3 to fit on the plot. For guidance, the horizontal dashed line indicates the abundance where all H is in H 2 . As seen in Fig. 19, the CO abundance at the last zone (indicated by a red triangle) for this density is comparable for both radiation fields, but the distribution throughout the column is very different.
As CO forms mainly deep into the shielding column (see e.g. right panel of Fig. 20), the emissivities calculated at the last zone ( Vol , solid lines) are not expected to be identical to the shielding column averaged emissivity ( Col , dotted lines) but they converge to the same values for temperatures log T[K] 4, where the shielding column becomes small. 6, photo-ionization dominates over collisional ionization and the line emissivity therefore depends on density (compare different panels). The differences between Vol and Col for log n H [ cm −3 ] = −6 (right panel) are unexpected as the gas is unshielded at this density and the variations within the small shielding column are negligible. This is likely an artefact in cloudy (see text for details).

SUMMARY
We use cloudy version 17.01 to tabulate the properties of unshielded and self-shielded gas as it is exposed to metagalactic and interstellar radiation fields. The tables cover a sufficiently large range in redshift (z = 0 -9), temperature (log T[K] = 1 -9.5), metallicity (log Z/Z = −4 -+0.5, Z = 0), and density (log n H [ cm −3 ] = −8 -+6) to enable users to use one set of tables for all gas phases (ionised, atomic, molecular) both during a simulation (i.e. cooling and heating rates) as well as for analysing the simulation output in post-processing (i.e. species fractions, line emissivities).
The gas column density is an important factor to estimate both the local star formation rate as well as the attenuation of the radiation field. The reference column density N ref in this work is based 9 on the Jeans column density N J , which is a typical scale for self-gravitating gas. As gas surface density is observed to correlate with the star formation rate surface density (Eq. 10), we add an interstellar radiation field (ISRF) and cosmic rays proportional to N 1.4 ref to the redshift-dependent UV/X-ray background. The UV (and X-ray) background is based on the model from FG20 but before hydrogen and helium reionization are complete, we attenuate the FG20 spectra to better match the effective rates 9 N ref equals N J for ISM gas but in addition includes a smooth transition to optically thin gas (see Eq. 6). As molecules form in self-shielded gas and their abundance and level population depend on the shielding column density, the emissivity at the last cloudy zone Vol (i.e. the most shielded part of the gas column, solid lines) of the selected CO lines is typically larger than the average line emissivity of the shielding column Col (dotted lines). For high densities (e.g. log n H [ cm −3 ] = 4, right panel) the emissivity of C ii ("C 2") is lower at high shielding column densities and therefore the column averaged emissivity is larger than the emissivity at the last cloudy zone ( Col > Vol ). inferred from observations by FG20 and therefore the observed electron scattering optical depth of τ e = 0.054 (Planck Collaboration et al. 2018).
Dust grains are an important catalyst for the formation of H 2 and other molecules and also contribute to selfshielding. A constant dust-to-metals ratio is assumed for star formation rate surface densities above the solar neighbourhood value, while for lower N ref , the dust-to-metal ratio scales like the ISRF (∝ N 1.4 ref , Eq. 18). The gas phase abundances of individual metals that are depleted onto dust grains are reduced accordingly.
At the centre of a gas cloud, the radiation field is attenuated by gas with a shielding column of N sh = 0.5N ref and the dust within. This is modelled by passing the incident radiation field through a gas shielding column, N sh . A large number of gas properties (see Table 6 for a full list) is stored in hdf5 datasets at the unattenuated (unshielded; tables with shield0) edge of the gas clump, as well as at the gas cloud centre (self-shielded; tables with shield1). Selected properties (e.g CO and H species fractions, line emissivities) are also available as averages over the shielding column. In the tables for self-shielded gas, all elements (also all metals) contribute to the self-shielding and in turn respond to the attenuated radiation field. Some datasets (e.g. cooling and heating rates in Figs. 11 and 9, ion fractions in Figs. 18) are presented for individual slices through the multi-dimensional hyper cuboid and we provide an easy to use graphical user interface to explore the full range.
As an example application, we showed characteristic Hdf5 files with cooling and heating rates, ion fractions, free electron fractions and molecule abundances for each model listed in Table 5. The content of each datafile is listed in Table 6 and described in Sec. 4.1. <modelname> lines.hdf5 Hd5 files for models UVB dust1 CR1 G1 shield0 and UVB dust1 CR1 G1 shield1 containing the emissivities for the emission lines listed in Table 7 and described in Sec. 4.2. modFG20 cloudy.ascii The modified FG20 spectra in the cloudy UVB format. See Appendix C for how to use this UVB with cloudy.
densities for the most important phase transitions (H 2 to H i, Fig. 16, and H i to H ii, Fig. 17) and compared the results to values from the literature. For the transition from ionised to neutral hydrogen, we compared our results to the fitting function from R13. Fig. 16 illustrates that self-shielding becomes important at very similar densities, although the increase in f H i with density is steeper in this work than in R13. The predicted transition between neutral and molecular hydrogen was compared to observations of individual sight lines within the MW disk and towards the MW halo, measuring both the H i and H 2 column densities (Fig. 17). For our fiducial ISRF and CR rates (UVB dust1 CR1 G1 shield1) the observations are very well matched.
Two tables (UVB dust1 CR1 G1 shield0 and UVB dust1 CR1 G1 shield1) include line emissivities for 183 selected lines. A full overview of the data files released with this work is in Table 8. In addition to the hdf5 data files we also provide a set of C routines that interpolate the tables and return the total net cooling rate for a given redshift, gas density, temperature, and metallicity or abundance ratio, which can be implemented in a hydrodynamics code, and a set of python3 routines that read the tables and produce plots via an easy to use graphical user interface. All data files and links to the github repositories for the python and C packages can be found at http://radcool.strw.leidenuniv.nl/ and https://www.sylviaploeckinger.com/radcool.
Finally, practical information about how to use the tables in simulations, either alone or coupled to a nonequilibrium network, such as chimes (Richings et al. 2014a,b) or grackle (Smith et al. 2017), as well as how to reproduce our results with the public code cloudyv17.01 are provided in the appendix.
If the metal cooling rates are scaled as Λ comb = Λ table (n e,comb /n e,table ) the cooling rates Λ comb , while calculated assuming ionization equilibrium (as in Λ table ), capture many of the non-equilibrium effects for metal cooling. For the heating rates this additional scaling is not necessary as Λ heat scales with n H .

APPENDIX B: INCLUDING THE EFFECTIVE FG20 BACKGROUND RADIATION FIELD
Before and during reionization, the photo-ionization rates calculated from the redshift-dependent UV / X-ray spectra from FG20 differ from their tabulated "effective" photoionization rates that are modelled to match the observed electron scattering optical depth of τ e = 0.054 (Planck Collaboration et al. 2018). As cloudy uses a full spectrum as input radiation field, we modify the spectra from FG20 to match their "effective" photo-ionization rates. Assuming that before H i reionization, the H i ionizing radiation is attenuated by neutral hydrogen gas with a column density N H i (and before He ii reionization, the He ii ionizing radiation is attenuated by singly ionized helium gas with column density N He ii ), the new spectrum J modFG20 are present in each run (e.g. "table ISM" is only included for tables with the ISRF).
There is no direct analogue in cloudy v17.01 to scale the metal depletion, as is done in this work. But the individual depletion factors can be set manually at the bottom of file abund.cpp in the cloudy source code to f gas,i = 1− f dust,i , with f dust,i from the dataset Depletion.
For the UVB used here (modified FG20), replace the hm12_galaxy.ascii file in cloudy's data folder with the provided modFG20_cloudy.ascii file, but keep the original filename (hm12_galaxy.ascii). As explained in Sec. 2.3 the dissociation rate of H 2 by CRs is re-normalized to better match the UMIST values. This has been done by changing one line in the cloudy file mole_reactions.cpp from newreact("H2,CRPHOT=>H,H","h2crphh",1.,0.,0.) to newreact("H2,CRPHOT=>H,H","h2crphh",1.746e-2,0.,0.). This paper has been typeset from a T E X/L A T E X file prepared by the author.