Modelling Mg II During Solar Flares. I. Partial Frequency Redistribution, Opacity, and Coronal Irradiation

TL;DR

This paper investigates the physics necessary for accurate Mg II line modeling during solar flares, emphasizing the importance of partial frequency redistribution, non-LTE effects, and coronal irradiation in radiation transfer simulations.

Contribution

It systematically compares modeling approaches for Mg II lines in flares, highlighting the physics required for realistic forward modeling and modifications to the RH code for including irradiation effects.

Findings

PRD remains essential during flare simulations.

Full angle-dependent PRD influences solutions but is computationally intensive.

Including Mg I NLTE has negligible effect on Mg II lines but affects the NUV continuum.

Abstract

The Interface Region Imaging Spectrograph (IRIS) has routinely observed the flaring Mg II NUV spectrum, offering excellent diagnostic potential and a window into the location of energy deposition. A number of studies have forward modelled both the general properties of these lines and specific flare observations. Generally these have forward modelled radiation via post-processing of snapshots from hydrodynamic flare simulations through radiation transfer codes. There has, however, not been a study of how the physics included in these radiation transport codes affects the solution. A baseline setup for forward modelling MgII in flares is presented and contrasted with approaches that add or remove complexity. It is shown for Mg II: (1) PRD is still required during flare simulations despite the increased densities, (2) using full angle-dependent PRD affects the solution but takes significantly longer to process a snapshot, (3) including Mg I in NLTE results in negligible differences to the Mg II lines but does affect the NUV quasi-continuum, (4) only hydrogen and Mg II need to be included in NLTE, (5) ideally the non-equilibrium hydrogen populations, with non-thermal collisional rates, should be used rather than the statistical equilibrium populations, (6) an atom consisting of only the ground state, h & k upper levels, and continuum level is insufficient to model the resonance lines, and (7) irradiation from a hot, dense flaring transition region can affect the formation of Mg II. We discuss modifications to the RH code allowing straightforward inclusion of transition region and coronal irradiation in flares.

Equations

$$I_{\lambda \mu }=\int \frac{1}{\mu }{S}_{\lambda }(z)e^{-{\tau }_{\lambda }\nu }{\chi }_{\lambda }\mathrm{d}z=\int C_{\lambda ,\mu }(z)\mathrm{d}z,$$

Full text

Sun: chromosphere - Sun: flares - Sun: UV radiation - radiative transfer - methods: numerical

1 Introduction

Solar flares are transient yet dramatic events in the solar atmosphere that release tremendous amounts of energy ($>10^{32}$ ergs per event) and drive space weather. It is thought that energy released in the corona during magnetic reconnection is transported to the transition region and chromosphere via directed beams of non-thermal particles (typically electrons), which lose energy via Coulomb interactions. This results in intense plasma heating and ionisation, mass flows, and consequently the broadband enhancement of the solar radiative output (Brown, 1971; Holman et al., 2011; Fletcher et al., 2011, and references therein). Recently it has also been suggested that alternative forms of energy transport can potentially play a role in transporting energy to the lower atmosphere, such as high frequency Alfvénic waves (Fletcher & Hudson, 2008; Reep et al., 2018; Kerr et al., 2016).

The particle acceleration mechanism is still uncertain, but the properties of the accelerated electrons can be inferred from non-thermal X-ray bremsstrahlung. Observations from the Reuven Ramaty High Energy Solar Spectroscopic Imager (RHESSI; Lin et al., 2002) have been extensively studied, and can be used to drive models of solar flares. Simulations of solar flares using non-thermal electron beams to transport energy can then be compared to observations to determine how consistent this model is to reality.

Spectroscopic observations of solar flares are essential if we wish to extract information about the state of the flaring plasma, which can be compared to advanced flare simulations. Both the spectral lines themselves, and the diagnostics they provide, are observables with which we can critically attack models of flare energy transport. Modelling the spectral lines under different physical conditions can also help us build diagnostics. The combination of observations and numerical modelling is powerful as a means to shed light on the physical processes at play in flares.

Spectral lines that originate from the solar chromosphere are of particular interest, as the chromosphere is the location of energy deposition, and the origin of the bulk of the radiative output in flares (Fletcher et al., 2011; Milligan et al., 2014). Since it’s launch in 2013 the Interface Region Imaging Spectrograph (IRIS; De Pontieu et al., 2014) has provided high spectral- and spatial- resolution observations of the solar chromosphere and transition region. The strongest chromospheric lines observed by IRIS are the Mg ii h & k resonance lines and the Mg ii subordinate triplet which are emitted in the near-ultraviolet (NUV). Indeed, these lines are amongst the strongest lines in the solar spectrum. These lines offer the potential to study a large swathe of the chromosphere, but are optically thick and require forward modelling with radiation transport codes to fully appreciate the complex formation properties, and, ultimately, to extract the information they carry.

The h & k lines are transitions $3p\,^{2}\!P^{\rm o}_{1/2,3/2}-3s\,^{2}\!S_{1/2}$ ($\lambda=2802.70$ Å & $\lambda=2795.52$ Å), and the subordinate triplet transitions are $3d\,^{2}\!D_{3/2}-3p\,^{2}\!P^{\rm o}_{1/2,3/2}$, and $3d\,^{2}\!D_{5/2}-3p\,^{2}\!P^{\rm o}_{3/2}$ ($\lambda=2790.78$ Å, $\lambda=2797.92$ Å & $\lambda=2798.00$ Å; the latter two are blended).

Early observations using balloon-borne and rocket experiments revealed that these lines are strongly affected by opacity, with central reversal features almost ubiquitous among the quiescent Sun profiles, though with single peaked profiles in sunspots and shallower reversals in plage and active regions, and intensity ratios that indicate optically thick formation (Lemaire & Skumanich, 1973; Kohl & Parkinson, 1976). Variations of the line width, intensity, emission peak separation and other features were compared to source types (Staath & Lemaire, 1995), but only one flare observation existed prior to the IRIS era (Lemaire et al., 1984). The subordinate lines were relatively little studied but were known to be in absorption on disk and in emission off-limb (Feldman & Doschek, 1977).

Modelling efforts used Mg ii lines to help build semi-empirical atmospheres and to compare models of temperature stratification to observations (e.g. Lemaire & Skumanich, 1973). It was demonstrated that partial frequency distribution (PRD) is required to model the formation of the resonance lines (e.g. Milkey & Mihalas, 1974), rather than the computationally more tractable problem of complete frequency distribution (CRD). In outline, PRD can be described as follows. Atomic species that are strongly scattering and which form in a relatively low density environment may not experience a sufficient number of elastic collisions to destroy coherency between an incident photon and the scattered or re-emitted photon. Photons absorbed in the line wings are re-emitted in the wing in PRD, where it is easier to escape, whereas in CRD they can be redistributed to the core. In PRD the absorption profile is not equal to the emission profile, and the source function is dependent on frequency. See Uitenbroek (2001, 2002) and Hubený (1982) for clear discussions of PRD and the form of redistribution functions

In the IRIS era, the Mg ii resonance and subordinate lines have been studied in many solar features, including in hundreds of flares or similar transient heating events. In flares, it has been noted that the profiles mostly appear single peaked in flare ribbons, are significantly enhanced and broadened over the quiet Sun, and show various Doppler shifts of a few $\times 10$ km s*-1*, and asymmetries, that have been attributed to mass flows (e.g. Kerr et al., 2015; Liu et al., 2015; Graham & Cauzzi, 2015; Tei et al., 2018; Tian & Chen, 2018). Despite appearing single peaked the lines remain optically thick, with a line intensity ratio of $R\sim 1.2$ ($R=2$ in the optically thin scenario) (Kerr et al., 2015). The subordinate lines undergo similar changes, and are almost exclusively in emission during flares. Panos et al. (2018) found that certain profiles (extremely wide with a blueshifted central reversal) appear at the leading front of flare ribbons, and related the line ratio to opacity to infer a higher opacity at flare peak.

Detailed and accurate modelling is required to understand the behaviour of these lines in flares, to confirm the origin of observational findings, and to develop diagnostics. Recent modelling has largely focussed on the quiet Sun (e.g. Leenaarts et al., 2013a, b; Pereira et al., 2015; Carlsson et al., 2015), which has related atmospheric properties to line features. Some flare modelling has been performed. Kerr et al. (2016) noted differences in the Mg ii line profiles that originated from atmospheres heated by different energy transport mechanisms, due to the resulting mass motions, similar to observations but much too narrow. In one of those cases the mass motions shifted the absorption profile so that the line appeared single peaked and asymmetric. Liu et al. (2015) found that a smaller electron beam spectral index (i.e. deeper heating) resulted in more intense line wings, and that the line core was affected by coronal pressure, with larger pressure acting to fill in the central reversal. A parametric study was performed by Rubio da Costa & Kleint (2017), who noted the effect of varying temperature, electron density, bulk velocity, and microturbulence on the emergent profiles, though not in a self-consistent manner. They noted that increasing the electron density would result in a single peaked profile, which was also noted by Zhu et al. (2019), who modelled a strong flare that produced a very high electron density. Zhu et al. (2019) also used more up to date Stark broadening data to explain the narrower-than-observed synthetic spectra, though were unable to fully reach the observed line widths.

When modelling Mg ii and other species that either require advanced radiation transport (e.g. PRD) or are not included in radiation hydrodynamic simulations, which generally only consider a select number of transitions important for energy balance, post-processing is required. The procedure has been to combine snapshots of dynamic atmospheres (from either a hydrodynamic or radiation hydrodynamic simulation) with static radiation transport codes. This means that statistical equilibrium is imposed on the solution, though non-equilibrium effects may be important (Carlsson & Stein, 2002). Some studies using semi-empirical flare atmospheres have also been performed. Clearly a dynamic flare simulation capable of non-equilibrium radiation transport that considers blends and PRD is the ideal, but such a resource is not currently available, in large part due to computational demands.

There can be numerous features of these radiation transport codes that can be variously used, such as employing CRD or PRD, or the number of additional atomic species to solve alongside the one of interest. Given the importance of the Mg ii as a window on the flaring chromosphere we investigate how these features may or may not affect the solution in flares to assess if the typical approaches used to model Mg ii spectra in flares is appropriate and sufficient.

Specifically, we look at if PRD is required, if the treatment of PRD affects the results, if a small model atom can be used, if Mg i must be included in the solution, and the impact of re-solving the hydrogen populations in statistical equilibrium. We also investigate transition region and coronal irradiation. All Mg ii modelling discussed here employs statistical equilibrium to obtain the atomic level populations. Paper 2 in this series will discuss non-equilibrium effects (Kerr et al., 2019a).

2 Numerical Resources and Flare Simulations

2.1 RADYN

RADYN (Carlsson & Stein, 1992, 1997) is a radiation hydrodynamics numerical code that models the solar atmosphere, with a particular focus on the chromosphere, by considering the coupled non-linear equations of hydrodynamics, radiation transport, and non-equilibrium atomic level populations (for certain species). It has the facility to simulate solar flares by the injection of energy via a directed beam of non-thermal particles (Abbett & Hawley, 1999; Allred et al., 2005, 2015) or by downward propagating Alfvénic waves (Kerr et al., 2016). We focus on electron beams for the remainder of this study. For a recent description of the RADYN code we point the reader to Kerr et al. (2019b), and to Allred et al. (2019), and references therein. We note important details below.

The non local thermodynamic equilibrium (NLTE) and non-equilibrium (NEQ) radiation transfer problem is solved in detail for species important for chromospheric energy (H, He i, He ii, Ca ii). The bound-bound and bound-free transitions from these species act as sources and sinks of energy, with non-local effects included. Additional radiative cooling is included by summing the emissivities from all transitions in the CHIANTI atomic database (Dere et al., 1997; Landi et al., 2013), excluding the bound-bound transitions solved in detail. Backwarming and photoionisation/photoexcitation by X-ray, extreme-ultraviolet, and ultraviolet emission from the flare-heated corona is included by injecting a donward-directed incident radiation from the loop apex. Thermal conduction is a modified form of Spitzer conductivity, saturating at the free-streaming limit (Smith & Auer, 1980). Radiation transport assumes complete frequency redistribution, but we truncate the Lyman lines at 10 Doppler widths in an effort to avoid overestimating radiative losses, which can happen if CRD is employed when PRD should really be used (e.g. Uitenbroek, 2002).

Flares are simulated by injecting a non-thermal particle distribution at the loop apex, with a power-law energy spectrum defined by the instantaneous energy flux $F$ [erg cm*-2* s*-1*], above a low-energy cutoff $E_{\mathrm{c}}$ [keV], and with a spectral index $\delta$. The energy losses and propagation of the beam of particles are modelled using the Fokker-Planck equations, which include the important scattering terms. The Fokker-Planck solver has recently been updated from the implementation described in Allred et al. (2015), and will be described more fully in Allred et al. (2019). Non-thermal collisional ionisation and excitation from the ground state of hydrogen is included using the methods of Fang et al. (1993).

2.2 Flare Simulations

Three flares of different strengths were simulated using RADYN, driven by non-thermal electron beams that were injected into a pre-flare atmosphere shown as the shaded portion in Figure 1. The $z=10$ Mm pre-flare atmosphere spanned the sub-photosphere through to the corona, with an apex temperature of $T_{\mathrm{cor}}=3.2$ MK and apex hydrogen density $n_{\mathrm{H,cor}}=6.6\times 10^{9}$ cm*-3*. The photosphere ($z=0$ where $\tau_{5000}=1$) had a temperature $T=5800$ K. These atmospheres were initially in radiative equilibrium, with non-radiative heating applied to maintain the corona and photosphere (column masses $c_{\mathrm{mass}}<1\times 10^{-6}$ g cm*-2* and $c_{\mathrm{mass}}>7.6$ g cm*-2*). A reflecting boundary condition was employed at the loop apex to mimic incoming disturbances from the other leg of the flare loop, and a transmitting boundary condition was used at the base of the loop. The electron beam fluxes were $F=[1\times 10^{9},1\times 10^{10},1\times 10^{11}]$ erg cm*-2* s*-1* (hereafter F9, F10, & F11) which were modelled as power law distributions with spectral index $\delta=5$ and a low-energy cutoff $E_{c}=20$ keV. Energy was injected at a constant rate for $t=10$ s, and the simulations 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, where it is clear that the atmosphere is disturbed more strongly and rapidly in the highest flux simulation. Electron densities and macroscopic flows were notably larger in the F11 simulation.

2.3 RH

RH (Uitenbroek, 2001) is a NLTE radiation transfer code that solves the coupled equations of statistical equilibrium (SE) and radiation transfer, using the MALI (Multi-level Approximate Lambda Iteration) formalism of Rybicki & Hummer (1991, 1992). Overlapping wavelengths (blends), and the effects of partial frequency redistribution are included where appropriate. A model atmosphere is provided to RH, with the temperature, electron density, macroscopic velocity, and microturbulent velocity defined on a given depth scale (in our case a 1D plane-parallel geometry, though other geometries are available). Atomic level populations are solved assuming statistical equilibrium for species specified by the user, with input atomic models containing details of bound-bound and bound-free transitions, collisional rates, photoionisation cross-sections, recombination rates and other relevant data. Multiple species are solved simultaneously. As well as opacities and emissivities from species selected to be treated in NLTE, background sources are included in LTE. These include those atomic species specified by the user, H*-* bound-bound & bound-free processes, Thomson and Rayleigh scattering, H free-free processes, and bound-free processes in OH and CH molecules.

When forward modelling radiation from dynamic atmospheres it is necessary to solve each time as an independent atmosphere. Unfortunately this means that any NEQ effects are not accounted for. In essence we neglect the ‘history’ of the simulation that led to that instantaneous snapshot, which can lead to an erroneous ionisation fraction (e.g. Carlsson & Stein, 2002, discuss this in the case of hydrogen). This is somewhat mitigated by using the non-equilibrium electron density provided in the model atmosphere fixed in the RH solution. Non-equilibrium effects are necessarily sacrificed in order to utilise the more advanced radiation transfer than is typically included in dynamic simulations.

To include hydrogen level populations with non-equilibrium effects and non-thermal ionisations/excitations when computing opacity sources in the flaring atmospheres, the populations from RADYN can be provided using the ‘OLD_POPULATIONS’ starting condition and setting hydrogen to ‘PASSIVE’ (‘PASSIVE’ species are treated only as a background opacity source, whereas the NLTE radiation transfer is solved for ‘ACTIVE’ species). This method allows the user provided populations to be used to compute the background opacity, but the NLTE lines and continua will not be computed by RH. We have instead modified RH to use user provided populations (in our case obtained from RADYN) for an ‘ACTIVE’ species, whilst skipping the level population recalculation step in the iterative loop so that populations remained fixed. RH will therefore use the provided populations directly to solve the NLTE radiation transfer, and for processes such as charge exchange that affect other species. This was implemented by introducing a new starting condition for the initial guess of atomic level populations before the iterations begin. During the iterative loop this ‘guess’ is unchanged. Other species not using this keyword are solved in the usual way, where the populations are recalculated in the iterative loop until convergence is achieved. Our modifications allow us to use NEQ populations from RADYN with the more advanced radiation transfer offered by RH, such as PRD and overlapping transitions. Note that we are forced for the time being to use the SE populations for Mg ii, as this species was not included in the RADYN flare simulations. A follow up work will investigate those effects, and for now we restrict this study to the other physics involved in the modelling of the Mg ii spectra. While the focus of this work is Mg ii and not hydrogen, we introduce our modifications here because the Mg ii results presented here were computed using this version of the code. Future research will use the hydrogen lines computed alongside the Mg ii profiles, and will employ this mechanism for other species.

Snapshots of our RADYN flares atmospheres were post-processed through RH, with a cadence of $0.5$ s in the heating and initial cooling phase ($t=0-15$ s), and with a cadence of $1$ s for the remainder of the simulation ($t=15-50$ s), giving 65 snapshots in total per simulation.

Our ‘standard’ setup, which serves as the benchmark, is as follows. As well as a 10-level-plus-continuum singly ionised Mg atom (we use the model atom from Leenaarts et al. 2013a), the NLTE atomic level populations were solved for: a 23-level-plus-continuum Silicon atom (that included transitions of Si i & Si ii, with the continuum level being the ground state of Si iii); an eight-level-plus-continuum carbon atom (that included transitions of C i & C ii, with the continuum level being the ground state of C iii); a thirteen-level-with-continuum oxygen atom (that included transitions of O i with the continuum level being O ii). Additionally we provided the NLTE, non-equilibrium atomic level populations of a five-level-plus-continuum hydrogen atom, and a five-level-plus-continuum singly ionised calcium atom, obtained from the RADYN simulations. The hydrogen populations that we provided also included non-thermal collisional excitations and ionisations from the ground state of hydrogen according to Fang et al. (1993).

For the lines in which coherent scattering was modelled we used the fast approximation of angle-dependent PRD treatment developed by Leenaarts et al. (2012), also known as ‘hybrid’ PRD (H-PRD). This scheme approximates full angle-dependent PRD (AD-PRD) for atmospheres with macroscopic flows, with substantially lower computational expense. It works by transforming from the observer’s rest frame to the rest frame of the moving gas parcel, allowing the assumption of radiation field isotropy to be used. In our simulations we employ H-PRD for the h & k lines and subordinate triplet. Since two of the subordinate triplet lines share an upper level ($3d\,^{2}\!D_{3/2}$) we also modelled cross redistribution effects (XRD), also known as Raman scattering. A comparison to AD-PRD is given in Section 4.

2.4 Note on H in CRD During Flares

While RADYN allows the NEQ H populations to be obtained (with non-thermal effects also included), it uses CRD when modelling the hydrogen lines, potentially leading to inaccuracies in the hydrogen ionisation fraction and electron density. The hydrogen ionisation fraction stratification, and therefore electron density, is affected by the choice of CRD or PRD (e.g. Hubeny & Lites, 1995). This will affect the ionisation fraction of other species. In addition to the impact of the electron density, other processes such as charge exchange with neutral hydrogen or protons, and optical pumping, can be important for certain species. For example, charge exchange is important for setting the O i/O ii fraction (Carlsson & Judge, 1993).

This is somewhat mitigated by truncating the Lyman lines at 10 Doppler widths, limiting the radiation losses from the wings. Furthermore, in flares densities are significantly enhanced, reducing the impact of PRD on the Lyman lines.

We forward modelled hydrogen using RH, with RADYN flare atmospheres as input, where Ly$\alpha$ and Ly$\beta$ were treated using both CRD and PRD in SE (with the caveat that the electron density was fixed). It can be concluded from our simulations that the CRD assumption employed by RADYN is, in fact, not so incorrect for moderate-larger flares (the F10 and F11 simulations). Using the NEQ RADYN hydrogen populations computed in CRD will not lead to overly spurious results when forward modelling other species. In fact the inclusion of NEQ and non-thermal effects seems to have a larger impact than the choice of CRD or PRD. The SE populations in PRD and CRD, the NEQ populations in CRD, and the resulting Ly $\alpha$ line profiles are are shown in Appendix A. Further discussion of the treatment of the hydrogen populations is found in Section 7.

3 Partial or Complete Frequency Redistribution

3.1 Emergent Intensities and Profiles

It is well known that PRD effects are important for the formation of Mg ii in the quiet Sun, but in flares the density is significantly enhanced, increasing collisions and potentially driving the formation closer to CRD (c.f. Section 2.4). If it transpired that CRD could be used for Mg ii in larger flares, in which the densities are higher, then the computational problem becomes more tractable.

The three flare simulations were processed through RH using the standard setup, and again with the Mg ii lines computed in CRD. Figures 2,3,4 & 5 show the profiles of the Mg ii k line in the F11, F10, & F9 simulations, and the Mg ii 2791 Å line in the F11 simulation, respectively. In each of those figures the PRD and CRD solutions are compared. Panels (a,b) show the profiles as functions of time, where the colour represents radiation temperature. Panels (c-f) show individual snapshots, and the percentage change between the two methods of solution, $(I_{\mathrm{CRD}}-I_{\mathrm{PRD}})/I_{\mathrm{PRD}}\times 100$.

It is immediately clear that unlike the Lyman lines, PRD is still required even in the strongest flare simulation, despite the increased densities. As flare strength (and consequently the electron densities) increases, PRD effects become less important, with the F9 simulation showing up to $\sim 1000$ % intensity changes compared to $\sim 200$ % in the F11 simulation. The resonance line cores and emission peaks were far less affected, if at all, while the inner wings showed the largest differences. While some features in the line wings were preserved in the CRD solution, for example the shoulder in the red wing in the F11 simulation at $t>7$ s (Figure 2(e,f)), not all were. In the F11 simulation a feature appears in the blue wing that is only present in the PRD solution, similarly in the F9 simulation a red wing feature forms only in PRD. We will comment on these asymmetries in Section 3.2. Given that both line intensities and, perhaps more importantly, profile shapes are not fully preserved in CRD we recommend that PRD is always used to model the Mg ii resonance lines even in large flares.

The subordinate lines are typically modelled in CRD, but our results suggest that PRD (with XRD) is required for these lines also. Figure 5 shows that while differences over much of the line are modest relative to the h & k lines, there are locations were a few $\times 10$ % differences are present, and that, similar to the resonance lines, certain features only appear in the PRD solution.

3.2 Line Formation

To understand the formation of line features we can use the contribution function to the emergent intensity, $C_{\lambda\mu}(z)$, which effectively describes where in the atmosphere the radiation originates (e.g. Magain, 1986; Carlsson, 1998). Integrating through a chosen depth scale yields the emergent intensity $I_{\lambda\mu}$. Defining on a height scale $z$, we can write the specific intensity in a plane-parallel semi-infinite atmosphere as:

[TABLE]

where $\mu=\cos\theta$ describes the viewing angle ($\theta$ being the angle between the line of sight and the normal), $S_{\lambda}(z)$ is the source function (with the line source function being independent of $\lambda$ in the CRD case), $\tau_{\lambda}$ is the optical depth, and $\chi_{\lambda}$ is the monochromatic opacity. Hereafter the $\mu$ subscript is dropped and we focus on near-disk centre results.

Figure 6 illustrates the differences between the PRD and CRD line formation for the Mg ii k line in the F11 simulation at $t=7$ s. In that figure the background image in the top panels (a,c) is the source function. The frequency dependence of the source function in the PRD solution is clear, with a larger source function in the core compared to the wings, since it is dependent on the redistribution function and the radiation field. In the CRD case the line source function is uniform in frequency, with redistribution from the wings to core. The source function follows the local temperature at all wavelengths, decoupling from the Planck function at the same height in the chromosphere. In the PRD case, however, different frequencies decouple from the Planck function at different heights due to coherency. Here we define the line core as the wavelength at which the height of $\tau_{\lambda}=1$ is maximum (i.e. the part of the optically thick component that forms highest in the atmosphere), and at this time the core is located in the redshifted component ($\sim\lambda_{\mathrm{rest,k}}+25$ Å. Even though in PRD the source function has decoupled from the Planck function much lower in the atmosphere compared to the CRD case which remains coupled until $\sim 0.9$ Mm (the orange dashed line shows the Planck function), the source functions at the emission height are identical, and have increased in response to the flare.

Perhaps the largest difference is that there is an increase in the source function in the blue wing, extending over a few hundred km, present in the PRD solution but absent in the CRD solution. This grows in strength with time, along with the bump in the redwing that is shared in both the PRD and CRD solutions.

The cause of these two bumps in the wings is different, as shown in the bottom panels Figure 6(b,d), which show the contribution function. Emission is considered optically thick if it forms close to the $\tau_{\lambda}=1$ layer, and optically thin if it forms far above that layer. In both CRD and PRD there are two optically thick line components, one nearly stationary, and the other redshifted by approximately $15-30$ km s*-1*. This redshifted component is caused by a dense condensation that is propagating deeper with time, emitting Doppler shifted Mg ii photons. At this time ($t=7$ s) it has reached a depth sufficient for the opacity in the red wing of Mg ii to increase. Previously the red wing feature was weaker and optically thin, still caused by emission from the condensation but from a greater altitude so that the Mg ii red wing opacity was not so increased. However, the blue wing feature in the PRD solution is not caused by an upflowing source. Instead it is a PRD effect. The absorption profile has shifted to the red so that there is a preference for absorbing redder wavelengths. In PRD, bluer photons are re-emitted coherently without scattering to other wavelengths, where they would be trapped by the much larger opacity. This feature grows in strength because the opacity in the core and the red wing increases when the condensation propagates deeper. Note also that this component is a mix of optically thick and thin emission.

Figure 7 shows the contribution function, opacity, and source functions for selected wavelengths in the k line at $t=7$ s in the F11 simulation, showing in a more quantitative sense the differences between the PRD and CRD solutions. At all wavelengths the PRD and CRD source functions converge between $z=1-1.1$ Mm. The larger temperatures and densities drive more collisions reducing the coherency fraction there, so that the CRD and PRD solutions are similar.

The Mg ii 2791 Å subordinate line formation is shown for comparison in Figure 8. Similar features to the resonance lines are present, with the same difference between the PRD and CRD solution. The redshifted component caused by the condensation is optically thin at this time because the subordinate line forms lower in the atmosphere, so that the condensation has to propagate deeper to sufficiently increase the opacity structure for this line. Despite this, the increased source function and contribution function in the blue wing can be seen. Note here that the subordinate lines form in the upper chromosphere only a couple of hundred km deeper than the resonance lines. This is quite different from the quiet Sun scenario, where they form in the lower atmosphere (Pereira et al., 2015). This means that quiet Sun diagnostics cannot be applied in the flaring scenario (see also Figure 11 in Kerr et al., 2019b).

A similar effect was found in the Ca ii K line, comparing the CRD and PRD solutions for that line revealed a blue wing bump that grew in strength with an optically thick red wing component. As with the Mg ii subordinate line this happened somewhat later in the Ca ii K line since it forms lower in the chromosphere, so the condensation took slightly longer ($\sim 1.5$ s) to reach heights where the effect became noticeable.

This highlights the ambiguity resulting from interpreting asymmetries and Doppler motions in optically thick lines, with the blue wing feature being a result of a condensation not evaporation (see also Kuridze et al., 2015; Kerr et al., 2016; Brown et al., 2018).

It is important to note here that these results used a value of microturbulence of $2$ km s*-1*, the nominal value used in RADYN. Defining the magnitude of microturbulence present is ambiguous, particularly in flares. Increasing microturbulence will increase the width of the Doppler core, having a consequential effect on the amount of redistribution. Some authors have attempted to use increased microturbulence to explain the narrower-than-observed line widths, though this alone is unlikely to be the main cause (Rubio da Costa & Kleint, 2017; Zhu et al., 2019). However, as an experiment we increased microturbulence to a value of $10$ km s*-1*. This smeared the stationary and redshifted optically thick components, and smoothed out the blue-wing bump so that the profile shapes appeared more similar in CRD and PRD. The CRD solution still overestimated the wing intensity. Thus, our conclusion that PRD is necessary stands, regardless of the level of microturbulence chosen.

4 Treatment of Partial Frequency Redistribution

Previously we adopted the fast angle-approximated PRD treatment (H-PRD) of Leenaarts et al. (2012), which saves considerably on computational time and was shown to be a sufficiently accurate approximation to full angle-dependent PRD (AD-PRD). The H-PRD scheme is typically employed when simulating Mg ii in flares using RH (e.g. Kerr et al., 2016; Rubio da Costa & Kleint, 2017; Reep et al., 2019; Zhu et al., 2019).

This H-PRD scheme was tested by Leenaarts et al. (2012) in the scenario of relatively small velocity fields compared to those present in flares. The mass motions driven during solar flares can be significant, with downflows of several tens of km s*-1*, and upflows of a few hundred km s*-1*. Therefore, this approximated scheme may not be applicable during flares, a prospect that we test in this section.

Using the F11 simulation we ran several snapshots with the Mg ii h & k and subordinate lines computed using AD-PRD. The snapshots chosen were $t=[2.5,3.5,4.5,5.5,6.5,7.5,8,9.5]$ s, which covered several points during the development of mass flows in the flare, sampling a range of velocity gradients. To save somewhat on computational time we used the Mg ii atomic level populations from the previously computed H-PRD simulation as the starting solution so that the initial guess was closer to the final result. Even with this effort to speed up convergence, the difference in the time taken to complete these simulations was substantial in comparison to the hybrid scheme. Using AD-PRD the time to complete each snapshot was 3-4 days, compared to only a few minutes for the H-PRD case. The AD-PRD solution took approximately $400$ times longer to complete than the H-PRD solution. Lowering the main RH convergence threshold to $1\times 10^{-2}$ in the AD-PRD experiments reduced convergence time to around one day, with negligible difference in the AD-PRD emergent Mg ii profiles (at most a $\sim 0.5$ % change). Performing high cadence studies flare simulations using AD-PRD will obviously be a computationally expensive and time-consuming endeavour, particularly if multiple flare simulations are involved.

Figure 9 shows the Mg ii k line computed for the F11 simulation using H-PRD, and AD-PRD. There is excellent agreement in line shape, with all features reproduced. There are, however, intensity differences. This can vary from only a few percent to up to $15-20$ % in narrow wavelength regions across the line where redistribution effects are maximised (recall the blue wing feature discussed previously). Balancing the computational cost to the accuracy of the results, we believe that, for most purposes, H-PRD is acceptable for use in flare simulations, despite the localised intensity differences that may result. It is, of course, up to an individual to decide if these intensity differences are tolerable for their specific study.

5 Using a Simpler Model Atom

It was demonstrated by Leenaarts et al. (2013a) that to adequately model Mg ii in quiescent atmospheres, an atom that included excited upper states was required, in order to support an ionisation/recombination loop that populated the h & k upper levels via recombinations to higher lying excited states, followed by cascades through to the h & k upper levels. To show the impact of using too simple a model atom in the flaring scenario, where atmospheric extremes could potentially lead to even more spurious results if too simple an atomic model is used, we show here the results from a three-level-with-continuum Mg ii atom, consisting of the Mg ii ground state, the h & k line upper levels, and the ground state of Mg iii.

As shown in Figure 10(a) there are non-negligible differences in both line intensity and in the shape of the profiles. This is a direct result of smaller population densities of the resonance lines upper levels, due to the lack of recombinations to higher lying excited states and subsequent cascades to the $3p\,^{2}\!P^{\rm o}_{3/2,1/2}$ levels. The lack of these recombinations changes the Mg ii ion fraction and opacity stratification, resulting in a somewhat lower $\tau_{\lambda}=1$ height, compared to the larger model atom. The lower formation height means that the source function is smaller, due to both the lower temperature at that height and the smaller upper level populations. Figure 10(b) shows the formation heights of the lines, along with the temperatures at those heights, and Figure 10(c) shows the lower population density when using the 3+1 level model atom.

Impacts of velocity gradients consequently appear at different times, with different magnitudes, so that the line profiles do not always share common features during the flare simulations. This is demonstrated in Figure 11, which shows the source and contribution functions for the smaller model atom (compare to the same in Figure 6(a,b) for the standard 10+1 level atom). For the larger model atom, the fraction of Mg ii at the height of the condensation is sufficient to produce a substantially larger opacity compared to the smaller model atom, resulting in a redshifted optically thick component to the line. At this same time for the smaller model atom, the condensation results in a weaker, optically thin, feature to appear in the red wing.

We recommend that the model atom of Leenaarts et al. (2013a) be used, or at least a model atom capable of supporting cascades through higher lying excited levels, as described by Leenaarts et al. (2013a).

6 Including Mg I

It has been shown that for modelling quiet Sun h & k profiles, the inclusion of Mg i mainly affected the very far wings with no effect on the line cores (Leenaarts et al., 2013a), and consequently if one is interested primarily in the h, k, or subordinate lines then omitting Mg i is safe to do. We assess here if this result also applies in the more extreme environment of flares. A 75-level-with-continuum model of Mg was used, with 65 levels of Mg i added to the 10 level model Mg ii atom plus the ground state of Mg iii.

Even with the inclusion of Mg i in our model atom, the majority of chromospheric/photospheric Mg exists as Mg ii (Figure 12(a)). Through the temperature minimum region $X_{Mg\textsc{i}}\sim 10$ %, dropping through the chromosphere to $X_{Mg\textsc{i}}\sim 1$ % or smaller. Figure 12(a) shows the Mg ion fractions at $t=5$ s into the F11 flare simulation. In that figure we also show the percentage change of $X_{Mg\textsc{ii}}$ resulting from including Mg i. Despite this small fraction, the amount of Mg ii is still somewhat reduced, which changes the opacity stratification, and the intensity of the Mg ii lines. This has only a small impact on the core of the line, resulting in a tolerable intensity difference in the line cores (on the order of a few percent). Further into the line wings, however, the effect becomes noticeable with a $15-20$ % difference as shown in Figure 12(b). The reduced fraction of Mg ii in the lower atmosphere reduces the opacity, so that the wings form somewhat deeper (Figure 12(c))

Including Mg i in our experiments increased the time to reach convergence by a factor of $\sim 2.5$. Given the minor intensity differences we feel it is safe to ignore Mg i if one is only interested in the line cores.

IRIS also observes the (quasi-) continuum near $\lambda=[2810-2834]$ Å, the intensity of which is affected by the opacity of the h & k absorption profiles (Kowalski et al., 2017). If one is concerned with modelling the h & k line far wings or the continuum in that region Mg i should probably be included.

7 Including Other Species

Hydrogen is an important source of opacity in the solar atmosphere, dominating the background and continuum at NUV and optical wavelengths. In our benchmark, and in the other experiments, we used the hydrogen atomic level populations obtained directly from the RADYN simulations, which account for NEQ ionisation and non-thermal effects. These populations may not always be available to a user wishing to simulate Mg ii so we tested the impact on Mg ii of allowing RH to solve the SE equations for hydrogen.

The F11 RADYN simulation was re-processed with hydrogen (and calcium, though this would not really affect the Mg ii result) solved in SE. The hydrogen populations can differ substantially between the solutions, shown in Figure 13, which has a resulting effect on the NUV/optical opacity. With a larger amount of hydrogen in the n = 2 state, the atmosphere is more opaque to Balmer wavelength photons. The LTE populations are also shown for comparison.

The Mg ii lines were not greatly affected by the hydrogen solution. Differences in the core intensity were negligible, though with increasing $\Delta\lambda$ from the line core into the line wings, discrepancies began to increase. Figure 13 shows the k line computed with hydrogen in SE and with the NEQ, non-thermal hydrogen populations, at $t=4$ s in the F11 simulation (the 2791Å subordinate line shows similar differences). As can be seen in that figure, differences range from almost zero to a few percent in the core, increasing to $>10$ % into the wings. This is to be expected as the Mg ii line opacity dominates the core, but in the wings the underlying hydrogen Balmer continuum becomes more important.

Figure 13(c) shows various sources of opacity in the Mg ii k line. The solid black line shows the full opacity (Mg ii, H NEQ plus the other NLTE and LTE species included as described previously) at $\lambda_{\mathrm{rest,k}}$, and the dashed black line shows $\lambda_{\mathrm{rest,k}}-2$ Å. The background H NEQ opacity (green line) is many orders of magnitude smaller than the core Mg ii opacity, though does contribute more to the wing opacity. The gold dashed line shows the opacity of NEQ H plus the other NLTE and LTE species excluding the Mg opacity. There is an increase in the temperature minimum region opacity compared to the NEQ H-only case, but the core of the Mg ii line is still not affected by the inclusion of these other species. Mg ii and hydrogen are the main contributors.

Given the relative un-importance of other species as sources of opacity at these wavelengths, we conclude that it is not necessary to include species other than hydrogen when forward modelling the Mg ii NUV lines observed by IRIS. Our tests showed no impact on the emergent profiles if the other species were included or not. Ideally one would use the hydrogen NEQ populations from RADYN, or a similar code, but these effects are confined to the wings and so the hydrogen SE populations are probably satisfactory for most purposes. However, the computational expense of including the other species was not prohibitive and it is useful to have lines from multiple species with which to compare against observations.

Another reason one may wish to include the more accurate NEQ hydrogen populations is their role in processes such as charge exchange. We did not study the importance of charge exchange with neutral hydrogen or protons on Mg ii, but it is known to be important for O i (Lin & Carlsson, 2015) and for Si iv (Kerr et al., 2019b).

8 Photoionization from coronal and transition region radiation

The corona in flares can reach temperatures in excess of $10$ MK, with ablation of chromospheric material increasing the density sufficiently so that there is a substantial increase in the flux of X-ray, extreme-ultraviolet, and ultraviolet emission (together: XEUV). At the same time the TR density increases, strongly enhancing radiation produced there. This optically thin emission is directed both outward to the observer, and downward towards the lower solar atmosphere, where it can heat the plasma and photoionise. For example, the emission of the He i $10830$ Å line is strongly influenced by coronal radiation, with photoionisation of helium and subsequent recombination to the $2s~{}^{3}S_{1}$ level (e.g Golding et al., 2014). The plasma heating effects of XEUV backwarming is included in the RADYN simulation (as are photoionisations) so that the atmospheric structure being solved by RH also includes this source of additional heating. However, the potential impacts of an increased photoionisation rate is not considered, as downward directed coronal radiation is not typically included in post-processing of flare simulations with radiation transport codes. An exception to this was an investigation into the origin of H2 emission, including using a semi-empirical flare atmosphere. The fluorescence from strong transition region lines was included as a downward incident radiation that resulted in enhanced molecular transitions (Jaeggli et al., 2018).

There is already functionality in RH to introduce a coronal source of radiation, but this requires a bespoke input file for each snapshot, with specific intensity of the injected radiation defined as a function of wavelength, covering all wavelengths that RH models. Producing a bespoke input file for hundreds of snapshots is a laborious task, especially if the atomic models are changed (which changes the wavelengths on which the injected intensity must be defined). It is preferable to provide a single common input file, and have RH internally compute the intensity to inject based on the atmospheric stratification. Using a similar method to that employed by RADYN we modified RH to permit a straightforward inclusion of coronal irradiance when post-processing flare simulations.

Contribution functions, $G({n_{e},T,\lambda})$ from all transitions in the CHIANTI atomic database were tabulated on a grid that spanned $\lambda=[1,40500]$ Å with $\Delta\lambda=1$ Å bins, $\mathrm{log}T=[4,9]$ with $\Delta\mathrm{log}T=0.05$ bins, and $\mathrm{log}N_{e}=[6,16]$ with $\Delta\mathrm{log}N_{e}=0.5$ bins. Bound-bound transitions already modelled by RH were omitted, as was Helium since those transitions are optically thick. The emissivity at each temperature and electron density in the grid was computed by $j_{\lambda}=A_{b}G(\lambda,ne,T)n_{e}n_{H}$, where $A_{b}$ is the elemental abundance, and $n_{H}$ is the hydrogen density. The abundance table sun$\_$coronal$\_$2012$\_$schmelz$\_$ext.abund from the CHIANTI database was used where low FIP elements are enhanced , while the high FIP elements are not (Schmelz et al., 2012; Caffau et al., 2011; Lodders et al., 2009). For $n_{H}$ we use the typical assumption for a hot optically thin plasma that $n_{H}=0.83n_{e}$. When computing $j_{\lambda}$ we modelled the lines as Gaussians with thermal broadening according to the temperature in each cell of our grid, and a bin size of $1$ Å. The units of this input grid were then erg cm*-3* s*-1* sr*-1* Å*-1*.

This grid of emissivities was input to RH, and the resulting emissivity in each grid cell of the atmospheric model was obtained by performing a trilinear interpolation to the local values of temperature and electron density, for each wavelength bin required by the RH active set of atoms and molecules. The wavelength bins are irregular in RH with some exceeding the bin size of the input grid. When that occurred, we created a new temporary grid that had finer resolution than the input grid. The input emissivities were interpolated at each point on this new wavelength grid, which was then integrated through wavelength to obtain the total emission within the original RH bin (erg cm*-3* s*-1* sr*-1*). This was divided by the RH bin size to return the averaged emissivity in that bin (erg cm*-3* s*-1* sr*-1* Å*-1*). In this manner the emission was conserved through the interpolation, and the potential of a transition falling through a missing wavelength bin was removed. Finally, the emissivity was converted to SI units (W m*-3* sr*-1* Hz*-1*) and integrated through height, providing the specific intensity spectrum (W m*-2* sr*-1* Hz*-1*) that was injected at the apex of the loop.

Several snapshots of the injected spectrum from the F11 simulations are shown in Figure 14 showing how the irradiation varies with time and with atmospheric state. Panels (a,c,e) show the case where emission from grid cells with $T>50$ kK, and (b,d,f) from cells with $T>100$ kK, illustrating that the lower-mid TR dominates the irradiating spectrum. This assumes optically thin emission, though it is likely that opacity effects will be present for certain species. We have already removed Helium for this reason, and future work will investigate in detail the opacity of certain species that are likely to show opacity effects, such as C iii, refining the emissivity grid accordingly. For this present work we performed two additional experiments. In the first we removed transitions of Si i - iv based on the finding in Kerr et al. (2019b) that Si iv lines can become optically thick in flares. In the second we also removed transitions of C i - iii. Figure 14(g,h) show the injected spectra for these two scenarios (both used a $50$ kK temperature threshold.

Emission from grid cells with $T>50$ kK only was initially considered for the irradiation. Coronal/TR irradiaton did affect the Mg ii formation at certain times during the simulation, the significance of which depended on the transitions included in the irradiating spectrum. Increased photoionisations at these times lowered the Mg ii fraction, decreasing the upper level populations of the resonance and subordinate lines, and decreasing the line intensity. These effects were only present at $t<5$ s, when the coronal irradiation was maximum. After $t=5$ s the injected spectrum became much weaker, due to the compression of the mid-lower transition region, the layer which contributed most to the irradiation. Both the resonance and subordinate line intensities were affected, though the profile shapes were preserved. Subordinate line cores were depressed by $\sim 20$ %, and the inner wings by $\sim 40$ %. The resonance line inner cores were little affected but moving outward from the inner cores the lines were depressed by up to $\sim 20$ % into the inner wings. Figure 15(a,b) shows the Mg ii k line and $2791$ Å line profiles at a time when irradiation had non-negligible effects. Removing Si i-iv transitions and then additionally C i-iii reduced the magnitude of the intensity decrease, suggesting that these transitions contributed a large amount of Mg ii photoionisations.

If the temperature threshold at which TR/coronal emissivity begins to be considered is increased to $T>100$ kK then the irradiating intensity drops significantly and there is no impact on the Mg ii formation. This tells us that it is the enhanced radiation originating in the lower-mid TR that has the potential of photoionising Mg ii.

These experiment were repeated for the F10 simulation, which showed no impact on the Mg ii formation even with the lower $50$ kK threshold. In that simulation there is not sufficient heating at higher densities to enhance radiation from the mid-lower TR to the point that it irradiates the chromosphere to the same degree as found for the F11 simulation.

An in-depth investigation of the appropriate temperatures and species to consider for irradiation will be the focus of future work, as this process is also of interest in the formation of the He i $10830$ Å and He i $\mathrm{D}_{3}$ lines, which have recently shown interesting dimming effects at the leading edge of flare ribbons (e.g. Xu et al., 2016; Libbrecht et al., 2019). For now, we note that lower-mid TR and coronal irradiation does have the potential to impact the formation of Mg ii, should ideally be included, but that further research is required to determine the most accurate means to do so. Those modelling Mg ii in flares should be aware of this additional source of ionisation.

9 Summary & Conclusions

We have compared various approaches of modelling the Mg ii NUV spectra during solar flares, using three radiation hydrodynamic flare simulations from the RADYN code, that were post-processed through the RH radiation transport code. This is an effort to ensure that we are including sufficient physics to properly model these lines, which have been routinely observed in flares by the IRIS spacecraft, and offer intriguing diagnostic potential of the flare chromosphere and a crucial observable with which to critically attack flare models. These model approaches either added or removed complexity, with the resulting impact on the Mg ii lines discussed both conceptually and quantitatively.

In the course of experimenting with the physics included in the models we made two modifications to the RH code. Those were:

1. (A)

Fixing input populations for specified atomic, or molecular, species, so that the radiation transport problem is solved using populations obtained from another source. We used this to include the RADYN NEQ hydrogen populations with non-thermal collisional excitation and ionisations also considered. With this method the NEQ hydrogen opacity was used by RH and the resulting radiation was output (as opposed to being treated as background only); 2. (B)

Inclusion of TR and coronal irradiation. RH can now accept a grid of emissivities as a function of temperature, electron density, and wavelength. In each grid cell above a defined temperature threshold ($50$ kK in our case) the local temperature and density is used to interpolate the emissivity at each wavelength modelled by RH’s atomic models. This is is then integrated through depth to obtain the intensity under the assumption that it is optically thin, and is subsequently injected as a downward-directed source of radiation from the apex of the loop. We provided every transition in the CHIANTI atomic database, with the transitions already solved by RH and certain other species removed, but the module is written such that a user can provide any grid they wish.

In summary, our experiments showed that:

1. (i)

PRD is required for modelling Mg ii in flares. The enhanced density in stronger flares does reduce, somewhat, the magnitude of the differences between CRD and PRD, but the wing intensities are still significantly overestimated with CRD. Features in the line wings due to redistribution effects also result in differences to the profile shapes that may be misinterpreted as being due to mass flows. 2. (ii)

Angle-dependent PRD is more accurate than the hybrid-PRD treatment when velocities are large. The profile shapes are unchanged, just the intensity in localised parts of the line. Given the substantial increase in computational time, we conclude that H-PRD is appropriate for most flare studies, so long as the caveats are understood. 3. (iii)

It is not necessary to model Mg i if one is only interested in the h & k lines, and subordinate lines, but there are intensity and formation differences in the far wings if Mg i is omitted. Therefore, if the (quasi-) continuum observed by IRIS is of interest, then Mg i should be included with Mg ii. 4. (iv)

The dominant sources of opacity at the Mg ii resonance line cores and near wings are are Mg ii and hydrogen, with hydrogen only really impacting the line wings. Only those species need be included in NLTE, with other sources safely treated in LTE. If the far wings or quasi- continuum are of interest then other species (e.g. Mg i, Si i) should be included. 5. (v)

Using the NEQ hydrogen populations did result in a more accurate modelling of the opacity in the Mg ii line wings, so that if possible those populations should be used. However, the differences were not very large. If it is not possible to include NEQ hydrogen populations for some reason then the results will not be very different. 6. (vi)

Using a three-level-with-continuum atom was found to be insufficient. A model atom containing excited states capable of supporting recombinations with cascades down to the h & k upper levels must be used. 7. (vii)

Lower-mid transition region irradiance does have an impact on the Mg ii formation, albeit with the effect seemingly very dependent on the transitions and temperatures considered. Irradiation from locations $T>50$ kK resulted in Mg ii line intensity differences of $\sim 10-30$ %. However, if irradiation only from locations $T>100$ kK is included then there is negligible impact on Mg ii. Removing Si i-iv and C i-iii from the irradiating spectrum reduces the magnitude of the differences. So, irradiation from the lower-mid TR should potentially be included, but this requires further investigation.

When modelling Mg ii in flares the above findings should be considered, and a suitable setup used. These experiments all assumed SE for the formation of Mg ii. Paper II in this series will discuss NEQ effects, which we found do have an impact during the initial heating and cooling phase of the flares (Kerr et al., 2019a).

Acknowledgments: GSK was funded by an appointment to the NASA Postdoctoral Program at Goddard Space Flight Center, administered by USRA through a contract with NASA. JCA acknowledges funding support from the Heliophysics Supporting Research and Heliophysics Innovation Fund programs. 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. We thank Dr. Han Uitenbroek for making the RH code publicly available, and Dr. Adrian Daw and Dr. Don Schmit for useful discussions regarding coronal irradiation.

Appendix A Ly $\alpha$ PRD

Here we show Ly $\alpha$ line profiles and hydrogen populations, for RADYN flare snapshots forward modelled using RH, comparing the PRD and CRD solutions. This was done for all three flares (F9, F10, and F11). In the weaker flare PRD effects are still clear, with CRD overestimating the line wing intensity and underestimating the emission peaks. The atomic level populations and proton density also vary somewhat. For the moderate and strong flares, where the electron density was significantly more enhanced, the PRD solution approaches the CRD solution, and atomic level populations show only modest changes.