EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 620 (December 2018) A&A, 620 (2018) A81 Full HTML
Free Access
Issue
A&A
Volume 620, December 2018
Article Number A81
Number of page(s) 27
Section Numerical methods and codes
DOI https://doi.org/10.1051/0004-6361/201834033
Published online 04 December 2018

© 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 CH2 and (hydro-)oxygen species like OH and H2O are consolidated in the form of (artificial) composite species CHx and OHx. 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 H2 and HD but ignores other molecular coolants like CO and H2O (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 × 104 K), and combining it with the ionization and cooling model of Oppenheimer & Schaye (2013a) appropriate in the high-temperature regime (T ≳ 2 × 104 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 H2 formation through the well-studied H−/ 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 ns of different species Xs, s = 1…Ns, of the form

(1)

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

between reactants Xs ∈ R(r) leading to products Xs ∈ P(r) with R(r) and P(r) being subsets of .

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 H2, , , HD, C2, O2, , CH, CH+, CH2, , , OH, OH+, CO, CO+, H2O, H2O+, H3O+, 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 ≈ 1010 m−3, 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 H2 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

(2)

with specified temperature range [5.5 × 103 K, 105 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 = 109 K. Although k becomes inaccurate for T ≫ 105 K, it will not have a large impact on the overall chemical evolution since the total reaction rate, knC+nH, for this process is negligible at such high temperatures owing to small number densities nC+ and nH. 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),

radiative recombination (RR),

and di-electronic recombination (DR),

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-published1 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 Fe13+, 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

for hydrogen CT ionization of concerned species X = D, He, C, N, O, Mg, Mg+, Si, Si+, Fe, Fe+, H2, O2, OH, and H2O;

for helium CT ionization of concerned species X = D, Si, H2, O2, OH, and H2O;

for hydrogen CT recombination of concerned species X = D, He, C, C+, C2+, C3+, N, N+, N2+, N3+, O, O+, O2+, O3+, Ne+, Ne2+, Ne3+, Mg+, Mg2+, Mg3+, Si+, Si2+, Si3+, Fe+, Fe2+, Fe3+, H2, and CO; and

for helium CT recombination of concerned species X = D, C2+, C3+, N+, N2+, N3+, O+, O2+, O3+, Ne+, Ne2+, Ne3+, Mg2+, Mg3+, Si2+, Si3+, Fe2+, and Fe3+. 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 database2, 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. H2 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 H2 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 H2 fractions are produced by pure gas-phase reactions.

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

H2 is then formed by associative detachment with H according to

This two-step H2 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 H2 grain catalysis. In addition to the above fundamental reactions, a couple of further reactions are taken into account involving also the ions and . Some of the reactions contribute to H2 formation like charge transfer from H to (reaction 303) or the dissociative detachment of 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 H2 to H+ (reaction 304).

At sufficiently high gas densities (such as n ≳ 1016 m−3), a second pathway for H2 formation is by ter-molecular reactions. Here, the following three important reactions are taken into account:

where X = H, He, H2.

The H2 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 H2 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

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:H2 abundance ratio relative to the prevalent D:H ratio. At higher temperatures, the analog neutral-neutral reactions

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, D2, deuterated hydronium, D3O+, 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−/ channel of H2 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 H2O. 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 H2O in H2-dominated regions, intermediate molecules and its ionized versions such as OH, OH+, the hydrocarbons CHx, , and the ions HCO+ and HOC+ play a role. The major route to form H2O is initiated by the reactions

forming the hydroxyl ion. This is followed by the sequence

and subsequent dissociative recombinations

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

Important routes to form CO are by the reactions

followed by

and by the reactions

where the formation and destruction of hydrocarbon radicals up to is governed by numerous reactions not discussed here, but fully included in the NC​2M. Further minor sources of CO are the dissociative recombination

the ion-neutral reactions

and the charge transfer reaction

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, H2, 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 H2,

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,

A summary of all CR-induced reactions included in the NC​2M 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−17 s−1 from Williams et al. 1998 is adopted in order to allow comparison with other work. However, more recent observations of 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,

(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,

(4)

where χH and χ are the ionization energies of hydrogen and the species in mind, respectively, and denotes the outer shell effective electron number of that species given by

(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 H2–H–He gas mixture3 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 H2, the yields of H(2p) excitations (Lyα photons), and the yields of electronic excitations of H2 (Lyman–Werner band photons) are computed as functions of the number density ratio nH2/nH and ionization degree. The net (primary plus secondary) cosmic ray ionization rate for atomic and molecular hydrogen is therefore represented by

(6)

(7)

where the numbers of secondary ionizations are

(8)

(9)

with xH = nH/(2nH2 + nH) and xH2 = nH2/(2nH2 + nH). Dalgarno et al. (1999) state that the expressions for and are accurate for ionization degrees up to 10%. For a neutral H–He gas mixture (xH = 1, ne = 0) , whereas for a neutral H2-He gas mixture (xH2 = 0.5, ne = 0) . 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 , and 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–H2 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 H2 by a 30 eV electron. Accordingly, the photoreaction rate RX in units s−1 for a species X is given by

(10)

where pX 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 pX a pure neutral H2 gas with was assumed, and yields of 0.3 and 0.12 for the dominating excitations to the (including cascading contribution from higher excited states) and C1Πu electronic states of H2 were adopted. For the NC​2M, the Gredel et al. outcome has been recalibrated to the situation of a partially ionized H–H2–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

(11)

This suggest the following modified photoreaction rate:

(12)

where qX 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, , is given by

(13)

We note that RX now depends on the ionization degree because the yields YLyW and YLyα do so. For a neutral H2-He gas mixture (xH2 = 0.5, nHe = 0.2nH2, ne = 0), the rate reduces to

(14)

which is close to the Gredel et al. (1989) result where . The efficiencies pX and qX for all CR-induced photoreactions in the NC​2M are summarized in Table A.2.

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

(15)

This approaches the original Gredel et al. (1987) result, RCO = pCO1.55ζH2, if electron degradation is assumed to take place in a neutral H2 gas where YLyW ≈ 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​2M adheres to the fundamental process of H2 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 nd(a) (a: grain radius). According to Draine (2002) the total dust cross-sectional area per unit volume for this model is

(16)

with nHnuc being the number density of H nuclei4. For a dust-to-gas ratio other than the average value in the Milky Way, it is common to scale (ndσd)tot linearly with the metallicity Z, hence, (ndσd)tot = 6.7 × 10−25m2 ⋅ nHnuc(Z/Z⊙).

The rate of H2 formation according to reaction 700, , can be expressed as (see, e.g., Cazaux & Tielens 2002)

(17)

where

(18)

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

(19)

The function ϵH2 describes the surface recombination efficiency, which is the fraction of accreted H atoms that leave the surface as H2. 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

(20)

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

(21)

where vth, HnHπa2SH is the number of H atoms per second landing on a spherical grain with radius a and Nass = 4πa2/W2 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 = 1010m−3 and T = Td = 10 K gives a value of F ≈ 5 × 10−8 monolayers per second. We note that RH2 becomes negligible for higher gas temperatures since vth, HSH ∝ 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, kX, for reactions 701–706 use the fitting formula from Weingartner & Draine (2001b), Eq. (8),

(22)

with the coefficients C0 − C6 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 RH2, and accounts for grain charging expressed by the (dimensionless) charging parameter (see Weingartner & Draine 2001c)

(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 RH2. As stressed by Weingartner & Draine (2001b), the expression for kX likely becomes inaccurate when . It is interesting to note that the limit G → 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 kX 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 ϵ,

(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 H2 formation reactions and cosmic ray heating. Among the cooling processes are chemical cooling due to ionization, recombination, and H2 dissociation, atomic line cooling of all elements, and molecular line cooling of H2, CO, H2O, and OH.

3.1. Chemical heating and cooling

The NC​2M takes into account chemical heating of the gas through the exothermic reactions of H2 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

(25)

Following Omukai (2000) and Hollenbach & McKee (1979) the gas heating rates associated with H2 formation by the different processes, namely, H−/ (reactions 301,303), three-body reactions 313, 314, 317, and grain-catalytic reaction 700 are given by

(26)

(27)

(28)

(29)

(30)

(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 H2 binding energy of 4.48 eV, which is assumed to be stored initially in excitation as well. In the heating rate due to H2 formation on dust grains, Γ700, 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 H2 (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):

(32)

(33)

(34)

(35)

(36)

(37)

(38)

(39)

(40)

In the recombination cooling coefficients Λ2, Λ4, and Λ6, αkBT ≈ 0.75kBT is the mean kinetic energy of the captured electron5. The second term in Λ3 accounts for collisional ionization cooling from the metastable 23S1 He state, which is proportional to (Black 1981), and where C(T)=(1 + T/105 K)−1.

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 11S0, 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 23S1 state to n = 2, 3, 4 triplet states adopting the rate from Cen (1992). The rates can be written as

(41)

(42)

(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

(44)

where the sum runs over all species s = C ...C5+, N ...N6+, O ...O7+, Ne ...Ne9+, Mg ...Mg11+, Si ...Si13+, Fe ...Fe25+. The cooling contribution of an individual species or ion is (the index s is suppressed in the following)

(45)

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

In order to compute the fractions ( for an ion, the statistical equilibrium equations are solved for a prescribed maximum number of atomic levels L subject to the conservation condition . 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 Qij (i ≠ j), it follows for any state i that

(46)

or in vector notation, setting for convenience,

(47)

The constraint is directly incorporated by exchanging with it the first equation of Eq. (46). Equation (47) is then replaced by

(48)

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

(49)

where (, ) is the electron (H, H2) collisional (de-)excitation rate. The Einstein coefficient Aij is defined to be zero for i ≤ j. The excitation rate (i <  j) is related to the de-excitation rate (i >  j) according to

(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 H2-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

(51)

where Υij(T) is the thermally-averaged collision strength for the i → j transition. To obtain 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 Λ/nionne, 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 dependence6 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−4 eV (IR) for the transition 2s22p3 2​P3/2 → 2s22p3 2​P1/2 of NI to 8.26 keV (X-ray) for the transition 3p 2​P3/2 → 1s 2​S1/2 of FeXXVI.

thumbnail Fig. 1.

Influence of parameter L on the cooling rate (top panels) and on the per-ion cooling coefficient (bottom panels) for OIII and FeVI. The dashed line corresponds to the highest possible L with Chianti. The results were obtained for an equilibrium model with nH = 106 m−3.

Table 1.

Recommended numbers of atomic energy levels for metal cooling computations as adopted in the NC​2M (z: ion charge).

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 = 104 K for which theoretical estimates in the low-density limit are given in Osterbrock & Ferland (2005):

(52)

(53)

(54)

The line intensity I(i → j)∝piΔEijAij. 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.

thumbnail Fig. 2.

Ratio of spectral line intensities obtained with the NC​2M (solid lines) and theoretical estimates (dashed lines, Eqs. (52)–(54)) for nH = 104 m−3.

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 H2, CO, H2O, 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. H2

For densities that are not too high (n ≲ 109 m−3) 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 H2 cooling I closely follow the work of Glover & Abel (2008) who considered collisional excitation of rotovibrational levels for collision partners H, H+, He, H2, and e−. The ratio of ortho-H2 to para-H2 is fixed to 3:1 independent of temperature. The cooling rate in an optically thin environment can be written as

(55)

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

(56)

The individual contributions are

(57)

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

(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 H2 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),

(59)

In the low-density limit, ΛH2, LTE ≫ ΛH2, 0, the effective cooling rate is ΛH2 → ΛH2, 0nH2 whereas in the high-density limit ΛH2 approaches ΛH2, LTEnH2. At high density (n ≳ 1016 m−3) 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

(60)

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

At even higher densities (n ≳ 1020 m−3) collision-induced emission (CIE) becomes a further important cooling effect. In the CIE process, deep collisions between H2 molecules or between H2 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

(61)

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

(62)

to approximately account for cooling in the optically thick regime when n ≳ 1022 m−3 (Hirano & Yoshida 2013). Under the assumption of LTE, ΛCIE is expressed by

(63)

summing up cooling contributions for the collision pairs H2–H2, H2–He, and H2–H. The quantity Bν(T) in units Js−1 m−2 Hz−1 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, nH2, and the number density ns 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 H2 molar density and per molar density of collision partner, (in units7 cm−1 amg−2), in the frequency interval and temperature range as specified in Table 2 for the different collision pairs. The specific does not depend on the number densities of collision partners and is related to αs, ν (in units m−1) in Eq. (63) according to

Table 2.

Sources of .

(64)

Equation (63) can then be rewritten as

(65)

with frequency-integrated cooling rate coefficients , 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:

(66)

The fitting coefficients al, b, c, Tu 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.

thumbnail Fig. 3.

CIE cooling fitting curves (lines) for the different collision pairs according to Eq. (66), and data points (symbols) derived from the sources in Table 2.

Table 3.

Fitting coefficients for expression 66.

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 ≳ 1010 m−3: a regime where H2 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

(67)

with the fitting function

(68)

The coefficients Dl, m are tabulated in Lipovka et al. (2005), Table 1. Since the excitation rate coefficients for HD-He collisions and HD–H2 collisions are comparable to the excitation rate for HD–H collisions, as pointed out by Flower & Le Bourlot (2000), an effective number density neff = nH + nH2 + nHe is used in Eq. (68). Lipovka et al. (2005) emphasize the high accuracy of the above approximation for 100 K ≲ T ≲ 2 × 104 K and 106 m−3 ≲ neff ≲ 1014 m−3. 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

(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 H2O

The NC​2M implementation of CO cooling and H2O 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 H2O cooling rates are given by

(70)

where Λs, rot and Λs, vib are cooling rate coefficients due to rotational and vibrational transitions, respectively, and neff, rot, s, neff, vib, s are effective number densities for collisions with H2, H, and e−. In the Neufeld et al. papers, neff, rot, s = neff, vib, s = nH2 assuming H2 dominates collisional excitation of CO and H2O. The cooling rate coefficients are determined from

(71)

for rotational transitions and from

(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 n1/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, n1/2, s, and exponent αs depend, in addition, on an optical depth variable . This variable is defined by

(73)

meaning an effective column density per unit velocity. The functions Λs, rot, LTE, Λs, vib, LTE, n1/2, s, and αs are sampled in Neufeld & Kaufman (1993), Tables 2, 3, and 5, and Neufeld et al. (1995), Tables 1–3, in -space for the range

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

(74)

whereas for vibrational cooling

(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

(76)

Following Hollenbach & McKee (1979), the cooling rate per OH molecule, , is given by

(77)

where

(78)

(79)

and

(80)

The quantity ncr 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),

(81)

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

(82)

Bremsstrahlung losses become significant at T ≈ 106 K where they compare with H and He atomic line cooling, and they dominate all metal cooling processes for T ≳ 107 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 H2 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-H2 abundance ratio. According to Dalgarno et al. (1999) the heating efficiency in a partly ionized H–H2–He gas mixture can be represented by

(83)

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

(84)

(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, ne/nH ≲ 0.1 or ne/nH2 ≲ 0.1, respectively.

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

(86)

To give some numbers, for a neutral H–He (H2-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 (H2-He) gas mixture with ne = 0.1nH (ne = 0.1nH2), 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−6. The equilibrium species abundances are calculated as a function of temperature in the range 10 K ≤ T ≤ 108 K at 801 nodes with ΔlogT = 8.75 × 10−3 for three different densities nH = 104 m−3, 106 m−3, 1010 m−3 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−17 s−1 is adopted. For grain-surface reactions, the dust temperature, Td, and the parameter G appearing in expression (23) for grain charging need to be defined. As in RSO14 we adopt a value of Td = 10 K since the rate coefficient given by Eq. (17) is not very sensitive to the dust temperature in the regime 6 K ≲ Td ≲ 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 H2, CO, and H2O and for the further molecules O2 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 H2, C nuclei density in the case of CO, and O nuclei density in the case of H2O, O2, and OH. Only the low-temperature regime is shown since molecules are effectively absent at T ≳ 2 × 104 K irrespective of density.

thumbnail Fig. 4.

Fractions of H2 (first row), CO (second row), H2O (third row), O2 (fourth row), and OH (fifth row) for different densities (from left to right). Equilibrium solutions for the full network (solid), dust-free network (dashed), CR-free network (dash-dotted), and pure gas-phase network (dotted) are shown.

Molecular abundances for the full network case are depicted as solid curves in Fig. 4. It is first noted that the obtained H2 fraction is larger than in RSO14 because dust-catalytic H2 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 assumption8. At moderate densities the curves of the O2 and OH fractions closely follow that of H2O. The H2, CO, and H2O fractions generally grow with increasing density.

At the highest density, nH = 1010 m−3, most hydrogen is in molecular form and the H2 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 H2. This cut-off is also striking as sharp transitions in the other molecules, since H2 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−2 associated with the increase in H2O fraction. Evidently, more oxygen becomes bound in H2O, 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 H2 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 nH = 1010 m−3) hydrogen is predominantly bound in H2 followed by a smaller fraction bound in H2O. The increased amount of H2 compared to the full network case is the result of missing cosmic ray ionization and CR-induced photodissociation of H2. The fraction of atomic hydrogen is very low in this temperature regime: nH/nHnuc <  10−6 (< 10−8, < 10−12) at nH = 104 m−3 (nH = 106 m−3, nH = 1010 m−3) compared to nH/nHnuc >  0.8 (> 0.18, > 5 × 10−5) 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 nH = 1010 m−3). Other C-bearing molecules like CH2, for example, exist only in small amounts with a fraction of nCH2/nCnuc <  10−10 for this molecule. There are three distinct temperature ranges where the H2O fraction is nearly constant. At a density nH = 104 m−3 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 O2 fraction drops by many order of magnitudes in the third temperature regime where the H2O fraction reaches its largest value. This O2 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 H2, has a substantial impact on the H2 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 H2 formation. Therefore, despite the lack of dust, the H2 fraction at nH = 1010 m−3 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 H2O (O2) 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 H2 fraction is up to two orders of magnitude smaller than in the full network case. At the same time the H2O 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 H2O fraction. Below a certain temperature both O2 and OH get substantially depleted independent of density because most oxygen is bound in H2O.

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).

thumbnail Fig. 5.

Ionization structure of C for different densities. Equilibrium solution for the full network (solid), dust-free network (dashed), and CR-free network (dash-dotted).

thumbnail Fig. 6.

Analog to Fig. 5 for O.

thumbnail Fig. 7.

Analog to Fig. 5 for Si.

thumbnail Fig. 8.

Analog to Fig. 5 for Fe.

At temperatures T ≲ 104 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 nH = 1010 m−3. For instance, the fraction of FeII grows by three orders of magnitude in the absence of dust (dashed lines) from ∼4 × 10−4 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 nH = 1010 m−3 where a significant amount of carbon is bound in CO. Molecule formation of O2 and H2O 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, ( for heating terms, for net cooling), where Λ* has units J m3 s−1 (=1013 erg cm3 s−1). Calculations are shown for the full network at a density of nH = 106 m−3 if not otherwise stated.

Figure 9 shows equilibrium cooling functions for the molecules H2, CO, H2O, and OH resulting from line cooling. The H2 cooling function, in addition, contains contributions from reaction kinetics (dissociation cooling and formation heating). Also, CIE cooling is included in the H2 function which is, however, unimportant at the gas densities explored here. The H2 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 H2 formation on dust grains. The second heating window is a consequence of fixing the ortho-to-para H2 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. H2 cooling dominates molecular cooling except in the vicinity of T ≈ 1000 K where cooling due to H2O 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 H2 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.

thumbnail Fig. 9.

Molecular cooling in chemical equilibrium for the full network.

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 × 104 K ≲ T ≲ 3.5 × 105 K, and by Fe for T ≳ 6 × 105 K. Around T ≈ 5 × 105 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 106 K (dashed line in Fig. 10, bottom panel).

thumbnail Fig. 10.

Metal cooling in chemical equilibrium for the full network. Top panel: contribution of individual oxygen ionization stages. Bottom panel: contribution of different metals. The dashed line represents the total metal cooling when the number of energy levels is restricted to 5 (L = 5) in all ion computations.

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 nH = 106 m−3. 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 H2 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 × 104 K, which is related to cooling processes of hydrogen and helium. Metal cooling then dominates over the large interval from T ≈ 3 × 104 K to T ≈ 107 K. Cooling due to Bremsstrahlung becomes unimportant before ∼106 K but dominates above ∼107 K.

thumbnail Fig. 11.

Total cooling function and its contributions from different cooling or heating agents in chemical equilibrium for the full network. Dash-dotted lines indicate heating.

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 H2O) and above T ≈ 1000 K (H2). 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 nH = 104 m−3 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 nH = 1010 m−3 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 nH = 106 m−3 case. The reduction comes about because, with increasing gas density, molecular cooling approaches the LTE regime, which lowers the cooling rate per , that is, Λ*.

thumbnail Fig. 12.

Impact of different gas densities and the absence of dust or cosmic rays on the net cooling function in chemical equilibrium. Dash-dotted lines indicate net heating, Λ* <  0.

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

(87)

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

(88)

which approaches 5/3 for atomic gas and 7/5 for H2 molecular gas9 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 nHnuc = 106 m−3 starts from a temperature T = 108 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.

thumbnail Fig. 13.

Ionization fractions of HeIII, CIV, OVII, and FeV obtained for the full network in equilibrium (dashed lines) and non-equilibrium (solid lines) in the 108 K cooling gas experiment.

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 × 104 K. The occurrence of the double peak for CIV at T ≈ 2.5 × 104 K and T ≈ 1.1 × 105 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 × 105 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 × 104 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 × 104 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.

thumbnail Fig. 14.

Cooling curve in equilibrium (dashed lines) and in non-equilibrium (solid lines) for the full network (black), dust-free network (yellow), and CR-free network (blue). For T ≳ 100 K the solid lines representing the different network cases nearly coincide. Crosses: Non-equilibrium cooling curve for a network without the metals N, Mg and Ne.

On the other hand, non-equilibrium cooling is enhanced for T ≲ 1.6 × 104 K. At T ≈ 104 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 × 104 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 H2 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 H2O cooling peaks for the full network under equilibrium conditions (see Fig. 9).

thumbnail Fig. 15.

Electron fraction (top panel) and molecular hydrogen fraction (bottom panel) as a function of temperature in equilibrium (dashed lines) and in non-equilibrium (solid lines) for the full network (black), dust-free network (yellow), and CR-free network (blue).

The T-dependence of the H2 fraction in non-equilibrium is illustrated in Fig. 15 (bottom panel). At any temperature the non-equilibrium H2 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 H2 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 × 105 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 × 105 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 H2 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–104 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 ne/nHnuc ≈ 0.3 whereas in equilibrium, ne/nHnuc ≈ 3 × 10−11. 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 H2 cooling.

In a second experiment the cooling down of high-density gas with nHnuc = 1010 m−3 and initial temperature of T = 105 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 105 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 H2 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 × 105 d. For instance, when T = 100 K at t ≈ 6 × 104 d, the equilibrium solution value at the same time is T ≈ 104 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 108 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 × 104 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−2 down to T ≈ 200 K whereas in the equilibium case it is of the order of 10−6 for T ≲ 104 K. The H2 fraction is diminished (Fig. 16, bottom panel), therefore the contribution of molecular cooling to the net cooling is reduced compared to equilibrium.

thumbnail Fig. 16.

Time dependence of the temperature (top panel) and the net cooling rate (panel 2), electron fraction (panel 3), and molecular hydrogen fraction (bottom panel) as a function of temperature for case A (black solid lines), case B (red lines), and for the equilibrium solution (dashed lines).

The results for case B are depicted as red lines in Fig. 16. Starting from 105 K the temperature first falls below case A because the initially high H2 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 × 104 K), which is striking in Fig. 16 as the sudden drop in the cooling rate and the sharp decline in the H2 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 × 105 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 × 104 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 H2 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 H2, CO, H2O, 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 (H2 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 × 105 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 H2 or by simply setting an upper threshold for the optical depth parameter in the cooling models for CO and H2O. 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 H2 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 H2 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.

3

With a fixed helium fraction of nHe = 0.1(2nH2 + nH).

4

This is a factor 6.7 larger than the total cross section used in RSO14, which was based on the dust UV absorption cross section per H nucleus mentioned in Krumholz et al. (2011).

5

Generally, considering Spitzer (1978), the quantity α (the quotient χ1/ϕ1 of the functions χ1 and ϕ1 in Spitzer (1978)) is a slowly varying function of T taking on values between 0.88 and 0.69 for T between 100 K and 105 K. Here, following Cen (1992), a value of α = 0.75 is adopted.

6

The temperature range has been restricted to the regime where the ion fraction of the element is larger than 10−5. This is an optimization feature of the code and can be switched off. It utilizes the fact that for a small ion fraction neighboring ionization stages dominate the cooling of the element (see also Fig. 10, top panel).

7

1 amg (amagat) equals 44.615 mol m−3.

8

There are still minor discrepancies probably caused by some differences in the rate coefficients. In particular, if we drop the T−1/2-dependence from k506 of the important reaction 506 (this dependence is somewhat hidden in the reference Petuchowski et al. 1989), excellent agreement is found as well.

9

Generally, the adiabatic index of H2 gas shows a dependence on temperature for T ≲ 400 K where γ >  7/5, and the caloric equation of state is more complex in this temperature regime (see, e.g., Boley et al. 2007, Fig. 1).

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

  1. Abel, T., Bryan, G. L., & Norman, M. L. 2002, Science, 295, 93 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  2. Adams, N. G., & Smith, D. 1976a, Int. J. Mass Spectrom. Ion Phys., 21, 349 [NASA ADS] [CrossRef] [Google Scholar]
  3. Adams, N. G., & Smith, D. 1976b, J. Phys. B, 9, 1439 [NASA ADS] [CrossRef] [Google Scholar]
  4. Adams, N. G., Smith, D., & Grief, D. 1978, Int. J. Mass Spectrom. Ion Phys., 26, 405 [NASA ADS] [CrossRef] [Google Scholar]
  5. Adams, N. G., Smith, D., & Paulson, J. F. 1980, J. Chem. Phys., 72, 288 [Google Scholar]
  6. Adams, N. G., Smith, D., & Millar, T. J. 1984, MNRAS, 211, 857 [NASA ADS] [Google Scholar]
  7. Aldrovandi, S. M. V., & Pequignot, D. 1973, A&A, 25, 137 [NASA ADS] [Google Scholar]
  8. Alge, E., Adams, N. G., & Smith, D. 1983, J. Phys. B, 16, 1433 [Google Scholar]
  9. Andreazza, C. M., & Singh, P. D. 1997, MNRAS, 287, 287 [NASA ADS] [CrossRef] [Google Scholar]
  10. Anicich, V. G., Futrell, J. H., Huntress, W. T., et al. 1975, Int. J. Mass Spectrom. Ion Phys., 18, 63 [NASA ADS] [CrossRef] [Google Scholar]
  11. Anicich, V. G., Huntress, W. T., & Futrell, J. H. 1976, Chem. Phys. Lett., 40, 233 [NASA ADS] [CrossRef] [Google Scholar]
  12. Arnaud, M., & Raymond, J. 1992, ApJ, 398, 394 [NASA ADS] [CrossRef] [Google Scholar]
  13. Azatyan, V. V., Aleksandrov, E. N., & Troshin, A. F. 1975, Kinet. Cat., 16, 306 [Google Scholar]
  14. Badnell, N. R. 2006, ApJS, 167, 334 [NASA ADS] [CrossRef] [Google Scholar]
  15. Barinovs, G., & van Hemert, M. C. 2006, ApJ, 636, 923 [NASA ADS] [CrossRef] [Google Scholar]
  16. Barlow, S. G. 1984, PhD Thesis, Univ. Colorado [Google Scholar]
  17. Baulch, D. L., Cobos, C. J., Cox, R. A., et al. 1992, J. Phys. Chem. Ref. Data, 21, 411 [NASA ADS] [CrossRef] [Google Scholar]
  18. Black, J. H. 1981, MNRAS, 197, 553 [NASA ADS] [Google Scholar]
  19. Black, J. H. 1987, Interstellar Processes, (Berlin: Springer Verlag), Astrophysics and Space Science Library, 134, 731 [NASA ADS] [CrossRef] [Google Scholar]
  20. Boley, A. C., Hartquist, T. W., Durisen, R. H., et al. 2007, ApJ, 660, L175 [NASA ADS] [CrossRef] [Google Scholar]
  21. Borysow, A. 2002, A&A, 390, 779 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. Borysow, A., Jorgensen, U. G., & Fu, Y. 2001, JQSRT, 68, 235 [Google Scholar]
  23. Bray, I., Burgess, A., Fursa, D. V., et al. 2000, A&AS, 146, 481 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Bruhns, H., Kreckel, H., Miller, K. A., et al. 2010, Phys. Rev. A, 82, 042708 [NASA ADS] [CrossRef] [Google Scholar]
  25. Burgess, A., & Tully, J. A. 1992, A&A, 254, 436 [NASA ADS] [Google Scholar]
  26. Cazaux, S., & Tielens, A. G. G. M. 2002, ApJ, 575, L29 [NASA ADS] [CrossRef] [Google Scholar]
  27. Cazaux, S., & Tielens, A. G. G. M. 2004, ApJ, 604, 222 [NASA ADS] [CrossRef] [Google Scholar]
  28. Cazaux, S., & Spaans, M. 2004, ApJ, 611, 40 [NASA ADS] [CrossRef] [Google Scholar]
  29. Cen, R. 1992, ApJS, 78, 341 [NASA ADS] [CrossRef] [Google Scholar]
  30. Cohen, N., & Westberg, K. R. 1979, J. Phys. Chem., 83, 46 [CrossRef] [Google Scholar]
  31. Cohen, N., & Westberg, K. R. 1983, J. Phys. Chem. Ref. Data, 12, 531 [NASA ADS] [CrossRef] [Google Scholar]
  32. Cravens, T. E., & Dalgarno, A. 1978, ApJ, 219, 750 [NASA ADS] [CrossRef] [Google Scholar]
  33. Croft, H., Dickinson, A. S., & Gadea, F. X. 1999, MNRAS, 304, 327 [NASA ADS] [CrossRef] [Google Scholar]
  34. Dalgarno, A., & Black, J. H. 1976, Rep. Progr. Phys., 39, 573 [NASA ADS] [CrossRef] [Google Scholar]
  35. Dalgarno, A., & Lepp, S. 1987, Astrochemistry (Dordrecht: Reidel) [Google Scholar]
  36. Dalgarno, A., Du, M. L., & You, J. H. 1990, ApJ, 349, 675 [NASA ADS] [CrossRef] [Google Scholar]
  37. Dalgarno, A., Yan, M., & Liu, W. 1999, ApJS, 125, 237 [NASA ADS] [CrossRef] [Google Scholar]
  38. Dean, A. J., Davidson, D. F., & Hanson, R. K. 1991, J. Phys. Chem., 95, 183 [CrossRef] [Google Scholar]
  39. Del Zanna, G., Dere, K. P., Young, P. R., et al. 2015, A&A, 582, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  40. Donahue, M., & Shull, J. M. 1991, ApJ, 383, 511 [NASA ADS] [CrossRef] [Google Scholar]
  41. Dove, J. E., Rusk, A. C. M., Cribb, P. H., et al. 1987, ApJ, 318, 379 [Google Scholar]
  42. Dubernet, M. L., Gargaud, M., & McCarroll, R. 1992, A&A, 259, 373 [NASA ADS] [Google Scholar]
  43. Draine, B. T. 2002, in The Cold Universe: Saas-Fee Advanced Course, eds. D. Pfenniger, & Y. Revaz (Springer Verlag), 32, 213 [Google Scholar]
  44. Draine, B. T., & Salpeter, E. E. 1979, ApJ, 231, 77 [NASA ADS] [CrossRef] [Google Scholar]
  45. Fairbairn, A. R. 1969, Proc. R. Soc. A, 312, 207 [NASA ADS] [CrossRef] [Google Scholar]
  46. Faure, A., Gorfinkiel, J. D., & Tennyson, J. 2004, MNRAS, 347, 323 [NASA ADS] [CrossRef] [Google Scholar]
  47. Federer, W., Villinger, H., Howorka, F., et al. 1984, Phys. Rev. Lett., 52, 2084 [CrossRef] [Google Scholar]
  48. Fehsenfeld, F. C. 1976, ApJ, 209, 638 [NASA ADS] [CrossRef] [Google Scholar]
  49. Ferland, G. J., Peterson, B. M., Horne, K., et al. 1992, ApJ, 387, 95 [NASA ADS] [CrossRef] [Google Scholar]
  50. Ferland, G. J., Porter, R. L., van Hoof, P. A. M., et al. 2013, Rev. Mex. Astron. Astrofis., 49, 137 [Google Scholar]
  51. Field, R. J., Körõs, E., Noyes, R. M. J., et al. 1972, J. Am. Soc., 94, 8649 [CrossRef] [Google Scholar]
  52. Field, D., Adams, N. G., & Smith, D. 1980, MNRAS, 192, 1 [NASA ADS] [Google Scholar]
  53. Flower, D. R., Le Bourlot, J., Pineau des Forêts, G., et al. 2000, MNRAS, 314, 753 [NASA ADS] [CrossRef] [Google Scholar]
  54. Frank, P. 1986, Proc. Int. Symp. Rarefied Gas Dyn., 2, 422 [Google Scholar]
  55. Fryxell, B., Olson, K., Ricker, P., et al. 2000, ApJS, 131, 273 [NASA ADS] [CrossRef] [Google Scholar]
  56. Furlanetto, S. R., & Stoever, S. J. 2010, MNRAS, 404, 1869 [NASA ADS] [Google Scholar]
  57. Geppert, W. D., Thomas, R. D., Ehlerding, A., et al. 2005, J. Phys. Conf. Ser., 4, 26 [NASA ADS] [CrossRef] [Google Scholar]
  58. Gerlich, D. 1982, in Symp. on Atomic and Surface Physics, eds. W. Lindinger, F. Howorka, & T. D. Märk (Dordrecht: Kluwer), 304 [Google Scholar]
  59. Gerlich, D., & Horning, S. 1992, Chem. Rev., 92, 1509 [NASA ADS] [CrossRef] [Google Scholar]
  60. Glover, S. C. O., & Jappsen, A.-K. 2007, ApJ, 666, 1 [NASA ADS] [CrossRef] [Google Scholar]
  61. Glover, S. C. O., & Abel, T. 2008, MNRAS, 388, 1627 [NASA ADS] [CrossRef] [Google Scholar]
  62. Glover, S. C. O., & Clark, P. C. 2012, MNRAS, 421, 116 [NASA ADS] [Google Scholar]
  63. Glover, S. C. O., Federrath, C., Mac Low, M.-M., et al. 2010, MNRAS, 404, 2 [NASA ADS] [Google Scholar]
  64. Gnat, O., & Sternberg, A. 2007, ApJS, 168, 213 [NASA ADS] [CrossRef] [Google Scholar]
  65. Gnat, O., & Ferland, G. J. 2012, ApJS, 199, 20 [NASA ADS] [CrossRef] [Google Scholar]
  66. Grassi, T., Bovino, S., Schleicher, D. R. G., et al. 2014, MNRAS, 439, 2386 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  67. Gray, C. 2002, Rose-Hulman Undergraduate Math. J., 3 [Google Scholar]
  68. Gray, W. J., Scannapieco, E., & Kasen, D. 2015, ApJ, 801, 107 [NASA ADS] [CrossRef] [Google Scholar]
  69. Gredel, R., Lepp, S., & Dalgarno, A. 1987, ApJ, 323, L137 [NASA ADS] [CrossRef] [Google Scholar]
  70. Gredel, R., Lepp, S., Dalgarno, A., et al. 1989, ApJ, 347, 289 [NASA ADS] [CrossRef] [Google Scholar]
  71. Grevesse, N., & Sauval, A. J. 1998, Space Sci. Rev., 85, 161 [NASA ADS] [CrossRef] [Google Scholar]
  72. Guberman, S. L. 1995, J. Phys. Chem., 102, 1699 [CrossRef] [Google Scholar]
  73. Gustafsson, M., & Frommhold, L. 2003, A&A, 400, 1161 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  74. Harding, L. B., Guadagnini, R., & Schatz, G. C. 1993, J. Phys. Chem., 97, 5472 [CrossRef] [Google Scholar]
  75. Herbst, E., & Klemperer, W. 1973, ApJ, 185, 505 [NASA ADS] [CrossRef] [Google Scholar]
  76. Hirano, S., & Yoshida, N. 2013, ApJ, 763, 52 [NASA ADS] [CrossRef] [Google Scholar]
  77. Hollenbach, D., & McKee, C. F. 1979, ApJS, 41, 555 [NASA ADS] [CrossRef] [Google Scholar]
  78. Hummer, D. G., & Storey, P. J. 1998, MNRAS, 297, 1073 [NASA ADS] [CrossRef] [Google Scholar]
  79. Huq, M. S., Doverspike, L. D., Champion, R. L., et al. 1982, J. Phys. B, 15, 951 [NASA ADS] [CrossRef] [Google Scholar]
  80. Indriolo, N., & McCall, B. J. 2012, ApJ, 745, 91 [NASA ADS] [CrossRef] [Google Scholar]
  81. Janev, R. K., Langer, W. D., Evans, K. J., et al. 1987, Elementary Processes in Hydrogen-Helium Plasmas (Berlin: Springer Verlag) [CrossRef] [Google Scholar]
  82. Jenkins, E. B. 2009, ApJ, 700, 1299 [NASA ADS] [CrossRef] [Google Scholar]
  83. Jensen, M. J., Bilodeau, R. C., Safvan, C. P., et al. 2000, ApJ, 543, 764 [NASA ADS] [CrossRef] [Google Scholar]
  84. Jones, J. D. C., Birkinshaw, K., & Twiddy, N. D. 1981, Chem. Phys. Lett., 77, 484 [NASA ADS] [CrossRef] [Google Scholar]
  85. Jorgensen, U. G., Hammer, D., Borysow, A., et al. 2000, A&A, 361, 283 [NASA ADS] [Google Scholar]
  86. Karpas, Z., Anicich, V., & Huntress, W. T. 1979, J. Chem. Phys., 70, 2877 [NASA ADS] [CrossRef] [Google Scholar]
  87. Kim, J. K., Theard, L. P., & Huntress, W. T. 1974, Int. J. Mass Spectrom. Ion Phys., 15, 223 [NASA ADS] [CrossRef] [Google Scholar]
  88. Kimura, M., Lane, N. F., Dalgarno, A., et al. 1993, ApJ, 405, 801 [NASA ADS] [CrossRef] [Google Scholar]
  89. Kingdon, J. B., & Ferland, G. J. 1996, ApJS, 106, 205 [NASA ADS] [CrossRef] [Google Scholar]
  90. Klippenstein, S. J., & Georgievskii, Y. 2010, J. Phys. Chem. A, 114, 278 [CrossRef] [Google Scholar]
  91. Krumholz, M. R., Leroy, A. K., & McKee, C. F. 2011, ApJ, 731, 25 [NASA ADS] [CrossRef] [Google Scholar]
  92. Langer, W. D. 1978, ApJ, 225, 860 [NASA ADS] [CrossRef] [Google Scholar]
  93. Lotz, W. 1967, ApJS, 14, 207 [NASA ADS] [CrossRef] [Google Scholar]
  94. Larson, Å., Le Padellec, A., Semaniak, J., et al. 1998, ApJ, 505, 459 [NASA ADS] [CrossRef] [Google Scholar]
  95. Lesaffre, P., Chieze, J.-P., Cabrit, S., et al. 2004, A&A, 427, 147 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  96. Lesaffre, P., Gerin, M., Hennebelle, P., et al. 2007, A&A, 469, 949 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  97. Le Petit, F., Nehmé, C., Le Bourlot, J., & Roueff, E. 2006, ApJS, 164, 506 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  98. Le Teuff, Y. H., Millar, T. J., & Marwick, A. J. 2000, A&AS, 146, 157 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [PubMed] [Google Scholar]
  99. Lepp, S., & Shull, J. M. 1983, ApJ, 270, 578 [NASA ADS] [CrossRef] [Google Scholar]
  100. Lim, A. J., Rawlings, J. M. C., & Williams, D. A. 1999, MNRAS, 308, 1126 [NASA ADS] [CrossRef] [Google Scholar]
  101. Linder, F., Janev, R. K., & Botero, J. 1995, Atomic and Molecular Processes in Fusion Edge Plasmas (New York: Plenum Press), 397 [CrossRef] [Google Scholar]
  102. Lipovka, A., Núñez-López, R., & Avila-Reese, V. 2005, MNRAS, 361, 850 [NASA ADS] [CrossRef] [Google Scholar]
  103. MacGregor, M., & Berry, R. S. 1973, J. Phys. B, 6, 181 [NASA ADS] [CrossRef] [Google Scholar]
  104. Maloney, P. R., Hollenbach, D. J., & Tielens, A. G. G. M. 1996, ApJ, 466, 561 [NASA ADS] [CrossRef] [Google Scholar]
  105. Massaglia, S., Mignone, A., & Bodo, G. 2005, A&A, 442, 549 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  106. Martin, P. G., Keogh, W. J., & Mandy, M. E. 1998, ApJ, 499, 793 [NASA ADS] [CrossRef] [Google Scholar]
  107. Martinez, Jr., O., Betts, N. B., Villano, S. M., et al. 2008, ApJ, 686, 1486 [NASA ADS] [CrossRef] [Google Scholar]
  108. Mauclaire, G., Derai, R., & Marx, R. 1978a, Int. J. Mass Spectrom. Ion Phys., 26, 284 [NASA ADS] [CrossRef] [Google Scholar]
  109. Mauclaire, G., Derai, R., & Marx, R. 1978b, Dyn. Mass Spectrom., 5, 139 [Google Scholar]
  110. McElroy, D., Walsh, C., Markwick, A. J., et al. 2013, A&A, 550, A36 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  111. MacLow, M. M., & Shull, J. M. 1986, ApJ, 302, 585 [NASA ADS] [CrossRef] [Google Scholar]
  112. McEwan, M. J., Scott, G. B. I., Adams, N. G., et al. 1999, ApJ, 513, 287 [NASA ADS] [CrossRef] [Google Scholar]
  113. Meijerink, R., & Spaans, M. 2005, A&A, 436, 397 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  114. Mielke, S. L., Perterson, K. A., Schwenke, D. W., et al. 2003, Phys. Rev. Lett., 91, 063201 [NASA ADS] [CrossRef] [Google Scholar]
  115. Milligan, D. B., & McEwan, M. J. 2000, Chem. Phys. Lett., 319, 482 [NASA ADS] [CrossRef] [Google Scholar]
  116. Mitchell, G. F. 1984, ApJS, 54, 81 [NASA ADS] [CrossRef] [Google Scholar]
  117. Mitchell, J. B. A. 1990, Phys. Rep., 186, 215 [NASA ADS] [CrossRef] [Google Scholar]
  118. Mitchell, G. F., & Deveau, T. J. 1983, ApJ, 266, 646 [NASA ADS] [CrossRef] [Google Scholar]
  119. Murrell, J. N., & Rodriguez, J. A. 1986, J. Mol. Struct. Theochem., 139, 267 [CrossRef] [Google Scholar]
  120. Natarajan, K., & Roth, P. 1987, Combust. Flame, 70, 267 [CrossRef] [Google Scholar]
  121. Nelson, R. P., & Langer, W. D. 1999, ApJ, 524, 923 [NASA ADS] [CrossRef] [Google Scholar]
  122. Neufeld, D. A., & Kaufman, M. J. 1993, ApJ, 418, 263 [NASA ADS] [CrossRef] [Google Scholar]
  123. Neufeld, D. A., Lepp, S., & Melnick, G. J. 1995, ApJS, 100, 132 [NASA ADS] [CrossRef] [Google Scholar]
  124. Novotny, O., Berg, M. H., Buhr, H., et al. 2010, J. Phys. Chem. A, 114, 4870 [CrossRef] [PubMed] [Google Scholar]
  125. Oldenborg, R. C., Loge, G. W., Harradine, D. M., et al. 1992, J. Phys. Chem., 96, 8426 [CrossRef] [Google Scholar]
  126. Omukai, K. 2000, ApJ, 534, 809 [NASA ADS] [CrossRef] [Google Scholar]
  127. Omukai, K., Tsuribe, T., Schneider, R., et al. 2005, ApJ, 626, 627 [Google Scholar]
  128. Oppenheimer, B. D., & Schaye, J. 2013a, MNRAS, 434, 1043 [NASA ADS] [CrossRef] [Google Scholar]
  129. Oppenheimer, B. D., & Schaye, J. 2013b, MNRAS, 434, 1063 [NASA ADS] [CrossRef] [Google Scholar]
  130. Orel, A. E. 1987, J. Chem. Phys., 87, 314 [NASA ADS] [CrossRef] [Google Scholar]
  131. Osterbrock, D. E., & Ferland, G. J. 2005, Astrophysics of Gaseous Nebulae and Active Galactic Nuclei, (University Science Books) [Google Scholar]
  132. Petuchowski, S. J., Dwek, E., Allen, Jr., J. E., et al. 1989, ApJ, 342, 406 [NASA ADS] [CrossRef] [Google Scholar]
  133. Poulaert, G., Brouillard, F., Claeys, W., et al. 1978, J. Phys. B, 11, L671 [NASA ADS] [CrossRef] [Google Scholar]
  134. Prasad, S. S., & Huntress, W. T. 1980, ApJS, 43, 1 [NASA ADS] [CrossRef] [Google Scholar]
  135. Raksit, A. B., & Warneck, P. 1980, J. Chem. Soc. Faraday Trans., 76, 1084 [CrossRef] [Google Scholar]
  136. Ramaker, D. E., & Peek, J. M. 1976, Phys. Rev. A, 13, 58 [NASA ADS] [CrossRef] [Google Scholar]
  137. Richings, A. J., Schaye, J., & Oppenheimer, B. D. 2014, MNRAS, 440, 3349 [NASA ADS] [CrossRef] [Google Scholar]
  138. Ripamonti, E., & Abel, T. 2004, MNRAS, 348, 1019 [NASA ADS] [CrossRef] [Google Scholar]
  139. Rodgers, S. D., & Millar, T. J. 1996, MNRAS, 280, 1046 [NASA ADS] [CrossRef] [Google Scholar]
  140. Rosén, S., Peverall, R., Larsson, M., et al. 1998, Phys. Rev., 57, 4462 [Google Scholar]
  141. Rosén, S., Derkatch, A., Semaniak, J., et al. 2000, Faraday Discuss., 115, 295 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  142. Shavitt, I. 1959, J. Chem. Phys., 31, 1359 [NASA ADS] [CrossRef] [Google Scholar]
  143. Savin, D. W. 2002, ApJ, 566, 599 [NASA ADS] [CrossRef] [Google Scholar]
  144. Savin, D. W., Krstic, P. S., Haiman, Z., et al. 2004, ApJ, 606, L167 [NASA ADS] [CrossRef] [Google Scholar]
  145. Schmutzler, T., & Tscharnuter, W. M. 1993, A&A, 273, 318 [NASA ADS] [Google Scholar]
  146. Schneider, I. F., Dulieu, O., Giusti-Suzor, A., et al. 1994, ApJ, 424, 983 [NASA ADS] [CrossRef] [Google Scholar]
  147. Schulz, G. J., & Asundi, R. K. 1967, Phys. Rev., 158, 25 [NASA ADS] [CrossRef] [Google Scholar]
  148. Scott, G. B. I., Fairley, D. A., Freeman, C. G., et al. 1997, J. Chem. Phys., 106, 3982 [NASA ADS] [CrossRef] [Google Scholar]
  149. Semenov, D., Hersant, F., Wakelam, V., et al. 2010, A&A, 522, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  150. Shapiro, P. R., & Kang, H. 1987, ApJ, 318, 32 [NASA ADS] [CrossRef] [Google Scholar]
  151. Shull, J. M., & van Steenberg, M. 1982, ApJS, 48, 95 [NASA ADS] [CrossRef] [Google Scholar]
  152. Sidhu, K. S., Miller, S., & Tennyson, J. 1992, A&A, 255, 453 [NASA ADS] [Google Scholar]
  153. Silk, J. 1970, Astrophys. Lett., 5, 283 [NASA ADS] [Google Scholar]
  154. Singh, P. D., Sanzovo, G. C., Borin, A. C., et al. 1999, MNRAS, 303, 235 [NASA ADS] [CrossRef] [Google Scholar]
  155. Smith, D., Spanel, P., & Mayhew, C. A. 1992, Int. J. Mass spectrom. Ion Proc., 117, 457 [NASA ADS] [CrossRef] [Google Scholar]
  156. Smith, D., & Adams, N. G. 1977a, Int. J. Mass Spectrom. Ion Phys., 23, 123 [NASA ADS] [CrossRef] [Google Scholar]
  157. Smith, D., & Adams, N. G. 1977b, Chem. Phys. Lett., 47, 383 [NASA ADS] [CrossRef] [Google Scholar]
  158. Smith, D., Adams, N. G., & Miller, T. M. 1978, J. Chem. Phys., 69, 308 [NASA ADS] [CrossRef] [Google Scholar]
  159. Smith, M. A., Schlemmer, S., von Richthofen, J., et al. 2002, ApJ, 578, L87 [NASA ADS] [CrossRef] [Google Scholar]
  160. Smith, I. W. M., Herbst, E., & Chang, Q. 2004, MNRAS, 350, 323 [NASA ADS] [CrossRef] [Google Scholar]
  161. Spitzer, Jr., L. 1978, Physical Processes in the Interstellar Medium (New York: Wiley) [Google Scholar]
  162. Stancil, P. C., & Dalgarno, A. 1998, Faraday Discuss., 109, 61 [NASA ADS] [CrossRef] [Google Scholar]
  163. Stancil, P. C., Schultz, D. R., Raković, M., et al. 2015, Charge Transfer Database http://www.physast.uga.edu/ugacxdb/ [Google Scholar]
  164. Stone, J. M., & Norman, M. L. 1992, ApJS, 80, 791 [Google Scholar]
  165. Sutherland, R. S., & Dopita, M. A. 1993, ApJS, 88, 253 [NASA ADS] [CrossRef] [Google Scholar]
  166. Takagi, H., Kosugi, N., & Le Dourneuf, M. 1991, J. Phys. B, 24, 711 [NASA ADS] [CrossRef] [Google Scholar]
  167. Teşileanu, O., Mignone, A., & Massaglia, S. 2008, A&A, 488, 429 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  168. Teyssier, R. 2002, A&A, 385, 337 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  169. Trevisan, C. S., & Tennyson, J. 2002, Plasma Phys. Control. Fusion, 44, 1263 [NASA ADS] [CrossRef] [Google Scholar]
  170. Tsang, W., & Hampson, R. F. 1986, J. Phys. Chem. Ref. Data, 15, 1087 [NASA ADS] [CrossRef] [Google Scholar]
  171. Vasiliev, E. O. 2013, MNRAS, 431, 638 [NASA ADS] [CrossRef] [Google Scholar]
  172. Verner, D. A., & Ferland, G. J. 1996, ApJS, 103, 467 [NASA ADS] [CrossRef] [Google Scholar]
  173. Viggiano, A. A., Howorka, F., Albritton, D. L., et al. 1980, ApJ, 236, 492 [NASA ADS] [CrossRef] [Google Scholar]
  174. Voronov, G. S. 1997, At. Data Nucl. Data Tables, 65, 1 [Google Scholar]
  175. Wagner-Redeker, W., Kemper, P. R., Jarrold, M. F., et al. 1985, J. Chem. Phys., 83, 1121 [NASA ADS] [CrossRef] [Google Scholar]
  176. Walkauskas, L. P., & Kaufman, F. 1975, Symp. Int. Combust. Proc., 15, 691 [CrossRef] [Google Scholar]
  177. Wakelam, V., Selsis, F., Herbst, E., et al. 2005, A&A, 444, 883 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  178. Wakelam, V., Herbst, E., Loison, J.-C., et al. 2012, ApJS, 199, 21 [NASA ADS] [CrossRef] [Google Scholar]
  179. Warnatz, J. 1984, in Combustion Chemistry, ed. W. C. Gardiner, Jr. (NY: Springer-Verlag), 197 [CrossRef] [Google Scholar]
  180. Watson, W. D., Anicich, V. G., & Huntress, W. T. 1976, ApJ, 205, L165 [NASA ADS] [CrossRef] [Google Scholar]
  181. Weingartner, J. C., & Draine, B. T. 2001a, ApJ, 548, 296 [NASA ADS] [CrossRef] [Google Scholar]
  182. Weingartner, J. C., & Draine, B. T. 2001b, ApJ, 563, 842 [NASA ADS] [CrossRef] [Google Scholar]
  183. Weingartner, J. C., & Draine, B. T. 2001c, ApJS, 134, 263 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  184. Williams, J. P., Bergin, E. A., & Caselli, P. 1998, ApJ, 503, 689 [NASA ADS] [CrossRef] [Google Scholar]
  185. Wishart, A. W. 1979, MNRAS, 187, 59P [NASA ADS] [CrossRef] [Google Scholar]
  186. Woodall, J., Agúndez, M., Markwick-Kemper, A. L., et al. 2007, A&A, 466, 1197 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  187. Yan, M. 1997, PhD Thesis, Center for Astrophysics, Harvard University [Google Scholar]
  188. Ziegler, U. 2008, CPC, 179, 227 [NASA ADS] [CrossRef] [Google Scholar]
  189. Ziegler, U. 2011, JCP, 230, 1035 [Google Scholar]
  190. Ziegler, U. 2012, SIAM J. Sci. Comput., 34, C102 [CrossRef] [Google Scholar]
  191. Ziegler, U. 2016, A&A, 586, A82 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  192. Zygelman, B., Dalgarno, A., Kimura, M., et al. 1989, Phys. Rev. A, 40, 2340 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: The chemical network

Table A.1.

NC​2M – gas-phase reactions.

Table A.2.

NC​2M – CR-induced reactions.

Table A.3.

NC​2M – dust-catalytic reactions.

All Tables

Table 1.

Recommended numbers of atomic energy levels for metal cooling computations as adopted in the NC​2M (z: ion charge).

In the text
Table 2.

Sources of .

In the text
Table 3.

Fitting coefficients for expression 66.

In the text
Table A.1.

NC​2M – gas-phase reactions.

In the text
Table A.2.

NC​2M – CR-induced reactions.

In the text
Table A.3.

NC​2M – dust-catalytic reactions.

In the text

All Figures

thumbnail Fig. 1.

Influence of parameter L on the cooling rate (top panels) and on the per-ion cooling coefficient (bottom panels) for OIII and FeVI. The dashed line corresponds to the highest possible L with Chianti. The results were obtained for an equilibrium model with nH = 106 m−3.
In the text
thumbnail Fig. 2.

Ratio of spectral line intensities obtained with the NC​2M (solid lines) and theoretical estimates (dashed lines, Eqs. (52)–(54)) for nH = 104 m−3.
In the text
thumbnail Fig. 3.

CIE cooling fitting curves (lines) for the different collision pairs according to Eq. (66), and data points (symbols) derived from the sources in Table 2.
In the text
thumbnail Fig. 4.

Fractions of H2 (first row), CO (second row), H2O (third row), O2 (fourth row), and OH (fifth row) for different densities (from left to right). Equilibrium solutions for the full network (solid), dust-free network (dashed), CR-free network (dash-dotted), and pure gas-phase network (dotted) are shown.
In the text
thumbnail Fig. 5.

Ionization structure of C for different densities. Equilibrium solution for the full network (solid), dust-free network (dashed), and CR-free network (dash-dotted).
In the text
thumbnail Fig. 6.

Analog to Fig. 5 for O.
In the text
thumbnail Fig. 7.

Analog to Fig. 5 for Si.
In the text
thumbnail Fig. 8.

Analog to Fig. 5 for Fe.
In the text
thumbnail Fig. 9.

Molecular cooling in chemical equilibrium for the full network.
In the text
thumbnail Fig. 10.

Metal cooling in chemical equilibrium for the full network. Top panel: contribution of individual oxygen ionization stages. Bottom panel: contribution of different metals. The dashed line represents the total metal cooling when the number of energy levels is restricted to 5 (L = 5) in all ion computations.
In the text
thumbnail Fig. 11.

Total cooling function and its contributions from different cooling or heating agents in chemical equilibrium for the full network. Dash-dotted lines indicate heating.
In the text
thumbnail Fig. 12.

Impact of different gas densities and the absence of dust or cosmic rays on the net cooling function in chemical equilibrium. Dash-dotted lines indicate net heating, Λ* <  0.
In the text
thumbnail Fig. 13.

Ionization fractions of HeIII, CIV, OVII, and FeV obtained for the full network in equilibrium (dashed lines) and non-equilibrium (solid lines) in the 108 K cooling gas experiment.
In the text
thumbnail Fig. 14.

Cooling curve in equilibrium (dashed lines) and in non-equilibrium (solid lines) for the full network (black), dust-free network (yellow), and CR-free network (blue). For T ≳ 100 K the solid lines representing the different network cases nearly coincide. Crosses: Non-equilibrium cooling curve for a network without the metals N, Mg and Ne.
In the text
thumbnail Fig. 15.

Electron fraction (top panel) and molecular hydrogen fraction (bottom panel) as a function of temperature in equilibrium (dashed lines) and in non-equilibrium (solid lines) for the full network (black), dust-free network (yellow), and CR-free network (blue).
In the text
thumbnail Fig. 16.

Time dependence of the temperature (top panel) and the net cooling rate (panel 2), electron fraction (panel 3), and molecular hydrogen fraction (bottom panel) as a function of temperature for case A (black solid lines), case B (red lines), and for the equilibrium solution (dashed lines).
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.