© ESO, 2016

1. Introduction

Astrochemistry deals with chemical reactions in astrophysical gases and, taking a broader view, embraces ionizing processes, dust-molecule interactions, and photo-chemical effects. Its significance in certain environments results from the fact that thermal and magnetic properties are determined by species abundances and the ionization state which, both, are regulated by chemical processes. Through diverse cooling and heating mechanisms entering the energy budget, and by controlling the coupling between magnetic field and matter, chemistry thus also influences the dynamical evolution of the system. The interstellar medium with its different prevailing phases of diffuse gas, molecular clouds, and dense gas cores is an environment where rich chemical processes take place. In the densest gas aggregates chemistry potentially manipulates fragmentation and the following gravitational collapse of pre-stellar cores to new stars. Further importance of interstellar chemistry comes from the role certain molecules play as tracers in observations such as the CO molecule for molecular gas detection or the NH₃ molecule as an indicator of the gas temperature.

In order to bring magneto-gasdynamical simulations and observations closer, astrochemistry must find its way into numerical models. However, embedding chemistry in time-dependent, multidimensional computations is a very challenging task. A major difficulty arises from a possible large disparity between dynamical and reactive timescales in the flow. Since the overall complexity of chemo-magneto-gasdynamical simulations requires split integration techniques a sufficiently good coupling between dynamics and chemistry is necessary for accuracy and stability reasons. In practice, this forces a reduction in the size of the time-step, which can be used in advancing the solution since chemical kinetics may involve smaller timescales than the dynamical processes. In addition, the computational cost per numerical cell for solving the chemical rate equations quickly dominates the magneto-gasdynamical part for larger chemical networks. In consequence, this often prohibits the use of complex chemical networks in high-resolution simulations. What makes chemical kinetics so expensive is the need for sophisticated ordinary differential equation (ODE) solvers that can account for the stiffness inherent to the system. This comes about because individual reactions can proceed with very different speed. Cheap explicit schemes are ruled out because of their strict stability requirements which would severely limit the time-step according to the fastest reaction speed. A further complication arises because chemical kinetics and their associated thermal processes are often closely connected with regard to their time evolution. A treatment of both effects within a common numerical framework seems therefore reasonable.

The problem of incorporating chemistry and related thermal processes in astrophysical simulations has attracted substantial attention in the past. Early attempts to combine chemistry with time-dependent gasdynamics were undertaken by Lim et al. (1999) in 2D and by Abel et al. (2002) in 3D using the ENZO code (Bryan et al. 2014). In a work by Glover & Mac Low (2007), who adopt the ZEUS code (Stone & Norman 1992), a rudimentary hydrogen network was evolved consistently with the equations of magneto-gasdynamics in order to study molecular cloud formation. Nonequilibrium models describing the ionization of an atomic multispecies gas with cooling in optically thin plasma have also been implemented into magneto-gasdynamical codes like ASTROBEAR (Poludnenko et al. 2005) or PLUTO (Teşileanu et al. 2008). Moreover, there are powerful 1D time-dependent codes that integrate magneto-gasdynamics alongside more complex chemical networks used, e.g., in high-resolution studies of magnetic shocks (Lesaffre et al. 2004) or diffusion fronts (Lesaffre et al. 2007). Stand-alone astrochemistry packages (not coupled to dynamics) have been designed for special purposes like the ALCHEMIC code (Semenov et al. 2010) for the chemistry in proto-planetary discs, the NAHOON code (Wakelam et al. 2012) for the study of the chemical evolution in molecular clouds, or the MEUDONPDR code (Le Petit et al. 2006) that couples chemistry with radiation transfer suitable for radiation-dominated regions. Recently, a publicly available astrochemistry package called KROME has been developed (Grassi et al. 2014) based on the DLSODES solver from the Hindmarsh ODE tool (Hindmarsh 1983). This code aims at a more general-purpose chemistry framework for the astrophysics community, and claims to be easily includable in any fluid dynamics simulation code which provides a multispecies functionality. KROME comes with a variety of predefined chemical networks and is prepared for a treatment of chemistry-related thermal processes, photochemistry and dust evolution.

This paper describes the implementation of a chemical reaction network solver into the astrophysics code NIRVANA¹ (Ziegler 2008, 2011, 2012). The chemistry module is constructed from scratch and adapted to the given code infrastructure. This proprietary coding effort allows full control over all aspects of implementation, and this was the main reason for not incorporating external software like KROME.

2. Chemical reaction network

A set of N_r chemical reactions of the form $\begin{array}{c}\sum \\ \mathit{s}{\mathit{X}}_{\mathit{s}}\hspace{0.17em}\mathrm{\in }\hspace{0.17em}\mathit{R}\mathrm{\left(}\mathit{r}\mathrm{\right)}\end{array}{\mathit{\alpha }}_{\mathit{s,r}}{\mathit{X}}_{\mathit{s}}\stackrel{{\mathit{k}}_{\mathit{r}}}{\mathrm{-}0\mathrm{\to }}\begin{array}{c}\sum \\ \mathit{s}{\mathit{X}}_{\mathit{s}}\hspace{0.17em}\mathrm{\in }\hspace{0.17em}\mathit{P}\mathrm{\left(}\mathit{r}\mathrm{\right)}\end{array}{\mathit{\beta }}_{\mathit{s,r}}{\mathit{X}}_{\mathit{s}}\mathit{,} \mathit{r}\mathrm{=}\mathrm{1}\mathit{,N}r$(1)between species $\mathrm{\{}{\mathit{X}}_{\mathit{s}}{\mathrm{\}}}_{\mathit{s}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{1}}^{{\mathit{N}}_{\mathit{s}}\mathrm{+}\mathrm{1}}$ is considered where, in each reaction r, a number of reactants X_s ∈ R(r) are transformed into a number of products X_s ∈ P(r) with R(r) and P(r) being subsets of { X_s }. The reaction speed is determined by the kinetic rate coefficient k_r(T), which is assumed to be a function of temperature T. The α’s and β’s are stoichiometric coefficients. Chemical species can be atoms, molecules, positively or negatively charged ions, and free electrons. The total number of species is N_s + 1 with index s = N_s + 1 reserved for electrons. A species X_s is associated with a number density n_s (expressed in m^-3), which is a function of space and time, a molecular weight μ_s and a charge q_s. At the current stage of development reactions including dust particles and photochemical processes are not considered. In principle, both effects could be treated within the framework presented here but would lead to a significant increase in complexity. It is naturally assumed that changes in the electron density are exclusively caused by changes in the ionic state of other species through ionization/recombination processes, electron attachment/detachment or charge transfer mechanisms. In this case, the electron density n_e (=n_{Ns + 1}) need not be treated as an independent variable but is given by ${\mathit{n}}_{\mathrm{e}}\mathrm{=}\sum _{\mathit{s}}{\mathit{q}}_{\mathit{s}}\mathit{n}s$(2)expressing total charge conservation². Here and in the following, an unspecified summation index s is assumed to run from 1 to N_s (i.e., excluding free electrons).

The chemical network module of NIRVANA is designed to handle an arbitrary number of species with an arbitrary number of reactions involved allowing an arbitrary number of reactants and products. Technically, individual reactions are easily specified as a single line in a text file obeying a generic input format of the form ID S R [& S R [& ...]] > S P [& S P [& ...]], where ID (a running number) uniquely identifies the reaction, S are stoichiometric coefficients (the α’s and β’s), R stands for reactant, and P stands for product. For example, for the reaction H₂+ e⁻→ 2 H + e⁻ being number twelve in a network the format would be 12 1 H2 & 1 e- > 2 H & 1 e-. By parsing the “reactions file” the code checks for the mass and charge balance in each reaction, and automatically derives and saves the structure of the reaction term (4) below which enters the ODE solver for chemical kinetics. A by-hand definition of the ODE system is redundant. Since reaction rates k_r can be complicated functions of temperature not necessarily fitting a prescribed generic expression, those rates are specified separately in order to permit the highest flexibility.

As a basis in this work, a simple network for a hydrogen–helium plasma was constructed to be used in some of the presented test problems. It consists of the 19 reactions cataloged in Table A.1 for the nine species H, H⁺, He, He⁺, He²⁺, H⁻, H₂, H$\begin{array}{c}\mathrm{+}\\ \mathrm{2}\end{array}$, and free electrons. The network includes ionization/recombination of H, He and a channel for H₂ formation via the H⁻, H$\begin{array}{c}\mathrm{+}\\ \mathrm{2}\end{array}$ ions.

2.1. Advection-reaction system

In gasdynamical simulations the time evolution of species number densities is described by a coupled system of advection–reaction equations, ${\mathit{\partial }}_{\mathit{t}}{\mathit{n}}_{\mathit{s}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{\left(}{\mathit{n}}_{\mathit{s}}{v}\mathrm{\right)}\mathrm{=}{\mathit{R}}_{\mathit{s}}\hspace{0.17em} \mathit{s}\mathrm{=}\mathrm{1}\mathit{,}{\mathit{N}}_{\mathit{s}}\mathit{,}$(3)with reaction term ${\mathit{R}}_{\mathit{s}}\mathrm{=}\sum _{\mathit{r}}\mathrm{(}{\mathit{\beta }}_{\mathit{s,r}}\mathrm{-}{\mathit{\alpha }}_{\mathit{s,r}}\mathrm{)}{\mathit{k}}_{\mathit{r}}\mathrm{\left(}\mathit{T}\mathrm{\right)}\begin{array}{c}\mathrm{􏽙}\\ \mathit{l}{\mathit{X}}_{\mathit{l}}\hspace{0.17em}\mathrm{\in }\hspace{0.17em}\mathit{R}\mathrm{\left(}\mathit{r}\mathrm{\right)}\end{array}{\mathit{n}}_{\mathit{l}}^{{\mathit{\alpha }}_{\mathit{l,r}}}$(4)subject to constraint (2). The product on the right-hand side running over reactants in a reaction makes the system nonlinear. Equation (3) has to be supplemented, of course, by other fluid equations describing the evolution of the flow speed v and pressure p (respective energy density). Maxwell’s equations add in magneto-gasdynamical problems involving the evolution of a magnetic field B. Pressure, temperature and chemical composition of the gas are connected via the ideal equation of state

$\mathit{p}\mathrm{=}{\mathit{k}}_{\mathrm{B}}\mathit{T}{\mathit{n}}_{\mathrm{tot}}\mathrm{=}{\mathit{k}}_{\mathrm{B}}\mathit{T}\frac{\mathit{\varrho }}{{\mathit{m}}_{\mathit{u}}\overline{\mathit{\mu }}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{Y}\mathrm{)}$with n_tot = ∑ _sn_s + n_e the total number density, k_B the Boltzmann constant and m_u the atomic mass unit. Taking constraint (2) into account, the gas density ϱ is linked to species densities by the relation $\mathit{\varrho }\mathit{/}{\mathit{m}}_{\mathit{u}}\mathrm{=}\sum _{\mathit{s}}{\mathit{\mu }}_{\mathit{s}}{\mathit{n}}_{\mathit{s}}\mathrm{+}{\mathit{\mu }}_{\mathrm{e}}{\mathit{n}}_{\mathrm{e}}\mathrm{=}\sum _{\mathit{s}}\begin{array}{c}\mathrm{˜}\\ {\mathit{\mu }}_{\mathit{s}}\end{array}\mathit{n}s$(5)with $\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{\mu }}_{\mathit{s}}\end{array}\mathrm{=}{\mathit{\mu }}_{\mathit{s}}\mathrm{+}{\mathit{q}}_{\mathit{s}}{\mathit{\mu }}_{\mathrm{e}}$ and the electron weight μ_e = 5.49 × 10^-4. The mean molecular weight $\overline{\mathit{\mu }}$ and electron number fraction Y is defined by

$\overline{\mathit{\mu }}\mathrm{=}\frac{{\sum }_{\mathit{s}}\begin{array}{c}\mathrm{˜}\\ {\mathit{\mu }}_{\mathit{s}}\end{array}{\mathit{n}}_{\mathit{s}}}{{\sum }_{\mathit{s}}{\mathit{n}}_{\mathit{s}}}\mathit{,} \mathit{Y}\mathrm{=}\frac{{\mathit{n}}_{\mathrm{e}}}{{\sum }_{\mathit{s}}{\mathit{n}}_{\mathit{s}}}\mathrm{·}$Usually, in gasdynamical simulations a continuity equation for ϱ expressing mass conservation is evolved in time: ${\mathit{\partial }}_{\mathit{t}}\mathit{\varrho }\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{\left(}\mathit{\varrho }{v}\mathrm{\right)}\mathrm{=}\mathrm{0.}$(6)Multiplying Eq. (3) by $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{\mu }}_{\mathit{s}}\end{array}$ and summing over species results in $\begin{array}{ccc}\sum _{\mathit{s}}\begin{array}{c}\mathrm{˜}\\ {\mathit{\mu }}_{\mathit{s}}\end{array}{\mathit{\partial }}_{\mathit{t}}{\mathit{n}}_{\mathit{s}}\mathrm{+}\sum _{\mathit{s}}\begin{array}{c}\mathrm{˜}\\ {\mathit{\mu }}_{\mathit{s}}\end{array}\mathrm{\nabla }\mathrm{·}\mathrm{\left(}{\mathit{n}}_{\mathit{s}}{v}\mathrm{\right)}& \mathrm{=}& \sum _{\mathit{s}}\begin{array}{c}\mathrm{˜}\\ {\mathit{\mu }}_{\mathit{s}}\end{array}{\mathit{R}}_{\mathit{s}}\\ & \mathrm{=}& \sum _{\mathit{r}}{\mathit{k}}_{\mathit{r}}\begin{array}{c}\mathrm{􏽙}\\ \mathit{l}\mathit{l}\hspace{0.17em}\mathrm{\in }\hspace{0.17em}\mathit{R}\mathrm{\left(}\mathit{r}\mathrm{\right)}\end{array}{\mathit{n}}_{\mathit{l}}^{{\mathit{\alpha }}_{\mathit{l,r}}}\sum _{\mathit{s}}\begin{array}{c}\mathrm{˜}\\ {\mathit{\mu }}_{\mathit{s}}\end{array}\mathrm{(}{\mathit{\beta }}_{\mathit{s,r}}\mathrm{-}{\mathit{\alpha }}_{\mathit{s,r}}\mathrm{)}\\ & \mathrm{=}& \frac{\mathrm{1}}{{\mathit{m}}_{\mathit{u}}}\left[{\mathit{\partial }}_{\mathit{t}}\mathit{\varrho }\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{\left(}\mathit{\varrho }{v}\mathrm{\right)}\right]\\ & \mathrm{=}& \mathrm{0.}\end{array}$(7)Hence, consistency with the continuity equation requires $\sum _{\mathit{s}}\begin{array}{c}\mathrm{˜}\\ {\mathit{\mu }}_{\mathit{s}}\end{array}{\mathit{\alpha }}_{\mathit{s,r}}\mathrm{=}\sum _{\mathit{s}}\begin{array}{c}\mathrm{˜}\\ {\mathit{\mu }}_{\mathit{s}}\end{array}{\mathit{\beta }}_{\mathit{s,r}} \mathrm{\forall }\mathit{r}$(8)expressing mass conservation in each reaction.

2.2. Thermal processes

Chemical reactions in astrophysical gases are often accompanied by special thermal processes. For instance, collisional ionization of hydrogen leads to a cooling of the gas because the energy needed to pass the hydrogen atom from its neutral to ionized stage is extracted from the thermal energy. Radiative recombination of ionized hydrogen to its neutral stage also leads to cooling of the gas because a free electron is captured, which removes part of its kinetic energy. The formation of H₂, on the other hand, is often an exothermic process that heats the gas. Macroscopically, cooling and heating processes are taken into account by a source term S_e in the thermal energy equation (and total energy equation). Ignoring dissipative processes ${\mathit{\partial }}_{\mathit{t}}\mathit{e}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{\left(}\mathit{e}{v}\mathrm{\right)}\mathrm{=}\mathrm{-}\mathit{p}\mathrm{\nabla }\mathrm{·}{v}\mathrm{+}\mathit{S}e$(9)holds for the thermal energy density e = p/ (γ−1) with adiabatic index γ assumed to be a constant value³. Since the full system of chemo-magneto-gasdynamical equations is solved in an operator-split fashion, from Eqs. (3) and (9) the thermo-reactive part in the solution procedure (ϱ, v, and B are assumed to be constant during this update) is governed by

$\begin{array}{ccc}& & \dot{{\mathit{n}}_{\mathit{s}}}\mathrm{=}{\mathit{R}}_{\mathit{s}}\hspace{0.17em} \mathit{s}\mathrm{=}\mathrm{1}\mathit{,}{\mathit{N}}_{\mathit{s}}\\ & & \frac{\mathrm{1}}{\mathit{\gamma }\mathrm{-}\mathrm{1}}\mathit{ṗ}\mathrm{=}{\mathit{S}}_{\mathrm{e}}\mathit{.}\end{array}$Using the equation of state, the expression for the total number density and Eq. (10) one has

$\mathit{ṗ}\mathrm{=}{\mathit{k}}_{\mathrm{B}}\mathit{ṅ}\mathrm{tot}\mathit{T}\mathrm{+}{\mathit{k}}_{\mathrm{B}}{\mathit{n}}_{\mathrm{tot}}\mathit{Ṫ}\mathrm{=}{\mathit{k}}_{\mathrm{B}}\mathit{T}\sum _{\mathit{s}}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{q}}_{\mathit{s}}\mathrm{)}{\mathit{R}}_{\mathit{s}}\mathrm{+}{\mathit{k}}_{\mathrm{B}}{\mathit{n}}_{\mathrm{tot}}\mathit{Ṫ}\hspace{0.17em}\mathit{.}$It follows for the evolution of temperature $\mathit{Ṫ}\mathrm{=}\frac{\mathit{\gamma }\mathrm{-}\mathrm{1}}{{\mathit{k}}_{\mathrm{B}}{\mathit{n}}_{\mathrm{tot}}}{\mathit{S}}_{\mathrm{e}}\mathrm{-}\frac{\mathit{T}}{{\mathit{n}}_{\mathrm{tot}}}\sum _{\mathit{s}}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{q}}_{\mathit{s}}\mathrm{)}{\mathit{R}}_{\mathit{s}}\mathrm{=}\mathrm{:}{\mathit{S}}_{\mathit{T}}\mathrm{\left(}{n}\mathit{,T}\mathrm{\right)}$(12)which is used instead of (11) in the operator-split formalism presented below.

The source term S_e, which generally depends on the species densities and temperature, is decomposed here into two parts according to ${\mathit{S}}_{\mathrm{e}}\mathrm{\left(}{n}\mathit{,T}\mathrm{\right)}\mathrm{=}\sum _{\mathit{r}}{\mathit{l}}_{\mathit{r}}\mathrm{\left(}\mathit{T}\mathrm{\right)}\underset{\mathit{l}\hspace{0.17em}\mathrm{\in }\hspace{0.17em}\mathit{R}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\mathrm{􏽙}}{\mathit{n}}_{\mathit{l}}^{{\mathit{\alpha }}_{\mathit{l,r}}}\mathrm{+}{\mathit{L}}_{\mathrm{e}}\mathrm{\left(}{n}\mathit{,T}\mathrm{\right)}\mathit{.}$(13)The first (canonical) term on the right-hand side has the same structural form as the source term R_s in the advection-reaction equation for n_s with a thermal rate coefficient l_r(T) instead of kinetic rate k_r(T) associated with a reaction r. Many thermal processes fit this canonical form with l_r ∝ k_r in some cases. In order to offer a high degree of flexibility other processes not describable within the first term of the energy source are subsumed in the second term L_e(n,T). Whenever possible, however, thermal terms should be treated within the canonical source term simply for reasons of computational efficiency: derivatives of products ${{}^{\mathrm{􏽑}}}_{\mathit{l}\mathrm{\in }\mathit{R}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\mathit{n}}_{\mathit{l}}^{{\mathit{\alpha }}_{\mathit{l,r}}}$ with respect to n_s (needed to compute the Jacobian matrix) must already be evaluated for the source terms R_s and so can be reused.

A summary of the thermal processes that have been considered in combination with the chemical network in Table A.1 is presented in Table B.1. A more detailed description of individual thermal processes can also be found in Appendix B, which includes examples for L_e-terms.

3. Numerical implementation

The following presentation concentrates on implementation details related to the introduction of reactive species. Other solution parts of the full system of chemo-magneto-gasdynamical equations are not described here. In this coverage numerical integration of the coupled system of advection-reaction Eq. (3) and temperature Eq. (12) must be carried out with respect to co-existing different timescales inherent to flow advection on the one hand and chemical kinetics on the other hand. Moreover, chemical reactions and related thermal processes may, on their own, proceed with disparate speeds introducing stiffness to the system and requiring robust and efficient solvers. Overall, the system is divided into an advection part described by the PDE ${\mathit{\partial }}_{\mathit{t}}{\mathit{n}}_{\mathit{s}}\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{\left(}{\mathit{n}}_{\mathit{s}}{v}\mathrm{\right)}\mathrm{=}\mathrm{0}\mathit{,} \mathit{s}\mathrm{=}\mathrm{1}\mathit{,N}s$(14)which is solved simultaneously with the advection of density, momentum equation, and energy equation with S_e = 0, and a source step given by the ODE system $\frac{\mathrm{d}}{\mathrm{d}\mathit{t}}{u}\mathrm{=}{S}\mathrm{\left(}{u}\mathrm{\right)}$(15)where ${u}\mathrm{=}\mathrm{\left(}\mathrm{\right\{}{\mathit{n}}_{\mathit{s}}{\mathrm{\}}}_{\mathit{s}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathit{s}}}\mathit{,T}\mathrm{)}\mathit{,}{S}\mathrm{=}\mathrm{(}\mathrm{\left\{}{\mathit{R}}_{\mathit{s}}{\mathrm{\}}}_{\mathit{s}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathit{s}}}\mathit{,}{\mathit{S}}_{\mathit{T}}\mathrm{\right)}$. In order to advance the solution one time-step Δtⁿ from tⁿ to t^{n + 1} = tⁿ + Δtⁿ a second-order accurate Strang-type operator splitting is used. Denoting the advection and source step operators as and the update formally reads

$\begin{array}{ccc}& & {{u}}^{\mathrm{n}\mathrm{+}\mathrm{1}\mathit{/}\mathrm{2}\mathit{,}\mathit{\star }}\mathrm{=}{\mathrm{𝒮}}_{\mathrm{\Delta }\mathit{t}\mathit{/}\mathrm{2}}\mathrm{\left(}{{u}}^{\mathrm{n}}\mathrm{\right)}\\ & & {\mathit{n}}_{\mathit{s}}^{\mathrm{n}\mathrm{+}\mathrm{1}\mathit{,}\mathit{\star }}\mathrm{=}{\mathrm{𝒜}}_{\mathrm{\Delta }\mathit{t}}\mathrm{\left(}{\mathit{n}}_{\mathit{s}}^{\mathrm{n}\mathrm{+}\mathrm{1}\mathit{/}\mathrm{2}\mathit{,}\mathit{\star }}\mathrm{\right)}\mathit{,} \mathit{s}\mathrm{=}\mathrm{1}\mathit{,}{\mathit{N}}_{\mathit{s}}\\ & & {{u}}^{\mathrm{n}\mathrm{+}\mathrm{1}}\mathrm{=}{\mathrm{𝒮}}_{\mathrm{\Delta }\mathit{t}\mathit{/}\mathrm{2}}\mathrm{\left(}{{u}}^{\mathrm{n}\mathrm{+}\mathrm{1}\mathit{,}\mathit{\star }}\mathrm{\right)}\mathit{.}\end{array}$The actual size of time-step is given by

$\mathrm{\Delta }{\mathit{t}}^{\mathrm{n}}\mathrm{=}\begin{array}{c}\min \end{array}\mathrm{(}\mathrm{\Delta }{\mathit{t}}_{\mathrm{dyn}}^{\mathrm{n}}\mathit{,}\mathrm{\Delta }{{\mathit{t}}_{\mathrm{reac}}^{\mathrm{n}}}^{\mathrm{)}}$where $\mathrm{\Delta }{\mathit{t}}_{\mathrm{dyn}}^{\mathrm{n}}$ is the advection/dynamical time-step dictated by the usual CFL condition for numerical stability of the discretized gasdynamical equations and $\mathrm{\Delta }{\mathit{t}}_{\mathrm{reac}}^{\mathrm{n}}$ is a reaction time-step defined in more detail later in this section. A proper choice of the latter serves to control the coupling between advection and reaction. Generally, adopting $\mathrm{\Delta }{\mathit{t}}_{\mathrm{dyn}}^{\mathrm{n}}$ as sole time-step may be insufficient because large splitting errors can lead to unphysical behavior or even destabilize numerics in situations where $\mathrm{\Delta }{\mathit{t}}_{\mathrm{dyn}}^{\mathrm{n}}\mathrm{\gg }\mathrm{\Delta }{\mathit{t}}_{\mathrm{reac}}^{\mathrm{n}}$.

The advection of individual species is performed within the Finite-Volume (FV) framework of the NIRVANA code (Ziegler 2011). The numerical procedure is analogous to the advection of the gas density ϱ. In the FV scheme the update of volume-averaged number densities requires computing fluxes ${\mathit{F}}^{\mathrm{\left(}{\mathit{n}}_{\mathit{s}}\mathrm{\right)}}\mathrm{\left(}{\mathit{n}}_{\mathit{s}}^{\mathit{\star }}\mathrm{\right)}$ at each grid cell surface, which depend on a properly upwinded state ${\mathit{n}}_{\mathit{s}}^{\mathit{\star }}\mathrm{=}{\mathit{n}}_{\mathit{s}}^{\mathit{\star }}\mathrm{\left(}{\mathit{n}}_{\mathit{s}}^{\mathit{L}}\mathit{,}{\mathit{n}}_{\mathit{s}}^{\mathit{R}}\mathrm{\right)}$. This is derived in accordance with other fluid variables from the approximate solution of a Riemann problem with left-sided/right-sided states ${\mathit{n}}_{\mathit{s}}^{\mathit{L,R}}$ at grid cell surfaces. These are obtained by a piecewise linear reconstruction procedure based on volume-averaged values in adjacent cells with slopes subject to the TVD limiter of van Leer (1977). However, since slope limiting is a nonlinear process and species profiles are reconstructed independently from each other, the consistency condition

${\mathit{\varrho }}^{\mathit{L,R}}\mathrm{=}\begin{array}{c}\sum \end{array}{\mathit{\varrho }}_{\mathit{s}}^{\mathit{L,R}}\mathrm{=}{\mathit{m}}_{\mathit{u}}\begin{array}{c}\sum \end{array}\begin{array}{c}\mathrm{˜}\\ {\mathit{\mu }}_{\mathit{s}}\end{array}{\mathit{n}}_{\mathit{s}}^{\mathit{L,R}}$is per se not fulfilled. If left uncorrelated, deviations between the total density and sum of partial densities ϱ_s can accumulate with time and may produce unphysical solutions. Hence, some modification in the advection algorithm is advisable in order to meet the conservation principle⁴. The simplest way to deal with the problem is to overwrite ϱ^L,R and to replace the corresponding total mass flux by the sum of partial mass fluxes. Although $\mathrm{|}{\mathit{\varrho }}^{\mathit{L,R}}\mathrm{-}\sum {\mathit{\varrho }}_{\mathit{s}}^{\mathit{L,R}}\mathrm{|}$ may be small in size, with this solution there is the danger of violating monotonicity and relaxing TVD constraints in the evolution of ϱ. This, in turn, may result in overshooting effects near local extrema or at discontinuities. An alternative, which is applied here, is profile flattening by reducing species density slopes, σ_s, while leaving the evolution of total density almost untouched. Slope modification should be carried out in a clever way to diminish overshooting, to maintain TVD properties where possible, and to retain second-order accuracy. Concretely, the individual species slopes in each spatial direction and in each numerical cell are modified according to

$\begin{array}{c}\mathrm{˜}\\ {\mathit{\sigma }}_{\mathit{s}}\end{array}\mathrm{=}{\mathit{\chi }}_{\mathit{s}}{\mathit{\sigma }}_{\mathit{s}}$where

${\mathit{\chi }}_{\mathit{s}}\mathrm{=}\mathrm{1}\mathrm{-}\begin{array}{c}\max \end{array}\left(\mathrm{0}\mathit{,}\begin{array}{c}\min \end{array}\left(\mathrm{1}\mathit{,}\frac{\mathrm{(}{\mathit{\sigma }}_{\sum }\mathrm{-}{\mathit{\sigma }}_{\mathit{\varrho }}\mathrm{)}\mathrm{|}\mathrm{sgn}\hspace{0.17em}{\mathit{\sigma }}_{\mathit{s}}\mathrm{+}\mathrm{sgn}\hspace{0.17em}{\mathit{\sigma }}_{\sum }\mathrm{|}}{\sum {\mathit{\sigma }}_{\mathit{s}}\mathrm{|}\mathrm{sgn}\hspace{0.17em}{\mathit{\sigma }}_{\mathit{s}}\mathrm{+}\mathrm{sgn}\hspace{0.17em}{\mathit{\sigma }}_{\sum }\mathrm{|}}\right)\right)\mathit{,}$0 ≤ χ_s ≤ 1, with the slope sum σ_∑ = ∑ σ_s and with σ_ϱ the slope in the ϱ-reconstruction. This choice of χ_s allows for an easy interpretation: First we note that flattening is selective, i.e., if σ_∑> 0 (σ_∑< 0) only components with positive (negative) slope are changed by the same factor, whereas species components with opposite slope are unaffected (χ_s = 1). Moreover, if sgn σ_∑ = sgn σ_ϱ and | σ_∑ | > | σ_ϱ | applying selective flattening recovers consistency in the sense that $\sum \begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{\sigma }}_{\mathit{s}}\end{array}\mathrm{=}{\mathit{\sigma }}_{\mathit{\varrho }}$. If | σ_∑ | < | σ_ϱ | no flattening is necessary and χ_s = 1. On the other hand if sgn σ_∑ ≠ sgn σ_ϱ, the total density slope can be expected to be close to zero. Then, the slope is set σ_ϱ ≡ 0 followed by species flattening as before, which also makes $\sum \begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{\sigma }}_{\mathit{s}}\end{array}\mathrm{=}\mathrm{0}$. We note that this does not generally mean setting $\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{\sigma }}_{\mathit{s}}\end{array}\mathrm{\equiv }\mathrm{0}$, which would imply reduction to a first-order algorithm. The choice of χ_s may be thought of as a best choice introducing a minimal amount of dissipation where needed in order to let | σ_∑ | converge to | σ_ϱ | from above. By contrast, convergence to | σ_ϱ | from below is suppressed, which would require profile steepening. Finally, number density fluxes $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{F}}^{\mathrm{\left(}{\mathit{n}}_{\mathit{s}}\mathrm{\right)}}\end{array}$ obtained on the basis of the modified reconstruction with slopes $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{\sigma }}_{\mathit{s}}\end{array}$ are rescaled to new values according to

${\mathit{F}}^{\mathrm{\left(}{\mathit{n}}_{\mathit{s}}\mathrm{\right)}}\mathrm{=}\frac{{\mathit{F}}^{\mathrm{\left(}\mathit{\varrho }\mathrm{\right)}}}{{\mathit{m}}_{\mathit{u}}\sum \begin{array}{c}\mathrm{˜}\\ {\mathit{\mu }}_{\mathit{s}}\end{array}\begin{array}{c}\mathrm{˜}\\ {\mathit{F}}^{\mathrm{\left(}{\mathit{n}}_{\mathit{s}}\mathrm{\right)}}\end{array}}\stackrel{}{{{\mathit{F}}^{\mathrm{\left(}{\mathit{n}}_{\mathit{s}}\mathrm{\right)}}}_{\mathrm{˜}}}$which ensures conservation, i.e., $\sum {\mathit{m}}_{\mathit{u}}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{\mu }}_{\mathit{s}}\end{array}{\mathit{F}}^{\mathrm{\left(}{\mathit{n}}_{\mathit{s}}\mathrm{\right)}}\mathrm{=}{\mathit{F}}^{\mathrm{\left(}\mathit{\varrho }\mathrm{\right)}}$ where F^(ϱ) is the total density mass flux. As a side note, in the similarly conservative multispecies FV advection scheme of Plewa & Müller (1999) the principle of selective flattening is applied to the upwinded cell interface values followed by a rescale of species fluxes. This is an equally valid approach applied to their third-order PPM scheme.

In summary, with updating ϱ mass is conserved up to machine round-off error accumulation as usual, whereas consistency of the gas species with ϱ in a FV sense is achieved, i.e., Eq. (5) – when volume-averaged over a grid cell – holds at t^{n + 1} when fulfilled at tⁿ.

For the source step integration a stiff ODE solver is needed. There is a wide variety of such solvers available in the mathematical literature which vary in accuracy and computational complexity from the simple first-order implicit Euler scheme to very precise high-order Peer methods. I decided to implement the Rosenbrock-Wanner method named GRK4(3)T developed in Kaps & Rentrop (1979). It is a fourth-order generalized Runge-Kutta method belonging to the class of semi-implicit schemes. This choice of solver was guided by the following thoughts:

• (1)

  Since the ODE solver is not envisaged as a stand-alone tool butis used operator-split together with second-order FVgasdynamics moderate accuracy requirements are sufficient,i.e., a fourth-order scheme like GRK4(3)T is convenient.

• (2)

  Since the time-step delivered to the ODE solver module can be much larger than the smallest timescale in the reaction system, the solver should have an efficient internal step-size control to meet accuracy requirements. GRK4(3)T embeds a third-order solution approximation almost for free, which can be used in a standard error estimator with user-prescribed relative (RTOL) and absolute (ATOL) error tolerances (see, e.g., Strehmel et al. 2012). In short, if u⁽⁴⁾ and u⁽³⁾ denote the fourth-order and third-order approximations (see Table 1), an estimate for the normalized error is given by $\mathit{ϵ}\mathrm{=}\underset{\mathit{s}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{1}\mathit{,}{\mathit{N}}_{\mathit{s}}\mathrm{+}\mathrm{1}}{\max }\begin{array}{c}\mathrm{⎧}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎨}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎩}\end{array}\frac{\mathrm{|}{\mathit{u}}_{\mathit{s}}^{\mathrm{\left(}\mathrm{4}\mathrm{\right)}}\mathrm{-}{\mathit{u}}_{\mathit{s}}^{\mathrm{\left(}\mathrm{3}\mathrm{\right)}}\mathrm{|}}{{\mathrm{ATOL}}_{\mathit{s}}\mathrm{+}\mathrm{RTOL}\mathrm{·}\max \mathrm{\left\{}\mathrm{\right|}{\mathit{u}}_{\mathit{s}}^{\mathrm{\left(}\mathrm{4}\mathrm{\right)}}\mathrm{\left|}\mathit{,}\mathrm{\right|}{\mathit{u}}_{\mathit{s}}^{\mathrm{\left(}\mathrm{3}\mathrm{\right)}}\mathrm{\left|}\mathrm{\right\}}}\begin{array}{c}\mathrm{⎫}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎬}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎭}\end{array}\mathrm{·}$If ϵ ≤ 1 the integration step with internal step-size h is accepted, whereas for ϵ> 1 it is rejected and is repeated with smaller step-size. In both cases the choice of a new step-size is error-controlled and adjusted according to $\mathit{h}\mathrm{-}\mathrm{\to }\mathit{h}\mathrm{·}\begin{array}{c}\min \end{array}\mathrm{\left\{}\mathrm{10}\mathit{,}\begin{array}{c}\max \end{array}\mathrm{\right\{}\mathrm{0.1}\mathit{,}\mathrm{0.9}\mathrm{·}{\mathit{ϵ}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathit{p}}\mathrm{\left\}}\mathrm{\right\}}$where p = 4 is the order of the scheme. All species components have the same value of ATOL whereas it is different, obviously, for the temperature component.

• (3)

  The species advection step was designed consistently with the gas density advection to ensure mass conservation. In order to maintain this property for the combined solver, the ODE method must conserve total species mass during reactions, i.e., recalling (7), $\frac{\mathit{d}}{\mathit{dt}}\sum _{\mathit{s}}\begin{array}{c}\mathrm{˜}\\ {\mathit{\mu }}_{\mathit{s}}\end{array}{\mathit{n}}_{\mathit{s}}\mathrm{=}\sum _{\mathit{s}}\begin{array}{c}\mathrm{˜}\\ {\mathit{\mu }}_{\mathit{s}}\end{array}{\mathit{R}}_{\mathit{s}}\mathrm{=}\mathrm{0}\hspace{0.17em}\mathit{.}$Our implementation of GRK4(3)T fulfills this condition to machine round-off⁵.

For interested readers, GRK4(3)T is sumarized in Table 1. It needs computation of the Jacobi matrix J = ∂S/∂u and the solution of a linear equation system of dimension N_s + 1 (number of species plus temperature) for four different right-hand side vectors corresponding to a four-stage procedure. At this stage of development, standard LU factorization is applied which is an (N_s + 1)³ operation and so becomes rapidly expensive for larger astrochemical networks. Costs may be reduced by utilizing a sparse structure of J inherent to astrochemistry. For several hundreds or thousands of species LU decomposition may still be much too expensive and, as a cheaper alternative, a Krylov-type solution approach could be applied. For the network in Table A.1 with N_s = 8 and N_r = 19, however, computation of J is the most expensive part and not the LU factorization. For efficiency reasons the Jacobian should be defined analytically, but the implementation also offers the possibility of an approximate evaluation by finite differencing if derivatives are too cumbersome.

In flow regions where the physical conditions are such that the ODE system is non-stiff, using an explicit method may be more efficient⁶. Therefore, a hybrid integration strategy was designed based on a measure of stiffness. The ODE system is advanced in time according to

$\mathit{h}\mathrm{·}\mathrm{\left|}\mathrm{\right|}\mathrm{J}\mathrm{\left|}\mathrm{\right|}\{\begin{array}{c}\\ \mathrm{\ge }\mathrm{2} \mathrm{useGRK}\mathrm{4}\mathrm{\left(}\mathrm{3}\mathrm{\right)}\mathrm{T}\\ \mathit{<}\mathrm{2} \mathrm{useDOPRI}\mathrm{5}\mathrm{\left(}\mathrm{4}\mathrm{\right)}\end{array}$where the norm⁷ of the Jacobi matrix is defined by

$\mathrm{\left|}\mathrm{\right|}\mathrm{J}\mathrm{\left|}\mathrm{\right|}\mathrm{=}\begin{array}{c}\max \end{array}\left\{\underset{\mathit{i}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{1}\mathit{,}{\mathit{N}}_{\mathit{s}}}{\max }\sum _{\mathit{j}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathit{s}}}\mathrm{\left|}{\mathit{J}}_{\mathit{i,j}}\mathrm{\right|}\mathit{,}\mathrm{\left|}{\mathit{J}}_{{\mathit{N}}_{\mathit{s}}\mathrm{+}\mathrm{1}\mathit{,}{\mathit{N}}_{\mathit{s}}\mathrm{+}\mathrm{1}}\mathrm{\right|}\right\}\hspace{0.17em}\mathit{.}$The explicit solver DOPRI5(4) is the well-known fifth-order accurate Runge-Kutta method of Dormand & Prince (1980). The value of two is roughly the stability limit of DOPRI5(4).

In order to better control the coupling between advection and reaction an a posteriori reaction time-step is introduced given by

$\mathrm{\Delta }{\mathit{t}}_{\mathrm{reac}}^{\mathrm{n}\mathrm{+}\mathrm{1}}\mathrm{=}\frac{\mathrm{\Delta }{\mathit{t}}^{\mathrm{n}}}{\mathrm{2}}\begin{array}{c}\min \end{array}\left({\mathit{ϵ}}_{\mathit{T}}\frac{{\mathit{T}}^{\mathrm{n}\mathrm{+}\mathrm{1}\mathit{,}\mathit{\star }}}{\mathrm{|}{\mathit{T}}^{\mathrm{n}\mathrm{+}\mathrm{1}}\mathrm{-}{\mathit{T}}^{\mathrm{n}\mathrm{+}\mathrm{1}\mathit{,}\mathit{\star }}\mathrm{|}}\mathit{,}{\mathit{ϵ}}_{\mathit{n}}\underset{\mathit{s}\mathrm{=}\mathrm{1}\mathit{,}{\mathit{N}}_{\mathit{s}}}{\min }\begin{array}{c}\mathrm{⎧}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎨}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎩}\end{array}\frac{{\mathit{n}}_{\mathit{s}}^{\mathrm{n}\mathrm{+}\mathrm{1}\mathit{,}\mathit{\star }}}{\mathrm{|}{\mathit{n}}_{\mathit{s}}^{\mathrm{n}\mathrm{+}\mathrm{1}}\mathrm{-}{\mathit{n}}_{\mathit{s}}^{\mathrm{n}\mathrm{+}\mathrm{1}\mathit{,}\mathit{\star }}\mathrm{|}}\begin{array}{c}\mathrm{⎫}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎬}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎭}\end{array}\right)\mathrm{·}$It means that the follow-up reaction time-step is limited by a maximum allowed relative change in species number densities (ϵ_n) due to chemical kinetics and relative change in temperature (ϵ_T) due to thermal processes computed from the last reaction source update in the Strang-splitting procedure. The coupling parameters ϵ_n and ϵ_T are user-prescribed quantities. Naturally, the values ϵ_n< 1,ϵ_T< 1 should be chosen with smaller values giving a better coupling at the expense of possibly reducing the overall time-step and, hence, increasing computational costs.

4. Validation and tests

4.1. Advection scheme validation

In order to validate the multispecies advection scheme proposed in this work an extended interacting blast wave problem was studied using the setup of Plewa & Müller (1999). The initial state is given by

$\mathrm{\left(}\mathit{\varrho ,p,v}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ \mathrm{\left(}\mathrm{1}\mathit{,}\mathrm{1000}\mathit{,}\mathrm{0}\mathrm{\right)}& \mathit{x}\mathrm{\le }\mathrm{0.1}\\ \mathrm{\left(}\mathrm{1}\mathit{,}\mathrm{0.01}\mathit{,}\mathrm{0}\mathrm{\right)}& \mathrm{0.1}\mathit{<}\mathit{x}\mathrm{\le }\mathrm{0.9}\\ \mathrm{\left(}\mathrm{1}\mathit{,}\mathrm{100}\mathit{,}\mathrm{0}\mathrm{\right)}& \mathit{x}\mathit{>}\mathrm{0.9}\end{array}$on the interval x = [0,1] discretized by 400 equidistant cells and adopting reflecting boundary conditions on both ends. Additionally, three nonreactive species components are introduced that have initial mass fractions of $\begin{array}{ccc}{\mathit{f}}_{\mathrm{1}}& \mathrm{=}& \mathrm{0.5}{\mathit{x}}^{\mathrm{2}}\\ {\mathit{f}}_{\mathrm{2}}& \mathrm{=}& \mathrm{0.5}{\sin }^{\mathrm{2}}\mathrm{\left(}\mathrm{20}\mathit{x}\mathrm{\right)}\\ {\mathit{f}}_{\mathrm{3}}& \mathrm{=}& \mathrm{1}\mathrm{-}{\mathit{f}}_{\mathrm{1}}\mathrm{-}{\mathit{f}}_{\mathrm{2}}\end{array}$as illustrated in Fig. 1 (left plot). An ideal equation of state with adiabatic index γ = 1.4 is used.

Figures 1 (right plot) and 2 illustrate the results at the final time t = 0.038. Figure 2 shows the solution for the hydrodynamic variables ϱ, p, and v. The spatial distribution of mass fractions in Fig. 1 agree well with those in Plewa & Müller (1999) (cf. Fig. 5) although their third-order and less dissipative PPM algorithm produces somewhat higher-amplitude variations in the region ${\mathrm{0.6}}_{\stackrel{\mathit{<}}{\mathrm{~}}}{\mathit{x}}_{\stackrel{\mathit{<}}{\mathrm{~}}}\mathrm{0.75}$. Moreover, no artificial discontinuity is created in the components f₂ and f₃ at x ≈ 0.6, which is visible in their solution obtained with the preferred sCMA method. As expected, the deviation from Eq. (5) over the integration time remains well bounded. Expressed in numbers for the mass fractions $\begin{array}{ccc}\mathrm{|}{\mathit{f}}_{\mathrm{1}}\mathrm{+}{\mathit{f}}_{\mathrm{2}}\mathrm{+}{\mathit{f}}_{\mathrm{3}}\mathrm{-}\mathrm{1}{\mathrm{|}}_{\max }& \mathit{<}& \mathrm{1.2}\mathrm{\times }{\mathrm{10}}^{-13}\mathit{,}\\ \mathrm{|}{\mathit{f}}_{\mathrm{1}}\mathrm{+}{\mathit{f}}_{\mathrm{2}}\mathrm{+}{\mathit{f}}_{\mathrm{3}}\mathrm{-}\mathrm{1}{\mathrm{|}}_{\mathrm{avg}}& \mathit{<}& \mathrm{2}\mathrm{\times }{\mathrm{10}}^{-14}\end{array}$measured on the grid for the maximum and spatially averaged error, and clearly demonstrating the conservative nature of the multispecies advection scheme.

Fig. 1 Spatial distribution of initial (left) and final (right) mass fractions (f1 – solid line; f2 – dashed line; f3 – dash-dotted line) in the interacting blast wave problem.

Fig. 2 Spatial distribution of density (top), pressure (middle), and velocity (bottom) at t = 0.038 for the multispecies interacting blast wave simulation.

4.2. Stiff solver validation

The accuracy of the GRK4(3)T solver and its correct implementation into the NIRVANA code was checked by performing several stiff problems. Here, the well-known problem named OREGO is presented. The OREGO problem is described in Hairer & Wanner (1996) and was originally developed by Field et al. (1972) as a simplified mathematical model for the oscillating behavior in the Belousov–Zhabotinskii chemical reaction (Gray 2002).

The problem is written in the form $\begin{array}{ccc}\mathit{ṅ}\mathrm{1}& \mathrm{=}& \mathrm{-}{\mathit{k}}_{\mathrm{1}}{\mathit{n}}_{\mathrm{1}}{\mathit{n}}_{\mathrm{2}}\mathrm{+}{\mathit{k}}_{\mathrm{2}}{\mathit{n}}_{\mathrm{1}}\mathrm{-}{\mathit{k}}_{\mathrm{3}}{\mathit{n}}_{\mathrm{1}}^{\mathrm{2}}\mathrm{+}{\mathit{k}}_{\mathrm{4}}{\mathit{n}}_{\mathrm{2}}\\ \mathit{ṅ}\mathrm{2}& \mathrm{=}& \mathrm{-}{\mathit{k}}_{\mathrm{1}}{\mathit{n}}_{\mathrm{1}}{\mathit{n}}_{\mathrm{2}}\mathrm{-}{\mathit{k}}_{\mathrm{4}}{\mathit{n}}_{\mathrm{2}}\mathrm{+}{\mathit{k}}_{\mathrm{5}}{\mathit{n}}_{\mathrm{3}}\\ \mathit{ṅ}\mathrm{3}& \mathrm{=}& {\mathit{k}}_{\mathrm{2}}{\mathit{n}}_{\mathrm{1}}\mathrm{-}{\mathit{k}}_{\mathrm{5}}{\mathit{n}}_{\mathrm{3}}\mathit{,}\end{array}$for concentrations n₁,n₂,n₃ with rate coefficients $\begin{array}{ccc}{\mathit{k}}_{\mathrm{1}}& \mathrm{=}& \mathrm{2.57555802}\mathrm{\times }{\mathrm{10}}^{\mathrm{8}}\\ {\mathit{k}}_{\mathrm{2}}& \mathrm{=}& \mathrm{7.72667406}\mathrm{\times }{\mathrm{10}}^{\mathrm{1}}\\ {\mathit{k}}_{\mathrm{3}}& \mathrm{=}& \mathrm{1.28777901}\mathrm{\times }{\mathrm{10}}^{\mathrm{7}}\\ {\mathit{k}}_{\mathrm{4}}& \mathrm{=}& \mathrm{1.29421790}\mathrm{\times }{\mathrm{10}}^{-2}\\ {\mathit{k}}_{\mathrm{5}}& \mathrm{=}& \mathrm{1.60972376}\mathrm{\times }{\mathrm{10}}^{-1}\mathit{.}\end{array}$Initial concentrations, obtained concentrations with GRK4(3)T at the final integration time t = 360, and high-precision reference values taken from the literature are given in Table 2. The numerical solution is illustrated in Fig. 3 demonstrating its oscillating character and the occurrence of very sharp variations over short time intervals. The maximum relative error obtained with GRK4(3)T is on the order of 5 × 10^-4 in the variable n₃ adopting numerical parameters RTOL= 10^-4 and ATOL= 10^-14.

Fig. 3 Numerical solution of the OREGO problem (with logarithmic scaling for n2).

4.3. Advective OREGO

This challenging test problem combines the previous OREGO problem with an underlying advection in a 1D setup. Choosing a constant velocity of v = 1 the OREGO equations in Sect. 4.2 are solved together with Eq. (14) on a domain x = [0,400] with initial conditions n₀(x,t = 0) = n_OREGO(t = 0) from Table 2. Outflow conditions are adopted on the right side of the domain and inflow conditions on the left side of the domain where concentrations are fixed at their initial values. The system is evolved up to time t = 360 corresponding to the integration time in the pure OREGO problem. The exact solution at the final time is given by n(x,t = 360) = n_OREGO(x/v) for x ≤ 360 and a constant continuation for x> 360, i.e., the spatial profiles of concentrations look identical to the temporal variation in the pure OREGO problem. We note that since the OREGO equations represent a nonconservative system (∑ _sṅ_s ≠ 0), this problem cannot serve as a consistency check of the multispecies advection scheme because the sum of concentrations does not obey a continuity equation.

The obtained numerical solution within the second-order operator-split framework with second-order advection scheme and GRK4(3)T solver is shown in Fig. 4 for two different grids with 400 and 1200 grid cells. A significant dependence on the numerical resolution is observed due to the presence of very sharp, thin structures advected with the flow. Numerical dissipation tends to broaden features and therefore to lower peak values. The largest of these errors are visible at the minimum peaks in variable n₂ in the lower resolved run. In this case insufficient resolution also causes a dispersion-like error noticeable as a moderate leftward shift of the right-hand feature. This dispersion error results from the interaction of numerical dissipation errors with the update of rate equations in virtue of the operator-split approximation. The more highly resolved case produces results much closer to the reference solution. The outcome of both low- and high-resolution runs is not sensitive to the choice of coupling parameter ϵ_n in the expression for Δt_reac because there is no back-reaction of the species evolution on the flow.

Fig. 4 Numerical solution of the OREGO problem with underlying advection for different spatial resolution (solid: 1200 grid cells; dashed: 400 grid cells).

4.4. Ionization equilibrium

The relaxation of a hydrogen-helium plasma with species H, H⁺, He, He⁺, He²⁺ in ionization equilibrium is considered. Assuming a fixed temperature of 10⁵ K and initial number densities far from equilibrium for this temperature, namely,

${\mathit{n}}_{\mathrm{H}}^{\mathrm{0}}\mathrm{=}{\mathit{n}}_{\mathrm{He}}^{\mathrm{0}}\mathrm{=}{\mathrm{10}}^{\mathrm{6}}\mathit{,}\hspace{0.17em}{\mathit{n}}_{{\mathrm{H}}^{\mathrm{+}}}^{\mathrm{0}}\mathrm{=}{\mathrm{10}}^{\mathrm{3}}\mathit{,}\hspace{0.17em}{\mathit{n}}_{{\mathrm{He}}^{\mathrm{+}}}^{\mathrm{0}}\mathrm{=}{\mathit{n}}_{{\mathrm{He}}^{\mathrm{2}\mathrm{+}}}^{\mathrm{0}}\mathrm{=}\mathrm{0}$measured in units m^-3, the ionization/recombination network consisting of reactions 1–6 in Table A.1 is solved. The evolution of the network is followed up to a maximum time of t_max = 10⁶ yr which is more than two orders of magnitude larger then the relaxation time for the ionization of He⁺.

In equilibrium (t → ∞), where ṅ_H = ṅ_He = ṅ_H⁺ = ṅ_He⁺ = ṅ_{He^{2 +}} = 0, the number densities are given by

$\begin{array}{ccc}& & {\mathit{n}}_{\mathrm{H}\mathit{,}\mathrm{eq}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}{\mathit{k}}_{\mathrm{1}}\mathit{/}{\mathit{k}}_{\mathrm{2}}}\mathrm{(}{\mathit{n}}_{\mathrm{H}}^{\mathrm{0}}\mathrm{+}{\mathit{n}}_{{\mathrm{H}}^{\mathrm{+}}}^{\mathrm{0}}\mathrm{)}\\ & & {\mathit{n}}_{{\mathrm{H}}^{\mathrm{+}}\mathit{,}\mathrm{eq}}\mathrm{=}\frac{{\mathit{k}}_{\mathrm{1}}}{{\mathit{k}}_{\mathrm{2}}}{\mathit{n}}_{\mathrm{H}\mathit{,}\mathrm{eq}}\\ & & {\mathit{n}}_{\mathrm{He}\mathit{,}\mathrm{eq}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}{\mathit{k}}_{\mathrm{3}}\mathit{/}{\mathit{k}}_{\mathrm{4}}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{k}}_{\mathrm{5}}\mathit{/}{\mathit{k}}_{\mathrm{6}}\mathrm{)}}{\mathit{n}}_{\mathrm{He}}^{\mathrm{0}}\\ & & {\mathit{n}}_{{\mathrm{He}}^{\mathrm{+}}\mathit{,}\mathrm{eq}}\mathrm{=}\frac{{\mathit{k}}_{\mathrm{3}}}{{\mathit{k}}_{\mathrm{4}}}{\mathit{n}}_{\mathrm{He}\mathit{,}\mathrm{eq}}\\ & & {\mathit{n}}_{{\mathrm{He}}^{\mathrm{2}\mathrm{+}}\mathit{,}\mathrm{eq}}\mathrm{=}\frac{{\mathit{k}}_{\mathrm{5}}}{{\mathit{k}}_{\mathrm{6}}}{\mathit{n}}_{{\mathrm{He}}^{\mathrm{+}}\mathit{,}\mathrm{eq}}\end{array}$with rate coefficients k_i as specified in Table A.1. Equilibrium values together with initial values and obtained numerical values are summarized in Table 3. The numerical solution of the relaxation process is shown in Fig. 5. Ionization of H takes place on a timescale of ≈50 yr followed by ionization of He on a timescale of ≈100 yr. Finally, built-up He⁺ is transformed into He²⁺ in a few thousand years. At t_max, equilibrium balance is reached with relative deviation less than 10^-8 from exact equilibrium values.

Fig. 5 Numerical solution of the ionization/recombination relaxation problem. Number densities of all species are shown as a function of time.

4.5. One-dimensional Riemann problem

As a first test in combining gasdynamics with chemistry, a 1D shock tube problem – a modified Sod problem – is considered. The chemical network as cataloged in Table A.1 is adopted accompanied by all the thermal processes summarized in Table B.1. The chemistry model is that used in Grassi et al. (2014) who performed the same Riemann problem to check their implementation of KROME in conjunction with the astrophysics codes RAMSES (Teyssier 2002) and FLASH (Fryxell et al. 2000).

The initial conditions are set up on a computational domain x = [0,1 pc] discretized by 1000 equidistant grid cells. Values on the regions left and right of the interface at x = 0.5 pc are given by

$\begin{array}{ccc}& & \mathrm{\left(}\mathit{\varrho ,T}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ & \\ & \end{array}\\ & & {v}\mathrm{=}\mathrm{0}\\ & & {\mathit{n}}_{\mathit{X}}\mathrm{=}{\mathit{f}}_{\mathit{X}}^{\mathrm{0}}\mathit{\varrho }\mathit{/}{\mathit{\mu }}_{\mathit{X}}{\mathit{m}}_{\mathit{u}}\end{array}$where ${\mathit{f}}_{\mathit{X}}^{\mathrm{0}}$ is the initial mass fraction of the different species X (see Table 4).

The outcome of this test is illustrated in Figs. 6 and 7. Figure 6 shows gasdynamical quantities and Fig. 7 shows the spatial variation of the mass fraction of species H₂, H⁻, and e⁻. In order to allow a qualitative comparison⁸ with the results of Grassi et al. (2014) the simulation was stopped at canonical time t = 2.45 × 10⁴ yr (see their Figs. 19–22). Differences from the pure gasdynamical case without chemistry are mainly due to the chemical evolution of H₂ and H₂-related cooling mechanisms which are the dominant thermal processes. It leads to a significant temperature reduction in regions of relatively high H₂ mass fraction. Finally, the overall characteristics of chemo-thermal evolution obtained with NIRVANA agree well with the behavior observed in Grassi et al. (2014)⁹.

Fig. 6 Gas density (top), temperature (middle), and velocity (bottom) as a function of distance (unit pc) for the modified Sod problem with chemistry. The dashed line corresponds to the solution of a standard Sod problem without chemistry.

Fig. 7 Mass fraction of species H2 (top), H− (middle), and e− (bottom) as a function of distance (unit pc) for the modified Sod problem with chemistry.

4.6. One-dimensional propagating shock

As a second example combining magneto-gasdynamics with chemistry a propagating shock problem is presented (see Teşileanu et al. 2008 and Massaglia et al. 2005). The background medium is a stratified gas with constant temperature T₀ = 1000 K, zero velocity v₀ = 0, and constant (perpendicular) magnetic field B_y0 = 10^-8 T. The background density profile is given by

${\mathit{\varrho }}_{\mathrm{0}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}{\mathit{\varrho }}_{\mathrm{00}}\frac{{\mathit{x}}_{\mathrm{0}}^{\mathrm{2}}}{{\mathit{x}}^{\mathrm{2}}\mathrm{+}{\mathit{x}}_{\mathrm{0}}^{\mathrm{2}}}$where ϱ₀₀ = 10¹¹m_u m^-3 and x₀ = 1.8 × 10¹² m. The computational domain ranges from x = 0 to x = 4 × 10¹³ m assuming outflow boundary conditions on both sides and adopting an effective spatial resolution¹⁰ of δx = 1.25 × 10⁹ m. The background medium is superimposed on a localized adiabatic perturbation that steepens into a single shock propagating down the stratification. In order to eliminate the reverse shock, the disturbance is constructed such that the non-magnetic Riemann invariant

$\mathit{v}\mathrm{-}\frac{\mathrm{2}}{\mathit{\gamma }\mathrm{-}\mathrm{1}}\sqrt{\mathit{\gamma p}\mathit{/}\mathit{\varrho }}\mathrm{=}\mathrm{const}\mathit{.}$is fulfilled with γ assumed 5/3. At t = 0 this results in a structure (background plus perturbation) given by $\begin{array}{ccc}\mathit{\varrho }\mathrm{\left(}\mathit{x}\mathrm{\right)}& \mathrm{=}& {\left[\frac{\mathit{\gamma }\mathrm{-}\mathrm{1}}{\mathrm{2}\sqrt{{\mathit{K}}_{\mathrm{0}}\mathit{\gamma }}}{\mathit{v}}^{\mathrm{\prime }}\mathrm{+}{\mathit{\varrho }}_{\mathrm{0}}^{\frac{\mathit{\gamma }\mathrm{-}\mathrm{1}}{\mathrm{2}}}\right]}^{\frac{\mathrm{2}}{\mathit{\gamma }\mathrm{-}\mathrm{1}}}\\ \mathit{p}\mathrm{\left(}\mathit{x}\mathrm{\right)}& \mathrm{=}& {\mathit{K}}_{\mathrm{0}}{\mathit{\varrho }}^{\mathit{\gamma }}\\ \mathit{v}\mathrm{\left(}\mathit{x}\mathrm{\right)}& \mathrm{=}& {\mathit{v}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{x}\mathrm{\right)}\end{array}$with ${\mathit{K}}_{\mathrm{0}}\mathrm{=}{\mathit{k}}_{\mathrm{B}}{\mathit{n}}_{\mathrm{tot}\mathit{,}\mathrm{0}}{\mathit{T}}_{\mathrm{0}}{\mathit{\varrho }}_{\mathrm{0}}^{\mathrm{-}\mathit{\gamma }}$ computed for the background medium. The velocity perturbation v′ is prescribed as

${\mathit{v}}^{\mathrm{\prime }}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ \frac{\widehat{{\mathit{v}}^{\mathrm{\prime }}}}{{\mathit{\sigma }}^{\mathrm{2}}}\mathrm{[}\mathrm{-}\mathrm{(}\mathit{x}\mathrm{-}{\mathit{x}}_{\mathrm{1}}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{2}\mathit{\sigma }\mathrm{(}\mathit{x}\mathrm{-}{\mathit{x}}_{\mathrm{1}}{\mathrm{)}}^{\mathrm{]}}& {\mathit{x}}_{\mathrm{1}}\mathit{<}\mathit{x}\mathit{<}{\mathit{x}}_{\mathrm{1}}\mathrm{+}\mathrm{2}\mathit{\sigma }\\ \mathrm{0}& \mathrm{otherwise}\end{array}$where the peak value $\widehat{{\mathit{v}}^{\mathrm{\prime }}}\mathrm{=}\mathrm{3}\mathrm{\times }{\mathrm{10}}^{\mathrm{4}}\hspace{0.17em}\mathrm{m}\mathit{/}\mathrm{s}$, σ = 2 × 10¹¹ m and x₁ = 10¹² m. We note that exterior to the interval (x₁,x₁ + 2σ) the configuration reduces to the isothermal background.

The gas is composed of species X with initial mass fractions ${\mathit{f}}_{\mathit{X}}^{\mathrm{0}}$ as summarized in Table 5. The chemical network in Table A.1 with thermal processes in Table B.1 is employed. Three different simulations are performed. The first includes the full chemistry model and all thermal processes, i.e., energy losses by ionization/recombination and line emission of H and He as well as H₂-related cooling/heating mechanisms. In the second simulation H₂ is treated as a passive scalar, H₂ chemistry is switched off, and there is no H₂ cooling/heating. It serves to isolate the effects of H₂ chemistry which are expected to be important in virtue of the relatively high initial H₂ mass fraction of ${\mathit{f}}_{{\mathrm{H}}_{\mathrm{2}}}^{\mathrm{0}}\mathrm{=}\mathrm{0.1}$ and chosen temperature T₀. This simulation is, at least qualitatively, comparable with the SNEq model in Teşileanu et al. (2008) despite the absence of metal cooling here. The third simulation is adiabatic with all species passively advected and no chemical reactions at all. In the simulations with a chemical reaction network included, the reaction time–step parameters are chosen ϵ_n = ϵ_T = 0.1, which adjusts it to any rate that is on the order of the dynamical time–step.

The results are shown in Fig. 8 for selected quantities at a time t = 13 yr where the simulations have been stopped. Evidently, the evolution of the shock and the structure of the post-shock region is crucially affected by the chemistry and cooling processes when compared with the adiabatic case. The shock speed is substantially reduced and the temperature drops faster behind the shock. Shock heating raises the temperature well above 10⁴ K, which partially ionizes the region with extent ≈3 × 10¹² m right behind the shock front. Also, the magnetic field reaches its maximum in this region. The influence of H₂ chemistry is to provide further cooling in regions where the number density of molecular hydrogen is sufficiently high. This is all of the region up to x ≈ 7 × 10¹² m and around the peak in molecular hydrogen density at x ≈ 1.3 × 10¹³ m where a thin shell has formed in the close vicinity of the left edge of ionized region. The shell lies within a broader region where the H₂ fraction otherwise was substantially reduced by collisional dissociation in passing the shock front – up to three orders of magnitude compared to the initial mass fraction. This feature is also clearly seen in the temperature that has dropped to ≈100 K there. The pre-shock region is at most marginally affected by cooling because of low effectiveness. Overall, the striking similarity in the dynamical evolution and spatial structure of the simulation without H₂ chemistry/cooling with the SNEq model in Teşileanu et al. (2008) further strengthens confidence in the validity of the presented algorithmic approach.

Fig. 8 From top to bottom: gas density, electron number fraction, H2 mass fraction, temperature, velocity and magnetic field at t = 13 yr (solid line: full chemistry; dashed line: no H2 chemistry; dotted line: no chemistry, adiabatic). The thin solid line is the initial profile.

5. Conclusions

A chemical reaction network solver was successfully implemented into the astrophysics code NIRVANA. For that purpose a numerical framework was devised that accounts for the accurate treatment of the set of advection-reaction equations for chemical species. These equations add to the already complex equations of magneto-gasdynamics. In particular, the ODE system built from the rate equations representing chemical kinetics and the temporal temperature change due to chemistry-connected thermal processes is solved coherently with a high-order stiff integrator of linear-implicit type. Another feature of the implementation is that it evolves chemical species consistently with the total gas density. This nice property is achieved by modifications in the advection part of the FV method and by a proper choice of ODE solver for the reaction part. It means that the mass conservation principle, which is usually fulfilled in a FV discretization of the gasdynamics equations, remains valid here for the case of reactive flows.

A series of tests have provided confidence in future applications of the code involving chemistry and complex thermal processes. However, computational costs for multidimensional simulations are a limiting factor in practice and, this time, prohibit the use of chemical networks that are too large. There is otherwise some room for optimizations of the employed ODE solver, for example by utilizing the sparse matrix structure to accelerate LU decomposition and/or by applying a Krylov-type approximation in order to reduce costs further. Such possible improvements are the subject of future work.

Complementary to the numerical framework an easy-to-use interface was developed which allows straightforward upgrades of the H-He chemical network that was constructed as a standard network in this work. Adding new species and reactions is performed just by editing a specifically formatted text file. By parsing that file the ODE system is generated automatically, i.e, the ODE system need not be given explicitly. There is, in principle, no restriction to the number of species, the number of reactions, and the number of reactants in a reaction. On the other hand, kinetic and thermal rate coefficients for each chemical reaction as well as special thermal processes have to be defined in separate modules. This allows the highest degree of flexibility regarding their dependence on temperature and, possibly, gas composition.

All in all, with the numerical and technical framework installed in this work the NIRVANA code has become a very powerful tool for time-dependent chemo-magneto-gasdynamical simulations. In combination with the code’s self-gravity solver, its efficient solver for dissipation processes and with its adaptive mesh refinement capability, challenging astrophysical problems seem attackable in the future limited more by computational costs than by physics.

Acknowledgments

I would like to thank the referee, Pierre Lesaffre, for valuable comments useful in improving the paper.

References

Appendix A: Chemical network

Table A.1 catalogs the chemical network of a hydrogen-helium plasma used in this work. The temperature-dependent rate coefficients for collisional ionization, radiative, and di-electronic recombination processes for H and He are collected from Janev et al. (1987) and Cen (1992). Other rates are mostly adopted from Abel et al. (1997) where the same network was used in a cosmological context to reflect a primordial gas composition.

Appendix B: Cooling and heating processes

All thermal processes that have been considered in conjunction with the chemical network in Table A.1 are summarized in Table B.1.

Appendix B.1: Hydrogen and helium cooling processes

Following Cen (1992) a standard description for cooling processes in a hydrogen-helium plasma is adopted. It includes ionization cooling of H, He, and He⁺ by collisions with electrons; radiative recombination cooling of H⁺, He⁺, and He²⁺; di-electronic recombination cooling of He⁺; collisional ionization cooling of the metastable He(2³S); and cooling due to collisional excitation of H (all atomic levels n), He (n = 2,3,4 triplets), and He⁺ (n = 2). Moreover, Bremsstrahlung for the ions H⁺, He⁺, and He²⁺ according to Black (1981) is taken into account. We note that the effective cooling rates for He(2³S) and He (n = 2,3,4 triplets) cooling are ∝${\mathit{n}}_{\mathrm{e}}^{\mathrm{2}}{\mathit{n}}_{{\mathrm{He}}^{\mathrm{+}}}$, hence, they cannot be treated as canonical source terms in (13).

Appendix B.2: Molecular hydrogen cooling and formation heating

Cooling due to collisional excitation of ro-vibrational levels of molecular hydrogen is considered. The net energy loss per volume per time is written in the form

${\mathit{L}}_{{\mathrm{H}}_{\mathrm{2}}}\mathrm{=}\mathrm{-}{\mathrm{\Lambda }}_{{\mathrm{H}}_{\mathrm{2}}}{\mathit{n}}_{{\mathrm{H}}_{\mathrm{2}}}\mathit{\tau }\mathrm{\left(}{n}\mathrm{\right)}$where Λ_H2 denotes the cooling function for an optically thin environment, and τ(n) is an opacity correction term once the gas enters the optically thick regime. Here, according to Ripamonti & Abel (2004)

$\mathit{\tau }\mathrm{=}\begin{array}{c}\min \end{array}\left[\mathrm{1}\mathit{,}{\left(\frac{{\mathit{n}}_{\mathrm{tot}}}{\mathrm{8}\mathrm{\times }{\mathrm{10}}^{\mathrm{15}}{\mathrm{m}}^{-3}}\right)}^{-0.45}\right]\mathrm{·}$The cooling function Λ_H2 measured in units [J/s] is given by

${\mathrm{\Lambda }}_{{\mathrm{H}}_{\mathrm{2}}}\mathrm{=}\frac{{\mathrm{\Lambda }}_{{\mathrm{H}}_{\mathrm{2}}\mathit{,}\mathrm{LTE}}}{\mathrm{1}\mathrm{+}{\mathrm{\Lambda }}_{{\mathrm{H}}_{\mathrm{2}}\mathit{,}\mathrm{LTE}}\mathit{/}{\mathrm{\Lambda }}_{{\mathrm{H}}_{\mathrm{2}}\mathit{,}\mathrm{CE}}}$which, in the high-density limit (Λ_H2,CE ≫ Λ_H2,LTE), reduces to the LTE cooling function $\begin{array}{ccc}{\mathrm{\Lambda }}_{{\mathrm{H}}_{\mathrm{2}}\mathit{,}\mathrm{LTE}}& \mathrm{=}& \mathrm{9.5}\mathrm{\times }{\mathrm{10}}^{-29}\frac{{\mathit{T}}_{\mathrm{3}}^{\mathrm{3.76}}}{\mathrm{1}\mathrm{+}\mathrm{0.12}{\mathit{T}}_{\mathrm{3}}^{\mathrm{2.1}}}\exp \mathrm{[}\mathrm{-}\mathrm{(}\mathrm{0.13}\mathit{/}{\mathit{T}}_{\mathrm{3}}{\mathrm{)}}^{\mathrm{3}}\mathrm{]}\\ & \mathrm{+}& \mathrm{3}\mathrm{\times }{\mathrm{10}}^{-31}\exp \mathrm{(}\mathrm{-}\mathrm{0.51}\mathit{/}{\mathit{T}}_{\mathrm{3}}\mathrm{)}\mathrm{+}\mathrm{6.7}\mathrm{\times }{\mathrm{10}}^{-26}\exp \mathrm{(}\mathrm{-}\mathrm{5.86}\mathit{/}{\mathit{T}}_{\mathrm{3}}\mathrm{)}\\ & \mathrm{+}& \mathrm{1.6}\mathrm{\times }{\mathrm{10}}^{\mathrm{11}}\exp \mathrm{(}\mathrm{-}\mathrm{11.7}\mathit{/}{\mathit{T}}_{\mathrm{3}}\mathrm{)}\mathit{,} {\mathit{T}}_{\mathrm{3}}\mathrm{=}\mathit{T}\mathit{/}{\mathrm{10}}^{\mathrm{3}}\mathrm{K}\mathit{,}\end{array}$following the work of Hollenbach & McKee (1979). In the low-density limit, Λ_H2,CE ≪ Λ_H2,LTE, the cooling function of Glover & Abel (2008) is adopted. It represents a sum of contributions from collisional excitation of H₂with collision partners X_k ={H, H⁺, He, H₂,e⁻}:

${\mathrm{\Lambda }}_{{\mathrm{H}}_{\mathrm{2}}\mathit{,}\mathrm{CE}}\mathrm{=}\sum _{\mathit{k}}{\mathrm{\Lambda }}_{{\mathrm{H}}_{\mathrm{2}}\mathit{,k}}{\mathit{n}}_{\mathit{k}}\mathit{.}$The logarithms of individual cooling rates are polynomial functions in log T, i.e., log Λ_H2,k = ∑ _lc_l(log T)^l, with fitting coefficients c_l taken from Table 8 in Glover & Abel (2008). In this fitting approach no distinction is made between excitation of ortho–H₂and para–H₂; instead, an ortho–para ratio of 3:1 is assumed.

Collisional dissociation of molecular hydrogen described by reactions 12 and 13 in Table A.1 is an endothermic process that leads to gas cooling. According to Omukai (2000) the energy removed from the gas per dissociated molecule is 4.48eV for both reactions. The thermal rate coefficients for these reactions are in canonical form and are given by

${\mathit{\ell }}_{\mathrm{12}\mathit{,}\mathrm{13}}\mathrm{=}\mathrm{-}\mathrm{4.48}\mathrm{eV}\mathrm{\times }{\mathit{k}}_{\mathrm{12}\mathit{,}\mathrm{13}}$with k_12,13 the kinetic rate coefficients.

In contrast, molecular hydrogen formation via associative detachment (reaction 8) and charge transfer (reaction 10) are exothermic processes that partly heat the gas. According to Omukai (2000), the binding energy released per produced molecule is 3.53 eV (1.83 eV) for reaction 8 (10). The thermal rate coefficients are given by

$\begin{array}{ccc}& & {\mathit{\ell }}_{\mathrm{8}}\mathrm{=}\mathrm{3.53}\mathrm{eV}\mathrm{\times }{\mathit{f}}_{\mathrm{cr}}{\mathit{k}}_{\mathrm{8}}\\ & & {\mathit{\ell }}_{\mathrm{10}}\mathrm{=}\mathrm{1.83}\mathrm{eV}\mathrm{\times }{\mathit{f}}_{\mathrm{cr}}{\mathit{k}}_{\mathrm{10}}\end{array}$where

${\mathit{f}}_{\mathrm{cr}}\mathrm{=}{\left[\mathrm{1}\mathrm{+}\frac{{\mathrm{10}}^{\mathrm{6}}{\mathit{T}}^{-0.5}}{\mathrm{1.6}{\mathit{n}}_{\mathrm{H}}\exp \left(\mathrm{-}{\mathrm{(}{\frac{\mathrm{400}}{\mathit{T}}}^{\mathrm{)}}}^{\mathrm{2}}\right)\mathrm{+}\mathrm{1.4}{\mathit{n}}_{{\mathrm{H}}_{\mathrm{2}}}\exp \mathrm{(}\mathrm{-}{\frac{\mathrm{12000}}{\mathit{T}\mathrm{+}\mathrm{1200}}}^{\mathrm{)}}}\right]}^{-1}\mathit{,}$according to Hollenbach & McKee (1979), approximates the fraction of liberated binding energy going into heating the gas. The factor f_cr arises from the assumption that all of the binding energy is first put into rotational/vibrational excitation of H₂which is then partly converted into heat via collisional de-excitation rather than radiated away. However, these heat sources only become important at high densities (n ≳ 10¹⁴ m^-3, say) at which f_cr approaches 1.

All Tables

All Figures

	Fig. 1 Spatial distribution of initial (left) and final (right) mass fractions (f1 – solid line; f2 – dashed line; f3 – dash-dotted line) in the interacting blast wave problem.
In the text

	Fig. 2 Spatial distribution of density (top), pressure (middle), and velocity (bottom) at t = 0.038 for the multispecies interacting blast wave simulation.
In the text

	Fig. 3 Numerical solution of the OREGO problem (with logarithmic scaling for n2).
In the text

	Fig. 4 Numerical solution of the OREGO problem with underlying advection for different spatial resolution (solid: 1200 grid cells; dashed: 400 grid cells).
In the text

	Fig. 5 Numerical solution of the ionization/recombination relaxation problem. Number densities of all species are shown as a function of time.
In the text

	Fig. 6 Gas density (top), temperature (middle), and velocity (bottom) as a function of distance (unit pc) for the modified Sod problem with chemistry. The dashed line corresponds to the solution of a standard Sod problem without chemistry.
In the text

	Fig. 7 Mass fraction of species H2 (top), H− (middle), and e− (bottom) as a function of distance (unit pc) for the modified Sod problem with chemistry.
In the text

	Fig. 8 From top to bottom: gas density, electron number fraction, H2 mass fraction, temperature, velocity and magnetic field at t = 13 yr (solid line: full chemistry; dashed line: no H2 chemistry; dotted line: no chemistry, adiabatic). The thin solid line is the initial profile.
In the text