© ESO 2018

1. Introduction

The composition of gas in the interstellar medium (ISM) can be subject to changes via chemical processes of great diversity. Knowledge of the underlying chemistry is potentially important because species abundances regulate the thermal properties of the gas. Thermal processes like radiation cooling can have, in turn, a strong impact on the gas dynamics. When the dynamics and thermal effects of gas get strongly coupled with its chemical changes, non-equilibrium behavior can occur. The assumption of instantaneous chemical equilibrium then becomes disputable and a time-dependent treatment of chemical processes seems imperative. For example, ionization equilibrium is quite assumed. This assumption breaks down, however, when gas cooling occurs more rapidly than ion recombination. In this case, ionization lags result, which lead to an enhancement or reduction of radiation cooling depending on the temperature regime (see, e.g., Gnat & Sternberg 2007; Oppenheimer & Schaye 2013a). In general, cooling rates can become rather inaccurate if chemical equilibrium is predisposed. The assumption of instantaneous chemical equilibrium can also break down in a flow that is sufficiently turbulent. If the eddy turnover time is short enough, physical conditions experienced by a gas parcel may change before chemical equilibrium can be established (Gray et al. 2015). Furthermore, in the presence of an ionizing radiation field that varies on a timescale shorter than the recombination time, the ionization balance is expected to be out of equilibrium and non-equilibrium effects come into play (Oppenheimer & Schaye 2013b).

Atomic and molecular line cooling is central to many ISM problems. For a reasonable modeling of non-equilibrium cooling efficiencies, the chemical state of the gas must be known. This requires us to solve a chemical network in accordance with the gas-dynamical equations. The most comprehensive chemical networks for the ISM are found in astrochemical databases like the UMIST database (Le Teuff et al. 2000) or the Kinetic Database for Astrochemistry (KIDA, Wakelam et al. 2012), which contain thousands of reactions between hundreds of species. Combining such sophisticated networks with large-scale dynamical simulations is far out of reach, however, due to prohibitive computational costs. One must therefore resort to reduced networks that are optimally constructed in such a way that relevant features of larger networks are subsumed. Of particular interest in this respect is to identify the most important coolants and to capture those reactions that have a significant influence on the abundance of those coolants. The selection of a proper subset of reactions from a larger network is also guided by the temperature and density range of the problem or diagnostic purposes. The size of a chemical network that can be tackled within gas-dynamical simulations is severely limited in any case by the available computing power to solve its rateequations.

An example of a reduced network is that of Nelson & Langer (1999) focusing on the formation of carbon monoxide, CO. In this simplified network, the effects of intermediate hydro-carbon species like CH or CH₂ and (hydro-)oxygen species like OH and H₂O are consolidated in the form of (artificial) composite species CH_x and OH_x. This greatly reduces the network complexity at the expense of not tracking distributions of such intermediate species. The benefit of the Nelson & Langer (1999) approach has been proved by Glover & Clark (2012) who compared different networks for modeling CO chemistry in molecular clouds. Another example of a useful simplification is the chemical network of Glover & Jappsen (2007) designed to study low-density, low-temperature gas with low metal content. This network consists of 74 reactions that evolve 18 species. It includes the molecular coolants H₂ and HD but ignores other molecular coolants like CO and H₂O (relevant at higher densities), and it neglects higher ionization states of metals like C, O, and Si (relevant at higher temperatures). This network has been later extended by adding carbon and oxygen molecular chemistry in the work of Glover et al. (2010; GFM10). Recently, a chemistry and cooling model for the diffuse ISM has been constructed by Richings et al. (2014; RSO14) taking the GFM10 chemical network as a basis to represent the low-temperature regime (T ≲ 2 × 10⁴ K), and combining it with the ionization and cooling model of Oppenheimer & Schaye (2013a) appropriate in the high-temperature regime (T ≳ 2 × 10⁴ K).

In an earlier paper (Ziegler 2016), a chemical reaction network solver for the magneto-hydrodynamics code NIRVANA (Ziegler 2008, 2011, 2012) was developed. For the solution of the ordinary differential equation (ODE) system resulting from the rate equations of chemical kinetics and thermal processes, a fourth-order stiff integrator of Rosenbrock-Wanner type was implemented. The thermo-chemistry solver was coupled to the magneto-hydrodynamics scheme in a Strang-type operator-splitting fashion. In conjunction with this development, a very rudimentary chemical model for the ISM was introduced describing H/He ionization and gas-phase H₂ formation through the well-studied H⁻/${\text{H}}_2^{+}$ channel.

In this paper a substantially extended NIRVANA chemistry and cooling module (NC^​2M) is presented inspired by the work of RSO14. To begin with, the NC^​2M comes with two simplifications compared to RSO14, which significantly reduce the network size. First, the less abundant elements S and Ca are neglected because they are minor coolants in interstellar gas. Second, any photoionizing background is ignored at the current stage of development. This allows us to cut down photochemical reactions and to dispense with the process of Auger ionization by energetic photons. The second assumption is a critical one and, to be valid in reality, requires the gas to be somehow shielded from the background ultraviolet radiation field of the ISM. This is the case, for instance, in the interior regions of dense clouds where chemical processes are controlled by gas-phase reactions, cosmic ray ionization, and by reactions on grain surfaces, and where photoionization and photodissociation effects play only a minor role. In such shielded regimes, which generally tend to be cold, the high-temperature chemistry included in the NC^​2M becomes non-essential and can be switched off as an option. Photochemical processes will be added in a future version of the NC^​2M to account for the effects of an ionizing background.

There are further differences to the work of RSO14. The RSO14 network does not consider deuterium chemistry. Following Omukai et al. (2005) and Glover & Abel (2008), it may become non-negligible in dense, low-metallicity gas. Although the NC^​2M aims at applications for metal-enriched gases in the ISM, a mini-network of deuterium chemistry has been included as an option. Fine-structure metal line cooling is treated in the NC^​2M from first principles by solving the statistical equilibrium equation for the ion level populations, thereby making use of the Chianti atomic database (Del Zanna et al. 2015). In RSO14 ion-by-ion cooling tables are applied that were precomputed with help of the spectral synthesis code Cloudy (Ferland et al. 2013). Additional minor differences affect the used kinetic- and thermal rate coefficients, the number of charge transfer reactions between metal species, details in the dust model, and the more elaborate treatment of cosmic-ray heating in the NC^​2M.

2. The chemical model

Following Ziegler (2016), the chemical evolution extends the system of gas-dynamical equations by adding advection-reaction equations for the number densities n_s of different species X_s, s = 1…N_s, of the form

$$\begin{array}{c}\hfill {\partial }_{\mathrm{t}}{n}_s+\mathrm{\nabla }·(n_s\mathbf{v})=R_s={\displaystyle \sum _{r=1}^{N_r}({\beta }_{s,r}-{\alpha }_{s,r})k_r\prod _{\begin{array}{c}l\\ X_l\in R(r)\end{array}}{n}_l^{{\alpha }_{l,r}},}\end{array}$$(1)

where v is the flow field, k_r are the kinetic rate coefficients for the chemical reactions, r = 1…N_r, and α and β are stoichiometric coefficients. The reaction source terms R_s in Eq. (1) derive from the set of chemical reactions

$$\begin{array}{c}\hfill {\displaystyle \sum _{\begin{array}{c}s\\ X_s\in R(r)\end{array}}{\alpha }_{s,r}{X}_s\stackrel{k_r}{\to }\sum _{\begin{array}{c}s\\ X_s\in P(r)\end{array}}{\beta }_{s,r}{X}_s,r=1\dots N_r}\end{array}$$

between reactants X_s ∈ R(r) leading to products X_s ∈ P(r) with R(r) and P(r) being subsets of ${\left\{X_s\right\}}_{s=1}^{N_s}$.

The NC^​2M evolves 121 species. Referring to standard solar abundances as cataloged in Grevesse & Sauval (1998), the nine most abundant elements H, He, C, N, O, Ne, Mg, Si, and Fe are considered. It is supplemented by deuterium, D, which is believed to play a role in primordial or low-metallicity contexts. The total set of species consists of all ionization states of these elements, the negatively charged species H⁻, C⁻, and O⁻, and the molecules H₂, ${\text{H}}_2^{+}$, ${\text{H}}_3^{+}$, HD, C₂, O₂, ${\text{O}}_2^{+}$, CH, CH⁺, CH₂, ${\text{CH}}_2^{+}$, ${\text{CH}}_3^{+}$, OH, OH⁺, CO, CO⁺, H₂O, H₂O⁺, H₃O⁺, HCO⁺, and HOC⁺. The chemical network has a total of 426 reactions, which is roughly half the size of the RSO14 network with 907 reactions. As already mentioned, the saving is due to the neglect of S and Ca and because of the absence of photochemical reactions. The NC^​2M includes gas-phase (atomic, ionic, molecular) reactions, cosmic ray reactions, and dust-assisted reactions. A compilation of all chemical reactions with references to sources for the rate coefficients is given in the appendix. The NC^​2M within its assessed limitations represents cooling gases with some fidelity for temperatures T ≳ 50 K and up to number densities of n ≈ 10¹⁰ m⁻³, say. These thresholds in density and temperature are not to be understood as sharp values but rather as plausible estimates that remind us of the fact that at even higher density and lower temperatures, more complex molecules may be formed contributing to cooling, but which are not covered by the NC^​2M. Despite such an upper density limit of validity, the NC^​2M optionally includes higher density cooling processes for the important H₂ molecule, namely, a correction for optically thick line cooling and collision-induced emission cooling. The neglect of photoreactions in the NC^​2M would, strictly speaking, also require us to set a lower density limit above which the assumption of shielded gas is reasonable.

Technically, the NC^​2M is applicable (and will be applied) for any temperature and density. In doing so, however, one has to ensure that the implementation does not suffer from unphysical behavior or unboundedness of the rate coefficients in the asymptotic limits T, n → 0, ∞. All rate coefficients are checked in this respect and, if necessary, are reasonably modified. An example is the rate coefficient from Kingdon & Ferland (1996) for the charge transfer reaction C⁺ + H → C + H⁺ given by

$$\begin{array}{c}\hfill k=4.88\times {10}^{-22}\frac{{\mathrm{m}}^3}{\mathrm{s}}{\left(\frac{T}{{10}^4\mathrm{K}}\right)}^{3.25}[1-1.12exp(-\frac{0.21T}{{10}^4\mathrm{K}})],\end{array}$$(2)

with specified temperature range [5.5 × 10³ K, 10⁵ K]. Below T ≈ 5390 K the rate becomes negative, which is unphysical. Therefore, k is set to zero below this threshold. Likewise, an upper floor is set at T = 10⁹ K. Although k becomes inaccurate for T ≫ 10⁵ K, it will not have a large impact on the overall chemical evolution since the total reaction rate, kn_C⁺n_H, for this process is negligible at such high temperatures owing to small number densities n_C⁺ and n_H. Generally, reaction rate coefficients are often not very precisely known. Uncertainties of 100% or even larger sometimes appear in the literature, which reduces confidence in the outcome of any reaction network (see Wakelam et al. 2005).

2.1. Gas-phase reactions

2.1.1. Ionization, recombination, and charge transfer

In calculating the ionization state of the gas, the chemical network of the NC^​2M includes the processes of collisional ionization (CI),

$$\begin{array}{c}\hfill {\mathrm{X}}^{n+}+{\mathrm{e}}^{-}\to {\mathrm{X}}^{(n+1)+}+2{\mathrm{e}}^{-},\end{array}$$

radiative recombination (RR),

$$\begin{array}{c}\hfill {\mathrm{X}}^{(n+1)+}+{\mathrm{e}}^{-}\to {\mathrm{X}}^{n+}(+h\nu ),\end{array}$$

and di-electronic recombination (DR),

$$\begin{array}{c}\hfill {\mathrm{X}}^{(n+1)+}+{\mathrm{e}}^{-}\to {\mathrm{X}}^{n+,\ast }\to {\mathrm{X}}^{n+}(+h\nu ),\end{array}$$

for all elements X = H, D, He, C, N, O, Ne, Mg, Si, and Fe. In the DR process, X^* denotes an intermediate excited state. Case A recombination is assumed for hydrogen and helium. CI, RR, and DR rate coefficients for the lighter elements H, D, and He are taken from various sources as specified in Table A.1, reactions 1–6, 350, 351. CI rate coefficients for metal species are described by the analytical fits in Voronov (1997), which encompass a large temperature range. RR and DR coefficients for almost all metal species rely on web-published¹ atomic data and fitting functions based on detailed calculations with the code AUTOSTRUCTURE (Badnell 2006). Unfortunately, for higher isoelectronic sequences, which concerns all the ionization states of iron below Fe¹³⁺, no data is available. Therefore, for those ions the older RR coefficients from Shull & van Steenberg (1982) are used, and DR coefficients are taken from Arnaud & Raymond (1992).

The network also covers charge transfer (CT) interactions of various species with hydrogen and helium. Either a neutral H (or He) atom transfers an electron to a species reducing its ionization state and H (or He) becomes ionized or, inversely, an electron is transferred from a species, increasing its ionization state, to a H⁺ (or He⁺) ion, which becomes neutral. CT reactions between metal species are neglected because they are of minor importance compared to H and He CT processes. Contrary to CI, RR, and DR processes, CT is generally regarded as a secondary effect leading to only moderate changes in the ionization balance. The four different CT reaction types are

$$\begin{array}{c}\hfill \mathrm{X}+{\mathrm{H}}^{+}\to {\mathrm{X}}^{+}+\mathrm{H},\end{array}$$

for hydrogen CT ionization of concerned species X = D, He, C, N, O, Mg, Mg⁺, Si, Si⁺, Fe, Fe⁺, H₂, O₂, OH, and H₂O;

$$\begin{array}{c}\hfill \mathrm{X}+{\mathrm{He}}^{+}\to {\mathrm{X}}^{+}+\mathrm{He},\end{array}$$

for helium CT ionization of concerned species X = D, Si, H₂, O₂, OH, and H₂O;

$$\begin{array}{c}\hfill {\mathrm{X}}^{+}+\mathrm{H}\to \mathrm{X}+{\mathrm{H}}^{+},\end{array}$$

for hydrogen CT recombination of concerned species X = D, He, C, C⁺, C²⁺, C³⁺, N, N⁺, N²⁺, N³⁺, O, O⁺, O²⁺, O³⁺, Ne⁺, Ne²⁺, Ne³⁺, Mg⁺, Mg²⁺, Mg³⁺, Si⁺, Si²⁺, Si³⁺, Fe⁺, Fe²⁺, Fe³⁺, H₂, and CO; and

$$\begin{array}{c}\hfill {\mathrm{X}}^{+}+\mathrm{He}\to \mathrm{X}+{\mathrm{He}}^{+},\end{array}$$

for helium CT recombination of concerned species X = D, C²⁺, C³⁺, N⁺, N²⁺, N³⁺, O⁺, O²⁺, O³⁺, Ne⁺, Ne²⁺, Ne³⁺, Mg²⁺, Mg³⁺, Si²⁺, Si³⁺, Fe²⁺, and Fe³⁺. CT rate coefficients come from different sources as specified in Table A.1. In particular, rates for the reactions involving lower ionized metal species are taken from the Stancil et al. (2015) web database², which catalogs data from previous works (e.g., Kingdon & Ferland 1996).

Atomic data and fitting functions for metal species in the NC^​2M are probably very similar, and in parts identical, to those used in other investigations. For example, the ionization model of Oppenheimer & Schaye (2013a) and the code Cloudy version 10 both apply the same Voronov (1997) fits for the CI rate coefficients and the same Badnell (2006) fits for RR and DR rate coefficients. Moreover, CT rates for many reactions in the NC^​2M with non-molecular species are likely identical to those in Oppenheimer & Schaye (2013a) since coefficients originate more or less from the same literature sources. Based on data common ground, the equilibrium ionization structure of metals in the high-temperature regime is therefore expected to be very similar. The ionization structure in the low-temperature regime, on the other hand, is more difficult since many processes involving molecules play a role, and calculations depend much more on model details.

2.1.2. H₂ chemistry

Molecular hydrogen is the most abundant molecule in the ISM, and it controls much of the chemistry running at a lower gas temperature. Large H₂ fractions are observed in dense and dusty regions where its formation is dictated by reactions on the surface of interstellar dust grains, which act as catalysts. Such dust-surface-catalytic reactions are considered in more detail in Sect. 2.3. In the absence of dust, generally smaller H₂ fractions are produced by pure gas-phase reactions.

Since the two-body association process 2H → H₂ is excluded, there remain two principle H₂ formation channels in the gas phase. The first pathway is initiated by the H⁻ ion, which is produced by the radiative association reaction

$$\begin{array}{c}\hfill \mathrm{H}+{\mathrm{e}}^{-}\to {\mathrm{H}}^{-}(+h\nu ).\end{array}$$

H₂ is then formed by associative detachment with H according to

$$\begin{array}{c}\hfill {\mathrm{H}}^{-}+\mathrm{H}\to {\mathrm{H}}_2+{\mathrm{e}}^{-}.\end{array}$$

This two-step H₂ formation sequence is probably the most important gas-phase mechanism in warm interstellar gas with sufficient fractional ionization. However, the sequence is a relatively slow process hampered by the usually low abundance of H⁻ and in general cannot compete with H₂ grain catalysis. In addition to the above fundamental reactions, a couple of further reactions are taken into account involving also the ions ${\text{H}}_2^{+}$ and ${\text{H}}_3^{+}$. Some of the reactions contribute to H₂ formation like charge transfer from H to ${\text{H}}_2^{+}$ (reaction 303) or the dissociative detachment of ${\text{H}}_3^{+}$ with free electrons and H (reactions 323, 325). Other reactions are counterproductive, like the mutual neutralization of H⁻ with H⁺ (reaction 309), which competes with associative detachment above, or charge transfer from H₂ to H⁺ (reaction 304).

At sufficiently high gas densities (such as n ≳ 10¹⁶ m⁻³), a second pathway for H₂ formation is by ter-molecular reactions. Here, the following three important reactions are taken into account:

$$\begin{array}{c}\hfill 2\mathrm{H}+\mathrm{X}\to {\mathrm{H}}_2+\mathrm{X},\end{array}$$

where X = H, He, H₂.

The H₂ subnetwork in the NC^​2M has been compiled from the chemical networks in Glover & Abel (2008) and GFM10. Reaction rates are adopted from those works except for some reactions that apply more recent rates: reaction 301 uses the rate from Bruhns et al. (2010), reaction 305 uses the rate from Trevisan & Tennyson (2002), and reaction 316 uses the UMIST RATE12 value (McElroy et al. 2013). RSO14 employs a similar set of reactions to treat gas-phase H₂ formation and destruction processes, except that reaction 328 of the NC^​2M, describing detachment of H⁻ with He, is missing in RSO14.

2.1.3. HD chemistry

The hydrogen reaction family is supplemented by a mini-network for HD formation as suggested by Glover & Jappsen (2007). Once molecular hydrogen is present it has been argued that the proton interchange reactions

$$\begin{array}{c}\hfill {\mathrm{H}}_2+{\mathrm{D}}^{+}\rightleftarrows \mathrm{HD}+{\mathrm{H}}^{+}\end{array}$$

are the most significant reactions in regulating the amount of HD, if the gas is sufficiently ionized. The importance of this reaction pair stems from the fact that below a few hundred Kelvin chemical fractionation can occur meaning a substantial enhancement of the HD:H₂ abundance ratio relative to the prevalent D:H ratio. At higher temperatures, the analog neutral-neutral reactions

$$\begin{array}{c}\hfill {\mathrm{H}}_2+\mathrm{D}\rightleftarrows \mathrm{HD}+\mathrm{H}\end{array}$$

may also influence the HD abundance. These two pairs of principle reactions together with ionization, recombination, and CT reactions for D build the subnetwork to handle deuterium chemistry in the NC^​2M. This network obviously neglects many D-bearing analogs to H-bearing species like molecular deuterium, D₂, deuterated hydronium, D₃O⁺, or the ions D⁻ and HD⁺. Reactions of D⁻ and HD⁺ with hydrogen, in particular, would provide an additional, albeit minor, source of HD formation analog to the H⁻/${\text{H}}_2^{+}$ channel of H₂ formation. Such a more detailed treatment of deuterium chemistry, which is not an intended goal here, can be found in Rodgers & Millar (1996), for example.

2.1.4. Carbon and oxygen chemistry

The gas-phase carbon and oxygen chemistry in the NC^​2M is based on the RSO14 compilation of reactions. The RSO14 compilation, in turn, is built upon the network presented in GFM10. RSO14 added three reactions to the GFM10 model, which they claimed to be of significance for the CO formation process discovered by comparison with the Cloudy software. The NC^​2M accounts for these reactions (483, 510, and 511 in Table A.1) using the rates given in RSO14. Other reaction rates are adopted from the GFM10 compilation with some exceptions: reactions 400, 404, 419, 437, 449, 454, 477, 479, 480-482, 488, 489, 491, 492, and 515 use updated rates from the UMIST RATE12 compilation. Reaction 493 applies the rate given in Martinez et al. (2008). The dissociative electron recombination reactions 498–500 with hydronium make use of novel rates from Novotny et al. (2010). Finally, reactions 504 and 505 use more recent rates from Klippenstein & Georgievskii (2010).

The applied network has been designed to accurately follow the abundances of the major molecular coolants CO and H₂O. It includes up to four-atom molecular species involved in their formation and destruction processes. Nevertheless, the network represents a considerable simplification with respect to comprehensive astrochemistry databases like UMIST. However, GFM10 demonstrated that this network is able to adequately model the abundances of C-bearing and O-bearing species when compared to a more extensive network derived from the UMIST RATE99 database, if equal reaction rates are assumed. Moreover, Glover & Clark (2012) have examined several approximate networks with different levels of detail for modeling CO chemistry among which the GFM10 network was the most detailed.

The complexity of chemical processes in diffuse and dense interstellar gas has been discussed in Herbst & Klemperer (1973) and Dalgarno & Black (1976) including C-,N-,O-,S-,Cl-,Si-bearing molecules. The network in the NC^​2M circumscribes the dominant chemical cycles for C-bearing and O-bearing molecules, ignoring chemical schemes of the other elements. For the formation of CO and H₂O in H₂-dominated regions, intermediate molecules and its ionized versions such as OH, OH⁺, the hydrocarbons CH_x, ${\text{CH}}_x^{+}$, and the ions HCO⁺ and HOC⁺ play a role. The major route to form H₂O is initiated by the reactions

$$\begin{array}{cc}& {\mathrm{O}}^{+}+{\mathrm{H}}_2\to {\mathrm{OH}}^{+}+\mathrm{H}\hfill \\ \hfill & \mathrm{O}+{\mathrm{H}}_3^{+}\to {\mathrm{OH}}^{+}+{\mathrm{H}}_2\hfill \end{array}$$

forming the hydroxyl ion. This is followed by the sequence

$$\begin{array}{cc}& {\mathrm{OH}}^{+}+{\mathrm{H}}_2\to {\mathrm{H}}_2{\mathrm{O}}^{+}+\mathrm{H}\hfill \\ \hfill & {\mathrm{H}}_2{\mathrm{O}}^{+}+{\mathrm{H}}_2\to {\mathrm{H}}_3{\mathrm{O}}^{+}+\mathrm{H}\hfill \end{array}$$

and subsequent dissociative recombinations

$$\begin{array}{cc}\hfill {\mathrm{OH}}^{+}+{\mathrm{e}}^{-}& \to \mathrm{O}+\mathrm{H}\hfill \\ \hfill {\mathrm{H}}_2{\mathrm{O}}^{+}+{\mathrm{e}}^{-}& \to \mathrm{OH}+\mathrm{H}\hfill \\ \hfill & \to \mathrm{O}+{\mathrm{H}}_2\hfill \\ \hfill & \to \mathrm{O}+2\mathrm{H}\hfill \\ \hfill {\mathrm{H}}_3{\mathrm{O}}^{+}+{\mathrm{e}}^{-}& \to {\mathrm{H}}_2\mathrm{O}+\mathrm{H}\hfill \\ \hfill & \to \mathrm{OH}+{\mathrm{H}}_2\hfill \\ \hfill & \to \mathrm{OH}+2\mathrm{H}\hfill \\ \hfill & \to \mathrm{O}+{\mathrm{H}}_2+\mathrm{H}.\hfill \end{array}$$

Where the level of ionization is low, the following neutral-neutral reactions contribute to H₂O formation:

$$\begin{array}{cc}& \mathrm{OH}+\mathrm{H}\to {\mathrm{H}}_2\mathrm{O}\hfill \\ \hfill & \mathrm{OH}+{\mathrm{H}}_2\to {\mathrm{H}}_2\mathrm{O}+\mathrm{H}\hfill \\ \hfill & 2\mathrm{OH}\to {\mathrm{H}}_2\mathrm{O}+\mathrm{O}.\hfill \end{array}$$

Important routes to form CO are by the reactions

$$\begin{array}{cc}\hfill {\mathrm{C}}^{+}+\mathrm{OH}& \to \mathrm{CO}+{\mathrm{H}}^{+}\hfill \\ \hfill & \to {\mathrm{CO}}^{+}+\mathrm{H},\hfill \end{array}$$

followed by

$$\begin{array}{cc}& {\mathrm{CO}}^{+}+{\mathrm{H}}_2\to {\mathrm{HCO}}^{+}+\mathrm{H}\hfill \\ \hfill & {\mathrm{HCO}}^{+}+{\mathrm{e}}^{-}\to \mathrm{CO}+\mathrm{H},\hfill \end{array}$$

and by the reactions

$$\begin{array}{cc}& \mathrm{O}+{\mathrm{CH}}^{+}\to \mathrm{CO}+{\mathrm{H}}^{+}\hfill \\ \hfill & \mathrm{O}+{\mathrm{C}}_2^{+}\to \mathrm{CO}+{\mathrm{C}}^{+},\hfill \end{array}$$

where the formation and destruction of hydrocarbon radicals up to ${\text{CH}}_3^{+}$ is governed by numerous reactions not discussed here, but fully included in the NC^​2M. Further minor sources of CO are the dissociative recombination

$$\begin{array}{c}\hfill {\mathrm{HOC}}^{+}+{\mathrm{e}}^{-}\to \mathrm{CO}+\mathrm{H},\end{array}$$

the ion-neutral reactions

$$\begin{array}{cc}& {\mathrm{HCO}}^{+}+\mathrm{C}\to \mathrm{CO}+{\mathrm{CH}}^{+},\hfill \\ \hfill & {\mathrm{HCO}}^{+}+{\mathrm{H}}_2\mathrm{O}\to \mathrm{CO}+{\mathrm{H}}_3{\mathrm{O}}^{+},\hfill \end{array}$$

and the charge transfer reaction

$$\begin{array}{c}\hfill {\mathrm{CO}}^{+}+\mathrm{H}\to \mathrm{CO}+{\mathrm{H}}^{+}.\end{array}$$

In very dense regions, three-body molecular associations with neutral helium acting as catalyst may become important. The NC^​2M includes a couple of reactions of this type that are numbered 519–525 in Table A.1.

2.2. Cosmic ray reactions

Cosmic rays (CR) are such highly energetic particles that they can deeply penetrate dense gas regions that are opaque to ionizing ultraviolet (UV) radiation otherwise. CR-induced ionizations of H, H₂, and He maintain a sufficient level of free electrons and so act as a driver for ion-molecule chemistry at low gas temperatures. The NC^​2M takes into account two distinct types of cosmic ray reactions: primary and secondary. Primary reactions (labeled “cr”-reactions) denote the direct impact of a cosmic ray particle resulting in the ionization of a specie X or the dissociation of H₂,

$$\begin{array}{c}\hfill \mathrm{X}\stackrel{\mathrm{cr}}{\to }{\mathrm{X}}^{+}+{\mathrm{e}}^{-},\end{array}$$

$$\begin{array}{c}\hfill {\mathrm{H}}_2\stackrel{\mathrm{cr}}{\to }2\mathrm{H}.\end{array}$$

Secondary reactions (labeled “γ_cr”-reactions) denote photoreactions due to ultraviolet photons generated internally by the interaction of energetic electrons ejected into primary cosmic ray ionizations with the gas. Secondary reactions are the photoionization of species X and the photodissociation of molecules XY,

$$\begin{array}{c}\hfill \mathrm{X}\stackrel{{\gamma }_{\mathrm{cr}}}{\to }{\mathrm{X}}^{+}+{\mathrm{e}}^{-},\end{array}$$

$$\begin{array}{c}\hfill \mathrm{XY}\stackrel{{\gamma }_{\mathrm{cr}}}{\to }\mathrm{X}+\mathrm{Y}.\end{array}$$

A summary of all CR-induced reactions included in the NC​²M is given in Table A.2. The subset of cosmic ray reactions is nearly equivalent to RSO14 but the rate coefficients differ for some reactions because of the dissimilar handling of secondary effects described below.

The primary cosmic ray ionization rate, ζ, is dominated by the energy spectrum of cosmic ray particles in the lower energy regime (≲100 MeV). Unfortunately, this part of the energy spectrum is observationally not well determined. For the cosmic ray ionization rate of hydrogen, the value of ζ_H = 2.5 × 10⁻¹⁷ s⁻¹ from Williams et al. 1998 is adopted in order to allow comparison with other work. However, more recent observations of ${\text{H}}_3^{+}$ in the Galactic diffuse ISM indicate that the rate may be higher by one order of magnitude (Indriolo & McCall 2012). The ionization rates for other species are assumed to be linearly correlated,

$$\begin{array}{c}\hfill \zeta ={\zeta }^{\mathrm{rel}}{\zeta }_{\mathrm{H}}.\end{array}$$(3)

The numbers ζ^rel are extracted from the UMIST RATE12 database insofar available. For all other “cr”-reactions, the factors are estimated on the basis of the theoretical considerations in Silk (1970) and Lotz (1967). Namely,

$$\begin{array}{c}\hfill {\zeta }^{\mathrm{rel}}=\overline{\xi }\frac{{\chi }_{\mathrm{H}}}{\chi },\end{array}$$(4)

where χ_H and χ are the ionization energies of hydrogen and the species in mind, respectively, and $\overline{\xi }$ denotes the outer shell effective electron number of that species given by

$$\begin{array}{c}\hfill \overline{\xi }=\overline{\chi }\sum _{j=1}^2\frac{{\xi }_j}{{\chi }_j},ln\overline{\chi }=\frac{{\sum }_{j=1}^2{\xi }_j/{\chi }_{j}ln{\chi }_j}{{\sum }_{j=1}^2{\xi }_j/{\chi }_j}·\end{array}$$(5)

Here, ξ_j and χ_j are the number of electrons and ionization potentials of the outermost (j = 1) and next inner (j = 2) subshell of the species. The various values can be found in Table 1 in Lotz (1967).

The ejected electron in a primary ionization has an average energy of about 30 eV for cosmic ray particles in the energy range 10…100 MeV (Cravens & Dalgarno 1978). Such energetic electrons interacting with the gas can give rise to secondary effects, which are of potential significance especially in low-temperature, high-density environments: further ionization, collisional excitation of atoms and molecules, dissociation of molecules, and Coloumb scattering in partly ionized gas. Dalgarno et al. (1999) have made detailed investigations of the energy degradation process of energetic electrons in a partly ionized H₂–H–He gas mixture³ for incident electron energies up to 1 keV. To account for such secondary effects in the NC^​2M, their fitting formulae for a 30 eV electron have been elaborated. This way the number of secondary ionizations of H and H₂, the yields of H(2p) excitations (Lyα photons), and the yields of electronic excitations of H₂ (Lyman–Werner band photons) are computed as functions of the number density ratio n_H2/n_H and ionization degree. The net (primary plus secondary) cosmic ray ionization rate for atomic and molecular hydrogen is therefore represented by

$$\begin{array}{cc}\hfill {\zeta }_{\mathrm{H}}^{\mathrm{tot}}& =(1+N_{{\mathrm{H}}^{+}}^{\sec }){\zeta }_{\mathrm{H}},\hfill \end{array}$$(6)

$$\begin{array}{cc}\hfill {\zeta }_{{\mathrm{H}}_2}^{\mathrm{tot}}& =(1+N_{{\mathrm{H}}_2^{+}}^{\sec }){\zeta }_{{\mathrm{H}}_2},\hfill \end{array}$$(7)

where the numbers of secondary ionizations are

$$\begin{array}{cc}\hfill N_{{\mathrm{H}}^{+}}^{\sec }& \approx 0.5x_{\mathrm{H}}{[1+221{\left(\frac{n_{\mathrm{e}}}{0.2n_{{\mathrm{H}}_2}+1.1n_{\mathrm{H}}}\right)}^{1.06}]}^{-1},\hfill \end{array}$$(8)

$$\begin{array}{cc}\hfill N_{{\mathrm{H}}_2^{+}}^{\sec }& \approx 0.76x_{{\mathrm{H}}_2}{[1+202{\left(\frac{n_{\mathrm{e}}}{1.2n_{{\mathrm{H}}_2}+0.1n_{\mathrm{H}}}\right)}^{1.15}]}^{-1},\hfill \end{array}$$(9)

with x_H = n_H/(2n_H2 + n_H) and x_H2 = n_H2/(2n_H2 + n_H). Dalgarno et al. (1999) state that the expressions for $N_{{\mathrm{H}}^{+}}^{\sec }$ and $N_{{\mathrm{H}}_2^{+}}^{\sec }$ are accurate for ionization degrees up to 10%. For a neutral H–He gas mixture (x_H = 1, n_e = 0) $N_{{\mathrm{H}}^{+}}^{\sec }\approx 0.5$, whereas for a neutral H₂-He gas mixture (x_H2 = 0.5, n_e = 0) $N_{{\mathrm{H}}_2^{+}}^{\sec }\approx 0.38$. The number of secondary ionizations decreases with increasing ionization degree because a larger fraction of the ejected electron energy goes into gas heating via Coloumb scattering lowering the yields of ionizations and excitations. At a 10% ionization level the number of secondary ionizations is reduced by more than one order of magnitude compared to neutral gas. Secondary ionizations of He (and also metals) are ignored because their contribution to the ionization rate is small. In the case of the above neutral H–He gas mixture, we find $N_{{\mathrm{He}}^{+}}^{\sec }\approx 5\times {10}^{-3}$, and $N_{{\mathrm{He}}^{++}}^{\sec }=0$ since the necessary ionization energy of 54.4 eV is above the ejected electron energy of 30 eV. In RSO14 the hydrogen secondary ionization rate is computed from data tables in Furlanetto & Stoever (2010) obtained by Monte-Carlo simulations of the energy degradation problem. However, these simulations were done for atomic gas and are therefore not suited to estimate secondary ionization rates in a H–H₂ gas mixture.

Gredel et al. (1987, 1989) have studied CR-induced photoionization and photodissociation for a variety of interstellar molecules. The authors determined photoreaction rates on the basis of detailed calculations of the emission spectrum resulting from the impact excitation of H₂ by a 30 eV electron. Accordingly, the photoreaction rate R_X in units s⁻¹ for a species X is given by

$$\begin{array}{c}\hfill R_{\mathrm{X}}=\frac{p_{\mathrm{X}}}{1-\omega }{\zeta }_{{\mathrm{H}}_2}^{\mathrm{tot}},\end{array}$$(10)

where p_X is the efficiency for reactions with Lyman–Werner photons per total number of ionizations and ω is the grain albedo (ω = 0.5 is used throughout). In computing p_X a pure neutral H₂ gas with $N_{{\mathrm{H}}_2^{+}}^{\sec }=0.55$ was assumed, and yields of 0.3 and 0.12 for the dominating excitations to the $B^1{\mathrm{\Sigma }}_{\text{u}}^{+}$ (including cascading contribution from higher excited states) and C¹Π_u electronic states of H₂ were adopted. For the NC^​2M, the Gredel et al. outcome has been recalibrated to the situation of a partially ionized H–H₂–He gas mixture making use of the energy degradation results in Dalgarno et al. (1999) a second time. If atomic hydrogen is present, excitations to the H(2p) state followed by Lyα emission contribute to the generated UV flux. From the fitting formulae in Dalgarno et al. (1999), the yields for Lyman–Werner photons and Lyα photons are approximately given by

$$\begin{array}{c}\hfill Y_{\mathrm{LyW}}\approx 1.09N_{{\mathrm{H}}_2^{+}}^{\sec },Y_{\mathrm{Ly}\alpha }\approx 1.8N_{{\mathrm{H}}^{+}}^{\sec }.\end{array}$$(11)

This suggest the following modified photoreaction rate:

$$\begin{array}{c}\hfill R_{\mathrm{X}}=(\frac{1.55}{0.42}{Y}_{\mathrm{LyW}}{p}_{\mathrm{X}}+Y_{\mathrm{Ly}\alpha }{q}_{\mathrm{X}}){(1-\omega )}^{-1}\overline{\zeta },\end{array}$$(12)

where q_X is the efficiency of a species X for reactions with Lyα photons. Values for various species can be found in Maloney et al. (1996). The averaged primary ionization rate, $\overline{\zeta }$, is given by

$$\begin{array}{c}\hfill \overline{\zeta }=\frac{n_{\mathrm{H}}{\zeta }_{\mathrm{H}}+n_{{\mathrm{H}}_2}{\zeta }_{{\mathrm{H}}_2}}{n_{\mathrm{H}}+n_{{\mathrm{H}}_2}}·\end{array}$$(13)

We note that R_X now depends on the ionization degree because the yields Y_LyW and Y_Lyα do so. For a neutral H₂-He gas mixture (x_H2 = 0.5, n_He = 0.2n_H2, n_e = 0), the rate reduces to

$$\begin{array}{c}\hfill R_X=\frac{p_{\mathrm{X}}}{1-\omega }1.53{\zeta }_{{\mathrm{H}}_2},\end{array}$$(14)

which is close to the Gredel et al. (1989) result where ${\zeta }_{{\mathrm{H}}_2}^{\mathrm{tot}}=1.55{\zeta }_{{\mathrm{H}}_2}$. The efficiencies p_X and q_X for all CR-induced photoreactions in the NC^​2M are summarized in Table A.2.

As discussed in Gredel et al. (1987), the efficiency p_CO of CO for reactions with Lyman–Werner photons must be treated more carefully since efficient destruction of CO occurs only if close coincidences between H₂ emission lines and CO absorption lines exist. The efficiency p_CO is tabulated in Gredel et al. (1987) as a function of temperature, CO abundance x_CO, and grain albedo. For a grain albedo of 0.5 a crude interpolation is p_CO ≈ 1.35(T/x_CO)^0.5. The CR-induced CO photodissociation rate is then approximated by

$$\begin{array}{c}\hfill R_{\mathrm{CO}}=\frac{1.55}{0.42}{Y}_{\mathrm{LyW}}{p}_{\mathrm{CO}}\overline{\zeta }.\end{array}$$(15)

This approaches the original Gredel et al. (1987) result, R_CO = p_CO1.55ζ_H2, if electron degradation is assumed to take place in a neutral H₂ gas where Y_LyW ≈ 0.41 according to Dalgarno et al. (1999).

2.3. Dust-catalytic reactions

To understand the high fraction of molecular hydrogen observed in various interstellar environments, dust grain chemistry seems indispensable. The NC​²M adheres to the fundamental process of H₂ formation on the surface of dust particles. Moreover, several relevant grain-surface ion recombination reactions are considered. A compilation of all dust-catalytic reactions is given in Table A.3. The efficacy of surface reactions depends on the total cross section of the dust mixture. A static dust model is predisposed, that is, the current NC^​2M implementation does not account for dust destruction (sputtering, sublimation) and formation (accretion, coagulation) phenomena. Therefore, the total dust cross section does not change with time. Also, the implementation of ion-dust recombination assumes that neutralized ions re-enter the gas phase after charge transfer and are not adsorbed onto the grain. Treatment of such element depletion processes is currently beyond the capabilities of the NC^​2M.

The NC^​2M applies the dust model of Weingartner & Draine (2001a), which reproduces the observed extinction curve and dust emission spectrum of the Milky Way by a mixture of graphite, PAHs, and silicate spherical grains with distinct size distributions n_d(a) (a: grain radius). According to Draine (2002) the total dust cross-sectional area per unit volume for this model is

$$\begin{array}{c}\hfill {(n_{\mathrm{d}}{\sigma }_{\mathrm{d}})}_{\mathrm{tot}}={\displaystyle \int \mathrm{d}a\frac{\mathrm{d}{n}_{\mathrm{d}}}{\mathrm{d}a}\pi a^2=6.7\times {10}^{-25}{\mathrm{m}}^2·n_{{\mathrm{H}}_{\mathrm{nuc}}},}\end{array}$$(16)

with n_Hnuc being the number density of H nuclei⁴. For a dust-to-gas ratio other than the average value in the Milky Way, it is common to scale (n_dσ_d)_tot linearly with the metallicity Z, hence, (n_dσ_d)_tot = 6.7 × 10⁻²⁵m² ⋅ n_Hnuc(Z/Z_⊙).

The rate of H₂ formation according to reaction 700, ${\dot{n}}_{{\mathrm{H}}_2}=R_{{\mathrm{H}}_2}=k_{700}{n}_{\mathrm{H}}^2$, can be expressed as (see, e.g., Cazaux & Tielens 2002)

$$\begin{array}{c}\hfill R_{{\mathrm{H}}_2}=\frac{1}{2}{n}_{\mathrm{H}}{\upsilon }_{\mathrm{th},\mathrm{H}}{(n_{\mathrm{d}}{\sigma }_{\mathrm{d}})}_{\mathrm{tot}}{S}_{\mathrm{H}}{ϵ}_{{\mathrm{H}}_2},\end{array}$$(17)

where

$$\begin{array}{c}\hfill {\upsilon }_{\mathrm{th},\mathrm{H}}=\sqrt{\frac{8kT}{\pi m_{\mathrm{H}}}}\end{array}$$(18)

is the mean thermal speed of the H gas component and S_H(T, T_d) is the sticking probability for adsorbing H atoms, which generally depends on both the gas temperature T and the dust temperature T_d. The sticking coefficient is taken from Hollenbach & McKee (1979)

$$\begin{array}{c}\hfill S_{\mathrm{H}}(T,T_{\mathrm{d}})={[1+0.04\sqrt{\frac{T+T_{\mathrm{d}}}{\mathrm{K}}}+0.002\frac{T}{\mathrm{K}}+8\times {10}^{-6}{\left(\frac{T}{\mathrm{K}}\right)}^2]}^{-1}.\end{array}$$(19)

The function ϵ_H2 describes the surface recombination efficiency, which is the fraction of accreted H atoms that leave the surface as H₂. The NC^​2M adopts the general expression in Cazaux & Tielens (2002), Eq. (15). This expression follows from the rate equation model of Cazaux & Tielens (2002, 2004) describing grain surface abundances in physisorbed and chemisorbed sites subject to different surface processes (desorption, migration: thermal diffusion and quantum mechanical tunneling, reaction). Using values for the surface characteristics parameters between those for silicate and carbonaceous grain materials as cataloged in Cazaux & Spaans (2004), one finds

$$\begin{array}{cc}\hfill {ϵ}_{{\mathrm{H}}_2}& ={[1+7\times {10}^{-14}\mathrm{s}Fexp\left(\frac{440\mathrm{K}}{T_{\mathrm{d}}}\right)+19exp(-\frac{225\mathrm{K}}{T_{\mathrm{d}}})]}^{-1}\hfill \\ \hfill & \times {[1+\frac{5.7\times {10}^{14}{\mathrm{s}}^{-1}}{F}exp(-\frac{45000\mathrm{K}}{T_{\mathrm{d}}})]}^{-1},\hfill \end{array}$$(20)

where F is the incoming H atom flux expressed in monolayers per second. This flux is approximated by

$$\begin{array}{c}\hfill F=\frac{1}{4}{\upsilon }_{\mathrm{th},\mathrm{H}}{n}_{\mathrm{H}}4\pi a^2{S}_{\mathrm{H}}\frac{1}{N_{\mathrm{ass}}}=\frac{1}{4}{\upsilon }_{\mathrm{th},\mathrm{H}}{n}_{\mathrm{H}}{S}_{\mathrm{H}}{W}^2,\end{array}$$(21)

where v_th, Hn_Hπa²S_H is the number of H atoms per second landing on a spherical grain with radius a and N_ass = 4πa²/W² denotes the number of adsorption sites on the surface. A typical value for the horizontal barrier width between adjacent adsorption sites is W ≈ 2Å (Cazaux & Tielens 2004). For an order of magnitude estimate in cold cores, for instance, adopting n = 10¹⁰m⁻³ and T = T_d = 10 K gives a value of F ≈ 5 × 10⁻⁸ monolayers per second. We note that R_H2 becomes negligible for higher gas temperatures since v_th, HS_H ∝ T^−1.5 for T ≫ 1000 K.

A second relevant grain-surface process that codetermines the ionization balance in the low-temperature regime is the neutralization of ions. This process competes with radiative recombination for gas temperatures below a few thousand Kelvin and counteracts cosmic ray ionization in cold, UV-shielded regions. The NC^​2M takes into account ion recombination of the elements X = H⁺, He⁺, C⁺, Mg⁺, Si⁺, and Fe⁺. The rate coefficients, k_X, for reactions 701–706 use the fitting formula from Weingartner & Draine (2001b), Eq. (8),

$$\begin{array}{c}\hfill {\dot{n}}_{\mathrm{X}}=-k_{\mathrm{X}}{n}_{\mathrm{X}}\approx -\frac{{10}^{-20}{C}_0{n}_{{\mathrm{H}}_{\mathit{nuc}}}{\mathrm{m}}^3{\mathrm{s}}^{-1}}{1+C_1{\stackrel{\sim }{\mathrm{\Psi }}}^{C_2}(1+C_3{\left(\frac{T}{\mathrm{K}}\right)}^{C_4}{\stackrel{\sim }{\mathrm{\Psi }}}^{-C_5-C_{6}ln(T/\mathrm{K})})}{n}_{\mathrm{X}},\end{array}$$(22)

with the coefficients C₀ − C₆ given in Table 2 of Weingartner & Draine (2001b). This formula has been derived for the Milky Way dust model of Weingartner & Draine (2001a), already used for R_H2, and accounts for grain charging expressed by the (dimensionless) charging parameter (see Weingartner & Draine 2001c)

$$\begin{array}{c}\hfill \stackrel{\sim }{\mathrm{\Psi }}=G\frac{{10}^6{\mathrm{m}}^{-3}}{n_{\mathrm{e}}}\sqrt{\frac{T}{\mathrm{K}}}·\end{array}$$(23)

The quantity G is a measure of the radiation intensity (in units of the Habing radiation field energy density) responsible for photoelectric grain charging. Equation (22) is multiplied by an additional factor Z/Z_⊙ to account for different dust-to-gas ratios as has been done for R_H2. As stressed by Weingartner & Draine (2001b), the expression for k_X likely becomes inaccurate when $\stackrel{\sim }{\mathrm{\Psi }}\lesssim {10}^2$. It is interesting to note that the limit G → 0 ($\stackrel{\sim }{\mathrm{\Psi }}\to 0$) would correspond to a completely radiation-free environment where grain charging is driven solely by electron collisional charging (ion collisional charging was neglected in Weingartner & Draine 2001b) resulting in over-balanced negatively charged grains, and k_X becomes a maximum.

3. Thermal processes

In time-dependent, gas-dynamical simulations thermal processes enter the energy equation as a net volumetric cooling rate Λ_net changing the energy density ϵ,

$$\begin{array}{c}\hfill \frac{\mathrm{d}ϵ}{\mathrm{d}t}=-{\mathrm{\Lambda }}_{\mathrm{net}},\end{array}$$(24)

where Λ_net = Λ − Γ is the difference between total cooling (Λ = ∑Λ_i ≥ 0) and total heating (Γ = ∑Γ_i ≥ 0) rates. The individual cooling and heating contributions, Λ_i and Γ_i, are described in more detail below. The heating processes include chemical heating associated with H₂ formation reactions and cosmic ray heating. Among the cooling processes are chemical cooling due to ionization, recombination, and H₂ dissociation, atomic line cooling of all elements, and molecular line cooling of H₂, CO, H₂O, and OH.

3.1. Chemical heating and cooling

The NC^​2M takes into account chemical heating of the gas through the exothermic reactions of H₂ formation both in the gas phase and dust-assisted. A substantial part of released energy in these reactions goes into the excitation of rotovibrational states of the newly formed hydrogen molecule, which is then partly transformed into gas heating via collisional de-excitation. The energy fraction converted to heat by collisional de-excitation (rather than radiated away) depends on the gas density and is approximated in Hollenbach & McKee (1979) by

$$\begin{array}{c}\hfill f={[1+\frac{{10}^{12}{(T/\mathrm{K})}^{-0.5}{\mathrm{m}}^{-3}}{1.6n_{\mathrm{H}}exp[-{\left(\frac{400\mathrm{K}}{T}\right)}^2]+1.4n_{{\mathrm{H}}_2}exp(-\frac{12000\mathrm{K}}{T+1200\mathrm{K}})}]}^{-1}.\end{array}$$(25)

Following Omukai (2000) and Hollenbach & McKee (1979) the gas heating rates associated with H₂ formation by the different processes, namely, H⁻/${\text{H}}_2^{+}$ (reactions 301,303), three-body reactions 313, 314, 317, and grain-catalytic reaction 700 are given by

$$\begin{array}{cc}\hfill {\mathrm{\Gamma }}_{301}& =(3.53\mathrm{eV})f{k}_{301}{n}_{{\mathrm{H}}^{-}}{n}_{\mathrm{H}}\hfill \end{array}$$(26)

$$\begin{array}{cc}\hfill {\mathrm{\Gamma }}_{303}& =(1.83\mathrm{eV})f{k}_{303}{n}_{{\mathrm{H}}_2^{+}}{n}_{\mathrm{H}}\hfill \end{array}$$(27)

$$\begin{array}{cc}\hfill {\mathrm{\Gamma }}_{313}& =(4.48\mathrm{eV})f{k}_{313}{n}_{\mathrm{H}}^3\hfill \end{array}$$(28)

$$\begin{array}{cc}\hfill {\mathrm{\Gamma }}_{314}& =(4.48\mathrm{eV})f{k}_{314}{n}_{\mathrm{H}}^2{n}_{{\mathrm{H}}_2}\hfill \end{array}$$(29)

$$\begin{array}{cc}\hfill {\mathrm{\Gamma }}_{317}& =(4.48\mathrm{eV})f{k}_{317}{n}_{\mathrm{H}}^2{n}_{\mathrm{He}}\hfill \end{array}$$(30)

$$\begin{array}{cc}\hfill {\mathrm{\Gamma }}_{700}& =(0.2\mathrm{eV}+4.2\mathrm{eV}f)k_{700}{n}_{\mathrm{H}}^2.\hfill \end{array}$$(31)

Reaction 301, 303 liberates 3.53 eV (1.83 eV) of energy, and it is assumed that all go into rotovibrational excitation. The three-body reactions release the H₂ binding energy of 4.48 eV, which is assumed to be stored initially in excitation as well. In the heating rate due to H₂ formation on dust grains, Γ₇₀₀, it is assumed that 4.2 eV of liberated binding energy goes into excitation, 0.2 eV is thermalized kinetic energy of the desorbing molecule, and the remaining amount of energy is absorbed by the grain.

Collisional ionization of H, He, and He⁺ (reactions 1, 3, 5), radiative recombination of H⁺ and He⁺⁺ (reactions 2, 6), radiative and di-electronic recombination of He⁺ (reaction 4), and collisional dissociation of H₂ (reactions 305, 306, 315) are the most relevant chemical cooling processes. The associated cooling rates are compiled from different sources (Spitzer 1978; Janev et al. 1987; Cen 1992; Verner & Ferland 1996; Omukai 2000):

$$\begin{array}{cc}& {\mathrm{\Lambda }}_1=(13.60\mathrm{eV})k_1{n}_{\mathrm{H}}{n}_{\mathrm{e}}\hfill \end{array}$$(32)

$$\begin{array}{cc}& {\mathrm{\Lambda }}_2=\alpha k_{\mathrm{B}}T{k}_2{n}_{{\mathrm{H}}^{+}}{n}_{\mathrm{e}}\hfill \end{array}$$(33)

$$\begin{array}{cc}& {\mathrm{\Lambda }}_3=(24.59\mathrm{eV})k_3{n}_{\mathrm{He}}{n}_{\mathrm{e}}+5.01\times {10}^{-46}\frac{{\mathrm{J}\mathrm{m}}^6}{\mathrm{s}}C(T)\hfill \\ \hfill & \times {(T/\mathrm{K})}^{-0.1687}exp(-55338\mathrm{K}/T)n_{{\mathrm{He}}^{+}}{n}_{\mathrm{e}}^2\hfill \end{array}$$(34)

$$\begin{array}{cc}& {\mathrm{\Lambda }}_4=(\alpha k_{\mathrm{B}}T{k}_{4,RR}+(40.8\mathrm{eV})k_{4,DR})n_{{\mathrm{He}}^{+}}{n}_{\mathrm{e}}\hfill \end{array}$$(35)

$$\begin{array}{cc}& {\mathrm{\Lambda }}_5=(54.42\mathrm{eV})k_5{n}_{{\mathrm{He}}^{+}}{n}_{\mathrm{e}}\hfill \end{array}$$(36)

$$\begin{array}{cc}& {\mathrm{\Lambda }}_6=\alpha k_{\mathrm{B}}T{k}_6{n}_{{\mathrm{He}}^{++}}{n}_{\mathrm{e}}\hfill \end{array}$$(37)

$$\begin{array}{cc}& {\mathrm{\Lambda }}_{305}=(4.48\mathrm{eV})k_{305}{n}_{{\mathrm{H}}_2}{n}_{\mathrm{e}}\hfill \end{array}$$(38)

$$\begin{array}{cc}& {\mathrm{\Lambda }}_{306}=(4.48\mathrm{eV})k_{306}{n}_{{\mathrm{H}}_2}{n}_{\mathrm{H}}\hfill \end{array}$$(39)

$$\begin{array}{cc}& {\mathrm{\Lambda }}_{315}=(4.48\mathrm{eV})k_{315}{n}_{{\mathrm{H}}_2}^2.\hfill \end{array}$$(40)

In the recombination cooling coefficients Λ₂, Λ₄, and Λ₆, αk_BT ≈ 0.75k_BT is the mean kinetic energy of the captured electron⁵. The second term in Λ₃ accounts for collisional ionization cooling from the metastable 2³S₁ He state, which is proportional to $n_{{\text{He}}^{+}}{n}_{\text{e}}^2$ (Black 1981), and where C(T)=(1 + T/10⁵ K)⁻¹.

3.2. Atomic line cooling

3.2.1. H and He

The NC^​2M takes into account radiation cooling from electron impact excitation of H, He, and He⁺. Rates for H and He⁺ are taken from Cen (1992) based on the work of Black (1981). The cooling rate for He⁺ only covers excitation to the n = 2 atomic level, whereas the rate for H covers all n. The rate for neutral He distinguishes between excitation from the ground state 1¹S₀, which is approximated by the function in Glover & Jappsen (2007), Table 7, based on data from Bray et al. (2000), and excitation from the metastable 2³S₁ state to n = 2, 3, 4 triplet states adopting the rate from Cen (1992). The rates can be written as

$$\begin{array}{cc}& {\mathrm{\Lambda }}_{\mathrm{H}}=7.5\times {10}^{-32}\frac{{\mathrm{J}\mathrm{m}}^3}{\mathrm{s}}C(T)exp(-\frac{118348\mathrm{K}}{T})n_{\mathrm{H}}{n}_{\mathrm{e}}\hfill \end{array}$$(41)

$$\begin{array}{cc}& {\mathrm{\Lambda }}_{\mathrm{He}}=1.1\times {10}^{-32}\frac{{\mathrm{J}\mathrm{m}}^3}{\mathrm{s}}{\left(\frac{T}{\mathrm{K}}\right)}^{0.082}exp(-\frac{230000\mathrm{K}}{T})n_{\mathrm{He}}{n}_{\mathrm{e}}\hfill \\ \hfill & +9.1\times {10}^{-46}\frac{{\mathrm{J}\mathrm{m}}^6}{\mathrm{s}}C(T){\left(\frac{T}{\mathrm{K}}\right)}^{-0.1687}exp(-\frac{13179\mathrm{K}}{T})n_{{\mathrm{He}}^{+}}{n}_{\mathrm{e}}^2\hfill \\ \hfill \end{array}$$(42)

$$\begin{array}{cc}& {\mathrm{\Lambda }}_{{\mathrm{He}}^{+}}=5.54\times {10}^{-30}\frac{{\mathrm{J}\mathrm{m}}^3}{\mathrm{s}}C(T){\left(\frac{T}{\mathrm{K}}\right)}^{-0.397}exp(-\frac{473638\mathrm{K}}{T})n_{{\mathrm{He}}^{+}}{n}_{\mathrm{e}}.\hfill \\ \hfill \end{array}$$(43)

3.2.2. Metals

The NC^​2M incorporates fine structure line cooling of the elements C, N, O, Ne, Mg, Si, and Fe in the optically thin approximation. The total cooling due to metals reads

$$\begin{array}{c}\hfill {\mathrm{\Lambda }}_{\mathrm{Z}}={\displaystyle \sum _s{n}_s{\mathrm{\Lambda }}_s,}\end{array}$$(44)

where the sum runs over all species s = C ...C⁵⁺, N ...N⁶⁺, O ...O⁷⁺, Ne ...Ne⁹⁺, Mg ...Mg¹¹⁺, Si ...Si¹³⁺, Fe ...Fe²⁵⁺. The cooling contribution of an individual species or ion is (the index s is suppressed in the following)

$$\begin{array}{c}\hfill \mathrm{\Lambda }={\displaystyle \sum _{i=1}^L{p}_i\sum _{j=1}^{i-1}\mathrm{\Delta }{E}_{\mathit{ij}}{A}_{\mathit{ij}}}\end{array}$$(45)

for transitions i → j of the upper atomic energy level i to the lower atomic energy level j with energy separation ΔE_ij and with the Einstein coefficients A_ij for spontaneous emission. The quantity p_i is the population fraction, that is, the relative number of ions populating level i.

In order to compute the fractions (${(p_i)}_{i=1}^L$ for an ion, the statistical equilibrium equations are solved for a prescribed maximum number of atomic levels L subject to the conservation condition $\sum _{i=1}^L{p}_i=1$. Different values for L are allowed for different ions. It is assumed that excitation processes populating the atomic states can be separated from those regulating the ionization state of an element, that is, that the calculation of the population fractions does not interfere with processes changing the ionization state. The principle of detailed balancing then reveals that the rate of ions leaving an atomic state is equal the rate of ions arriving at that state by (collisional) excitation and de-excitation processes and spontaneous emission (photoexcitation and stimulated emission are absent). Denoting the transition rate i → j as Q_ij (i ≠ j), it follows for any state i that

$$\begin{array}{c}\hfill p_i{\displaystyle \sum _{j=1;j\ne i}^L{Q}_{\mathit{ij}}=\sum _{j=1;j\ne i}^L{p}_j{Q}_{\mathit{ji}},i=1\dots L,}\end{array}$$(46)

or in vector notation, setting $Q_{\mathit{ii}}=-\sum _{j=1;j\ne i}^L{Q}_{\mathit{ij}}$ for convenience,

$$\begin{array}{c}\hfill {\mathrm{Q}}^{\top }\mathbf{p}=0.\end{array}$$(47)

The constraint $\sum _{i=1}^L{p}_i=1$ is directly incorporated by exchanging with it the first equation of Eq. (46). Equation (47) is then replaced by

$$\begin{array}{c}\hfill \stackrel{\sim }{\mathrm{Q}}\mathbf{p}={(1,0,\dots ,0)}^{\top },{\stackrel{\sim }{Q}}_{\mathit{ij}}=\{\begin{array}{cc}1\hfill & i=1\hfill \\ {({\mathrm{Q}}^{\top })}_{\mathit{ij}}\hfill & i\ne 1\hfill \end{array},\end{array}$$(48)

which can be solved for p with standard techniques. The transition rate is given by

$$\begin{array}{c}\hfill Q_{\mathit{ij}}=n_e{Q}_{\mathit{ij}}^{\mathrm{e}}+n_{\mathrm{H}}{Q}_{\mathit{ij}}^{\mathrm{H}}+n_{{\mathrm{H}}_2}{Q}_{\mathit{ij}}^{{\mathrm{H}}_2}+A_{\mathit{ij}},\end{array}$$(49)

where $Q_{\mathit{ij}}^{\mathrm{e}}$ ($Q_{\mathit{ij}}^{\mathrm{H}}$, $Q_{ij}^{{\text{H}}_2}$) is the electron (H, H₂) collisional (de-)excitation rate. The Einstein coefficient A_ij is defined to be zero for i ≤ j. The excitation rate (i <  j) is related to the de-excitation rate (i >  j) according to

$$\begin{array}{c}\hfill Q_{\mathit{ij}}^{{\mathrm{e},\mathrm{H},\mathrm{H}}_2}=\{\begin{array}{cc}{q}_{\mathit{ij}}^{{\mathrm{e},\mathrm{H},\mathrm{H}}_2}\hfill & i<j(\mathrm{excitation})\hfill \\ \frac{{\omega }_j}{{\omega }_i}exp(\mathrm{\Delta }{E}_{\mathit{ij}}/kT)q_{\mathit{ji}}^{{\mathrm{e},\mathrm{H},\mathrm{H}}_2}\hfill & i>j(\mathrm{de}-\mathrm{excitation})\hfill \end{array},\end{array}$$(50)

where ω_i (ω_j) is the statistical weight of level i (j). Collisions with neutral and molecular hydrogen is a relevant excitation mechanism only in the low-temperature regime when the electron density becomes small. The NC^​2M includes H-excitation of lower lying fine-structure levels of the important coolants CI, CII, OI, SiII, and FeII and H₂-excitation of CI and OI. The set of atomic transitions and the corresponding de-excitation rates are taken from Table 6 in Glover & Jappsen (2007) and from Table B1 in Grassi et al. (2014).

Collisions by electrons is the dominant excitation process at higher temperatures. The electron excitation rate is given by

$$\begin{array}{c}\hfill q_{\mathit{ij}}^{\mathrm{e}}=2.172\times {10}^{-14}\frac{{\mathrm{m}}^3}{\mathrm{s}}{\left(\frac{13.61\mathrm{eV}}{\mathit{kT}}\right)}^{1/2}exp(-\frac{\mathrm{\Delta }{E}_{\mathit{ij}}}{\mathit{kT}})\frac{{\mathrm{{\rm Y}}}_{\mathit{ij}}(T)}{{\omega }_i},\end{array}$$(51)

where Υ_ij(T) is the thermally-averaged collision strength for the i → j transition. To obtain $q_{\mathit{ij}}^{\mathrm{e}}$ the Chianti atomic database version 8 is consulted (Del Zanna et al. 2015). This database comprises critically evaluated atomic data like energy levels, radiation transition probabilities, statistical weights, and collision strength data for all metal ions considered here. In particular, Chianti samples Υ_ij in temperature applying the transition-type-dependent scaling of Burgess & Tully (1992). This scaling allows us to accurately reproduce Υ_ij(T) over a large temperature range and guarantees a well-defined behavior in the limits T → 0, ∞. In order to enable data access, an interface between the NIRVANA code and Chianti database has been designed. For each transition i → j a spline function is generated to represent Υ_ij(T) for all T. The spline-defining parameters are precomputed and stored for each transition, which permits a computationally cheap on-the-fly evaluation of Υ_ij(T).

The Chianti database is very comprehensive and contains many atomic fine-structure levels with a large number of radiative transitions. The volume of Chianti is utilized to study first the impact of the choice of uppermost atomic level L, which is used in the calculation of the ion population fractions, on the cooling rate. Computations are performed separately for each metal ion. Often, a fiducial value of L = 5 is adopted, which might be insufficient to represent ion cooling rates accurately enough. The L-influence on the ion cooling rates, Λ, is illustrated in Fig. 1 as an example for the two ions OIII and FeVI. Figure 1 also shows this influence on the corresponding per-ion cooling coefficients, defined as Λ/n_ionn_e, which in the case of electron-collision-dominated excitation are functions of temperature only and, thus, are particularly suited for tabulation (e.g., Oppenheimer & Schaye 2013a). The temperature dependence⁶ is displayed for L = 5, for the maximum value, which is possible with Chianti (L = 46 for OIII, L = 96 for FeVI) and for the default value used in the NC^​2M (L = 11 for OIII, L = 30 for FeVI). Using five atomic levels, the cooling rate of both ions is off by an order of magnitude compared to the case exploiting the full Chianti database (dashed lines in Fig. 1). Deviations are much less (at most 40%) for the NC^​2M default values, which were adopted as a trade-off between accuracy and computational costs. In summary, for many ions L >  5 is required to ensure sufficient convergence of their cooling rates at all T. The recommended default values of L are cataloged in Table 1. This set represents about 2400 emission lines with photon energies ranging from 1.24 × 10⁻⁴ eV (IR) for the transition 2s²²p^3 2​P_3/2 → 2s²²p^3 2​P_1/2 of NI to 8.26 keV (X-ray) for the transition 3p ^2​P_3/2 → 1s ^2​S_1/2 of FeXXVI.

Metal line cooling is a cost-intensive and complex part of the NC^​2M. To further check its implementation certain spectral line intensity ratios for NII, OIII, and NeIII have been computed as a function of temperature in the vicinity of T = 10⁴ K for which theoretical estimates in the low-density limit are given in Osterbrock & Ferland (2005):

$$\begin{array}{cc}& [\mathrm{NII}]2s^{2}2p^2:\frac{I{(}^1{D}_2{\to }^3{P}_1)+I{(}^1{D}_2{\to }^2{P}_1)}{I{(}^1{S}_0{\to }^1{D}_2)}\simeq \frac{8.23exp(2.5\times {10}^4\mathrm{K}/T)}{1+4.4\times {10}^{-9}\frac{n_{\mathrm{e}}}{{\mathrm{m}}^{-3}}{\left(\frac{T}{\mathrm{K}}\right)}^{-1/2}}\hfill \end{array}$$(52)

$$\begin{array}{cc}& [\mathrm{OIII}]2s^{2}2p^2:\frac{I{(}^1{D}_2{\to }^3{P}_1)+I{(}^1{D}_2{\to }^2{P}_1)}{I{(}^1{S}_0{\to }^1{D}_2)}\simeq \frac{7.9exp(3.29\times {10}^4\mathrm{K}/T)}{1+4.5\times {10}^{-10}\frac{n_{\mathrm{e}}}{{\mathrm{m}}^{-3}}{\left(\frac{T}{\mathrm{K}}\right)}^{-1/2}}\hfill \end{array}$$(53)

$$\begin{array}{cc}& [\mathrm{NeIII}]2s^{2}2p^4:\frac{I{(}^1{D}_2{\to }^3{P}_1)+I{(}^1{D}_2{\to }^2{P}_1)}{I{(}^1{S}_0{\to }^1{D}_2)}\simeq \frac{13.7exp(4.3\times {10}^4\mathrm{K}/T)}{1+3.8\times {10}^{-11}\frac{n_{\mathrm{e}}}{{\mathrm{m}}^{-3}}{\left(\frac{T}{\mathrm{K}}\right)}^{-1/2}}.\hfill \end{array}$$(54)

The line intensity I(i → j)∝p_iΔE_ijA_ij. The outcome is depicted in Fig. 2. Despite the overall good agreement, the NC^​2M calculation is expected to be more precise since it uses temperature-dependent collision strengths whereas the theoretical expressions assume constant values.

3.3. Molecular cooling

Molecular line radiation becomes an important gas coolant for temperatures below a few thousand Kelvin besides atomic fine structure cooling of C, O, Si, and Fe. Among the most relevant molecular coolants in the ISM are the molecules H₂, CO, H₂O, and OH. In primordial or low-metallicity gas the hydrogen deuteride molecule, HD, also contributes to gas cooling if sufficient amounts can be produced by chemical fractionation. The NC^​2M includes these five molecular coolants.

3.3.1. H₂

For densities that are not too high (n ≲ 10⁹ m⁻³) and for temperatures T ≈ 1000 K − 4000 K, molecular hydrogen cooling dominates other molecular cooling agents in the ISM. For higher densities and lower temperatures, CO cooling becomes equally important. To model H₂ cooling I closely follow the work of Glover & Abel (2008) who considered collisional excitation of rotovibrational levels for collision partners H, H⁺, He, H₂, and e⁻. The ratio of ortho-H₂ to para-H₂ is fixed to 3:1 independent of temperature. The cooling rate in an optically thin environment can be written as

$$\begin{array}{c}\hfill {\mathrm{\Lambda }}_{{\mathrm{H}}_2}=\frac{{\mathrm{\Lambda }}_{{\mathrm{H}}_2,\mathrm{LTE}}}{1+{\mathrm{\Lambda }}_{{\mathrm{H}}_2,\mathrm{LTE}}/{\mathrm{\Lambda }}_{{\mathrm{H}}_2,0}}{n}_{{\mathrm{H}}_2},\end{array}$$(55)

where Λ_H2, 0 denotes the sum of cooling rates per H₂ molecule in the low-density limit due to the different collision partners:

$$\begin{array}{c}\hfill {\mathrm{\Lambda }}_{{\mathrm{H}}_2,0}={\displaystyle \sum _s{\mathrm{\Lambda }}_{{\mathrm{H}}_2,s}{n}_s+{\mathrm{\Lambda }}_{{\mathrm{H}}_2,{\mathrm{H}}^{+},J=0\leftrightarrow 1}{n}_{{\mathrm{H}}^{+}}s=\mathrm{H},{\mathrm{H}}^{+},\mathrm{He},{\mathrm{H}}_2,{\mathrm{e}}^{-}.}\end{array}$$(56)

The individual contributions are

$$\begin{array}{c}\hfill log\left(\frac{{\mathrm{\Lambda }}_{{\mathrm{H}}_2,s}}{{10}^{-7}{\mathrm{Jm}}^3{\mathrm{s}}^{-1}}\right)={\displaystyle \sum _{l=0}^5{a}_{s,l}{(log{T}_3)}^l,}\end{array}$$(57)

with the fitting coefficients a_s, l given in Table 8 of Glover & Abel (2008) and T₃ = T/10³ K. Eq. (56) contains the special rate

$$\begin{array}{c}\hfill {\mathrm{\Lambda }}_{{\mathrm{H}}_2,{\mathrm{H}}^{+},J=0\leftrightarrow 1}=4.76\times {10}^{-37}\frac{{\mathrm{Jm}}^3}{\mathrm{s}}[2.25exp(-\frac{170.5\mathrm{K}}{T})-0.75],\end{array}$$(58)

which separately accounts for the cooling due to proton collisional excitation of the rotational transition J = 0 ↔ 1 associated with ortho-to-para state conversion, and this contribution is not implicit in the underlying fitting formula for Λ_H2, H⁺. Equation (58) leads to cooling or heating if the ortho-to-para ratio deviates from its thermal equilibrium value, which is the case here assuming a fixed 3:1 mix.

The function Λ_H2, LTE represents the cooling rate per H₂ molecule when rotovibrational level populations are in local thermodynamic equilibrium (LTE) as expected at high gas densities. The temperature dependence of Λ_H2, LTE is approximated by the formula from Hollenbach & McKee (1979),

$$\begin{array}{cc}\hfill \frac{{\mathrm{\Lambda }}_{{\mathrm{H}}_2,\mathrm{LTE}}}{{\mathrm{J}\mathrm{s}}^{-1}}& =9.5\times {10}^{-29}\frac{T_3^{3.76}}{1+0.12T_3^{2.1}}exp(-\frac{2.197\times {10}^{-3}}{T_3^3})\hfill \\ \hfill & +3\times {10}^{-31}exp(-\frac{0.51}{T_3})+6.7\times {10}^{-26}exp(-\frac{5.86}{T_3})\hfill \\ \hfill & +1.6\times {10}^{-25}exp(-\frac{11.7}{T_3}).\hfill \end{array}$$(59)

In the low-density limit, Λ_H2, LTE ≫ Λ_H2, 0, the effective cooling rate is Λ_H2 → Λ_H2, 0n_H2 whereas in the high-density limit Λ_H2 approaches Λ_H2, LTEn_H2. At high density (n ≳ 10¹⁶ m⁻³) the gas starts to become opaque in the molecular line radiation. This is to the lowest order accounted for by multiplying Λ_H2 with the attenuation factor

$$\begin{array}{c}\hfill {\tau }_{{\mathrm{H}}_2}=min[1,{\left(\frac{n}{8\times {10}^{15}{\mathrm{m}}^{-3}}\right)}^{-0.45}]\end{array}$$(60)

proposed by Ripamonti & Abel (2004), and which depends only on the local density.

At even higher densities (n ≳ 10²⁰ m⁻³) collision-induced emission (CIE) becomes a further important cooling effect. In the CIE process, deep collisions between H₂ molecules or between H₂ and H or He form a short-term supramolecular complex with an induced non-zero electric dipole. Such a configuration has a high probability of an energy transition followed by photon emission or absorption. According to Ripamonti & Abel (2004), the CIE cooling rate is written as

$$\begin{array}{c}\hfill {\mathrm{\Lambda }}_{\mathrm{CIE},\mathrm{thick}}={\mathrm{\Lambda }}_{\mathrm{CIE}}min[1,\frac{1-e^{-{\tau }_{\mathrm{CIE}}}}{{\tau }_{\mathrm{CIE}}}],\end{array}$$(61)

where Λ_CIE is the cooling rate in the optically thin limit. The optical depth correction parameter τ_CIE is given by

$$\begin{array}{c}\hfill {\tau }_{\mathrm{CIE}}={\left(\frac{n_{{\mathrm{H}}_2}}{7\times {10}^{21}{\mathrm{m}}^{-3}}\right)}^{2.8},\end{array}$$(62)

to approximately account for cooling in the optically thick regime when n ≳ 10²² m⁻³ (Hirano & Yoshida 2013). Under the assumption of LTE, Λ_CIE is expressed by

$$\begin{array}{c}\hfill {\mathrm{\Lambda }}_{\mathrm{CIE}}={\displaystyle \sum _{s}4\pi {\int }_0^{\infty }{\alpha }_{s,\nu }(T,n_{{\mathrm{H}}_2},n_s)B_{\nu }(T)\mathrm{d}\nu ,s={\mathrm{H}}_2,\mathrm{He},\mathrm{H},}\end{array}$$(63)

summing up cooling contributions for the collision pairs H₂–H₂, H₂–He, and H₂–H. The quantity B_ν(T) in units Js⁻¹ m⁻² Hz⁻¹ denotes the black-body spectral radiance density (Planck’s function) and α_s, ν, for each collision pair s, is the spectral collision-induced absorption coefficient that depends on T, n_H2, and the number density n_s of the respective collision partner. In order to perform the integral over frequency, ν data for α_s, ν from various sources is used as listed in Table 2. These sources sample absorption coefficients per H₂ molar density and per molar density of collision partner, ${\stackrel{\sim }{\alpha }}_{s,\nu }$ (in units⁷ cm⁻¹ amg⁻²), in the frequency interval and temperature range as specified in Table 2 for the different collision pairs. The specific ${\stackrel{\sim }{\alpha }}_{s,\nu }$ does not depend on the number densities of collision partners and is related to α_s, ν (in units m⁻¹) in Eq. (63) according to

$$\begin{array}{c}\hfill {\alpha }_{s,\nu }(T,n_{{\mathrm{H}}_2},n_s)=1.3854\times {10}^{-49}{n}_{{\mathrm{H}}_2}{n}_s{\stackrel{\sim }{\alpha }}_{s,\nu }(T).\end{array}$$(64)

Equation (63) can then be rewritten as

$$\begin{array}{c}\hfill {\mathrm{\Lambda }}_{\mathrm{CIE}}={\displaystyle \sum _s{\mathrm{\Lambda }}_{\mathrm{CIE},s}{n}_s{n}_{{\mathrm{H}}_2},s={\mathrm{H}}_2,\mathrm{He},\mathrm{H},}\end{array}$$(65)

with frequency-integrated cooling rate coefficients ${\mathrm{\Lambda }}_{\mathrm{CIE},s}(T)\propto \int \mathrm{d}\nu {\stackrel{\sim }{\alpha }}_{s,\nu }(T)B_{\nu }(T)$, which are only functions of temperature. The integration in ν is performed with the trapezoidal rule. Polynomial interpolation in logT is used to accurately represent those integrated rates as a function of temperature within the ranges given in Table 2, and a power-law extrapolation is adopted at higher temperatures for which no data exists:

$$\begin{array}{cc}\hfill log\left(\frac{{\mathrm{\Lambda }}_{\mathrm{CIE},s}}{{\mathrm{Js}}^{-1}{\mathrm{m}}^3}\right)& =\{\begin{array}{cc}{\sum }_l{a}_l{[log(T/\mathrm{K})]}^l\hfill & T<T_{\mathrm{u}}\hfill \\ b+clog(T/T_{\mathrm{u}})\hfill & T\ge T_{\mathrm{u}}\hfill \end{array}.\hfill \end{array}$$(66)

The fitting coefficients a_l, b, c, T_u are summarized in Table 3 for the various collision pairs. The temperature dependence is shown in Fig. 3 demonstrating the quality of Λ_CIE, s-fitting in the regimes where data is available. The approximation becomes uncertain for very low temperatures but the fitting formula can safely be applied.

3.3.2. HD

In primordial gas or low-metallicity gas, the HD molecule can be an important coolant at temperatures below a few hundred Kelvin and densities n ≳ 10¹⁰ m⁻³: a regime where H₂ cooling becomes increasingly ineffective. HD cooling is represented in the NC^​2M by the cooling function of Lipovka et al. (2005). This cooling function is based on a detailed level population balance calculation covering rotovibrational transitions with ν ≤ 3 and J ≤ 8, and applying excitation rates for collisions with H.

Accordingly, the HD cooling rate is given by

$$\begin{array}{c}\hfill {\mathrm{\Lambda }}_{\mathrm{HD}}={\stackrel{\sim }{\mathrm{\Lambda }}}_{\mathrm{HD}}{n}_{\mathrm{HD}},\end{array}$$(67)

with the fitting function

$$\begin{array}{c}\hfill log\left(\frac{{\stackrel{\sim }{\mathrm{\Lambda }}}_{\mathrm{HD}}}{{10}^{-7}{\mathrm{Js}}^{-1}}\right)={\displaystyle \sum _{l,m=0}^4{D}_{l,m}{[log\left(\frac{n_{\mathrm{eff}}}{{10}^6{\mathrm{m}}^{-3}}\right)]}^l{[log\left(\frac{T}{\mathrm{K}}\right)]}^m.}\end{array}$$(68)

The coefficients D_l, m are tabulated in Lipovka et al. (2005), Table 1. Since the excitation rate coefficients for HD-He collisions and HD–H₂ collisions are comparable to the excitation rate for HD–H collisions, as pointed out by Flower & Le Bourlot (2000), an effective number density n_eff = n_H + n_H2 + n_He is used in Eq. (68). Lipovka et al. (2005) emphasize the high accuracy of the above approximation for 100 K ≲ T ≲ 2 × 10⁴ K and 10⁶ m⁻³ ≲ n_eff ≲ 10¹⁴ m⁻³. Within this density range the fitting function behaves well with respect to T and can be safely applied outside of the prescribed temperature range. Outside the given density range, however, the fitting function is ill-conditioned. To circumvent mathematical problems, the HD cooling rate is therefore extended according to

$$\begin{array}{c}\hfill {\mathrm{\Lambda }}_{\mathrm{HD}}⟶\{\begin{array}{cc}\frac{n_{\mathrm{eff}}}{{10}^6{\mathrm{m}}^{-3}}{\mathrm{\Lambda }}_{\mathrm{HD}}\hfill & n_{\mathrm{eff}}<{10}^6{\mathrm{m}}^{-3}\hfill \\ {\mathrm{\Lambda }}_{\mathrm{HD}}(n_{\mathrm{eff}}={10}^{14}{\mathrm{m}}^{-3})\hfill & n_{\mathrm{eff}}>{10}^{14}{\mathrm{m}}^{-3},\hfill \end{array}\end{array}$$(69)

that is, it is linearly extrapolated down to zero density from the lower density bound, and it is assumed to be density-independent above the upper bound.

3.3.3. CO and H₂O

The NC^​2M implementation of CO cooling and H₂O cooling is based on the work of Neufeld & Kaufman (1993) and Neufeld et al. (1995) with some modifications as addressed in GFM10. The CO and H₂O cooling rates are given by

$$\begin{array}{c}\hfill {\mathrm{\Lambda }}_s=({\mathrm{\Lambda }}_{s,\mathrm{rot}}{n}_{\mathrm{eff},\mathrm{rot},s}+{\mathrm{\Lambda }}_{s,\mathrm{vib}}{n}_{\mathrm{eff},\mathrm{vib},s})n_s,s={\mathrm{CO},\mathrm{H}}_2\mathrm{O},\end{array}$$(70)

where Λ_s, rot and Λ_s, vib are cooling rate coefficients due to rotational and vibrational transitions, respectively, and n_{eff, rot, s}, n_{eff, vib, s} are effective number densities for collisions with H₂, H, and e⁻. In the Neufeld et al. papers, n_{eff, rot, s} = n_{eff, vib, s} = n_H2 assuming H₂ dominates collisional excitation of CO and H₂O. The cooling rate coefficients are determined from

$$\begin{array}{c}\hfill \frac{1}{{\mathrm{\Lambda }}_{s,\mathrm{rot}}}=\frac{1}{{\mathrm{\Lambda }}_{s,\mathrm{rot},0}}[1+{\left(\frac{n_{\mathrm{eff},\mathrm{rot},s}}{n_{1/2,s}}\right)}^{{\alpha }_s}(1-\frac{n_{1/2,s}{\mathrm{\Lambda }}_{s,\mathrm{rot},0}}{{\mathrm{\Lambda }}_{s,\mathrm{rot},\mathrm{LTE}}})]+\frac{n_{\mathrm{eff},\mathrm{rot},s}}{{\mathrm{\Lambda }}_{s,\mathrm{rot},\mathrm{LTE}}},\end{array}$$(71)

for rotational transitions and from

$$\begin{array}{c}\hfill \frac{1}{{\mathrm{\Lambda }}_{s,\mathrm{vib}}}=\frac{1}{{\mathrm{\Lambda }}_{s,\mathrm{vib},0}}+\frac{n_{\mathrm{eff},\mathrm{vib},s}}{{\mathrm{\Lambda }}_{s,\mathrm{vib},\mathrm{LTE}}},\end{array}$$(72)

for vibrational transitions. Here, Λ_{s, rot, 0} and Λ_{s, vib, 0} denote the cooling coefficients in the low-density limit, Λ_{s, rot, LTE} and Λ_{s, vib, LTE} represent cooling rates per s molecule in case the level population is in LTE, and n_1/2, s is the effective density at which Λ_s, rot = Λ_{s, rot, 0}/2. The low-density rate coefficients are pure functions of the temperature whereas Λ_{s, rot, LTE}, Λ_{s, vib, LTE}, n_1/2, s, and exponent α_s depend, in addition, on an optical depth variable ${\stackrel{\sim }{N}}_s$. This variable is defined by

$$\begin{array}{c}\hfill {\stackrel{\sim }{N}}_s=\frac{n_s}{|\mathrm{\nabla }·\mathbf{v}|},\end{array}$$(73)

meaning an effective column density per unit velocity. The functions Λ_{s, rot, LTE}, Λ_{s, vib, LTE}, n_1/2, s, and α_s are sampled in Neufeld & Kaufman (1993), Tables 2, 3, and 5, and Neufeld et al. (1995), Tables 1–3, in $\{T,log{\stackrel{\sim }{N}}_s\}$-space for the range

$$\begin{array}{c}\hfill \left(\frac{T}{\mathrm{K}}\right)\times log\left(\frac{{\stackrel{\sim }{N}}_s}{\frac{{\mathrm{cm}}^{-2}}{\mathrm{km}{\mathrm{s}}^{-1}}}\right)=\{\begin{array}{cc}[10,4000]\times [10,19]\hfill & {\mathrm{H}}_2\mathrm{O},\mathrm{rot}\hfill \\ [10,100]\times [14.5,19]\hfill & \mathrm{CO},\mathrm{rot}\hfill \\ [100,2000]\times [14,19]\hfill & \mathrm{CO},\mathrm{rot}\hfill \\ [100,4000]\times [13,20]\hfill & {\mathrm{H}}_2\mathrm{O}/\mathrm{CO},\mathrm{vib}.\hfill \end{array}\end{array}$$

Linear interpolation within these tables is applied for function evaluation, and constant function continuation is adopted outside of the validated range.

The expressions for the effective number densities are taken from GFM10 referring to the work of Yan (1997), Faure et al. (2004), and Meijerink & Spaans (2005). Accordingly, for rotational cooling

$$\begin{array}{c}\hfill n_{\mathrm{eff},\mathrm{rot},s}=\{\begin{array}{c}{n}_{{\mathrm{H}}_2}+10n_{\mathrm{H}}+1291n_{\mathrm{e}}{\left(\frac{T}{\mathrm{K}}\right)}^{-1/2}\hfill \\ \times {10}^{15.749{\left(\frac{T}{\mathrm{K}}\right)}^{-1/6}-47.137{\left(\frac{T}{\mathrm{K}}\right)}^{-1/3}+76.648{\left(\frac{T}{\mathrm{K}}\right)}^{-1/2}-60.191{\left(\frac{T}{\mathrm{K}}\right)}^{-2/3}}\hfill & \\ {\mathrm{H}}_2\mathrm{O}\hfill & \\ n_{{\mathrm{H}}_2}+1.75n_{\mathrm{H}}{\left(\frac{T}{\mathrm{K}}\right)}^{1/4}+680n_{\mathrm{e}}{\left(\frac{T}{\mathrm{K}}\right)}^{-1/4}\mathrm{CO},\hfill \end{array}\end{array}$$(74)

whereas for vibrational cooling

$$\begin{array}{c}\hfill n_{\mathrm{eff},\mathrm{vib},s}=\{\begin{array}{cc}{n}_{{\mathrm{H}}_2}+10n_{\mathrm{H}}+4.06\times {10}^8{n}_{\mathrm{e}}\frac{exp\left(\frac{47.5}{{(T/\mathrm{K})}^{1/3}}\right)}{\sqrt{T/\mathrm{K}}}\hfill & {\mathrm{H}}_2\mathrm{O}\hfill \\ n_{{\mathrm{H}}_2}+50n_{\mathrm{H}}+43n_{\mathrm{e}}{\left(\frac{T}{\mathrm{K}}\right)}^{0.938}exp\left(\frac{68}{{(T/\mathrm{K})}^{1/3}}\right)\hfill & \mathrm{CO},\hfill \end{array}\end{array}$$(75)

3.3.4. OH

Cooling of the OH molecule in the NC^​2M is based on the model developed in Hollenbach & McKee (1979) for rotational cooling by molecules with dipole moment including optical depth effects and collisional de-excitation.

The OH cooling rate is written as

$$\begin{array}{c}\hfill {\mathrm{\Lambda }}_{\mathrm{OH}}={\stackrel{\sim }{\mathrm{\Lambda }}}_{\mathrm{OH}}{n}_{\mathrm{OH}}.\end{array}$$(76)

Following Hollenbach & McKee (1979), the cooling rate per OH molecule, ${\stackrel{\sim }{\mathrm{\Lambda }}}_{\mathrm{OH}}$, is given by

$$\begin{array}{c}\hfill {\stackrel{\sim }{\mathrm{\Lambda }}}_{\mathrm{OH}}=2.35\times {10}^{-26}\frac{\mathrm{J}}{\mathrm{s}}{\left(\frac{T}{\mathrm{K}}\right)}^2\left[\frac{2+y_{\mathrm{m}}+0.6y_{\mathrm{m}}^2}{1+c_{\tau }+n_{\mathrm{cr}}/n_{{\mathrm{H}}_{\mathrm{nuc}}}+1.5{(n_{\mathrm{cr}}/n_{{\mathrm{H}}_{\mathrm{nuc}}})}^{0.5}}\right]\end{array}$$(77)

where

$$\begin{array}{c}\hfill y_{\mathrm{m}}=ln[1+\frac{c_{\tau }}{1+10n_{\mathrm{cr}}/n_{{\mathrm{H}}_{\mathrm{nuc}}}}],\end{array}$$(78)

$$\begin{array}{c}\hfill n_{\mathrm{cr}}=1.5\times {10}^{16}{\mathrm{m}}^{-3}\sqrt{\frac{T}{{10}^3\mathrm{K}}},\end{array}$$(79)

and

$$\begin{array}{c}\hfill c_{\tau }=\frac{{\tau }_T{[2\pi ln(2.13+{({\tau }_T/e)}^2)]}^{0.5}}{e^{-{\tau }_{\mathrm{d}}}/(1+{\tau }_{\mathrm{d}}^2)+2{\tau }_{\mathrm{d}}{[ln(1+{\tau }_T/e)ln({\tau }_T/{\tau }_{\mathrm{d}}e)]}^{0.5}}·\end{array}$$(80)

The quantity n_cr denotes a critical density above which collisional de-excitation starts to become important and c_τ represents the optically thick regime. In the expression for c_τ, the quantity e is the Euler constant, and τ_T and τ_d are the optical depths in the lines and due to dust absorption, respectively. In the one-cell numerical experiments below, OH cooling is considered in the optically thin limit, which is obtained by setting τ_T = τ_d = 0.

3.4. Bremsstrahlung

The NC^​2M takes into account plasma radiation losses in free-free transitions (Bremsstrahlung). A formula for the frequency-integrated rate is given in Shapiro & Kang (1987),

$$\begin{array}{c}\hfill {\mathrm{\Lambda }}_{\mathrm{ff}}=1.426\times {10}^{-40}\frac{{\mathrm{J}\mathrm{m}}^3}{\mathrm{s}}\sqrt{\frac{T}{\mathrm{K}}}{n}_{\mathrm{e}}{\displaystyle \sum _i{z}_i^2{g}_{\mathrm{ff}}(T,z_i)n_i,}\end{array}$$(81)

where the sum is over relevant ions with charge number z_i. The mean Gaunt factor, g_ff, is given by

$$\begin{array}{c}\hfill g_{\mathrm{ff}}=\{\begin{array}{cc}0.79464+0.1243log\left(\frac{T}{\mathrm{K}}{z}_i^{-2}\right)\hfill & T{z}_i^{-2}<3.2\times {10}^5\mathrm{K}\hfill \\ 2.13164-0.1240log\left(\frac{T}{\mathrm{K}}{z}_i^{-2}\right)\hfill & T{z}_i^{-2}>3.2\times {10}^5\mathrm{K}\hfill \end{array}.\end{array}$$(82)

Bremsstrahlung losses become significant at T ≈ 10⁶ K where they compare with H and He atomic line cooling, and they dominate all metal cooling processes for T ≳ 10⁷ K. The main contribution comes from the ions H⁺, He⁺, and He⁺⁺. In a fully ionized, solar abundance plasma, the high-z metal ions of O, Ne, and Fe contribute to about another 10%. Hence, those ions were included in the sum of Eq. (81).

3.5. Cosmic ray heating

The kinetic energy of the ejected electron in cosmic ray ionizations is partly converted into gas heating by Coloumb scattering, momentum transfer, and collisional de-excitation of excited H₂ rotational states. The amount of energy deposited in the gas as heat is often assumed to be a constant fraction of the ejected electron energy. For instance, the models in GFM10 and RSO14 use a canonical value of 20 eV per primary ionization. More generally, the heating fraction must be considered taking into account the ionization degree of the gas and its H-to-H₂ abundance ratio. According to Dalgarno et al. (1999) the heating efficiency in a partly ionized H–H₂–He gas mixture can be represented by

$$\begin{array}{c}\hfill \eta =\frac{10n_{{\mathrm{H}}_2}{\eta }_{{\mathrm{H}}_2-\mathrm{He}}+n_{\mathrm{H}}{\eta }_{\mathrm{H}-\mathrm{He}}}{10n_{{\mathrm{H}}_2}+n_{\mathrm{H}}},\end{array}$$(83)

where η_H − He and η_{H2 − He} are the efficiencies for separate H–He and H₂-He gas mixtures each with a 10% helium fraction in particle numbers. The individual efficiencies are given by the fitting formulae

$$\begin{array}{cc}\hfill {\eta }_{\mathrm{H}--\mathrm{He}}& =1-0.74{[1+164.8{(n_{\mathrm{e}}/n_{\mathrm{H}})}^{0.984}]}^{-1}\hfill \end{array}$$(84)

$$\begin{array}{cc}\hfill {\eta }_{{\mathrm{H}}_2-\mathrm{He}}& =1-0.9{[1+16.11{(n_{\mathrm{e}}/n_{{\mathrm{H}}_2})}^{0.515}]}^{-1},\hfill \end{array}$$(85)

and they are accurate to within 10% to the Monte-Carlo simulation results of Dalgarno et al. (1999) up to moderate gas ionization levels, n_e/n_H ≲ 0.1 or n_e/n_H2 ≲ 0.1, respectively.

In the NC^​2M the total heating rate for a mean ejected electron energy of 30 eV due to H and H₂ primary cosmic ray ionizations is then expressed by

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathit{cr}}=(30\mathrm{eV})\eta ({\xi }_{\mathrm{H}}{n}_{\mathrm{H}}+{\xi }_{{\mathrm{H}}_2}{n}_{{\mathrm{H}}_2}).\end{array}$$(86)

To give some numbers, for a neutral H–He (H₂-He) gas mixture one obtains a heating efficiency of η = 0.26 (η = 0.1) corresponding to a mean energy per cosmic ray ionization of 7.8 eV (3 eV). The presence of free electrons substantially increases the energy loss of the ejected electron by Coloumb scattering. For a partly ionized H–He (H₂-He) gas mixture with n_e = 0.1n_H (n_e = 0.1n_H2), the heating efficiency grows to η = 0.96 (η = 0.85) corresponding to a mean energy per cosmic ray ionization of 28.8 eV (25.5 eV).

4. Equilibrium chemistry and cooling

The NC^​2M has been developed with the ultimate aim of following non-equilibrium chemistry in multidimensional magnetohydrodynamics applications. However, it is essential to first check the implementation by considering chemical equilibrium models. To find the chemical equilibrium for the network specified in Tables A.1–A.3 at a fixed temperature, density, and elemental composition the numerical solver as described in Ziegler (2016) is used. The system of kinetic rate equations is evolved for a sufficiently long time until the solution has converged close to chemical equilibrium, which is defined by a relative residual of 10⁻⁶. The equilibrium species abundances are calculated as a function of temperature in the range 10 K ≤ T ≤ 10⁸ K at 801 nodes with ΔlogT = 8.75 × 10⁻³ for three different densities n_H = 10⁴ m⁻³, 10⁶ m⁻³, 10¹⁰ m⁻³ in order to cover ISM conditions from dilute to dense gas. At this point we recall that in the presence of an UV background, dilute gas is not shielded from photoionizing or photodissociating radiation, which certainly would have a profound influence on the low-temperature equilibrium balance obtained here. The calculations assume an elemental composition given by default solar abundances as cataloged in Grevesse & Sauval 1998. This results in a total metallicity of Z = 0.96 Z_⊙. Other metals missing in the NC^​2M make a contribution of 0.04 Z_⊙.

In order to illustrate the potential importance of grain-surface catalysis and cosmic ray ionization effects, equilibrium computations are performed for the full chemical network including dust and cosmic rays as well as for reduced networks neglecting either the dust component or the cosmic ray component. A pure gas-phase network is also considered ignoring both dust and cosmic rays. Variations in metallicity are not investigated. In the models including cosmic ray reactions, a canonical cosmic ray ionization rate of ξ_H = 2.5 × 10⁻¹⁷ s⁻¹ is adopted. For grain-surface reactions, the dust temperature, T_d, and the parameter G appearing in expression (23) for grain charging need to be defined. As in RSO14 we adopt a value of T_d = 10 K since the rate coefficient given by Eq. (17) is not very sensitive to the dust temperature in the regime 6 K ≲ T_d ≲ 50 K. As far as concerns grain charging, we set G ≡ 0 consistent with the limiting case of no radiation, that is, no photoelectric charging effect.

4.1. Molecule abundances

Equilibrium abundances for the three important molecular coolants H₂, CO, and H₂O and for the further molecules O₂ and OH are shown in Fig. 4 as a function of temperature for different gas densities (different panels in a row) and for the considered chemical network variants (different line styles of curves). Molecular fractions have been normalized to the H nuclei density in the case of H₂, C nuclei density in the case of CO, and O nuclei density in the case of H₂O, O₂, and OH. Only the low-temperature regime is shown since molecules are effectively absent at T ≳ 2 × 10⁴ K irrespective of density.

Molecular abundances for the full network case are depicted as solid curves in Fig. 4. It is first noted that the obtained H₂ fraction is larger than in RSO14 because dust-catalytic H₂ formation is more effective due to the higher dust total cross section used in our standard model. However, when assuming identical dust parameters, excellent agreement can be found between RSO14 and the full network case. The resulting CO distribution also matches well with that in RSO14 under this assumption⁸. At moderate densities the curves of the O₂ and OH fractions closely follow that of H₂O. The H₂, CO, and H₂O fractions generally grow with increasing density.

At the highest density, n_H = 10¹⁰ m⁻³, most hydrogen is in molecular form and the H₂ fraction approaches unity up to a temperature of ∼2500 K, and rapidly drops hereafter. The weakly visible cut-off near this temperature is attributed to the functional form of coefficient Eq. (17) for dust-surface catalysis of H₂. This cut-off is also striking as sharp transitions in the other molecules, since H₂ largely controls their chemistry. The CO fraction is close to unity for T ≲ 500 K and in the window 2500 K ≲ T ≲ 4000 K. In between, 500 K ≲ T ≲ 2500 K, the CO fraction decreases to ∼10⁻² associated with the increase in H₂O fraction. Evidently, more oxygen becomes bound in H₂O, which is no longer available to form CO. Also, the formation of molecular oxygen is largely suppressed in the temperature regime 300 K ≲ T ≲ 2500 K in favor of the synthesis of water.

The absence of cosmic ray reactions has a paramount influence on the chemical equilibrium abundances since cosmic rays as the dominant driver for ionization at lower temperature ceases, and the ionization balance is determined by collisional effects and neutralization of ions on dust grains only. This results in a very low level of ionization and the chemistry is controlled by neutral-neutral reactions and H₂ formation on grains. Abundances for the CR-free network are depicted as dash-dotted curves in Fig. 4. Up to a temperature of ∼4000 K (∼3300 K for n_H = 10¹⁰ m⁻³) hydrogen is predominantly bound in H₂ followed by a smaller fraction bound in H₂O. The increased amount of H₂ compared to the full network case is the result of missing cosmic ray ionization and CR-induced photodissociation of H₂. The fraction of atomic hydrogen is very low in this temperature regime: n_H/n_Hnuc <  10⁻⁶ (< 10⁻⁸, < 10⁻¹²) at n_H = 10⁴ m⁻³ (n_H = 10⁶ m⁻³, n_H = 10¹⁰ m⁻³) compared to n_H/n_Hnuc >  0.8 (> 0.18, > 5 × 10⁻⁵) in the full network case. Likewise, carbon is almost completely bound in CO with a fraction close to unity below T ≈ 5000 K (≈4500 K for n_H = 10¹⁰ m⁻³). Other C-bearing molecules like CH₂, for example, exist only in small amounts with a fraction of n_CH2/n_Cnuc <  10⁻¹⁰ for this molecule. There are three distinct temperature ranges where the H₂O fraction is nearly constant. At a density n_H = 10⁴ m⁻³ the intervals are T ≲ 220 K, 250 K ≲ T ≲ 1000 K, and 1100 K ≲ T ≲ 4000 K, respectively. The intervals slightly change with density, and for the highest density the first two temperature regimes seem to merge. The O₂ fraction drops by many order of magnitudes in the third temperature regime where the H₂O fraction reaches its largest value. This O₂ depletion effect is also seen in the high-density full network case but in the CR-free network it occurs at all densities.

Equilibrium abundances for the dust-free network are depicted as dashed curves in Fig. 4. The absence of dust-assisted reactions, in particular the catalysis of H₂, has a substantial impact on the H₂ abundance, and therefore on the overall molecular chemistry. Molecular fractions are substantially reduced at lower temperatures compared to the full network case except at the highest density where three-body reactions contribute to H₂ formation. Therefore, despite the lack of dust, the H₂ fraction at n_H = 10¹⁰ m⁻³ is roughly one order of magnitude larger than for the other densities for T ≲ 2500 K. Compared with the full network case, the buckle in CO around 1000 K is absent, and a substantial amount of H₂O (O₂) is present in the vicinity of 1000 K (400 K).

Abundances for the pure gas-phase network, which ignores both dust and cosmic rays, are depicted as dotted curves in Fig. 4. At moderate densities and below a temperature of ∼1000 K the H₂ fraction is up to two orders of magnitude smaller than in the full network case. At the same time the H₂O fraction is larger by up to several orders of magnitude, and is nearly constant. The CO fraction is close to one up to T ≈ 5000 K and a somewhat smaller temperature at the highest density. The buckle in the CO distribution most prominent in the high-density full network case vanishes since its appearance is tightly correlated with a simultaneous high H₂O fraction. Below a certain temperature both O₂ and OH get substantially depleted independent of density because most oxygen is bound in H₂O.

4.2. Ionization structure

The ionization structure is illustrated in Figs. 5–8 as examples for the metal elements C, O, Si, and Fe. As expected, the ionization fractions in the high-temperature regime are pure functions of the temperature and are independent of density. Ion recombination on dust and cosmic ray processes has a negligible effect at higher temperatures, and the ionization fraction is solely determined by the balance between electron collisional ionization, electron-ion recombination, and, to a lesser extent, charge transfer process. The high-temperature ionization structure, in particular the shapes and peak positions of the individual ionization state curves, is in excellent agreement with previous literature results, such as Sutherland & Dopita (1993) or Gnat & Sternberg (2007).

At temperatures T ≲ 10⁴ K the lower ionization stages of metal elements show an explicit dependence of their fractions on the gas density, dust content, and cosmic ray strength. In the absence of cosmic rays (dash-dotted lines) the fractions of SiII and FeII, for instance, are substantially reduced compared to the full network case (solid lines) where they are close to unity except for the highest density. The overall impact of dust is less pronounced at lower densities but becomes apparent for n_H = 10¹⁰ m⁻³. For instance, the fraction of FeII grows by three orders of magnitude in the absence of dust (dashed lines) from ∼4 × 10⁻⁴ to a value of ∼0.3. Generally, the CI and OI abundance is heavily influenced by molecular reactions taking place in the low-temperature regime. This is prominent in CI at n_H = 10¹⁰ m⁻³ where a significant amount of carbon is bound in CO. Molecule formation of O₂ and H₂O leads to the broad gap in the abundance of OI around T = 1000 K especially when cosmic rays are absent.

4.3. Cooling

In this section cooling functions computed under the assumption of chemical equilibrium are presented. For convenience, the functions are normalized to the H nuclei number density squared, ${\mathrm{\Lambda }}_{*}=\mathrm{\Lambda }/n_{{\text{H}}_{\text{nuc}}}^2$ ($\mathrm{\Gamma }/n_{{\text{H}}_{\text{nuc}}}^2$ for heating terms, ${\mathrm{\Lambda }}_{\text{net}}/n_{{\text{H}}_{\text{nuc}}}^2$ for net cooling), where Λ_* has units J m³ s⁻¹ (=10¹³ erg cm³ s⁻¹). Calculations are shown for the full network at a density of n_H = 10⁶ m⁻³ if not otherwise stated.

Figure 9 shows equilibrium cooling functions for the molecules H₂, CO, H₂O, and OH resulting from line cooling. The H₂ cooling function, in addition, contains contributions from reaction kinetics (dissociation cooling and formation heating). Also, CIE cooling is included in the H₂ function which is, however, unimportant at the gas densities explored here. The H₂ cooling function is negative for T ≲ 40 K and between T ≈ 60 K and T ≈ 100 K, that is, there is net heating indicated by dash-dotted lines in Fig. 9. The heating below 40 K is caused by chemical heating mainly due to H₂ formation on dust grains. The second heating window is a consequence of fixing the ortho-to-para H₂ ratio to 3:1, which results in heating (or cooling) below (or above) T ≈ 155 K according to Eq. (58). This process interleaves with the other thermal processes and leads to net heating in this temperature interval. If Eq. (58) were skipped, the second heating window would vanish. H₂ cooling dominates molecular cooling except in the vicinity of T ≈ 1000 K where cooling due to H₂O peaks, and above T ≈ 6000 K where CO significantly contributes to molecular cooling. The cooling due to OH is more than one order of magnitude below that of H₂ and CO over the whole temperature range at this density. In general, the relative cooling contributions of these molecules are expected to depend on the gas density, the amount of dust present and cosmic ray effects regulating the molecular chemistry.

The fine-structure line cooling produced by metal species is illustrated in Fig. 10. The top panel example shows the cooling function of oxygen expanded in its contributions from individual ionization stages OI–OVIII. The bottom panel of Fig. 10 gives the contributions from each metal element. In the low-temperature regime the primary metal coolants are C, O, Si, and Fe. This is consistent with the findings of RSO14 for their radiation-free network. However, the importance of these metal elements for cooling at low temperature may be differently qualified when the effect of metal depletion onto dust grains is considered. As shown by Jenkins (2009) metal depletion can substantially diminish the fraction of metals in the gas phase, hence reducing their cooling contribution. The NC^​2M does not yet account for such depletion effects. In the high-temperature regime metal cooling is dominated by C and O between 2 × 10⁴ K ≲ T ≲ 3.5 × 10⁵ K, and by Fe for T ≳ 6 × 10⁵ K. Around T ≈ 5 × 10⁵ K Ne provides an important contribution. The computation of metal line cooling was based on the recommendations summarized in Table 1. When assuming only a five level atomic fine-structure model (L = 5) for each ion, the amount of cooling is underestimated by up to a factor of five, especially in the vicinity of 10⁶ K (dashed line in Fig. 10, bottom panel).

The total cooling function representing the sum of the distinctive cooling and heating agents is shown in Fig. 11. We recall that the following results are obtained for the full network at n_H = 10⁶ m⁻³. There is net heating up to 50 K caused by the dominance of cosmic ray heating. There is a further narrow window of net heating between T ≈ 56 K and T ≈ 69 K, which is the fingerprint of proton-collisional H₂ excitations involving nuclear spin reversals (Eq. (58)). The corresponding cooling function for the radiation-free network in RSO14, Fig. 3, does not show this complicated structure below 100 K but states a unique thermal equilibrium temperature of about 100 K. Ignoring Eq. (58) our model would likewise achieve a unique thermal equilibrium at ∼76 K. The temperature regime from a few hundred Kelvin to ∼8000 K is governed by molecular cooling. The total cooling function reaches a maximum at T ≈ 1.6 × 10⁴ K, which is related to cooling processes of hydrogen and helium. Metal cooling then dominates over the large interval from T ≈ 3 × 10⁴ K to T ≈ 10⁷ K. Cooling due to Bremsstrahlung becomes unimportant before ∼10⁶ K but dominates above ∼10⁷ K.

The impact of different densities, absence of dust, and absence of cosmic rays on the total cooling function is presented in Fig. 12. As can be seen the cooling behavior is insensitive to such influence in the high-temperature regime but reveals a strong dependence in the low-temperature regime. The absence of dust (but presence of cosmic rays) reduces molecule fractions and thus diminishes molecular cooling and therefore lowers total cooling (yellow curve in Fig. 12). In the absence of cosmic rays (but presence of dust) molecule fractions are higher and molecular cooling is increased, in particular below T ≈ 100 K (CO and H₂O) and above T ≈ 1000 K (H₂). Because of this and because the low degree of ionization renders Eq. (58) unimportant, there is no net heating at T ≥ 10 K (light blue curve in Fig. 12). The cooling function for the full network at the lower density n_H = 10⁴ m⁻³ is depicted as a red line in Fig. 12. Since cosmic ray heating is roughly proportional to the gas density it becomes relatively more important at lower densities than cooling processes, being proportional to the density squared in the low-density limit. Consequently, the thermal equilibrium temperature lies much higher at T ≈ 700 K. For the full network with n_H = 10¹⁰ m⁻³ there is no heating at T ≥ 10 K, and the cooling efficiency is reduced by roughly one order of magnitude at T ≲ 1000 K (dark blue curve in Fig. 12) compared with the n_H = 10⁶ m⁻³ case. The reduction comes about because, with increasing gas density, molecular cooling approaches the LTE regime, which lowers the cooling rate per $n_{{\text{H}}_{\text{nuc}}}^2$, that is, Λ_*.

5. Non-equilibrium chemistry and cooling

The cooling functions computed for chemical equilibrium in the last section are now compared with cooling functions in time-dependent situations where chemical equilibrium is not prevailing. To demonstrate the impact of non-equilibrium effects, the chemical rate equations and thermal rate equation are solved simultaneously under isochoric conditions using a modified version of the NIRVANA code operating on a single numerical cell. We focus on the isochoric case here since constant density is usually assumed during the chemistry and cooling update in a hydrodynamic simulation time step. Under isochoric assumption the temperature evolves according to

$$\begin{array}{c}\hfill \frac{\mathrm{d}T}{\mathrm{d}t}=-\frac{(\gamma -1){\mathrm{\Lambda }}_{\mathrm{net}}}{k{n}_{\mathrm{tot}}}-\frac{T}{n_{\mathrm{tot}}}\frac{\mathrm{d}{n}_{\mathrm{tot}}}{\mathrm{d}t},\end{array}$$(87)

where the time derivative of total number density n_tot can be replaced by making use of the reaction source terms as shown in Ziegler (2016). The adiabatic index is approximated by

$$\begin{array}{c}\hfill \gamma =\frac{5(n_{\mathrm{H}}+n_{\mathrm{He}}+n_{\mathrm{e}})+7n_{{\mathrm{H}}_2}}{3(n_{\mathrm{H}}+n_{\mathrm{He}}+n_{\mathrm{e}})+5n_{{\mathrm{H}}_2}},\end{array}$$(88)

which approaches 5/3 for atomic gas and 7/5 for H₂ molecular gas⁹ Integration of Eq. (87) allows us to represent chemical abundances and cooling rates in non-equilibrium as a function of time-dependent temperature and enables us to compare those with ones obtained by assuming instantaneous chemical equilibrium.

Isochoric simulations of cooling gas are performed for two different densities. The first experiment with n_Hnuc = 10⁶ m⁻³ starts from a temperature T = 10⁸ K and assumes that the gas is in chemical equilibrium initially. Over the course of evolution, cooling drives the gas towards lower temperatures up to the point where thermal equilibrium is reached and a new chemical equilibrium is established. Results of the computations for the full network and for the reduced networks (no dust and no cosmic rays) are displayed in the Figs. 13–15. Simulations have been tracked to the final equilibrium state with associated thermal equilibrium temperatures of ∼69 K and ∼74 K for the full network and dust-free network, respectively. The final equilibrium temperature for the CR-free network lies below 10 K outside of the presented temperature regime.

The most significant non-equilibrium effect is the development of ionization lags because the recombination of many ions drags behind cooling. It means that higher ionization stages of species exist down to lower temperatures compared to equilibrium. Such ionization lags are not an exclusive phenomenon of metal species but also appear for hydrogen and helium. This is illustrated in Fig. 13, which compares the computed ionization fractions for the selected ions HeIII, CIV, OVII, and FeV with their corresponding equilibrium fractions for the full network. All of the ions shown have comparatively long recombination times and therefore persist down to lower temperatures in non-equilibrium than predicted for equilibrium. In the case of HeIII, for instance, the ionization fraction does not fall below 1% before ∼1000 K whereas in equilibrium this temperature is ∼5 × 10⁴ K. The occurrence of the double peak for CIV at T ≈ 2.5 × 10⁴ K and T ≈ 1.1 × 10⁵ K is a further characteristic feature that occurs for metal ions where di-electronic recombination dominates radiative recombination in the temperature range between the two appearing maxima. In general, all the computed non-equilibrium ionization curves are qualitatively similar to those obtained in previous studies (e.g., Schmutzler & Tscharnuter 1993; Gnat & Sternberg 2007). The results were obtained for a fixed metallicity of Z = 0.96 Z_⊙. For higher (lower) metallicities, cooling times are shorter (longer), and non-equilibrium ionization lags are expected to be greater (smaller), which has implications for the cooling efficiency (Oppenheimer & Schaye 2013a).

The non-equilibrium cooling function is illustrated in Fig. 14. Results for the full network are discussed first, which are depicted as black lines in Fig. 14. Departures from the equilibrium become apparent for temperatures T ≲ 5 × 10⁵ K. This is the point when cooling proceeds more rapidly than recombination. At higher temperatures the cooling timescale is generally large so that the gas has sufficient time to approximately adjust chemical equilibrium yielding a close-to-equilibrium cooling rate. Down to T ≈ 1.6 × 10⁴ K non-equilibrium cooling (solid line) is suppressed compared to equilibrium cooling (dashed line) by a maximum factor of three appearing near the hydrogen Lyα peak at about 2 × 10⁴ K. The distinct H Lyα bump and the unobtrusive bumps of the C, O, and Ne metal contributions present in the equilibrium cooling function are smeared out in non-equilibrium because of the larger overlap of generally broader ion distributions. The broadening of ion distributions due to ionization lags as seen in Fig. 13 leads to the observed suppression of cooling in this temperature regime. This is because the most effective ion coolants at a given temperature are less abundant than in ionization equilibrium.

On the other hand, non-equilibrium cooling is enhanced for T ≲ 1.6 × 10⁴ K. At T ≈ 10⁴ K the non-equilibrium rate is nearly two orders of magnitude higher than the equilibrium rate. In the equilibrium case metal cooling drops significantly around 1.6 × 10⁴ K owing to a rapid decline of the density of free electrons before molecular cooling becomes important at even lower temperatures. In the non-equilibrium case, due to the recombination delay, a high electron number density exists down to ∼100 K as demonstrated in Fig. 15 (top panel). As a consequence metal cooling is enhanced by roughly two orders of magnitude compared with equilibrium. At the same time, molecular cooling is substantially diminished in non-equilibrium since the formation of H₂ and other important molecular coolants takes too long and thus cannot keep up with cooling. The increase in metal cooling over-compensates the diminished molecular cooling resulting in an overall enhancement compared to equilibrium except in the close vicinity of ∼1000 K where H₂O cooling peaks for the full network under equilibrium conditions (see Fig. 9).

The T-dependence of the H₂ fraction in non-equilibrium is illustrated in Fig. 15 (bottom panel). At any temperature the non-equilibrium H₂ fraction (solid line) remains below the equilibrium fraction (dashed line) with up to three orders of magnitude difference in the interval 100 K–1000 K. At late times, when thermal equilibrium has been established, the H₂ fraction takes on its equilibrium value at T ≈ 69 K (in Fig. 15 this is when the black solid line meets the dashed line). Just before equilibrium settles, the final equilibrium temperature is slightly undershot (in Fig. 15 this is when the curve has a vertical tangent for the first time) as a result of non-equilibrium cooling effects.

In Sect. 4.3 we observed that N and Mg give only minor contributions to the total metal cooling, and that Ne has some significant contribution only near T ≈ 5 × 10⁵ K. To study their neglect on the non-equilibrium cooling curve, the gas cooling experiment is repeated for a network without the elements N, Mg, and Ne. The result is represented by the crosses in Fig. 14. The deviations to the full network are largest near T ≈ 5 × 10⁵ K due to missing Ne cooling but, in general, may be on an acceptable level. When including Ne the discrepancies become rather small, which means that N and Mg coolants may be safely ignored.

The non-equilibrium cooling curves obtained for the reduced networks are virtually indistinguishable from the respective cooling curve for the full network. Neither the absence of dust nor the absence of cosmic rays has a remarkable impact on the non-equilibrium cooling efficiency. This is because of the dominant role played by ionization lags and the relative unimportance of dust and cosmic rays in boostin molecule formation in cooling gas starting from a highly ionized state. The non-equilibrium H₂ fraction for both the dust-free network (yellow solid line in Fig. 15) and the CR-free network (blue solid line in Fig. 15) is some orders of magnitude below the corresponding equilibrium fraction (dashed lines in Fig. 15) down to 100 K. Like for the full network, metal cooling is substantially enhanced in the temperature range 100 K–10⁴ K due to an increased number density of free electrons. In particular, the non-equilibrium electron fraction for the CR-free network (blue solid line in Fig. 15) is up to ten orders of magnitude larger compared to equilibrium (blue dashed line in Fig. 15). For example, at T = 1000 K the non-equilibrium value is n_e/n_Hnuc ≈ 0.3 whereas in equilibrium, n_e/n_Hnuc ≈ 3 × 10⁻¹¹. The T-dependence of the non-equilibrium electron fraction is very alike for all three networks down to 100 K with slightly larger values for the dust-free network because of the missing ion neutralization on dust grains. Near ∼3000 K the non-equilibrium cooling function for the CR-free network falls below the equilibrium function, which is dominated near this temperature by relatively strong H₂ cooling.

In a second experiment the cooling down of high-density gas with n_Hnuc = 10¹⁰ m⁻³ and initial temperature of T = 10⁵ K was studied. Simulations were performed only for the full chemical network including dust and cosmic ray effects. Two different cases are considered that differ in their initial conditions. Case A starts from chemical equilibrium for 10⁵ K. Case B starts with a chemical composition that corresponds to chemical equilibrium at a lower temperature of 100 K. One may think here of a situation where gas is impulsively heated to a high temperature as it occurs in the interaction with high Mach number shocks. Because the gas in case B is predominantly in molecular form at the beginning, the initial phase of evolution is expected to differ significantly from case A. Both cases are contrasted with a solution where equilibrium abundances are imposed at all times. The simulations are stopped when the time-dependent temperature falls below 10 K, which is larger than the thermal equilibrium value of ∼4.47 K.

Figure 16 subsumes the results of these calculations displaying the time evolution of temperature (top panel) in units of days, and shows the net cooling rate (panel 2), electron fraction (panel 3), and H₂ fraction (bottom panel) as a function of temperature for all cases. There is a one-to-one relation between time and temperature, which can be converted into each other. In case A the temperature starts to deviate significantly from the equilibrium solution after one hundred days. Subsequently, the non-equilibrium temperature (black solid line) remains substantially below the equilibrium temperature (dashed line) down to ∼80 K, which is reached after ∼6 × 10⁵ d. For instance, when T = 100 K at t ≈ 6 × 10⁴ d, the equilibrium solution value at the same time is T ≈ 10⁴ K. The different temperature histories are explained with the different cooling rates in the course of evolution. The cooling functions are illustrated in panel 2 of Fig. 16. As in the previous 10⁸ K experiment, the cooling in case A is affected by ionization lags, which increase the electron density at lower temperatures (see panel 3 of Fig. 16) and lead to stronger metal (and net) cooling between T ≈ 4500 K and T ≈ 1.6 × 10⁴ K compared with equilibrium. Consequently, the temperature decreases more rapidly in this temperature regime. Figure 16 (panel 3) demonstrates that the electron fraction stays above 10⁻² down to T ≈ 200 K whereas in the equilibium case it is of the order of 10⁻⁶ for T ≲ 10⁴ K. The H₂ fraction is diminished (Fig. 16, bottom panel), therefore the contribution of molecular cooling to the net cooling is reduced compared to equilibrium.

The results for case B are depicted as red lines in Fig. 16. Starting from 10⁵ K the temperature first falls below case A because the initially high H₂ fraction provides more molecular cooling and net cooling compared to both case A and the equilibrium case. This initial phase of enhanced net cooling is a rather short-term phenomenon and lasts for less than ∼16 d because rapid dissociation of molecules sets in. The dissociation boost halts at t ≈ 16.9 d (T ≈ 4 × 10⁴ K), which is striking in Fig. 16 as the sudden drop in the cooling rate and the sharp decline in the H₂ fraction. The liberated atomic hydrogen is rapidlyionized and this produces a similarly sharp increase in the electron fraction. From t ≈ 300 d to t ≈ 6 × 10⁵ d the temperature in case B exceeds the value for case A because of a lower net cooling rate between T ≈ 100 K and T ≈ 2 × 10⁴ K. The weaker cooling for case B in this temperature regime reflects the fact that ionization lags are less important, and atomic line cooling is closer to equilibrium. The H₂ fraction is higher than in case A but remains relatively unimportant except in the very initial phase ofevolution.

6. Conclusions

In this paper a chemistry and cooling module, NC^​2M, for the multiphysics, adaptive-grid, and MPI-parallelized code NIRVANA has been presented. The NC^​2M has been checked successfully both by computing equilibrium solutions and performing time-dependent, non-equilibrium cooling gas simulations with an isochoric equation of state. NIRVANA now allows for more realistic simulations of interstellar gas by explicitly following non-equilibrium chemistry and thermal processes. The network of the NC^​2M includes chemical schemes for the important coolants H₂, CO, H₂O, and OH and covers a fully-fledged ionization model for the elements H, D, He, C, N, O, Ne, Mg, Si, and Fe. It accounts for dust catalysis (H₂ formation, ion neutralization) and treats cosmic ray effects (ionization, molecule dissociation). The implementation of cooling processes is state of the art. For example, atomic data and electron excitation coefficients for metal line cooling are based on the up-to-date Chianti database version 8. The NC^​2M is seriously applicable within a wide range of temperatures and densities and can be used without restriction from a technical point of view. In extended mode it evolves 121 species linked by 426 chemical reactions. However, exploitation of the full power of the NC^​2M in multidimensional, gas-dynamical simulations is still extremely challenging in terms of computational costs. A reduction of complexity seems possible by neglecting the elements N, Mg, and Ne since N and Mg are minor coolants relative to other metals and Ne essentially contributes to cooling only near T = 5 × 10⁵ K.

Nevertheless, NC^​2M suffers from several weaknesses. The first important aspect concerns the opaqueness of the gas-dust mixture. At sufficiently high gas densities, atomic fine-structure lines and molecular rotovibrational lines become optically thick and are subject to dust absorption as well. In the NC^​2M atomic line cooling is treated in the optically thin approximation. Rotovibrational line cooling crudely respects the optically thick regime by either adopting a local-density-dependent attenuation factor in the case of H₂ or by simply setting an upper threshold for the optical depth parameter in the cooling models for CO and H₂O. A more proper treatment of an optically thick, dusty medium would require the solution of a complex radiative transfer problem, which is currently beyond our capabilities. Second, in the present stage of development the NC^​2M ignores photoionization and photodissociation by incident radiation. The interstellar background radiation field as cataloged by Black (1987) possesses an ionizing UV component from starlight. In regions of unshielded gas, photochemical reactions can therefore have a strong impact on the ionization level and molecular chemistry in the low-temperature regime which, in turn, affects the cooling properties. Photoionization of hydrogen, helium, and metals cause a reduction of their cooling efficiencies because fewer bound electrons mean fewer line transitions. Third, dust in the NC^​2M is considered as a passive constituent. There is no time-dependent treatment of dust destruction (thermal sputtering, sublimation) and dust formation processes (accretion, coagulation). Since H₂ production and ion neutralization are potentially influenced by dust, such destruction and formation processes become important if their timescales are short compared to other relevant timescales. For instance, in the gas cooling simulations starting with a high temperature (Sect. 5), dust would be destroyed by thermal sputtering in a time shorter than the cooling time. At the same time, dust would not reform at an equal rate since dust growth by accretion is rather inefficient at high temperatures (Draine & Salpeter 1979). Consequently, the influence of dust-catalytic reactions like H₂ formation are likely to be overestimated in these experiments by assuming a fixed dust component. Future development efforts aim to remove at least some of these weaknesses.

---

Acknowledgments

I would like to thank the unknown referee for her or his dedicated report and many useful comments, which helped to improve the paper.

References

Appendix A: The chemical network

All Tables

All Figures