11institutetext: Max-Planck Institute for Solar System Research, 37077
Göttingen, Germany 22institutetext: School of Space Research, Kyung Hee
University, Yongin, Gyeonggi 446-701, Republic of Korea 33institutetext: High
Altitude Observatory, NCAR, P.O. Box 3000, Boulder, Colorado 80307, USA
44institutetext: Institute for Solar Physics, Dept. of Astronomy, Stockholm
University, AlbaNova University Centre, SE-10691 Stockholm, Sweden

# Chromospheric Extension of the MURaM Code

D. Przybylski 11 R. Cameron 11 S.K. Solanki 1122 M. Rempel 33 J. Leenaarts 44
L. S. Anusha 11 V. Witzke 11 A.I. Shapiro 11<EMAIL_ADDRESS>
###### Abstract

Context. Detailed numerical models of chromosphere and corona are required to
understand the heating of the solar atmosphere. An accurate treatment of the
solar chromosphere is complicated by the effects arising from Non Local
Thermodynamic Equilibrium (NLTE) radiative transfer. A small number of strong,
highly scattering lines dominate the cooling and heating in the chromosphere.
Additionally, the recombination times of ionised hydrogen are longer than the
dynamical timescales, requiring a non-equilibrium (NE) treatment of hydrogen
ionisation.

Aims. We describe a set of necessary additions to the MURaM code so that it
might handle some of the important NLTE effects. We investigate the impact on
models of the solar chromosphere caused by NLTE and NE effects in radiation
magnetohydrodynamic (rMHD) simulations of the solar atmosphere.

Methods. The MURaM code is extended to include the physical process required
for accurate simulation of the solar chromosphere, as implemented in the
Bifrost code. This includes a time-dependent treatment of hydrogen ionisation,
a scattering multi-group radiation transfer scheme and approximations for NLTE
radiative cooling.

Results. The inclusion of NE and NLTE physics has a large impact on the
structure of the chromosphere; the NE treatment of hydrogen ionisation leads
to a higher ionisation fraction and enhanced populations in the first excited
state throughout cold inter-shock regions of the chromosphere. Additionally
this prevents hydrogen ioniation from buffering energy fluctuations, leading
to hotter shocks and cooler inter-shock regions. The hydrogen populations in
the ground and first excited state are enhanced by $10^{2}-10^{3}$ in the
upper chromosphere and up to $10^{9}$ near the transition region.

Conclusions. Including the necessary NLTE physics leads to significant
differences in chromospheric structure and dynamics. The thermodynamics and
hydrogen populations calculated using the extended version of the MURaM code
are consistent with previous non-equilibrium simulations. The electron number
and temperature calculated using the non-equilibrium treatment of the
chromosphere are required to accurately synthesise chromospheric spectral
lines.

###### Key Words.:

magnetohydrodynamics, radiative transfer, Sun:chromosphere

##  1 Introduction

The importance of the solar chromosphere for resolving a number of the large
open questions in solar physics is undoubted. E.g., 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 and the local
plasma heated, as well as the solar wind 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 is necessary to model the
structure of the chromosphere.

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

A multi-dimensional treatment of the chromosphere is complicated by the non-
locality of 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 3D spectral synthesis of a single
chromospheric 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 radiation transfer
($23.23\leavevmode\nobreak\ \mathrm{s}$ per grid point for hydrogen in PRD)
compared to an update of the multi-group LTE rMHD code
($15.80\leavevmode\nobreak\ \mu\mathrm{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 multi-dimensional 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/recombination of hydrogen and
helium in the solar chromosphere. For multi-dimensional simulations, we turn
to approximations of these effects in order to make 3D simulations
computationally tractable.

The 3D computational modelling of the solar choromosphere has been so far 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
3D simulations of the chromosphere currently available. The approximations
include tabulated recipes for computationally efficient NLTE chromospheric
line losses, based on detailed synthesis of the radiation field including PRD
effects (Carlsson & Leenaarts, 2012), scattering multi-group radiation
transfer (Skartlien, 2000; Hayek et al., 2010), and a NE hydrogen Equation of
State (EoS) (Leenaarts et al., 2007). To more accurately simulate the
transition region, Bifrost includes a 3D 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) which 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 extensions of Golding et
al. (2016).

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), magnetic flux
emergence (Cheung et al., 2007; Chen et al., 2017), etc. and has been used to
investigate the effect of ambipolar-diffusion in 3D rMHD simulation (Cheung &
Cameron, 2012; Danilovic, 2017). The code has 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). It has also 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 close one of the large remaining gaps in its
ability to treat phenomena in the solar atmosphere, namely covering the solar
chromosphere.

In Sect. 2 we describe the numerical methods used, including the diffusion
scheme, NE equation of state 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 analyze the radiative
cooling in the chromosphere, and in Sect. 5 we discuss the effect of the NE
treatment of hydrogen in the chromosphere. 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 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 4-stage explicit time-update scheme. The code’s original
hyperdiffusion 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
Equation of State and the implementation of a non-equilibrium 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 $\Delta t$. First the radiative heating and cooling (Sect 2.5) are
calculated using the system state at the previous time $t_{0}$. In each sub-
stage $n$ the steps followed are: calculate the right-hand-side (RHS) of the
MHD equations using a directionally unsplit approach (Sect. 2.6), integrate in
time to the new state $t_{*}=t_{0}+1/(5-n)\Delta t$, then apply the boundary
conditions, advect the populations, and evaluate the EoS (Sect. 2.2). No
hyperdiffusion 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 $\nabla\cdot B$ cleaner (Dedner et al., 2002)
is applied. After the MHD variables have been updated to the next timestep,
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 which is based around slope-
limiters and includes higher order hyper-diffusion 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 undershoot or overshoot
that occurs near discontinuities when a high-order (HO) scheme is used.

To apply the limiter to a quantity $\Phi$, the flux through the interface
$i+1/2$ between cell $i$ and $i+1$ must be calculated. We use a 4th order
centred reconstruction to calculate the high-order value at the cell interface
$\Phi_{i+1/2}^{\mathrm{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.

$\displaystyle\Phi^{l}_{i+1/2}$ $\displaystyle=$
$\displaystyle\Phi^{HO}_{i+1/2}-\mathrm{sign}\left(\nabla^{+}\right)\mathrm{max}\left[0,|\Phi^{HO}_{i+1/2}-\Phi_{i}|\right.$
(1) $\displaystyle-$
$\displaystyle\left.C_{\mathrm{pdm}}\left(\nabla^{-}\nabla^{+}-|\nabla^{-}\nabla^{+}|\right)|\nabla^{-}|\right],\leavevmode\nobreak\
\mathrm{and}$ $\displaystyle\Phi^{r}_{i+1/2}$ $\displaystyle=$
$\displaystyle\Phi^{HO}_{i+1/2}+\mathrm{sign}\left(\nabla^{+}\right)\mathrm{max}\left[0,|\Phi_{i+1}-\Phi^{HO}_{i+1/2}|\right.$
(2) $\displaystyle-$
$\displaystyle\left.C_{\mathrm{pdm}}\left(\nabla^{++}\nabla^{+}-|\nabla^{++}\nabla^{+}|\right)|\nabla^{++}|\right],$

where $\nabla^{+}=\Phi_{i+1}-\Phi_{i}$, $\nabla^{-}=\Phi_{i}-\Phi_{i-1}$,
$\nabla^{++}=\Phi_{i+2}-\Phi_{i+1}$. $C_{\mathrm{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

$\displaystyle f_{i+1/2}$ $\displaystyle=$
$\displaystyle-\frac{c_{i+1/2}}{2}|\mathrm{sign}\left(\Phi^{r}_{i+1/2}-\Phi^{l}_{i+1/2}\right)$
(3) $\displaystyle+$
$\displaystyle\mathrm{sign}\left(\nabla^{+}\right)|\left(\Phi^{r}_{i+1/2}-\Phi^{l}_{i+1/2}\right),$

where $c_{i+1/2}=c_{\mathrm{s}}+v_{\mathrm{A}}+|\mathbf{v}|$ the
characteristic velocity in terms of the sound speed $c_{\mathrm{s}}$, Alfvén
velocity $v_{\mathrm{A}}$ and velocity vector $\mathbf{v}$. Once the fluxes
are calculated, the remainder of the scheme is the same as that described in
Rempel (2017). The diffusion scheme is applied to the logarithm of the density
$\rho$, internal energy per gram $\epsilon$, and population fractions
$\mathbf{n}/\rho$, as well as the velocity components $v_{x},v_{y},v_{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 are 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 is added in the
vertical direction to
$\log\rho,\log\epsilon,\log\mathbf{n}\leavevmode\nobreak\
\mathrm{and}\leavevmode\nobreak\ v_{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 gridpoints which
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 in the simulation is
$C_{\mathrm{pdm}}=2$. These switches include; To reduce the errors in
$\nabla\cdot B$ the diffusive flux of $B_{z}$ in the z-direction is set to
zero at the vertical boundaries and the numerical diffusivity of $\mathbf{B}$
in the direction of $\mathbf{B}$ is 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}$ is limited to a maximum of
$1\leavevmode\nobreak\ \mathrm{km\;s}^{-1}$.

Two hard switches are included to prevent the formation of instabilities which
can occur in the solar atmosphere. Firstly, if the maximum density contrast
between grid-points exceeds 10, then the mass diffusivity is increased
($C_{\mathrm{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 ($\gamma<1.12$) is increased
($C_{\mathrm{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, use 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 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 require an
equation of state (EoS). We combine two approaches; firstly, we use a pre-
tabulated LTE equation of state to model the interior based on the
prescription introduced by Vardya (1965), and extended by Mihalas (1967) and
Wittmann (1974) (VMW), see also Vitas & Khomenko (2015). Secondly, an equation
of state including a non-equilibrium treatment of hydrogen is used for the
solar chromosphere and corona. The pre-tabulated LTE EoS is used for pressure
greater than $2\times 10^{5}\;\mathrm{dyn\;cm}^{-2}$, and an EoS state
including the non-equilibrium treatment of hydrogen is used for lower
pressures. As the non-ideal gas and non-equilibrium ionisation effects are
negligible at this pressure the two methods join smoothly. We use a mixture of
the 15 most abundant elements of the Sun and include the molecules $H_{2}$ in
non-equilibrium, and $H_{2}^{+}$ and $H^{-}$ in chemical equilibrium. We
include up to three ionisation states for all LTE elements and a 5-level plus
continuum model of hydrogen in NE. In this Section we focus on the NE EoS, 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_{\mathrm{int}}$ and density $\rho$ of the plasma. From the
density and the atomic abundances, the total number density of hydrogen nuclei
$n_{\mathrm{H,tot}}$ and non-hydrogen nuclei $n_{\mathrm{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_{\mathrm{e}}$) and population number densities
($\mathbf{n}=\left\\{n_{a,i,j}\right\\}$). 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_{\mathrm{e,nonH}}$) and
energies of excitation and ionisation per hydrogen nuclei
($\epsilon_{\mathrm{nonH}}$) for non-hydrogen atoms are tabulated as a
function of temperature and electron number. These are calculated as;

$n_{\mathrm{\mathrm{e,nonH}}}\left(n_{\mathrm{e}},T\right)=\sum_{a=1}^{15}\sum_{i=0}^{3}i\frac{n_{a,i}}{n_{a,\mathrm{tot}}}\frac{n_{a,\mathrm{tot}}}{n_{\mathrm{H,tot}}},$
(4)
$\epsilon_{\mathrm{nonH}}\left(n_{\mathrm{e}},T\right)=\sum_{a=1}^{15}\sum_{i=0}^{3}\left[\chi_{a,i}+k_{\mathrm{B}}T^{2}\frac{\partial\ln\left(U_{a,i}\right)}{\partial
T}\right]\frac{n_{a,i}}{n_{a,\mathrm{tot}}}\frac{n_{a,\mathrm{tot}}}{n_{\mathrm{H,tot}}},$
(5)

Where $n_{a,\mathrm{tot}}/n_{H,\mathrm{tot}}$ is the number fraction of
element $a$ relative to hydrogen, $\chi_{a,i}$ is the ionisation energy,
$k_{\mathrm{B}}$ the Boltzmann constant, $U$ is the partition function and
$n_{a,i}/n_{a,\mathrm{tot}}$ is the fraction of element $a$ in ionisation
stage $i$ calculated using Saha-Boltzmann (Eqn. 37). 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

$\displaystyle f_{0}$ $\displaystyle=$ $\displaystyle
1-\frac{1}{E_{\mathrm{int}}}\left(\frac{3k_{\mathrm{B}}T}{2}\left[n_{\mathrm{e}}+n_{\mathrm{nonH}}+n_{\mathrm{H}_{2}}+n_{\mathrm{H}_{2},1}+n_{\mathrm{H},-1}\right.\right.$
(6) $\displaystyle+$
$\displaystyle\left.\left.\sum_{i,j}n_{{}_{\mathrm{H},i,j}}\right]+n_{\mathrm{H,tot}}\epsilon_{\mathrm{nonH}}+n_{\mathrm{H}_{2}}E_{\mathrm{H}_{2}}+n_{\mathrm{H}_{2},1}E_{\mathrm{H}_{2},1}\right.$
$\displaystyle+$
$\displaystyle\left.n_{\mathrm{H},-1}E_{\mathrm{H},-1}+\sum_{i,j}n_{{}_{\mathrm{H},i,j}}E_{{}_{\mathrm{H},i,j}}\right)=0,$

where the $E_{a,i,j}$ gives the ionisation, excitation and dissociation
energies of the atom/molecule. The derivatives of
$n_{\mathrm{\mathrm{e,nonH}}}$ and $\epsilon_{\mathrm{nonH}}$ with respect to
$n_{\mathrm{e}}$ and $T$ are calculated numerically from the table
interpolants. The energies of the hydrogen species are given by Eqns. 40-43.
The equation of charge conservation is

$f_{1}=1-\frac{1}{n_{\mathrm{e}}}\left(n_{\mathrm{H},1}+n_{\mathrm{H}_{2},1}-n_{\mathrm{H},-1}+n_{\mathrm{H,tot}}n_{\mathrm{e,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_{\mathrm{H,tot}}}\left(\sum_{i,j}n_{{}_{\mathrm{H},i,j}}+2n_{\mathrm{H}_{2}}+2n_{\mathrm{H}_{2},1}+n_{\mathrm{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 $\mathbf{v}$, the rate of collisional
($C_{ij,kl}$) and radiative ($R_{ij,kl}$) transitions between level
$\left(ij\right)$ and level $\left(kl\neq ij\right)$, 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 form ($r^{-}$) the species $\left(aij\right)$.
Defining $P_{ij,kl}=C_{ij,kl}+R_{ij,kl}$, the rate equation of a species
$n_{aij}$ is described by

$\displaystyle\frac{\partial n_{aij}}{\partial t}$ $\displaystyle=$
$\displaystyle-\nabla\cdot\left(n_{aij}\mathbf{v}\right)$ (9) $\displaystyle+$
$\displaystyle\sum_{kl\neq ij}n_{ak}P_{kl,ij}-n_{aij}\sum_{kl\neq
ij}P_{ij,kl}$ $\displaystyle+$
$\displaystyle\sum_{r^{+}}n_{r^{+},A}n_{r^{+},B}n_{r^{+},C}K_{r^{+}}$
$\displaystyle-$
$\displaystyle\sum_{r^{-}}n_{r^{-},A}n_{r^{-},B}n_{r^{-},C}K_{r^{-}}.$

In this work we treat only hydrogen and $H_{2}$ in non-equilibrium. To solve
Eqn. 9 we solve separately for the continuity and the rates. The advection of
the species by the macroscopic fluid velocity $\mathbf{v}$ is then given by

$\frac{\partial n_{aij}}{\partial
t}=-\nabla\cdot\left(n_{aij}\mathbf{v}\right).$ (10)

We found that advecting the populations with the 4th order central differences
and 4-stage temporal integration scheme used for MHD variables leads to
frequent negative values. Instead, we use a flux-limited unsplit donor cell
method. Each sub-stage of the temporal integration scheme advances the fluid
variables $\rho,\epsilon,\mathbf{v}$ in time, with a timestep $\Delta t$, from
$t_{0}$ to $t_{*}=t_{0}+\Delta t$. To advect the populations we use the
average velocity
$\mathbf{v}^{1/2}=0.5\left(\mathbf{v}^{t_{0}}+\mathbf{v}^{t_{*}}\right)$. For
each direction, the velocities $\mathbf{v}^{1/2}$ are 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 are
calculated using the PDM limiter (Eqn. 1 & 2), as applied in the diffusion
scheme (Sect. 2.1). The flux of the quantity $\Phi$ through the cell interface
$i+1/2$, in direction $d=x,y,z$, is calculated as

$\displaystyle f_{d,i+1/2}$ $\displaystyle=$
$\displaystyle\frac{1}{2}(v^{1/2}_{d,i+1/2}+|v^{1/2}_{d,i+1/2}|)\Phi^{l}_{i+1/2}$
(11) $\displaystyle+$
$\displaystyle\frac{1}{2}(v_{d,i+1/2}^{1/2}-|v_{d,i+1/2}^{1/2}|)\Phi^{r}_{i+1/2},$

and the populations are updated using the sum of the fluxes through all faces
of the cell,

$\Phi^{t_{*}}=\Phi^{t_{0}}+\sum_{d=x,y,z}\frac{\Delta t}{\Delta
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 Eqns. 6-8 ensures consistency between the EoS variables,
populations, and magnetohydrodynamic energy $E_{\mathrm{int}}^{*}$ and density
$\rho^{*}$. After all sub-stages are complete the directionally split
diffusion scheme is run. The populations are diffused using the scheme
described in Section 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,

$\displaystyle f_{3+i,j}$ $\displaystyle=$
$\displaystyle\frac{n_{aij}^{t_{0}+\Delta t}}{n_{aij}^{t_{0}}}-\frac{\Delta
t}{n_{aij}^{t_{0}}}\left(\sum_{kl\neq ij}n_{akl}P_{kl,ij}-n_{aij}\sum_{kl\neq
ij}P_{ij,kl}\right.$ (13) $\displaystyle+$
$\displaystyle\left.\sum_{r^{+}}n_{r^{+},A}n_{r^{+},B}n_{r^{+},C}K_{r^{+}}\right.$
$\displaystyle-$
$\displaystyle\left.\sum_{r^{-}}n_{r^{-},A}n_{r^{-},B}n_{r^{-},C}K_{r^{-}}\right)-1=0,$

where we include 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

$\displaystyle f_{8}$ $\displaystyle=$
$\displaystyle\frac{n_{H_{2}}^{t_{0}+\Delta
t}}{n_{H_{2}}^{t_{0}}}-\frac{\Delta
t}{n_{H_{2}}^{t_{0}}}\left(\sum_{r^{+}}n_{r^{+},A}n_{r^{+},B}n_{r^{+},C}K_{r^{+}}\right.$
(14) $\displaystyle-$
$\displaystyle\left.\sum_{r^{-}}n_{r^{-},A}n_{r^{-},B}n_{r^{-},C}K_{r^{-}}\right)-1=0,$

where the rates for the formation and destruction of molecular hydrogen are
listed in Appendix D. We use 5 out of 6 of the atomic hydrogen rate equations
and the rate equation of molecular hydrogen. We include the nucleus
conservation equation, and discard 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 non-equilibrium 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_{2}$ fraction becomes dominant. It is
therefore necessary to include additional mechanisms to prevent temperatures
becoming too low. We include three mechanisms. Firstly, an additional time-
step constraint is included. The time-step $\Delta t$ is limited such that
$|Q_{\mathrm{rad}}|\Delta t/E_{\mathrm{int}}\leq 0.25$, $Q_{rad}$ is the total
radiative cooling/heating. This damps large decreases in energy due to
radiative cooling. Secondly, a minimum temperature threshold is set. Rather
than including a parameterised heating term, as in Leenaarts et al. (2011), a
temperature floor is implemented. For a minimum temperature
$T_{\mathrm{min}}$, a minimum value $E_{\mathrm{min}}$ for the instantaneous
parts of the internal energy equation is calculated. This includes the kinetic
$nk_{\mathrm{B}}T_{\mathrm{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 is 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_{\mathrm{min}}=2500\leavevmode\nobreak\ \mathrm{K}$. Finally, if the
populations will not converge to the required tolerance, due to the
temperature dropping below a threshold value of $1000\leavevmode\nobreak\
\mathrm{K}$ during a solver call, we allow a small amount of $H^{+}$ to be
recombined to ensure convergence. This allows faster convergence of grid
points which will be limited anyway by the $2500\leavevmode\nobreak\
\mathrm{K}$ floor.

###  2.4 Scattering Multi-group Radiation Transfer

The MURaM code uses a multigroup 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 is extended to include a scattering term. We follow 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 solution of the time-independent radiative transfer
equation

$\frac{dI_{\nu}}{d\tau_{\nu}}=S_{\nu}-I_{\nu},$ (15)

where $\nu$ is the frequency, $S_{\nu}$ the source function, $I_{\nu}$ is the
specific intensity, and $d\tau_{\nu}=\chi_{\nu}ds$ is the optical thickness
over a path length $ds$, where $\chi_{\nu}$ is the plasma opacity. Eqn. 15
must be solved for a number of ray directions $\mathbf{\hat{n}}$. The source
function has been expanded to include a scattering term,

$S_{\nu}=\left(1-\epsilon_{\nu}\right)J_{\nu}+\epsilon B_{\nu},$ (16)

where $\epsilon_{\nu}$ is the photon destruction probability, $B_{\nu}$ the
LTE (Planck) source function, and the mean intensity $J_{\nu}$ is calculated
as the integral over all angles

$J_{\nu}=\frac{1}{4\pi}\int_{4\pi}I_{\nu}\left(\mathbf{\hat{n}}\right)d\omega.$
(17)

The type A quadrature of Carlson (1963), including three points per quadrant,
is used to perform the angular integration. The radiative energy flux
$\mathbf{F_{\nu}}$ is

$\mathbf{F_{\nu}}=\int_{4\pi}I_{\nu}\left(\mathbf{\hat{n}}\right)\mathbf{\hat{n}}d\omega.$
(18)

For the formal solution to the radiative transfer Eqn. 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}}+\Phi_{U}S_{U}+\Phi_{O}S_{O},$ (19)

where $(U)$ is the upwind point, $\tau_{U}$ is the optical distance on the
segment $UO$ and the $\Phi$ quantities are the weights. The formal solution,
Eqns. 16 & 19, can be written as

$J_{\nu}=\Lambda\left[S_{\nu}\right]+\mathcal{J_{\nu}},$ (20)

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

Direct solution of Eqn. 20 is expensive, and simply updating Eqn. 16 with the
new $J_{\nu}$ leads to slow convergence. We employ the Approximate Lambda
Iteration (ALI) method (Cannon, 1973), and the diagonal operator as the
approximate operator. The scheme is iterated until a tolerance of $10^{-4}$ is
reached on the relative correction of the source function. Once the source
function has converged the radiative cooling/heating $Q_{\nu}$ is calculated
from two equivalent expressions

$Q_{RT}=-\int_{\nu}\left(\nabla\cdot\mathbf{F_{\nu}}\right)d\nu=4\pi\int_{\nu}\chi_{\nu}\left(J_{\nu}-S_{\nu}\right)d\nu.$
(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^{3}-10^{5}$ 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 $4-5$ bands (Vögler et al.,
2004). To calculate the radiation field we require three band-integrated
quantities, the extinction coefficient $\chi_{j}$, the mean scattering albedo
$\left(1-\epsilon\right)_{j}$ and the band-integrated emissivity
$\left(\epsilon B\right)_{j}$. The binning process is discussed in detail in
Appendix E. The integral over wavelength in Eqn 21 is then solved as

$\displaystyle Q_{RT}$ $\displaystyle=$
$\displaystyle-\sum_{j}\left[\left(\nabla\cdot\mathbf{F_{j}}\right)\left(1-e^{-\tau_{j}/\tau_{0}}\right)\right.$
(22) $\displaystyle+$
$\displaystyle\left.4\pi\chi_{j}\left(J_{j}-S_{j}\right)e^{-\tau_{j}/\tau_{0}}\right]$

where $\tau_{j}$ is the band mean optical depth, $\tau_{0}=0.1$ and the term
$e^{-\tau_{j}/\tau_{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\approx S$, which are amplified
as $\chi$ grows exponentially with depth (Bruls et al., 1999).

###  2.5 Radiative Cooling/Heating

To perform radiation MHD simulations from the convection zone to the corona a
radiation scheme is required which can accurately model a range of physical
regimes. These include the deep interior, which can be treated using the
diffusion approximation. The photosphere where spectral line formation becomes
important, requiring a 3D multi-group approach (Sect. 2.4). Finally, in the
chromosphere the radiation transport must include NLTE effects, such as PRD
and scattering. Detailed NLTE radiation transfer is too computationally
expensive to be performed in a 3D time-dependent simulation. To include
accurate and fast radiative cooling/heating in the chromosphere and corona, we
use 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 use 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 take 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 continuua are the
Lyman-$\alpha$, H-$\alpha$ and the Lyman continuum of hydrogen, the Mg ii h &
k and the Ca ii H & K lines. These lines are modelled using a simplified
description of the heating/losses, for element $a$ in ionisation stage $i$;

$Q_{X}=-f(T)q_{a,i}(\tau)n_{a,\mathrm{tot}}n_{\mathrm{e}},$ (23)

where $f(T)=L_{a,i}(T)F_{a,i}(T)$. Here three pre-tabulated quantities111The
tables for the chromospheric lines are available as part of an IRIS data
release http://iris.lmsal.com/bf/code.tar.bz2. 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}(\tau)$ is the escape
probability as a function of some optical depth proxy $\tau$. For the escape
probability of calcium and magnesium the column mass is used as a proxy for
$\tau$, and for hydrogen the neutral hydrogen column density is used. The
optically thin radiative losses $Q_{\mathrm{thin}}$ are calculated as;

$Q_{\mathrm{thin}}=-\Lambda(T)n_{\mathrm{H,tot}}n_{\mathrm{e}},$ (24)

where $\Lambda(T)$ is given as a table in terms of temperature.

Additionally, we include the back-heating of chromospheric plasma
($Q_{back}$). This is performed using the 3D radiation transport scheme,
following Carlsson & Leenaarts (2012). The emissivity is given in terms of the
optically thin coronal losses given by:

$\eta_{back}=-\frac{Q_{\mathrm{thin}}}{4\pi},$ (25)

and the opacity at the ionisation edge of helium is used;

$\chi_{\mathrm{back}}=\alpha\frac{n_{\mathrm{He\;I}}}{n_{\mathrm{He}}}\left(T,p_{e}\right)\frac{n_{\mathrm{He}}}{\rho},$
(26)

where $\alpha$ is the opacity at the ionisation edge of helium,
$n_{\mathrm{He\;I}}/n_{\mathrm{He}}$ is the neutral helium fraction and is
pre-tabulated in LTE in terms of electron pressure and temperature, and
$n_{\mathrm{He}}/\rho$ is the number of particles of helium per gram of
stellar material.

The full radiative cooling/heating prescription is then:

$Q_{\mathrm{rad}}=Q_{\mathrm{RT}}+Q_{\mathrm{H}}+Q_{\mathrm{Mg}}+Q_{\mathrm{Ca}}+Q_{\mathrm{thin}}+Q_{\mathrm{back}},$
(27)

where $Q_{\mathrm{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 is switched off based on the band-averaged optical depth
$\tau$. This is performed using a function of the form
$\tau^{2}/\left(\tau^{2}+\tau_{\mathrm{cutoff}}^{2}\right)$. A value of
$\tau_{\mathrm{cutoff}}=1.0\times 10^{-4}$ is used in these simulations. The
multi-group RT scheme groups frequencies into bands. Which band a particular
frequency $\nu$ goes into depends on the height where $\tau_{\nu}=1$. We use
as a reference $\tau_{500}$, the optical depth at $500\leavevmode\nobreak\
\mathrm{nm}$. Previous work has found 4-bands are sufficient to capture back-
heating and line-blanketing in the photosphere and temperature minimum (Vögler
et al., 2005). We use a 4-band setup similar to Carlsson et al. (2016), with
boundaries at the heights where $\tau_{500}=10^{-1/2}$, $10^{-3/2}$, and
$10^{-5/2}$.

###  2.6 MHD Equations

The set of equations solved by the MURaM code, are

$\displaystyle\frac{\partial\rho}{\partial t}$ $\displaystyle=$
$\displaystyle-\nabla\cdot\left(\rho\mathbf{v}\right)$ (28)
$\displaystyle\frac{\partial\rho\mathbf{v}}{\partial t}$ $\displaystyle=$
$\displaystyle-\nabla\cdot\left(\rho\mathbf{v}\mathbf{v}\right)-\nabla
p+\rho\mathbf{g}+\mathbf{F}_{\mathrm{L}}+\mathbf{F}_{\mathrm{SR}}$ (29)
$\displaystyle\frac{\partial E_{\mathrm{HD}}}{\partial t}$ $\displaystyle=$
$\displaystyle-\nabla\cdot\left[\mathbf{v}\left(E_{HD}+p\right)+q\mathbf{\hat{b}}\right]+\rho\mathbf{v}\cdot\mathbf{g}+\mathbf{v}\cdot\mathbf{F_{L}}$
(30) $\displaystyle+$
$\displaystyle\mathbf{v}\cdot\mathbf{F_{\mathrm{SR}}}+Q_{\mathrm{rad}}+Q_{\mathrm{res}}$
$\displaystyle\frac{\partial\mathbf{B}}{\partial t}$ $\displaystyle=$
$\displaystyle\nabla\times\left(\mathbf{v}\times\mathbf{B}\right)$ (31)

where $\rho$ is the plasma density, $\mathbf{v}$ the velocity vector, $p$ the
gas pressure, $\mathbf{g}$ the gravitational acceleration,
$\mathbf{F}_{\mathrm{L}}$ the Lorentz force, $\mathbf{F}_{\mathrm{SR}}$ the
semi-relativistic (Boris) correction, $E_{\mathrm{HD}}$ the hydrodynamic
energy, q the Spitzer heat flux, $\mathbf{\hat{b}}=\mathbf{B}/|\mathbf{B}|$
the unit vector in the direction of the magnetic field vector $\mathbf{B}$,
$Q_{\mathrm{rad}}$ the radiative cooling/heating, $Q_{\mathrm{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 heating and cooling effects are given by
$Q_{\mathrm{rad}}$ and described in detail in Section 2.5 above. The
hydrodynamic energy
$E_{\mathrm{HD}}=E_{\mathrm{int}}+\frac{1}{2}\rho\mathbf{v}^{2}$ is used
instead of the total energy. This prevents numerical errors in calculating
$E_{\mathrm{int}}$ from $E_{\mathrm{tot}}$ in low-$\beta$ 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_{\mathrm{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_{\mathrm{\mathrm{sat}}}\sigma
T^{\frac{5}{2}}\left(\mathbf{\hat{b}}\cdot\nabla\right)T-q\right)\\\ $ (32)

where $\sigma=1.1\times 10^{-6}\leavevmode\nobreak\
\mathrm{erg\;cm^{-1}\;s^{-1}\;K^{-7/2}}$ is the constant of Spitzer heat
conductivity, $f_{\mathrm{sat}}$ controls the saturation of thermal conduction
and $\tau$, which is used to control the transition between parabolic and
hyperbolic solutions to the heat conduction equation, has the form

$\tau=\left(f_{\mathrm{CFL}}\frac{\Delta x_{\mathrm{min}}}{\Delta
t}-|v|\right)^{-2}\frac{f_{\mathrm{sat}}\sigma T^{7/2}}{E_{\mathrm{int}}},$
(33)

where $f_{\mathrm{CFL}}\frac{\Delta x_{\mathrm{min}}}{\Delta t}-|v|$ is used
as a maximum propagation speed in order to avoid violations of the
Courant–Friedrichs–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
$(c=f_{\mathrm{CFL}}\frac{\Delta x_{\mathrm{min}}}{\Delta t})$, in terms of
the minimum spatial grid-scale $x_{\mathrm{min}}$ and time-step $\Delta t$. In
order to explicitly integrate the system of equations we set a lower limit on
$\tau$ of $4\Delta t$.

The treatment of the Lorentz force $\mathbf{F}_{\mathrm{L}}$ follows that of
Rempel (2017),

$\displaystyle\mathbf{F}_{\mathrm{L}}$ $\displaystyle=$
$\displaystyle\frac{f_{A}}{4\pi}\nabla\cdot\left(\mathbf{B}\mathbf{B}-\frac{1}{2}\mathbf{I}\mathbf{B}\right)$
(34) $\displaystyle+$
$\displaystyle\frac{1-f_{A}}{4\pi}\left(\nabla\times\mathbf{B}\right)\times\mathbf{B},$

where $f_{A}=\frac{1}{\sqrt{1+\left(v_{a}/c_{\mathrm{max}}\right)^{4}}}$ is
the Alfvén limit factor in terms of the reduced speed of light
$c_{\mathrm{max}}$ and Alfvén speed $v_{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+(v_{a}/c_{\mathrm{max}})^{2})$. 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 $\mathbf{F}_{\mathrm{SR}}$
given by

$\displaystyle\mathbf{F}_{\mathrm{SR}}$ $\displaystyle=$
$\displaystyle-\left(1-f_{A}\right)\left[\mathcal{I}-\mathbf{\hat{b}}\mathbf{\hat{b}}\right]$
(35)
$\displaystyle\left(-\rho\left(\mathbf{v}\cdot\nabla\right)\mathbf{v}-\nabla
p+\rho\mathbf{g}+\mathbf{F_{L}}+\nabla\cdot\boldsymbol{\tau}^{\mathrm{diff}}\right).$

The components of the numerical viscous stress tensor
($\boldsymbol{\tau}^{\mathrm{diff}}$) are calculated

$\left(\nabla\cdot\boldsymbol{\tau}^{\mathrm{diff}}\right)_{i}=\sum_{j=x,y,z}\frac{1}{\Delta_{j}}\left(f^{j+1/2}_{v_{i}}-f^{j-1/2}_{v_{i}}\right),$
(36)

where $\Delta_{j}$ is the grid resolution, and $f^{j}_{v_{i}}$ are the
diffusive fluxes of velocity component $i$ in the $j$ direction, calculated
using Eqn. 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 has been run in non-equilibrium for 3850 seconds, starting from
the first publicly available snapshot ”285”.

Figure 1: A snapshot of the simulation after 3850s, showing: panel a) The
vertical velocity at the photosphere, b) the normalised intensity at
$500\leavevmode\nobreak\ \mathrm{nm}$, panel c) the vertical magnetic field at
the photosphere, and panel d) the vertical magnetic field at a height of
$2\leavevmode\nobreak\ \mathrm{Mm}$ in the chromosphere. The green dashed
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 $\tau_{500}=1$.

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 is performed using log-linear
interpolations for density, energy, electron and population numbers, and
linear interpolations for velocities and magnetic fields.

A different equation of state 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 calculate an equation of state using the same
abundances. There remain inconsistencies between the two EoS’s, coming from
the non-ideal formulation used by Bifrost, and differences in partition
functions for atoms and molecules. The internal energy (Eqn 6) was
recalculated with the EoS described in Section 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\times 10^{-3}$, with a maximum of $6\times 10^{-2}$ near the transition
region.

The resulting simulation has $504\times 504\times 840$ grid-points, spanning
$24\leavevmode\nobreak\ \mathrm{Mm}$ in the horizontal direction and
$16.8\leavevmode\nobreak\ \mathrm{Mm}$ in the vertical, with the lower
boundary at $-2.44\leavevmode\nobreak\ \mathrm{Mm}$, where 0 is the averaged
solar surface. This corresponds to a horizontal resolution of
$47.6\leavevmode\nobreak\ \mathrm{km}$ and a vertical resolution of
$20\leavevmode\nobreak\ \mathrm{km}$. The diffusion scheme is described in
Section 2.1, and a PDM coefficient of $C_{\mathrm{pdm}}=2$ is used for both
the diffusion scheme and the population advection. The viscous and resistive
heating is determined from the momentum and magnetic field fluxes as described
in (Rempel, 2017). We enforce 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 (OSb) condition described
by Rempel (2014).

We limit 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_{\mathrm{max}}=\max\left(2c_{\mathrm{s},\;\mathrm{max}},3v_{max}\right)$ in
terms of the maximum sound speed ($c_{\mathrm{s},\;\mathrm{max}}$) and
velocity ($v_{max}$) in the simulation domain. In addition, we impose a
dynamic ceiling on velocity ($v_{\mathrm{max}}$) and internal energy
($\epsilon_{\mathrm{max}}$) in order to prevent extreme values, only realised
in a few grid points, dominating the simulation time step. The value of the
ceiling $v_{\mathrm{max}}$ is chosen dynamically so that it affects fewer than
one in one million grid points. If more than $9.4\times 10^{-7}$ of the
simulation grid-points are above 0.95 of the limited $v_{\mathrm{max}}$ (or
$\epsilon_{\mathrm{max}}$), then the value is increased by 1%. If fewer than
$4.7\times 10^{-7}$ points are above 0.95 of the limit then is lowered by 1%.
The simulation presented has a maximum velocity of
$300-400\leavevmode\nobreak\ \mathrm{km\;s}^{-1}$. The maximum speed of light
used in the Boris correction is then calculated
$c_{\mathrm{max}}=\max\left(2c_{\mathrm{s,max}},3v_{\mathrm{max}}\right)$. We
do not allow the speed of light in the box to decrease below
$2000\leavevmode\nobreak\ \mathrm{km\;s}^{-1}$. 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\leavevmode\nobreak\
\mathrm{km\;s}^{-1}$ at the start of the simulation. The simulation was run
until this transient passed and the RMS velocity stabilised, which took about
600 seconds. The additional diffusion and damping is then slowly removed and
the simulation is run for an hour until the corona reheated.

Figure 2: Slices through the chromosphere after 3850s 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 & b) and the upper chromosphere (panels g & h). The middle rows show
vertical slices through a quieter inter-network region (panels c & d) and
through the centre of the network field (panels e & f). Animation available
online. Figure 3: A time distance diagram of temperature spanning 14 minutes
of the simulation. The slice is taken at $y=12\leavevmode\nobreak\
\mathrm{Mm}$, and at two heights $z=1\leavevmode\nobreak\ \mathrm{Mm}$ (panel
a) and $z=2\leavevmode\nobreak\ \mathrm{Mm}$ (panel b).

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 222 Animation Available online,
in the mid chromosphere (panel a & b, $z=1\leavevmode\nobreak\ \mathrm{Mm}$)
is dominated by shocks. The shock-fronts show velocities over
$20\leavevmode\nobreak\ \mathrm{km\;s^{-1}}$ 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\leavevmode\nobreak\ \mathrm{kK}$ in the shock rarefactions to above
$10\leavevmode\nobreak\ \mathrm{kK}$ in the shock fronts. Above the network
fields the temperatures are higher, reaching nearly $100\leavevmode\nobreak\
\mathrm{kK}$ in regions where the transition region is depressed. The
chromospheric velocity field at $2\leavevmode\nobreak\ \mathrm{Mm}$ shows
strong shocks, with velocities above $25\leavevmode\nobreak\
\mathrm{km\;s}^{-1}$. At $3\leavevmode\nobreak\ \mathrm{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 $\approx 5\leavevmode\nobreak\ \mathrm{min}$. In
the network field regions periodic brightenings are seen at a
$x=16\leavevmode\nobreak\ \mathrm{Mm}$ with a period of $\approx
3\leavevmode\nobreak\ \mathrm{min}$. Due to the treatment of helium in LTE, a
large fraction of the plasma at $3\leavevmode\nobreak\ \mathrm{Mm}$ height
plasma sits at around $10\leavevmode\nobreak\ \mathrm{kK}$, the preferred
temperature of the first ionisation stage of helium.

##  4 Radiative Cooling and Heating

Figure 4: Radiative cooling in the photosphere and chromosphere for a slice
through the simulation. Panel a) shows the temperature, b) shows the losses
through the 3D multi-group radiation transport scheme, the chromospheric line
losses are shown for magnesium c), calcium d) and hydrogen e), and panel f)
shows the optically thin losses in the corona. Figure 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.

The radiative losses and heating in the simulation can be split into the
multi-group scheme in and below the temperature minimum and low chromosphere,
the chromospheric line losses, and the optically thin 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 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^{10}\leavevmode\nobreak\
\mathrm{erg\;g}^{-1}\mathrm{s}^{-1}$, are strongly localised near areas of the
transition region where the optically thin losses are high.

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

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
$\tau=E_{\mathrm{inst}}/Q$ using the terms in the internal energy that can
instantaneously change
$E_{\mathrm{inst}}=3/2Nk_{\mathrm{B}}T+E_{\mathrm{nonH}}+E_{\mathrm{mol,CE}}$.
This includes the microscopic kinetic energy, the ionisation of non-hydrogen
species and the formation of $H^{-}$ and $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 $\approx
10\leavevmode\nobreak\ kK$, and hydrogen dominates radiative losses from
$\approx 10kK$ to the transition region, reaching timescales lower than $10$
seconds. These results are similar to those presented in Carlsson & Leenaarts
(2012), with hydrogen dominating the cooling above the mid-chromosphere
($\approx 1.5\leavevmode\nobreak\ \mathrm{Mm}$) and being marginally lower
than magnesium in the low chromosphere.

##  5 Hydrogen populations

To investigate the effects of the non-equilibrium 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 non-equilibrium equation
of state 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 prevent hydrogen ionisation from buffering the
temperature fluctuations. The higher ionisation fraction is also seen in the
electron number density, which is higher than LTE in the low-to-mid
chromosphere.

Figure 7: Histograms of the temperature (top) and electron number (bottom)
throughout the photosphere and chromosphere of the simulation. Comparing the
non-equilibrium values (left columns) to those obtained from the pre-tabulated
LTE EoS (right columns). Figure 8: Properties of the non-equilibrium 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 $b_{T}=T_{\mathrm{NE}}/T_{\mathrm{LTE}}$, c) the NE ionisation
fraction $F_{i}=n_{\mathrm{H},1}/n_{\mathrm{H,tot}}$, d) the number density of
the first excited level of hydrogen $n_{\mathrm{H},0,1}$, and the departure
coefficients of e) the ground state, f) the first excited state, g) protons
and h) molecular hydrogen.

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 non-equilibrium, to the value
calculated using the LTE equation of state,
$b_{X}=X_{\mathrm{NE}}/X_{\mathrm{LTE}}$. The departure coefficient of
temperature ($b_{T}$, panel b) shows up to $35\%$ higher temperature in shocks
and the transition region in the non-equilibrium 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_{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_{2}$ departure
coefficient occurs largely in locations where the temperature departure
coefficient is decreased (or increased). The departure coefficients of the
hydrogen ground state ($b_{1}$, panel e), and first excited state ($b_{2}$,
panel f) in the corona are $b_{1}=1.7\times 10^{4}$ and $b_{2}=4\times
10^{4}$, similar to the values of Leenaarts et al. (2007). Although small
regions with $b_{1}>10^{9}$ are seen, we do not see large regions with
extremely high departure coefficients just above the transition region
($b_{1}>10^{10}$), as observed in Leenaarts et al. (2007). The departure
coefficients of the Bifrost public simulation (Carlsson 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. 11 of Appendix
F.

##  6 Numerical performance

Table 1: Simulation Setups

Simulation Setup | Timing | Computational Cost
---|---|---
RTaa$a$Radiation transfer (RT) calculations are either multi-band or single frequency (grey), and include scattering effects (S) or do not (A). | EoS |  $c_{\mathrm{max}}$bb$b$The minimum value of the dynamically adjusted reduced-speed of light. | RTcc$c$The frequency (in iterations) at which the radiation field is updated. | Back | MHDdd$d$MHD routines include the calculation of the right-hand-side of the MHD equations, the div-B cleaner, and the time integration. | RTee$e$RT routines include the interpolation of the opacity tables, calculation of radiation intensities and calculation of the radiative heating/cooling source term. | Diffff$f$Calculation of the diffusive fluxes. When the non-equilibrium module is turned on, this includes calls to the EoS in between directional sweeps. | EoS | Totalgg$g$Other routines that contribute to the total time include boundary updates, time-step synchronisation and grid exchanges. | dt |  $\mu$s perhh$h$Seconds taken for one core to update one grid-point. | Wall timeii$i$Million CPU hours required to calculate one solar hour.
type | type | ($\mathrm{km}\;\mathrm{s}^{-1}$) | freq. | heating | (s) | (s) | (s) | (s) | (s) | (ms) | update | (Mcore-h)
4-band S | NE | 2000 | 5 | on | 0.996 | 1.486 | 1.250 | 4.697 | 8.490 | 8.89 | 57.30 | 1.38
Grey A | LTE | 2000 | 5 | off | 0.645 | 0.138 | 0.048 | 0.075 | 0.980 | 8.89 | 6.61 | 0.16
Grey A | LTE | 2000 | 5 | on | 0.645 | 0.393 | 0.049 | 0.075 | 1.234 | 8.89 | 8.33 | 0.20
4-band S | LTE | 2000 | 5 | on | 0.648 | 1.494 | 0.048 | 0.075 | 2.341 | 8.89 | 15.80 | 0.38
4-band S | NE | 5000 | 5 | on | 0.997 | 1.162 | 1.565 | 4.483 | 8.270 | 3.75 | 55.81 | 3.18
4-band S | NE | 1000 | 5 | on | 0.991 | 1.619 | 1.190 | 4.616 | 8.476 | 14.0 | 57.20 | 0.87
4-band S | NE | 2000 | 1 | on | 1.024 | 4.52 | 1.206 | 4.675 | 11.50 | 8.92 | 77.61 | 1.86
4-band S | NE | 2000 | 10 | on | 0.995 | 0.913 | 1.210 | 4.588 | 7.786 | 8.89 | 52.55 | 1.261

333

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 use 20 nodes, or 1440 cores. The MURaM code is written with MPI
communication, and does not support hybrid shared memory calculations.

The simulation presented in this work, row 1 of table 1, uses an Alfvén speed
limit of $2000\leavevmode\nobreak\ \mathrm{km\;s^{-1}}$, has a typical time-
step of $8.89\leavevmode\nobreak\ \mathrm{ms}$. We calculate the average time
per iteration from 200 timestep updates, the expected computational time per
grid cell update and the cost of simulating one hour 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\leavevmode\nobreak\ \mu\mathrm{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 perform a number of test simulations. The simplest of these is a
simulation with a LTE equation of state, grey LTE multi-group radiation
transfer, 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\leavevmode\nobreak\ \mu\mathrm{s}$ per grid-point per core). This
LTE, grey simulation spends $14\%$ of the computational time on the radiation
transfer (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 radiation transfer
modules. Including the EUV back heating of the chromosphere, row 3 of table 1,
increases the cost of radiation transfer by almost $300\%$. This corresponds
to a $26\%$ increase of 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 3D
multi-group radiation scheme further increases the computational cost of the
simulation by $90\%$. Including the more realistic treatment of radiation
transfer makes radiation transfer 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 2 or 3 iterations, similar to the grey
radiation bin. However, the optically thin chromospheric lines bin can take up
to 7 iterations to converge.

The greatest computational cost is the inclusion of the non-equilibrium
ionisation of hydrogen in the equation of state. 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 latter 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.6x more computational power per gridpoint than the LTE
coronal simulations presented in Rempel (2017).

In addition, we vary two approximations that have a significant effect on the
computational cost. Firstly, we perform two simulations where the minimum
limit on the speed of light is changed. Increasing the speed of light, to
$5000\leavevmode\nobreak\ \mathrm{km\;s}^{-1}$, 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\leavevmode\nobreak\ \mathrm{km\;s}^{-1}$, row 6 of
table 1, allows a larger time-step and lowers the computational cost by 37% .

A second choice which 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
radiation transfer 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.

##  7 Discussion and Conclusion

The current work improves on the original LTE MURaM code through the
implementation of three main modules; a NE treatment of hydrogen in the EoS,
NLTE tabulated losses in the chromosphere, and a scattering multi-group
radiation transfer 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 equation of state. A detailed comparison between the results from Bifrost
and MuRAM will be the subject of a separate paper.

The simulation shows the importance of a non-equilibrium 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_{1}$ and $b_{2}$ are approximately $10^{3}$ to
$10^{6}$ in the upper chromosphere. The departure populations calculated 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 to 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 non-
equilibrium equation of state are based upon a number of simplifying
assumptions.

  1. 1.

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

  2. 2.

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

  3. 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. 4.

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

  5. 5.

Helium is treated in LTE.

The above approximations are necessary for the simulation of large 3D models
including a non-equilibrium 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 non-equilibrium 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 3D cooling
and heating effects while ignoring scattering and PRD effects. Additionally,
the optically thin losses have been split into 6 EUV frequency bands to
incorporate the Lyman continuum and a simplified 3-level helium atom. These
improvements greatly increase the realism of the method in the upper
chromospheric layers.

Another proposal for the fast non-equilibrium 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 1D plane-parallel
nature may lead to unrealistic variations between neighbouring horizontal
pixels. It is also suggested to decouple the equation of state from atomic
populations, solving the problem 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 equation of state, potentially
offering a significant speed up when solving the non-equilibrium problem.

Finally, a new time-implicit numerical method for solving the detailed NE
radiative MHD problem, built as a 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 (Martínez-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 in non-equilibrium. Including NE ionisation reduces
the impact ambipolar diffusion on the simulation, however a decreased
efficiency of shock heating, and heating of cool low-lying loops are seen
(Martínez-Sykora et al., 2020). Ambipolar diffusion and the Hall effect are
included in the MURaM code (Cheung & Cameron, 2012; Rempel & Przybylski, 2021)
and simulations including both effects will soon be performed.

###### 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. JL 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. LSA, VW, and AIS acknowledge support from the European
Research Council (ERC) under the European Union’s Horizon 2020 research and
innovation program (grant no. 715947).

## References

  * Anusha et al. (2021) Anusha, L. S., van Noort, M., & Cameron, R. H. 2021, ApJ, 911, 71
  * Beeck et al. (2013) Beeck, B., Cameron, R. H., Reiners, A., & Schüssler, M. 2013, A&A, 558, A48
  * Beeck et al. (2012) Beeck, B., Collet, R., Steffen, M., et al. 2012, A&A, 539, A121
  * Beeck et al. (2015) Beeck, B., Schüssler, M., Cameron, R. H., & Reiners, A. 2015, A&A, 581, A42
  * Boris (1970) Boris, J. P. 1970, NRL Memorandum Report 2167
  * Bruls et al. (1999) Bruls, J. H. M. J., Vollmöller, P., & Schüssler, M. 1999, A&A, 348, 233
  * Cannon (1973) Cannon, C. J. 1973, ApJ, 185, 621
  * Cardona et al. (2005) Cardona, O., Simonneau, E., & Crivellari, L. 2005, Revista Mexicana de Fisica, 51, 476
  * Carlson (1963) Carlson, B. G. 1963, in Methods in Computational Physics, Vol. 1, ed. B. Alder & S. Fernbach, 1
  * Carlsson et al. (2016) Carlsson, M., Hansteen, V. H., Gudiksen, B. V., Leenaarts, J., & De Pontieu, B. 2016, A&A, 585, A4
  * Carlsson & Leenaarts (2012) Carlsson, M. & Leenaarts, J. 2012, A&A, 539, A39
  * Carlsson & Stein (1992) Carlsson, M. & Stein, R. F. 1992, ApJ, 397, L59
  * Carlsson & Stein (2002) Carlsson, M. & Stein, R. F. 2002, ApJ, 572, 626
  * Chen et al. (2017) Chen, F., Rempel, M., & Fan, Y. 2017, ApJ, 846, 149
  * Cheung & Cameron (2012) Cheung, M. C. M. & Cameron, R. H. 2012, ApJ, 750, 6
  * Cheung et al. (2019) Cheung, M. C. M., Rempel, M., Chintzoglou, G., et al. 2019, Nature Astronomy, 3, 160
  * Cheung et al. (2007) Cheung, M. C. M., Schüssler, M., & Moreno-Insertis, F. 2007, A&A, 467, 703
  * Collet et al. (2011) Collet, R., Hayek, W., Asplund, M., et al. 2011, A&A, 528, A32
  * Danilovic (2017) Danilovic, S. 2017, A&A, 601, A122
  * de la Cruz Rodríguez & Piskunov (2013) de la Cruz Rodríguez, J. & Piskunov, N. 2013, ApJ, 764, 33
  * Dedner et al. (2002) Dedner, A., Kemm, F., Kröner, D., et al. 2002, Journal of Computational Physics, 175, 645
  * Eggleton et al. (1973) Eggleton, P. P., Faulkner, J., & Flannery, B. P. 1973, A&A, 23, 325
  * Forrey (2013) Forrey, R. C. 2013, ApJ, 773, L25
  * Galli & Palla (1998) Galli, D. & Palla, F. 1998, A&A, 335, 403
  * Golding et al. (2016) Golding, T. P., Leenaarts, J., & Carlsson, M. 2016, ApJ, 817, 125
  * Gombosi et al. (2002) Gombosi, T. I., Tóth, G., De Zeeuw, D. L., et al. 2002, Journal of Computational Physics, 177, 176
  * Gudiksen et al. (2011) Gudiksen, B. V., Carlsson, M., Hansteen, V. H., et al. 2011, A&A, 531, A154
  * Gustafsson et al. (1975) Gustafsson, B., Bell, R. A., Eriksson, K., & Nordlund, A. 1975, A&A, 500, 67
  * Halenka & Madej (2002) Halenka, J. & Madej, J. 2002, Acta Astron., 52, 195
  * Halenka et al. (2001) Halenka, J., Madej, J., Langer, K., & Mamok, A. 2001, Acta Astron., 51, 347
  * Hayek et al. (2010) Hayek, W., Asplund, M., Carlsson, M., et al. 2010, A&A, 517, A49
  * Hirasawa (1969) Hirasawa, T. 1969, Progress of Theoretical Physics, 42, 523
  * Jameson (2017) Jameson, A. 2017, AIAA Journal, 55, 1487
  * Judge (2005) Judge, P. G. 2005, J. Quant. Spec. Radiat. Transf., 92, 479
  * Judge (2017) Judge, P. G. 2017, ApJ, 851, 5
  * Kurucz (1970) Kurucz, R. L. 1970, SAO Special Report, 309
  * Leenaarts et al. (2011) Leenaarts, J., Carlsson, M., Hansteen, V., & Gudiksen, B. V. 2011, A&A, 530, A124
  * Leenaarts et al. (2007) Leenaarts, J., Carlsson, M., Hansteen, V., & Rutten, R. J. 2007, A&A, 473, 625
  * Leenaarts & Wedemeyer-Böhm (2006) Leenaarts, J. & Wedemeyer-Böhm, S. 2006, A&A, 460, 301
  * Ludwig (1992) Ludwig, H. 1992, PhD thesis
  * Martínez-Sykora et al. (2012) Martínez-Sykora, J., De Pontieu, B., & Hansteen, V. 2012, ApJ, 753, 161
  * Martínez-Sykora et al. (2020) Martínez-Sykora, J., Leenaarts, J., De Pontieu, B., et al. 2020, ApJ, 889, 95
  * McElroy et al. (2013) McElroy, D., Walsh, C., Markwick, A. J., et al. 2013, A&A, 550, A36
  * Mihalas (1967) Mihalas, D. 1967, Methods in Computational Physics, 7, 1
  * Mihalas et al. (1988) Mihalas, D., Dappen, W., & Hummer, D. G. 1988, ApJ, 331, 815
  * Nordlund (1982) Nordlund, A. 1982, A&A, 107, 1
  * Palla et al. (1983) Palla, F., Salpeter, E. E., & Stahler, S. W. 1983, ApJ, 271, 632
  * Panja et al. (2020) Panja, M., Cameron, R., & Solanki, S. K. 2020, ApJ, 893, 113
  * Pereira & Uitenbroek (2015) Pereira, T. M. D. & Uitenbroek, H. 2015, A&A, 574, A3
  * Popovas & Jørgensen (2016) Popovas, A. & Jørgensen, U. G. 2016, A&A, 595, A130
  * Rempel (2014) Rempel, M. 2014, ApJ, 789, 132
  * Rempel (2017) Rempel, M. 2017, ApJ, 834, 10
  * Rempel & Przybylski (2021) Rempel, M. & Przybylski, D. 2021, ApJ, 923, 79
  * Rempel et al. (2009) Rempel, M., Schüssler, M., Cameron, R. H., & Knölker, M. 2009, Science, 325, 171
  * Rogers et al. (1996) Rogers, F. J., Swenson, F. J., & Iglesias, C. A. 1996, ApJ, 456, 902
  * Schüssler et al. (2003) Schüssler, M., Shelyag, S., Berdyugina, S., Vögler, A., & Solanki, S. K. 2003, ApJ, 597, L173
  * Schüssler & Vögler (2006) Schüssler, M. & Vögler, A. 2006, ApJ, 641, L73
  * Shapiro et al. (2017) Shapiro, A. I., Solanki, S. K., Krivova, N. A., et al. 2017, Nature Astronomy, 1, 612
  * Shelyag et al. (2016) Shelyag, S., Khomenko, E., de Vicente, A., & Przybylski, D. 2016, ApJ, 819, L11
  * Shelyag et al. (2007) Shelyag, S., Schüssler, M., Solanki, S. K., & Vögler, A. 2007, A&A, 469, 731
  * Skartlien (2000) Skartlien, R. 2000, ApJ, 536, 465
  * Sollum (1999) Sollum, E. 1999, Master’s thesis
  * Stancil (1994) Stancil, P. C. 1994, ApJ, 430, 360
  * Sukhorukov & Leenaarts (2017) Sukhorukov, A. V. & Leenaarts, J. 2017, A&A, 597, A46
  * Uitenbroek (2001) Uitenbroek, H. 2001, ApJ, 557, 389
  * van Regemorter (1962) van Regemorter, H. 1962, ApJ, 136, 906
  * Vardya (1965) Vardya, M. S. 1965, MNRAS, 129, 205
  * Vitas & Khomenko (2015) Vitas, N. & Khomenko, E. 2015, Annales Geophysicae, 33, 703
  * Vögler et al. (2004) Vögler, A., Bruls, J. H. M. J., & Schüssler, M. 2004, A&A, 421, 741
  * Vögler & Schüssler (2007) Vögler, A. & Schüssler, M. 2007, A&A, 465, L43
  * Vögler et al. (2005) Vögler, A., Shelyag, S., Schüssler, M., et al. 2005, A&A, 429, 335
  * Warnecke & Bingert (2020) Warnecke, J. & Bingert, S. 2020, Geophysical and Astrophysical Fluid Dynamics, 114, 261
  * Wittmann (1974) Wittmann, A. 1974, Sol. Phys., 35, 11
  * Witzke et al. (2021) Witzke, V., Shapiro, A. I., Cernetic, M., et al. 2021, A&A, 653, A65
  * Yeo et al. (2017) Yeo, K. L., Solanki, S. K., Norris, C. M., et al. 2017, Phys. Rev. Lett., 119, 9.1102
  * Zhang et al. (2019) Zhang, B., Sorathia, K. A., Lyon, J. G., et al. 2019, ApJS, 244, 20

##  Appendix A LTE Equation of State

The ideal EoS is pre-tabulated to calculate the thermodynamic variables
$n_{\mathrm{e}}$, $T$, and $p$ in terms of density $\rho$ and internal energy
$\epsilon$. To start, we calculate the LTE ionisation fractions for a given
density $\rho$, electron number $n_{\mathrm{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_{\mathrm{e}}}\left(\frac{2\pi
m_{e}k_{\mathrm{B}}T}{h^{2}}\right)^{1.5}\frac{U_{a,i+1}}{U_{a,i}}e^{-\chi_{a,i}/\left(k_{\mathrm{B}}T\right)+\left(Z_{a,i}+1\right)\Delta\mu},$
(37)

where $k_{\mathrm{B}}$ is Boltzmann’s constant, $m_{e}$ the electron mass, h
is Planck’s constant and $\Delta\mu$ the ionisation lowering due to
interactions with surrounding particles. $U_{a,i}$ is the partition function,
$\chi_{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 use 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) are used. The H- fraction is calculated from Eqn. 37 and the H2 and
H${}_{2}^{+}$ fractions are

$\frac{n_{H,0}n_{H,0}}{n_{\mathrm{H}_{2,0}}}=\left(\frac{2\pi
mk_{\mathrm{B}}T}{h^{2}}\right)^{1.5}\frac{U_{H,0}^{2}}{U_{\mathrm{H}_{2},0}}e^{-D_{\mathrm{H}_{2},0}/\left(k_{\mathrm{B}}T\right)},\leavevmode\nobreak\
\mathrm{and}$ (38)
$\frac{n_{H,0}n_{\mathrm{H},1}}{n_{H_{2},1}}=\left(\frac{2\pi
mk_{\mathrm{B}}T}{h^{2}}\right)^{1.5}\frac{U_{H,0}}{U_{\mathrm{H}_{2},1}}e^{-D_{\mathrm{H}_{2,1}}/\left(k_{\mathrm{B}}T\right)},$
(39)

where $m$ is the reduced mass and
$D_{\mathrm{H}_{2}}=4.478007\leavevmode\nobreak\ \mathrm{ev}$ and
$D_{\mathrm{H}_{2}^{+}}=2.650639\leavevmode\nobreak\ \mathrm{ev}$ are the
dissociation energies. The partition functions of $\mathrm{H}_{2}$ are taken
from (Popovas & Jørgensen 2016) and the partition functions of
$\mathrm{H}_{2}^{+}$ are from a polynomial fit to the table of (Stancil 1994).
Following Mihalas et al. (1988), the energy of the hydrogen species are

$\displaystyle E_{\mathrm{H}_{2}}$ $\displaystyle=$ $\displaystyle
k_{\mathrm{B}}T^{2}\frac{\partial\ln U_{\mathrm{H}_{2},0}}{\partial T},$ (40)
$\displaystyle E_{\mathrm{H}_{2},1}$ $\displaystyle=$ $\displaystyle
D_{\mathrm{H}_{2},0}-D_{\mathrm{H}_{2},1}+\chi_{\mathrm{H},1}+k_{\mathrm{B}}T^{2}\frac{\partial\ln
U_{\mathrm{H}_{2},1}}{\partial T},$ (41) $\displaystyle E_{\mathrm{H},-1}$
$\displaystyle=$ $\displaystyle
0.5D_{\mathrm{H}_{2}}-\chi_{\mathrm{H},-1},\leavevmode\nobreak\ \mathrm{and}$
(42) $\displaystyle E_{\mathrm{H},i,e}$ $\displaystyle=$ $\displaystyle
0.5D_{\mathrm{H}_{2}}+\chi_{\mathrm{H},i,e},$ (43)

where $\chi_{\mathrm{H},i,e}$ is the excitation or ionisation energy of the
hydrogen level $(i,e)$. The ionisation energy of neutral hydrogen is
$\chi_{\mathrm{H},1}=13.59844\leavevmode\nobreak\ \mathrm{ev}$, and for
$\mathrm{H}^{-}$ is $\chi_{\mathrm{H},-1}=0.754\leavevmode\nobreak\
\mathrm{ev}$. Once the fractions of all elements and of the hydrogen molecules
are determined we follow the formulation of the VMW equation of state. This
process involves iterating the ionisation and molecular fractions, electron
number density and temperature until a convergence criteria is reached. We
iterate this procedure until the electron number density converges to a
tolerance of $1.0\times 10^{-8}$. Once converged, the energy and pressure are
calculated

$E_{\mathrm{int}}=E_{\mathrm{exi}}+E_{\mathrm{trans}}+E_{\mathrm{rad}}+E_{\mathrm{C}}+E_{\mathrm{PI}},\leavevmode\nobreak\
\mathrm{and}$ (44)
$p=p_{\mathrm{trans}}+p_{\mathrm{rad}}+p_{\mathrm{C}}+p_{\mathrm{PI}},$ (45)

in terms of excitation and ionisation $(\mathrm{exi})$, translational
$(\mathrm{trans})$, radiation $(\mathrm{rad})$, Coloumb $(\mathrm{C})$ and
pressure ionisation $(\mathrm{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_{\mathrm{trans}}=\left(\sum_{a,i}n_{a,i}+n_{\mathrm{e}}\right)k_{\mathrm{B}}T$
and $E_{\mathrm{trans}}=3/2\;p_{\mathrm{trans}}$.

The contribution to the energy from excitation and ionisation is calculated as

$E_{\mathrm{exi}}=\sum_{a,i}\left(\chi_{a,i}+k_{\mathrm{B}}T^{2}\frac{\partial\ln
U_{a,i}}{\partial T}\right)n_{a,i}.$ (46)

For the Coulomb correction we follow 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_{\mathrm{C}}=-\frac{2\pi^{1/2}e^{3}}{3k_{\mathrm{B}}^{1/2}}\frac{1}{\left(VT\right)^{1/2}}\left(\sum_{a,i}N_{a,i}Z^{2}_{a,i}+N_{\mathrm{e}}\right)^{3/2}\tau\left(x\right),$
(47)

where the function $\tau$ is

$\tau\left(x\right)=3x^{-3}\left(\ln\left(1+x\right)-x+\frac{1}{2}x^{2}\right),$
(48)

and $x$ is given by

$\displaystyle x$ $\displaystyle=$
$\displaystyle\frac{4\pi^{1/2}e^{3}}{3k_{\mathrm{B}}^{3/2}}\frac{1}{V^{1/2}T^{3/2}}$
(49)
$\displaystyle\frac{\left(\sum_{a,i}N_{a,i}Z_{a,i}\right)}{\left(\sum_{a,i}N_{a,i}\right)}\left(\sum_{a,i}N_{a,i}Z^{2}_{a,i}+N_{\mathrm{e}}\right)^{1/2}.$

The required derivatives of $\tau$ and $x$ are

$\frac{\partial\tau}{\partial x}=-\frac{3}{x}\left(\tau-\frac{1}{1+x}\right),$
(50) $\frac{\partial x}{\partial
N_{\mathrm{e}}}=x\frac{1}{2\left(\sum_{a,i}N_{a,i}Z^{2}_{a,i}+N_{\mathrm{e}}\right)},$
(51) $\frac{\partial x}{\partial V}=-x\frac{1}{2V},\mathrm{and}$ (52)
$\frac{\partial x}{\partial T}=-x\frac{3}{2T}.$ (53)

The reduction of the ionisation potential $\Delta\mu_{\mathrm{C}}$, used in
the Saha-Boltzmann equation (37), is

$\displaystyle\Delta\mu_{\mathrm{C}}$ $\displaystyle=$
$\displaystyle-\frac{1}{k_{\mathrm{B}}T}\frac{\partial
F_{\mathrm{C}}}{\partial N_{\mathrm{e}}}$ $\displaystyle=$
$\displaystyle-\frac{F_{\mathrm{C}}}{k_{\mathrm{B}}T}\left(\frac{3}{2}\frac{1}{\sum_{a,i}N_{a,i}Z^{2}_{a,i}+N_{\mathrm{e}}}+\frac{1}{\tau}\frac{\partial\tau}{\partial
x}\frac{\partial x}{\partial N_{\mathrm{e}}}\right),$

giving a pressure correction;

$P_{\mathrm{C}}=-\frac{\partial F_{\mathrm{C}}}{\partial
V}=F_{\mathrm{C}}\left(\frac{1}{2V}-\frac{1}{\tau}\frac{\partial\tau}{\partial
x}\frac{\partial x}{\partial V}\right),$ (55)

and the internal energy correction;

$\epsilon_{\mathrm{C}}=-T^{2}\frac{\partial}{\partial
T}\left(\frac{F_{\mathrm{C}}}{T}\right)=F_{\mathrm{C}}\left(\frac{3}{2}-\frac{1}{\tau}\frac{\partial\tau}{\partial
x}T\frac{\partial x}{\partial T}\right).$ (56)

Simplifying Eqns. 55-56, and taking $V=1.0\;\mathrm{cm}^{-3}$ gives

$\displaystyle\Delta\mu_{C}$ $\displaystyle=$
$\displaystyle\frac{1}{2k_{\mathrm{B}}T}(3\tau+\frac{\partial\tau}{\partial
x}x)$ (57)
$\displaystyle\frac{2\pi^{1/2}e^{3}}{3k_{\mathrm{B}}^{1/2}}\frac{1}{T^{1/2}}\left(\sum_{a,i}n_{a,i}Z^{2}_{a,i}+n_{e}\right)^{1/2},$
$\displaystyle P_{C}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}(\tau+\frac{\partial\tau}{\partial x}x)$ (58)
$\displaystyle\frac{2\pi^{1/2}e^{3}}{3k_{\mathrm{B}}^{1/2}}\frac{1}{T^{1/2}}\left(\sum_{a,i}n_{a,i}Z^{2}_{a,i}+n_{\mathrm{e}}\right)^{3/2},\mathrm{and}$
$\displaystyle\epsilon_{C}$ $\displaystyle=$
$\displaystyle-\frac{3}{2}\left(\tau+\frac{\partial\tau}{\partial x}x\right)$
(59)
$\displaystyle\frac{2\pi^{1/2}e^{3}}{3k_{\mathrm{B}}^{1/2}}\frac{1}{T^{1/2}}\left(\sum_{a,i}n_{a,i}Z^{2}_{a,i}+n_{\mathrm{e}}\right)^{3/2}.$

Additionally, the pressure ionisation device described in the Eggleton,
Faulkner and Flannery (EFF) EoS (Eggleton et al. 1973) is included. This
method provides thermodynamically consistent result that gives a qualitatively
correct pressure ionisation as density increases.

$F_{PI}=\frac{\Omega\left(T\right)}{V}\left(N_{\mathrm{e},0}^{2}-N_{\mathrm{e}}^{2}\right),$
(60)

where $N_{e,0}$ is the electron number of the gas when it is fully ionised,
and $\Omega$ is

$\Omega=\frac{a_{0}^{3}}{2}\left(k_{\mathrm{B}}T+20\chi_{0}\right),$ (61)

where $a_{0}=5.23e^{-9}/\left<Z\right>\leavevmode\nobreak\ \mathrm{cm}$ and
$\chi_{0}=2.16e^{-11}\left<Z\right>^{2}\leavevmode\nobreak\ \mathrm{erg}$ and
$\left<Z\right>$ is the mean charge per nucleus, and the required derivative

$\frac{\partial\Omega}{\partial T}=\frac{a_{0}^{3}}{2}k_{\mathrm{B}}.$ (62)

The change in the potential can then be calculated,

$\Delta\mu_{PI}=-\frac{1}{k_{\mathrm{B}}T}\frac{\partial F_{PI}}{\partial
N_{\mathrm{e}}}=\frac{2\Omega}{k_{\mathrm{B}}T}\left(T\right)n_{\mathrm{e}},$
(63)

the pressure correction;

$P_{PI}=-\frac{\partial F_{PI}}{\partial
V}=\Omega\left(T\right)\left(n_{e,0}^{2}-n_{\mathrm{e}}^{2}\right),$ (64)

and finally the internal energy correction;

$\epsilon_{PI}=-T^{2}\frac{\partial}{\partial
T}\left(\frac{F_{PI}}{T}\right)=\left(\Omega\left(T\right)-VT\frac{\partial\Omega\left(T\right)}{\partial
T}\right)\left(n_{e,0}^{2}-n_{\mathrm{e}}^{2}\right).$ (65)

Once a complete solution is obtained for all required density and energy
values, the entropy $s$ is calculated by integrating over the table

$s=\int_{\epsilon}\int_{\rho}\left(\frac{\epsilon}{T}\partial\ln\epsilon-\frac{p}{T\rho}\partial\ln\rho\right).$
(66)

##  Appendix B Derivatives of the Non-Equilibrium Equation of State

EoS derivatives for chemical equilibrium of $\mathrm{H}_{2}^{+}$;

$\displaystyle\frac{\partial n_{\mathrm{H}_{2,1}}}{\partial T}$
$\displaystyle=$ $\displaystyle-
n_{\mathrm{H}_{2,1}}\left(\frac{3}{2T}+\frac{D_{\mathrm{H}_{2,1}}}{k_{\mathrm{B}}T^{2}}\right.$
(67) $\displaystyle+$
$\displaystyle\left.\frac{1}{U_{\mathrm{H},0}}\frac{\partial U_{H,0}}{\partial
T}-\frac{1}{U_{\mathrm{H}_{2,1}}}\frac{\partial U_{\mathrm{H}_{2,1}}}{\partial
T}\right),$ $\displaystyle\frac{\partial n_{\mathrm{H}_{2,1}}}{\partial
n_{\mathrm{H},0,j}}$ $\displaystyle=$
$\displaystyle\frac{n_{\mathrm{H}_{2,1}}}{\sum_{j}n_{\mathrm{H},0,j}},\leavevmode\nobreak\
\mathrm{and}$ (68) $\displaystyle\frac{\partial n_{\mathrm{H}_{2,1}}}{\partial
n_{\mathrm{H},1}}$ $\displaystyle=$
$\displaystyle\frac{n_{\mathrm{H}_{2,1}}}{n_{\mathrm{H},1}},$ (69)

And for $\mathrm{H}^{-1}$;

$\displaystyle\frac{\partial n_{{}_{\mathrm{H},-1}}}{\partial T}$
$\displaystyle=$ $\displaystyle-
n_{{}_{\mathrm{H},-1}}\left(\frac{3}{2T}+\frac{\chi_{{}_{\mathrm{H},-1}}}{k_{\mathrm{B}}T^{2}}+\frac{1}{U_{H,0}}\frac{\partial
U_{\mathrm{H},0}}{\partial T}\right),$ (70) $\displaystyle\frac{\partial
n_{{}_{\mathrm{H},-1}}}{\partial n_{\mathrm{e}}}$ $\displaystyle=$
$\displaystyle\frac{n_{{}_{\mathrm{H},-1}}}{n_{\mathrm{e}}},\leavevmode\nobreak\
\mathrm{and}$ (71) $\displaystyle\frac{\partial
n_{{}_{\mathrm{H},-1}}}{\partial n_{\mathrm{H},0,j}}$ $\displaystyle=$
$\displaystyle\frac{n_{{}_{\mathrm{H},-1}}}{\sum_{j}n_{\mathrm{H},0,j}}.$ (72)

The derivatives for the energy conservation equation $f_{0}$ are

$\displaystyle\frac{\partial f_{0}}{\partial T}$ $\displaystyle=$
$\displaystyle-\frac{1}{E_{\mathrm{int}}}\left(\frac{3k_{\mathrm{B}}}{2}\left[n_{\mathrm{e}}+n_{\mathrm{nonH}}+n_{\mathrm{H}_{2}}+n_{\mathrm{H}_{2,1}}\right.\right.$
(73) $\displaystyle+$
$\displaystyle\left.n_{{}_{\mathrm{H},-1}}+\sum_{i,j}n_{{}_{\mathrm{H},i,j}}\right]+\frac{3k_{\mathrm{B}}T}{2}\left[\frac{\partial
n_{noH}}{\partial T}+\frac{\partial n_{{}_{\mathrm{H},-1}}}{\partial
T}\right.$ $\displaystyle+$ $\displaystyle\left.\frac{\partial
n_{{}_{\mathrm{H_{2}},1}}}{\partial T}\right]+n_{\mathrm{H,tot}}\frac{\partial
E_{\mathrm{nonH}}}{\partial T}+\frac{\partial E_{H_{2}}}{\partial
T}n_{\mathrm{H}_{2}}$ $\displaystyle+$ $\displaystyle\left.\frac{\partial
E_{{}_{\mathrm{H},-1}}n_{{}_{\mathrm{H},-1}}}{\partial T}+\frac{\partial
E_{H_{2,1}}n_{\mathrm{H}_{2,1}}}{\partial T}\right),$
$\displaystyle\frac{\partial f_{0}}{\partial n_{\mathrm{e}}}$ $\displaystyle=$
$\displaystyle-\frac{1}{E_{\mathrm{int}}}\left(\frac{3k_{\mathrm{B}}T}{2}\left[1.0+\frac{\partial
n_{{}_{\mathrm{H},-1}}}{\partial
n_{\mathrm{e}}}\right]+n_{\mathrm{H,tot}}\frac{\partial
E_{\mathrm{nonH}}}{\partial n_{\mathrm{e}}}\right.$ (74) $\displaystyle+$
$\displaystyle\left.\frac{\partial n_{{}_{\mathrm{H},-1}}}{\partial
n_{\mathrm{e}}}E_{{}_{\mathrm{H},-1}}\right),$ $\displaystyle\frac{\partial
f_{0}}{\partial n_{\mathrm{H},0,j}}$ $\displaystyle=$
$\displaystyle-\frac{1}{E_{\mathrm{int}}}\left(\frac{3k_{\mathrm{B}}T}{2}+E_{\mathrm{H},0,j}\right.$
(75) $\displaystyle+$ $\displaystyle\left.\frac{\partial
n_{\mathrm{H}_{2,1}}}{\partial n_{\mathrm{H},0,j}}E_{H_{2,1}}+\frac{\partial
n_{{}_{\mathrm{H},-1}}}{\partial
n_{\mathrm{H},0,j}}E_{{}_{\mathrm{H},-1}}\right),\leavevmode\nobreak\
\mathrm{and}$ $\displaystyle\frac{\partial f_{0}}{\partial n_{\mathrm{H},1}}$
$\displaystyle=$
$\displaystyle-\frac{1}{E_{\mathrm{int}}}\left(\frac{3k_{\mathrm{B}}T}{2}+E_{H,1}+\frac{\partial
n_{\mathrm{H}_{2,1}}}{\partial n_{\mathrm{H},1}}E_{H_{2,1}}\right)$ (76)

The derivatives for the charge conservation equation $f_{1}$ are

$\displaystyle\frac{\partial f_{1}}{\partial T}$ $\displaystyle=$
$\displaystyle-\frac{1}{n_{\mathrm{e}}}\left(-\frac{\partial
n_{{}_{\mathrm{H},-1}}}{\partial T}+n_{\mathrm{H,tot}}\frac{\partial
n_{\mathrm{\mathrm{e,nonH}}}}{\partial T}\right),$ (77)
$\displaystyle\frac{\partial f_{1}}{\partial n_{\mathrm{e}}}$ $\displaystyle=$
$\displaystyle-\frac{1}{n_{\mathrm{e}}}\left(1-\frac{\partial
n_{{}_{\mathrm{H},-1}}}{\partial
n_{\mathrm{e}}}+n_{\mathrm{H,tot}}\frac{\partial
n_{\mathrm{\mathrm{e,nonH}}}}{\partial n_{\mathrm{e}}}\right)$
$\displaystyle+$
$\displaystyle\frac{1}{n_{\mathrm{e}}^{2}}\left(n_{\mathrm{H},1}+n_{\mathrm{H}_{2,1}}-n_{{}_{\mathrm{H},-1}}+n_{\mathrm{H,tot}}n_{\mathrm{\mathrm{e,nonH}}}\right)$
$\displaystyle\frac{\partial f_{1}}{\partial n_{\mathrm{H},0,j}}$
$\displaystyle=$ $\displaystyle-\frac{1}{n_{\mathrm{e}}}\left(\frac{\partial
n_{\mathrm{H}_{2,1}}}{\partial n_{\mathrm{H},0,j}}+\frac{\partial
n_{{}_{\mathrm{H},-1}}}{\partial n_{\mathrm{H},0,j}}\right),$ (79)
$\displaystyle\frac{\partial f_{1}}{\partial n_{\mathrm{H},1}}$
$\displaystyle=$ $\displaystyle-\frac{1}{n_{\mathrm{e}}}\left(1-\frac{\partial
n_{\mathrm{H}_{2,1}}}{\partial n_{\mathrm{H},1}}\right).$ (80)

The derivatives for the hydrogen nucleus conservation equation $f_{2}$ are

$\displaystyle\frac{\partial f_{2}}{\partial T}$ $\displaystyle=$
$\displaystyle-\frac{1}{n_{H,tot}}\left(1+\frac{\partial
n_{\mathrm{H}_{2,1}}}{\partial T}+\frac{\partial
n_{{}_{\mathrm{H},-1}}}{\partial T}\right),$ (81) $\displaystyle\frac{\partial
f_{2}}{\partial n_{\mathrm{e}}}$ $\displaystyle=$
$\displaystyle-\frac{1}{n_{H,tot}}\frac{\partial
n_{{}_{\mathrm{H},-1}}}{\partial n_{\mathrm{e}}},$ (82)
$\displaystyle\frac{\partial f_{2}}{\partial n_{\mathrm{H},0,j}}$
$\displaystyle=$ $\displaystyle-\frac{1}{n_{H,tot}}\left(1+\frac{\partial
n_{\mathrm{H}_{2,1}}}{\partial n_{\mathrm{H},0,j}}+\frac{\partial
n_{{}_{\mathrm{H},-1}}}{\partial n_{\mathrm{H},0,j}}\right),$ (83)
$\displaystyle\frac{\partial f_{2}}{\partial n_{\mathrm{H},1}}$
$\displaystyle=$ $\displaystyle-\frac{1}{n_{H,tot}}\frac{\partial
n_{\mathrm{H}_{2,1}}}{\partial n_{\mathrm{H},1}}.$ (84)

The derivatives of the rate-equations $f_{2+i}$ are

$\displaystyle f_{3+ij}$ $\displaystyle=$ $\displaystyle\frac{\Delta
t}{n_{aij}^{t_{0}}}\left(\sum_{kl\neq ij}n_{akl}\frac{\partial
P_{kl,ij}}{\partial T}-n_{aij}\sum_{kl\neq ij}\frac{\partial
P_{ij,kl}}{\partial T}\right.$ (85) $\displaystyle+$
$\displaystyle\left.\sum_{r^{+}}n_{r^{+},A}n_{r^{+},B}n_{r^{+},C}\frac{\partial
K_{r^{+}}}{\partial T}\right.$ $\displaystyle-$
$\displaystyle\left.\sum_{r^{-}}n_{r^{-},A}n_{r^{-},B}n_{r^{-},C}\frac{\partial
K_{r^{-}}}{\partial T}\right),$ $\displaystyle\frac{\partial
f_{3+ij}}{\partial n_{\mathrm{e}}}$ $\displaystyle=$
$\displaystyle-\frac{\Delta t}{n_{aij}^{o}}\left(\sum_{kl\neq
ij}n_{akl}\frac{\partial P_{kl,ij}}{\partial n_{\mathrm{e}}}\right.$ (86)
$\displaystyle-$ $\displaystyle\left.n_{aij}\sum_{kl\neq ij}\frac{\partial
P_{ij,kl}}{\partial n_{\mathrm{e}}}\right),$ $\displaystyle\frac{\partial
f_{3+ij}}{\partial n_{aij}}$ $\displaystyle=$
$\displaystyle\frac{1}{n_{aij}^{o}}+\frac{\Delta t}{n_{aij}^{o}}\sum_{kl\neq
ij}P_{ij,kl},$ (87) $\displaystyle\frac{\partial f_{3+ij}}{\partial n_{akl}}$
$\displaystyle=$ $\displaystyle-\frac{\Delta t}{n_{aij}}P_{kl,ij},$ (88)
$\displaystyle\frac{\partial f_{3+ij}}{\partial n_{r,A}}$ $\displaystyle=$
$\displaystyle-\frac{\Delta t}{n_{aij}}n_{r,B}n_{r,C}K_{r}.$ (89)

##  Appendix C Rate equations for solution of non-equilibrium hydrogen

The radiative rates used in this work are calculated using the prescription of
Sollum (1999), see also Leenaarts & Wedemeyer-Böhm (2006) for a description of
their implementation in a 3D simulation. The angle averaged radiation field
$J_{\nu}$ for each transition is given in terms of a height-dependent
radiation temperature $T_{\mathrm{rad}}$ by setting
$J_{\nu}=B_{\nu}\left(T_{\mathrm{rad}}\right)$, where

$B_{\nu}\left(T_{\mathrm{rad}}\right)=\frac{2h\nu_{0}^{3}}{c^{2}}\frac{1}{e^{h\nu/k_{\mathrm{B}}T_{\mathrm{rad}}}-1}.$
(90)

In the upper atmosphere $T_{\mathrm{rad}}$ is constant, using the values
$(T_{\mathrm{Sol}})$ prescribed for each transition. These values were chosen
to match comprehensive 1D RADYN simulations (Sollum 1999). Below the
photosphere $T_{\mathrm{rad}}$ is equal to the local gas temperature $T$,
these are smoothly joined by setting

$J_{T_{\mathrm{rad}}}\left(z\right)=B_{\nu}\left(T_{\mathrm{sol}}\right)+\left[B_{\nu}\left(T_{\mathrm{crit}}\right)-B_{\nu}\left(T_{\mathrm{sol}}\right)\right]\left(\frac{m_{c}\left(z\right)}{m_{\mathrm{c,crit}}}\right)^{H_{\mathrm{Sol}}}$
(91)

in terms of column mass $m_{c}$, and as well as a cutoff temperature
$T_{\mathrm{crit}}$ and column mass $m_{\mathrm{c,crit}}$. The parameter
$H_{\mathrm{Sol}}$ is defined for each transition, and fit to a Radyn
simulation by Sollum (1999). The critical values of temperature and column
mass are determined by finding the lowest point
$\left(z_{\mathrm{crit}}\right)$ for which
$B_{\nu}\left(\mathrm{T}\right)=2B_{\nu}\left(T_{\mathrm{Sol}}\right)$. Below
this point $T_{\mathrm{rad}}=T$ and above it we use Eqns. 90 & 91 to calculate
$T_{\mathrm{rad}}$. This method decouples the radiation field from the
thermodynamic properties, allowing rapid calculation and fast convergence of
the hydrogen populations.

The radiative rates can then be calculated from $T_{\mathrm{rad}}$ and
$J_{\nu}$, see Sollum (1999) for a detailed derivation. For a lower level $l$
and an upper level $u$ the upwards $R_{l,u}$ and downwards $R_{u,l}$ radiative
rates are

$\displaystyle R_{l,u}$ $\displaystyle=$
$\displaystyle\frac{4\pi^{2}e^{2}}{h\nu_{0}m_{e}c}f_{l,u}J_{\nu},$ (92)
$\displaystyle R_{u,l}$ $\displaystyle=$
$\displaystyle\frac{g_{l}}{g_{u}}e^{\frac{h\nu_{0}}{k_{\mathrm{B}}T_{\mathrm{rad}}}}R_{l,u},$
(93)

where $e$ is the electron charge, $g_{i}$ the statistical weight of level i,
$f_{l,u}$ is the oscillator strength, and $\nu_{0}$ is the line centre
frequency. The radiation temperature $T_{\mathrm{rad}}=T$ when $z<z_{crit}$
and is constant when $z\geq z_{\mathrm{crit}}$. The temperature derivative of
$J_{\nu}$ is

$\displaystyle\frac{\partial J_{\nu}}{\partial T}$ $\displaystyle=$
$\displaystyle\frac{h\nu_{0}}{k_{\mathrm{B}}T^{2}}\frac{e^{\frac{h\nu_{0}}{k_{\mathrm{B}}T}}}{e^{\frac{h\nu_{0}}{k_{\mathrm{B}}T}}-1}J_{\nu},\leavevmode\nobreak\
\mathrm{for}\leavevmode\nobreak\ z<z_{\mathrm{crit}}$ $\displaystyle=$
$\displaystyle 0,\leavevmode\nobreak\ \mathrm{for}\leavevmode\nobreak\ z\geq
z_{\mathrm{crit}},$

and the derivatives of the bound-bound rates are

$\displaystyle\frac{\partial R_{l,u}}{\partial T}$ $\displaystyle=$
$\displaystyle\frac{4\pi^{2}e^{2}}{h\nu_{0}m_{e}c}f_{l,u}\frac{\partial
J_{\nu}}{\partial T},\leavevmode\nobreak\ \mathrm{for}\leavevmode\nobreak\
z<z_{\mathrm{crit}}$ $\displaystyle=$ $\displaystyle 0,\leavevmode\nobreak\
\mathrm{for}\leavevmode\nobreak\ z\geq z_{\mathrm{crit}},$
$\displaystyle\frac{\partial R_{u,l}}{\partial T}$ $\displaystyle=$
$\displaystyle\frac{g_{l}}{g_{u}}\frac{\partial R_{l,u}}{\partial
T},\leavevmode\nobreak\ \mathrm{for}\leavevmode\nobreak\ z<z_{\mathrm{crit}}$
$\displaystyle=$ $\displaystyle 0,\leavevmode\nobreak\
\mathrm{for}\leavevmode\nobreak\ z\geq z_{\mathrm{crit}}.$

The bound-free rates between a lower level $l$ and the continuum $6$

$\displaystyle R_{l,6}$ $\displaystyle=$
$\displaystyle\frac{8\pi}{c^{2}}\alpha_{0}\nu_{0}^{3}\sum_{n=0}^{\infty}E_{1}\left(n\frac{h\nu_{0}}{k_{\mathrm{B}}T}\right),$
(97) $\displaystyle R_{6,l}$ $\displaystyle=$ $\displaystyle
R_{l,6}\left[\frac{n_{l}}{n_{6}}\right]_{LTE},\leavevmode\nobreak\
\mathrm{for}\leavevmode\nobreak\ z<z_{\mathrm{crit}}$ $\displaystyle=$
$\displaystyle\frac{8\pi}{c^{2}}\alpha_{0}\nu_{0}^{3}\left[\frac{n_{l}}{n_{6}}\right]_{LTE}$
$\displaystyle\sum_{n=0}^{\infty}E_{1}\left(\left[n\frac{T}{T_{\mathrm{rad}}}+1\right]\frac{h\nu_{0}}{k_{\mathrm{B}}T}\right),\leavevmode\nobreak\
\mathrm{for}\leavevmode\nobreak\ z\geq z_{\mathrm{crit}}$

where $\left[\frac{n_{l}}{n_{6}}\right]_{LTE}$ is the LTE population ratio,
and $\alpha_{0}$ is the radiative absorption cross-section at the ionisation
edge frequency $\nu_{0}$. The derivatives of the bound-free rates are

$\displaystyle\frac{\partial R_{l,6}}{\partial T}$ $\displaystyle=$
$\displaystyle\frac{8\pi}{c^{2}}\alpha_{0}\nu_{0}^{3}\frac{1}{T\left(e^{\frac{h\nu_{0}}{k_{\mathrm{B}}T}}\right)},\leavevmode\nobreak\
\mathrm{for}\leavevmode\nobreak\ z<z_{\mathrm{crit}}$ $\displaystyle=$
$\displaystyle 0,\leavevmode\nobreak\ \mathrm{for}\leavevmode\nobreak\ z\geq
z_{\mathrm{crit}}$ $\displaystyle\frac{\partial R_{6,l}}{\partial T}$
$\displaystyle=$ $\displaystyle\frac{\partial R_{l,6}}{\partial
T}\left[\frac{n_{l}}{n_{6}}\right]_{LTE}$ $\displaystyle-$ $\displaystyle
R_{6,l}\left(\frac{3}{2T}+\frac{h\nu_{0}}{k_{\mathrm{B}}T^{2}}\right),\leavevmode\nobreak\
\mathrm{for}\leavevmode\nobreak\ z<z_{\mathrm{crit}}$ $\displaystyle=$
$\displaystyle\frac{8\pi}{c^{2}}\alpha_{0}\nu_{0}^{3}\left[\frac{n_{l}}{n_{6}}\right]_{LTE}\sum_{n=0}^{\infty}\frac{\exp{\left(-\frac{nh\nu_{0}}{k_{\mathrm{B}}T_{\mathrm{rad}}}-\frac{h\nu_{0}}{k_{\mathrm{B}}T}\right)}}{nT+T_{\mathrm{rad}}}\frac{T}{T_{\mathrm{rad}}}$
$\displaystyle-$ $\displaystyle
R_{6,l}\left(\frac{3}{2T}+\frac{h\nu_{0}}{k_{\mathrm{B}}T^{2}}\right),\leavevmode\nobreak\
\mathrm{for}\leavevmode\nobreak\ z\geq z_{crit},$ $\displaystyle\frac{\partial
R_{l,6}}{\partial n_{\mathrm{e}}}$ $\displaystyle=$ $\displaystyle 0,$ (101)
$\displaystyle\frac{\partial R_{6,l}}{\partial n_{\mathrm{e}}}$
$\displaystyle=$ $\displaystyle\frac{R_{6,l}}{n_{\mathrm{e}}}.$ (102)

In Fig. 9 we compare detailed SE calculations made with the RH code
(Uitenbroek 2001; Pereira & Uitenbroek 2015) to those using the Sollum
radiative rates. The statistical equilibrium solution calculated using the
Sollum rates closely match the detailed solution for much of the photosphere
up to the mid-chromosphere. In the upper chromosphere, the Sollum treatment of
Lyman lines in detailed balance leads to a higher ionisation fraction and a
lower population in the first excited state. When Lyman lines are treated in
detailed balance in the RH calculation, the result closely matches that from
the Sollum rates.

Figure 9: Statistical equilibrium populations of hydrogen, comparing the
Sollum radiative rates with a detailed calculation using the RH code. Panel a)
shows the temperature and electron density of the FALC atmosphere. Panel b)
compares the ionisation fraction, and panel c) the $n_{2}$ populations. These
are shown for LTE (red-dash), the Sollum radiative rates (blue-dash), and from
two detailed calculations with the RH code, one of which treats Lyman-lines in
detailed balance (green-dash), and a second which calculates the radiative
rates for the Lyman-lines (solid-black), see the legend in Panel c).

##  Appendix D Molecular hydrogen rates

For the time dependent solution of the molecular $\mathrm{H}_{2}$ we use a set
of rates described in table 2, these are often described in terms of the
Arrhenius equation

$\mathrm{arr}\left(\alpha,\beta,\gamma\right)=\alpha\left(\frac{T}{300}\right)^{\beta}e^{-\frac{\gamma}{T}},$
(103)

and its derivative

$\frac{\partial\mathrm{arr}}{\partial
T}\left(\alpha,\beta,\gamma\right)=\left(\frac{\beta}{T}+\frac{\gamma}{T^{2}}\right)\mathrm{arr}\left(\alpha,\beta,\gamma\right).$
(104)

The rate coefficients, $K$, are used in the equations of molecular hydrogen
and the hydrogen ground state. In the case of $K_{6}$ & $K_{7}$ the
$\mathrm{H}^{-}$ and $\mathrm{H}_{2}^{+}$ molecules are assumed to
automatically dissociate.

Table 2: Molecular hydrogen formation and dissociation rates

$K$ | Reaction | Formula | Reference
---|---|---|---
$K_{1}$ | $H+H+H\rightarrow H_{2}+H$ | $5.0e^{-32}T^{-0.25}+2.0e^{-31}T^{-0.5}$ | 11$1$Forrey (2013)
$K_{2}$ | $H_{2}+H+H\rightarrow H_{2}+H_{2}$ | $\frac{1}{8}K_{HHH}$ | 22$2$Palla et al. (1983)
$K_{3}$ | $H_{2}+H\rightarrow H+H+H$ | $\mathrm{arr}(4.67e^{-8},-1.0,5.5e^{4})$ | UMIST 142 33$3$McElroy et al. (2013)
$K_{4}$ | $H_{2}+H_{2}\rightarrow H_{2}+H+H$ | $\mathrm{arr}(1.0e^{-8},0,8.41e^{4})$ | UMIST 135 33$3$McElroy et al. (2013)
$K_{5}$ | $H_{2}+e^{-}\rightarrow H+H+e^{-}$ | $\mathrm{arr}(3.22e^{-9},0.35,102000)$ | UMIST 140 33$3$McElroy et al. (2013)
$K_{6}$ | $H_{2}+e^{-}\rightarrow H+H^{-}\rightarrow H+H+e^{-}$ | $\mathrm{arr}(1.92934e^{-11},-1.27,43000)$ | 44$4$Hirasawa (1969)
$K_{7}$ | $H_{2}+H^{+}\rightarrow H_{2}^{+}+H\rightarrow H+H+H^{+}$ | $\mathrm{arr}(1.5178e-9,-0.4563,2.1812e4)$ | 55$5$Galli & Palla (1998)

##  Appendix E Opacity Binning

To create the group-integrated opacities for the multi-group scattering scheme
we require the absorption ($\kappa$), scattering ($\sigma$) and total
($\chi=\sigma+\kappa$) opacities. The opacity contains a contribution from the
continuum c and from lines l. Detailed Opacity Distribution Functions (ODFs)
are taken from the Merged Parallelised Simplified ATLAS code (MPS-ATLAS)
(Witzke et al. 2021), an upgraded version of ATLAS9 (Kurucz 1970). For the
continuum, the total and scattering opacities are available through the ATLAS
package, tabulated with frequency, temperature and pressure. For the spectral
line contribution, ODFs are used. The ODFs are tabulated on the same
frequency, temperature and pressure grid as the continuum values, with an
additional 12 sub-bins per frequency point. To extract the absorption and
scattering opacity from the total opacity in the ODFs we follow Skartlien
(2000), using the approximation of van Regemorter (1962). The photon
destruction probability of a frequency $\nu$ is given by the probability for
collisional de-excitation from the upper level $j$ to the lower level $i$:

$\epsilon^{l}_{\nu}\approx\frac{C_{j,i}/A_{j,i}}{1+C_{j,i}/A_{j,i}}.$ (105)

Here $A_{j,i}$ is the Einstein coefficient for spontaneous radiative de-
excitation and $C_{j,i}$ is the collisional de-excitation parameter. Using van
Regemorter’s approximation

$\frac{C_{j,i}}{A_{j,i}}=20.6\lambda^{3}n_{\mathrm{e}}T^{-1/2}P\left(\frac{\Delta
E_{\nu}}{kT}\right),$ (106)

where the function $P\left(\frac{\Delta E_{\nu}}{kT}\right)$ is pre-tabulated
by van Regemorter (1962). This allows the scattering and absorption opacities
to be determined using $\epsilon^{l}_{\nu}$

$\kappa^{l}_{\nu}=\epsilon^{l}_{\nu}\chi^{l}_{\nu},\leavevmode\nobreak\
\mathrm{and}$ (107) $\sigma^{l}_{\nu}=(1-\epsilon^{l}_{\nu})\chi^{l}_{\nu}.$
(108)

The total opacities are then calculated

$\kappa_{\nu}=\kappa^{l}_{\nu}+\kappa^{c}_{\nu},\leavevmode\nobreak\
\mathrm{and}$ (109) $\sigma_{\nu}=\sigma^{l}_{\nu}+\sigma^{c}_{\nu}.$ (110)

For each band j of the multigroup scheme, consisting of a set of frequencies
$\Omega_{j}=\\{\nu\\}$, different averages are used to calculate the group-
integrated opacities, see Skartlien (2000) for a detailed discussion. These
are the Rosseland mean opacity $\chi_{j}^{R}$

$\chi_{j}^{R}=\left.\int_{\Omega_{j}}\frac{dB_{\nu}}{dT}d\nu\middle/\int_{\Omega_{j}}\frac{1}{\chi_{\nu}}\frac{dB_{\nu}}{dT}d\nu\right.,$
(111)

the Planck mean opacity $\chi_{\nu}^{P}$

$\chi_{j}^{P}=\left.\int_{\Omega_{j}}\chi_{\nu}B_{\nu}d\nu\middle/\int_{\Omega_{j}}B_{\nu}d\nu\right.,\leavevmode\nobreak\
\mathrm{and}$ (112)

the intensity-weighted mean $\chi_{\nu}^{J}$

$\chi_{j}^{J}=\left.\int_{\Omega_{j}}\chi_{\nu}J_{\nu}^{\mathrm{pp}}d\nu\middle/\int_{\Omega_{j}}J_{\nu}^{\mathrm{pp}}d\nu\right..$
(113)

in terms of a mean intensity $J^{\mathrm{pp}}$, calculated for each bin of the
ODF. This is performed using a 1D plane-parallel reference atmosphere and a
short-characteristics scheme, similar to that used in the MURaM code, with
$\mu=\pm 1/\sqrt{3}$ and arbitrary azimuths. In this work we use an atmosphere
calculated using a column mass average of the full time-series of the Bifrost
public release (Carlsson et al. 2016), shown in Fig. 10. To extrapolate from
the 1D $J^{\mathrm{pp}}$ to the full range of $p,T$ a reference optical depth
is used, we use $\tau_{500}$.

The extrapolation is performed following a procedure similar to that described
by (Collet et al. 2011). For each $x,y$ pixel of each snapshot of the
reference simulation time series we bin $\tau_{500}$ values for temperature
and pressure corresponding to the ODF grid. Additionally we calculate a column
mass averaged atmosphere of the 3D simulation. The log-mean
$\tau_{500}^{\mathrm{3D}}$ is calculated for each temperature and pressure
point in the table. This $\tau_{500}^{\mathrm{3D}}$ is interpolated to fill
any gaps and then extrapolated to the full range of temperature and pressure
used in the ODF table. The column-mass averaged background model is used to
calculate $\tau^{\mathrm{ref}}_{500}$ and $J^{\mathrm{pp}}$. Using the
assumption that $J^{\mathrm{pp}}$ is constant over points with the same
optical depth ($\tau_{500}$) the intensity mean is calculated for the full
table.

The resulting background model, extrapolated $\tau_{500}^{\mathrm{3D}}$,
intensity $J^{\mathrm{pp}}$ and photon destruction probability $\epsilon$ can
be seen in Fig. 10.

Figure 10: The atmospheric properties required to calculate the multi-group
scattering opacities a) The density and temperature of the reference
atmosphere b) Binned optical depth $\tau_{500}^{\mathrm{3D}}$, the white
contour shows the range of $T,p$ values encountered in the 3D simulation, c)
Intensity $J^{\mathrm{pp}}$ used for averaging in the streaming regime and d)
photon destruction probability calculated using the van-Regemorter
approximation.

We then follow the formalism of Ludwig (1992) to join the diffusion (optically
thick) and streaming (optically thin) domains. Firstly the extinction
coefficient is given as

$\chi_{j}=\chi_{j}^{J}e^{\tau_{j}/\tau_{0}}+\chi_{j}^{R}\left(1-e^{\tau_{j}/\tau_{0}}\right),$
(114)

where $\tau_{0}=0.1$ and the group optical depth $\tau_{j}$ is approximated
using the expression

$\tau_{j}\approx\frac{\kappa^{R}_{j}p}{g_{\bigodot}},$ (115)

where $g_{\odot}$ is the gravity at the photosphere. Similarly, the scattering
albedo is

$\left(1-\epsilon\right)_{j}=\frac{\sigma_{j}^{J}e^{\tau_{j}/\tau_{0}}+\sigma_{j}^{P}\left(1-e^{\tau_{j}/\tau_{0}}\right)}{\chi_{j}^{J}e^{\tau_{j}/\tau_{0}}+\chi_{j}^{P}\left(1-e^{\tau_{j}/\tau_{0}}\right)},$
(116)

and the integrated emissivity is

$\left(\epsilon
B\right)_{j}=\frac{\int_{\Omega_{j}}\kappa_{\nu}B_{\nu}d\nu}{\chi_{j}^{J}e^{\tau_{j}/\tau_{0}}+\chi_{j}^{P}\left(1-e^{\tau_{j}/\tau_{0}}\right)}.$
(117)

##  Appendix F Initial hydrogen populations

In Fig 8 we show the departure coefficients from the initial snapshot, which
was computed in Bifrost code (Carlsson et al. 2016). This snapshot is used as
the initial condition for the simulation presented in this paper. The
departure coefficients are seen in Fig. 11. The $\mathrm{H}_{2}$ populations
are very small in the upper chromosphere, and the large departure coefficients
are energetically and dynamically insignificant.

Figure 11: Properties of the non-equilibrium hydrogen populations in the
initial snapshot, through the centre of the enhanced network region. The
panels show a) temperature in the NE simulations, b) the departure coefficient
of temperature comparing $T_{\mathrm{LTE}}/T_{\mathrm{NE}}$, c) the NE
ionisation fraction $n_{\mathrm{H},1}/n_{\mathrm{H,tot}}$, d) the first
excited level of hydrogen $n_{\mathrm{H},0,1}$, and the departure coefficients
of e) the ground state, f) the first excited state, g) protons and h)
molecular hydrogen.

##  Appendix G Comparison to RH

In order to test if the code correctly reproduces the limiting LTE and coronal
equilibrium cases we reproduce the test from Sect. 11.6 of Gudiksen et al.
(2011). In Fig. 12 we compare the hydrogen populations of the MURaM NE module
with those calculated in statistical equilibrium in RH (Uitenbroek 2001;
Pereira & Uitenbroek 2015). The ground state population densities match the
statistical equilibrium calculation for the interior and low chromosphere. In
the upper-chromosphere the ground state population varies strongly from the
statistical equilibrium solution, and in the corona the proton number differs,
although variations remain within a factor of 10. The proton densities match
the statistical equilibrium calculation for the interior, transition region
and corona. Throughout the chromosphere the proton number densities differ
strongly from the statistical equilibrium solution.

Figure 12: Comparison of the ground state (panel a) and protons (panel b) from
the MURaM NE treatment of hydrogen and RH.