Issue 
A&A
Volume 586, February 2016



Article Number  A82  
Number of page(s)  13  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/201527262  
Published online  28 January 2016 
A chemical reaction network solver for the astrophysics code NIRVANA
LeibnizInstitut für Astrophysik Potsdam, An der Sternwarte 16, 14482 Potsdam, Germany
email: uziegler@aip.de
Received: 27 August 2015
Accepted: 23 November 2015
Context. Chemistry often plays an important role in astrophysical gases. It regulates thermal properties by changing species abundances and via ionization processes. This way, timedependent cooling mechanisms and other chemistry–related energy sources can have a profound influence on the dynamical evolution of an astrophysical system. Modeling those effects with the underlying chemical kinetics in realistic magnetogasdynamical simulations provide the basis for a better link to observations.
Aims. The present work describes the implementation of a chemical reaction network solver into the magnetogasdynamical code NIRVANA. For this purpose a multispecies structure is installed, and a new module for evolving the rate equations of chemical kinetics is developed and coupled to the dynamical part of the code. A small chemical network for a hydrogen–helium plasma was constructed including associated thermal processes which is used in test problems.
Methods. Evolving a chemical network within timedependent simulations requires the additional solution of a set of coupled advectionreaction equations for species and gas temperature. Secondorder Strangsplitting is used to separate the advection part from the reaction part. The ordinary differential equation (ODE) system representing the reaction part is solved with a fourthorder generalized RungeKutta method applicable for stiff systems inherent to astrochemistry.
Results. A series of tests was performed in order to check the correctness of numerical and technical implementation. Tests include wellknown stiff ODE problems from the mathematical literature in order to confirm accuracy properties of the solver used as well as problems combining gasdynamics and chemistry. Overall, very satisfactory results are achieved.
Conclusions. The NIRVANA code is now ready to handle astrochemical processes in timedependent simulations. An easytouse interface allows implementation of complex networks including thermal processes. In combination with NIRVANA’s selfgravity solver, its efficient solver for dissipation terms and its adaptive mesh refinement capability challenging astrophysical problems come into reach with the code.
Key words: magnetohydrodynamics (MHD) / astrochemistry / methods: numerical
© ESO, 2016
1. Introduction
Astrochemistry deals with chemical reactions in astrophysical gases and, taking a broader view, embraces ionizing processes, dustmolecule interactions, and photochemical 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 prestellar 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_{3} molecule as an indicator of the gas temperature.
In order to bring magnetogasdynamical simulations and observations closer, astrochemistry must find its way into numerical models. However, embedding chemistry in timedependent, 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 chemomagnetogasdynamical 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 timestep, 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 magnetogasdynamical part for larger chemical networks. In consequence, this often prohibits the use of complex chemical networks in highresolution 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 timestep 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 timedependent 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 magnetogasdynamics 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 magnetogasdynamical codes like ASTROBEAR (Poludnenko et al. 2005) or PLUTO (Teşileanu et al. 2008). Moreover, there are powerful 1D timedependent codes that integrate magnetogasdynamics alongside more complex chemical networks used, e.g., in highresolution studies of magnetic shocks (Lesaffre et al. 2004) or diffusion fronts (Lesaffre et al. 2007). Standalone astrochemistry packages (not coupled to dynamics) have been designed for special purposes like the ALCHEMIC code (Semenov et al. 2010) for the chemistry in protoplanetary 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 radiationdominated 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 generalpurpose 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 chemistryrelated thermal processes, photochemistry and dust evolution.
This paper describes the implementation of a chemical reaction network solver into the astrophysics code NIRVANA^{1} (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 (1)between species 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 (2)expressing total charge conservation^{2}. 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_{2}+ 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 byhand 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^{2+}, H^{−}, H_{2}, H, and free electrons. The network includes ionization/recombination of H, He and a channel for H_{2} formation via the H^{−}, H ions.
2.1. Advectionreaction system
In gasdynamical simulations the time evolution of species number densities is described by a coupled system of advection–reaction equations, (3)with reaction term (4)subject to constraint (2). The product on the righthand 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 magnetogasdynamical 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
with n_{tot} = ∑ _{s}n_{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 (5)with and the electron weight μ_{e} = 5.49 × 10^{4}. The mean molecular weight and electron number fraction Y is defined by
Usually, in gasdynamical simulations a continuity equation for ϱ expressing mass conservation is evolved in time: (6)Multiplying Eq. (3) by and summing over species results in (7)Hence, consistency with the continuity equation requires (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_{2}, 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 (9)holds for the thermal energy density e = p/ (γ−1) with adiabatic index γ assumed to be a constant value^{3}. Since the full system of chemomagnetogasdynamical equations is solved in an operatorsplit fashion, from Eqs. (3) and (9) the thermoreactive part in the solution procedure (ϱ, v, and B are assumed to be constant during this update) is governed by
Using the equation of state, the expression for the total number density and Eq. (10) one has
It follows for the evolution of temperature (12)which is used instead of (11) in the operatorsplit 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 (13)The first (canonical) term on the righthand side has the same structural form as the source term R_{s} in the advectionreaction 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 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 chemomagnetogasdynamical equations are not described here. In this coverage numerical integration of the coupled system of advectionreaction Eq. (3) and temperature Eq. (12) must be carried out with respect to coexisting 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 (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 (15)where . In order to advance the solution one timestep Δt^{n} from t^{n} to t^{n + 1} = t^{n} + Δt^{n} a secondorder accurate Strangtype operator splitting is used. Denoting the advection and source step operators as and the update formally reads
The actual size of timestep is given by
where is the advection/dynamical timestep dictated by the usual CFL condition for numerical stability of the discretized gasdynamical equations and is a reaction timestep 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 as sole timestep may be insufficient because large splitting errors can lead to unphysical behavior or even destabilize numerics in situations where .
The advection of individual species is performed within the FiniteVolume (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 volumeaveraged number densities requires computing fluxes at each grid cell surface, which depend on a properly upwinded state . This is derived in accordance with other fluid variables from the approximate solution of a Riemann problem with leftsided/rightsided states at grid cell surfaces. These are obtained by a piecewise linear reconstruction procedure based on volumeaveraged 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
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^{4}. 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 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 secondorder accuracy. Concretely, the individual species slopes in each spatial direction and in each numerical cell are modified according to
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 . 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 . We note that this does not generally mean setting , which would imply reduction to a firstorder 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 obtained on the basis of the modified reconstruction with slopes are rescaled to new values according to
which ensures conservation, i.e., 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 thirdorder PPM scheme.
In summary, with updating ϱ mass is conserved up to machine roundoff error accumulation as usual, whereas consistency of the gas species with ϱ in a FV sense is achieved, i.e., Eq. (5) – when volumeaveraged over a grid cell – holds at t^{n + 1} when fulfilled at t^{n}.
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 firstorder implicit Euler scheme to very precise highorder Peer methods. I decided to implement the RosenbrockWanner method named GRK4(3)T developed in Kaps & Rentrop (1979). It is a fourthorder generalized RungeKutta method belonging to the class of semiimplicit schemes. This choice of solver was guided by the following thoughts:

(1)
Since the ODE solver is not envisaged as a standalone tool butis used operatorsplit together with secondorder FVgasdynamics moderate accuracy requirements are sufficient,i.e., a fourthorder scheme like GRK4(3)T is convenient.

(2)
Since the timestep 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 stepsize control to meet accuracy requirements. GRK4(3)T embeds a thirdorder solution approximation almost for free, which can be used in a standard error estimator with userprescribed relative (RTOL) and absolute (ATOL) error tolerances (see, e.g., Strehmel et al. 2012). In short, if u^{(4)} and u^{(3)} denote the fourthorder and thirdorder approximations (see Table 1), an estimate for the normalized error is given by If ϵ ≤ 1 the integration step with internal stepsize h is accepted, whereas for ϵ> 1 it is rejected and is repeated with smaller stepsize. In both cases the choice of a new stepsize is errorcontrolled and adjusted according to 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), Our implementation of GRK4(3)T fulfills this condition to machine roundoff^{5}.
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 righthand side vectors corresponding to a fourstage procedure. At this stage of development, standard LU factorization is applied which is an (N_{s} + 1)^{3} 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 Krylovtype 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 nonstiff, using an explicit method may be more efficient^{6}. Therefore, a hybrid integration strategy was designed based on a measure of stiffness. The ODE system is advanced in time according to
where the norm^{7} of the Jacobi matrix is defined by
The explicit solver DOPRI5(4) is the wellknown fifthorder accurate RungeKutta 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 timestep is introduced given by
It means that the followup reaction timestep 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 Strangsplitting procedure. The coupling parameters ϵ_{n} and ϵ_{T} are userprescribed 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 timestep and, hence, increasing computational costs.
Fourstage GRK4(3)T method.
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
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 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 thirdorder and less dissipative PPM algorithm produces somewhat higheramplitude variations in the region . Moreover, no artificial discontinuity is created in the components f_{2} and f_{3} 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 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 (f_{1} – solid line; f_{2} – dashed line; f_{3} – dashdotted 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 wellknown 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 for concentrations n_{1},n_{2},n_{3} with rate coefficients Initial concentrations, obtained concentrations with GRK4(3)T at the final integration time t = 360, and highprecision 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_{3} adopting numerical parameters RTOL= 10^{4} and ATOL= 10^{14}.
Concentrations in the OREGO problem.
Fig. 3
Numerical solution of the OREGO problem (with logarithmic scaling for n_{2}). 
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_{0}(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 secondorder operatorsplit framework with secondorder 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_{2} in the lower resolved run. In this case insufficient resolution also causes a dispersionlike error noticeable as a moderate leftward shift of the righthand feature. This dispersion error results from the interaction of numerical dissipation errors with the update of rate equations in virtue of the operatorsplit approximation. The more highly resolved case produces results much closer to the reference solution. The outcome of both low and highresolution runs is not sensitive to the choice of coupling parameter ϵ_{n} in the expression for Δt_{reac} because there is no backreaction 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 hydrogenhelium plasma with species H, H^{+}, He, He^{+}, He^{2+} in ionization equilibrium is considered. Assuming a fixed temperature of 10^{5} K and initial number densities far from equilibrium for this temperature, namely,
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^{6} 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+} = ṅ_{He2 +} = 0, the number densities are given by
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, builtup He^{+} is transformed into He^{2+} in a few thousand years. At t_{max}, equilibrium balance is reached with relative deviation less than 10^{8} from exact equilibrium values.
Number densities in the relaxation problem.
Fig. 5
Numerical solution of the ionization/recombination relaxation problem. Number densities of all species are shown as a function of time. 
4.5. Onedimensional 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
where is the initial mass fraction of the different species X (see Table 4).
Initial mass fractions of species in the Riemann problem.
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_{2}, H^{−}, and e^{−}. In order to allow a qualitative comparison^{8} with the results of Grassi et al. (2014) the simulation was stopped at canonical time t = 2.45 × 10^{4} yr (see their Figs. 19–22). Differences from the pure gasdynamical case without chemistry are mainly due to the chemical evolution of H_{2} and H_{2}related cooling mechanisms which are the dominant thermal processes. It leads to a significant temperature reduction in regions of relatively high H_{2} mass fraction. Finally, the overall characteristics of chemothermal evolution obtained with NIRVANA agree well with the behavior observed in Grassi et al. (2014)^{9}.
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 H_{2} (top), H^{−} (middle), and e^{−} (bottom) as a function of distance (unit pc) for the modified Sod problem with chemistry. 
4.6. Onedimensional propagating shock
As a second example combining magnetogasdynamics 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_{0} = 1000 K, zero velocity v_{0} = 0, and constant (perpendicular) magnetic field B_{y0} = 10^{8} T. The background density profile is given by
where ϱ_{00} = 10^{11}m_{u} m^{3} and x_{0} = 1.8 × 10^{12} m. The computational domain ranges from x = 0 to x = 4 × 10^{13} m assuming outflow boundary conditions on both sides and adopting an effective spatial resolution^{10} of δx = 1.25 × 10^{9} 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 nonmagnetic Riemann invariant
is fulfilled with γ assumed 5/3. At t = 0 this results in a structure (background plus perturbation) given by with computed for the background medium. The velocity perturbation v′ is prescribed as
where the peak value , σ = 2 × 10^{11} m and x_{1} = 10^{12} m. We note that exterior to the interval (x_{1},x_{1} + 2σ) the configuration reduces to the isothermal background.
The gas is composed of species X with initial mass fractions 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_{2}related cooling/heating mechanisms. In the second simulation H_{2} is treated as a passive scalar, H_{2} chemistry is switched off, and there is no H_{2} cooling/heating. It serves to isolate the effects of H_{2} chemistry which are expected to be important in virtue of the relatively high initial H_{2} mass fraction of and chosen temperature T_{0}. 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 postshock 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^{4} K, which partially ionizes the region with extent ≈3 × 10^{12} m right behind the shock front. Also, the magnetic field reaches its maximum in this region. The influence of H_{2} 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^{12} m and around the peak in molecular hydrogen density at x ≈ 1.3 × 10^{13} 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_{2} 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 preshock 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_{2} chemistry/cooling with the SNEq model in Teşileanu et al. (2008) further strengthens confidence in the validity of the presented algorithmic approach.
Initial mass fractions of species in the propagating shock problem.
Fig. 8
From top to bottom: gas density, electron number fraction, H_{2} mass fraction, temperature, velocity and magnetic field at t = 13 yr (solid line: full chemistry; dashed line: no H_{2} 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 advectionreaction equations for chemical species. These equations add to the already complex equations of magnetogasdynamics. In particular, the ODE system built from the rate equations representing chemical kinetics and the temporal temperature change due to chemistryconnected thermal processes is solved coherently with a highorder stiff integrator of linearimplicit 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 Krylovtype approximation in order to reduce costs further. Such possible improvements are the subject of future work.
Complementary to the numerical framework an easytouse interface was developed which allows straightforward upgrades of the HHe 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 timedependent chemomagnetogasdynamical simulations. In combination with the code’s selfgravity 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.
This is an incomplete maximum norm composed of the maximum norm of the chemical kinetics subsystem and an inverse temperature timescale due to thermal processes (J_{Ns + 1,Ns + 1}element). Use of the full maximum norm is not suitable here because matrix elements do not all have the same physical dimension.
In Grassi et al. (2014) their temperature Eq. (61) misses the term ∝ṅ_{tot}/n_{tot} which is, however, very small in this test problem.
In order to speed up the simulation an adaptive mesh with a 4000 equidistantcell base grid and three additional levels of refinement was used (32 000 effective grid). This is a slightly better resolution than in Teşileanu et al. (2008), but not as good as in Massaglia et al. (2005).
Acknowledgments
I would like to thank the referee, Pierre Lesaffre, for valuable comments useful in improving the paper.
References
 Abel, T., Anninos, P., Zhang, Y., & Norman, M. L. 1997, New Astron., 2, 181 [Google Scholar]
 Abel, T., Bryan, G. L., & Norman, M. L. 2002, Science, 295, 93 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Black, J. H. 1981, MNRAS, 197, 553 [NASA ADS] [Google Scholar]
 Bryan, G. L., Norman, M. L., O’Shea, B. W., et al. 2014, ApJS, 211, 19 [NASA ADS] [CrossRef] [Google Scholar]
 Cen, R. 1992, ApJS, 78, 341 [NASA ADS] [CrossRef] [Google Scholar]
 Dalgarno, A., & Lepp, S. 1987, Astrochemistry (Dordrecht: Reidel) [Google Scholar]
 Donahue, M., & Shull, J. M. 1991, ApJ, 383, 511 [NASA ADS] [CrossRef] [Google Scholar]
 Dormand, J. R., & Prince, P. J. 1980, J. Comput. Appl. Math., 6, 19 [CrossRef] [MathSciNet] [Google Scholar]
 Dove, J. E., & Mandy, M. E. 1986, ApJ, 311, L93 [NASA ADS] [CrossRef] [Google Scholar]
 Field, R. J., Körös, E., & Noyes, R. M. 1972, J. Am. Soc., 94, 8649 [CrossRef] [Google Scholar]
 Fryxell, B., Olson, K., Ricker, P., et al. 2000, ApJS, 131, 273 [NASA ADS] [CrossRef] [Google Scholar]
 Glover, S. C. O., & Abel, T. 2008, MNRAS, 388, 1627 [NASA ADS] [CrossRef] [Google Scholar]
 Glover, S. C. O., & Mac Low, M.M. 2007, ApJ, 659, 1317 [NASA ADS] [CrossRef] [Google Scholar]
 Grassi, T., Bovino, S., Schleicher, D. R. G., et al. 2014, MNRAS, 439, 2386 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Gray, C. 2002, RoseHulman Undergraduate Mathematics Journal, 3, 1 [Google Scholar]
 Hairer, E., & Wanner, G. 1996, Solving Ordinary Differential Equations II. Stiff and Differentialalgebraic Problems (Springer Verlag) [Google Scholar]
 Hindmarsch, A. C. 1983, IMACS Trans. Sci. Comput., 1, 55 [Google Scholar]
 Hollenbach, D., & McKee, C. F. 1979, ApJS, 41, 555 [NASA ADS] [CrossRef] [Google Scholar]
 Janev, R. K., Langer, W. D., & Evans, K. 1987, Elementary Processes in Hydrogen–Helium Plasmas (Berlin: Springer Verlag) [Google Scholar]
 Kaps, P., & Rentrop, P. 1978, Numer. Math., 33, 55 [CrossRef] [Google Scholar]
 Karpas, Z., Anicich, V., & Huntress, W. T. 1979, J. Chem. Phys., 70, 2877 [NASA ADS] [CrossRef] [Google Scholar]
 Lesaffre, P., Chièze, J.P., Cabrit, S., & Pineau des Forêts, G. 2004, A&A, 427, 147 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lesaffre, P., Gerin, M., & Hennebelle, P. 2007, A&A, 469, 949 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Le Petit, F., Nehme, C., Le Bourlot, J., & Roueff, E. 2006, ApJS, 164, 506 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lim, A. J., Rawlings, J. M. C., & Williams, D. A. 1999, MNRAS, 308, 1126 [NASA ADS] [CrossRef] [Google Scholar]
 Massaglia, S., Mignone, A., & Bodo, G. 2005, A&A, 442, 549 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Omukai, K. 2000, ApJ, 534, 809 [NASA ADS] [CrossRef] [Google Scholar]
 Plewa, T., & Müller, E. 1999, A&A, 342, 179 [NASA ADS] [Google Scholar]
 Poludnenko, A., Varniere, P., Frank, A., & Mitran, S. 2005, Springer’s Lecture Notes in Computational Science and Engineering, 41 [Google Scholar]
 Poulaert, G., Brouillard, F., Claeys, W., McGowan, J. W., & Van Wassenhove, G. 1978, J. Phys. B, 11, L671 [NASA ADS] [CrossRef] [Google Scholar]
 Ripamonti, E., & Abel, T. 2004, MNRAS, 348, 1019 [NASA ADS] [CrossRef] [Google Scholar]
 Semenov, D., Hersant, F., Wakelam, V., et al. 2010, A&A, 522, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Shapiro, P. R., & Kang, H. 1987, ApJ, 318, 32 [NASA ADS] [CrossRef] [Google Scholar]
 Stone, J. M., & Norman, M. L. 1992, ApJS, 80, 791 [Google Scholar]
 Strehmel, K., Weiner, R., & Podhaisky, H. 2012, Numerik gewöhnliche Differentialgleichungen (Wiesbaden: Vieweg+Teubner Verlag) [Google Scholar]
 Teşileanu, O., Mignone, A., & Massaglia, S. 2008, A&A, 488, 429 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Teyssier, R. 2002, A&A, 385, 337 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 van Leer, B. 1977, J. Comput. Phys., 23, 276 [Google Scholar]
 Wakelam, V., Herbst, E., Loison, J.C., et al. 2012, ApJS, 199, 21 [NASA ADS] [CrossRef] [Google Scholar]
 Ziegler, U. 2008, Comput. Phys. Common., 179, 227 [Google Scholar]
 Ziegler, U. 2011, J. Comput. Phys., 230, 1035 [Google Scholar]
 Ziegler, U. 2012, SIAM J. Sci. Comput., 34, C102 [CrossRef] [Google Scholar]
Appendix A: Chemical network
Table A.1 catalogs the chemical network of a hydrogenhelium plasma used in this work. The temperaturedependent rate coefficients for collisional ionization, radiative, and dielectronic 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.
Hydrogenhelium network.
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 hydrogenhelium plasma is adopted. It includes ionization cooling of H, He, and He^{+} by collisions with electrons; radiative recombination cooling of H^{+}, He^{+}, and He^{2+}; dielectronic recombination cooling of He^{+}; collisional ionization cooling of the metastable He(2^{3}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^{2+} according to Black (1981) is taken into account. We note that the effective cooling rates for He(2^{3}S) and He (n = 2,3,4 triplets) cooling are ∝, 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 rovibrational levels of molecular hydrogen is considered. The net energy loss per volume per time is written in the form
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)
The cooling function Λ_{H2} measured in units [J/s] is given by
which, in the highdensity limit (Λ_{H2,CE} ≫ Λ_{H2,LTE}), reduces to the LTE cooling function following the work of Hollenbach & McKee (1979). In the lowdensity 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_{2}with collision partners X_{k} ={H, H^{+}, He, H_{2},e^{−}}:
The logarithms of individual cooling rates are polynomial functions in log T, i.e., log Λ_{H2,k} = ∑ _{l}c_{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_{2}and para–H_{2}; 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
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
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_{2}which is then partly converted into heat via collisional deexcitation rather than radiated away. However, these heat sources only become important at high densities (n ≳ 10^{14} 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 (f_{1} – solid line; f_{2} – dashed line; f_{3} – dashdotted 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 n_{2}). 

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 H_{2} (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, H_{2} mass fraction, temperature, velocity and magnetic field at t = 13 yr (solid line: full chemistry; dashed line: no H_{2} chemistry; dotted line: no chemistry, adiabatic). The thin solid line is the initial profile. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.