© D. Przybylski et al. 2022

1 Introduction

The importance of the solar chromosphere for resolving a number of the large open questions in solar physics is undoubted. For example, it provides the connection between the solar surface, the source of the energy for the upper solar atmosphere, and the corona, where this energy is deposited, the local plasma is heated, and the solar wind is accelerated. A detailed understanding of the dynamics and structure of the solar chromosphere has consequently been a major goal of solar physics for decades. Energy transfer in the chromosphere is strongly affected by the interplay between radiation and the chromospheric plasma. An accurate treatment of radiation transfer (RT) is necessary to model the structure of the chromosphere.

The theoretical treatment of RT in the chromosphere is difficult because it cannot be treated as optically thin, unlike in the corona, nor can it be treated in local thermodynamic equilibrium (LTE), unlike in the photosphere. Additionally, the large recombination timescale of ionised hydrogen and helium means the problem cannot be treated in statistical equilibrium (SE). The dynamics and radiation effects must be solved together in non-equilibrium (NE; Carlsson & Stein 2002; Judge 2005). Additionally, the low ionisation fraction and low collisional frequencies may lead to a drift between the ionised and neutral components of the plasma. The weak coupling between ions and neutrals can lead to ambipolar diffusion and Hall drift becoming significant, or may even require a multi-fluid treatment.

A multi-dimensional treatment of the chromosphere is complicated by the non-locality of non-local thermodynamic equilibrium (NLTE) radiation transport. The important chromospheric spectral lines are strongly scattering and should in principle be treated with partial frequency redistribution (PRD). For NLTE RT simulations, the computational time scales proportionally to $n_z^2$, where n_z is the number of points in the vertical direction. The three-dimensional spectral synthesis of a single chromo-spheric spectral line, such as Ca H&K, Mg H&K, or hydrogen Lyman-Alpha, can cost 50–200 kcore-H per million grid points (Sukhorukov & Leenaarts 2017). The computational cost of detailed NLTE RT (23.23 s per grid point for hydrogen in PRD) compared to an update of the multi-group LTE radiation magnetohydrodynamic (rMHD) code (15.80 µs per grid point) leaves large multi-dimensional simulations, including detailed radiative NLTE transfer, out of reach. One-dimensional codes exist, for example the RADYN code (Carlsson & Stein 1992). In multidimensional rMHD simulations of the solar chromosphere, two NLTE effects are critical: the radiative cooling due to NLTE chromospheric spectral lines and the rate of ionisation (or recombination) of hydrogen and helium in the solar chromosphere. For multi-dimensional simulations, we turn to approximations of these effects in order to make three-dimensional simulations computationally tractable.

The three-dimensional computational modelling of the solar chromosphere has so far been led by the Bifrost group (Gudiksen et al. 2011). The Bifrost code includes a number of approximations to NLTE and NE physics, producing the most realistic three-dimensional simulations of the chromosphere currently available. The approximations include tabulated recipes for computationally efficient NLTE chromospheric line losses based on a detailed synthesis of the radiation field, including PRD effects (Carlsson & Leenaarts 2012), scattering multi-group RT (Skartlien 2000; Hayek et al. 2010), and an NE hydrogen equation of state (EoS; Leenaarts et al. 2007). To more accurately simulate the transition region, Bifrost includes a three-dimensional treatment of hydrogen Lyman lines and the addition of a computationally efficient helium model atom (Golding et al. 2016).

In this paper we introduce an updated version of the MURaM¹ code (Vögler et al. 2005; Rempel 2014, 2017) that includes prescriptions for NLTE and NE effects. The prescriptions used are the same as those described in the Bifrost code (Gudiksen et al. 2011) but do not include the Golding et al. (2016) extensions.

The MURaM code has been employed to treat many phenomena in the solar photosphere, such as umbral dots (Schüssler & Vögler 2006), sunspots (Rempel et al. 2009), small-scale dynamo (Vögler & Schüssler 2007), and magnetic flux emergence (Cheung et al. 2007; Chen et al. 2017), and has been used to investigate the effect of ambipolar diffusion in three-dimensional rMHD simulations (Cheung & Cameron 2012; Danilovic 2017). The code has also been used to compute molecular (Schüssler et al. 2003) and atomic (Shelyag et al. 2007) diagnostics, as well as for modelling of the solar irradiance variability (Shapiro et al. 2017; Yeo et al. 2017). Furthermore, it has successfully treated phenomena in the solar corona (Rempel 2017; Cheung et al. 2019) and in stellar photospheres (Beeck et al. 2013, 2015; Panja et al. 2020). With the extension described here, this versatile code will be able to cover the solar chromosphere and thereby close one of the large remaining gaps in its ability to treat phenomena in the solar atmosphere.

In Sect. 2, we describe the numerical methods used, including the diffusion scheme, the NE EoS, and radiative losses. In Sect. 3, we outline the experimental setup and present the first simulations of the solar chromosphere using this updated version of MURaM. In Sect. 4, we analyse the radiative cooling in the chromosphere, in Sect. 5 we discuss the effect of the NE treatment of hydrogen in the chromosphere. In Sect. 6, we discuss the performance impact of the new routines. Finally, in Sect. 7 we discuss the results and present our conclusions.

2 Numerical Approach

The MURaM code (Vögler et al. 2005) solves the conservative magnetohydrodynamic (MHD) equations on a Cartesian grid in one-, two- or three-dimensions. Spatial derivatives are calculated using a fourth-order central difference scheme. Temporal integration is performed with the Jameson-Schmidt-Turkel scheme (Jameson 2017), a four-stage explicit time-update scheme. The code's original hyper-diffusion scheme was replaced by a hybrid scheme based around slope-limiters and higher order hyperdiffusion scheme (Rempel 2014). Further enhancements were made by Rempel (2017) to allow simulations of the solar corona, including; optically thin losses, point-implicit heat conduction and a semi-relativistic ‘Boris correction’ to circumvent the time-step restrictions due to the high Alfvén velocity (Boris 1970). A pre-tabulated EoS using the Opal (Rogers et al. 1996) or Uppsala (Gudiksen et al. 2011) packages is used to calculate temperature, pressure, and electron number from the plasma density and internal energy.

In this section, we describe further extensions to realistically simulate the solar chromosphere, including modifications to the diffusion scheme, the EoS and the implementation of an NE treatment of hydrogen populations, and radiative cooling and heating.

We first summarise the procedure used to update the system of equations over a time-step, ∆t. First, the radiative cooling and heating (Sect. 2.5) are calculated using the system state at the previous time, t₀. In each sub-stage, n, the steps followed are: calculate the right-hand side (RHS) of the MHD equations using a directionally un-split approach (Sect. 2.6), integrate in time to the new state t_⋆ = t₀ + 1/(5 − n)∆t, then apply the boundary conditions, advect the populations, and evaluate the EoS (Sect. 2.2). No hyper-diffusion or explicit diffusivities are included in the code, instead once all the sub-stages have been completed the directionally split diffusion scheme (Sect. 2.1) is run. After each directional sweep the EoS equations must be solved. Next the hyperbolic ∇·B cleaner (Dedner et al. 2002) is applied. After the MHD variables have been updated to the next time step, the EoS is solved once more, including a set of hydrogen rate equations (Sect. 2.3).

2.1 Numerical Diffusion Scheme

The MURaM code includes a hybrid diffusion scheme that is based around slope-limiters and includes higher order hyperdiffusion terms (Rempel 2014, 2017). For the simulations presented in this work, we introduce an additional scheme based around the partial donor cell method (PDM) as described in detail in Zhang et al. (2019). The PDM limiter is used to avoid the undershoot or overshoot that occurs near discontinuities when a high-order (HO) scheme is used.

To apply the limiter to a quantity Φ, the flux through the interface i + 1/2 between cell i and i + 1 must be calculated. We used a fourth-order centred reconstruction to calculate the HO value at the cell interface ${\Phi }_{i+1/2}^{\text{HO}}$. When a discontinuity is detected, the PDM limiter decides whether the value at the cell interfaces needs to be ‘limited’. The left (l) and right (r) interface values are calculated: $$\begin{array}{ll}{\text{Φ}}_{i+1/2}^l=\hfill & {\text{Φ}}_{i+1/2}^{HO}-\text{sign\hspace{0.17em}}\left({\nabla }^{+}\right)\max \left[0,\left|{\text{Φ}}_{i+1/2}^{HO}-{\text{Φ}}_i\right|\right]\hfill \\ \hfill & -C_{\text{pdm}}\left({\nabla }^{-}{\nabla }^{+}-\left|{\nabla }^{-}{\nabla }^{+}\right|\right)\left|{\nabla }^{-}\right|,\text{\hspace{0.17em}and}\hfill \end{array}$$(1) $$\begin{array}{ll}{\text{Φ}}_{i+1/2}^r=\hfill & {\text{Φ}}_{i+1/2}^{HO}+\text{sign\hspace{0.17em}}\left({\nabla }^{+}\right)\max \left[0,\left|{\text{Φ}}_{i+1}-{\text{Φ}}_{i+1/2}^{HO}\right|\right]\hfill \\ \hfill & -C_{\text{pdm}}\left({\nabla }^{++}{\nabla }^{+}-\left|{\nabla }^{++}{\nabla }^{+}\right|\right)\left|{\nabla }^{++}\right|,\hfill \end{array}$$(2)

where ∇⁺ = Φ_i+1 − Φ_i, ∇⁻ = Φ_i − Φ_i−1, ∇⁺⁺ = Φ_i+2 − Φ_i+1. C_pdm is a parameter that controls the amount of diffusion, with C_pdm = 0 equivalent to a first order donor cell scheme, and C_pdm > 0 corresponding to a lower diffusivity. The diffusive flux f_i+1/2 across the cell interface is calculated from the limited values, $$\begin{array}{rr}\hfill f_{i+1/2}=& \hfill -\frac{c_{i+1/2}}{2}|\text{sign\hspace{0.17em}}\left({\text{Φ}}_{i+1/2}^r-{\text{Φ}}_{i+1/2}^l\right)\\ \hfill & \hfill +\text{sign\hspace{0.17em}}\left({\nabla }^{+}\right)|\left({\text{Φ}}_{i+1/2}^r-{\text{Φ}}_{i+1/2}^l\right),\end{array}$$(3)

where C_i+1/2 = C_s + υ_A + |v| is the characteristic velocity in terms of the sound speed, C_s, Alfvén velocity, υ_A, and velocity vector, v. Once the fluxes are calculated, the remainder of the scheme is the same as that described in Rempel (2017). The diffusion scheme was applied to the logarithm of the density ρ, internal energy per gram ϵ, and population fractions n/ρ, as well as the velocity components υ_x, υ_y, υ_z, and magnetic field vector B_x, B_y, B_z. Applying the diffusion to the logarithm in a stratified atmosphere reduces the systematic vertical diffusive fluxes. The energy, momentum and population numbers were corrected for mass diffusion.

To increase stability and minimise diffusion, a number of enhancements were made, following those in Rempel (2017). In order to remove wiggles superposed on the stratified atmosphere, fourth-order hyper diffusion was added in the vertical direction to log ρ, log ϵ, log n, and υ_z.

Additionally, a number of switches exist to allow the code to run at a less diffusive setting while taking care of the few localised grid points that require higher diffusion. One can operate the scheme without these switches, but a more diffusive setting would be required everywhere in order to keep the code stable. The default value of diffusion coefficient used in the simulation was C_pdm = 2. These switches include; To reduce the errors in ∇·B the diffusive flux of B_z in the z-direction was set to zero at the vertical boundaries and the numerical diffusivity of B in the direction of B was reduced by a factor of 0.2. Secondly, to reduce the diffusion in the convection zone the sound-speed contribution to the characteristic velocity C_i+1/2 was limited to a maximum of 1 km s⁻¹.

Two hard switches are included to prevent the formation of instabilities that can occur in the solar atmosphere. Firstly, if the maximum density contrast between grid points exceeds 10, then the mass diffusivity is increased (C_pdm = 0). This setting ensures stability of the code, providing extra diffusion at sharp shock-fronts and other extreme events. Secondly, diffusivity in regions with a low adiabatic index (γ < 1.12) is increased (C_pdm = 0.5). The increased diffusivity at low gamma was included due to an instability in the chromosphere. When a strong flow is present, and the plasma is around the ionisation temperature of hydrogen (the dominant species), recombination of protons to neutral hydrogen can cause a sharp drop in pressure. The pressure gradient, combined with the diffusion scheme, was found to drive an instability at the grid scale. The diffusion scheme used is not a smooth Laplacian, the numerical diffusivities are highly intermittent. The addition of these two extra switches does not greatly change the behaviour of the diffusion scheme, but allows for the use of an overall lower diffusivity in the numerical domain.

The MURaM code, like many solar and stellar rMHD codes, uses non-isotropic grid spacing in the horizontal and vertical direction, many codes (e.g. Bifrost; Gudiksen et al. 2011) additionally use a non-uniform spacing in the vertical direction. The ratio of horizontal to vertical grid spacing is typically around 1.5–3. This ratio will lead to anisotropies in the numerical diffusion. In simulations of the convection zone and atmosphere, the structures modelled are highly anisotropic. Even the use of an isotropic numerical grid will not give isotropic diffusivities. A higher vertical resolution is suitable for convection simulations, as the vertical-to-horizontal ratio of convective cells is approximately 1/3. Although this argument breaks down in the atmosphere, the increased vertical resolution is also advantageous for radiative transport. The grid anisotropy allows a better resolution of sharp vertical gradients, such as those present at the photosphere and transition region. Simulations with the MURaM code and a variety of different resolutions have found no significant effects due to the grid anisotropy (Vögler & Schüssler 2007; Rempel 2014). Comparison of these simulations and others with non-uniform meshes (Beeck et al. 2012) have found no significant systematic differences with using the anisotropic grid. An additional systematic effect of the diffusivity can exist due to a tendency for vertical diffusive fluxes in the stratified atmosphere. The systematic effects in the horizontally averaged diffusive fluxes of mass, energy and vertical momentum remain small at the photosphere. We find that the horizontally averaged hydrostatic balance is conserved to within a percent in the convection zone and lower atmosphere. In the corona, large cross-field gradients of thermodynamic quantities can exist, which can lead to enhanced diffusive fluxes.

2.2 Non-Equilibrium Equation of State

To perform self-consistent simulations of the solar atmosphere, we required an EoS. We combined two approaches; firstly, we used a pre-tabulated LTE EoS to model the interior based on the prescription introduced by Vardya (1965) and extended by Mihalas (1967) and Wittmann (1974) (the VMW prescription; see also Vitas & Khomenko 2015). Secondly, an EoS that includes an NE treatment of hydrogen was used for the solar chromosphere and corona. The pre-tabulated LTE EoS was used for pressures greater than 2 × 10⁵ dyncm⁻², and an EoS state that includes the NE treatment of hydrogen was used for lower pressures. As the non-ideal gas and NE ionisation effects are negligible at this pressure, the two methods join smoothly. We used a mixture of the 15 most abundant elements of the Sun and included the molecules H₂ in NE, and ${\text{H}}_2^{+}$ and H^– in chemical equilibrium. We included up to three ionisation states for all LTE elements and a five-level plus continuum model of hydrogen in NE. In this section we focus on the NE EoS, and an overview of the LTE EoS is given in Appendix A.

At each sub-stage of the iteration scheme, or after each directional-sweep of the diffusion scheme, the MHD solver provides updated values of the internal energy density E_int and density ρ of the plasma. From the density and the atomic abundances, the total number density of hydrogen nuclei n_H,tot and non-hydrogen nuclei n_nonH,tOt are calculated. We then find a solution to the equations of energy conservation, charge conservation and nuclei conservation, in terms of temperature (T), electron number density (n_e) and population number densities (n = {n_a,i,j}). We use the notation (a, i, j) to represent a species (atom or molecule) a of ionisation stage i and energy level j.

Two tables are required for the NE EoS to include the contribution from non-hydrogen atoms. The thermodynamics of non-hydrogen atoms are treated in LTE. The electron number density per hydrogen nuclei (n_e,nonH) and energies of excitation and ionisation per hydrogen nuclei (ϵ_nonH) for non-hydrogen atoms are tabulated as a function of temperature and electron number. These are calculated as $$n_{\text{e,nonH}}\left(n_{\text{e}},T\right)={\displaystyle \sum _{a=1}^{15}{\displaystyle \sum _{i=0}^{3}i\frac{n_{a,i}}{n_{a\text{,tot}}}\frac{n_{a,\text{tot}}}{n_{\text{H,tot}}},}}$$(4) $${є}_{\text{nonH}}\left(n_{\text{e}},T\right)={\displaystyle \sum _{a=1}^{15}{\displaystyle \sum _{i=0}^3\left[{\chi }_{a,i}+k_{\text{B}}{T}^2\frac{\partial \ln \left(U_{a,i}\right)}{\partial T}\right]}}\frac{n_{a,i}}{n_{a\text{,tot}}}\frac{n_{a\text{,tot}}}{n_{\text{H,tot}}},$$(5)

where n_a,tot/n_H,tot is the number fraction of element a relative to hydrogen,χ_a,i is the ionisation energy, k_B the Boltzmann constant, U is the partition function and n_a,i/n_a,tot is the fraction of element a in ionisation stage i calculated using Saha–Boltzmann (Eq. (A.1)). The partition functions used are described in Appendix A.

Following Leenaarts et al. (2007) the NE EoS is evaluated by solving a system of equations using a Newton–Raphson method. The equation of energy conservation is $$\begin{array}{ll}{f}_0=\hfill & 1-\frac{1}{E_{\mathrm{int}}}(\left(\frac{3k_{\text{B}}T}{2}\right)[n_{\text{e}}+n_{\text{nonH}}+n_{{\text{H}}_2}+n_{{\text{H}}_2,1}+n_{\text{H}}{,}_{-1}\hfill \\ \hfill & +{\displaystyle \sum _{i,j}{n}_{\text{H,}i,j}}]+n_{\text{H,TOT}}{є}_{\text{nonH}}+n_{{\text{H}}_2}+n_{{\text{H}}_2,}{}_1{E}_{{\text{H}}_2,1}\hfill \\ \hfill & +n_{\text{H,}-1}{E}_{\text{H,}-1}+{\displaystyle \sum _{i,j}{}_{\text{H,}i,j}})=0,\hfill \end{array}$$(6)

where the E_a,i,j gives the ionisation, excitation and dissociation energies of the atom/molecule. The derivatives of n_e,nonH and ϵ_nonH with respect to n_e and T were calculated numerically from the table interpolants. The energies of the hydrogen species are given by Eqs. (A.4)–(A.7). The equation of charge conservation is $$f_1=1-\frac{1}{n_{\text{e}}}\left(n_{\text{H,1}}+n_{{\text{H}}_2\text{,1}}-n_{\text{H,tot}}{n}_{\text{e,}}{}_{\text{nonH}}\right)=0.$$(7)

Finally, nucleus conservation must be maintained. All non-hydrogen elements are considered in LTE. The equation of hydrogen nucleus conservation is $$f_2=1-\frac{1}{n_{\text{H,tot}}}\left(n_{\text{H,}i,j}+2n_{{\text{H}}_2}+n_{{\text{H}}_2,1}+n_{\text{H,-1}}\right)=0.$$(8)

A full set of the derivatives are described in Appendix B.

2.3 Non-Equilibrium Hydrogen Populations

The time evolution of a species n_a,i,j depends on the advection of the populations with the bulk fluid velocity v, the rate of collisional (C_ij,kl) and radiative (R_ij,kl) transitions between level (ij) and level (kl ≠ ij), and the rate of molecule formation or destruction by gas-phase reactions (r). The gas-phase reactions are described in terms of a set of reactants (A, B, C) and rate coefficient, K_r, which form (r⁺) or destroy (r⁻) the species (aij). Defining P_ij,kl = C_ij,kl + R_ij,kl, the rate equation of a species n_aij is described by $$\begin{array}{ll}\frac{\partial n_{ai}{}_j}{\partial t}=\hfill & -\nabla \cdot \left(n_{aij}\text{v}\right)\hfill \\ \hfill & +{\displaystyle \sum _{kl\ne ij}{n}_{ak}{P}_{kl,ij}-{\displaystyle \sum _{kl\ne ij}{P}_{ij,kl}}}\hfill \\ \hfill & +{\displaystyle \sum _{r^{+}}{n}_{r^{+},A^n{r}^{+},B^n{r}^{-}{,}_C{K}_{r^{+}}}}\hfill \\ \hfill & -{\displaystyle \sum _{r^{-}}{n}_{r^{-},A^n{r}^{-},B^n{r}^{-}{,}_C{K}_{r^{-}}.}}\hfill \end{array}$$(9)

In this work, we treated only hydrogen and H2 in NE. To solve Eq. (9), we separated it into continuity and rate components. The advection of the species by the macroscopic fluid velocity v is then given by $$\frac{\partial n_{aij}}{\partial t}=-\nabla \cdot \left(n_{aij}\text{v}\right).$$(10)

We found that advecting the populations with the fourth-order central differences and four-stage temporal integration scheme used for MHD variables leads to frequent negative values. Instead, we used a flux-limited un-split donor cell method. Each sub-stage of the temporal integration scheme advances the fluid variables ρ, ϵ, v in time, with a time step ∆t, from t₀ to t_* = t₀ + ∆t. To advect the populations we used the average velocity v^1/2 = 0.5(v^t0 + v^t*). For each direction, the velocities v^1/2 were interpolated to the cell interfaces using quadratic Bezier interpolation (de la Cruz Rodríguez & Piskunov 2013). The limited left(l)- and right(r)-interface values were calculated using the PDM limiter (Eqs. (1) and (2)), as applied in the diffusion scheme (Sect. 2.1). The flux of the quantity Φ through the cell interface i + 1/2, in direction d = x,y, z, is calculated as $$\begin{array}{ll}{f}_{d,i+1/2}=\hfill & \frac{1}{2}\left({\upsilon }_{d,i+1/2}^{1/2}+\left|{\upsilon }_{d,i+1/2}^{1/2}\right|\right){\text{Φ}}_{i+1/2}^l\hfill \\ \hfill & +\frac{1}{2}\left({\upsilon }_{d,i+1/2}^{1/2}-\left|{\upsilon }_{d,i+1/2}^{1/2}\right|\right){\text{Φ}}_{i+1/2}^r,\hfill \end{array}$$(11)

and the populations are updated using the sum of the fluxes through all faces of the cell, $${\text{Φ}}^{t_{*}}={\text{Φ}}^{t_0}+{\displaystyle \sum _{d=x,y,z}\frac{\text{Δ}t}{\text{Δ}d}\left(f_{d,i-1/2}-f_{d,i+1/2}\right)}.$$(12)

Once the advected populations $n_{a,i,j}^{*}$ have been calculated, an update of the EoS Eqs. (6)–(8) ensures consistency between the EoS variables, populations, and MHD energy $E_{\mathrm{int}}^{*}$ and density ρ^*. After all sub-stages are complete the directionally split diffusion scheme is run. The populations are diffused using the scheme described in Sect. 2.1. The populations are also corrected for any mass diffusion.

In order to evaluate the system of hydrogen rates concurrently with the EoS, the rate equations are written in a form suitable for solution with the Newton–Raphson method, $$\begin{array}{ll}{f}_{3+i,j}=\hfill & \frac{n_{aij}^{t_0+\text{Δ}t}}{n_{aij}^{t_0}}-\frac{\text{Δ}t}{n_{aij}^{t_0}}({\displaystyle \sum _{kl\ne ij}{n}_{akl}{P}_{kl,ij}-n_{aij}}{\displaystyle \sum _{kl\ne ij}{P}_{ij,kl}}\hfill \\ \hfill & +{\displaystyle \sum _{r^{+}}{n}_{r^{+},A}{n}_{r^{+},B}{n}_{r^{+},C}{K}_{r^{+}}}\hfill \\ \hfill & -{\displaystyle \sum _{r^{-}}{n}_{r^{-},A}{n}_{r^{-},B}{n}_{r^{-},C}{K}_{r^{-}}})-1=0,\hfill \end{array}$$(13)

where we included the ground level, four excited states and the continuum. The radiative and collisional rates used for atomic hydrogen follow the method of Sollum (1999), and are described in Appendix C. The rate equation describing molecular hydrogen is $$\begin{array}{ll}{f}_8=\hfill & \frac{n_{{\text{H}}_2}^{t_0+\text{Δ}t}}{n_{{\text{H}}_2}^{t_0}}-\frac{\text{Δ}t}{n_{{\text{H}}_2}^{t_0}}({\displaystyle \sum _{r^{+}}{n}_{r^{+},A}{n}_{r^{+}{,}_C{K}_r{}^{+}}}\hfill \\ \hfill & -{\displaystyle \sum _{r^{-}}{n}_{r^{-},A}{n}_{r^{-},B}{n}_{r^{-},C}{K}_{r^{-}}})-1=0,\hfill \end{array}$$(14)

where the rates for the formation and destruction of molecular hydrogen are listed in Appendix D. We used five out of six of the atomic hydrogen rate equations and the rate equation of molecular hydrogen. We included the nucleus conservation equation, and discarded the rate equation of the level with the highest population. In principle this method can be extended to include an arbitrary choice of atoms and molecules.

In the wake of strong chromospheric shocks, the internal energy density of the plasma can become very low. In NE simulations, recombination is too slow for ionisation energy to be released as heat. Very low temperatures may lead to the EoS and opacity tables becoming inaccurate, or to poor convergence of the solver when the H₂ fraction becomes dominant. It is therefore necessary to include additional mechanisms to prevent temperatures from becoming too low. We included three mechanisms. Firstly, an additional time-step constraint was included. The time step ∆t was limited such that |Q_rad|∆t/E_int ≤ 0.25, where Q_rad is the total radiative cooling or heating. This damps large decreases in energy due to radiative cooling. Secondly, a minimum temperature threshold was set. Rather than including a parameterised heating term, as in Leenaarts et al. (2011), a temperature floor was implemented. For a minimum temperature T_min, a minimum value E_min for the instantaneous parts of the internal energy equation is calculated. This includes the kinetic nk_BT_min term, an increase in the ionisation and excitation energies of non-hydrogen atoms and the dissociation of molecules treated in chemical equilibrium. The internal energy was then limited to this minimum value after a call to the EoS routine. In order to match the simulation of Carlsson et al. (2016) we use T_min = 2500 K. Finally, if the populations would not converge to the required tolerance, due to the temperature dropping below a threshold value of 1000 K during a solver call, we allowed a small amount of H⁺ to be recombined to ensure convergence. This allows a faster convergence of grid points, which will be limited anyway by the 2500 K floor.

2.4 Scattering Multi-Group Radiation Transfer

The MURaM code uses a multi-group method (Nordlund 1982) to accurately and efficiently compute the frequency-dependent photospheric radiation field (Vögler et al. 2004). To accurately simulate the low chromosphere, the treatment of radiation was extended to include a scattering term. We followed the prescription of Skartlien (2000) and the extension of Hayek et al. (2010) to short characteristics (see also Collet et al. 2011). Calculation of the radiation field requires a solution of the time-independent radiative transfer equation $$\frac{\text{d}{I}_v}{\text{d}{\tau }_v}=S_v-I_v,$$(15)

where v is the frequency, S _v the source function, I_v the specific intensity, and dτ_v = χ_vds the optical thickness over a path length ds, where χ_v is the plasma opacity. Equation (15) must be solved for a number of ray directions $\widehat{n}$. The source function has been expanded to include a scattering term, $$S_v=\left(1-{\text{є}}_v\right)J_v+єB_v,$$(16)

where ϵ_v is the photon destruction probability, B_v the LTE (Planck) source function, and the mean intensity J_v is calculated as the integral over all angles, $$J_v=\frac{1}{4\pi }{\displaystyle {\int }_{4\pi }{I}_v\left(\widehat{n}\right)\text{\hspace{0.17em}d}\omega .}$$(17)

The type A quadrature of Carlson (1963), including three points per quadrant, was used to perform the angular integration. The radiative energy flux, F_v, is $${\text{F}}_v={\displaystyle {\int }_{4\pi }{I}_v\left(\widehat{n}\right)\widehat{n}d\omega .}$$(18)

For the formal solution to the radiative transfer Eq. (15) we use a short-characteristics scheme with linear interpolation. The intensity at a given point (O) is calculated as $$I_O=I_U{e}^{-\delta {\tau }_U}+{\text{Φ}}_U{S}_U+{\text{Φ}}_O{S}_O,$$(19)

where (U) is the upwind point, τ_U is the optical distance on the segment UO and the Φ quantities are the weights. The formal solution, Eqs. (16) and (19), can be written as $$J_v=\Lambda \left[S_v\right]+{\mathcal{J}}_v,$$(20)

where ${\Im }_v$ represents the transmitted contribution to J_v due to the given incident radiation at the boundaries of the computational domain, and Λ is the angle-averaged Lambda operator.

Direct solution of Eq. (20) is expensive, and simply updating Eq. (16) with the new J_v leads to slow convergence. We employed the approximate lambda iteration method (Cannon 1973) and the diagonal operator as the approximate operator. The scheme was iterated until a tolerance of 10⁻⁴ is reached on the relative correction of the source function. Once the source function has converged, the radiative cooling (or heating) Q_v is calculated from two equivalent expressions, $$Q_{\text{RT}}=-{\displaystyle {\int }_v\left(\nabla \cdot F_v\right)}\text{\hspace{0.17em}d}v=4\pi {\displaystyle \underset{v}{\int }{\chi }_v\left(J_v-S_v\right)\text{\hspace{0.17em}d}v.}$$(21)

The multi-group scheme is used to simplify the frequency spectrum into a number of subsets j, known as bands. Instead of detailed calculations incorporating 10³–10⁵ frequency points, many of the important processes for a radiative MHD magneto-convection simulation, such as line blanketing, can be captured by using as few as four to five bands (Vögler et al. 2004). To calculate the radiation field, we require three band-integrated quantities, the extinction coefficient, χj, the mean scattering albedo, (1 – ϵ)_j, and the band-integrated emissivity, (ϵB)_j. The binning process is discussed in detail in Appendix E. The integral over wavelength in Eq. (21) is then solved as $$\begin{array}{ll}{Q}_{\text{RT}}=\hfill & -{\displaystyle \sum _j[\left(\nabla \cdot {\text{F}}_{\text{j}}\right)\left(1-e^{-{\tau }_j/{\tau }_0}\right)}\hfill \\ \hfill & +4\pi {\chi }_j\left(J_j-S_j\right){\text{e}}^{-{\tau }_j/{\tau }_0}],\hfill \end{array}$$(22)

where τ_j is the band mean optical depth, τ₀ = 0.1 and the term e^–τj/τ0 provides a transition between the radiative energy and flux-divergence form of the equation. This prevents numerical round-off errors in the optically thick regime where J ≈ S, which are amplified as χ grows exponentially with depth (Bruls et al. 1999).

2.5 Radiative Cooling and Heating

To perform rMHD simulations from the convection zone to the corona, a radiation scheme that can accurately model a range of physical regimes is required. These include the deep interior, which can be treated using the diffusion approximation. The photosphere where spectral line formation becomes important, requiring a three-dimensional multi-group approach (Sect. 2.4). Finally, in the chromosphere the radiation transport must include NLTE effects, such as PRD and scattering. Detailed NLTE RT is too computationally expensive to be performed in a three-dimensional time-dependent simulation. To include accurate and fast radiative cooling and heating of the chromosphere and corona, we used the pre-tabulated radiative losses calculated by Carlsson & Leenaarts (2012).

The prescription of Carlsson & Leenaarts (2012) consists of two parts; NLTE line losses in the chromosphere from hydrogen, calcium and magnesium, and optically thin coronal losses. We used the overlap interval approach of Rempel (2017). In the original implementation the overlap interval is calculated only in the vertical direction. In this work, we took the isotropic average of the overlap interval in all three directions to better model the cooling around the irregularly shaped transition region.

In the chromosphere the most important spectral lines and continua are the Lyman-α, H-α, and the Lyman continuum of hydrogen and the Mg ii H & K and Ca ii H & K lines. These lines are modelled using a simplified description of the cooling and heating, for element a in ionisation stage i, $$Q_X=-f\left(T\right)q_{a,i}\left(\tau \right)n_{a,\text{tot}}{n}_{\text{e}},$$(23)

where f(T) = L_a,i(T)F_a,i(T). Here three pre-tabulated quantities² are used; L_a,i(T) is the optically thin radiative loss function, F_a,i(T) is the fraction of element a in ionisation stage i, and q_a,i(τ) is the escape probability as a function of some optical depth proxy τ. For the escape probability of calcium and magnesium the column mass is used as a proxy for τ, and for hydrogen the neutral hydrogen column density is used. The optically thin radiative losses Q_thin are calculated as $$Q_{\text{thin}}=-\text{Λ}\left(T\right)n_{\text{H,tot}}{n}_{\text{e}},$$(24)

where Λ(T) is given as a table in terms of temperature.

Additionally, we included the back-heating of chromospheric plasma (Q_back). This was performed using the three-dimensional radiation transport scheme, following Carlsson & Leenaarts (2012). The emissivity is given in terms of the optically thin coronal losses given by $${\eta }_{\text{back}}=-\frac{Q_{\text{thin}}}{4\pi },$$(25)

and the opacity at the ionisation edge of helium is used, $${\chi }_{\text{back}}=\alpha \frac{n_{\text{He\hspace{0.17em}I}}}{n_{\text{He}}}\left(T,p_e\right)\frac{n_{\text{He}}}{\rho },$$(26)

where α is the opacity at the ionisation edge of helium, n_{He I}/n_He is the neutral helium fraction and is pre-tabulated in LTE in terms of electron pressure and temperature, and n_He/ρ is the number of particles of helium per gram of stellar material.

The full radiative cooling/heating prescription is then $$Q_{\text{rad}}=Q_{\text{RT}}+Q_{\text{H}}+Q_{\text{Mg}}+Q_{\text{Ca}}+Q_{\text{thin}}+Q_{\text{back}},$$(27)

where Q_RT is the heating/cooling from the multi-group radiative transport scheme described in Sect. 2.4 above. To prevent over-cooling from both the LTE and NLTE losses in the upper chromosphere, the cooling in each radiation band was switched off based on the band-averaged optical depth τ. This was performed using a function of the form ${\tau }^2/\left({\tau }^2+{\tau }_{\text{cutoff}}^2\right)$. A value of τ_cutoff = 1.0 × 10^–4 was used in these simulations. The multi-group RT scheme groups frequencies into bands. Which band a particular frequency v goes into depends on the height where τ_v = 1. We used as a reference τ₅₀₀, the optical depth at 500 nm. Previous work has found that four bands are sufficient to capture back-heating and line-blanketing in the photosphere and temperature minimum (Vögler et al. 2005). We used a 4-band setup similar to Carlsson et al. (2016), with boundaries at the heights where τ₅₀₀ = 10^−1/2, 10^–3/2, and 10^−5/2.

2.6 MHD Equations

The set of equations solved by the MURaM code is $$\frac{\partial \rho }{\partial t}=-\nabla \cdot \left(\rho \text{v}\right)$$(28) $$\frac{\partial \rho \text{v}}{\partial t}=-\nabla \cdot \left(\rho \text{vv}\right)-{\nabla }_p+\rho g+F_{\text{L}}+F_{\text{SR}}$$(29) $$\begin{array}{ll}\frac{\partial E_{\text{HD}}}{\partial t}=\hfill & -\nabla \cdot \left[v\left(E_{\text{HD}}+p\right)+q\widehat{\text{b}}\right]+\rho v\cdot \text{g+}v\cdot F_{\text{L}}\hfill \\ \hfill & +v\cdot F_{\text{SR}}+Q_{\text{rad}}+Q_{\text{res}}\hfill \end{array}$$(30) $$\frac{\partial B}{\partial t}=\nabla \times \left(v\times B\right),$$(31)

where ρ is the plasma density, v the velocity vector, p the gas pressure, g the gravitational acceleration, F_L the Lorentz force, F_SR the semi-relativistic (Boris) correction, E_HD the hydrodynamic energy, q the Spitzer heat flux, $\widehat{b}=B/\left|B\right|$ the unit vector in the direction of the magnetic field vector B, Q_rad the radiative cooling/heating, Q_res the resistive heating, and T the gas temperature. Additional diffusive terms are applied to each equation, based on the scheme described in Sect. 2.1.

The various radiative cooling and heating effects are given by Q_rad and described in detail in Sect. 2.5 above. The hydrodynamic energy $E_{\text{HD}}=E_{\mathrm{int}}+\frac{1}{2}\rho v^2$ is used instead of the total energy. This prevents numerical errors when calculating E_int from E_tot in low-β regions where the magnetic energy dominates the total energy. To conserve the total energy the heating from the diffusion scheme is then added as Q_res. The Spitzer heat flux, q, is solved using the hyperbolic method (see Rempel 2017 for a full derivation), $$\frac{\partial q}{\partial t}=\frac{1}{\tau }\left(-f_{\text{sat}}\sigma T^{\frac{5}{2}}\left(\widehat{\text{b}}\cdot \nabla \right)T-q\right),$$(32)

where σ = 1.1 × 10⁻⁶ erg cm⁻¹ s⁻¹ K^−7/2 is the constant of Spitzer heat conductivity, f_sat controls the saturation of thermal conduction and τ, which is used to control the transition between parabolic and hyperbolic solutions to the heat conduction equation, has the form $$\tau ={\left(f_{\text{CFL}}\frac{\text{Δ}{x}_{\min }}{\text{Δ}t}-\left|\upsilon \right|\right)}^{-2}\frac{f_{\text{sat}}\sigma T^{7/2}}{E{}_{\mathrm{int}}},$$(33)

where $f_{\text{CFL}}\frac{\Delta x_{\min }}{\Delta t}-\left|\upsilon \right|$ is used as a maximum propagation speed in order to avoid violations of the CourantFriedrichs–Lewy (CFL) condition. The first term is chosen so that the maximum wave speed of the hyperbolic heat conduction is comparable to the maximum MHD wave speed with the limited Alfvén velocity $\left(c=f_{\text{CFL}}\frac{\Delta x_{\min }}{\Delta t}\right)$, in terms of the minimum spatial grid scale, x_min, and time step, Δt. In order to explicitly integrate the system of equations, we set a lower limit on τ of 4Δt.

The treatment of the Lorentz force F_L follows that of Rempel (2017), $$\begin{array}{ll}{F}_{\text{L}}=\hfill & \frac{f_A}{4\pi }\nabla \cdot \left(BB-\frac{1}{2}\text{I}B\right)\hfill \\ \hfill & +\frac{1-f_A}{4\pi }\left(\nabla \times B\right)\times B,\hfill \end{array}$$(34)

where $f_A=\frac{1}{\sqrt{1+{\left({\upsilon }_a/c_{\max }\right)}^4}}$ is the Alfvén limit factor in terms of the reduced speed of light C_max and Alfvén speed υ_a. This form is used to have a sharper transition between the limited and non-limited regime, the semi-relativistic form would give 1/(1 + (υ_a/C_max)²). The reduction of the Alfvén velocity is achieved through a semi-relativistic treatment with reduced speed of light (Boris correction), which can be implemented through a projection operator in the momentum equation (Gombosi et al. 2002; Rempel 2017) by adding the force term F_SR given by $$\begin{array}{ll}{F}_{\text{SR}}=\hfill & -\left(1-f_A\right)\left[ℐ-\widehat{b}\widehat{b}\right]\hfill \\ \hfill & \left(-\rho \left(\text{v}\cdot \nabla \right)\text{\hspace{0.17em}v}-\nabla p+\rho \text{g}+F_L+\nabla \cdot {\tau }^{\text{diff}}\right).\hfill \end{array}$$(35)

The components of the numerical viscous stress tensor (τ^dlff) are calculated as $${\left(\nabla \cdot {\tau }^{\text{diff}}\right)}_i={\displaystyle \sum _{j=x,y,z}\frac{1}{{\text{Δ}}_j}\left(f_{{\upsilon }_i}^{j+1/2}-f_{{\upsilon }_i}^{j-1/2}\right)},$$(36)

where Δ_j is the grid resolution, and $f_{v_i}^j$ are the diffusive fluxes of velocity component i in the jdirection, calculated using Eq. (3).

3 Simulation Setup

We present results of a simulation continued from the publicly available Bifrost simulation (Carlsson et al. 2016). The initial condition consists of a bipolar magnetic field region modelling an enhanced network magnetic field. The simulation was run in NE for 3850 s, starting from the first publicly available snapshot, ‘285’.

One difference between the two simulations comes from the use of a stretched grid in the Bifrost code. This must be interpolated onto a constantly spaced vertical grid suitable for the MURaM code. This was performed using log-linear interpolations for density, energy, electron and population numbers, and linear interpolations for velocities and magnetic fields.

A different EoS is used by the Bifrost and MURaM codes. The Bifrost code uses a non-ideal EoS based on the free-energy minimisation method (Gustafsson et al. 1975). We calculated an EoS using the same abundances. There remain inconsistencies between the two EoSs, coming from the non-ideal formulation used by Bifrost, and differences in partition functions for atoms and molecules. The internal energy (Eq. (6)) was recalculated with the EoS described in Sect. 2.2, using the hydrogen population levels, temperature and electron number densities of the original simulation. The relative difference in the internal energy has a rms value of 3 × 10⁻³, with a maximum of 6 × 10⁻² near the transition region.

The resulting simulation has 504 × 504 × 840 grid points, spanning 24 Mm in the horizontal direction and 16.8 Mm in the vertical, with the lower boundary at −2.44 Mm, where 0 is the averaged solar surface. This corresponds to a horizontal resolution of 47.6 km and a vertical resolution of 20 km. The diffusion scheme is described in Sect. 2.1, and a PDM coefficient of C_pdm = 2 was used for both the diffusion scheme and the population advection. The viscous and resistive heating was determined from the momentum and magnetic field fluxes as described in Rempel (2017). We enforced a minimum temperature of 2500 K to prevent over-cooling in the post-shock rarefactions. The upper boundary condition imposes a potential field and is open to outflows, but closed to inflows. The lower boundary condition is the open symmetric-field condition described by Rempel (2014).

We limited the Alfvén speed through the use of the Boris correction (semi-relativistic MHD with an artificially reduced speed of light), in combination with the dynamic limiting scheme described by Rempel (2017). The maximum speed of light used in the Boris correction is calculated C_max = max (2C_{s, max}, 3υ_max) in terms of the maximum sound speed (C_{s, max}) and velocity (υ_max) in the simulation domain. In addition, we imposed a dynamic ceiling on velocity (υ_max) and internal energy (ϵ_max) in order to prevent extreme values, only realised in a few grid points, dominating the simulation time step. The value of the ceiling υ_max was chosen dynamically so that it affected fewer than one in 1 million grid points. If more than 9.4 × 10⁻⁷ of the simulation grid points are above 0.95 of the limited υ_max (or ϵ_max), then the value is increased by 1%. If fewer than 4.7 × 10⁻⁷ points are above 0.95 of the limit then the value is lowered by 1%. The simulation presented has a maximum velocity of 300–400 km s⁻¹. The maximum speed of light used in the Boris correction is then calculated C_max = max 2C_s,max, 3υ_max). We do not allow the speed of light in the box to decrease below 2000 km s⁻¹. The chosen limits on velocity and speed of light in the box ensure a minimal effect on the chromospheric structure and dynamics. The impact of the choice of maximum speed of light has been studied for strong field active region simulations (Rempel 2017; Warnecke & Bingert 2020). The simulation presented in this work has a significantly weaker magnetic field than the active region simulations, and far fewer grid points will require limiting.

When the simulation was started the potential field upper boundary condition causes the coronal field to become more vertical. The transition region then produces a large transient. Differences resulting from the interpolation, and slight differences in the way the EoS was constructed, likely contribute to this transient. To reduce the timescale of the transient the diffusion was temporarily increased on velocities over 100 km s⁻¹ at the start of the simulation. The simulation was run until this transient passed and the RMS velocity stabilised, which took about 600 s. The additional diffusion and damping was then slowly removed and the simulation was run for an hour until the corona reheated.

This simulation cannot be directly compared to the original Bifrost public snapshot due to the differences resulting from the large transient and the lower viscosity and resistivity. The resulting model is shown in Fig. 1. The photosphere shows a bipolar enhanced network regions contain strong field concentrations of 1–3 kG. The magnetic field in the chromosphere (panel d) is dominated by the large-scale bipolar fields, with finely structured strands, including regions of opposite polarity.

The chromospheric dynamics are shown in Fig. 2³, in the mid chromosphere (panels a and b, z = 1 Mm) they are dominated by shocks. The shock-fronts show velocities over 20 km s⁻¹ in the inter-network regions. Above the strong network fields the shocks are suppressed and the dynamics follows the magnetic field. The temperatures range from the minimum value of 2.5 kK in the shock rarefactions to above 10 kK in the shock fronts. Above the network fields the temperatures are higher, reaching nearly 100 kK in regions where the transition region is depressed. The chromospheric velocity field at 2 Mm shows strong shocks, with velocities above 25 km s⁻¹. At 3 Mm loops are seen between the bipolar fields, while the inter-network regions show the shock canopy with temperatures around 1 kK.

The time-evolution of the chromospheric shocks is shown in a time-distance diagram in Fig. 3. In the low chromosphere (panel a), the quiet regions show shock fronts with a period of ≈5 min. In the network field regions periodic brightenings are seen at a x = 16 Mm with a period of ≈3 min. Due to the treatment of helium in LTE, a large fraction of the plasma at 3 Mm height plasma sits at around 10 kK, the preferred temperature of the first ionisation stage of helium.

Fig. 1 Snapshot of the simulation after 3850 s, showing: the vertical velocity at the photosphere (panel a) the normalised intensity at 500 nm (panel b), the vertical magnetic field at the photosphere (panel c), and the vertical magnetic field at a height of 2 Mm in the chromosphere (panel d). The dashed green lines show the slices taken in Fig. 2. Panels c and d have been saturated in order to show the fine structure of the magnetic field. Slices at the photosphere are taken at the contour where τ500 = 1.

4 Radiative Cooling and Heating

The radiative losses and heating in the simulation can be split into contributions from the multi-group scheme in and below the temperature minimum and low chromosphere, the chromospheric line losses, and the optically thin coronal losses. Figure 4 shows these different components for a slice through the model. The multi-group RT scheme cools and heats the photosphere and shocks near the temperature minimum. In the low chromosphere, the calcium and magnesium losses are strongest. The hydrogen losses are strong throughout the chromosphere, dominating the upper chromosphere and peaking in the lower transition region. The optically thin losses dominate above the transition region and are strongest in a narrow region immediately above the transition region.

The prescription for Lyman-alpha and the Lyman-continuum do not provide any heating in the upper chromosphere. The prescriptions for calcium and magnesium can provide a small amount of heating when the temperature decreases below 3.5 kK. The extreme ultraviolet (EUV) back-heating is strongest in the upper chromosphere, below the transition region, where neutral helium can form. The angle averaged intensity and heating rate of the EUV bin are shown in Fig. 5 for a slice through the simulation. High heating rates, above 10¹⁰ ergg⁻¹ S⁻¹, are strongly localised near areas of the transition region where the optically thin losses are high.

To investigate the relative importance of the different lines we plot the cooling/heating timescales in Fig. 6. The ionisation and recombination times in the chromosphere are long, preventing fast recombination of hydrogen as the gas is cooled. We therefore calculate the timescales τ = E_inst/Q using the terms in the internal energy that can instantaneously change E_inst = 3/2Nk_BT+E_nonH+E_mol,CE. This includes the microscopic kinetic energy, the ionisation of non-hydrogen species and the formation of H⁻ and ${\text{H}}_2^{+}$ molecules in chemical equilibrium. Figure 6 shows histograms of the timescale with temperature, for the chromospheric line losses and back-heating. Calcium cooling extends to lower temperatures, affecting shocks down to the temperature minimum, but it is lower than magnesium and hydrogen through the mid-to-high chromosphere. Magnesium cooling is strongest in regions below ≈10 kK, and hydrogen dominates radiative losses from ≈10 kK to the transition region, reaching timescales lower than 10 s. These results are similar to those presented in Carlsson & Leenaarts (2012), with hydrogen dominating the cooling above the mid chromosphere (≈ 1.5 Mm) and being marginally lower than magnesium in the low chromosphere.

Fig. 2 Slices through the chromosphere after 3850 s of the simulation, showing temperature in the left column and vertical velocity in the right column. Slices in the horizontal plane are taken through the low chromosphere (panels a and b) and the upper chromosphere (panels ɡ and h). The middle rows show vertical slices through a quieter inter-network region (panels c and d) and through the centre of the network field (panels e and f). An animation is available online.

Fig. 3 Time-distance diagram of temperature spanning 14 min of the simulation. The slice is taken at y = 12 Mm and at two heights, z = 1 Mm (panel a) and z = 2 Mm (panel b).

Fig. 4 Radiative cooling in the photosphere and chromosphere for a slice through the simulation. Panel a shows the temperature, and panel b shows the losses through the three-dimensional multi-group radiation transport scheme. The chromospheric line losses are shown for magnesium, calcium, and hydrogen in panels c, d, and e, respectively. Panel f shows the optically thin losses in the corona.

Fig. 5 Radiative back-heating of the upper chromosphere due to optically thin coronal losses. Panel a shows the intensity due to the coronal losses, and panel b shows the heating in the chromosphere.

Fig. 6 Histograms of radiative timescales in the chromosphere with temperature, showing the recipes for calcium (panel a), magnesium (panel b), and hydrogen losses (panel c). Panel d shows the back-heating from optically thin losses in the corona.

5 Hydrogen Populations

To investigate the effects of the NE treatment of hydrogen on the thermodynamics of the simulation we plot histograms of the temperature and electron number density in Fig. 7. The results of the NE EoS are compared to those calculated in LTE. The preferred temperature of hydrogen ionisation, around 8 kK is a prominent feature in the LTE results, but it is inconspicuous in the NLTE case. Two preferred temperature bands remain at 10 kK and 30 kK, caused by the first and second ionisation stages of helium, which is treated in LTE (Leenaarts et al. 2011; Golding et al. 2016). The spread of temperatures in the chromosphere is wider due to the long recombination timescales preventing hydrogen ionisation from buffering the temperature fluctuations. The effect of the higher ionisation fraction is also seen in the electron number density, which is higher than LTE in the low to mid chromosphere.

A detailed look at the hydrogen populations can be seen in Fig. 8. In order to compare the NE simulation with LTE we calculate the departure coefficient b of a quantity X as the ratio of the value in NE, to the value calculated using the LTE EoS, b_X = X_NE/X_LTE. The departure coefficient of temperature (b_T, panel b) shows up to 35% higher temperature in shocks and the transition region in the NE simulation, while behind the shocks the temperature is reduced by up to 25%. The ionisation fraction (panel c) is smooth throughout the chromosphere, as the long recombination timescales prevent neutral formation in the inter-shock regions. The departure coefficient for molecular hydrogen ($b_{{\text{H}}_2}$, panel h) is 1 in the photosphere and the temperature minimum, and around unity in cold chromospheric pockets. In warmer regions of the mid chromosphere the departure coefficient can be temporarily enhanced, and in hot shocks and the upper chromosphere it is reduced. The increase (or decrease) in H₂ departure coefficient occurs largely in locations where the temperature departure coefficient is decreased (or increased). The departure coefficients of the hydrogen ground state (b₁, panel e), and first excited state (b₂, panel f) in the corona are b₁ = 1.7 × 10⁴ and b₂ = 4 × 10⁴, similar to the values of Leenaarts et al. (2007). Although small regions with b₁ > 10⁹ are seen, we do not see large regions with extremely high departure coefficients just above the transition region (b₁ > 10¹⁰), as observed in Leenaarts et al. (2007). The departure coefficients of the Bifrost public simulation (Carls son et al. 2016) show similar magnitudes and behaviour as those shown in Fig. 8. They are shown for the initial snapshot, calculated from the Bifrost code, in Fig. F.1.

Fig. 7 Histograms of the temperature (top) and electron number (bottom) throughout the photosphere and chromosphere of the simulation. We compare the NE values (left columns) with those obtained from the pre-tabulated LTE EoS (right columns).

6 Numerical Performance

In this section we investigate the numerical cost of the newly implemented routines. These simulations were performed on the Max-Planck Computation data facilities ‘Raven’ cluster. This cluster contains 1592 compute nodes, each consisting of Intel Xeon IceLake-SP processors (Platinum 8360Y) processors, with 72 cores run at 2.4 GHz and connected with Mellanox HDR InfiniBand network (100 Gbit s⁻¹) interconnects. For the results presented in this section we used 20 nodes, or 1440 cores. The MURaM code is written with message passing interface (MPI) communication and does not support hybrid shared memory calculations.

The simulation presented in this work, row 1 of Table 1, was performed with an Alfvén speed limit of 2000 km s⁻¹, has a typical time step of 8.89 ms. We calculated the average time per iteration from 200 time-step updates, the expected computational time per grid cell update and the cost of simulating 1 h of solar evolution. A summary of the timing, the time step and the wall time are presented in Table 1. The simulation setup shown in the paper will take 57.30 µs per grid point per core. This gives a wall time cost of approximately 1.38 million CPU-hours (Mcore-h) per hour of simulated time.

To determine the computational costs of the new physics implemented in this work we performed a number of test simulations. The simplest of these is a simulation with a LTE EoS, grey LTE multi-group RT, and no back-heating due to coronal EUV radiation. This setup, row 2 of Table 1, is similar to that presented by Rempel (2017) utilising the coronal extension to the MURaM code and costs 0.16 Mcore-h per hour of solar time (6.61 µs per grid point per core). This grey LTE simulation spends 14% of the computational time on the RT modules, 66% on MHD, and 5% and 7.5% on the EoS and diffusion treatments, respectively. By comparison, in the chromospheric simulations presented in this work, the computational cost is dominated by the EoS, and to a lesser extent the RT.

First, we consider the effects of including the extended RT modules. Including the EUV back heating of the chromosphere, row 3 of Table 1, increases the cost of RT by almost 300%. This corresponds to a 26% increase in the total computational time. This large increase is due to the optically thin nature of the EUV radiation in the corona, rays can cross many computational sub-domains and take more iterations to converge. Including a four-band scattering formulation, row 4 of Table 1, for the three-dimensional multi-group radiation scheme further increases the computational cost of the simulation by 90%. Including the more realistic treatment of RT makes RT the most expensive component of the simulation, requiring 64% of the computational time. This large increase is from additional iterations of the strongly scattering sub-bins. Most of the radiation groups converge quickly, in two or three iterations, similar to the grey radiation bin. However, the optically thin chromospheric lines bin can take up to seven iterations to converge.

The greatest computational cost is the inclusion of the NE ionisation of hydrogen in the EoS. The new module increases the total runtime by 360%. This includes the solution of the hydrogen rate equations, as well as overhead in the MHD and diffusion modules. The increase in the MHD and diffusion modules is caused by advection and diffusion of the atomic populations, and calls to the EoS to maintain consistency of the solution between directional sweeps of the diffusion routine. The simulation setup presented in this paper requires 8.6 times more computational power per grid point than the LTE coronal simulations presented in Rempel (2017).

In addition, we varied two approximations that have a significant effect on the computational cost. Firstly, we performed two simulations where the minimum limit on the speed of light is changed. Increasing the speed of light, to 5000 km s⁻¹, row 5 of Table 1, reduces the time step leading to a 230% increase in the computational cost. Reducing the speed of light to 1000 km s⁻¹, row 6 of Table 1, allows a larger time step and lowers the computational cost by 37%.

A second choice that affects the computational cost is the frequency (in iterations) at which the radiation field is calculated. Updating the radiation field every iteration, row 7 of Table 1, increases the computational cost of RT by 300%. Reducing the frequency to every 10 iterations, row 8 of Table 1, reduces it by 40%. This change does not scale linearly as more iterations are needed to converge the radiation field to the required tolerance when it is calculated less regularly.

Fig. 8 Properties of the NE hydrogen populations through the centre of the enhanced network region. The panels show (a) temperature in the NE simulations, (b) the departure coefficient of temperature, bT = Tne/Tlte, (c) the NE ionisation fraction, Fi = nH,1/nH,tot, (d) the number density of the first excited level of hydrogen, nH,0,1, and departure coefficients of (e) the ground state, (f) the first excited state, (g) protons, and (h) molecular hydrogen.

7 Discussion and Conclusion

The current work improves on the original LTE MURaM code through the implementation of three new modules: an NE treatment of hydrogen in the EoS, NLTE tabulated losses in the chromosphere, and a scattering multi-group RT scheme. An initial simulation has been performed, beginning from the publicly available Bifrost snapshot. The simulation differs significantly from the original Bifrost model due to differences in the diffusion scheme, the potential boundary condition, and small differences in the EoS. A detailed comparison between the results from Bifrost and MURaM will be the subject of a separate paper.

The simulation shows the importance of an NE treatment of hydrogen in the chromosphere. The upper chromosphere is highly dynamic, with strong shocks and large departure coefficients of the ground state and n = 2 energy level. Despite the strong gradients and fine structure in velocity and temperature, the hydrogen populations in the upper atmosphere are smooth due to the long recombination times, relative to the dynamical timescales. The departure coefficients b₁ and b₂ are approximately 10³ to 10⁶ in the upper chromosphere. The calculated departure populations match those in the Bifrost code (Leenaarts et al. 2007; Carlsson et al. 2016). The temperatures are around 25% different from the LTE values, and the electron number density remains higher in cold shock expansions. These differences occur due to the inability of protons to recombine into neutral hydrogen before a new shock passes through the chromosphere. These differences will be important for the accurate synthesis of chromospheric spectral lines.

The current implementation of chromospheric radiative losses and the NE EoS are based upon a number of simplifying assumptions:

1. The tabulated chromospheric line losses ignore significant scatter around pre-tabulated values of escape probability, the ionisation fraction, and the radiative loss function. The optically thin formalism cannot simulate three-dimensional heating effects due to shocks and explosive events;

2. The radiation field used for the NE treatment of hydrogen is isotropic in the chromosphere;

3. The treatment of Lyman-alpha in radiative equilibrium for the NE treatment of hydrogen is inaccurate near the transition region (Carlsson & Stein 2002; Golding et al. 2016);

4. Atomic and molecular populations in the multi-group radiation transport scheme are treated in LTE;

5. Helium is treated in LTE.

The above approximations are necessary for the simulation of large three-dimensional models that include an NE chromosphere. They will also have a significant effect on the physics and chemistry acting in the chromosphere. It is important to investigate new methods to relax these assumptions. The work by Golding et al. (2016) has extended the NE EoS and chromospheric line cooling to be more accurate in the upper chromosphere and transition region. This includes an approximate Lyman-alpha bin, allowing for three-dimensional cooling and heating effects while ignoring scattering and PRD effects. Additionally, the optically thin losses have been split into six EUV frequency bands to incorporate the Lyman continuum and a simplified three-level helium atom. These improvements greatly increase the realism of the method in the upper chromospheric layers.

Another proposal for the fast NE treatment of atoms in a radiative MHD simulation is described by Judge (2017). This method uses the escape probability approach to quickly converge the populations. The method is easily extendable to elements other than hydrogen. The one-dimensional plane-parallel nature may lead to unrealistic variations between neighbouring horizontal pixels. It is also suggested that the EoS be decoupled from atomic populations and the problem solved in stages, where the populations are updated and then used to calculate the new temperature, electron number density, and pressure. This would reduce the complexity of the system of equations that are solved and allow the use of a pre-tabulated EoS, potentially offering a significant increase in speed when solving the NE problem.

Finally, a new time-implicit numerical method for solving the detailed NE radiative MHD problem, built as an extension to the MURaM code, is presented in Anusha et al. (2021). The formulation allows larger time steps to be taken when solving the system of rates, making the method promising for time-dependent multi-dimensional simulations.

Recent studies (Martinez-Sykora et al. 2012; Shelyag et al. 2016) have shown the importance of including ion-neutral interactions, in particular ambipolar diffusion, in simulations of the chromosphere. The collisional rates and electron number density are strongly tied to the hydrogen ionisation fraction, and this varies greatly between LTE and NE. Including NE ionisation reduces the impact of ambipolar diffusion on the simulation; however, a decreased efficiency of shock heating is observed, as is the heating of cool low-lying loops (Martinez-Sykora et al. 2020). Ambipolar diffusion and the Hall effect are included in the MURaM code (Cheung & Cameron 2012; Rempel & Przybylski 2021), and simulations that include both effects will soon be performed.

Movie

Movie 1 associated with Fig. 2 (fig2_movie) Access here

Acknowledgements

We thank the anonymous referee for suggestions that improved the paper. We would like to thank S. Danilovic for assisting with understanding of the MURaM code and facilitating collaboration with the Stockholm group. D.P. would like to acknowledge the help of I. Milic, K. Sowmya, H.N. Smitha, M. van Noort, R. Collet, and P. Judge for helpful discussion relating to NLTE and NE physics. D.P. would also like to thank the Bifrost group for helpful explanations of the Bifrost code, especially V. Hansteen. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 695075). We gratefully acknowledge the computational resources provided by the Cobra & Raven supercomputer systems of the Max Planck Computing and Data Facility (MPCDF) in Garching, Germany. J.L. was supported by a grant from the Knut and Alice Wallenberg foundation (2016.0019). This material is based upon work supported by the National Center for Atmospheric Research, which is a major facility sponsored by the National Science Foundation under Cooperative Agreement No. 1852977. L.S.A., V.W., and A.I.S. acknowledge support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant no. 715947).

Appendix A LTE Equation of State

The ideal EoS is pre-tabulated to calculate the thermodynamic variables n_e, T, and p in terms of density ρ and internal energy ϵ. To start, we calculate the LTE ionisation fractions for a given density ρ, electron number n_e, temperature T and set of abundances A_a, using the Saha-Boltzmann equation: $$\frac{n_{a,i+1}}{n_{a,i}}=\frac{2}{n_{\text{e}}}{\left(\frac{2\pi m_e{k}_{\text{B}}T}{h^2}\right)}^{1.5}\frac{U_{a,i+1}}{U_{a,i}}{e}^{-{\chi }_{a,i}/\left(k_{\text{B}}T\right)+\left(Z_{a,i}+1\right)\text{Δ}\mu },$$(A.1)

where k_B is Boltzmann’s constant, m_e the electron mass, h is Planck’s constant and Δμ the ionisation lowering due to interactions with surrounding particles. U_a,i is the partition function, χ_a,i the ionisation energy, Z_a,i the charge, and and n_a,i the number density for element a and ionisation stage i. Where available, we used the polynomial partition functions of Cardona et al. (2005), except for iron and nickel where the tables of Halenka et al. (2001) and Halenka & Madej (2002) were used. The H⁻ fraction was calculated from Eq. (A.1) and the H₂ and ${\text{H}}_2^{+}$ fractions are $$\frac{n_{H,0}{n}_{H,0}}{n_{H_{2,0}}}={\left(\frac{2\pi m_e{k}_{\text{B}}T}{h^2}\right)}^{1.5}\frac{U_{H,0}}{U_{H_2,0}}{e}^{-D_{{\text{H}}_2,0}/\left(k_{\text{B}}T\right)}\text{and}$$(A.2) $$\frac{n_{H,0}{n}_{H,1}}{n_{H_{2,1}}}={\left(\frac{2\pi m_e{k}_{\text{B}}T}{h^2}\right)}^{1.5}\frac{U_{H,0}}{U_{H_2,1}}{e}^{-D_{{\text{H}}_2,1}/\left(k_{\text{B}}T\right)},$$(A.3)

where m is the reduced mass and $D_{{\text{H}}_2}=4.478007\text{\hspace{0.17em}ev}$ and $D_{{\text{H}}_2^{+}}=2.650639\text{\hspace{0.17em}ev}$ are the dissociation energies. The partition functions of H₂ were taken from (Popovas & Jørgensen 2016) and the partition functions of ${\text{H}}_2^{+}$ are from a polynomial fit to the Stancil ((1994) table. Following Mihalas et al. (1988), the energy of the hydrogen species are $$E_{{\text{H}}_2}=k_{\text{B}}{T}^2\frac{\partial \ln U_{{\text{H}}_2,0}}{\partial T},$$(A.4) $$E_{{\text{H}}_2,1}=D_{{\text{H}}_2\text{,0}}-D_{{\text{H}}_2\text{,1}}+{\chi }_{\text{H,1}}+k_{\text{B}}{T}^2\frac{\partial \ln U_{{\text{H}}_2,1}}{\partial T},$$(A.5) $$E_{\text{H,}-1}=0.5D_{{\text{H}}_2}-{\chi }_{\text{H,}-1},and$$(A.6) $$E_{\text{H,}i,e}=0.5D_{{\text{H}}_2}+{\chi }_{\text{H,}i,e}\text{,}$$(A.7)

where χ_H,i,e is the excitation or ionisation energy of the hydrogen level (i, e). The ionisation energy of neutral hydrogen is χ_H,1 = 13.59844 ev, and for H⁻ is χ_H,−1 = 0.754 ev. Once the fractions of all elements and of the hydrogen molecules are determined we followed the formulation of the VMW EoS. This process involves iterating the ionisation and molecular fractions, electron number density and temperature until a convergence criteria is reached. We iterated this procedure until the electron number density converged to a tolerance of 1.0 × 10⁻⁸. Once converged, the energy and pressure are calculated $$E_{\mathrm{int}}=E_{\text{exi}}+E_{\text{trans}}+E_{\text{rad}}+E_{\text{C}}+E_{\text{PI}}\text{\hspace{0.17em}and}$$(A.8) $$p=p_{\text{trans}}+p_{\text{rad}}+p_{\text{C}}+p_{\text{PI}},$$(A.9)

in terms of excitation and ionisation (exi), translational (trans), radiation (rad), Coloumb (C), and pressure ionisation (PI) components. We ignore the effects of degenerate and relativistic electrons as they are small within the physical regime of these simulations. The radiation energy and pressure terms are also small and not included. The translational components are calculated as p_trans = (Σ_a,i n_a,i + n_e) k_BT and E_trans = 3/2 p_trans.

The contribution to the energy from excitation and ionisation is calculated as $$E_{\text{exi}}={\displaystyle \sum _{a,i}\left({\chi }_{a,i}+k_{\text{B}}{T}^2\frac{\partial \ln U_{a,i}}{\partial T}\right)}{n}_{a,i}.$$(A.10)

For the Coulomb correction we followed the prescription of Mihalas et al. (1988). The free energy, ignoring electron density, is calculated for a parcel of gas of volume V and a particle number N = nV as $$F_{\text{C}}=-\frac{2{\pi }^{1/2}{e}^3}{3k_{\text{B}}^{{}^{1/2}}}\frac{1}{{\left(VT\right)}^{{}^{1/2}}}{\left({\displaystyle \sum _{a,i}{N}_{a,i}{Z}_{a,i}^2+N_{\text{e}}}\right)}^{3/2}\tau \left(x\right),$$(A.11)

where the function τ is $$\tau \left(x\right)=3x^{-3}\left(\ln \left(1+x\right)-x+\frac{1}{2}{x}^2\right),$$(A.12)

and x is given by $$\begin{array}{ll}x=\hfill & \frac{4{\pi }^{1/2}{e}^3}{3k_{\text{B}}^{3/2}}\frac{1}{V^{1/2}{T}^{3/2}}\hfill \\ \hfill & \frac{\left({\displaystyle {\sum }_{a,i}{N}_{a,i}{Z}_{a,i}}\right)}{\left({\displaystyle {\sum }_{a,i}{N}_{a,i}}\right)}{\left({\displaystyle \sum _{a,i}{N}_{a,i}{Z}_{a,i}^2+N_{\text{e}}}\right)}^{1/2}.\hfill \end{array}$$(A.13)

The required derivatives of τ and x are $$\frac{\partial \tau }{\partial x}=-\frac{3}{2}\left(\tau -\frac{1}{1+x}\right),$$(A.14) $$\frac{\partial x}{\partial N_{\text{e}}}=x\frac{1}{2\left({\displaystyle {\sum }_{a,i}{N}_{a,i}{Z}_{a,i}^2+N_{\text{e}}}\right)},$$(A.15) $$\frac{\partial x}{\partial V}=-x\frac{1}{2V},and$$(A.16) $$\frac{\partial x}{\partial V}=-x\frac{3}{2T}.$$(A.17)

The reduction of the ionisation potential Δµ_C, used in the Saha-Boltzmann equation (Eq. (A.1)), is $$\begin{array}{ll}\text{Δ}{\mu }_{\text{C}}\hfill & =-\frac{1}{k_{\text{B}}T}\frac{\partial F_{\text{C}}}{\partial N_{\text{e}}}\hfill \\ \hfill & =-\frac{F_{\text{C}}}{k_{\text{B}}T}\left(\frac{3}{2}\frac{1}{{\displaystyle \sum _{a,i}{N}_{a,i}{Z}_{a,i}^2}+N_{\text{e}}}+\frac{1}{\tau }\frac{\partial \tau }{\partial x}\frac{\partial x}{\partial N_{\text{e}}}\right),\hfill \end{array}$$(A.18)

giving a pressure correction of $$P_{\text{C}}=-\frac{\partial F_{\text{C}}}{\partial V}=F_{\text{C}}\left(\frac{1}{2V}-\frac{1}{\tau }\frac{\partial \tau }{\partial x}\frac{\partial x}{\partial V}\right)$$(A.19)

and an internal energy correction of $${є}_{\text{C}}=-T^2\frac{\partial }{\partial T}\left(\frac{F_{\text{C}}}{T}\right)=F_{\text{C}}\left(\frac{3}{2}-\frac{1}{\tau }\frac{\partial \tau }{\partial x}T\frac{\partial x}{\partial T}\right).$$(A.20)

Simplifying Eqs. A.19–A.20 and taking V = 1.0 cm⁻³ gives $$\begin{array}{ll}\text{Δ}{\mu }_C=\hfill & \frac{1}{2k_{\text{B}}T}\left(3\tau +\frac{\partial \tau }{\partial x}x\right)\hfill \\ \hfill & \frac{2{\pi }^{1/2}{e}^3}{3k_{\text{B}}^{1/2}}\frac{1}{T^{1/2}}{\left({\displaystyle \sum _{a,i}{n}_{a,i}{z}_{a,i}^2+n_{\text{e}}}\right)}^{1/2},\hfill \end{array}$$(A.21) $$\begin{array}{ll}{p}_C=\hfill & \frac{1}{2}\left(\tau +\frac{\partial \tau }{\partial x}x\right)\hfill \\ \hfill & \frac{2{\pi }^{1/2}{e}^3}{3k_{\text{B}}^{1/2}}\frac{1}{T^{1/2}}{\left({\displaystyle \sum _{a,i}{n}_{a,i}{z}_{a,i}^2+n_{\text{e}}}\right)}^{3/2},\text{\hspace{0.17em}and}\hfill \end{array}$$(A.22) $$\begin{array}{ll}{є}_C=\hfill & \frac{3}{2}\left(\tau +\frac{\partial \tau }{\partial x}x\right)\hfill \\ \hfill & \frac{2{\pi }^{1/2}{e}^3}{3k_{\text{B}}^{1/2}}\frac{1}{T^{1/2}}{\left({\displaystyle \sum _{a,i}{n}_{a,i}{z}_{a,i}^2+n_{\text{e}}}\right)}^{3/2}.\hfill \end{array}$$(A.23)

Additionally, the pressure ionisation device described in the Eggleton, Faulkner, and Flannery EoS (Eggleton et al. 1973) was included. This method provides a thermodynamically consistent result that gives a qualitatively correct pressure ionisation as density increases: $$F{}_{PI}=\frac{\text{Ω}\left(T\right)}{V}\left(N_{\text{e,0}}^2-N_{\text{e}}^2\right),$$(A.24)

where N_e,0 is the electron number of the gas when it is fully ionised and Ω is $$\text{Ω\hspace{0.17em}\hspace{0.17em}}=\frac{a_0^3}{2}\left(k_{\text{B}}T+{20}_{\chi 0}\right),$$(A.25)

where a₀ = 5.23e⁻⁹/〈Z〉 cm, χ0 = 2.16e⁻¹¹ 〈Z〉² erg, and 〈Z〉 is the mean charge per nucleus and the required derivative, $$\frac{\partial \text{Ω}}{\partial T}=\frac{a_0^3}{2}{k}_{\text{B}}.$$(A.26)

The change in the potential can then be calculated as $$\text{Δ}{\mu }_{PI}=-\frac{1}{k_{\text{B}}T}\frac{\partial F_{PI}}{\partial N_{\text{e}}}=\frac{2\text{Ω}}{k_{\text{B}}T}\left(T\right)n_{\text{e}},$$(A.27)

the pressure correction as $$P_{PI}=-\frac{\partial F_{PI}}{\partial V}=\text{Ω\hspace{0.17em}}\left(T\right)\left(n_{e,0}^2-n_e^2\right),$$(A.28)

and finally the internal energy correction as $${є}_{PI}=-T^2\frac{\partial }{\partial T}\left(\frac{F_{PI}}{T}\right)=\left(\text{Ω\hspace{0.17em}}\left(T\right)-VT\frac{\partial \text{Ω\hspace{0.17em}}\left(T\right)}{\partial T}\right)\left(n_{e,0}^2-n_{\text{e}}^2\right).$$(A.29)

Once a complete solution is obtained for all required density and energy values, the entropy s is calculated by integrating over the table, $$s={\displaystyle {\int }_{є}{\displaystyle {\int }_{\rho }\left(\frac{є}{T}\partial \ln є-\frac{p}{T\rho }\partial \ln \rho \right).}}$$(A.30)

Appendix B Derivatives of the Non-Equilibrium Equation of State

EoS derivatives for chemical equilibrium of ${\text{H}}_2^{+}$, $$\begin{array}{ll}\frac{\partial n_{{\text{H}}_{2,1}}}{\partial T}\hfill & =-n_{{\text{H}}_{2,1}}(\frac{3}{2T}+\frac{D_{{\text{H}}_{2,1}}}{k_{\text{B}}{T}^2}\hfill \\ \hfill & +\frac{1}{U_{\text{H,0}}}\frac{\partial U_{H,0}}{\partial T}-\frac{1}{U_{{\text{H}}_{\text{2,1}}}}\frac{\partial U_{{\text{H}}_{2,1}}}{\partial T}),\hfill \end{array}$$(B.1) $$\frac{\partial n_{{\text{H}}_{2,1}}}{\partial n_{\text{H,0,}j}}=\frac{n_{{\text{H}}_{2,1}}}{{\Sigma }_j{n}_{\text{H,0,}j}},\text{and}$$(B.2) $$\frac{\partial n_{{\text{H}}_{2,1}}}{\partial n_{\text{H,}1}}=\frac{n_{{\text{H}}_{2,1}}}{n_{\text{H,}1}},$$(B.3)

and for H⁻¹, $$\frac{\partial n_{\text{H,-1}}}{\partial T}=-n_{\text{H,-1}}\left(\frac{3}{2T}+\frac{{\chi }_{\text{H,-1}}}{k_{\text{B}}{T}^2}+\frac{1}{U_{H,0}}\frac{\partial U_{\text{H,0}}}{\partial T}\right),$$(B.4) $$\frac{\partial n_{\text{H,-1}}}{\partial n_{\text{e}}}=\frac{n_{\text{H,-1}}}{n_{\text{e}}},and$$(B.5) $$\frac{\partial n_{\text{H,-1}}}{\partial n_{\text{H,0,}j}}=\frac{n_{\text{H,-1}}}{{\Sigma }_j{n}_{\text{H,0,}j}}.$$(B.6)

The derivatives for the energy conservation equation, ƒ₀, are $$\begin{array}{ll}\frac{\partial f_0}{\partial T}\hfill & =-\frac{1}{E_{\mathrm{int}}}(\frac{3k_{\text{B}}}{2}[n_{\text{e}}+n_{\text{nonH}}+n_{{\text{H}}_2}+n_{{\text{H}}_{\text{2,1}}}\hfill \\ \hfill & +n_{\text{H,-1}}+{\displaystyle \sum _{i,j}{n}_{\text{H},i,j}}\left]+\frac{3k_{\text{B}}T}{2}\right[\frac{\partial n_{noH}}{\partial T}+\frac{\partial n_{\text{H,-1}}}{\partial T}\hfill \\ \hfill & +\frac{\partial n_{{\text{H}}_2,1}}{\partial T}]+n_{\text{H,tot}}\frac{\partial E_{\text{nonH}}}{\partial T}+\frac{\partial E_{{\text{H}}_2}}{\partial T}{n}_{{\text{H}}_2}\hfill \\ \hfill & +\frac{\partial E_{\text{H,-1}}{n}_{\text{H,-1}}}{\partial T}+\frac{\partial E_{H_{2,1}}{n}_{{\text{H}}_{2,1}}}{\partial T}),\hfill \end{array}$$(B.7) $$\begin{array}{ll}\frac{\partial f_0}{\partial n_{\text{e}}}\hfill & =-\frac{1}{E_{\mathrm{int}}}(\frac{3k_{\text{B}}T}{2}\left[1.0+\frac{\partial n_{\text{H,}-1}}{\partial n_{\text{e}}}\right]+n_{\text{H,tot}}\frac{\partial E_{\text{nonH}}}{\partial n_{\text{e}}}\hfill \\ \hfill & +\frac{\partial n_{\text{H,}-1}}{\partial n_{\text{e}}}{E}_{\text{H,-1}}),\hfill \end{array}$$(B.8) $$\begin{array}{ll}\frac{\partial f_0}{\partial n_{\text{H,0,}j}}\hfill & =-\frac{1}{E_{\mathrm{int}}}(\frac{3k_{\text{B}}T}{2}+E_{\text{H,0,}j}\hfill \\ \hfill & +\frac{\partial n_{{\text{H}}_{2,1}}}{\partial n_{\text{H,0,}j}}{E}_{{\text{H}}_{2,1}}+\frac{\partial n_{\text{H,-1}}}{\partial n_{\text{H,0,}j}}{E}_{\text{H,-1}}),and\hfill \end{array}$$(B.9) $$\frac{\partial f_0}{\partial n_{\text{H,1}}}=-\frac{1}{E_{\mathrm{int}}}\left(\frac{3k_{\text{B}}T}{2}{E}_{\text{H,1}}+\frac{\partial n_{{\text{H}}_{2,1}}}{\partial n_{\text{H,1}}}{E}_{{\text{H}}_{2,1}}\right).$$(B.10)

The derivatives for the charge conservation equation, f₁ , are $$\frac{\partial f_1}{\partial T}=-\frac{1}{n_{\text{e}}}\left(-\frac{\partial n_{\text{H,-1}}}{\partial T}+n_{\text{H,tot}}\frac{\partial n_{\text{e,nonH}}}{\partial T}\right),$$(B.11) $$\begin{array}{ll}\frac{\partial f_1}{\partial n_{\text{e}}}\hfill & =-\frac{1}{n_{\text{e}}}\left(1-\frac{\partial n_{\text{H,-1}}}{\partial n_{\text{e}}}+n_{\text{H,tot}}\frac{\partial n_{\text{e,nonH}}}{\partial n_{\text{e}}}\right)\hfill \\ \hfill & +\frac{1}{n_{\text{e}}^2}\left(n_{\text{H,}1}+n_{{\text{H}}_{2,1}}-n_{\text{H,}-1}+n_{\text{H,tot}}{n}_{\text{e,nonH}}\right)\hfill \end{array}$$(B.12) $$\frac{\partial f_1}{\partial n_{\text{H,0},j}}=-\frac{1}{n_{\text{e}}}\left(\frac{\partial n_{{\text{H}}_{2,1}}}{\partial n_{\text{H,0,}j}}+\frac{\partial n_{\text{H,-1}}}{\partial n_{\text{H,0,}j}}\right),$$(B.13) $$\frac{\partial f_1}{\partial n_{\text{H,1}}}=-\frac{1}{n_{\text{e}}}\left(1-\frac{\partial n_{{\text{H}}_{2,1}}}{\partial n_{\text{H,1}}}\right).$$(B.14)

The derivatives for the hydrogen nucleus conservation equation, f₂ , are $$\frac{\partial f_2}{\partial T}=-\frac{1}{n_{H,tot}}\left(1+\frac{\partial n_{{\text{H}}_{2,1}}}{\partial T}+\frac{\partial n_{\text{H,-1}}}{\partial T}\right),$$(B.15) $$\frac{\partial f_2}{\partial n_{\text{e}}}=-\frac{1}{n_{H,tot}}\frac{\partial n_{\text{H,-1}}}{\partial n_{\text{e}}},$$(B.16) $$\frac{\partial f_2}{\partial n_{\text{H},0,j}}=-\frac{1}{n_{H,tot}}\left(1+\frac{\partial n_{{\text{H}}_{2,1}}}{\partial n_{\text{H},0,j}}+\frac{\partial n_{\text{H,-1}}}{\partial n_{\text{H},0,j}}\right),$$(B.17) $$\frac{\partial f_2}{\partial n_{\text{H},1}}=-\frac{1}{n_{H,tot}}\frac{\partial n_{{\text{H}}_{2,1}}}{\partial n_{\text{H,1}}}.$$(B.18)

The derivatives of the rate-equations, f_2+i, are $$\begin{array}{ll}{f}_{3+ij}\hfill & =\frac{\text{Δ}t}{n_{aij}^{t_0}}({\displaystyle \sum _{kl\ne ij}{n}_{akl}}\frac{\partial P_{kl,ij}}{\partial T}-n_{aij}{\displaystyle \sum _{kl\ne ij}\frac{\partial P_{ij,kl}}{\partial T}}\hfill \\ \hfill & +{\displaystyle \sum _{r^{+}}{n}_{r^{+},A}}{n}_{r^{+},B}{n}_{r^{+},C}\frac{\partial K_{r^{+}}}{\partial T}\hfill \\ \hfill & -{\displaystyle \sum _{r^{-}}{n}_{r^{-},A}{n}_{r^{-},B}{n}_{r^{-},C}\frac{\partial K_{r^{-}}}{\partial T}}),\hfill \end{array}$$(B.19) $$\begin{array}{ll}\frac{\partial f_{3+ij}}{\partial n_e}\hfill & =-\frac{\text{Δ}t}{n_{aij}^o}({\displaystyle \sum _{kl\ne ij}{n}_{akl}}\frac{\partial P_{kl,ij}}{\partial n_{\text{e}}}\hfill \\ \hfill & -n_{aij}{\displaystyle \sum _{kl\ne ij}\frac{\partial P_{ij,kl}}{\partial n_{\text{e}}}}),\hfill \end{array}$$(B.20) $$\frac{\partial f_{3+ij}}{\partial n_{aij}}=\frac{1}{n_{aij}^o}+\frac{\text{Δ}t}{n_{aij}^o}{\displaystyle \sum _{kl\ne ij}{P}_{ij,kl},}$$(B.21) $$\frac{\partial f_{3+ij}}{\partial n_{akl}}=-\frac{\text{Δ}t}{n_{aij}}{P}_{kl,ij},$$(B.22) $$\frac{\partial f_{3+ij}}{\partial n_{r,A}}=-\frac{\text{Δ}t}{n_{aij}}{n}_{r,B}{n}_{r,C}{K}_r.$$(B.23)