# Non-equilibrium helium ionization in an MHD simulation of the solar atmosphere

###### Abstract

The ionization state of the gas in the dynamic solar chromosphere can depart strongly from the instantaneous statistical equilibrium commonly assumed in numerical modeling. We improve on earlier simulations of the solar atmosphere that only included non-equilbrium hydrogen ionization by performing a 2D radiation-magneto-hydrodynamics simulation featuring non-equilibrium ionization of both hydrogen and helium. The simulation includes the effect of hydrogen Lyman- and the EUV radiation from the corona on the ionization and heating of the atmosphere. Details on code implementation are given. We obtain helium ion fractions that are far from their equilibrium values. Comparison with models with LTE ionization shows that non-equilibrium helium ionization leads to higher temperatures in wave fronts and lower temperatures in the gas between shocks. Assuming LTE ionization results in a thermostat-like behaviour with matter accumulating around the temperatures where the LTE ionization fractions change rapidly. Comparison of DEM curves computed from our models shows that non-equilibrium ionization leads to more radiating material in the temperature range 11-18 kK compared to models with LTE helium ionization. We conclude that non-equilbrium helium ionization is important for the dynamics and thermal structure of the upper chromosphere and transition region. It might also help resolve the problem that intensities of chromospheric lines computed from current models are smaller than those observed.

###### Subject headings:

magnetohydrodynamics (MHD) — methods: numerical — radiative transfer — Sun: atmosphere — Sun: chromosphere## 1. Introduction

The solar atmosphere is the Sun’s dynamic outer layer, residing above the convection zone. Here, (partially) ionized gas interacts with magnetic fields, and we see different kinds of waves, jets, and other phenomena. Numerical modeling is a powerful tool for determining which physical processes are important in the solar atmosphere, either together with observations (e.g., Ortiz et al., 2014; De Pontieu et al., 2015), or in their own right (where Martínez-Sykora et al. (2012) and Bourdin et al. (2015) are examples).

Depending on the intended use, various levels of model sophistication are required. A model should include treatment of physics relevant for the regime under consideration.

We aim at modelling the solar atmosphere, all the way from the convection zone to the corona. One difficulty with this is the treatment of the chromosphere and transition region. Here the gas goes from nearly neutral to almost completely ionized, the temperature increases from a few thousand kelvin to a million kelvin, and the dynamics goes from being dominated by the gas pressure to being dominated by the magnetic pressure. A reasonable description of the chromosphere and transition region is important for instance for coronal heating, a long standing problem in solar physics (Klimchuk, 2006; Parnell & De Moortel, 2012), since all the energy that ultimately ends up heating the corona, at some point must have been transported through these regions.

In this paper we present a method for treating the time-dependent ionization state of the atmosphere. Hydrogen is the most abundant element in the Sun, and its non-equilibrium ionization has been found to result in higher temperature wave fronts and lower temperature gas between wave fronts, compared to assuming local thermodynamic equilibrium (LTE, Carlsson & Stein, 2002; Leenaarts et al., 2007). Although the number of helium particles is around ten times less than the number of hydrogen particles, its ionization state does have an effect on the energy balance of the upper chromosphere and transition region (Golding et al., 2014). In fact, for a parcel of neutral solar gas, ionizing half of the helium atoms requires an amount of energy equivalent to raising the temperature from 10 kK to 20 kK.

An account of the non-equilibrium ionization state involves solving a set of rate equations for the atomic population densities. The transition rate coefficients involve frequency integrals over the intensity, so a general description must necessarily take radiative transfer into account. The wave code presented in Rammacher & Ulmschneider (2003) and the code RADYN (Carlsson & Stein, 1992, 1995, 1997) are examples of codes that include this complexity. They solve the hydrodynamic equations and the rate equations together with the radiative transfer. One serious drawback of these codes is that they operate only in one dimension. In 2D or 3D a detailed treatment of radiative transfer, such as that featured in the mentioned codes, is challenging because the amount of computational work needed exceeds the capacity of current supercomputers.

Leenaarts et al. (2007) overcame this difficulty and used simplifying assumptions of the radiation
field that enabled them to carry out 2D radiation-magnetohydrodynamics (RMHD) simulations,
including the effects of non-equilibrium hydrogen ionization. These simulations were performed
using the Oslo Stagger Code (Hansteen, 2004; Hansteen et al., 2007), a predecessor of the stellar atmosphere code Bifrost (Gudiksen et al., 2011).
This method was later implemented in the Bifrost code (Golding, 2010). An
example of a 3D simulation with this package has been made available for
download^{1}^{1}1For download details see IRIS Technical Note 33 at
https://iris.lmsal.com/documents.html
(Carlsson et al., 2016). We expand on this work by including
also a description of non-equilibrium helium ionization. The paper is laid out in the following
way: we explain the developed method in Section 2, and present the
results in Section 3. Finally we summarize and draw conclusions in Section
4.

## 2. Method

The stellar atmosphere code Bifrost solves the equations of radiation-magnetohydrodynamics (RMHD):

(1) | |||||

(2) | |||||

(3) | |||||

(4) | |||||

(5) |

where is the mass density, is the velocity field, is the stress tensor, is the gas pressure, is the current density, is the magnetic field, is the gravitational acceleration, is the internal energy density per unit volume, is the heating due to radiation, is the heating due to heat conduction and viscous and ohmic dissipation, is the magnetic diffusivity, and is the vacuum permeability. We express radiative heating as the sum of different contributions,

(6) |

where is the radiative heating from the photosphere described in Gudiksen et al. (2011), is losses from the chromosphere due to strong lines, is heating from the Lyman- line of hydrogen, and is the heating from EUV photons, corresponding to the thin radiative losses from the transition region and corona (negative ) absorbed in the chromosphere (positive ). Recipes for the three latter contributions are described in Carlsson & Leenaarts (2012) (from now on CL12). These recipes are based on empirical fits. In this work we discuss the effects of non-equilibrium ionization and the absorption of coronal radiation in the chromosphere. These effects enter the RMHD equations through , and . In the remaining part of the method section we describe how we compute these quantities self-consistently with the population densities of hydrogen and helium. For further details on the Bifrost code we refer the reader to Gudiksen et al. (2011).

### 2.1. Equations of State

To close the RMHD equation set, we need to relate to and . is in general given by the relation,

(7) |

where , , and are the Boltzmann constant, the gas temperature, the electron density, and the population density of an atom or molecule in the ionization stage occupying the excitation state . The quantities , and are not present in the RMHD equations, and must be specified through extra equations. The MHD model requires a single temperature for all species, so we express the internal energy density as

(8) |

where is the dissociation, ionization or excitation energy of an element or molecule in the ionization stage occupying the excitation state . In addition, MHD assumes charge neutrality, so that the electron density is given by

(9) |

where denotes a neutral particle, a singly ionized particle, etc. The total number of atomic nuclei of each element (including atoms bound in molecules) is conserved, so we add conservation equations for each element. The additional equations required to constrain and close the set of equations depend on the physical system one wants to model.

We refer to the Equations 7, 8 and 9 and the extra equations used to constrain (see 2.1.1 and 2.1.2) collectively as the equations of state (EOS).

#### 2.1.1 The Local Thermodynamic Equilibrium EOS

The most common way of constraining the population densities is to assume local thermodynamic equilibrium (LTE). In that case the temperature, electron density, and population densities are related through a set of Saha-Boltzmann equations and molecular equilibrium equations (e.g., Mihalas, 1978). This combination of equations leads to an EOS that is local: once the internal energy, mass density, and the elemental abundances are given, , , and , and hence , can be computed directly. Based on these variables other relevant quantities, such as opacities, can be computed. This method is computationally fast because all relevant variables can be pre-computed and read from a table. Most stellar atmosphere codes use this assumption (e.g., Nordlund, 1982; Vögler et al., 2005; Gudiksen et al., 2011; Freytag et al., 2012; Wray et al., 2015). The Bifrost implementation employs tables that depend on and .

#### 2.1.2 Non-Equilibrium Ionization of Abundant Elements

The assumption of LTE breaks down in the layers above the solar photosphere. Radiative transition rates become dominant over collisional rates leading to non-local coupling of different parts of the atmosphere through radiation. The transition rates themselves become so small that the ionization- recombination and molecular association-dissociation timescale can become long compared to typical hydrodynamical timescales in the atmosphere (Joselyn et al., 1979).

In that case the population densities should be determined from a continuity equation:

(10) |

The gains and losses represent processes that add or remove particles out of state , for example collisions with electrons, radiative transitions, or processes that form or destroy molecules. This type of lack of any equilibrium in the population densities is generally referred to as non-equilibrium.

The terms that matter the most in the EOS are those that are associated with the most abundant elements. In the solar atmosphere these are hydrogen and helium (Asplund et al., 2009), which are both susceptible to non-equilbrium effects: Non-equilibrium hydrogen ionization leads to increased temperature fluctuations throughout the chromosphere as shown in detailed 1D wave simulations of the solar atmosphere (Carlsson & Stein, 2002). Leenaarts et al. (2011) investigated non-equilibrium formation of H molecules, which can lead to very low temperatures in the chromosphere because the exothermic H formation rate is too low to counteract the rapid adiabatic cooling between internetwork shocks. Golding et al. (2014) showed that non-equilibrium helium ionization leads to significant differences in temperature in the upper chromosphere and transition region.

Solving the continuity equations involves solving the complete radiative transfer problem because radiative rates enter into the gain and loss terms in Equation 10. Each radiative transition rate coefficient is found by integrating the mean intensity over frequency. The mean intensity for each frequency is essentially found by integrating the source function over the computational domain. This procedure is too costly to apply directly in a multi-dimensional RMHD code where it must be repeated every time step, and typical time steps are of the order of milliseconds. Simplifications that lead to a computationally tractable problem are fortunately possible.

### 2.2. Non-Equilibrium Hydrogen in Bifrost

The Bifrost code already includes an option to compute the EOS including non-equilibrium hydrogen ionization and H molecule formation using such simplifications. A detailed description of the method can be found in Leenaarts et al. (2007, 2011) and Gudiksen et al. (2011). Here we briefly restate the essential points as a foundation for the further extensions to the method that are the topic of this paper.

The non-equilibrium hydrogen ionization method is based on approximations by Sollum, E. (1999) who observed that the mean intensity in non-Lyman hydrogen transitions decouples from the gas temperature in the photosphere, and that the mean intensity above this decoupling point stays constant. This means that each radiative rate coefficient in the chromosphere and above can be described by one parameter - the radiation temperature above the photosphere. In the photosphere and below the radiation temperature is simply the gas temperature.

The continuity equations for atomic hydrogen read

(11) |

where is the sum of the collisional () and radiative () rate coefficients and are population densities in the excitation and ionization stages of hydrogen. The collisional rate coefficients depend on the local electron density and temperature only, and the non-Lyman radiative rate coefficients are computed from Sollum’s radiation temperatures, so that the equations are again local. For all other elements we assume LTE.

The assumption that the radiation field decouples from the temperature in the photosphere does not hold for the Lyman transitions. Instead, these are assumed to be in detailed radiative balance, i.e. , which is equivalent to setting in the hydrogen continuity equations (where subscript 1 denotes the ground state).

A continuity equation for H is also included, with three-body association and collisional dissociation with neutral hydrogen atoms as source and sink terms.

The continuity equations are solved using operator splitting. First the populations are advected using an explicit first order upwind scheme. The source and sink part are solved implicitly together with Equations 8 and 9, for the population densities, electron density and temperature. The pressure then follows from Equation 7.

### 2.3. Heating and cooling in Lyman-

The assumption that the Lyman transitions are in detailed balance breaks down in the upper chromosphere (see for example Vernazza et al., 1981), and the method described in Section 2.2 will thus not be accurate there.

Lyman- photons are, roughly speaking, released from the transition region or from shocks, and absorbed by surrounding cold chromospheric material (CL12). The upward rate coefficient for a bound-bound transition depends on the radiation field. Computing the frequency-dependent radiation field every time we advance the RMHD equations with a timestep is computationally too expensive. Therefore we use a simple one-frequency approach. Lyman- photons that contribute to the transport of energy are those that manage to escape from the location where they are emitted. Deep in the solar atmosphere most Lyman- photons are emitted and absorbed at the same location, i.e. the radiation temperature is equal to the gas temperature. These photons do not affect the energy balance of the gas, so we do not model them. We only consider the photons that are able to escape from where they are emitted.

CL12 use an escape probability to parameterize hydrogen losses. We adopt this escape probability to account only for the photons that make a difference in the energy budget. The escape probability at a specific location, , is modelled as a function of the column density of neutral hydrogen (). This escape probability is monotonically decreasing with increasing . In the corona and transition region there is very little neutral hydrogen in the column above, so . Further down in the atmosphere increases and the escape probability goes to zero. We let the downward radiative rate coefficient in the Lyman- transition be proportional to the escape probability,

(12) |

where is the Einstein coefficient for spontaneous deexcitation. From this rate coefficient we express the frequency-integrated Lyman- emissivity,

(13) |

where is Planck’s constant, is the line center frequency, and is the number density of the upper level of the line. We ignore stimulated emission and express the Lyman- opacity as

(14) |

where is the Einstein coefficient for radiative excitation and is a frequency-averaged profile function. We set Hz which is half of the maximum value of the corresponding frequency-dependent Doppler profile at 10 kK. From the emissivity and opacity we obtain the frequency-integrated mean intensity, , by solving the equation of radiative transfer using a short-characteristics method for decomposed domains (Hayek et al., 2010; Carlsson & Leenaarts, 2012). The upward radiative rate coefficient is then expressed as

(15) |

When these non-zero radiative rate coefficients in the Lyman- transition are included in the hydrogen rate equations, we compute the Lyman- heating,

(16) |

where and denote the population densities of the lower and upper level. All other Lyman lines are still assumed to be in detailed balance. The Lyman continuum is however taken into account, as described in Section 2.4.

Our handling of the Lyman- transfer represents an extreme simplification. Nevertheless, we find it worthwhile. Our description qualitatively produces what we want: cooling in the transition region and in wave fronts, and heating in the colder ambient gas.

### 2.4. Non-Equilibrium Helium Ionization

We obtain non-equilibrium helium population densities by solving the rate equations for helium (Eq. 10). We use a three level model atom, consisting of the ground states of He i and He ii, and the doubly ionized state He iii (see Fig. 1). Golding et al. (2014) showed that the model atom reproduces the correct ionization state remarkably well given its extreme simplification.

The model atom includes collisional ionization/recombination (Arnaud & Rothenflug, 1985, see Section B.2), photoionization/recombination and an extra recombination rate coefficient that mimics the net radiative recombination to excited states. Further details and the derivation of the extra recombination rate coefficient are given in Golding et al. (2014). The photoionization rate coefficient is a frequency integral over the mean intensity weighted by the photoionization cross section. As mentioned earlier, this is a quantity which is too computationally costly to be computed in detail. To make the problem computationally tractable we use a bin formulation, reducing the integral to a sum with only a few terms. It is the EUV photons emitted downward from the transition region and corona that are responsible for ionizing helium in the chromosphere. Photons with wavelengths shorter than 504 Å are prone to ionize neutral helium. However, they might just as well ionize hydrogen in the Lyman continuum. The Lyman continuum has its ionization edge at 911 Å. We include the hydrogen Lyman continuum transition and for that reason we choose the first bin to have its upper limit at 911 Å. We let the last bin have a lower limit at wavelength 20 Å, well below the photoionization edge wavelength for He ii at 228 Å.

Figure 2 shows the continuum opacity per hydrogen atom as a function of wavelength for a parcel of solar gas where hydrogen is 30% ionized, helium is 75% neutral and 25% singly ionized – values typical for the upper chromosphere. Here, and in the transition region, is where the ionization state of helium is most important for the energy balance. Radiation bins are indicated in the figure. They are chosen in the following way: bin boundaries are set at the photoinioniziation edges of He i and He ii, resulting in three main bins which are further split into a number of sub-bins.

In each of the main bins the continuum opacity falls off with decreasing wavelength. We set the boundary between sub-bins at the wavelengths where the relative change of opacity is equal in each of the sub-bins. The upper and lower frequency of the bin is denoted and . We adopt a formulation where x denotes the transition and is the bin number. and are the population densities of the ground and continuum states.

There are three continuum transitions, the ground state He i continuum, the ground state He ii continuum, and hydrogen Lyman continuum. That means that we in principle need 3 photoionization cross sections ( denotes the total number of bins). The photoionization cross section for transition x in bin is denoted . Some of these photoionization cross sections are zero. For instance, corresponding to the transition between ground states He i and He ii is zero in all bins that have their lower wavelength boundary at 504 Å or higher.

We find for a transition x in a specific bin by equating the upward radiative rate coefficients in the binned and continuum formulation,

(17) |

where and are the the frequency dependent mean intensity and the frequency dependent photoionization cross section for the transition x. is the mean intensity that is found using in the opacity. This equation is solved by iteration for a point in the chromosphere of a reference atmosphere. As reference atmosphere we use the initial snapshot of the simulations from Golding et al. (2014). Having all the atomic constants determined, we can compute the photoionization and radiative recombination rate coefficients,

(18) | |||||

(19) |

where is a bin-dependent constant, is the LTE ground state to continuum population density ratio for the transition x, given by the Saha relation (see for instance Mihalas, 1978), and is the mean of the Planck function corrected for simulated emission. We ignore stimulated emission, so that the radiative recombination coefficient is independent of the mean intensity. A detailed derivation of these expressions is given in Section B.1.

Test computations showed that the bin divisions drawn in Figure 2 reproduce the actual photoionization rate coefficients of the reference atmosphere fairly well. These six bins are the ones we use in the simulations presented in this paper.

#### 2.4.1 EUV Radiative Transfer

We compute the heating and cooling in the EUV spectrum based on the actual population densities of H and He. The bin opacity and bin emissivity are necessary to compute the mean intensity, . The opacity contribution for transition x in the th bin is

(20) |

where is the lower level population density of the transition. We include two contributors of photons to the radiation. First, the photons released when hydrogen and helium ions recombine. They give rise to an emissivity that is transition and bin dependent,

(21) |

where is the upper level population density of the transition.

Then, second, we include photons produced by collisional excitation followed by radiative deexcitation in lines in the transition region and corona. An exact account of the resulting emissivity is not possible, as it would require of us to solve the rate equations for all relevant ions. Rather, we assume ionization equilibrium and sum up the line losses from relevant ions in each radiation bin. Photons from lines with wavelengths larger than 911 Å or shorter that 20 Å are treated as photons which, after emission, do not interact with matter, i.e. we have no opacity for these wavelengths. We can express this coronal emissivity in the th bin as

(22) |

where is a temperature-dependent loss function per electron per hydrogen atom per steradian. The exponential factor ensures that we avoid including losses from hot regions that are dense. We set to a value typical of the upper chromosphere to make sure these losses only come from the transition region and corona. We use CHIANTI to compute with the abundances given in sun_coronal.abund and the ionization equilibrium values from chianti.ioneq (Dere et al., 1997, 2009). Figure 3 shows for the bins we use in the simulations presented here. Losses due to helium are dominated by the He ii 304 Å line, and therefore we include helium line losses using the non-equilibrium He ii population:

(23) |

Here the first factor is the photon energy per steradian, the middle factor is the number density of singly ionized helium per hydrogen atom and is the collisional excitation rate coefficient per electron taken from CHIANTI. is added to the appropriate . The total bin emissivity is then finally expressed as the sum of both the contributions,

(24) |

We obtain the bin mean intensity, , with the same formal solver that is used for the Lyman- radiation.

We now have all relevant quantities and can compute the heating rate in the upper chromosphere, transition region and corona caused by the EUV spectrum:

(25) |

### 2.5. Four Simulation Runs

In order to compare the various approximations to the ionization balance of hydrogen and helium we perform a differential study. Therefore we performed four two-dimensional simulation runs starting out from the same initial snapshot. The simulations are meant to be comparable to quiet sun conditions. The spatial domain spans a region 15.8 Mm 16.6 Mm with an equidistant horizontal resolution of 33 km and a vertical resolution of 28 km at Mm, continuously increasing to 150 km in the corona at Mm. Mm is set in the photosphere where the optical depth at 5000 Å is unity. The initial magnetic field has a mean absolute value of 65 G at Mm. The flux is concentrated in four regions separated roughly by 4 Mm, the strongest of which has a negative sign ( Mm) and the three remaining concentrations slightly weaker and with a positive sign. The two concentrations at Mm and Mm form footpoints for looplike structures, and we refer to them as network. In 2D models there is not enough magnetic dissipation to sustain coronal temperatures self consistently. We therefore use a hotplate boundary condition in the corona that will heat or cool the plasma towards a temperature of 1 MK on a timescale of around 400 seconds. The runs differ in the way the EOS, and thus also and , is modeled.

The first run (LTE) treats all elements, including H and He, in LTE. and are computed from the recipes in CL12.

The second run (HION) treats hydrogen in non-equilibrium as described in 2.2, with the Lyman transitions in detailed balance. All other elements are in LTE. and are computed from the recipes given in CL12.

The third run (LYA-HION) treats hydrogen in non-equilibrium as described in Section 2.2, but with the Lyman- transition computed using our one-frequency recipe (see Section 2.3). is computed as described in Equation 16 and is computed from the recipe in CL12.

Finally, the fourth run (HELIUM) treats both hydrogen and helium in non-equilibrium. is computed from Equation 16 and is computed from Equation 25.

All the runs include tabulated losses from helium in the coronal emissivity (see Figure 3), except the HELIUM-run which includes helium losses as given by Equation 23.

Compared to running with an LTE EOS, the computing time needed per time step is 2-3 times longer when using a non-equilibrium hydrogen EOS and 4-5 times longer when using the non-equilibrium hydrogen and helium EOS. The four simulations are run for at least 3000 solar seconds and snapshots are written every 10 seconds. The details of each run are summarized in Table 1.

Simulation | H | He | He304 | ||
---|---|---|---|---|---|

LTE | LTE | recipe | LTE | eq. | recipe |

HION | non-eq. | recipe | LTE | eq. | recipe |

LYA-HION | non-eq. | detailed | LTE | eq. | recipe |

HELIUM | non-eq. | detailed | non-eq. | non-eq. | non-eq. |

## 3. Results

### 3.1. Effect of Non-Equilibrium Hydrogen Ionization

Figure 4 shows temperatures and ion fractions at seconds for all four simulation runs. We first focus our attention on the two left columns, corresponding to the LTE-run and the HION-run. Comparing panels (a) and (b), we see that the HION-run has a hotter chromosphere and more contrasted structure in the temperature than the LTE-run. This is due to long hydrogen ionization recombination timescales, resulting in a slowly changing ionization state. The fraction of H as H ii in the chromosphere of the HION-run (panel (r)) is more stable and less responsive to waves (and other perturbations) than the ion fraction from the LTE-run (panel (q)). The change of internal energy associated with the waves will either go into ionizing atoms (i.e. bound in ions) or increase the temperature. We see this clearer in the vertical space-time cut shown in Figure 5. The location of the cut is indicated by the green vertical lines of Figure 4. We can see waves propagating in the chromosphere in panels (a) and (b) as the tilted brighter lines. The amplitude in temperature of the wave fronts is higher for the HION-run than for the LTE-run. Again, the hydrogen ion fraction of the HION-run (panel (r)) is less structured than the LTE-run ion fraction (panel (q)). This confirms what was reported in Leenaarts et al. (2007).

Helium is in LTE in both the LTE-run and in the HION-run. Panels (i) and (j) in Figures 4 and 5 show the fraction of He as He ii for the two simulations. For the region below the arched structures separating the chromosphere from the corona, this fraction is larger in the HION-run than in the LTE-run. This is simply because of the higher chromospheric temperatures featured in the HION-run. Neutral helium ionizes at 10 kK in LTE. The two ion fractions are different in value, but they correlate well with the patterns and structures seen in temperature.

#### 3.1.1 Lyman- Heating

All four simulations feature Lyman- cooling, but only the LYA-HION-run and the HELIUM-run include the Lyman- heating and cooling self-consistently with the rate equations. This self consistent treatment is the only thing separating the HION-run from the LYA-HION-run, so these two runs are used to identify the effects. We compare the temperatures in Figure 4 (panels (b) and (c)), and the temperature space-time cuts from Figures 5 and 7 (panels (b) and (c) in both figures). We do not see any systematic difference between the two runs. We do observe, however, when comparing the fraction of H as H ii, that including the Lyman- heating leads to a more extended region where the fraction of H as H ii is of order (panels (r) and (s) from Figures 4 and 5). This is expected for two reasons, which we illustrate in Figure 6. First, Lyman- cooling in the transition region will contribute to a stronger net downward rate making it easier for atoms to remain neutral, even at high temperatures. Second, absorption deeper down, in colder regions, will contribute to a stronger net upward rate and more hydrogen occupying the first excited state. This results in more ions since photoionization in the Balmer continuum is the most important ionization process for hydrogen in the chromosphere (Carlsson & Stein, 2002).

### 3.2. Effects of Non-Equilibrium Helium Ionization

We now investigate the effects of non-equilibrium helium ionization and focus first on the two right columns of Figure 4, corresponding to the LYA-HION run and the HELIUM-run. Panels (c) and (d) show the temperature. The difference is not dramatic, but there is a tendency towards structural features in the chromosphere standing more out in the HELIUM-run. These structures are less pronounced in the HELIUM-run ion fractions. See for instance the fraction of He as He ii in panel (l). It shows little or no correlation with the structures that can be seen in the chromosphere temperature in panel (d). Panel (k) displays this fraction in the LYA-HION-run and there is a solid correlation between it and the chromosphere structures showing in panel (c). The vertical and horizontal space-time temperature diagrams of Figures 5 and 7 (panels (c) and (d) in both figures) indicate that the gas in the wave fronts is hotter, and gas between wave fronts is colder in the HELIUM-run than what they are in the LYA-HION-run. We compare the vertical space-time diagrams of the fraction of He as He ii (panels (k) and (l) of Figure 5) and see that these waves are essentially not showing in the HELIUM-run, but they are in the LYA-HION-run.

We explain the effect of non-equilibrium helium ionization in the same way we explain the effect of non-equilibrium hydrogen ionization. The long ionization-recombination timescale of helium prevents the ionization state to respond to waves and other perturbations. The increased internal energy associated with the wave compression will, instead of ionizing the gas, lead to increased temperatures. Conversely, the expanding gas between shocks cools off the material instead of maintaining its temperature by releasing energy from recombining helium ions.

### 3.3. Preferred Temperatures

In the simulations where LTE ionization is used in the EOS, certain temperatures are more frequently occurring than others. We can see this in the horizontal temperature space-time diagrams shown in Figure 7. In the LTE-run (panel (a)) there are many points with a temperature around 6 kK corresponding to blue color. In the two HION-runs there are many points with a temperature of about 10 kK corresponding to orange/red color. Finally, in panel (d) we see that the temperature spans more of the color table range.

The preferred temperatures are also clearly visible as dark bands in Figure 8, where we map out the occurrence rate of points on a grid. Using an LTE EOS clusters grid cells around these temperatures. There is a cut in these panels at (2 kK). This is due to an artificial heating term that kicks in when the temperature drops below this value, effectively setting 2 kK as a lower limit temperature.

If we heat a parcel of gas with LTE hydrogen ionization, the temperature will not rise above 6 kK before all of the hydrogen atoms are ionized. It is similar for LTE helium ionization. As we heat the parcel of gas, the temperature will not exceed 10 kK before all of the He i atoms are ionized. Heating it even more, eventually the temperature rises until we reach 22 kK where all of the He ii ions are ionized before the temperature can increase to higher values. This happens because the Saha equation which governs the LTE ionization is particularly sensitive to temperature. Atoms or ions will ionize over a small temperature range, while the corresponding range of internal energy is large. In the chromosphere these small temperature ranges are centered at the preferred temperatures. There hydrogen and helium act as thermostats when their ionization state is described by LTE.

### 3.4. EOS and Radiative Capability

Using a non-equilibrium ionization EOS might change how a model atmosphere radiates. We use the differential emission measure (DEM) to get a qualitative understanding of this change. Assuming ionization equilibrium, the emergent intensity of a thin line from a column in our simulation can be expressed as,

(26) |

where is the line dependent contribution function determined by atomic data. The DEM, , is defined as . We compute this quantity on a temperature grid ranging from 10-25 kK. To get a statistically reasonable result, we average over all columns of all the snapshots covering the timespan 1000-3000 s.

Figure 9 shows the DEMs corresponding to the four simulation runs. The HELIUM-run DEM falls off smoother than the DEMs from the other runs. At temperatures above 11 kK it deviates from the two HION-run DEMs, whereas they are very similar at lower temperatures. In the temperature range 11-18 kK the HELIUM-run DEM has a higher value by a factor of around 2 compared to the LYA-HION run. The bumps in the DEM curves at 22 kK in all but the HELIUM-run are due to LTE ionization of helium and the resulting 22 kK preferred temperature. Since the LTE-run has less material at temperatures between the two preferred temperatures 6 kK and 10 kK, its DEM in that temperature range is down by an order of of the DEMs from the runs including non-equilbrium hydrogen ionization (not shown here). Although DEMs are used mostly for analysis concerning coronal lines formed at higher temperatures than the ones featured here, the quantity is correlated to the amount of mass and its ability to radiate at a given temperature - also for the low temperatures we are considering here.

### 3.5. Consequence for Modelling Helium Lines

Using a DEM to model helium resonance lines with Equation 26 is not ideal for two reasons. First, the lines are not thin, and second, the assumption of ionization equilibrium is not valid. We elaborate on the latter in Figure 10. Here we show the probability density functions of helium ions as a function of temperature, from the HELIUM-run. The time dependent fraction of He as He ii peaks at a temperature near 10 kK - well below the corresponding ionization equilibrium value at around 50 kK. This happens because the coronal radiation is photoionizing the ”cold” and neutral helium in the upper chromosphere.

#### 3.5.1 Helium Ionization-Recombination Timescales

We use the helium transition rates to compute , the ionization-recombination timescale (Judge, 2005). This timescale is shown in Figure 11 at s for the HELIUM-run (the same snapshot as in Figure 4). Around the two network regions at Mm and Mm the timescales in the chromosphere are the shortest, ranging from s to s. The transition region and corona above these regions are hot and dense, leading to high EUV emissivities () and strong EUV heating of the chromosphere, i.e. large transition rates and short timescales. In the more elevated parts of the chromosphere residing under the arched structures, the timescale is an order of magnitude higher at s. Here the gas has a low density and the incident radiation field is weaker than in the network regions, both of which lead to slow rates and long timescales.

## 4. Conclusion

We present a method for including non-equilibrium ionization of abundant elements in the EOS of RMHD simulations. The method is implemented in the stellar atmosphere code Bifrost. To assess the effects of the dynamics of the chromosphere on the ionization, we present four 2D numerical simulations of the solar atmosphere. They feature different setups for ionization in the EOS: it is either described by a set of Saha equations (LTE) or by a set of rate equations for hydrogen and/or helium resulting in a non-equilibrium ionization state where possible long ionization-recombination timescales are taken into account.

We find that the use of a non-equilibrium ionization EOS affects the state of the chromosphere in the following way: there is a larger variation in the chromospheric temperature which is caused by smaller variation in the ionization state than what is the case for LTE simulations. A slowly changing ionization state means that the internal energy bound in ions is also slowly changing. Any short lived chromospheric perturbation in the internal energy will therefore manifest itself as a perturbation in temperature. Both non-equilibrium hydrogen and helium ionization contribute to this effect. The ionization state of helium impacts the upper chromosphere the most because hydrogen is already highly ionized there.

A side effect of using LTE ionization in the EOS is that grid cells that make up the chromosphere and lower transition region, where the ionization state of the gas plays a decisive role in the energy balance, are likely to have one of three preferred temperatures, 6, 10 or 22 kK, associated with the ionization of H i, He i and He ii respectively. This leads to chromospheric configurations that are not necessarily very physical, and that, in turn, may effect significantly how certain chromospheric lines form. We illustrate this with a DEM calculation that shows an increase by a factor two between 11 kK and 18 kK when including non-equilibrium helium ionization. This change in atmospheric structure might resolve part of the problem that chromospheric spectral lines calculated from RMHD models typically show too little emission (e.g., Leenaarts et al., 2013; Lin & Carlsson, 2015). Changes in the electron density and temperature might affect the first ionization potential effect (FIP) observed in the corona (e.g., Laming, 2015). The non-equilibrium ionization of helium might influence the mechanisms that cause abundance differences of helium in the corona (Laming & Feldman, 2003), and on the anomalous helium line intensities derived from DEM models (Giunta et al., 2015).

The timescale of helium ionization and recombination in the chromosphere span a range to seconds. The lower end of this range is consistent with the timescale we computed in Golding et al. (2014). There we found and gave examples of non-equilibrium effects on the He i 10830 and He ii 304 lines. Since a timescale in the low end of the range is sufficient to cause non-equilibrium effects, these effects might be more pronounced in atmospheric regions where the timescale is in the higher end of the range. We thus expect a strong influence of non-equilibrium helium ionization on the formation of the He i 10830 Å line, which is a popular diagnostics of chromospheric magnetic fields (e.g., Lagg et al., 2004; Centeno et al., 2010; Schad et al., 2013).

Our conlusions are based on a differential study using 2D simulations. Since the physical mechanisms do not change when going to 3D, we expect our conclusions to hold then as well.

While our simulations with non-equilibrium hydrogen and helium ionization are an improvement over previous work, they have some limitations. Our treatment of the radiative transfer in the hydrogen Ly and He ii 304 lines is highly simplified. Also, ion-neutral interaction effects are not included, but they have been found to affect the thermodynamics of the chromosphere (Martínez-Sykora et al., 2012).

At the time of writing a 3D simulation including non-equilibrium ionization of hydrogen and helium is underway. The natural next step would be the inclusion of ion-neutral interactions.

## Appendix A Code Implementation of Non-Equilibrium Helium Ionization

The variables included in experiments with non-equilibrium helium ionization are the temperature, , the electron density, , the atomic hydrogen population densities, , for five bound states ( to ) and the continuum (), the molecular hydrogen population density , and, finally, the helium population densities for the ground states of He i, He ii and He iii (). This adds up to 12 variables, so we need to solve 12 equations. The equations are linearized and solved with the iterative Newton-Raphson scheme. The first 9 equations are the internal energy conservation equation (Eq. 8), charge conservation equation (Eq. 9), a hydrogen atom conservation equation, and six rate equations, one for molecular hydrogen and five for atomic hydrogen. For the three last equations we use two rate equations for atomic helium and one helium atom conservation equation.

The helium and hydrogen atom conservation equations have identical forms. For brevity we let refer to either or and express the atom conservation equations:

(A1) |

where is the total number density of either hydrogen or helium, which is proportional to the mass density with an abundance dependent coefficient, and is the number of atomic states. For hydrogen is 6 and for helium it is 3. The atom conservation equations only depend on the number density and the derivative is

(A2) |

The atomic rate equations of hydrogen and helium also have identical forms. We define

(A3) |

where is the sum of the collisional and radiative rate coefficient for the transition from atomic state to atomic state : . The atomic rate equations then take the form

(A4) |

These equations depend on electron density, temperature and number densities. The derivatives are:

(A5) | |||||

(A6) | |||||

(A7) | |||||

(A8) |

In the old non-equilibrium hydrogen description (the HION-run), all Lyman transitions are assumed to be in detailed radiative balance, i.e. the radiative rate coefficients to or from the ground state of hydrogen () are all set to zero. We have derived approximate expressions for the Lyman- and Lyman continuum radiative rate coefficients, , , and given in Equations 15, 12, 18 and 19, respectively. Of these rate coefficients, only has non-zero derivatives. These are given in the next section. For all other derivatives we refer the reader to the Appendix of Leenaarts et al. (2007). For the rate equation for molecular hydrogen and its derivatives we refer the reader to Appendix B of Gudiksen et al. (2011).

## Appendix B Transition Rate Coefficients

The total transition rate coefficient needed for the rate equations, , is the sum of a collisional part, and a radiative part . In this section we give a derivation of the radiative rate coefficients and provide expressions for the collisional rate coefficients.

### b.1. Derivation of Photoionization Rate Coefficients

The general expression for the photoionization rate coefficient (e.g., Mihalas, 1978) is given by

(B1) |

where is the frequency dependent photoionization cross section and is frequency dependent mean intensity, is Planck’s constant, and corresponds to the difference in energy between the ground and continuum states. The subscript x denotes the transition under consideration. We have adopted a bin formulation. In each bin the photoionization cross section is constant and determined by Equation 17. The mean intensity is also constant for each bin, . and denotes the lower and upper frequency boundaries of the th bin. We can then express the photoionzation rate coefficient as a weighted sum over the binned mean intensity,

(B2) |

where the weights are defined,

(B3) |

The radiative recombination coefficient, ignoring stimulated emission, is expressed

(B4) |

where the Planck function neglecting stimulated emission is

(B5) |

Here is the speed of light and is Boltzmann’s constant. is the LTE ground state to continuum ratio. Given the statistical weights, and , the electron mass, , and the difference in energy between the ground and continuum state, , the LTE ratio is expressed:

(B6) |

In the bin formulation we use the binned photoionization cross section, and a bin integrated . The radiative recombination coefficient then becomes a weighted sum, much like the photoionization coefficient,

(B7) |

where

(B8) |

All derivatives of are zero. depends on the electron density and the temperature. The derivatives are:

(B9) | |||||

(B10) | |||||

where

(B11) |

In the code we interpolate in pre computed tables of and as functions of .

### b.2. Collisional ionization/recombination rate coefficients

The collisional ionization and recombination rate coefficients are

(B12) | |||||

(B13) |

where is a temperature dependent function described in Arnaud & Rothenflug (1985). The derivatives are

(B14) | |||||

(B15) | |||||

(B16) | |||||

(B17) | |||||

In the code we interpolate pre computed tables containing and as functions of .

## References

- Arnaud & Rothenflug (1985) Arnaud, M., & Rothenflug, R. 1985, A&AS, 60, 425, 425
- Asplund et al. (2009) Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481, 481
- Bourdin et al. (2015) Bourdin, P.-A., Bingert, S., & Peter, H. 2015, A&A, 580, A72, A72
- Carlsson et al. (2016) Carlsson, M., Hansteen, V. H., Gudiksen, B. V., Leenaarts, J., & De Pontieu, B. 2016, A&A, 585, A14, A14
- Carlsson & Leenaarts (2012) Carlsson, M., & Leenaarts, J. 2012, A&A, 539, A39, A39
- Carlsson & Stein (1992) Carlsson, M., & Stein, R. F. 1992, ApJ, 397, L59, L59
- Carlsson & Stein (1995) —. 1995, ApJ, 440, L29, L29
- Carlsson & Stein (1997) —. 1997, ApJ, 481, 500, 500
- Carlsson & Stein (2002) —. 2002, ApJ, 572, 626, 626
- Centeno et al. (2010) Centeno, R., Trujillo Bueno, J., & Asensio Ramos, A. 2010, ApJ, 708, 1579, 1579
- De Pontieu et al. (2015) De Pontieu, B., McIntosh, S., Martinez-Sykora, J., Peter, H., & Pereira, T. M. D. 2015, ApJ, 799, L12, L12
- Dere et al. (1997) Dere, K. P., Landi, E., Mason, H. E., Monsignori Fossi, B. C., & Young, P. R. 1997, A&AS, 125, 149, 149
- Dere et al. (2009) Dere, K. P., Landi, E., Young, P. R., et al. 2009, A&A, 498, 915, 915
- Freytag et al. (2012) Freytag, B., Steffen, M., Ludwig, H.-G., et al. 2012, Journal of Computational Physics, 231, 919, 919
- Giunta et al. (2015) Giunta, A. S., Fludra, A., Lanzafame, A. C., et al. 2015, ApJ, 803, 66, 66
- Golding (2010) Golding, T. P. 2010, Master Thesis, University of Oslo
- Golding et al. (2014) Golding, T. P., Carlsson, M., & Leenaarts, J. 2014, ApJ, 784, 30, 30
- Gudiksen et al. (2011) Gudiksen, B. V., Carlsson, M., Hansteen, V. H., et al. 2011, A&A, 531, A154, A154
- Hansteen (2004) Hansteen, V. H. 2004, in IAU Symp. 223: Multi Wavelength Investigations of Solar Activity, ed. A.V.Stepanov, E.E.Benevolenskaya, & E.G.Kosovichev, 385–386
- Hansteen et al. (2007) Hansteen, V. H., Carlsson, M., & Gudiksen, B. 2007, in Astronomical Society of the Pacific Conference Series, Vol. 368, The Physics of Chromospheric Plasmas, ed. P. Heinzel, I. Dorotovič, & R. J. Rutten, 107–114
- Hayek et al. (2010) Hayek, W., Asplund, M., Carlsson, M., et al. 2010, A&A, 517, A49, A49
- Joselyn et al. (1979) Joselyn, J. A., Munro, R. H., & Holzer, T. E. 1979, ApJS, 40, 793, 793
- Judge (2005) Judge, P. G. 2005, J. Quant. Spec. Radiat. Transf., 92, 479, 479
- Klimchuk (2006) Klimchuk, J. A. 2006, Sol. Phys., 234, 41, 41
- Lagg et al. (2004) Lagg, A., Woch, J., Krupp, N., & Solanki, S. K. 2004, A&A, 414, 1109, 1109
- Laming (2015) Laming, J. M. 2015, Living Reviews in Solar Physics, 12, 2, 2
- Laming & Feldman (2003) Laming, J. M., & Feldman, U. 2003, ApJ, 591, 1257, 1257
- Leenaarts et al. (2011) Leenaarts, J., Carlsson, M., Hansteen, V., & Gudiksen, B. V. 2011, A&A, 530, A124, A124
- Leenaarts et al. (2007) Leenaarts, J., Carlsson, M., Hansteen, V., & Rutten, R. J. 2007, A&A, 473, 625, 625
- Leenaarts et al. (2013) Leenaarts, J., Pereira, T. M. D., Carlsson, M., Uitenbroek, H., & De Pontieu, B. 2013, ApJ, 772, 90, 90
- Lin & Carlsson (2015) Lin, H.-H., & Carlsson, M. 2015, ApJ, 813, 34, 34
- Martínez-Sykora et al. (2012) Martínez-Sykora, J., De Pontieu, B., & Hansteen, V. 2012, ApJ, 753, 161, 161
- Mihalas (1978) Mihalas, D. 1978, Stellar atmospheres (2nd ed.; San Francisco, CA: W. H. Freeman and Co.)
- Nordlund (1982) Nordlund, A. 1982, A&A, 107, 1, 1
- Ortiz et al. (2014) Ortiz, A., Bellot Rubio, L. R., Hansteen, V. H., de la Cruz Rodríguez, J., & Rouppe van der Voort, L. 2014, ApJ, 781, 126, 126
- Parnell & De Moortel (2012) Parnell, C. E., & De Moortel, I. 2012, Royal Society of London Philosophical Transactions Series A, 370, 3217, 3217
- Rammacher & Ulmschneider (2003) Rammacher, W., & Ulmschneider, P. 2003, ApJ, 589, 988, 988
- Schad et al. (2013) Schad, T. A., Penn, M. J., & Lin, H. 2013, ApJ, 768, 111, 111
- Sollum, E. (1999) Sollum, E. 1999, Master Thesis, University of Oslo
- Vernazza et al. (1981) Vernazza, J. E., Avrett, E. H., & Loeser, R. 1981, ApJS, 45, 635, 635
- Vögler et al. (2005) Vögler, A., Shelyag, S., Schüssler, M., et al. 2005, A&A, 429, 335, 335
- Wray et al. (2015) Wray, A. A., Bensassi, K., Kitiashvili, I. N., Mansour, N. N., & Kosovichev, A. G. 2015, arXiv:1507.07999