Modeling Mg ii during Solar Flares. II. Nonequilibrium Effects

Solar flares, and other transient energy release events, can dramatically disturb the solar atmosphere, driving it out of equilibrium. The importance of nonequilibrium (NEQ) ionization and recombination in determining the atomic level populations (and therefore the radiative response to the disturbance) will likely vary from species to species and the charge states of those species. However, when modeling the radiative response of the atmosphere to flares it is often the case that the statistical equilibrium (SE) solution is sought, given the computational demands of performing time-dependent simulations that include physical processes, such as partial frequency redistribution (PRD), that are required for certain transitions (recent examples include Rubio da Costa & Kleint 2017; Kerr et al. 2019; Zhu et al. 2019). This is partly mitigated by using the NEQ electron density from dynamic simulations. We address the importance of NEQ ionization on the modeling of Mg ii h and k lines during solar flares, a routinely observed pair of strong spectral lines from the Sun.

In the standard flare model, magnetic reconnection occurs in the corona, releasing vast amounts of magnetic energy, up to 10³² erg (Fletcher et al. 2011). This results in in situ heating and particle acceleration. Electrons (and likely ions) are accelerated at relativistic speeds along flare loops, loosing energy via Coulomb interactions as they travel into denser regions (Brown 1971; Holman et al. 2011; Jeffrey et al. 2019). A broadband enhancement to the solar radiative output results from both thermal and nonthermal processes following plasma heating and ionization. Chromospheric ablations ("evaporations"), with bulk upflows reaching several hundreds of km s⁻¹, and condensations, with bulk downflows of dense material reaching a few tens of km s⁻¹, are also driven during flares, revealed by Doppler shifts of spectral lines (e.g., Fisher et al. 1985; Fisher 1989; Milligan & Dennis 2009; Graham & Cauzzi 2015).

The transition region and chromosphere are both the sites of energy deposition and the origin of the bulk of the radiative output during flares (Fletcher et al. 2011; Milligan et al. 2014). They are crucial locations for both diagnosing the flaring plasma, and as test grounds for models of energy transport during flares. For the former, forward modeling radiation from different flare atmospheres allows us to understand how to extract physical properties of the flaring plasma from observations. For the latter, model–data comparisons allow us to understand how well simulated flare atmospheres driven by a particular energy transport mechanism compare with the actual flaring chromosphere, and therefore how consistent our models of energy transport are with the reality.

Since the launch of the Interface Region Imaging Spectrograph (De Pontieu et al. 2014, IRIS), observations of the Mg ii h and k resonance and subordinate lines have become commonly used in observational studies (e.g., Graham & Cauzzi 2015; Kerr et al. 2015; Liu et al. 2015; Panos et al. 2018; Tei et al. 2018; Tian & Chen 2018), and are attractive as a means to diagnose or critically attack flare (and other) models (e.g., Kerr et al. 2016; Rubio da Costa et al. 2016; Zhu et al. 2019). These lines have been observed in hundreds of flare events since the launch of IRIS, show significant variation from the quiet Sun during flare heating, and form over a range of chromospheric layers meaning they have the potential to diagnose a large part of the lower solar atmosphere. Particularly if they are combined with other IRIS or ground-based observables.

These lines are optically thick, and as such are complex to interpret, requiring forward modeling to understand their formation (Leenaarts et al. 2013). Modeling of the Mg ii NUV spectrum has been a key complementary activity to the analysis of observations, and it is essential that we are modeling these lines in an accurate manner. This is both to provide an accurate physical interpretation of these complex optically thick lines, but also to facilitate high fidelity model–data comparison by which models of flare energy transport are critically interrogated.

In Paper I (Kerr et al. 2019), we used RADYN (Carlsson & Stein 1992, 1997; Allred et al. 2015) radiation hydrodynamic flare atmospheres with the radiation transport code RH (Uitenbroek 2001), a common approach to modeling these lines in flares. This was to investigate how various setups affect the Mg ii solution such as PRD versus complete frequency redistribution (CRD), inclusion of other species, and effects of coronal irradiation. In particular we noted that PRD is still required, resulting in substantial differences to the CRD solution. However, in that work the level populations were obtained under the assumption of SE, since RH is a stationary code. Each atmospheric snapshot was treated in isolation, and the SE solution obtained. This was necessary in order to include the more advanced radiation transfer from the RH code (in particular PRD), that would be computationally very demanding to include in a time-dependent flare simulation. The history of the atmosphere and any time-dependent effects were neglected. In this paper we will investigate the impact of NEQ processes on the formation of Mg ii during solar flares, with the aim of determining if the usual practice of omitting these effects is safe.

The chromosphere is dynamic, especially during solar flares. Carlsson & Stein (1999, 2002) demonstrated, from RADYN simulations of propagating acoustic waves, that the ionization and recombination timescale (τ ∼ 10³–10⁵ s) for hydrogen is long compared to the dynamical timescale, and that if SE is assumed then the ionization fraction is underestimated by several orders of magnitude in certain locations of the atmosphere. The electron density would consequently be very different. Similarly, the requirements for considering NEQ effects on other species have been noted, such as helium (Golding et al. 2014, 2016), O iv (Olluri et al. 2013), and Si iv (Martínez-Sykora et al. 2016). Leenaarts et al. (2013) investigated if NEQ processes were important for Mg ii in the quiet Sun, concluding that whenever the temperature was large enough to produce significant amounts of Mg iii, the relaxation time was shortened, and that using SE was sufficient. Of course flare chromospheres are very dynamic and the conclusions of Leenaarts et al. (2013) might not apply in flaring conditions. We investigate the flaring scenario in this work.

The NEQ atomic level population equation is:

$$\begin{eqnarray}&&\displaystyle \frac{\partial {n}_{i}}{\partial t}+\displaystyle \frac{\partial {n}_{i}v}{\partial z}-\left(\displaystyle \sum _{j\ne i}^{{N}^{{\prime} }}{n}_{j}{P}_{j,i}-{n}_{i}\displaystyle \sum _{i\ne j}^{{N}^{{\prime} }}{P}_{i,j}\right)=0,\end{eqnarray} \tag{ 1 }$$

where v is the atmospheric velocity, N' is the total number of states, t is time, z is the height in the atmosphere, ${P}_{i,j}$ describes the total rates (collisional plus radiative) from i to j, and ${P}_{j,i}$ is the total rate from j to i. The transition rates are functions of the local atmospheric conditions (including energy input) which vary with time in dynamical simulations. If the local thermodynamic state of the atmosphere or the radiation field vary faster than the timescale for ionization and recombination then there is not enough time for the atmosphere to reach equilibrium—the "history" of the atmosphere becomes important.

If the ionization and recombination timescales are sufficiently fast then the populations can be approximated by SE (setting $\partial {n}_{i}/\partial t$ and ∂n_iv/∂z to zero in Equation (1)):

$$\begin{eqnarray}&&\displaystyle \sum _{j\ne i}^{{N}^{{\prime} }}{n}_{j}{P}_{j,i}={n}_{i}\displaystyle \sum _{i\ne j}^{{N}^{{\prime} }}{P}_{i,j}.\end{eqnarray} \tag{ 2 }$$

In situations where the history of the atmosphere is important but SE is used, there may consequently be errors in the population densities of atomic states, and in the synthetic spectra.

Three flare simulations were produced using the RADYN radiation hydrodynamic code (Carlsson & Stein 1992, 1997; Abbett & Hawley 1999; Allred et al. 2005, 2015, J. C. Allred et al. 2019 in preparation), with the flares driven by nonthermal electron beams that were injected into a preflare atmosphere that spanned the subphotosphere through to the corona. Both the code and details of these simulations are described in Paper I (Kerr et al. 2019). The electron beam fluxes were F = [1 × 10⁹, 1 × 10¹⁰, 1 × 10¹¹] erg cm⁻² s⁻¹ (hereafter F9, F10, and F11) which were modeled as power-law distributions with spectra index δ = 5 and a low-energy cutoff E_c = 20 keV. Energy was injected at a constant rate for t = 10 s, and the simulations were allowed to continue to evolve until t = 50 s. Figure 1 shows the evolution of the atmospheres as a function of time in each simulation. The gray shaded portion is the preflare atmosphere.

Zoom In Zoom Out Reset image size

RADYN simulated the time-dependent response of the atmospheres to flare energy injection so that the atomic level populations, electron density, and temperature stratification are all NEQ, but with radiation computed using the simplifying assumption of CRD (necessary to make the dynamic flare problem computationally tractable). Species important for energy balance are considered in the main simulation (H, He, and Ca ii), with other species included in an optically thin radiation loss function. While Mg ii h and k are very strong chromospheric lines, PRD is required to accurately model the radiative losses. Omitting the Mg ii h and k lines is likely safe to do from an energetic balance standpoint, given that the Ca ii H and K lines are included, but in CRD. So, losses from Ca ii H and K are overestimated, but this is mitigated by ignoring Mg ii h and k.

In order to obtain the Mg ii NEQ atomic level populations and emergent profiles we must therefore turn to MS_RADYN, the minority species version of RADYN. This code uses each internal timestep of an existing RADYN RHD solution to solve the NEQ NLTE radiation transport problem for a desired minority species, such as Mg ii. It includes the time-dependent and advection terms when solving the atomic level populations. Since the variables required for this simulation are stored for each internal timestep (and not the output cadence of the main set of variables, which can be many times larger) MS_RADYN can capture changes to the atomic level populations on very short timescales that result from the changing atmospheric state. MS_RADYN was used in this fashion to study C ii emission (Judge et al. 2003) and more recently to investigate radiative transfer effects on Si iv during flares (Kerr et al. 2019).

Our three flare simulations were used as input in MS_RADYN to simulate the atomic level populations and synthetic spectra of a 10-level-with-continuum model of Mg ii, the same as that used in Leenaarts et al. (2013) and that we used in Kerr et al. (2019). This model atom included the h and k transitions, and the subordinate line transitions. To test the impact of NEQ effects we then repeated these simulations, but switched off the time-dependent and advection terms when computing the atomic level populations so that MS_RADYN then used SE (Equation (2)). In this latter series of simulations, the Mg ii problem was solved using the NEQ hydrogen and electron densities (similar to the more typical post-processing of RHD/HD snapshots through RH).

Using MS_RADYN in this fashion meant that the only differences were the terms included in the atomic level population equation (Equation (1) versus (2)), with the same flare atmospheres and background opacities used throughout. Any differences in the level populations and emergent spectra are then due to NEQ processes.

How the emergent profiles differ between each simulation is the most important consideration indicating if there is a negligible, tolerable, or significant difference between NEQ and SE when modeling Mg ii, given our goal of comparing to observations. Here we discuss these differences in each scenario, but we do not address in detail the line features themselves and how they vary with flare strength or how well they compare to observations, which will be the focus of other investigations.

To quantify the difference between the two results we compute the percentage change between the intensities across the resonance and subordinate line profiles, $\tfrac{{I}_{\mathrm{SE}}-{I}_{\mathrm{NEQ}}}{{I}_{\mathrm{NEQ}}}\times 100$, where I_SE is the intensity computed using SE, and I_NEQ is the intensity computed using NEQ. A positive percentage means that I_SE is more intense.

For the bulk of the duration of the flares the differences between the results is only minor or negligible, though there are exceptions. Figures 2–4 show a comparison between the NEQ and SE Mg ii k line profiles, and the percentage change between the two solutions, for the F11, F10, and F9 simulations respectively. The Mg ii 2791 Å line is shown for the F9 simulation in Figure 5. In each figure, panel (a) shows the line profile as a function of time (stacked on the y-axis), and panel (b) shows the percentage change. Panels (c)–(g) show profiles at selected times of interest.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It is immediately clear from these figures that while assuming SE results in very little change during the main heating phase, it does result in under- or overestimations of the line intensity across the whole k line during the initial energy injection and following cessation of the beam. The subordinate line differences are more limited to the line core and near wings. Line shapes and features are preserved, though in the decay phase the subordinate lines return from emission to absorption more rapidly in the NEQ simulations.

Within the first second of the F11 and F10 simulations, there is an intensity change −[5–50]%, with the strongest changes around the emission peaks. These differences rapidly subside in both F11 and F10, until the cessation of energy input after which a very narrow region around the line core shows approximately a 10%–20% change, again with NEQ more intense. In the F10 simulation there are several seconds following the cessation of energy input where the NEQ solution is significantly less intense, up to approximately 230% in the blue emission peak. The F9 simulation evolves more gradually, and shows close agreement between the NEQ and SE solutions in the initial stages. Unlike the other two simulations, however, differences appear during the heating phase, with a larger NEQ result, with a difference on the order of 5%–20%. In the cooling phase the F9 behaves similarly to the F10 simulation, with a stronger SE intensity. The h and k lines are not affected by the same magnitude meaning the k:h line ratio will vary between the NEQ and SE solutions, as discussed in Section 4.

The largest differences appear around the emission peaks in the decay phase, with the wings showing less significant changes. Integrating across the line profiles therefore reduces the discrepancy at these times substantially. The discrepancies in the initial heating phase for the F10 and F11 simulations affected the line wings also, so the percentage change is of a similar magnitude to the specific intensity case at those times. Figures 6(a)–(c) show the lightcurves of the Mg ii k line, integrated ±0.5 Å around the line core.

Zoom In Zoom Out Reset image size

In the moderate-to-strong flares the initial heating ($t\lt 1$ s) shows a local maximum in the NEQ solution, but not in the SE solution. In the F10 flare, the NEQ solution shows a steeper decrease following the cessation of the beam. At other times the lightcurves generally show similar behavior, with small intensity differences (note the logarithmic scales, so that while the lines do appear in close agreement they can actually differ by a few tens of percent).

The k:h line ratio, ${R}_{k:h}$, is a useful metric that can indicate if the lines are optically thick, if there are relative changes between the lines during the flare, or if radiative excitation processes are significant (e.g., Harra et al. 2014). Kerr et al. (2015) noted that in some areas of a flare ribbon observed by IRIS, ${R}_{k:h}$ changed in response to the flare. The ratio decreased slightly, and the relative spread in values decreased (there was a tighter correlation of ${R}_{k:h}$ than in areas outside of the ribbon). They speculated that this could be caused by the h and k lines forming closer together than in the quiet Sun, and/or sampling a more uniform chromosphere during the flare. In the optically thin limit this ratio is ${R}_{k:h}=2$, the ratio of the statistical weights of the k and h upper levels. It is typical that ${R}_{k:h}\approx 1.2$, both in the quiet Sun and in flares (e.g., Kerr et al. 2015; Liu et al. 2015), indicating optically thick line formation. Note, though, that ${R}_{k:h}$ can theoretically have a value of two even in the optically thick case as the ratio of the source functions can in effect take any value (Rathore & Carlsson 2015).

It is important to determine if this ratio is affected by NEQ effects that will confuse the interpretation of ${R}_{k:h}$ variation in flares if forward modeled using SE.

Figure 7 shows ${R}_{k:h}$ in each simulation where again we show the percentage change between NEQ and SE. The magnitude of the variations between the NEQ and SE solutions differs for the h and k lines, meaning that ${R}_{k:h}$ consequently shows differences over time. The magnitude of this difference is relatively small, on the order of <10%. Generally the temporal profile of ${R}_{k:h}$ is preserved, though in the weaker flares the rate of change at the start of the decay phase is smaller in the NEQ solution, and in the F10 and F11 simulations the NEQ ratios decrease somewhat whereas the SE ratios increase within the first second. Both NEQ and SE solutions have R_k:h < 2, indicating optically thick line formation, which we confirmed from inspecting the detailed line formation properties.

Zoom In Zoom Out Reset image size

The ratio is influenced by the relative separation of the h and k lines in the flare atmosphere, and the formation height of the lines is also useful in relating model results to observations of the flaring plasma. Figures 8 and 9 show the formation heights of the h, k, and 2791 Å lines in the NEQ solution (panel (a)), the formation height differences of the lines (Δz) between the NEQ and SE solutions (panel (b)), and the h line source function in the NEQ and SE solutions (panel (c)), for the F11 and F10 simulations respectively. Note that the scales vary between the two flares. The formation height here is defined as the height at which the τ_ν = 1 surface is maximal (that is, we are defining the line core to the part of the line forming highest in the atmosphere, with the greatest opacity).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In both cases during the main phase of the flare there is little difference in the formation heights between the NEQ or SE solutions. However, at the very start of the heating phase in both flares the formation heights can differ. In the F11 simulation this is on the order of kilometers, and in the F10 simulation this is on the order of tens to a few hundred kilometers. The Δz of k and h lines are of different magnitudes, and can be of different directions in the F11 simulation (the k line forms deeper in the NEQ solution whereas the h line forms higher, compared to the SE case). In both cases the subordinate line at 2791 Å actually shows a much greater formation height variation than the h and k lines, forming a few hundred kilometers deeper in the NEQ solution in the F10 flare during the initial decay phase. The h and k lines in that flare do reach Δz ∼ 100 km, but this rapidly decreases to Δz ∼ 5 km or smaller.

Given the generally small differences in the ratios and formation heights, and their short lived nature, we do not envisage using SE will result in any significant misinterpretations when relating line profile features to plasma properties. A possible exception is the subordinate line during the decay phase in weak or moderate flare simulations, or if extremely strong gradients are present near the formation heights of the Mg ii lines.

The differences in the line profile results obtained from the NEQ and SE solutions discussed previously can be understood by studying the ion fraction stratification and the atomic level populations.

Figures 10(a)–(c) show the fraction Mg ii/Mg, χ_Mgii, at the very start of the heating phase of the three simulations. In the two stronger flares, where the flare disturbs the atmosphere in a more impulsive and dramatic manner, χ_Mgii is larger in the NEQ solution than in the SE solution, and consequently there is a larger population in the resonance and subordinate line upper levels. This then results in more radiative decays, and more intense lines. In the F11 simulation there is a narrow region where there is up to an 80% change in χ_Mgii between SE and NEQ. The discrepancy rapidly reduces, both in magnitude and spatial extent. The F10 behaves in a similar manner. The weak flare, however, shows a difference only through the TR, slightly above the formation height of Mg ii, so that the line intensity is not really affected.

Zoom In Zoom Out Reset image size

During the main heating phase the NEQ and SE solutions give largely similar results with only small percentage changes in the F11 and F10 solutions. Figures 10(d)–(f) show t = 6 s, illustrating that in the main heating phase of the flares the NEQ and SE solutions have only marginal differences. In the F9 simulation, which evolved more slowly, differences have started to appear, explaining why there are line intensity differences at this stage in the flare.

Finally, Figure 10(g)–(i) shows the start of the cooling phase, shortly after the beam has ceased depositing energy, where the atmospheric temperature rapidly drops. Here the situation is reversed, with the SE solution predicting more χ_Mgii/Mg. The F11 simulation, in which the flare induced electron density and temperature enhancements were much larger, shows a smaller difference in comparison to the F10 and F9 simulations.

The ion fraction differences between the NEQ and SE solutions are due to the ionization/recombination timescales in each simulation, which vary with atmospheric state (temperature and electron density).

We determined the timescales for the ionization equilibrium of Mg ii/Mg iii following the methodology of Carlsson & Stein (2002) and Leenaarts et al. (2013): For a selected timestep of a flare simulation, the temperature was increased by 1% throughout the atmosphere and the rate equations for the full magnesium model atom were solved as a function of time (keeping the hydrodynamic state constant at the perturbed value), following the relaxation of the population densities from the initial state toward the new equilibrium. The relaxation timescale was calculated from a fit of the time evolution of the population density of Mg iii to the analytic solution for a two-level atom:

$$\begin{eqnarray}&&n(t)=n(\infty )+(n(0)-n(\infty )){e}^{-t/{\tau }_{\mathrm{relax}}},\end{eqnarray} \tag{ 3 }$$

where n(t) is the population density of Mg iii at time t, $n(\infty )$ is the SE population density in the perturbed atmosphere, n(0) is the equilibrium population density in the initial atmosphere, and τ_relax is the relaxation timescale.

Figures 11, 12, and 13 show the relaxation timescales (top panels) and temperatures (bottom panels) at various times in the three simulations. The relaxation timescales are initially on the order of τ_relax ∼ 1–2 s in the upper chromosphere. Through the first second of flare heating in the F11 and F10 simulations, the relaxation time decreases to the order of τ_relax < 0.1 s in response to the substantial temperature and electron density increase in the chromosphere. NEQ effects are prominent within this initial second of heating, because even though the relaxation timescale is decreasing, the atmosphere evolves very impulsively. In the main heating phase the atmosphere evolves somewhat more slowly, so that Mg ii ionization equilibrium keeps pace given the small τ_relax.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the F9 simulation τ_relax = 1 does not decrease as much as in the stronger flares, but the atmosphere also evolves more gradually so that NEQ effects only become apparent when the dynamics become faster. At these times 0.1 < τ_relax < 1 s, and Mg ii is somewhat out of equilibrium.

During the decay phase of the F9 and F10 simulations τ_relax rises back to the order of 1 s when the temperature drops, so that there is initially a large disagreement between NEQ and SE while the atmosphere cools rapidly. After the initial sharp decrease in temperature (and electron density), the rate of change of the atmosphere where Mg II forms is smaller so that the differences between NEQ and SE get smaller (but are still present to some degree), and Mg ii is only partially out of equilibrium.

In the F11 case the temperature is so much larger at the end of the heating phase that even through the decay phase the relaxation timescale is still sufficiently small that the discrepancies are reduced.

Nonequilibrium effects on the formation of Mg ii spectra during solar flares of different magnitudes have been investigated using RADYN and MS_RADYN simulations. The time-dependent NLTE NEQ atomic level populations and synthetic spectra were computed, and compared to the NLTE SE solution.

While line profile shapes are preserved in each solution, the intensities of lines can differ, sometimes substantially, if NEQ effects are taken into consideration. These largely appear in the initial heating phase (t < 1 s), and in the initial seconds of the decay phase. Investigation of Mg ii ionization equilibrium showed that changes in the atmospheric state can lower the relaxation timescale, meaning that Mg ii can be very close to equilibrium during the main heating phase, but that when the temperature rapidly falls Mg ii is again driven out of equilibrium.

In Kerr et al. (2019) we found that PRD is required to accurately forward model Mg ii in flares, with intensity changes of 200%–1000% between the CRD and PRD solutions, and certain features only appearing due to redistribution. With MS_RADYN we are able to include NEQ effects, but the line profiles are computed assuming CRD. PRD is seemingly more important to include than NEQ effects. While the magnitudes of differences between NEQ and SE were generally smaller than the differences between CRD and PRD, and shorter lived (PRD was required throughout the duration of the flare), if one is interested in the initial flare heating, or in the decay phase, then NEQ effects should ideally be considered.

A time-dependent code capable of including both NEQ and PRD effects in flares should be developed, but this is a computationally demanding endeavor. For the moment, if we wish to study the initial heating and cooling, we can capture both NEQ and PRD by using a multistep process. Radiative hydrodynamic flare atmospheres are produced by RADYN. The RADYN solutions are used in MS_RADYN to obtain the NEQ CRD Mg ii populations. Those NEQ populations are fed into RH which, using the modifications we described in Kerr et al. (2019) to fix the level populations, will solve the PRD Mg ii radiation transfer.

Another useful flare line observed by IRIS is Fe xxi 1354.1 Å, which forms in 10 MK plasma and is likely to experience NEQ effects (Bradshaw & Klimchuk 2011). It is commonplace to model this line using data from the CHIANTI atomic database (Dere et al. 1997, 2019), under the assumption of ionization equilibrium and optically thin emission (e.g., Young et al. 2015; Polito et al. 2019). While the latter assumption is safe, ignoring NEQ effects may not be correct. A similar comparison to that presented here for Mg ii should be performed, comparing the predicted level Fe xxi ion fraction from MS_RADYN to those from CHIANTI given model atmospheres from RADYN.

G.S.K. was funded by an appointment to the NASA Postdoctoral Program at Goddard Space Flight Center, administered by USRA through a contract with NASA. This research was supported by the Research Council of Norway through its Centres of Excellence scheme, project number 262622, and through grants of computing time from the Programme for Supercomputing. J.C.A. acknowledges funding support from the Heliophysics Supporting Research and Heliophysics Innovation Fund programs.