# The effect of non-equilibrium metal cooling on the interstellar medium

Pedro R. Capelo, Stefano Bovino, Alessandro Lupi, Dominik R. G. Schleicher, Tommaso Grassi

► To cite this version:

Pedro R. Capelo, Stefano Bovino, Alessandro Lupi, Dominik R. G. Schleicher, Tommaso Grassi. The effect of non-equilibrium metal cooling on the interstellar medium. Monthly Notices of the Royal Astronomical Society, 2018, 475 (3), pp.3283-3304. 10.1093/mnras/stx3355. hal-01807432

HAL Id: hal-01807432

https://hal.science/hal-01807432v1

Submitted on 4 Jul 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# The effect of non-equilibrium metal cooling on the interstellar medium

Pedro R. Capelo,<fnref n="1" />* Stefano Bovino,<fnref n="2" />,<fnref n="3" /> Alessandro Lupi,<fnref n="4" /> Dominik R. G. Schleicher<fnref n="2" /> and Tommaso Grassi<fnref n="5" />,<fnref n="6" />,<fnref n="7" />

Accepted 2017 December 23. Received 2017 December 7; in original form 2017 October 2

## ABSTRACT

By using a novel interface between the modern smoothed particle hydrodynamics code GASOLINE2 and the chemistry package KROME, we follow the hydrodynamical and chemical evolution of an isolated galaxy. In order to assess the relevance of different physical parameters and prescriptions, we constructed a suite of 10 simulations, in which we vary the chemical network (primordial and metal species), how metal cooling is modelled (non-equilibrium versus equilibrium; optically thin versus thick approximation), the initial gas metallicity (from 10 to 100 per cent solar), and how molecular hydrogen forms on dust. This is the first work in which metal injection from supernovae, turbulent metal diffusion, and a metal network with non-equilibrium metal cooling are self-consistently included in a galaxy simulation. We find that properly modelling the chemical evolution of several metal species and the corresponding non-equilibrium metal cooling has important effects on the thermodynamics of the gas, the chemical abundances, and the appearance of the galaxy: the gas is typically warmer, has a larger molecular-gas mass fraction, and has a smoother disc. We also conclude that, at relatively high metallicity, the choice of molecular-hydrogen formation rates on dust is not crucial. Moreover, we confirm that a higher initial metallicity produces a colder gas and a larger fraction of molecular gas, with the low-metallicity simulation best matching the observed molecular Kennicutt–Schmidt relation. Finally, our simulations agree quite well with observations that link star formation rate to metal emission lines.

**Key words:** astrochemistry – molecular processes – methods: numerical – ISM: molecules – galaxies: evolution – galaxies: ISM.

## 1 INTRODUCTION

Numerical simulations of galaxy formation and evolution, both cosmological and idealized, have been extremely important in shaping our understanding of structure formation for both the dark matter (DM) component and the baryonic collapsing objects within. In particular, DM-only cosmological simulations have been very successful at reproducing observations of the large-scale structure of the Universe, helping the cold DM model with a cosmological

constant to become the most widely accepted paradigm for structure formation (e.g. Springel et al. 2005b).

Simulating the baryonic content of the Universe, however, has proven a lot more difficult, as hydrodynamical simulations have, for a long time, failed to produce consistent results and to reproduce galaxies that match observations (see e.g. Schaye et al. 2010; Scannapieco et al. 2012; Schaye et al. 2015).

This was due to a variety of factors, linked in part to the use of different hydrodynamics algorithms [e.g. smoothed particle hydrodynamics (SPH) and adaptive mesh refinement (AMR) methods; see e.g. Agertz et al. 2007] and, most importantly, to how small-scale physical processes were modelled. Indeed, due to the huge dynamic range involved in the simulations and to limitations in

<footnote n="1">¹Center for Theoretical Astrophysics and Cosmology, Institute for Computational Science, University of Zurich, Winterthurerstrasse 190, CH-8057 Zürich, Switzerland</footnote>

<footnote n="2">²Departamento de Astronomía, Facultad Ciencias Físicas y Matemáticas, Universidad de Concepción, Av. Esteban Iturra s/n Barrio Universitario, Casilla 160, Concepción, Chile</footnote>

<footnote n="3">³Hamburger Sternwarte, Universität Hamburg, Gojenbergsweg 112, D-21029 Hamburg, Germany</footnote>

<footnote n="4">⁴Institut d'Astrophysique de Paris, Sorbonne Universités, UPMC Univ. Paris 06 et CNRS, UMR 7095, 98bis boulevard Arago, F-75014 Paris, France</footnote>

<footnote n="5">⁵University Observatory Munich, Scheinerstr. 1, D-81679 Munich, Germany</footnote>

<footnote n="6">⁶Excellence Cluster Origin and Structure of the Universe, Boltzmannstr. 2, D-85748 Garching, Germany</footnote>

<footnote n="7">⁷Centre for Star and Planet Formation, Niels Bohr Institute & Natural History Museum of Denmark, University of Copenhagen, Øster Voldgade 5-7, DK-1350 Copenhagen K, Denmark</footnote>

<footnote>* E-mail: pcapelo@physik.uzh.ch</footnote>
computer power, many processes have to be modelled with so-
called subgrid recipes, which heavily depend on resolution. Large
uncertainties come from, e.g. modelling stellar evolution and the
initial mass function (IMF) in stellar population synthesis models
(e.g. Conroy, Gunn & White 2009; Gonzalez-Perez et al. 2014),
supernova (SN) and supermassive black hole (BH) feedback (e.g.
Dalla Vecchia & Schaye 2008, 2012; Bellovary et al. 2010; Durier &
Dalla Vecchia 2012), star formation (SF; e.g. Hopkins, Narayanan &
Murray 2013; Braun & Schmidt 2015), and magnetic fields (e.g.
Pakmor & Springel 2013; Rodenbeck & Schleicher 2016).

In recent years, there has been much improvement, due in part to
the introduction of modern-SPH codes, which has alleviated several
of the issues raised in Agertz et al. (2007), and to a strong effort
in developing subgrid recipes (e.g. Keller et al. 2014; Kimm et al.
2015; Hopkins et al. 2017).

A crucial aspect of the modelling of galaxies is the proper mod-
elling of the gas chemistry and cooling which, so far, has received
relatively limited attention. The cooling rate of gas has important
dynamical effects, as it controls, e.g. the accretion mode of gas in
DM haloes, the amount of SF in galaxies, and accretion on to BHs.

Several galaxy simulations assume that gas is in (collisional and/or photo) ionization equilibrium and employ cooling tables or functions (often constructed with CLOUDY; Ferland et al. 2013) to compute gas cooling rates, based on gas density, temperature, and redshift. A few studies follow a few atomic primordial species to compute non-equilibrium cooling, but still employ the use of tables for molecular hydrogen and metals (e.g. Shen, Wadsley & Stinson 2010).

Only during the past decade, non-equilibrium chemical networks
with molecular hydrogen have been incorporated in galaxy simula-
tions. Gnedin, Tassis & Kravtsov (2009; see also Gnedin & Kravtsov
2011 for a revised version) were the first to self-consistently im-
plement an eight-species primordial chemical model (with H I, H II,
H2 I, H2 II, H −, He I, He II, and He III)<fnref n="1" /> to follow the rate of formation
and destruction of molecular hydrogen, applying it to cosmological
AMR simulations of galaxies down to redshift z = 4. This model, in
which SF is directly linked to the abundance of molecular hydrogen,
was successfully used to study in detail the Kennicutt–Schmidt (KS;
Schmidt 1959, 1963; Kennicutt 1989, 1998) relation (Gnedin &
Kravtsov 2010; Feldmann, Gnedin & Kravtsov 2011; Gnedin &
Kravtsov 2011).

Christensen et al. (2012) later implemented a model related to
that of Gnedin et al. (2009) to run isolated SPH simulations of
a Milky Way-sized galaxy and cosmological SPH simulations of
a dwarf galaxy to z = 0, finding a colder, denser, and clumpier
interstellar medium (ISM) with respect to ‘standard’ simulations.
Tomassetti et al. (2015) also implemented a similar model and ran
cosmological AMR simulations of a Milky Way-sized galaxy to
z = 2, to assess the effect of using different SF recipes.

In all the above studies, only primordial species are evolved.
Metal cooling is therefore still modelled assuming equilibrium.
Richings & Schaye (2016) perform SPH simulations of isolated
galaxies with a static DM potential, using the CHIMES chemistry
solver module described in Richings, Schaye & Oppenheimer
(2014a,b). They follow the evolution of 157 chemical species and
the non-equilibrium cooling from the ionization states of 11 el-
ements, to assess the effects of metallicity and ultraviolet (UV)

radiation in galaxies. In their implementation, however, they do
not model metal injection from SNae nor turbulent metal diffusion,
therefore keeping the gas metallicity fixed throughout the entire
simulation. This assumption, though, is likely to fail in star-forming
galaxies, wherein the gas metallicity can increase even by one order
of magnitude (see Section 3.2).

Hu et al. (2016; see also Hu et al. 2017) go a step further in that
direction, by modelling metal injection from SNae in their isolated
SPH simulations of dwarf galaxies. However, they track only six
species (H2 I, H II, CO, H I, C II, and O I), of which only the first
three are explicitly followed.

Further investigations are clearly needed to explore the treat-
ment of metal-line cooling (following the evolution of more metal
species), the comparison of equilibrium versus non-equilibrium
cooling, the effect of UV radiation models, and dust physics. Lupi
et al. (2018), using initial conditions and settings similar to ours,
study in detail the effect of having a local UV field from young
stars and several ways to model it (see also Nickerson, Teyssier &
Rosdahl 2018). For this reason, in this work, we focus more on the
role of metals as coolants and on how molecular hydrogen forms
on dust.

By using a novel interface between a modern-SPH code (GASOLINE2; Wadsley, Keller & Quinn 2017) and a chemistry package (KROME; Grassi et al. 2014), we explicitly follow nine primordial and seven metal species, model the corresponding non-equilibrium cooling, compute both metal injection from SNae and turbulent metal diffusion, and explore models of H2I formation on dust.

Properly modelling the gas chemistry is not only important to
accurately compute gas cooling, but is also relevant in the context
of observational efforts: computing the abundance of molecular hy-
drogen is important to compare the simulation results to observed
relations between SF rate (SFR) and H2I surface densities (the
so-called molecular KS relation; Bigiel et al. 2008, 2010). More-
over, evolving a chemical network that includes metals allows us
to use C I, C II, and O I as tracers of the ISM (e.g. Malhotra et al.
2001; Cormier et al. 2015; Michałowski et al. 2016; Pallottini et al.
2017a,b) and of SFR (e.g. De Looze et al. 2014), giving us a better
theoretical understanding on how all these species trace different
components of the ISM and of star-forming regions. This is par-
ticularly timely, given the wealth of far-infrared (FIR) data coming
from recent (e.g. the Herschel Space Observatory; Pilbratt et al.
2010; low redshift) and current (e.g. the Atacama Large Millime-
ter/submillimeter Array – ALMA; ALMA Partnership et al. 2015;
high redshift) facilities.

In Section 2, we describe in detail the new code interface, the
initial conditions, and the suite of simulations. In Section 3, we
present the results from our suite, focusing on the role of non-
equilibrium metal cooling, the dependence on metallicity and on
the H2I formation model, and on the role of metal species as tracers
of the ISM and SFR. We summarize and conclude in Section 4.

## 2 NUMERICAL SET-UP

For all the simulations presented in this work, we interfaced
two publicly available codes<fnref n="2" /> – the *N*-body SPH code GASOLINE2
(Wadsley et al. 2017; an extension of the *N*-body tree code PKDGRAV;
Stadel 2001) and the chemistry package KROME (Grassi et al. 2014)

<footnote n="1">¹ In this work, we adopt the astronomy notation in which, e.g. H, H I, and H2 I indicate total hydrogen, atomic neutral hydrogen, and molecular neutral hydrogen, respectively.</footnote>

<footnote n="2">² Unless otherwise stated, when we simply write ‘the code’, we refer to the ensemble of the two codes: GASOLINE2 (http://gasoline-code.com/) and KROME (http://kromepackage.org/).</footnote>
– to accurately follow the hydrodynamical and chemical evolution of an isolated galaxy.

The modern-SPH GASOLINE2 code differs from the original GASOLINE code (Wadsley, Stadel & Quinn 2004) in many respects (see Wadsley et al. 2017 for a thorough description). The main additions employed in the simulations of this work, compared to, e.g. those performed in Capelo et al. (2015) with a previous version of the code for galaxy simulations, are the use of an improved SPH smoothing kernel (the Wendland C² kernel; Wendland 1995; Dehnen & Aly 2012; Keller et al. 2014) that does not suffer from the pairing/clumping instability (Schuessler & Schmitt 1981), an increased number of neighbours (64), the geometric-density-average force expression (Keller et al. 2014; see also Monaghan 1992; Ritchie & Thomas 2001), found to minimize numerical surface tension, and the use of a pressure floor (Robertson & Kravtsov 2008; Brook et al. 2012; Roškar et al. 2015), to ensure the correct fragmentation behaviour of the gas.

The chemistry package KROME (Grassi et al. 2014) is a flexible code that can be embedded as a library in any hydrodynamics code and that, given any chemical network, generates an optimized FORTRAN code to solve the corresponding system of coupled ordinary differential equations (ODEs) and follow the time-dependent evolution of the chemical species as well as the gas temperature. Additionally, the package provides a variety of modules to model several physical processes, including radiative and thermochemical cooling/heating, photochemistry, and dust physics (e.g. formation of molecular hydrogen by catalysis on dust, gas cooling by dust, and dust evaporation; see also Grassi et al. 2017). KROME has already been successfully interfaced with a variety of hydrodynamics codes, such as the AMR codes ENZO (Bryan et al. 2014; see e.g. Bovino et al. 2014), RAMSES (Teyssier 2002; see e.g. Katz, Sijacki & Haehnelt 2015), and FLASH (Fryxell et al. 2000; see e.g. Seifried & Walch 2016; Körtgen et al. 2017), and, very recently, the mesh-less finite-mass code GIZMO (Hopkins 2015; see Lupi et al. 2018). An updated list of applications and interfaced codes can be found at the KROME website.

In this work, the first application of a working interface between KROME and an SPH implementation of the equations of hydrodynamics is presented.

## 2.1 Initial conditions

The initial conditions were built using the MAKEDISK code (Springel & White 1999; Springel, Di Matteo & Hernquist 2005a) and are identical for all the runs, except for the initial metallicity of the gas which, depending on the simulation, can be 0.1 Z⊙, 0.5 Z⊙, or Z⊙, where Z⊙ = 0.0134 (Asplund et al. 2009).

The galaxy is composed of a DM halo and of a baryonic disc and bulge. The DM halo is described by a spherical Navarro–Frenk–White (NFW; Navarro, Frenk & White 1996) profile up to the virial radius, $r_{vir}$, and by an exponentially decaying NFW profile outside $r_{vir}$ (Springel & White 1999). The DM spin and concentration parameters are 0.04 (Vitvitska et al. 2002) and 4 (Dutton & Macciò 2014; Diemer & Kravtsov 2015), respectively, with the concentration parameter setting both the DM scale radius and the exponential decay pace. The baryonic disc is composed of a stellar and a gaseous exponential disc with an isothermal sheet (Spitzer 1942; Camm 1950), which share the same initial disc scale radius, $r_d$, obtained from imposing conservation of specific angular momentum of the material that forms the disc (Mo, Mao & White 1998), and disc scale height $z_d = 0.1r_d$. The virial mass fraction of the baryonic disc is 0.02, of which 60 per cent is gaseous and 40 per cent is stellar,

typical of high-redshift galaxies (Tacconi et al. 2010). The vertical structure of the gas is governed by hydrostatic equilibrium, which sets the initial gas temperature to be between ~10² and ~10⁶ K, depending on the gas density, and we additionally impose a global gas temperature floor of 10 K. The stellar bulge is described by a spherical Hernquist (1990) profile, with a virial mass fraction of 0.004 and a scale radius $r_b = 0.2r_d$. The virial mass and radius of the galaxy are $2 \times 10^{11} M_{\odot}$ (Adelberger et al. 2005) and 45 kpc, respectively, and $r_d = 1.28$ kpc.

The (constant) gravitational softenings of the stellar, gas, and DM particles are $\epsilon_{star} = 10$, $\epsilon_{gas} = 20$, and $\epsilon_{DM} = 30$ pc, respectively. The smoothing length $h$, defined as half the distance to the furthest of the 64 neighbours, cannot be smaller than $0.1\epsilon_{gas}$. The initial particle mass for stars, gas, and DM is $3.3 \times 10^3$, $4.6 \times 10^3$, and $1.1 \times 10^5 M_{\odot}$, respectively, making the total initial number of particles ~4 × 10⁶.

The galaxy's structure is similar to that of the main galaxy in the gas-rich galaxy merger of Capelo et al. (2015), where only equilibrium chemistry was followed, and of the galaxy in Lupi et al. (2018), where the non-equilibrium chemistry of a primordial network and several radiative-transfer methods were studied. The structural parameters are typical of high redshift ($z=3$) galaxies, but the same system can be interpreted as a representative of relatively small, gas-rich low-redshift galaxies (see also Capelo et al. 2017).

## 2.2 Star formation, feedback, and metal injection

SF and stellar feedback (energy, mass, and metal injection) are modelled within GASOLINE2 and follow the recipes described in Stinson et al. (2006; see also Katz 1992; Christensen et al. 2010; Shen et al. 2010), which we summarize here.

At the beginning of the simulation, the existing stars are all set to be 2 Gyr old. New stars are allowed to form from gas with temperature $T_{gas} < T_{SF} = 10^3$ K and density $\rho_{gas} > \rho_{SF} = 100 m_H \text{ cm}^{-3}$ (Capelo et al. 2015), where $m_H$ is the hydrogen mass. When the above conditions are met, gas particles are stochastically selected to form stars according to the probability $p = (m_{gas}/m_{star})[1 - \exp(-(\epsilon_{*} \Delta t/t_{dyn}))]$, where $m_{gas}$ is the mass of the (potentially) star-forming gas particle, $m_{star}$ is the mass of the (potentially) formed star (set to half the initial gas particle mass), $\epsilon_{*}$ is the SF efficiency parameter, $\Delta t = 10^6$ yr is the SF time-scale (i.e. how often SF is computed), $t_{dyn} = (4\pi G \rho_{gas})^{-1/2}$ is the local dynamical time, and $G$ is the gravitational constant.

The above probability function implies that, on average, $dM_{*}/dt = \epsilon_{*}m_{gas}/t_{dyn}$, effectively ensuring that SF approximately follows the slope of the KS relation between SFR and gas surface densities, whereas the normalization can be matched by tuning the SF efficiency parameter. In our simulations, we first 'relax' the galaxy, i.e. we run the simulation for 0.05 Gyr using $\epsilon_{*} = 0.005$, followed by another 0.05 Gyr using $\epsilon_{*} = 0.01$, in order to counteract the fact that, at the beginning of the simulation, there is no heating from SNae which could prevent excessive gas cooling and consequent (unphysically high) bursts of SF. After the 'relaxation' period of 0.1 Gyr, we run the simulation with $\epsilon_{*} = 0.015$ (Krumholz, Dekel & McKee 2012), found to match the normalization of the KS law fairly well and to be consistent with the average SF efficiency computed as a function of $\rho_{SF}$ in Lupi et al. (2018).

Every stellar particle, being much more massive than real individual stars, represents a stellar population, assumed to cover the entire IMF described in Kroupa (2001). Stellar masses are then converted into stellar lifetimes according to Raiteri, Villata & Navarro (1996), in order to compute when and if a star explodes as an SN.
Stars with masses between 8 and 40 M⊙ can explode as a Type II SN (SNII), each injecting $E_{SN} = 10^{51}$ erg into the surrounding gas particles (using the SPH smoothing kernel) as thermal energy, according to the ‘blastwave model’ of Stinson et al. (2006), in which, in order to avoid the gas from radiating away the SN energy due to limited resolution, radiative cooling is disabled during the survival time of the hot low-density shell of the SN (McKee & Ostriker 1977). Together with energy, for each SN event, a given amount of total, iron, and oxygen mass, dependent on the progenitor mass $m_{star}$,

$$ 
\begin{aligned}
M_{\text{ejected,SNII}} &= 0.7682 (m_{\text{star}}/M_{\odot})^{1.056} M_{\odot}, \\
M_{\text{Fe,SNII}} &= 2.802 \times 10^{-4} (m_{\text{star}}/M_{\odot})^{1.864} M_{\odot}, \\
M_{\text{O,SNII}} &= 4.586 \times 10^{-4} (m_{\text{star}}/M_{\odot})^{2.721} M_{\odot},
\end{aligned} 
\quad (1) $$

is injected into the surrounding gas (Woosley & Weaver 1995; Raiteri et al. 1996).

Type Ia SNae (SNIa), occurring in evolved binary systems, are also modelled in the code, injecting into the surrounding gas $E_{SN} = 10^{51}$ erg and a fixed amount of mass and metals, independent of the progenitor mass (Thielemann, Nomoto & Yokoi 1986; Raiteri, Villata & Navarro 1996):

$$ 
\begin{aligned}
M_{\text{ejected,SNIa}} &= 1.4 M_{\odot}, \\
M_{\text{Fe,SNIa}} &= 0.63 M_{\odot}, \\
M_{\text{O,SNIa}} &= 0.13 M_{\odot}.
\end{aligned} 
\quad (2) $$

Since SNIa stars should not be clustered and should not produce large blastwaves, shutting the cooling in this case would produce too large an effect.

Stars with masses $1 < m_{star} < 8 M_{\odot}$ do not explode as SNae but release part of their mass as stellar winds, according to Kennicutt, Tamblyn & Congdon (1994) and Weidemann (1987), with the returned gas having the same metallicity of the low-mass stars.

The injection of metals from SNae changes the metallicity of the gas, which is computed as

$$ M_{\text{metals}} = 1.06 M_{\text{Fe}} + 2.09 M_{\odot}. \quad (3) $$

Another phenomenon that can change the metallicity of the gas is turbulent diffusion, modelled in the code by estimating the particle diffusion coefficients as $d = C_{\text{turb}}|S|h^2$, where S is the trace-free local shear tensor and $C_{\text{turb}} = 0.05$ (Wadsley, Veeravalli & Couchman 2008; Shen, Wadsley & Stinson 2010). The same model and coefficient are applied to thermal diffusion.

The change in metallicity, $f_{\text{metals}}$, naturally implies a change in the mass fraction of primordial elements ($f_{\text{H}}$ and $f_{\text{He}}$ for hydrogen and helium,<fnref n="3" /> respectively), and we follow the implementation of Shen et al. (2010; see also Jimenez et al. 2003), who adopt

$$ 
\begin{aligned}
f_{\text{He}} &= 0.236 + 2.1 f_{\text{metals}} \quad (\text{for } f_{\text{metals}} \le 0.1), \\
f_{\text{He}} &= 0.446[1 - 10(f_{\text{metals}} - 0.1)/9] \quad (\text{for } f_{\text{metals}} > 0.1), \\
f_{\text{H}} &= 1 - f_{\text{He}} - f_{\text{metals}}.
\end{aligned} 
\quad (4) $$

We note that, in our simulations, we never have $f_{\text{metals}} > 0.1$.

## 2.3 Thermodynamics, radiation, and chemistry

Chemistry and radiative cooling/heating are computed within KROME and follow the models described in Bovino et al. (2016). Depending on the simulation, we either follow a system of nine primordial species (‘eq-metals runs’: primordial species, i.e. H I, H II, H2 I, H2 II, H $^-$, He I, He II, He III, and e $^-$) and dust (to treat the formation of H2 I on dust grains), as well as metal-line cooling via an effective cooling function (see below), or a system with the same nine primordial species and dust plus seven metal species (‘non-eq-metals runs’: primordial species, dust, C I, C II, O I, O II, Si I, Si II, and Si III) and non-equilibrium metal cooling at low temperature. In the latter case, a linear system for the individual metal excitation levels is solved on-the-fly for the most important coolants in the ISM (C, O, and Si; see also Maio et al. 2007; Glover & Jappsen 2007). KROME employs the implicit high-order ODE solver DLSODES (Hindmarsh 1983) to integrate the system of differential equations corresponding to nine species and 40 reactions (eq-metals runs) or 16 species and 65 reactions (non-eq-metals runs). No three-body reactions were included, because of the relatively low gas density achieved in the simulations. The reactions are listed in Appendix B.

We assume a uniform extragalactic UV background (Haardt & Madau 2012) at $z = 3$, using $5 \times 10^3$ energy bins ranging from $1.237 \times 10^{-1}$ to $4.997 \times 10^7$ eV. We note that we did not include any radiative feedback from the local stellar population. However, we refer the reader to Lupi et al. (2018), in which the effects of the local radiation, together with different methods to account for radiative feedback, were studied, using a galaxy model almost identical to that of this work. In that case, the number of energy bins had to be reduced to 10, as the coupling to an on-the-fly radiative-transfer module is computationally very expensive.

To account for the attenuation of the radiation in high-density regions, where the radiation background is known to be shielded, we compute the optical depth $\tau = \sum_i \sigma_i N_i$, as the sum over every species i of the product of the photo cross-section $\sigma$ and the column density N. We approximate the column density as $N_i \sim n_i L_{\text{jeans}}$, where $n_i$ is the species number density and $L_{\text{jeans}}$ is the Jeans length (Jeans 1902),<fnref n="4" /> which was shown in Safranek-Shrader et al. (2017) to be a reasonable approximation for the shielding length.

In our simulations, we follow the formation of H2 I both in the gas phase and through catalysis on dust. In particular, H2 I is formed through gas-phase channels (e.g. reactions 7–8 and 9–10 in Appendix B) but, given the relatively high values of metallicity and dust-to-gas ratio considered in this work, the predominant formation channel (Hollenbach & McKee 1979) is by catalysis on dust grains (reaction 31). We assume that the dust-to-gas ratio, $D \equiv \rho_{\text{dust}}/\rho_{\text{gas}}$, scales linearly with metallicity as $D = D_{\odot}Z/Z_{\odot}$, where $D_{\odot} = 0.00934$ (Yamasawa et al. 2011).

We employ two different models of H2 I formation on dust. The first model (‘J75’; from Jura 1975) assumes fixed dust properties and no dependence on gas temperature, yielding the rate

$$ \frac{dn_{\text{H2I}}}{dt} = 3 \times 10^{-17} n_{\text{H}} n_{\text{H}} \frac{Z}{Z_{\odot}} C_{\rho} \text{ cm}^{-3} \text{ s}^{-1}, \quad (5) $$

<footnote n="3">³ In this work, we use $f_{\text{species}}$ and $n_{\text{species}}$ to signify the mass fraction and number density of such species, respectively. We translate from the astronomy notation, $f_{\text{H}} = X$ and $f_{\text{He}} = Y$, but use interchangeably Z and $f_{\text{metals}}$, where ‘metals’ includes all species that are not H and He.</footnote>

<footnote n="4">⁴ In the case of molecular-hydrogen self-shielding, we followed the model of Wolcott-Green, Haiman & Bryan (2011) but used as the characteristic length-scale a fixed length of 20 pc (i.e. the gravitational softening of the gas) instead of the Jeans length. We re-ran 0.2 Gyr of Run 01 using the model of Richings et al. (2014b) and the Jeans length and obtained indistinguishable results.</footnote>
where $n_H = n_{H1} + n_{HII} + n_{H^-} + 2n_{H2I} + 2n_{H2II}$ and $C_\rho$ is the clumping factor, which accounts for the missed H2I formation in the high-density regions due to limited resolution. In this work, we use a constant clumping factor which, depending on the simulation, can be either $C_\rho = 1$ (i.e. no clumping factor) or 10 (see also Lupi et al. 2018, where a variable clumping factor model is implemented).

The second model ('CS09'; from Cazaux & Spaans 2009) computes the formation rate as a function of dust type and size, dust and gas number density and temperature, and gas thermal velocity $v_{gas}$ (Hollenbach & McKee 1979; CS09). If one assumes only one type of grain for simplicity, the dust density can be written as

$$ \rho_{\text{dust}} = C_{\text{dust}} \int_{a_{\text{min}}}^{a_{\text{max}}} m(a)\varphi(a)da, \qquad (6) $$

where $m(a) = 4\pi\rho_0 a^3/3$ is the mass of a grain particle of size $a$, with $a_{\text{min}} \le a \le a_{\text{max}}$, $\varphi(a) = dn(a)/da$, where $n(a)$ is the number density of grains of size $a$, and $\rho_0$ is the bulk density of the grains (2.25 and 3.13 g cm$^{-3}$ for carbonaceous and silicates, respectively; Zhukovska, Gail & Trieloff 2008). The constant $C_{\text{dust}}$ is a normalization factor constrained by $\rho_{\text{dust}} = D\rho_{\text{gas}}$. The equation can be extended to a mix of grain types (see Grassi et al. 2017). For this work, we assumed two grain types (carbonaceous and silicates), $\varphi(a) = a^{-3.5}$ (Mathis, Rumpl & Nordsieck 1977; Draine & Lee 1984), $a_{\text{min}} = 5 \times 10^{-7}$ cm, and $a_{\text{max}} = 2.5 \times 10^{-5}$ cm; divided the size range into $N_{\text{dust}} = 20$ logarithmically spaced bins; and computed the rate of H2I formation on dust as

$$ \frac{dn_{\text{H2I}}}{dt} = \frac{\pi}{2} n_{\text{H1}} v_{\text{gas}} C_{\text{dust}} \sum_{j} \sum_{i} n_{\text{dust},ij} a_{ij}^2 \epsilon_{ij} S_i, \qquad (7) $$

where the summations are over the grain types $j = (C, Si)$ and size bins $i = (1-N_{\text{dust}})$, $n_{\text{dust},ij}$ is the number density of dust type $j$ in bin $i$, and the efficiency factors $\epsilon_{ij}$ and sticking coefficients $S_i$ both depend on $T_{\text{gas}}$ and on the size-dependent dust temperature $T_{\text{dust}}$ (Hollenbach & McKee 1979; CS09; Grassi et al. 2017).

Following the methodology proposed by Grassi et al. (2017), we tabulated the reaction rates of H2I formation on dust as a function of gas number density ($10^{-6} \le n_{\text{gas}} \le 10^6$ cm$^{-3}$) and temperature ($3 \le T_{\text{gas}} \le 10^9$ K), adopting a constant dust-evaporation temperature of $1.8 \times 10^3$ K, computing the dust opacity as in Omukai et al. (2005), and assuming zero gas opacity, since grains are the dominant contribution. Our results are rather insensitive to the specific radiation background employed for calculating the dust temperatures (Bovino et al. 2016); therefore, here we only considered the cosmic microwave background (CMB) radiation from $z=3$ for definiteness. Using the same methodology, we also created similar tables for the averaged dust temperature and for the gas cooling by dust.<fnref n="6" /> All tables were created assuming solar metallicity, therefore, cooling by dust and the rate of H2I formation on dust are then (linearly) rescaled with $Z/Z_\odot$ (and, in the case of the formation rate, multiplied by $C_\rho$) within KROME. We note that, when we employ the first method of H2I formation on dust (i.e. J75), we do not compute the dust temperature nor the gas cooling by dust. However, it was

shown in Grassi et al. (2017) that the effect of gas cooling by dust at low gas density ($\lesssim 10^7$ cm$^{-3}$) is not relevant (see also Bovino et al. 2016).

Gas cooling is an extremely important physical process that has profound effects on the evolution of a galaxy, since it has direct consequences on SF (and on other phenomena, like e.g. accretion on to BHs, not modelled here; e.g. Bellovary et al. 2013). In this work, we include several cooling processes (see Grassi et al. 2014; Bovino et al. 2016 for a thorough explanation): H1, He1, and HeII collisional excitation and ionization; HII, HeII, and HeIII recombination; HeI dielectronic recombination; Compton cooling from the CMB; Bremsstrahlung (Cen 1992); H2I roto-vibrational cooling (Bovino et al. 2016; Glover & Abel 2008; Glover 2015) and collisional dissociation (Omukai 2000); grain-surface recombination (Bakes & Tielens 1994); gas cooling by dust (Hollenbach & McKee 1979; except in the J75 runs); and metal-line cooling, as explained below.<fnref n="7" />

Metal-line cooling is dominant for gas of relatively high metallicity and low density, such as that simulated in this work. Ideally, it should be computed in non-equilibrium, as it was shown that gas cooling (isochorically or isobarically) from $10^{-7} \sim 8$ K departs from equilibrium already at $\sim 10^6$ K, either in the absence (Gnat & Sternberg 2007) or presence (Oppenheimer & Schaye 2013a) of a photoionizing background. This was confirmed also when the background is fluctuating, both in idealized (Oppenheimer & Schaye 2013b) and cosmological simulations (Segers et al. 2017). However, at high temperatures, given the large number of metal species, their ionization states, transitions, and collisional processes, it becomes expensive to properly follow a metal network (even just considering the most important metal coolants like, e.g. carbon, oxygen, and silicon) and its related non-equilibrium cooling. This is why most authors employ cooling metal tables, such as collisional ionization equilibrium tables (CIE; e.g. Sutherland & Dopita 1993) or, when assuming a given photoionizing radiation background, photoionization equilibrium (PIE) tables (see Bovino et al. 2016 for a summary of the types of tables used in the literature). The usage, however, is common also at low temperatures (i.e. $T_{\text{gas}} < 10^4$ K), where the departure from equilibrium is larger.

In this work, for $T_{\text{gas}} \ge 10^4$ K, we utilize the (appropriately translated to be used within KROME) PIE table computed by Shen et al. (2013; see also Shen et al. 2010) using CLOUDY (Ferland et al. 1998), under the assumption of an extragalactic radiation background by Haardt & Madau (2012). The tables are provided as a function of gas temperature ($10 \le T_{\text{gas}} \le 10^9$ K), density ($10^{-9} \le n_H \le 10^4$ cm$^{-3}$), and redshift ($0 \le z \le 15.1$), although in this work we consider a fixed redshift ($z=3$), and assume the default CLOUDY solar composition – containing the first 30 elements in the periodic table – except for H and He. We then linearly rescale the cooling rates with metallicity, as done for the cooling by dust. We caution that the PIE tables were constructed under the optically thin approximation (i.e. with no shielding), therefore, when using them, we overestimate metal cooling at $T_{\text{gas}} < 10^4$ K.

<footnote n="5">⁵ We compute the gas number density within GASOLINE2 as $n_{\text{gas}} = \rho_{\text{gas}}/(\mu m_{\text{H}})$, with $\mu$, the mean molecular weight of the gas, calculated as $\mu^{-1} = (f_e m_{\text{H}}/m_e) + f_{\text{H1}} + f_{\text{HII}} + 0.5(f_{\text{H2I}} + f_{\text{H2II}}) + 0.25(f_{\text{He1}} + f_{\text{HeII}} + f_{\text{HeIII}}) + Z/(A)_{\text{metals}}$, where $m_e$ is the electron mass and $\langle A \rangle_{\text{metals}} = 17.6003$ (see also Bovino et al. 2016; there is a typo in their equations 28 and 30).</footnote>

<footnote n="6">⁶ Since usually $T_{\text{dust}} < T_{\text{gas}}$, this is typically called cooling. However, one can have gas heating by dust when $T_{\text{dust}} > T_{\text{gas}}$.</footnote>

<footnote n="7">⁷ We also include the following heating processes: photoheating of atoms and molecules due to ionisations and photodissociations (Grassi et al. 2014), including H2I UV pumping and H2I direct photodissociation heating (Burton, Hollenbach & Tielens 1990); heating due to the formation of H2I on dust and in the gas phase (Hollenbach & McKee 1979; Omukai 2000); and photoelectric heating by dust grains (Bakes & Tielens 1994). In this work, we do not include any ionization and heating from cosmic rays, as they become important only at very high densities (e.g. Papadopoulos & Thi 2013).</footnote>
For $T_{gas} < 10^4$ K, we either use these same tables (eq-metals runs) or non-equilibrium cooling by C I, C II, O I, O II, Si I, and Si II (non-eq-metals runs). We chose these elements because they are the most important metal coolants in the ISM (e.g. Wolfire et al. 2003). In this case, we solved a linear system for the individual metal excitation levels of the most important atoms and ions time dependently. In all cases, we impose a temperature floor of 10 K.

## 2.4 Evolving the galaxy

From the initial metallicity of the gas (0.1Z⊙, 0.5Z⊙, or Z⊙), we first compute $f_{H}$ and $f_{He}$ based on equation (4), to be consistent with what done by GASOLINE2 during the subsequent evolution. In the case of the non-eq-metals runs, we also impose the total abundances $f_{C} = 0.18f_{metals}$, $f_{O} = 0.43f_{metals}$, and $f_{Si} = 0.05f_{metals}$ (Asplund et al. 2009). We make an initial guess of the abundances of the 9 or 16 species, assuming a mostly neutral atomic gas. We then run the KROME code on all the gas particles in the initial conditions, for a time that depends on the density of each gas particle (shorter for denser gas), to evolve the system towards chemical equilibrium, keeping $T_{gas}$ constant.

With these new initial conditions, we run the code for 0.1 Gyr (roughly three times the local dynamical time of the galaxy, when using a typical $\rho_{gas} = 10^{-24}$ g cm$^{-3}$) using low values of the SF efficiency parameter (as explained in Section 2.2), during which the galaxy is relaxed. The ‘real’ runs start at the end of this ‘chemical equilibrium and relaxation’ process, when we reset the time to $t=0$, and last for 0.4 Gyr. We chose this final time as a compromise between the wish to have a truly in-equilibrium system (free of any numerical artefacts stemming from the initial conditions) and the notion that a system can be realistically considered in isolation for no more than a few $10^8$ yr, before cosmological gas inflows (Dekel, Sari & Ceverino 2009) and/or mergers (Genel et al. 2009) become important. We also note that, even if there is no gas replenishment from cosmological gas inflows, this problem is made less severe by the stellar winds from low-mass stars (Stinson et al. 2006).

In order to ensure the correct fragmentation behaviour of the gas, it is important that the Jeans mass ($M_{Jeans}$) and length of the gas are always resolved. In Capelo et al. (2015), where the mass resolution and SF density threshold were identical to those of this work, a temperature floor of 500 K was imposed, so that the Jeans mass of gas with $\rho_{gas} \lesssim \rho_{SF}$ would be resolved by at least 64 particles (see Capelo et al. 2017 for a more accurate explanation). In this work, we adopt instead a pressure floor, as described in Roškar et al. (2015; see also Agertz, Teyssier & Moore 2009). In GASOLINE2, the minimum gas pressure is set to $P_{gas,min} = \alpha G [\max(\epsilon_{gas}, h)]^2 \rho_{gas}^2$, where $\alpha$ is a parameter that sets the minimum number $N_{Jeans}$ of resolution elements, $\max(\epsilon_{gas}, h)$, that must resolve the Jeans length. If we define $L_{Jeans} = [\pi\gamma P_{gas}/(G\rho_{gas}^2)]^{1/2}$, where $\gamma = 5/3$ is the gas adiabatic index, imposing $L_{Jeans} > N_{Jeans} \max(\epsilon_{gas}, h)$ is equivalent to imposing $P_{gas} > [\rho_{gas}N_{Jeans} \max(\epsilon_{gas}, h)]^2 G/(\pi\gamma)$. Therefore, the number of resolution elements that resolve $L_{Jeans}$ is related to $\alpha$ by $N_{Jeans} = (\alpha\pi\gamma)^{1/2}$. We impose $\alpha = 3$ and obtain $N_{Jeans} = 4$ (see Truelove et al. 1997). A similar floor, within RAMSES, has already been tested in Gabor et al. (2016) when modelling a galaxy with similar resolution and SF density threshold as in this work. As a sanity check, we also impose $L_{Jeans} \ge \epsilon_{gas}$. We note that the presence of a pressure floor effectively makes the gas non-ideal when $P_{gas,min}$ is reached, since the gas temperature is allowed to cool below the ideal-gas equivalent $T_{gas,min}$ (see Fig. 1), as done in, e.g. Richings & Schaye (2016).

<table><caption>**Table 1.** Main simulations parameters. (1) Run name. (2) Metallicity of all the gas at the beginning of the relaxation. (3) Model of H2I formation on dust: CS09 or J75. (4) Clumping factor. (5) A chemical network of nine primordial and seven metal species is (is not) solved, together with non-equilibrium metal cooling, for the gas at low temperatures ($T < 10^4$ K) for the non-eq-metals (eq-metals) runs.</caption><thead><tr><th>Run</th><th>Z/Z<sub>☉</sub></th><th>H2I</th><th>C<sub>ρ</sub></th><th>Non-eq. metals</th></tr></thead><tbody><tr><td>01</td><td>0.5</td><td>CS09</td><td>1</td><td>No</td></tr><tr><td>01m</td><td>0.5</td><td>CS09</td><td>1</td><td>Yes</td></tr><tr><td>02</td><td>0.5</td><td>CS09</td><td>10</td><td>No</td></tr><tr><td>02m</td><td>0.5</td><td>CS09</td><td>10</td><td>Yes</td></tr><tr><td>03</td><td>0.5</td><td>J75</td><td>1</td><td>No</td></tr><tr><td>03m</td><td>0.5</td><td>J75</td><td>1</td><td>Yes</td></tr><tr><td>04</td><td>1.0</td><td>CS09</td><td>1</td><td>No</td></tr><tr><td>04m</td><td>1.0</td><td>CS09</td><td>1</td><td>Yes</td></tr><tr><td>05</td><td>0.1</td><td>CS09</td><td>1</td><td>No</td></tr><tr><td>05m</td><td>0.1</td><td>CS09</td><td>1</td><td>Yes</td></tr></tbody></table>

Metal injection from stellar evolution and turbulent metal diffusion imply that $f_{H}$, $f_{He}$, and $f_{metals}$ all evolve with time (see equation 4). At each time-step, before computing the chemical evolution, we rescale all species abundances according to

$$f_{\text{Hspecies, new}} = f_{\text{Hspecies, old}}(f_{\text{H, new}}/f_{\text{H, old}}), \quad (8)$$

where $f_{\text{Hspecies, new}}$ and $f_{\text{H, new}}$ are the new/old mass fraction of a hydrogen species (e.g. H I or H II) and of hydrogen, respectively. The same is done for He and, in the non-eq-metals runs, for metals. This way, the abundance ratios are conserved. Moreover, the mass fractions of C and Si are always approximately consistent with the metal ejection from SNa (we do not follow the evolution of Fe).

We emphasize here that this is the first work in which a metal network, its related non-equilibrium metal cooling, metal injection from SNa, and turbulent metal diffusion are all self-consistently modelled in a galaxy simulation. We remind the reader that SNa can change the metallicity of the gas by more than an order of magnitude in less than half Gyr (see Section 3.2).

### 2.4.1 The suite of simulations

The suite of simulations is composed of two distinct sub-sets, the eq-metals and the non-eq-metals set, which differ only in how metals are modelled (see Section 2.3). Within each sub-set, we vary (a) the method of H2I formation on dust (and include/exclude gas cooling by dust), (b) the clumping factor of such formation mechanism, and (c) the initial metallicity. This brings the total number of simulations to 10.

In Run 01, a chemical network of nine primordial species (H I, H II, H2I, H2 II, H$^-$, He I, He II, He III, and e$^-$) is followed and several photo and thermal processes are modelled during the evolution of an isolated gas-rich galaxy in which the initial metallicity of the gas is 50 per cent solar; H2I is formed on dust following the model of CS09, with no clumping factor (i.e. $C_ρ = 1$); and a shielded Haardt & Madau (2012) radiation background at $z=3$ is assumed.

Run 01m is the same as Run 01, except for additional seven metal species (C I, C II, O I, O II, Si I, Si II, and Si III) in the chemical network; and C I, C II, O I, O II, Si I, and Si II non-equilibrium cooling at low temperatures (i.e. $T_{gas} < 10^4$ K).

The other eight runs are variations of the first two runs, in which we vary one parameter at a time. The main parameters of the simulations are listed in Table 1.
# 3 RESULTS

In this section, we describe the results of our 10 simulations, mostly by means of comparing different sub-sets of runs. In Section 3.1, we assess the effects of including the chemical evolution and non-equilibrium cooling (for $T_{\text{gas}} < 10^4$ K) of seven metal species. In Section 3.2, we study the dependence on the (initial) gas metallicity. In Section 3.3, we evaluate the dependence of the results on how the formation of H2 on dust is modelled. In Section 3.4, we focus on the role of metals as tracers of the ISM and SFR. Unless otherwise stated, all results are shown at 0.4 Gyr, but we note that, being the system an isolated galaxy, all results hold also at earlier times ($t \gtrsim 0.15$ Gyr).

## 3.1 Equilibrium versus non-equilibrium metal cooling

In this section, we compare the results of two simulations, which are identical in everything except in the treatment of metal chemistry and cooling. We focus here on the comparison between Runs 01 and 01m, but we note that similar comparisons (i.e. Runs 02 and 02m, Runs 03 and 03m, Runs 04 and 04m, and Runs 05 and 05m) yield similar results, unless otherwise stated. We remind the reader that we use eq-metal – PIE – cooling in both runs for $T_{\text{gas}} \ge 10^4$ K, since we assume that gas at those temperatures is in equilibrium. Therefore, we expect to see differences only for $T_{\text{gas}} \lesssim 10^4$ K.

We summarize the main differences in the modelling of metals between Runs 01 and 01m. In addition to the obvious fact that (i) metal cooling is computed assuming equilibrium in Run 01 and lack thereof in Run 01m, we note that (ii) the network utilized for the PIE tables contains the first 30 elements in the periodic table, whereas the non-eq-metal computations include seven metal species. Additionally, we stress that (iii) the PIE tables were constructed assuming an optically thin gas, whereas the non-equilibrium metal cooling explicitly assumes a shielded UV background.

Assuming an optically thin or thick gas has potentially considerable implications. As a reference, at $T_{\text{gas}} < 10^4$ K, the PIE (UV at $z=3$; Haardt & Madau 2012) cooling is two orders of magnitude stronger than the CIE (no UV) cooling (Bovino et al. 2016; see also Shen et al. 2010). Therefore, the PIE tables likely overestimate the metal cooling at low temperatures.

The different number of species can also in principle produce different values of metal cooling. For example, Richings et al. (2014a; see also Richings et al. 2014b) claim that Fe, which we have not included in our network, could be potentially relevant. On the other hand, Fe is believed to be strongly depleted into solid form (Jenkins 2009), and there is general agreement that C, O, and Si (the metal elements we do follow in our simulations) are by far the dominant metal coolants in the ISM (e.g. Tielens & Hollenbach 1985; Wolfire et al. 1995, 2003; Shen et al. 2010).

Using the same networks employed in this work, Bovino et al. (2016) performed isochoric one-zone tests with constant gas temperature to compute the ratio of equilibrium to non-equilibrium metal cooling (in both cases in the optically thin approximation; see their fig. 16), showing that gas at both low ($0.1 \text{ cm}^{-3}$) and high ($10^2 \text{ cm}^{-3}$) densities experiences two evolution regimes: during the first part of the evolution, the non-equilibrium metal cooling is stronger than the PIE cooling (even by orders of magnitude), because of the larger number of electrons (due to incomplete recombinations) and consequent enhanced collisions; during the second part, the system evolves towards an equilibrium state and the differences between equilibrium and non-equilibrium metal cooling tend to diminish, with the ratio between the eq-metal to non-eq-metal

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

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

**Figure 1.** Dependence on the modelling of metals. Gas $\rho$–$T$ diagrams at 0.4 Gyr, for the eq-metals Run 01 (upper panel) and the non-eq-metals Run 01m (lower panel). The horizontal and vertical lines show the SF temperature and density thresholds, respectively. The diagonal lines show the temperature corresponding to the pressure floor, assuming an ideal gas and $\mu = 0.60$ (lower line) or 2.30 (upper line). For all the gas phase diagrams in this work, we chose the limits of the colour bars such that a dark blue pixel means there is no gas in that phase-space region, since the (initial) mass of a gas particle is $4.6 \times 10^3 \text{ M}_\odot$. The number of logarithmic bins in the x- and y-axis is always $10^2$. The eq-metals gas is on average colder than the eq-metals gas.

cooling rates (both at equilibrium) varying between ~0.6 (for low $T_{\text{gas}}$) and ~3 (for $T_{\text{gas}} = 9 \times 10^3$ K).

In Fig. 1, we show the gas density–temperature diagrams for the two runs. In both diagrams, one can observe three distinct regions in the phase space: (1) the gas with $T_{\text{gas}} \gtrsim 10^4$ K and $10^{-25} \lesssim \rho_{\text{gas}} \lesssim 10^{-22} \text{ g cm}^{-3}$ is the ISM heated by SNae;<fnref n="8" /> (2) this same gas eventually cools down and fills a region of the phase diagram at $\rho_{\text{gas}} \lesssim 10^{-25} \text{ g cm}^{-3}$ and $T_{\text{gas}} \sim 10^4$ K (the so-called warm neutral medium, WNM; e.g. Ferrière 2001), under which temperature the

<footnote n="8">⁸ We note that this relatively dense, hot gas is kept artificially hot for some time because of the cooling shut-off SN model (see Section 2.2 and Stinson et al. 2006 for more details), hence it only shares the gas temperature with the observed hot ionized medium (Ferrière 2001).</footnote>
magnitude of the cooling function drops dramatically; (3) part of this gas, together with the gas that was never heated by SNaе, manages to cool down to $T_{gas} < 10^4$ K, with the shape of this third region heavily depending on the low-temperature cooling function.

The phase diagrams exhibit differences according to the above expectations. The eq-metals gas is on average colder than the non-eq-metals gas at all densities, with the difference in temperature increasing as the density decreases, ranging from a factor of a few at high densities to ~one order of magnitude at low densities. In particular, the gas mass fraction of WNM (defined here as the gas with $3.8 < \log_{10}(T_{gas}[K]) < 4.2$ and $\rho_{gas} < 10^{-24} \text{ g cm}^{-3}$) is much larger in the non-eq-metals case than in the eq-metals case (16 versus 2 per cent). Most importantly, the gas mass fraction of gas with $T_{gas} < 10^2$ K (the upper temperature limit for the cold neutral medium; e.g. Ferrière 2001) is 1.5 and 76.8 per cent for the non-eq-metals and eq-metals gas, respectively.

Because of the substantial differences between CIE and PIE cooling for $T_{gas} < 10^4$ K observed in, e.g. Shen et al. (2010) and Bovino et al. (2016), we attribute the differences between Runs 01 and 01m mostly to the optically thin/thick approximation in the metal cooling. To assert this, we ran a version of Runs 01 and 01m in which we assumed the UV background to not be shielded and found that the gas thermodynamics of the two no-shielding runs is very similar (see Appendix A), with a couple of minor differences: (a) in the density range $\rho_{gas} \lesssim 10^{-23} \text{ g cm}^{-3}$, the slope is slightly different (constant temperature versus marginally decreasing temperature); (b) in the density range $10^{-23} \lesssim \rho_{gas} \lesssim 10^{-22} \text{ g cm}^{-3}$, the eq-metals gas is slightly colder than the non-eq-metals gas. The differences in the two gas $\rho-T$ diagrams of the original Runs 01 and 01m are therefore likely due to the implementation of shielding (or lack thereof), rather than to the differences in the modelling of metals.

In order to understand the causes of the remaining, minor differences and disentangle non-equilibrium effects from those due to different metal networks, we once again consider the isochoric tests at constant temperature performed by Bovino et al. (2016). We note that the non-eq-metals models reach equilibrium after ~$10^3-10^5$ yr, depending on the density. These time-scales are much shorter than the dynamical time-scales at those densities, implying that, in our simulations of an isolated galaxy, the gas can be considered somewhat in equilibrium. We caution, however, that the comparison in Bovino et al. (2016) is on metal cooling rather than on net cooling. Therefore, we conservatively state that the minor differences that are not attributable to shielding are due to a combination of different metal networks and non-equilibrium effects.

The different $T_{gas}$ has direct consequences on the chemical composition of the gas. In Fig. 2, we show the dependence of the mean molecular weight on gas density. In both runs, the vast majority of the gas lies around (and in between) the two phase-space regions of neutral molecular ($\mu \sim 2.30$) and neutral atomic ($\mu \sim 1.24$) gas, with only a small fraction as ionized atomic ($\mu \sim 0.60$) gas. The difference between the two runs lies predominantly in the neutral molecular region: there is more H2I in the non-eq-metals gas (23 per cent of the gas mass has $\mu > 2.26$) than in the eq-metals gas (13 per cent). This is because formation of H2I on dust, according to the CS09 model, peaks at $\log_{10}(T_{gas}) \sim 2.2$ (owing to the opposing dependences on $T_{gas}$ of dust sticking and gas thermal velocity), and there is much more non-eq-metals gas at or around the peak temperature (36 per cent of the gas mass has $2 < \log_{10}(T_{gas}[K]) < 2.4$) than eq-metals gas (6 per cent). This interpretation is corroborated by the fact that the gas mass fractions of gas with $\mu > 2.26$ for the two runs using the J75 model of H2I formation on dust (in which there is no dependence on gas temperature) are closer to each other: 8

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Dependence on the modelling of metals. Gas $\rho-\mu$ diagrams at 0.4 Gyr, for the eq-metals Run 01 (upper panel) and the non-eq-metals Run 01m (lower panel). The horizontal lines show the values of $\mu$ for (from top to bottom) a fully neutral molecular gas ($f_{H} = f_{H2I}$, $f_{He} = f_{HeI}$), a fully neutral atomic gas ($f_{H} = f_{H1}$, $f_{He} = f_{HeI}$), and a fully ionized atomic gas ($f_{H} = f_{H1}$, $f_{He} = f_{HeII}$), computed for a gas with $Z = 0.5Z_{\odot}$. The mean molecular weight can get higher than the upper horizontal line because of the range in metallicity (see the upper middle panel of Fig. 8 for Run 01; the equivalent plot for Run 01m is very similar), since $\mu$ increases with Z. The non-eq-metals gas is more molecular than the eq-metals gas.</figcaption></figure>

and 13 per cent for Runs 03 and 03m, respectively. Moreover, even though they are less important than catalysis on dust, we note that the gas-phase reactions of H2I formation peak at $10^4-10^5$ K, where the non-eq-metals run has much more gas than in the eq-metals run (see also the difference in WNM gas mass stated above).

A colder gas is also prone to more fragmentation (Toomre 1964). In Fig. 3, we show the surface density maps, viewed face-on, for the total gas. The central ($\lesssim 5$ kpc), denser regions of the galaxy are very similar in the two runs. On the other hand, the outer regions – where the gas density is lower and therefore the temperature difference between the two runs is larger (see Fig. 1) – are quite different. The eq-metals gas map shows a clumpy outer region, whereas the outer gas in the non-eq-metals run is very smooth. The H2I maps in Fig. 4 show a slightly more extended disc in the eq-metals case. This is linked to the enhanced fragmentation in Run
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Dependence on the modelling of metals. Gas surface density maps at 0.4 Gyr for the eq-metals Run 01 (upper panel) and the non-eq-metals Run 01m (lower panel). The eq-metals gas is more fragmented than the non-eq-metals gas.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Dependence on the modelling of metals. H2I surface density maps at 0.4 Gyr for the eq-metals Run 01 (upper panel) and the non-eq-metals Run 01m (lower panel). The eq-metals molecular disc is slightly more extended than the non-eq-metals counterpart.</figcaption></figure>

01, where the local H2I fraction is increased within clumps. We note, however, that, despite this, the total amount of H2I is larger in the non-eq-metals run.

The large difference in the density-temperature diagrams does not translate into large differences in SF, as it would be naively expected (since SF occurs in cold dense gas). In Fig. 5, we show the SFR as a function of time for Runs 01 and 01m. In both simulations, the SFR varies around 0.7–0.9 M$_{\odot}$ yr$^{-1}$, with the non-eq-metals run forming slightly more stars than the eq-metals run. The difference, however, is minimal and approaches zero in the last ~100 Myr of the simulation, with the final stellar mass equal to $2.7 \times 10^9$ M$_{\odot}$ in both runs.

In Figs 6 and 7, we show the classic and molecular KS diagrams, respectively. To construct these diagrams, we positioned the galaxy face-on, divided it into 50 concentric cylindrical bins of thickness 200 pc between 0 and 10 kpc, and computed the surface density of the H1+H2I gas (for the classic KS diagram) or H2I gas (for the molecular KS diagram), and of the SFR, defined here as the mass of new stars formed in the past 100 Myr, divided by 100 Myr. In

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Dependence on the modelling of metals. SFR versus time for the eq-metals Run 01 (black solid line) and the non-eq-metals Run 01m (red dashed line). The two runs exhibit very similar SF histories.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Dependence on the modelling of metals. KS diagrams at 0.4 Gyr between H1+H2I gas and SFR (in the last 100 Myr) surface density, for the 50 innermost 200-pc-thick concentric cylindrical bins, for the eq-metals Run 01 (upper panel) and the non-eq-metals Run 01m (lower panel). The data from the simulation (white points) are compared to the contours by Bigiel et al. (2010). The two runs display very similar classic KS diagrams and both agree very well with the observations.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Dependence on the modelling of metals. Molecular KS diagrams at 0.4 Gyr between H2I gas and SFR (in the last 100 Myr) surface density, for the 50 innermost 200-pc-thick concentric cylindrical bins, for the eq-metals Run 01 (upper panel) and the non-eq-metals Run 01m (lower panel). The data from the simulation (white points) are compared to the contours by Bigiel et al. (2010). The two runs feature very similar molecular KS diagrams.</figcaption></figure>

all the diagrams, we overplotted the data by Bigiel et al. (2010).<fnref n="9" /> We note that our simulations agree very well with the observed classic KS relation, both in slope and in normalization (the only discrepancy occurring at high gas surface densities, due to small number statistics in the observed data). This was expected, since the SF recipe and SF efficiency parameter were constructed and tuned to match observations.

Our data also agree reasonably well with the normalization of the observed molecular KS relation – in the sense that the simulated and observed curves intersect in the middle of the gas surface density range ($10 \lesssim \Sigma_{H2I} \lesssim 30 M_{\odot} pc^{-2}$) observed by Bigiel et al. (2010) – although they have a steeper slope than the observed one. As noted in the next sections, however, the slope discrepancy becomes less severe when we decrease the gas metallicity (see Section 3.2).

or change the method of H2I formation on dust (see Section 3.3). We stress here that, in our SF recipe, we did not directly link the abundance of molecular hydrogen to SF, as done in other works (e.g. Gnedin et al. 2009; Christensen et al. 2012; Tomassetti et al. 2015; Pallottini et al. 2017b). The correlation between H2I and SFR surface densities emerges naturally from the simulations, hinting to the idea that SF occurs in cold, dense gas, where molecular hydrogen is generally more abundant (see also discussion in Lupi et al. 2018).

The reason why the SFR and KS diagrams are not that different despite the large differences in the gas density–temperature diagrams is because those differences mostly occur at relatively low densities. The gas structure in the SF region (defined by $T_{gas} < T_{SF}$ and $\rho_{gas} > \rho_{SF}$) is almost the same, as the gas mass fraction of potentially star-forming gas in Runs 01 and 01m is 7 and 8 per cent, respectively. Moreover, the difference in H2I abundance is not large enough to produce discernible discrepancies in the molecular KS diagrams. All of the above suggests that our model is relatively robust, when it comes to SF observables, in spite of some differences in the

<footnote n="9">⁹ From their online data, we combined the Spiral2 sample of table 2 (applying the filter $log_{10}(\Sigma_*[M_{\odot}yr^{-1}kpc^{-2}]) > -3.7$) and the Spirals sample of table 3.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Dependence on the initial gas metallicity. Gas phase diagrams at 0.4 Gyr for eq-metals runs of different initial metallicity [Run 05 (0.1Z⊙; left-hand panels), Run 01 (0.5Z⊙; central panels), and Run 04 (Z⊙; right-hand panels)]: metallicity (upper panels), temperature (middle panels), and H2I hydrogen mass fraction (lower panels), all as a function of gas density. The spread of metallicity and the gas ρ–T structure are similar amongst the three runs, whereas the amount of H2I increases with $f_{\text{metals}}$.</figcaption></figure>

gas thermodynamics. Indeed, the mass of potentially star-forming gas is very similar for a large range of SF density and temperature thresholds, down to $T_{\text{SF}} \sim 10^2$ K, when the gas thermodynamic differences become important and the eq-metals simulation forms more stars than the non-eq-metals run.

## 3.2 Dependence on metallicity

In this section, we study the dependence on the (initial) metallicity of the gas, by comparing Runs 01, 04, and 05. We remind the reader that, due to continuous metal injection from SNae and turbulent metal diffusion, the metallicity of the gas varies significantly with time, as it can be seen in the upper panels of Fig. 8. We note that the magnitude of the spread in metallicity is basically the same in all runs: the amount of metals deposited into the gas from SNae is more or less the same, regardless of the initial metallicity. This is because the amount of SF (and subsequent SNae) is very similar, yielding the same increase in average $Z/Z_⊙$ (equal to 0.3) for all three runs. For this reason, the final metallicity simply depends on the initial metallicity, which was varied to test its effects on the gas thermodynamics and abundances and not to reproduce the observed mass–metallicity relation. Nevertheless, given the wide

range of possible values of gas metallicity for a stellar mass of $2.7 \times 10^9 M_⊙ (Z/Z_⊙ \sim 0.2–1.3$; see e.g. Sánchez et al. 2017), all our galaxies are within the observed range.

In our model, there are several physical phenomena linked to the amount of metals in the gas. Metal cooling increases with metallicity, because of the enhanced density of colliders (ions, electrons), and this increase is commonly modelled in PIE runs with the metal cooling scaling linearly with $Z/Z_⊙$. A higher metallicity also implies a higher dust-to-gas ratio, which in this work is assumed to scale as $D = D_⊙ Z/Z_⊙$. Consequently, we also expect stronger gas cooling by dust, a larger abundance of molecular hydrogen, and stronger molecular-hydrogen cooling (which, however, at the densities and metallicities considered in this work, is sub-dominant with respect to metal cooling).

In agreement with the above expectations, the gas density–temperature diagrams (middle panels of Fig. 8) show that the gas gets colder as the metallicity increases. In particular, when comparing Runs 05 (of initial metallicity 0.1Z⊙) and 04 (Z⊙), we observe relatively more gas in the low-metallicity run at $\sim 10^4$ K (at low densities) and in the region of phase space that connects the $10^4$-K gas with the $10^2$-K gas (at $\sim 10^{-25}$ g cm$^{-3}$). Moreover, if we select gas at $10^{-26}$ g cm$^{-3}$, the coldest gas in Runs 05 and 04 has a
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Dependence on the initial gas metallicity. Face-on surface density maps of the H2I gas (upper panels) and KS diagrams for the H2I gas (middle panels) and H1+H2I gas (lower panels) at 0.4 Gyr, for eq-metals runs of different initial metallicity: Run 05 (0.1Z⊙; left-hand panels), Run 01 (0.5Z⊙; central panels), and Run 04 (Z⊙; right-hand panels). The decrease of H2I with decreasing metallicity implies a less steep slope in the molecular KS diagram, which better matches observations.</figcaption></figure>

temperature of ~10⁴ and 300 K, respectively. The dependence on $f_{\text{metals}}$, however, does not appear to be very strong: this is due to the spread in metallicity, which reduces the effective difference between runs (e.g. from the initial factor of 10 to only a few, when comparing Runs 05 and 04). This is the likely explanation why the effect is visibly stronger in Richings & Schaye (2016), who do not include any chemical enrichment from SNae (nor turbulent metal diffusion) and therefore have a constant metallicity throughout their simulations.

As expected from the direct link between metals, dust, and H2I formation, the amount of molecular gas increases with metallicity, as it can be seen in the lower panels of Fig. 8. The gas mass fraction of gas with $f_{\text{H2I}}/f_{\text{H}} > 0.9$ for Runs 05, 04, and 01 is 12, 22, and 30 per cent, respectively. We note that the density at which the gas starts transitioning from atomic to molecular is more or less the same in the three runs (at ~10⁻²⁴ g cm⁻³), in disagreement with previous studies (both assuming equilibrium – McKee & Krumholz 2010; Sternberg et al. 2014 – and not – Gnedin et al. 2009). This is caused by the fact that metal injection and diffusion smear out the differences in metallicity between the runs (e.g. the ratio between the average metallicities of Runs 04 and 05 starts at 10 and ends at 3). We also note a wide spread of densities in the region of gas transitioning from atomic to molecular.

In Fig. 9, we show gas maps and KS diagrams for the same runs. We note that there is almost no difference in the total gas surface density (not shown), but there is a clear increase of H2I with metallicity (upper panels; consistent with the lower panels of Fig. 8). This difference in H2I abundance translates into a visible change in the molecular KS law (middle panels): because the difference is larger at low densities (as at higher densities the gas saturates to fully molecular), the molecular KS law in Run 05 has a slope that matches quite well the data by Bigiel et al. (2010), even though the observed galaxies, being local, have a slightly higher metallicity (we remind the reader that the final average gas metallicity of Run 05 is Z/Z⊙ = 0.4). The classic KS diagrams (lower panels) do not show any difference, as expected.

When performing the same kind of comparison with the non-eq-metals runs (i.e. Runs 01m, 04m, and 05m), we obtain similar results.

### 3.3 H2I formation on dust: models and clumping

In this section, we briefly study the effects of changing the model of formation of H2I on dust (J75 versus CS09) and the value of the clumping factor ($C_\rho = 1$ or 10). We chose the value 10 because it is the same value used in Gnedin et al. (2009), Christensen et al.
<figure><img src="image_12.png" /><figcaption>Figure 10. Dependence on the model of H2 I formation on dust. We compare the H2 I hydrogen mass fraction for Run 01 (CS09 and $C_\rho = 1$; left-hand panel), Run 02 (CS09 and $C_\rho = 10$; central panel), and Run 03 (J75 and $C_\rho = 1$; right-hand panel) at 0.4 Gyr.</figcaption></figure>

(2012), and Tomassetti et al. (2015). Moreover, it is also consistent with the values found in the simulations with a variable clumping factor by Lupi et al. (2018).

The amount of H2 I varies significantly when we vary the formation model and/or the clumping factor (see Fig. 10). As expected, assuming a clumping factor of 10 increases the amount of H2 I: the gas mass fraction of gas with $f_{H2I} / f_H > 0.9$ for Runs 01 ($C_\rho = 1$) and 02 ($C_\rho = 10$) is 22 and 50 per cent, respectively. The choice of formation model is also important, with the J75 model (Run 03) yielding less H2 I than the CS09 model (15 versus 22 per cent), but not as crucial as the choice of $C_\rho$. This is because our simulations are all run at relatively high metallicity, where the two models of H2 I formation on dust are most similar. The differences are visible also in the H2 I surface density maps (not shown), even though they are quite clear only when comparing different clumping factors.

The different H2 I abundance does not have a large effect on the density-temperature structure of the gas (not shown). This is due to the fact that the most important H2 I-related phenomenon that could in principle affect the gas thermodynamics, molecular-hydrogen cooling, is negligible with respect to metal cooling, at the densities and metallicities considered in this work.

The lack of significant variations in the gas $\rho$-T diagram explains the absence of discrepancy in the classic KS diagrams (not shown). On the other hand, because of the disparity in H2 I abundances, the molecular KS diagrams are slightly different, with Runs 01 and 03 matching the observed slope marginally better than Run 02. In our simulations, we do not need a (constant) clumping factor to improve the match with the molecular KS relation. In fact, the simulations with $C_\rho = 1$ perform slightly better. We caution, however, that this is also the result of our choice of parameters (including the SF efficiency parameter).

We performed the same kind of comparison with the non-eq-metals runs (i.e. Runs 01m, 02m, and 03m) and obtained the same results.

## 3.4 Metal species as tracers of star formation and the interstellar medium

In this section, we study the role of metal species as tracers of the ISM and SFR, focusing on the non-eq-metals Run 01m.

In Fig. 11, we show the abundances of C I, C II, O I, and Si II as a function of gas density. C I and C II are somewhat complementary: C I is more abundant for $\rho_{gas} \gtrsim 10^{-23}$ g cm$^{-3}$, whereas most of C II is found at lower densities. O I and Si II, on the other hand, have a constant abundance at almost all densities. The same phase diagrams for O II, Si I, and Si III (not shown) reveal instead that these species are negligible. In Fig. 12, the gas surface density maps reinforce the

notion that C I traces higher-density gas (in the central regions of the galaxy). C II, on the other hand, traces the entire ISM, as do O I and Si II (not shown).

Metal lines have also been used as tracers of SFR (e.g. Stacey et al. 1991). De Looze et al. (2014) analysed the applicability of metal lines as tracers of the ISM in different types of galaxies, obtaining the general relation

$$ \log_{10}(\dot{\Sigma}_*) = \xi + \chi \log_{10}(\Sigma_{\text{line}}), \quad (9) $$

where $\dot{\Sigma}_*$ is given in $M_{\odot} \text{yr}^{-1} \text{kpc}^{-2}$, $\Sigma_{\text{line}}$ (with line = C II or O I) is given in $L_{\odot} \text{kpc}^{-2}$, and [$\chi, \xi$] = [1.01 ± 0.02, -6.99 ± 0.14] and [1.00 ± 0.03, -6.79 ± 0.22] for C II and O I, respectively. For simplicity, hereafter we assume $\chi = 1$ and $\xi = -7$.

In order to compare our results to observations, we first translate the metal-species gas mass surface densities from our simulations into metal-line surface brightnesses. To do so, we follow Pallottini et al. (2017b), who find (for $n_{gas} \lesssim 10^3 \text{ cm}^{-3}$, the critical density for C II emission)

$$ \Sigma_{\text{CII,lum}} = \epsilon_{\text{CII}} \Sigma_{\text{gas}}, \quad (10) $$

where the emissivity

$$ \epsilon_{\text{CII}} = 0.1 \left( \frac{n_{\text{gas}}}{100[\text{cm}^{-3}]} \right) \frac{Z}{Z_{\odot}} [L_{\odot}/M_{\odot}]. \quad (11) $$

Combining equations (9)–(11), considering typical values from Run 01m [$\mu = 2, Z/Z_{\odot} = 0.5$, and $f_{CII} = 10^{-3}$ (see the upper right-hand panel of Fig. 11)], and assuming for simplicity that $\Sigma_{\text{CII, mass}} = 2\rho_{CII}z_{\text{d}}$, we obtain

$$ \log_{10}(\dot{\Sigma}_*) = 0.6 + 2\log_{10}(\Sigma_{\text{CII, mass}}), \quad (12) $$

where $\dot{\Sigma}_*$ is given in $M_{\odot} \text{yr}^{-1} \text{kpc}^{-2}$ and $\Sigma_{\text{CII}}$ is given in $M_{\odot} \text{pc}^{-2}$. Equation (12) is plotted in the left-hand diagram of Fig. 13, where we show the relation between SFR and C II surface density for the non-eq-metals Run 01m.

The agreement between our data and observations is very good for $\Sigma_{\text{CII, mass}} \lesssim 0.04 M_{\odot} \text{pc}^{-2}$. The mismatch at high surface densities is due to the drop of $f_{CII}$ at high gas densities (where recombination is more efficient), as shown in the upper right-hand panel of Fig. 11.

Since we do not have an equivalent of equations (10) and (11) for O I, we made the simple assumption that a relation with the same slope of equation (12) holds. In order to match the simulation data, we then divided the normalization factor by 10. This is equivalent to assuming an O I emissivity 10 times smaller than that of C II, which is consistent with observations (e.g. Malhotra et al. 2001; Cormier et al. 2015; Michałowski et al. 2016).
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Metals as tracers of the ISM. Mass fractions at 0.4 Gyr, all as a function of gas density, for the non-eq-metals Run 01m: C I (upper left-hand panel), C II (upper right-hand panel), O I (lower left-hand panel), and Si II (lower right-hand panel). O I and Si II trace the entire ISM, whereas C I and C II trace different ISM regimes.</figcaption></figure>

## 4 CONCLUSIONS

We used a novel interface between the modern-SPH code GASOLINE2 and the chemistry package KROME, to accurately follow the hydrodynamical and chemical evolution of an isolated galaxy at z = 3. This is the first work in which metal injection from SNae, turbulent metal diffusion, and a metal network with its corresponding non-equilibrium metal cooling were self-consistently modelled in a galaxy simulation.

We built a suite of 10 high-resolution hydrodynamical simulations, in which we varied the chemical network, how metal cooling is modelled, the initial gas metallicity, and how H2I forms on dust, to assess the importance of each physical parameter/prescription and, in the case of the non-eq-metals runs, to study the role of metals as tracers of the ISM and SFR.

We itemize our findings below:

(i) Our default run with half-solar metallicity gas, a primordial network (H I, H II, H−, H2 I, H2 II, He I, He II, He III, and e−) with dust, optically thin equilibrium (PIE) metal cooling, and the CS09 model of H2I formation on dust with no clumping factor is able to create a realistic galaxy that follows the classic KS relation (Fig. 6). Indeed, all the simulations in this work recover the classic KS relation. The simulated molecular KS relation has a slightly steeper slope than

the observed one, but the simulated and observed curves intersect in the middle of the observed gas surface density range (Fig. 7).

(ii) Modelling the chemical evolution of seven metal species (C I, C II, O I, O II, Si I, Si II, and Si III) and optically thick non-equilibrium metal cooling has a substantial effect on the thermodynamics of the gas, the chemical abundances, and the appearance of the galaxy. With respect to the eq-metals gas, non-eq-metals gas is typically warmer (Fig. 1), has a larger molecular-gas mass fraction (Fig. 2), a smoother disc (Fig. 3), and a slightly less extended molecular disc (Fig. 4). The differences are due to the proper account of shielding when computing metal cooling and, to a lesser extent, to different metal networks and non-equilibrium effects.

(iii) However, the differences in the gas thermodynamics are typically small for cold, dense gas, where stars form, implying that both the SF history (Fig. 5) and KS diagrams (Figs 6 and 7) are similar for eq-metals and non-eq-metals simulations.

(iv) Varying the initial gas metallicity has direct consequences on gas cooling and H2I formation, with a higher metallicity producing a colder gas and a larger fraction of molecular gas (Fig. 8). The differences, however, are less conspicuous than previously thought because metal enrichment from SNae and turbulent metal diffusion dilute the initial metallicity discrepancy between simulations. The low-metallicity run recovers not only the classic but also the molecular KS relation.
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Metals as tracers of the ISM. Surface density maps at 0.4 Gyr for the non-eq-metals Run 01m: face-on C<sub>I</sub> (upper left-hand panel), face-on C<sub>II</sub> (upper right-hand panel), edge-on C<sub>I</sub> (lower left-hand panel), and edge-on C<sub>II</sub> (lower right-hand panel).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 13.** KS-like diagrams at 0.4 Gyr, for the non-eq-metals Run 01m, between C<sub>II</sub> (left-hand panel) or O<sub>II</sub> (right-hand panel) and SFR (in the last 100 Myr) surface densities, for the 50 innermost 200-pc-thick concentric cylindrical bins (black points). The solid line in the left-hand panel shows the relation from equation (12), whereas that in the right-hand panel shows the same relation with Σ<sub>CII</sub> replaced by Σ<sub>OII</sub> and with the normalization divided by 10.</figcaption></figure>

(v) Varying the model of H2 I formation on dust and the clumping factor has a direct effect on the abundance of molecular hydrogen (Fig. 10). Expectedly, a larger clumping factor yields more H2 I. We also find that using the J75 model results in a smaller fraction of H2 I than when using the CS09 model, hence producing a molecular

KS diagram with a slope slightly closer to the observed one. Using a clumping factor of 10 moves the simulated molecular KS relation away from the observed one.

(vi) Properly following the chemical evolution of metal species is not only useful to account for the non-equilibrium metal cooling:
the metal species themselves can be a valuable tracer for the ISM (Figs 11 and 12) and SFR (Fig. 13). We showed that our simulations agree quite well with observations that link the SFR to C II emission lines.

In summary, an improved and more accurate prescription of opti-
cally thick non-equilibrium metal cooling, together with the possi-
bility to model and follow individual metal species (in the context of
comparisons with observations), are the main qualities of this new
interface. We note here that improvements can and will be made in
future work. On the hydrodynamics code side, the most important is
the addition of local radiation from young stars, and we refer to Lupi
et al. (2018) for a thorough study performed on a galaxy identical to
ours. On the chemistry side, we wish to add more metal and molec-
ular species, to improve our modelling of non-equilibrium metal
cooling and be able to follow more metal and molecular lines as
tracers of the ISM and SFR. Additionally, we plan on constructing
synthetic emission maps from selected lines (e.g. the [C II] 157.74
µm line) to provide a better comparison with current and future FIR
observations.

The work presented here was the *first* application of an interface between KROME and an SPH implementation of the equations of hydrodynamics, and we expect this interface to be used in several applications, e.g. from isolated simulations of galaxies, mergers of galaxies, and circumnuclear discs, to cosmological zoom-in simulations of galaxy and BH formation.

ACKNOWLEDGEMENTS

We thank the reviewer for the useful comments that greatly im-
proved this work. P.R.C. thanks Robert Feldmann, Francesco
Haardt, Alexander P. Hobbs, Andrea Pallottini, Rafael Souza Lima,
and James W. Wadsley for useful exchanges, and Sijing Shen
for fruitful discussions and for having provided the metal cool-
ing table described in Shen et al. (2010, 2013). Numerical cal-
culations were performed on Piz Daint, at the Swiss National
Supercomputing Center (CSCS). Part of the analysis was per-
formed using the PYNBODY package (Pontzen et al. 2013). P.R.C.
acknowledges research funding by the Deutsche Forschungsge-
meinschaft (DFG) Sonderforschungsbereiche (SFB) 963 and sup-
port by the Tomalla Foundation. A.L. acknowledges support from
the European Research Council (Projects No. 267117, ‘DARK’,
and No. 614199, ‘BLACK’). S.B. and D.R.G.S. thank for funding
through the DFG priority programme ‘The Physics of the Interstel-
lar Medium’ (projects BO 4113/1-2 and SCHL 1964/1-2). D.R.G.S.
acknowledges funding via Fondecyt regular (project 1161247),
the ‘Concurso Proyectos Internacionales de Investigación, Con-
vocatoria 2015’ (project PII20150171), ALMA–Conicyt (project
31160001), and the BASAL Centro de Astrofísica y Tecnologías
Afines (CATA) PFB-06/2007. T.G. acknowledges the Centre for
Star and Planet Formation funded by the Danish National Research
Foundation.

REFERENCES

Abel T., Anninos P., Zhang Y., Norman M. L., 1997, New A, 2, 181

Adelberger K. L., Steidel C. C., Pettini M., Shapley A. E., Reddy N. A., Erb
D. K., 2005, ApJ, 619, 697

Agertz O. et al., 2007, MNRAS, 380, 963

Agertz O., Teyssier R., Moore B., 2009, MNRAS, 397, L64

Aldrovandi S. M. V., Pequignot D., 1973, A&A, 25, 137

ALMA Partnership et al., 2015, ApJ, 808, L1

Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A, 47, 481

Bakes E. L. O., Tielens A. G. G. M., 1994, ApJ, 427, 822

Barlow S. E., 1984, PhD thesis, Univ. Colorado Boulder

Bellovary J. M., Governato F., Quinn T. R., Wadsley J., Shen S., Volonteri
M., 2010, ApJ, 721, L148

Bellovary J., Brooks A., Volonteri M., Governato F., Quinn T., Wadsley J., 2013, ApJ, 779, 136

Bigiel F., Leroy A., Walter F., Brinks E., de Blok W. J. G., Madore B., Thornley M. D., 2008, AJ, 136, 2846

Bigiel F., Leroy A., Walter F., Blitz L., Brinks E., de Blok W. J. G., Madore B., 2010, AJ, 140, 1194

Bovino S., Grassi T., Schleicher D. R. G., Latif M. A., 2014, ApJ, 790, L35

Bovino S., Grassi T., Capelo P. R., Schleicher D. R. G., Banerjee R., 2016, A&A, 590, A15

Braun H., Schmidt W., 2015, MNRAS, 454, 1545

Brook C. B. et al., 2012, MNRAS, 426, 690

Bryan G. L. et al., 2014, ApJS, 211, 19

Burton M. G., Hollenbach D. J., Tielens A. G. G. M., 1990, ApJ, 365, 620

Camm G. L., 1950, MNRAS, 110, 305

Capelo P. R., Volonteri M., Dotti M., Bellovary J. M., Mayer L., Governato
F., 2015, MNRAS, 447, 2123

Capelo P. R., Dotti M., Volonteri M., Mayer L., Bellovary J. M., Shen S., 2017, MNRAS, 469, 4437

Capitelli M., Coppola C. M., Diomede P., Longo S., 2007, A&A, 470, 811

Cazaux S., Spaans M., 2009, A&A, 496, 365

Cen R., 1992, ApJS, 78, 341

Christensen C. R., Quinn T., Stinson G., Bellovary J., Wadsley J., 2010, ApJ, 717, 121

Christensen C., Quinn T., Governato F., Stilp A., Shen S., Wadsley J., 2012, MNRAS, 425, 3058

Conroy C., Gunn J. E., White M., 2009, ApJ, 699, 486

Coppola C. M., Longo S., Capitelli M., Palla F., Galli D., 2011, ApJS, 193, 7

Cormier D. et al., 2015, A&A, 578, A53

Corrigan S. J. B., 1965, J. Chem. Phys., 43, 4381

Dalgarno A., Lepp S., 1987, in Vardya M. S., Tarafdar S. P., eds, Proc. IAU Symp. Vol. 120, Astrochemistry, Reidel, Dordrecht, p. 109

Dalla Vecchia C., Schaye J., 2008, MNRAS, 387, 1431

Dalla Vecchia C., Schaye J., 2012, MNRAS, 426, 140

de Jong T., 1972, A&A, 20, 263

De Looze I. et al., 2014, A&A, 568, A62

Dehnen W., Aly H., 2012, MNRAS, 425, 1068

Dekel A., Sari R., Ceverino D., 2009, ApJ, 703, 785

Diemer B., Kravtsov A. V., 2015, ApJ, 799, 108

Dove J. E., Rusk A. C. M., Cribb P. H., Martin P. G., 1987, ApJ, 318, 379

Draine B. T., Lee H. M., 1984, ApJ, 285, 89

Durier F., Dalla Vecchia C., 2012, MNRAS, 419, 465

Dutton A. A., Macciò A. V., 2014, MNRAS, 441, 3359

Feldmann R., Gnedin N. Y., Kravtsov A. V., 2011, ApJ, 732, 115

Ferland G. J., Peterson B. M., Horne K., Welsh W. F., Nahar S. N., 1992, ApJ, 387, 95

Ferland G. J., Korista K. T., Verner D. A., Ferguson J. W., Kingdon J. B., Verner E. M., 1998, PASP, 110, 761

Ferland G. J. et al., 2013, Rev. Mex. Astron. Astrofis., 49, 137

Ferrière K. M., 2001, Rev. Mod. Phys., 73, 1031

Fryxell B. et al., 2000, ApJS, 131, 273

Gabor J. M., Capelo P. R., Volonteri M., Bournaud F., Bellovary J., Governato F., Quinn T., 2016, A&A, 592, A62

Genel S., Genzel R., Bouché N., Naab T., Sternberg A., 2009, ApJ, 701, 2002

Glover S. C. O., 2015, MNRAS, 451, 2082

Glover S. C. O., Abel T., 2008, MNRAS, 388, 1627

Glover S. C. O., Jappsen A.-K., 2007, ApJ, 666, 1

Gnat O., Sternberg A., 2007, ApJS, 168, 213

Gnedin N. Y., Kravtsov A. V., 2010, ApJ, 714, 287

Gnedin N. Y., Kravtsov A. V., 2011, ApJ, 728, 88

Gnedin N. Y., Tassis K., Kravtsov A. V., 2009, ApJ, 697, 55

Gonzalez-Perez V., Lacey C. G., Baugh C. M., Lagos C. D. P., Helly J., Campbell D. J. R., Mitchell P. D., 2014, MNRAS, 439, 264
Grassi T., Krstic P., Merlin E., Buonomo U., Piovan L., Chiosi C., 2011, A&A, 533, A123

Grassi T., Bovino S., Schleicher D. R. G., Prieto J., Seifried D., Simoncini E., Gianturco F. A., 2014, MNRAS, 439, 2386

Grassi T., Bovino S., Haugbølle T., Schleicher D. R. G., 2017, MNRAS, 466, 1259

Haardt F., Madau P., 2012, ApJ, 746, 125

Heays A. N., Bosman A. D., van, Dishoeck E. F., 2017, A&A, 602, A105

Hernquist L., 1990, ApJ, 356, 359

Hindmarsh A. C., 1983, IMACS Trans. Sci. Comput., 1, 55

Hollenbach D., McKee C. F., 1979, ApJS, 41, 555

Hopkins P. F., 2015, MNRAS, 450, 53

Hopkins P. F., Narayanan D., Murray N., 2013, MNRAS, 432, 2647

Hopkins P. F. et al., 2017, MNRAS, preprint (arXiv:1707.07010)

Hu C.-Y., Naab T., Walch S., Glover S. C. O., Clark P. C., 2016, MNRAS, 458, 3528

Hu C.-Y., Naab T., Glover S. C. O., Walch S., Clark P. C., 2017, MNRAS, 471, 2151

Huebner W. F., Mukherjee J., 2015, Planet. Space Sci., 106, 11

Janev R. K., Langer W. D., Evans K., 1987, Elementary Processes in Hydrogen-Helium plasmas – Cross-sections and Reaction Rate Coefficients. Springer-Verlag, Berlin

Jeans J. H., 1902, Phil. Trans. Roy. Soc., 199, 1

Jenkins E. B., 2009, ApJ, 700, 1299

Jimenez R., Flynn C., MacDonald J., Gibson B. K., 2003, Science, 299, 1552

Jura M., 1975, ApJ, 197, 575

Karpas Z., Anicich V., Huntress W. T., 1979, J. Chem. Phys., 70, 2877

Katz N., 1992, ApJ, 391, 502

Katz H., Sijacki D., Haehnelt M. G., 2015, MNRAS, 451, 2352

Keller B. W., Wadsley J., Benincasa S. M., Couchman H. M. P., 2014, MNRAS, 442, 3013

Kennicutt R. C., Jr, 1989, ApJ, 344, 685

Kennicutt R. C., Jr, 1998, ApJ, 498, 541

Kennicutt R. C., Jr, Tamblyn P., Congdon C. E., 1994, ApJ, 435, 22

Kimm T., Cen R., Devriendt J., Dubois Y., Slyz A., 2015, MNRAS, 451, 2900

Kimura M., Lane N. F., Dalgarno A., Dixson R. G., 1993, ApJ, 405, 801

Kingdon J. B., Ferland G. J., 1996, ApJS, 106, 205

Körtgen B., Bovino S., Schleicher D. R. G., Giannetti A., Banerjee R., 2017, MNRAS, 469, 2602

Kreckel H., Bruhns H., Čížek M., Glover S. C. O., Miller K. A., Urbain X., Savin D. W., 2010, Science, 329, 69

Kroupa P., 2001, MNRAS, 322, 231

Krumholz M. R., Dekel A., McKee C. F., 2012, ApJ, 745, 69

Le Teuff Y. H., Millar T. J., Markwick A. J., 2000, A&AS, 146, 157

Lenzuni P., Chernoff D. F., Salpeter E. E., 1991, ApJS, 76, 759

Lupi A., Bovino S., Capelo P. R., Volonteri M., Silk J., 2018, MNRAS, 474, 2884

Maio U., Dolag K., Giardi B., Tornatore L., 2007, MNRAS, 379, 963

Malhotra S. et al., 2001, ApJ, 561, 766

Mathis J. S., Rumpl W., Nordsieck K. H., 1977, ApJ, 217, 425

Mckee C. F., Krumholz M. R., 2010, ApJ, 709, 308

Mckee C. F., Ostriker J. P., 1977, ApJ, 218, 148

Michałowski M. J. et al., 2016, A&A, 595, A72

Mitchell G. F., Deveau T. J., 1983, ApJ, 266, 646

Mo H. J., Mao S., White S. D. M., 1998, MNRAS, 295, 319

Monaghan J. J., 1992, ARA&A, 30, 543

Nahar S. N., 1995, ApJS, 101, 423

Nahar S. N., 1996, ApJS, 106, 213

Nahar S. N., 1999, ApJS, 120, 131

Nahar S. N., 2000, ApJS, 126, 537

Nahar S. N., Pradhan A. K., 1997, ApJS, 111, 339

Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563

Nickerson S., Teyssier R., Rosdahl J., 2018, preprint (arXiv:1802.00445)

Omukai K., 2000, ApJ, 534, 809

Omukai K., Tsuribe T., Schneider R., Ferrara A., 2005, ApJ, 626, 627

Oppenheimer B. D., Schaye J., 2013a, MNRAS, 434, 1043

Oppenheimer B. D., Schaye J., 2013b, MNRAS, 434, 1063

Pakmor R., Springel V., 2013, MNRAS, 432, 176

Pallottini A., Ferrara A., Gallerani S., Vallini L., Maiolino R., Salvadori S., 2017a, MNRAS, 465, 2540

Pallottini A., Ferrara A., Bovino S., Vallini L., Gallerani S., Maiolino R., Salvadori S., 2017b, MNRAS, 471, 4128

Papadopoulos P. P., Thi W.-F., 2013, in Torres D. F., Reimer O., eds, Astrophysics and Space Science Proceedings Vol. 34, Cosmic Rays in Star-Forming Environments. Springer-Verlag Berlin, p. 41

Pilbratt G. L. et al., 2010, A&A, 518, L1

Pontzen A., Roškar R., Stinson G., Woods R., 2013, Astrophysics Source Code Library, record ascl:1305.002

Poulaert G., Brouillard F., Claeys W., McGowan J. W., Van Wassenhove G., 1978, J. Phys. B At. Mol. Phys., 11, L671

Raiteri C. M., Villata M., Navarro J. F., 1996, A&A, 315, 105

Richings A. J., Schaye J., 2016, MNRAS, 458, 270

Richings A. J., Schaye J., Oppenheimer B. D., 2014a, MNRAS, 440, 3349

Richings A. J., Schaye J., Oppenheimer B. D., 2014b, MNRAS, 442, 2780

Ritchie B. W., Thomas P. A., 2001, MNRAS, 323, 743

Robertson B. E., Kravtsov A. V., 2008, ApJ, 680, 1083

Rodenbeck K., Schleicher D. R. G., 2016, A&A, 593, A89

Roškar R., Fiacconi D., Mayer L., Kazantzidis S., Quinn T. R., Wadsley J., 2015, MNRAS, 449, 494

Safranek-Shrader C., Krumholz M. R., Kim C.-G., Ostriker E. C., Klein R. I., Li S., McKee C. F., Stone J. M., 2017, MNRAS, 465, 885

Sánchez S. F. et al., 2017, MNRAS, 469, 2121

Savin D. W., Krstić P. S., Haiman Z., Stancil P. C., 2004, ApJ, 606, L167

Scannapieco C. et al., 2012, MNRAS, 423, 1726

Schaye J. et al., 2010, MNRAS, 402, 1536

Schaye J. et al., 2015, MNRAS, 446, 521

Schmidt M., 1959, ApJ, 129, 243

Schmidt M., 1963, ApJ, 137, 758

Schneider I. F., Dulieu O., Giusti-Suzor A., Roueff E., 1994, ApJ, 424, 983

Schuessler I., Schmitt D., 1981, A&A, 97, 373

Segers M. C., Oppenheimer B. D., Schaye J., Richings A. J., 2017, MNRAS, 471, 1026

Seifried D., Walch S., 2016, MNRAS, 459, L11

Shapiro P. R., Kang H., 1987, ApJ, 318, 32

Shen S., Wadsley J., Stinson G., 2010, MNRAS, 407, 1581

Shen S., Madau P., Guedes J., Mayer L., Prochaska J. X., Wadsley J., 2013, ApJ, 765, 89

Spitzer L., Jr, 1942, ApJ, 95, 329

Springel V., White S. D. M., 1999, MNRAS, 307, 162

Springel V., Di Matteo T., Hernquist L., 2005a, MNRAS, 361, 776

Springel V. et al., 2005b, Nature, 435, 629

Stacey G. J., Geis N., Genzel R., Lugten J. B., Poglitsch A., Sternberg A., Townes C. H., 1991, ApJ, 373, 423

Stadel J. G., 2001, PhD thesis, Univ. Washington

Stancil P. C. et al., 1998, ApJ, 502, 1006

Stancil P. C., Schultz D. R., Kimura M., Gu J.-P., Hirsch G., Buenker R. J., 1999, A&AS, 140, 225

Stenrup M., Larson Å., Elander N., 2009, Phys. Rev. A, 79, 012713

Sternberg A., Le Petit F., Roueff E., Le Bourlot J., 2014, ApJ, 790, 10

Stinson G., Seth A., Katz N., Wadsley J., Governato F., Quinn T., 2006, MNRAS, 373, 1074

Sutherland R. S., Dopita M. A., 1993, ApJS, 88, 253

Tacconi L. J. et al., 2010, Nature, 463, 781

Teyssier R., 2002, A&A, 385, 337

Thielemann F.-K., Nomoto K., Yokoi K., 1986, A&A, 158, 17

Tielens A. G. G. M., Hollenbach D., 1985, ApJ, 291, 747

Tomassetti M., Porciani C., Romano-Díaz E., Ludlow A. D., 2015, MNRAS, 446, 3330

Toomre A., 1964, ApJ, 139, 1217

Truelove J. K., Klein R. I., McKee C. F., Holliman II J. H., Howell L. H., Greenough J. A., 1997, ApJ, 489, L179

Verner D. A., Ferland G. J., 1996, ApJS, 103, 467

Verner D. A., Ferland G. J., Korista K. T., Yakovlev D. G., 1996, ApJ, 465, 487
Vitvitska M., Klypin A. A., Kravtsov A. V., Wechsler R. H., Primack J. R., Bullock J. S., 2002, ApJ, 581, 799

Voronov G. S., 1997, At. Data Nucl. Data Tables, 65, 1

Wadsley J. W., Stadel J., Quinn T., 2004, New A, 9, 137

Wadsley J. W., Veeravalli G., Couchman H. M. P., 2008, MNRAS, 387, 427

Wadsley J. W., Keller B. W., Quinn T. R., 2017, MNRAS, 471, 2357

Weidemann V., 1987, A&A, 188, 74

Weingartner J. C., Draine B. T., 2001, ApJ, 563, 842

Wendland H., 1995, Adv. Comput. Math., 4, 389

Wolcott-Green J., Haiman Z., Bryan G. L., 2011, MNRAS, 418, 838

Wolfire M. G., Hollenbach D., McKee C. F., Tielens A. G. G. M., Bakes E. L. O., 1995, ApJ, 443, 152

Wolfire M. G., McKee C. F., Hollenbach D., Tielens A. G. G. M., 2003, ApJ, 587, 278

Woosley S. E., Weaver T. A., 1995, ApJS, 101, 181

Yamasawa D., Habe A., Kozasa T., Nozawa T., Hirashita H., Umeda H., Nomoto K., 2011, ApJ, 735, 44

Yoshida N., Omukai K., Hernquist L., Abel T., 2006, ApJ, 652, 6

Zhao L. B. et al., 2004, ApJ, 615, 1063

Zhukovska S., Gail H.-P., Trieloff M., 2008, A&A, 479, 453

## APPENDIX A: THE IMPORTANCE OF SHIELDING

In this section, we show the effect of not including any shielding from the photoionizing UV background in Runs 01 and 01m. In Fig. A1, we show the gas $\rho-T$ diagrams for the eq-metals Run 01 (upper panel) and the non-eq-metals Run 01m (lower panel) with no shielding. The two diagrams are very similar to each other, except mostly for the density range $10^{-23} \lesssim \rho_{\text{gas}} \lesssim 10^{-22} \text{ g cm}^{-3}$, in which the non-eq-metals is slightly warmer than the eq-metals gas. This implies that the major differences observed between the two panels of Fig. 1 are likely due to the implementation of shielding, with the minor remaining differences due to a combination of different metal networks and non-equilibrium effects (see the discussion in Section 3.1).

## APPENDIX B: REACTION RATES EMPLOYED IN THE PRESENT WORK

The chemical network used in this work is very similar to that described in Bovino et al. (2016; a slightly updated version of that network, named *react_galaxy_ism*, is publicly available with the KROME package under the folder KROME/NETWORKS). However, since the publication of Bovino et al. (2016), a few typos were found and there have been a number of changes, mostly involving the temperature limits of validity of the reactions (we note that, in our simulations, the gas mass fraction of particles with $T_{\text{gas}} > 10^7$ K is $\lesssim 0.01$ per cent and there are zero particles with $T_{\text{gas}} > 10^8$ K). Therefore, for definiteness, we list the reaction rates in Tables B1–B3 and describe the changes from Bovino et al. (2016) in the notes. In this section only, we use the usual chemical notation: $H = H^+$; $H^+ = H_2$, $H_2 = H2^+$; $H_2^+ = H2_+$, $He = He^+$; $He^+ = He_2$, $He^{++} = He_3$, $C = C^+$; $C^+ = C_2$, $O = O^+$; $O^+ = O_2$, $Si = Si^+$; $Si^+ = Si_2$, and $Si^{++} = Si_3$.

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

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

**Figure A1.** Dependence on the modelling of metals, when assuming no shielding. Gas $\rho-T$ diagrams at 0.4 Gyr, for the eq-metals Run 01 (upper panel) and the non-eq-metals Run 01m (lower panel) with no shielding. The lines are the same described in Fig. 1. The differences between the upper and lower diagrams are negligible, when compared to the discrepancies between the panels of Fig. 1.
<table><caption>**Table B1** List of reactions and rates included in our chemical network.</caption><thead><tr><th>Reaction</th><th>Rate coefficient (cm<sup>3</sup> s<sup>-1</sup>)</th><th>Temp. range</th><th>Ref.</th></tr></thead><tbody><tr><td>1. H + e<sup>-</sup> → H<sup>+</sup> + 2e<sup>-</sup></td><td>k<sub>1</sub> = exp[-32.713 967 86 + 13.536 5560 ln T<sub>e</sub><br/>-5.739 328 75 (ln T<sub>e</sub>)<sup>2</sup> + 1.563 154 98 (ln T<sub>e</sub>)<sup>3</sup><br/>-0.287 705 60 (ln T<sub>e</sub>)<sup>4</sup> + 3.482 559 77 × 10<sup>-2</sup>(ln T<sub>e</sub>)<sup>5</sup><br/>-2.631 976 17 × 10<sup>-3</sup>(ln T<sub>e</sub>)<sup>6</sup> + 1.119 543 95 × 10<sup>-4</sup>(ln T<sub>e</sub>)<sup>7</sup><br/>-2.039 149 85 × 10<sup>-6</sup>(ln T<sub>e</sub>)<sup>8</sup>]</td><td></td><td>1</td></tr><tr><td>2. H<sup>+</sup> + e<sup>-</sup> → H + γ</td><td>k<sub>2</sub> = 3.92 × 10<sup>-13</sup> T<sub>e</sub><sup>-0.6353</sup><br/>= exp[-28.613 033 806 892 32<br/>-7.241 125 657 826 851 × 10<sup>-1</sup> ln T<sub>e</sub><br/>-2.026 044 731 984 691 × 10<sup>-2</sup> (ln T<sub>e</sub>)<sup>2</sup><br/>-2.380 861 877 349 834 × 10<sup>-3</sup> (ln T<sub>e</sub>)<sup>3</sup><br/>-3.212 605 213 188 796 × 10<sup>-4</sup> (ln T<sub>e</sub>)<sup>4</sup><br/>-1.421 502 914 054 107 × 10<sup>-5</sup> (ln T<sub>e</sub>)<sup>5</sup><br/>+4.989 108 920 299 513 × 10<sup>-6</sup> (ln T<sub>e</sub>)<sup>6</sup><br/>+5.755 614 137 575 758 × 10<sup>-7</sup> (ln T<sub>e</sub>)<sup>7</sup><br/>-1.856 767 039 775 261 × 10<sup>-8</sup> (ln T<sub>e</sub>)<sup>8</sup><br/>-3.071 135 243 196 595 × 10<sup>-9</sup> (ln T<sub>e</sub>)<sup>9</sup>]</td><td>T ≤ 5500 K<br/>T > 5500 K</td><td>2</td></tr><tr><td>3. He + e<sup>-</sup> → He<sup>+</sup> + 2e<sup>-</sup></td><td>k<sub>3</sub> = exp[-44.098 648 86 + 23.915 965 63 ln T<sub>e</sub><br/>-10.753 2302 (ln T<sub>e</sub>)<sup>2</sup><br/>+3.058 038 75 (ln T<sub>e</sub>)<sup>3</sup><br/>-5.685 1189 × 10<sup>-1</sup> (ln T<sub>e</sub>)<sup>4</sup><br/>+6.795 391 23 × 10<sup>-2</sup> (ln T<sub>e</sub>)<sup>5</sup><br/>-5.009 056 10 × 10<sup>-3</sup> (ln T<sub>e</sub>)<sup>6</sup><br/>+2.067 236 16 × 10<sup>-4</sup> (ln T<sub>e</sub>)<sup>7</sup><br/>-3.649 161 41 × 10<sup>-6</sup> (ln T<sub>e</sub>)<sup>8</sup>]</td><td></td><td>1</td></tr><tr><td>4. He<sup>+</sup> + e<sup>-</sup> → He + γ</td><td>k<sub>4</sub> = 3.92 × 10<sup>-13</sup> T<sub>e</sub><sup>-0.6353</sup><br/>= 3.92 × 10<sup>-13</sup> T<sub>e</sub><sup>-0.6353</sup><br/>+1.54 × 10<sup>-9</sup> T<sub>e</sub><sup>-1.5</sup> [1.0 + 0.3 / exp (8.099 328 789 667/T<sub>e</sub>)]<br/>/[exp (40.496 643 948 336 62/T<sub>e</sub>)]</td><td>T<sub>e</sub> ≤ 0.8 eV<br/>T<sub>e</sub> > 0.8 eV</td><td>3</td></tr><tr><td>5. He<sup>+</sup> + e<sup>-</sup> → He<sup>++</sup> + 2e<sup>-</sup></td><td>k<sub>5</sub> = exp[-68.710 409 902 120 01 + 43.933 476 326 35 ln T<sub>e</sub><br/>-18.480 669 935 68 (ln T<sub>e</sub>)<sup>2</sup><br/>+4.701 626 486 759 002 (ln T<sub>e</sub>)<sup>3</sup><br/>-7.692 466 334 492 × 10<sup>-1</sup> (ln T<sub>e</sub>)<sup>4</sup><br/>+8.113 042 097 303 × 10<sup>-2</sup> (ln T<sub>e</sub>)<sup>5</sup><br/>-5.324 020 628 287 001 × 10<sup>-3</sup> (ln T<sub>e</sub>)<sup>6</sup><br/>+ 1.975 705 312 221 × 10<sup>-4</sup> (ln T<sub>e</sub>)<sup>7</sup><br/>-3.165 581 065 665 × 10<sup>-6</sup> (ln T<sub>e</sub>)<sup>8</sup>]</td><td></td><td>4</td></tr><tr><td>6. He<sup>++</sup> + e<sup>-</sup> → He<sup>+</sup> + γ</td><td>k<sub>6</sub> = 1.891 × 10<sup>-10</sup> [(1.0 + √T/9.37)<sup>0.2476</sup> × (1.0 + √T/(2.774 × 10<sup>6</sup>))<sup>1.7524</sup> (√T/9.37)]<sup>-1</sup></td><td></td><td>5</td></tr><tr><td>7. H + e<sup>-</sup> → H<sup>-</sup> + γ</td><td>k<sub>7</sub> = 1.4 × 10<sup>-18</sup>T<sup>0.928</sup>exp(-T/16 200)</td><td></td><td>6</td></tr><tr><td>8. H<sup>-</sup> + H → H<sub>2</sub> + e<sup>-</sup></td><td>k<sub>8</sub> = a<sub>1</sub>(T<sup>a2</sup> + a<sub>3</sub>T<sup>a4</sup> + a<sub>5</sub>T<sup>a6</sup>)/(1.0 + a<sub>7</sub>T<sup>a8</sup> + a<sub>9</sub>T<sup>a10</sup> + a<sub>11</sub>T<sup>a12</sup>)<br/>a<sub>1</sub> = 1.35 × 10<sup>-9</sup><br/>a<sub>2</sub> = 9.8493 × 10<sup>-2</sup><br/>a<sub>3</sub> = 3.2852 × 10<sup>-1</sup><br/>a<sub>4</sub> = 5.5610 × 10<sup>-1</sup><br/>a<sub>5</sub> = 2.7710 × 10<sup>-7</sup><br/>a<sub>6</sub> = 2.1826<br/>a<sub>7</sub> = 6.1910 × 10<sup>-3</sup><br/>a<sub>8</sub> = 1.0461<br/>a<sub>9</sub> = 8.9712 × 10<sup>-11</sup><br/>a<sub>10</sub> = 3.0424<br/>a<sub>11</sub> = 3.2576 × 10<sup>-14</sup><br/>a<sub>12</sub> = 3.7741</td><td></td><td>7</td></tr><tr><td>9. H + H<sup>+</sup> → H<sub>2</sub><sup>+</sup> + γ</td><td>k<sub>9</sub> = 2.10 × 10<sup>-20</sup>(T/30)<sup>-0.15</sup><br/>= dex[-18.20 - 3.194 log<sub>10</sub>T + 1.786(log<sub>10</sub>T)<sup>2</sup> - 0.2072(log<sub>10</sub>T)<sup>3</sup>]<br/>k<sub>10</sub> = 6.0 × 10<sup>-10</sup></td><td>T < 30 K<br/>T ≥ 30 K</td><td>8</td></tr><tr><td>10. H<sub>2</sub><sup>+</sup> + H → H<sub>2</sub> + H<sup>+</sup></td><td>k<sub>11</sub> = dex[∑<sup>7</sup><sub>i=0</sub> a<sub>i</sub> [log<sub>10</sub>(T)]<sup>i</sup> (see Table B2)]</td><td></td><td>9</td></tr><tr><td>11. H<sub>2</sub> + H<sup>+</sup> → H<sub>2</sub><sup>+</sup> + H</td><td>k<sub>12</sub> = exp[-18.018 493 342 73<br/>+ 2.360 852 208 681 ln T<sub>e</sub><br/>-2.827 443 061 704 × 10<sup>-1</sup> (ln T<sub>e</sub>)<sup>2</sup><br/>+1.623 316 639 567 × 10<sup>-2</sup> (ln T<sub>e</sub>)<sup>3</sup><br/>-3.365 012 031 362 999 × 10<sup>-2</sup> (ln T<sub>e</sub>)<sup>4</sup><br/>+1.178 329 782 711 × 10<sup>-2</sup> (ln T<sub>e</sub>)<sup>5</sup><br/>-1.656 194 699 504 × 10<sup>-3</sup> (ln T<sub>e</sub>)<sup>6</sup><br/>+1.068 275 202 678 × 10<sup>-4</sup> (ln T<sub>e</sub>)<sup>7</sup><br/>-2.631 285 809 207 × 10<sup>-6</sup> (ln T<sub>e</sub>)<sup>8</sup>]</td><td>T ≤ 10<sup>8</sup> K</td><td>10</td></tr></tbody></table>
<table><caption>**Table B1 – continued**</caption><thead><tr><th>Reaction</th><th>Rate coefficient (cm<sup>3</sup> s<sup>-1</sup>)</th><th>Temp. range</th><th>Ref.</th></tr></thead><tbody><tr><td>13. H<sup>-</sup> + H → 2H + e<sup>-</sup></td><td>k<sub>13</sub> = 2.56 × 10<sup>-9</sup> T<sub>e</sub><sup>1.78186</sup><br>= exp[-20.372 608 965 333 24<br>+1.139 449 335 841 631 ln T<sub>e</sub>)<sup>2</sup><br>+8.464 455 386 63 × 10<sup>-3</sup> (ln T<sub>e</sub>)<sup>3</sup><br>-1.432 764 121 2992 × 10<sup>-3</sup> (ln T<sub>e</sub>)<sup>4</sup><br>+2.012 250 284 791 × 10<sup>-4</sup> (ln T<sub>e</sub>)<sup>5</sup><br>+8.663 963 243 09 × 10<sup>-5</sup> (ln T<sub>e</sub>)<sup>6</sup><br>-2.585 009 680 264 × 10<sup>-5</sup> (ln T<sub>e</sub>)<sup>7</sup><br>+2.455 501 197 0392 × 10<sup>-6</sup> (ln T<sub>e</sub>)<sup>8</sup><br>-8.068 382 461 18 × 10<sup>-8</sup> (ln T<sub>e</sub>)<sup>9</sup></td><td>T<sub>e</sub> ≤ 0.1 eV<br>T<sub>e</sub> > 0.1 eV</td><td>11</td></tr><tr><td>14. H<sup>-</sup> + H<sup>+</sup> → 2H + γ</td><td>k<sub>14</sub> = 2.96 × 10<sup>-6</sup>/√T - 1.73 × 10<sup>-9</sup> + 2.50 × 10<sup>-10</sup>√T<br>-7.77 × 10<sup>-13</sup></td><td>T > 10 K</td><td>12</td></tr><tr><td>15. H<sup>-</sup> + H<sup>+</sup> → H<sub>2</sub><sup>+</sup> + e<sup>-</sup></td><td>k<sub>15</sub> = 10<sup>-8</sup> × T<sup>-0.4</sup></td><td>T ≤ 617 K</td><td>13</td></tr><tr><td>16. H<sub>2</sub><sup>+</sup> + e<sup>-</sup> → 2H + γ</td><td>k<sub>16</sub> = 1.0 × 10<sup>-8</sup><br>= 1.32 × 10<sup>-6</sup>T<sup>-0.76</sup></td><td>T > 617 K</td><td>14</td></tr><tr><td>17. H<sub>2</sub><sup>+</sup> + H<sup>-</sup> → H + H<sub>2</sub></td><td>k<sub>17</sub> = 5.0 × 10<sup>-7</sup> (10<sup>2</sup>/T)<sup>0.5</sup></td><td></td><td>15</td></tr><tr><td>18. H<sub>2</sub> + H → 3H</td><td>k<sub>18</sub> = 6.67 × 10<sup>-12</sup>√T exp[-(1 + 63 593/T)]</td><td></td><td>16</td></tr><tr><td>19. H<sub>2</sub> + H<sub>2</sub> → H<sub>2</sub> + H + H</td><td>k<sub>19</sub> = 5.996 × 10<sup>-30</sup>T<sup>4.1881</sup>(1 + 6.761 × 10<sup>-6</sup>T)<sup>-5.6881</sup><br>× exp(-54 657.4/T)</td><td></td><td>16</td></tr><tr><td>20. He<sup>+</sup> + H → He + H<sup>+</sup></td><td>k<sub>20</sub> = 1.20 × 10<sup>-15</sup>T(300)<sup>0.25</sup></td><td></td><td>17</td></tr><tr><td>21. He + H<sup>+</sup> → He<sup>+</sup> + H</td><td>k<sub>21</sub> = 1.26 × 10<sup>-9</sup>T<sup>-0.75</sup>exp(-1.275 × 10<sup>5</sup>/T)<br>= 4 × 10<sup>-37</sup>T<sup>4.74</sup></td><td>T ≤ 10<sup>4</sup> K<br>T > 10<sup>4</sup> K</td><td>18</td></tr><tr><td>22. H<sub>2</sub> + e<sup>-</sup> → H + H + e<sup>-</sup></td><td>k<sub>22</sub> = 4.38 × 10<sup>-10</sup>T<sup>0.35</sup>exp(-102 000/T)</td><td></td><td>19</td></tr><tr><td>23. H<sub>2</sub> + e<sup>-</sup> → H + H<sup>-</sup></td><td>k<sub>23</sub> = 35.5 × T<sup>-2.28</sup>exp(-46707/T)</td><td></td><td>20</td></tr><tr><td>24. H<sub>2</sub> + He → H + H + He</td><td>k<sub>24</sub> = dex[-27.029 + 3.801 log<sub>10</sub>(T) - 29487/T]</td><td></td><td>21</td></tr><tr><td>25. H<sub>2</sub> + He<sup>+</sup> → He + H + H<sup>+</sup></td><td>k<sub>25</sub> = 3.7 × 10<sup>-14</sup> exp(-35.0/T)</td><td></td><td>22</td></tr><tr><td>26. H<sub>2</sub> + He<sup>+</sup> → H<sub>2</sub><sup>+</sup> + He</td><td>k<sub>26</sub> = 7.2 × 10<sup>-15</sup></td><td></td><td>22</td></tr><tr><td>27. H + H → H + H<sup>+</sup> + e<sup>-</sup></td><td>k<sub>27</sub> = 1.2 × 10<sup>-17</sup>T<sup>1.2</sup>exp(-157800/T)</td><td></td><td>24</td></tr><tr><td>28. H + He → He + H<sup>+</sup> + e<sup>-</sup></td><td>k<sub>28</sub> = 1.75 × 10<sup>-17</sup>T<sup>1.3</sup>exp(-157800/T)</td><td></td><td>24</td></tr><tr><td>29. H<sup>+</sup> + e<sup>-</sup> → H<sup>dust</sup></td><td>k<sub>29</sub> = 1.225 × 10<sup>-13</sup>(Z/Z<sub>☉</sub>)[1.0 + 8.074 × 10<sup>-6</sup>ψ<sup>1.378</sup><br>× (1 + 5.087 × 10<sup>2</sup>T<sup>0.01586</sup>ψ<sup>-0.4723-1.102×10<sup>-5</sup>ln(T)</sup>]<sup>-1</sup></td><td></td><td>25</td></tr><tr><td>30. He<sup>+</sup> + e<sup>-</sup> → He<sup>dust</sup></td><td>k<sub>30</sub> = 5.572 × 10<sup>-14</sup>(Z/Z<sub>☉</sub>)[1.0 + 6.089 × 10<sup>-3</sup>ψ<sup>1.1728</sup><br>× (1 + 4.331 × 10<sup>2</sup>T<sup>0.04845</sup>ψ<sup>-0.8120-1.333×10<sup>-4</sup>ln(T)</sup>]<sup>-1</sup></td><td></td><td>25</td></tr><tr><td>31. H + H → H<sup>dust</sup></td><td>see text</td><td></td><td>26</td></tr><tr><td>32. C<sup>+</sup> + e<sup>-</sup> → C + γ</td><td>k<sub>32</sub> = 4.67 × 10<sup>-12</sup>(T/300)<sup>-0.6</sup><br>= 1.23 × 10<sup>-17</sup>(T/300)<sup>2.49</sup>exp((21 845.6)/T)<br>= 9.62 × 10<sup>-8</sup>(T/300)<sup>-1.37</sup>exp(-115 786.2/T)</td><td>T ≤ 7950 K<br>7950 &lt; T ≤ 21 140 K<br>T &gt; 21 140 K</td><td>27</td></tr><tr><td>33. Si<sup>+</sup> + e<sup>-</sup> → Si + γ</td><td>k<sub>33</sub> = 7.5 × 10<sup>-12</sup>(T/300)<sup>-0.55</sup><br>= 4.86 × 10<sup>-12</sup>(T/300)<sup>-0.32</sup><br>= 9.08 × 10<sup>-14</sup>(T/300)<sup>0.818</sup></td><td>T ≤ 2000 K<br>2000 &lt; T ≤ 10<sup>4</sup> K<br>T &gt; 10<sup>4</sup> K</td><td>28</td></tr><tr><td>34. O<sup>+</sup> + e<sup>-</sup> → O + γ</td><td>k<sub>34</sub> = 1.30 × 10<sup>-10</sup>T<sup>-0.64</sup><br>= 1.41 × 10<sup>-10</sup>T<sup>-0.66</sup> + 7.4 × 10<sup>-4</sup>T<sup>-1.5</sup><br>× exp(-175 000/T)[1.0 + 0.062 × exp(-145 000/T)]</td><td>T ≤ 400 K<br>T &gt; 400 K</td><td>29</td></tr><tr><td>35. C + e<sup>-</sup> → C<sup>+</sup> + 2e<sup>-</sup></td><td>k<sub>35</sub> = 6.85 × 10<sup>-8</sup>(0.193 + u)<sup>-1</sup>u<sup>0.25</sup>e<sup>-u</sup></td><td>u = 11.26/T<sub>e</sub></td><td>30</td></tr><tr><td>36. Si + e<sup>-</sup> → Si<sup>+</sup> + 2e<sup>-</sup></td><td>k<sub>36</sub> = 1.88 × 10<sup>-7</sup>(1.0 + u<sup>0.5</sup>)(0.376 + u)<sup>-1</sup>u<sup>0.25</sup>e<sup>-u</sup></td><td>u = 8.2/T<sub>e</sub></td><td>30</td></tr><tr><td>37. O + e<sup>-</sup> → O<sup>+</sup> + 2e<sup>-</sup></td><td>k<sub>37</sub> = 3.59 × 10<sup>-8</sup>(0.073 + u)<sup>-1</sup>u<sup>0.34</sup>e<sup>-u</sup></td><td>u = 13.6/T<sub>e</sub></td><td>30</td></tr><tr><td>38. O<sup>+</sup> + H → O + H<sup>+</sup></td><td>k<sub>38</sub> = 4.99 × 10<sup>-11</sup>T<sup>0.405</sup> + 7.54 × 10<sup>-10</sup>T<sup>-0.458</sup></td><td></td><td>31</td></tr><tr><td>39. O + H<sup>+</sup> → O<sup>+</sup> + H</td><td>k<sub>39</sub> = [1.08 × 10<sup>-11</sup>T<sup>0.517</sup> + 4.00 × 10<sup>-10</sup>T<sup>0.00669</sup>]<exp(-227/T)</td><td></td><td>32</td></tr><tr><td>40. O + He<sup>+</sup> → O<sup>+</sup> + He</td><td>k<sub>40</sub> = 4.991 × 10<sup>-15</sup>(T/1000)<sup>0.3794</sup>exp(-1121 000/T)<br>+ 2.780 × 10<sup>-15</sup>(T/1000)<sup>-0.2163</sup>exp((815 800)/T)</td><td>T &gt; 10 K</td><td>33</td></tr><tr><td>41. C + H<sup>+</sup> → C<sup>+</sup> + H</td><td>k<sub>41</sub> = 3.9 × 10<sup>-16</sup>T<sup>0.213</sup></td><td></td><td>32</td></tr><tr><td>42. C<sup>+</sup> + H → C + H<sup>+</sup></td><td>k<sub>42</sub> = 6.08 × 10<sup>-14</sup>(T/1000)<sup>1.96</sup>exp(-170 000/T)</td><td></td><td>32</td></tr><tr><td>43. C + He<sup>+</sup> → C<sup>+</sup> + He</td><td>k<sub>43</sub> = 8.58 × 10<sup>-17</sup>T<sup>0.757</sup><br>= 3.25 × 10<sup>-17</sup>T<sup>0.968</sup><br>= 2.77 × 10<sup>-19</sup>T<sup>1.597</sup></td><td>T ≤ 200 K<br>200 &lt; T ≤ 2000 K<br>T &gt; 2000 K</td><td>34</td></tr><tr><td>44. Si + H<sup>+</sup> → Si<sup>+</sup> + H</td><td>k<sub>44</sub> = 5.88 × 10<sup>-13</sup>T<sup>0.848</sup><br>= 1.45 × 10<sup>-13</sup>T</td><td>T ≤ 10<sup>4</sup> K<br>T &gt; 10<sup>4</sup> K</td><td>35</td></tr><tr><td>45. Si + He<sup>+</sup> → Si<sup>+</sup> + He</td><td>k<sub>45</sub> = 3.3 × 10<sup>-9</sup></td><td></td><td>36</td></tr><tr><td>46. C<sup>+</sup> + Si → C + Si<sup>+</sup></td><td>k<sub>46</sub> = 2.1 × 10<sup>-9</sup></td><td></td><td>36</td></tr></tbody></table>
<table><caption>**Table B1 – continued**</caption><thead><tr><th>Reaction</th><th>Rate coefficient (cm<sup>3</sup> s<sup>-1</sup>)</th><th>Temp. range</th><th>Ref.</th></tr></thead><tbody><tr><td>47. Si<sup>+</sup> + H<sup>+</sup> → Si<sup>++</sup> + H</td><td>k<sub>47</sub> = 4.10 × 10<sup>-10</sup> (<sup>T</sup>/<sub>10000</sub>)<sup>0.24</sup><br/>× [1.0 + 3.17 exp(<sup>-T</sup>/<sub>2.39×10<sup>6</sup></sub>)] exp(<sup>-3.178</sup>/<sub>T<sub>e</sub></sub>)</td><td></td><td>35</td></tr><tr><td>48. Si<sup>++</sup> + H → Si<sup>+</sup> + H<sup>+</sup></td><td>k<sub>48</sub> = 1.23 × 10<sup>-9</sup> (<sup>T</sup>/<sub>10000</sub>)<sup>0.24</sup><br/>× [1.0 + 3.17 exp(<sup>-T</sup>/<sub>2.39×10<sup>6</sup></sub>)]</td><td></td><td>35</td></tr><tr><td>49. Si<sup>++</sup> + e<sup>-</sup> → Si<sup>+</sup> + γ</td><td>k<sub>49</sub> = 1.75 × 10<sup>-12</sup> (<sup>T</sup>/<sub>10000</sub>)<sup>-0.6346</sup></td><td></td><td>37</td></tr><tr><td>50. C<sup>+</sup> + e<sup>-dust→C</sup></td><td>k<sub>50</sub> = 4.558 × 10<sup>-13</sup>(Z/Z<sub>☉</sub>)[1.0 + 6.089 × 10<sup>-3</sup>ψ<sup>1.128</sup><br/>× (1.0 + 4.331 × 10<sup>2</sup>T<sup>0.04845</sup>ψ<sup>-0.8120-1.333×10<sup>-4</sup> ln T</sup>]<sup>-1</sup></td><td></td><td>25</td></tr><tr><td>51. O<sup>+</sup> + e<sup>-dust→O</sup></td><td>k<sub>51</sub> = <sup>1</sup>/<sub>4</sub>k<sub>29</sub></td><td></td><td>25,38</td></tr><tr><td>52. Si<sup>+</sup> + e<sup>-dust→Si</sup></td><td>k<sub>52</sub> = 2.166 × 10<sup>-14</sup>(Z/Z<sub>☉</sub>)[1.0 + 5.678 × 10<sup>-8</sup>ψ<sup>1.874</sup><br/>× (1.0 + 4.375 × 10<sup>4</sup>T<sup>1.635×10<sup>-6</sup></sup>ψ<sup>-0.8964-7.538×10<sup>-5</sup> ln T</sup>]<sup>-1</sup></td><td></td><td>25</td></tr></tbody></table>

References – 1: Janev, Langer & Evans (1987); 2: Abel et al. (1997) fit by data from Ferland et al. (1992); 3: Cen (1992); Aldrovandi & Pequignot (1973); 4: Aladdin data base https://www-amdis.iaea.org/ALADDIN/, see Abel et al. (1997); 5: Verner & Ferland (1996); 6: de Jong (1972); 7: Kreckel et al. (2010); 8: Coppola et al. (2011); 9: Karpas, Anicich & Huntress (1979); 10: Savin et al. (2004) and Grassi et al. (2011), see table B2 for the coefficients; 11: Abel et al. (1997), based on Janev et al. (1987); 12: Stenrup, Larson & Elander (2009); 13: Poulaert et al. (1978); 14: Abel et al. (1997), fit from data by Schneider et al. (1994); 15: Dalgarno & Lepp (1987); 16: Glover & Abel (2008); 17: Yoshida et al. (2006); 18: Kimura et al. (1993); 19: Mitchell & Deveau (1983), fit of data by Corrigan (1965); 20: Capitelli et al. (2007); 21: Dove et al. (1987); 22: Barlow (1984); 23: Lenzuni, Chernoff & Salpeter (1991); 24: Weingartner & Draine (2001); 25: Cazaux & Spaans (2009), see also Grassi et al. (2014); 26: Nahar & Pradhan (1997); 27: Nahar (2000); 28: Nahar (1999); 29: Voronov (1997); 30: Stancil et al. (1999); 31: Stancil et al. (1998); 32: Zhao et al. (2004); 33: Kimura et al. (1993); 34: Kingdon & Ferland (1996); 35: Le Teuff, Millar & Markwick (2000); 36: Nahar (1995, 1996); 37: Glover & Jappsen (2007).

Notes – In this and the next tables, $T = T_{gas}$, $T_e = k_B T$, where $k_B$ is the Boltzmann constant, is the gas temperature in eV. Reaction 6 and one coefficient in Reaction 8 were correct in the code used by Bovino et al. (2016) but incorrect in the published paper. The parameter $\psi = G_0 T^{1/2}/n_e$, where $G_0$ is the Habing flux and $n_e$ is the electron number density, is the charging parameter. Reaction 17 and the temperature ranges of Reaction 13 were correct in the code used by Bovino et al. (2016) but incorrect in the published paper. Reactions 14 and 25 and the temperature ranges of Reactions 11 (see also Table B2) and 14 were changed from the code used by Bovino et al. (2016). Reactions 47 and 48 and the temperature ranges of Reactions 40, 47, and 48 were changed from the code used by Bovino et al. (2016). The functions and temperature ranges given for Reactions 40, 47, and 48 are slightly different from what given in the references.

<table><caption>**Table B2 Fitting coefficients for Reaction 11 (see Table B1). Coefficients are in the form $a(b) = a \times 10^b$**</caption><thead><tr><th>a<sub>i</sub></th><th>T &lt; 10<sup>5</sup> K</th><th>10<sup>5</sup> &le; T &le; 10<sup>8</sup> K</th></tr></thead><tbody><tr><td>a<sub>0</sub></td><td>-1.915 3214(+2)</td><td>-8.875 5774(+3)</td></tr><tr><td>a<sub>1</sub></td><td>4.012 9114(+2)</td><td>1.008 1246(+4)</td></tr><tr><td>a<sub>2</sub></td><td>-3.744 6991(+2)</td><td>-4.860 6622(+3)</td></tr><tr><td>a<sub>3</sub></td><td>1.907 8410(+2)</td><td>1.288 9659(+3)</td></tr><tr><td>a<sub>4</sub></td><td>-5.726 3467(+1)</td><td>-2.031 9575(+2)</td></tr><tr><td>a<sub>5</sub></td><td>1.013 3210(+1)</td><td>1.905 7493(+1)</td></tr><tr><td>a<sub>6</sub></td><td>-9.801 2853(-1)</td><td>-9.853 0668(-1)</td></tr><tr><td>a<sub>7</sub></td><td>4.002 3414(-2)</td><td>2.167 5387(-2)</td></tr></tbody></table>
<table><caption>**Table B3.** List of photoprocesses included in our chemical network.</caption><thead><tr><th>Reaction</th><th>Cross-section (cm<sup>2</sup>)</th><th>Energy range (eV)</th><th>Ref.</th></tr></thead><tbody><tr><td>P1. H + γ → H<sup>+</sup> + e<sup>-</sup></td><td>KROME data base</td><td></td><td>1</td></tr><tr><td>P2. He + γ → He<sup>+</sup> + e<sup>-</sup></td><td>KROME data base</td><td></td><td>1</td></tr><tr><td>P3. He<sup>+</sup> + γ → He<sup>++</sup> + e<sup>-</sup></td><td>KROME data base</td><td></td><td>1</td></tr><tr><td>P4. H<sup>-</sup> + γ → H + e<sup>-</sup></td><td>SWRI data base</td><td>E > 0.755</td><td>2</td></tr><tr><td>P5. H<sub>2</sub> + γ → H<sub>2</sub><sup>+</sup> + e<sup>-</sup></td><td>SWRI data base</td><td>E > 15.4</td><td>2</td></tr><tr><td>P6. H<sub>2</sub><sup>+</sup> + γ → H<sup>+</sup> + H</td><td>LEIDEN data base</td><td>E > 2.65</td><td>3</td></tr><tr><td>P7. H<sub>2</sub><sup>+</sup> + γ → H<sup>+</sup> + H<sup>+</sup> + e<sup>-</sup></td><td>dex[ - 16.926 - 4.528 × 10<sup>-2</sup>E + 2.238 × 10<sup>-4</sup>E<sup>2</sup> + 4.245 × 10<sup>-7</sup>E<sup>3</sup>]</td><td>30 < E < 90</td><td>4</td></tr><tr><td>P8. H<sub>2</sub> + γ → H + H</td><td>Direct</td><td>14.159 < E < 17.7</td><td>5</td></tr><tr><td>P9. H<sub>2</sub> + γ → H<sub>2</sub><sup>+</sup> → H + H</td><td>Solomon</td><td></td><td>6</td></tr><tr><td>P10. C + γ → C<sup>+</sup> + e<sup>-</sup></td><td>KROME data base</td><td></td><td>1</td></tr><tr><td>P11. O + γ → O<sup>+</sup> + e<sup>-</sup></td><td>KROME data base</td><td></td><td>1</td></tr><tr><td>P12. Si + γ → Si<sup>+</sup> + e<sup>-</sup></td><td>KROME data base</td><td></td><td>1</td></tr><tr><td>P13. Si<sup>+</sup> + γ → Si<sup>++</sup> + e<sup>-</sup></td><td>KROME data base</td><td></td><td>1</td></tr></tbody></table>

References. 1: Verner & Ferland (1996); Verner et al. (1996); 2: http://phidrates.space.swri.edu; 3: http://home.strw.leidenuniv.nl/~ewine/photo/; 4: Shapiro & Kang (1987); 5: Abel et al. (1997); 6: this work, based on the formula reported by Glover & Jappsen (2007); 7: Huebner & Mukherjee (2015); 8: Heays, Bosman & van Dishoeck (2017).

This paper has been typeset from a TeX/LATEX file prepared by the author.