Solar flares are transient, broadband brightenings to the Sun’s radiative output following the liberation of a tremendous amount of energy (up to 10³² erg, or larger: Emslie et al., 2012; Aschwanden et al., 2015) during magnetic reconnection (e.g., Priest and Forbes, 2002; Shibata and Magara, 2011; Janvier et al., 2013; Emslie et al., 2012). It is thought that this energy is subsequently transported predominately in the form of non-thermal particles. We primarily consider non-thermal electrons ¹ , accelerated during the reconnection process. Once they reach the denser lower solar atmosphere they thermalise via Coulomb collisions (e.g., Brown, 1971), heating and ionising the plasma and generating mass flows: chromospheric evaporation (upflowing material) and chromospheric condensations (downflowing material). Alternative mechanisms of energy transport in flares include non-thermal protons or heavier ions, thermal conduction following direct heating of the corona, and Alfvénic waves, discussed in more detail in Section 3.

There is unambiguous evidence for the presence of non-thermal particles in flares, due to the hard X-rays they produce via bremsstrahlung. Their ubiquitousness and the close spatial and temporal association with other flare emission (e.g., optical and UV) has bolstered the ‘electron-beam’ model of solar flares. Observations of hard X-rays, e.g., from the Reuven Ramaty High Energy Solar Spectroscopic Imager (RHESSI; Lin et al., 2002), can be used to infer the underlying non-thermal electron energy distribution, that itself can drive models of solar flares. There is a substantial body of literature describing the various characteristics of flares, and the means in which we observe them. I direct the reader towards the following reviews of flare observations, and flare particle acceleration and thermalisation: Benz (2008); Fletcher et al., 2011; Holman et al., 2011; Kontar et al., 2011; Zharkova et al., 2011; Milligan (2015). The bulk of the flare radiative output originates from the chromosphere and transition region, making those regions important areas of study for their diagnostic potential regarding the plasma response to energy injection, and the energy transport and release process themselves. However, speaking candidly, this potential has been somewhat squandered by the lack of routine high spatial-, temporal-, and spectral-observations of the chromosphere and transition region at UV wavelengths during flares (crucially, we lacked routine imaging spectroscopy of the flare chromosphere). That observational gap has fortunately been plugged following the launch of the Interface Region Imaging Spectrograph (IRIS; De Pontieu et al., 2014); in 2013, that now gives us an unprecedented view of the flaring chromosphere and transition region, yielding crucial new insights. Given the complex environment of these particular layers, parallel efforts to forward model the flaring lower atmosphere, and its impacts on the flaring corona, are required to make substantial progress in understanding the physics at play in flares.

This is the second paper in a two part review of how solar flare loop models in concert with IRIS observations have improved our understanding of solar flares. Between both parts I hope to emphasise that it is only by attacking the problem of flare physics via the combination of high quality observations and state-of-the-art models, that include the pertinent physical processes, that we can make rapid progress. Overall I aim to show: 1) how modelling has helped interpret the IRIS observations; 2) how IRIS observations have been used to interrogate and validate model predictions; and 3) how, when models fail to stand up to the stubborn reality of those observations, IRIS has led to model improvements. In Paper 1 of this review (Kerr, 2022) I provided a detailed overview of each numerical code, and discussed what we have learned from the study of Doppler motions from IRIS in the context of the non-thermal electron beam driven flare model. Also in Paper one is a more extensive introduction to solar flares. Here in Paper 2 I demonstrate how we have used the combination of IRIS and flare loop modelling to learn about plasma properties and flare energy transport mechanisms, and provide some thoughts on future directions.

IRIS is a NASA Small Explorer mission that has observed many hundreds of flares, including dozens of M and X class events. Both images (via the slit-jaw imager, SJI) and spectra (via the slit-scanning spectrograph, SG) are provided in the far- and near-UV (FUV and NUV), with a spatial resolution 0.33–0.4 arcseconds. High cadences are achievable, as low as 1 s but more generally a few seconds to tens of seconds. The strongest lines observed are Mg ii h 2,803 Å and k 2,796 Å lines (chromosphere), C ii 1,334 Å and 1,335 Å and Si iv 1,394 Å and 1,403 Å (transition region), and Fe xxi 1,354.1 Å ( ∼ 11 MK plasma), with numerous other lines observed in the three passbands [1,332–1,358] Å [1,389–1,407] Å, and [2,783–2,834] Å. The spectral resolution is ∼ 53 mÅ in the NUV and ∼ 26 mÅ in the FUV. The SJI observes at 2,796 ± 4 Å (Mg ii k), 2,832 ± 4 Å (Mg ii wing plus quasi-continuum), 1,330 ± 55 Å (C ii), and 1,400 ± 55 Å (Si iv). See De Pontieu et al., 2021 for a review of the various successes over the first near-decade of IRIS observations.

The models employed to study flares are generally field-aligned (1D) numerical codes (though there are some exceptions, e.g., the 3D radiative magnetohydrodynamic, MHD, model of Cheung et al., 2019). These codes are nimble enough to be run on timescales that make performing parameter studies of flare energy transport processes a tractable activity, and they allow us to include the relevant physical processes at the required spatial resolution (down to sub-metres) that is not yet feasible in 3D RMHD simulations. I focus on the RADYN, HYDRAD, FLARIX, and PREFT models here. A brief overview is presented below but see Paper 1 for a full description of each code.

The hydrodynamic field-aligned codes HYDRAD (Bradshaw and Mason, 2003a; Bradshaw and Mason, 2003b; Reep et al., 2013; Reep et al., 2019), RADYN (Carlsson and Stein, 1992; Carlsson and Stein, 1997; Carlsson and Stein, 2002; Abbett and Hawley, 1999; Allred et al., 2005; Allred et al., 2015), and FLARIX (Kas̆parová et al., 2009; Varady et al., 2010; Heinzel et al., 2016) are now well established and widely used by the solar flare community. These codes solve the equations describing the conservation of mass, momentum, charge, and energy in a single field-aligned magnetic strand rooted in the photosphere and stretching out to include the chromosphere, transition region, and corona. HYDRAD and RADYN use an adaptive grid where the size of the grid cells can vary to allow shocks and steep gradients in the atmosphere to be resolved as required (with HYDRAD varying the number of grid cells also), while FLARIX uses a fixed, but optimized, grid with ∼ 2000 points. The codes have various similarities and differences as regards treatment of radiation and flare energy transport, and with the numerical approaches themselves.

All three simulate the response of the atmosphere to injection of energy, typically via a beam of non-thermal electrons (but flare-accelerated ions can be included too). RADYN uses a Fokker-Plank treatment to model the evolution of the non-thermal electron distribution as a function of time (including return current effects), that was recently updated to use the standalone state-of-the-art non-thermal particle transport code FP ² (Allred et al., 2020). HYDRAD uses the analytic treatment of Emslie (1978) and Hawley and Fisher (1994), and FLARIX uses a test-particle module that provides the time-dependent beam propagation including scattering terms. Dissipation of Alfvénic waves has also been implemented in both HYDRAD and RADYN (Reep and Russell, 2016; Reep et al., 2018b; Kerr et al., 2016), and all codes can include ad hoc time dependent heating.

Each code has been conceived and developed to focus on particular details of the flaring plasma physics problem. RADYN and FLARIX are radiation hydrodynamic codes which couple the hydrodynamic equations to the non-LTE (NLTE) 1D radiative transfer and time-dependent non-equilibrium atomic level population equations, for elements important for chromospheric energy balance. RADYN considers H, He and Ca, with Mg also sometimes included, whereas FLARIX considers H, Ca, and Mg (with plans to update the code to include He). Continua from other species are treated in LTE as background metal opacities. Optically thin losses are included by summing all transitions from the CHIANTI atomic database (Dere et al., 1997; Del Zanna et al., 2015; Del Zanna et al., 2021) ³ apart from those transitions solved in detail. Additional backwarming and photoionisations by soft X-ray, extreme ultraviolet, and ultraviolet radiation is included. Both currently use the assumption of complete frequency redistribution (CRD) ⁴ when solving the radiation transport problem, so that post-processing via other radiation transport codes such as RH (Uitenbroek, 2001), RH15D (Pereira and Uitenbroek, 2015), or MALI (Heinzel, 1995) is required. In RADYN and FLARIX the loop is modelled as one leg of a symmetric flux tube. RADYN also allows to calculate aposteriori (i.e. with no feedback on the plasma equations of mass, momentum, and energy) the time-dependent non-equilibrium populations and radiation transport of a desired ion via the minority species version of that code, MS_RADYN (Judge et al., 2003; Kerr et al., 2019b; Kerr et al., 2019c).

HYDRAD does not solve the detailed optically thick radiation transport and atomic level population equations, instead employing approximations of chromospheric radiation losses. Losses from H, Ca and Mg are included via the approach of Carlsson and Leenaarts (2012). The code has also recently adopted a more accurate method for computing NLTE H populations following the prescription of Sollum (1999) which approximates the radiation field in the chromosphere (Reep et al., 2019). Ion population equations, however, are solved self-consistently in full non-equilibrium ionization (NEI) for any desired element, returning a more accurate calculation of the optically thin radiative losses and spectral synthesis of optically thin lines using those ion fractions. While the treatment of optically thick radiation is less robust than in RADYN or FLARIX, HYDRAD has the advantage of being significantly less computationally demanding. Other important differences are that HYDRAD features a multi-fluid plasma that treats the electron and hydrogen temperatures separately, it solves a full length flux tube (foot-point to foot-point) of arbitrary geometry (e.g., based on a magnetic field extrapolation) and includes effects due to cross-sectional area expansion (varying inversely with the magnetic field strength), which has been shown to play an important role in dynamics (Reep et al., 2022).

PREFT (Guidoni and Longcope, 2010; Longcope and Guidoni, 2011; Longcope and Klimchuk, 2015; Longcope et al., 2016) is a rather different code than the other three, and is a powerful tool to study the impact of magnetic loop dynamics during flares. It is a 1D MHD code that solves the thin flux tube (TFT) equations. The tube is initialized at the instant after a localized reconnection process within the current sheet has linked sections of equilibrium tubes from opposite sides of the current sheet. No further reconnection occurs, and any heating from the initializing event is neglected. In its subsequent evolution, the tube retracts under magnetic tension releasing magnetic energy and converting it to bulk kinetic energy in flows which include a component parallel to the tube. The collision between the parallel components generates a pair of propagating slow magnetosonic shocks, which resemble gas dynamic shocks as they must in the parallel limit. Radiative losses are optically thin, and normally an isothermal, but gravitationally stratified, chromosphere is included mostly as a mass reservoir. Solutions of the TFT equations show that thermal conduction carries heat away from the shocks, drastically altering the temperature and density of the post-flare plasma (Longcope and Guidoni, 2011; Longcope and Klimchuk, 2015; Longcope et al., 2016; Unverferth and Longcope, 2020).

The flare chromosphere has been studied extensively using optically thick lines, which while presenting challenges with respect to extracting useful information due to their complex formation properties, offer important diagnostics of how the plasma responds to flares. Most notably, the Hα line has been both observed and modelled in flare studies too numerous to exhaustively list or describe in detail here. Some examples include: Canfield and Ricchiazzi (1980), Canfield et al., 1984, Canfield and Gayley (1987), Canfield et al. (1991), Heinzel (1991), Gayley and Canfield (1991), Gan et al., 1992, Gan et al., 1993, Li et al., 2006, Kuridze et al., 2015, Rubio da Costa et al. (2015). Some important numerical results inform us about how conditions in the upper chromosphere result in varying characteristics of the line, such as the depth of central reversal, are related to plasma conditions. The coronal pressure, for example, was found to be an important factor in determining the depth of the central reversal, with a high coronal pressure ( > 100 dyne cm⁻²) required to markedly reduce the depth or fill in the reversal (e.g., Canfield et al., 1984). The pressure is closely related to density at the formation region of Hα, such that increasing the pressure forces the transition region, and hence Hα formation height, deeper. As we will see in this section the Mg ii central reversal depth is also seemingly related to conditions in the upper chromosphere. Additionally, exploring line widths and intensities could help constrain densities (via Stark wings) and temperatures (relative intensities of lines). Other lines have received significant attention in the past have included the Ca ii H and K resonance lines, and the infrared triplet (e.g., Machado and Linsky, 1975; Fang et al., 1992; Falchi and Mauas, 2002; Kuridze et al., 2018; Ding, 1999). Non-thermal processes have also been studied numerically using chromospheric spectral lines. Non-thermal collisions between the particle beam and hydrogen or calcium have been shown to be a significant process, affecting the intensity and shape the lines (e.g., Fang et al., 1993; Druett and Zharkova, 2018; Kas̆parová et al., 2009), and charge exchange between ion beams and hydrogen, producing highly redshifted Lyman line emission, has been suggested as a means to diagnose the non-thermal proton distributions in the lower atmosphere (e.g. Orrall and Zirker, 1976; Canfield and Chang, 1985).

In this section I introduce the Mg ii NUV spectra as observed by IRIS, which being one of the strongest lines in the IRIS passbands has become a workhorse for studying the chromosphere, including during flares. Like the spectral lines of hydrogen and calcium, their characteristics are very sensitive to atmospheric properties, with various flare effects changing both the formation location as well as local plasma conditions. Magnesium is 18 times more abundant than Calcium, consequently forming higher in altitude and sampling the upper chromosphere, a useful region for understanding energy deposition during flares. The Mg ii h and k Doppler width is much smaller than that of Hα, offering advantages in sensitivity to both Doppler motions and turbulence.

Prior EUV observations of hot flare lines have shown anomalously broad lines, of unknown origin (e.g., Milligan, 2011; Milligan, 2015). While several suggestions have been made, a definitive solution remains elusive (as you have no doubt realised by now, line widths are a sore spot for flare modellers). As discussed in Paper 1 in relation to probing the duration and magnitudes of chromospheric evaporation, the Fe xxi 1,354.1 Å line observed by IRIS offers a window at high spatial, temporal, and spectral resolution on hot flare plasma at ∼ 11 MK. Here I discuss a few studies that attempted to model Fe xxi 1,354.1 Å emission in flares, concentrating on line broadening. This is a different problem than the Mg ii missing width, as opacity effects play no role for Fe xxi 1,354.1 Å emission, which is a forbidden line.

The Fe xxi 1,354.1 Å line has been observed in numerous flares by IRIS (e.g., Graham and Cauzzi, 2015; Polito et al., 2015; Polito et al., 2016; Young et al., 2015; Tian et al., 2015). It is observed to be largely symmetric, with significantly enhanced line widths. It is initially weak and broad, and becomes more narrow and intense over time. The line widths during flares have ranged from the instrumental + thermal width W ∼ 0.43 Å (assuming ionisation equilibrium) at loop tops to W ∼ [0.5–1] Å or larger in ribbons. Kerr et al., 2020 performed a superposed epoch analysis similar to Graham and Cauzzi (2015)’s Doppler shift analysis, to understand the typical evolution of Fe xxi line widths over time in the 2014-September-10th X-class flare. That event showed a large amount of scatter during the impulsive phase of the flare, but with W ∼ [0.6–1.2] Å, peaking after t ∼ 200 s, before gradually returning to pre-flare values over the subsequent 500 s or so.

A popular suggestion for the origin of the broad line profiles of hot lines is the superposition of flows along the line of sight from numerous Doppler shifted line components. Using the RADYN code, Polito et al., 2019 explored this idea. They produced several field-aligned flare simulations, with a t = 60 s heating duration. Synthetic Fe xxi emission was produced, with Doppler shifts applied as appropriate as a function of height along the loop. From those simulations they constructed both single and multi-stranded loop models, which for the latter had 100 identical threads each with a randomly selected start time within 15 s of the first thread start time. Other loop lengths and energy fluxes was also tested, but no change to the overall conclusions was found. The threads were then orientated in several ways that investigated the effects of loops being co-spatial or along a ribbon-like structure within 20 IRIS pixels, either aligned in the same angle or at different orientations. Emission from each set-up was summed to mimic different scenarios of IRIS looking through different lines of sight. Polito et al., 2019 found that there was a non-negligible asymmetry, with an anti-correlation between broadening and asymmetry (broader lines were more asymmetric), in contrast to observations which showed largely symmetric lines no matter the width, in each of their scenarios. Narrow profiles were quite symmetric, and were largely due to superposition of several upflows that had similar magnitudes, but the superposition of loops was unable to characterise the broad, symmetric Fe xxi profiles in this case. Figure 4 shows a sample of these experiments, for a single loop model with different orientations. The synthetic spectra and resulting broadening are shown, where it is clear that asymmetries are present.

Kerr et al., 2020 also studied synthetic Fe xxi line widths using RADYN simulations. The model framework they developed, RADYN_Arcade, is described in Paper 1. From the field-aligned loops grafted onto the observed magnetic field skeleton, the superposition along the line of sight, and loop geometry was automatically taken into account. Qualitatively they produced synthetic Fe xxi emission that largely followed the observations. There was line broadening that was strongest at flare footpoints, and which narrowed along the loops towards looptops. However, while some profiles exceeded W ∼ 0.8 Å the majority of the Fe xxi lines only reached W ∼ [0.5–0.6] Å, much narrower than observed. There was not a very strong correlation between asymmetry and line width (broad profiles could be both symmetric or asymmetric), but the largest asymmetries were associated with the broadest profiles. The majority of the profiles were fit well with a single Gaussian, with only a subset requiring multiple components. There was also not a strong correlation between line width and Doppler shift, in contrast to some observations of hot flare lines studied by Milligan (2011). Finally, a synthetic superposed epoch analysis showed again a qualitative similarity to observations, but with profiles that were too narrow (with synthetic FWHM ∼ [0.4–0.8] Å, clustered around FWHM ∼ 0.5 Å, compared to observations with a range of FWHM ∼ [0.4–1.4] Å, clustered around FWHM ∼ 0.8 Å), and profiles that both peaked and narrowed in width too quickly (in observations the profiles took several hundred seconds to narrow towards pre-flare values, compared to t < 100 s in the model). Figure 5 shows a map of Fe xxi from the RADYN_Arcade model at various snapshots, alongside the synthetic superposed widths, and the asymmetries versus line widths, illustrating the discrepancies. These results agree with the earlier findings of Polito et al., 2019.

As set out nicely by Polito et al., 2019, there are several possible physical conditions that the RADYN modelling did not account for, which might explain how to obtain a closer match to the IRIS observations. Turbulence (including MHD wave turbulence) could broaden lines symmetrically. In fact, Allred et al., 2022 recently demonstrated that by suppressing thermal conduction in a RADYN simulation, via non-local effects or turbulence (e.g. Emslie and Bian, 2018), could lengthen the gradual phase of a flare, and produce a flow pattern more consistent with observations (e.g. Milligan and Dennis, 2009). Taking the turbulent mean free path of the best fit model, Allred et al., 2022 were able to estimate the broadening associated with turbulence for numerous lines (including hot lines from high charge states of Fe), finding that they were broadened substantially and symmetrically. There are plans to perform follow on studies to the modelling of Kerr et al., 2020, using these RADYN updates, and with the observed line widths from IRIS and the Hinode/EUV Imaging Spectrograph (EIS) observations as constraints on the degree of suppression to include. Using PREFT, Dr. William Ashfield and Dr. Dana Longcope are exploring the creation of MHD turbulence following loop retraction with added drag (private communication 2022). I eagerly await the application of their modelling to the formation of Fe xxi. Another source of broadening could be non-equilibrium effects, such that the ion temperature is very much larger than the equilibrium value of 11.2 MK. This would require an ion temperature on the order 40–60 MK, which also likely requires decoupling of the ion and electron temperatures (see also de Jager, 1985; Polito et al., 2018a). Given the high densities in flare footpoints, it is not clear if such extreme non-equilibrium processes apply, but the HYDRAD code is the ideal resource to study this in flares. Finally, Alfvénic waves propagating downwards from the magnetic reconnection site could broaden spectral lines via ion motions (see Section 3).

The extreme gradients through the transition region make it an important interface for mass and energy transport during flares. Strong lines produced in the transition region that are observed by IRIS are the Si iv and C ii resonance lines. Aside from study of their Doppler shifts, C ii has been relatively little studied in flare loop models. Si iv 1,394 Å and 1,402 Å, however, have been modelled in a few studies. I discuss their Doppler shifts in Paper 1, but here focus on their intensity ratio, and what that might tell us about temperatures and densities in the flaring transition region.

These lines have been used to infer flows in flares from their Doppler shifts, under the assumption that they are optically thin. They increase in intensity, broaden to a line width of similar magnitude to Mg ii, and exhibit red wing asymmetries indicative of mass flows on the order of a few ×10–100 km s⁻¹ (e.g., Tian et al., 2015; Li et al., 2017; Yu et al., 2020). Non-Gaussian line shapes have been observed in flare ribbons, but most observations suggest optically thin formation from the ratio of the 1,394 Å/1,402 Å intensities (R _1394∕1402 = 2; though it is often the case that only one of the lines is included in the IRIS lineslists), with exceptions discussed later in this section. Given the optically thin assumption, most flare modelling of this line was performed by computing the emissivity from CHIANTI atomic data alongside the stratification of physical variables from flare loop model atmospheres, and summing through height in some fashion to obtain the total intensity. However, some quiescent Sun studies suggested that effects of non-equilibrium radiation transfer, photoexcitation, or charge exchange may be important in setting the Si ion fraction stratification (e.g., Dudík et al., 2017; Dzifc̆áková et al., 2017; Dzifc̆áková and Dudík, 2018; Kerr et al., 2019c). Observations of line ratios in stellar flares have also suggested that the resonance lines of Si iv, and also of C iv, exhibit opacity effects (e.g., Bloomfield et al., 2002; Mathioudakis et al., 1999) and form under optically thick conditions, or, at least with some non-negligible optical depth τ > 0.1.

To determine the importance of radiation transfer effects on the formation of Si iv resonance lines during flares, Kerr et al., 2019c simulated a large number of electron beam driven flares using RADYN, and then generated the synthetic Si iv emission in two ways. The first was the standard optically thin synthesis using CHIANTI contribution functions and ionisation fractions, assuming equilibrium. The second was to use the minority species version of RADYN, MS_RADYN, to synthesise the two Si iv lines including the effects of photoionisation and photoexcitation, non-equilibrium ionisation, opacity, and charge exchange. MS_RADYN takes as input certain hydrodynamic variables at the cadence of RADYN′s internal timestep (i.e. the relevant timescales to capture dynamics, not simply the output cadence), and then solves just the NLTE non-equilibrium radiation transfer for a given minority species. That is, there is no feedback of the radiation from the minority species on the atmosphere itself. The line profiles from each method were quite different, even from the pre-flare where charge exchange broadened the Si iv ion fraction stratification, which peaked slightly cooler in temperature than it does in ionisation equilibrium (T ∼ 66 kK versus T ∼ 80 kK). Charge exchange is generally not considered in transition region line modelling, but both the results of Kerr et al., 2019c, and recent quiet Sun modelling of Dufresne et al. (2021a), Dufresne et al. (2021b) highlight its importance. Weaker simulated flares generally produced similar results from both synthesis methods. However, in stronger flares they differed. The peak intensities of MS_RADYN Si iv profiles were smaller, but the lines broader overall due to opacity effects so that the integrated line intensities were higher. Those profiles also showed stronger asymmetries, and self-absorption features. Crucially, the intensity ratio deviated from the optically thin ratio of R _1394∕1402 = 2. Opacity effects were present in all simulations with an energy injection F > 5 × 10¹⁰ erg s⁻¹ cm⁻², and for some weaker flares with softer electron distributions since they more easily heated the upper chromosphere and lower transition region. Some of these flares only exhibited opacity effects for a transient period, since the transition region and upper chromosphere compressed quickly, meaning there was not a sufficient column mass of Si iv to build up opacity. When there was an extended flaring lower-transition region (that is temperatures climbing through 30 kK < T < 100 kK over a large height range before sharply rising to MK coronal temperatures), opacity effects were present. Roughly, opacity effects were present when the temperatures were enhanced to 40 < T < 100 kK above a column mass 5 × 10^–6 g cm⁻².

There have since been a number of observations of R _1394∕1402 deviating from the optically thin limit R _1394∕1402 = 2 (e.g., Mulay and Fletcher, 2021; Zhou et al., 2022). Mulay and Fletcher (2021) found R _1394∕1402 ≠ 2 at several locations along flare ribbons in an M7.3 flare. Zhou et al., 2022 report similar results, noting also that the ratio varies across the line profile, with stronger opacity in the core so that photons scattered from an optically thick core can easily escape through optically thin line wings. In those observations, we might infer that the flaring atmosphere produced the extended region of 40 < T < 100 kK at sufficiently high density. It is important to note that if we do not see much observational evidence for these potentially short lived opacity effects, then our models may be predicting too much density at these intermediate temperatures. Further RT modelling, particularly of other transition region lines in conjunction with Si iv is sorely required, as are high cadence observations to catch potentially transient opacity effects. Another impact of potential opacity effects in transition region lines, and motivation for their further study, is that these are important contributors to the (assumed) optically thin radiative loss functions, which are a major component energy loss in the simulations, governing the plasma response in our models.

Panos et al., 2021 and Panos and Kleint (2021) explored, using machine learning techniques based around mutual information theory (MI; Li, 1990) ¹³ , the correlation between the various spectral lines of IRIS through the transition region and chromosphere. I encourage the reader to read their detailed analysis carefully, in particular the subtleties surrounding selecting flaring areas and how this might impact correlations, but note the headline results here. They find weak correlations between spectral line pairs during quiescent periods, but substantially enhanced correlations of those pairs during solar flares. Mg ii and C ii have the strongest correlation, followed by their correlations with Si iv. Other lines (e.g. O iv) are more weakly correlated, and others such as Fe ii only show strong correlations directly over flare ribbons. This coupling meant that Panos and Kleint (2021) were able to predict the most probable spectrum of a certain IRIS observable given an input Mg ii spectra, for example. The strong correlation of Si iv to the chromospheric lines, despite the weak correlation of other transition region lines such as O iv, could be due to the deeper formation height suggested during the Kerr et al., 2019c simulations. Further, the coherency that flares introduce could be a result of the strong compression of the chromosphere and transition that occurs in many flare simulations. For example, Figure 11 in Kerr et al., 2019c shows that over time the range of formation height of the IRIS line cores can shrink to a very small Δz. The “big data” studies of Dr. Panos and collaborators provide an excellent test bed against which models can be critiqued—our models should be able to produce similar coherency between the various lines observed by IRIS, and this should be a target of our efforts in the near future.

Finally, I note briefly that it is typical in flare simulations from HYDRAD, RADYN, and FLARIX to produce large enhancements in electron density through the chromosphere and into the corona. These can be in excess of n _e > 10^13–14 cm⁻³ in the chromosphere, and n _e > 10^10–12 cm⁻³ through the transition region and lower corona. Indeed, as discussed in the preceding sections, a very large electron density at the Mg ii formation temperatures is required to explain the single peaked profiles. IRIS and Hinode/EIS density sensitive lines from the corona and transition region can demonstrate if these densities are consistent with observations. Polito et al., 2016 measured the ratio of the O iv 1,399.77 Å and O iv 1,401.16 Å line pair, which form at T ∼ 158 kK, during the impulsive phase of the X2-class flare that occurred on 2014-October-27th. The ratio reached the high-density limit, indicating that the density of the flare transition region reached n _e > 10¹² cm⁻³. A caveat here is the assumption of ionisation equilibrium, so that the observed ratio may be in part due to non-equilibrium effects. Other assumptions are that the lines are free of unknown blends, and that the plasma is a Maxwellian, which may not be the case in solar flares, or even active regions, which have been seen to exhibit κ distributions (e,g, Jeffrey et al., 2016; Dzifc̆áková et al., 2018; Del Zanna et al., 2022). Similar analysis using EUV spectral lines from EIS indicated a coronal density at 2 MK of n _e > 10^10–11 cm⁻³. Polito et al., 2016 then modelled this flare using HYDRAD, finding that the electron density in the synthetic flaring atmosphere (both flare footpoints and the transition region/lower corona) were consistent with the observationally derived values.

Pivoting slightly to white light observations, the IRIS NUV Balmer continuum modelling and observations seem to suggest that there is likely some contribution to the optical continuum excess in flares from recombination radiation in the upper chromosphere (see Section 3.3). Given the dependence of bound-free (and also free-free) emission on electron density, NUV and optical continuum observations of the chromosphere can give us some means to investigate the density there, and off-limb observations allow us to isolate the chromospheric portion. Off-limb observations of white light flares have revealed both the typical footpoint sources at the base of flare loops as well as bright loop structures (also referred to as prominence loop systems in some literature). The former was discussed by Heinzel et al., 2017, who analysed SDO/HMI continuum data of off-limb flares that revealed co-spatial HMI 6173 Å and RHESSI hard X-ray emission, with a characteristic height of ∼ 1000 km (see also Krucker et al., 2015). Using an analytical argument of the relative strength of white-light continuum emission mechanisms, Heinzel et al., 2017 determined that for electron densities above n _e > 10¹² cm⁻³ Balmer bound-free recombination emission dominated over Thomson scattering of incident radiation from the solar disk, with some contribution from free-free emission. Further comparisons using FLARIX electron-beam driven flares confirmed their supposition, with the simulations containing a high electron density (n _e > 10^12–13 cm⁻³) at the height range of the observed HMI emission, and with an intensity as a function of height resembling the HMI observations (with some assumed loop thickness). These results are consistent with the IRIS NUV continuum flare footpoint observations. To my knowledge no off-limb IRIS NUV flare observations have been reported, and it is likely they would be fairly weak unless a long exposure time was used, but those would be very interesting to compare to the HMI sources. Looptop structures that are readily apparent in optical continuum observations pose a bit more of a challenge to models to reproduce, namely due to the very high coronal densities they imply. Several studies of very strong flares have inferred looptop electron densities between n _e = 10¹²–10¹³ cm⁻³, usually in the gradual phase of flares, presumably once loops have cooled. For example, Hiei et al., 1992 studied both footpoint and loop sources in the 16 August 1989 flare that was estimated to be an X20 class event. They predicted the intensity of emission from Thomson scattering, free-free, or recombination radiation for a range of temperatures given an assumed emitting volume, inferring from the observed intensity that n e = 5 × 1 0 11 , 2 × 1 0 12 , 1 × 1 0 13 cm⁻³ for T = [10⁴, 10⁵, 10⁶] K, respectively, in a loop source several hundred km above the white-light footpoint sources. A similar analysis was performed by Jejčič et al., 2018 using HMI observations from the X8.2 10 September 2017 event, finding n _e = 10¹²–10¹³ cm⁻³ were the most likely values in a large parameter space of temperatures and emitting thicknesses. Inverting Ca ii 8,542 Å and H β data taken by the Swedish Solar Telescope of that same flare, Koza et al., 2019 found consistent values in the cool loops. By studying polarisation of HMI data from the X2.8 flare on 13 May 2013, Saint-Hilaire et al., 2014 determined that the emission could not be solely due to Thomson scattering and estimated an electron density in the range n _e = 3.5 × 10¹¹–1.8 × 10¹² cm⁻³. While flare models can readily explain electron densities up to a few × 10¹¹ cm⁻³ in the upper portion of the corona, obtaining higher electron densities at looptops is less straightforward and demands an explanation.

First proposed as a means of heating the temperature minimum region where non-thermal electrons likely could not reach, but which observational evidence suggested experienced a modest temperature rise in flares, Emslie and Sturrock (1982) constructed a simple but informative model of energy transport via downward propagating, coronally-generated, Alfvénic waves. In this model, waves would be produced from the reconnection site, propagating through the corona into the lower atmosphere to the temperature minimum region where they were damped by resistivity. These simulations assumed Mono-chromatic (single frequency) waves, employed the WKB approximation (that is, waves were not reflected by density gradients), and assumed an instantaneous travel time. These assumptions allowed a straightforward formulation of a damping length to model the dissipation.

This notion was revisited by Fletcher and Hudson (2008) who investigated the possibility that Alfvénic waves could not only deliver the energy liberated by magnetic reconnection to the chromosphere, and accelerate electrons in the corona via field-aligned electric fields but could also potentially locally accelerate electrons in the chromosphere via mode-conversion to high wave-numbers resulting in turbulent acceleration. More work is required to understand the role of these waves in particle acceleration. The thought experiments of Fletcher and Hudson (2008) explored Alfvénic waves as an alternative to the electron beam model as a means to deliver flare energy and explain observations of both hard X-rays and broadband enhancements of the UV/optical/infrared. This was motivated by perceived issues with the coronal acceleration problem, namely the vast numbers of electrons required ( > 1 0 36 elec s⁻¹), which can quickly deplete the coronal volume of ambient electrons. Return currents can resupply the corona with electrons, however, mitigating this problem.

Even if they are not required as a complete replacement to electron beams (which is still a source of vigorous debate), it is important that we continue to properly consider the role of Alfvénic waves in flares. Flares are, fundamentally, a violent restructuring of the magnetic field, meaning that MHD waves are undoubtedly produced. The question is, do they carry sufficient energy to play a non-negligible role in transporting energy compared to coronally accelerated electrons, and can they efficiently heat the chromosphere (either alongside or instead of those electrons). Additionally, we do not see hard X-rays all along the flare ribbons. Perhaps different parts of ribbons are heated by different mechanisms. Some MHD simulations by Russell and Fletcher (2013) and Russell and Stackhouse (2013) revealed that Alfvénic waves could penetrate the transition region density boundary if they had a high enough frequency, f > 1 Hz, meaning the WKB approximation could be used within loop models to further investigate high-frequency Alfvénic waves. They also noted that ion-neutral interactions were important, alongside electron resistivity, in damping the Alfvén waves.

Inspired by these results Reep and Russell (2016) modified HYDRAD to model Alfvén waves using the WKB approach of Emslie and Sturrock (1982), but with an updated treatment of damping which included ambipolar effects. Thus, the waves were damped by ion-neutral, neutral-electron, and electron-ion collisions. Modelling a range of Alfvén wave parameters, including the injected Poynting flux, mono-chromatic frequency, and wave number they found that they could strongly heat the chromosphere, and that they could drive explosive chromospheric evaporation. This model was further improved in HYDRAD by Reep et al., 2018b, to include the wave travel time, via ray-tracing so that the waves propagate at the local Alfvén speed. They show that in addition to certain wave parameters being damped more effectively in the lower atmosphere than in the upper chromosphere, that leading waves can effectively bore a hole through the chromosphere allowing following rays to penetrate deeper into the lower atmosphere. This occurred due to ionisation by the leading waves, reducing the local damping.

Following the approach of Reep and Russell (2016), Kerr et al., 2016 included Alfvén waves as a mode of energy transport into RADYN. This initial work employed the instant-travel approximation where the wave propagation was ignored. I have since updated RADYN to include the travel time of the wave in the same manner as Reep et al., 2018b. For the remainder of this section I mostly discuss the results using the Kerr et al., 2016 model, since those are the published results relevant to IRIS, but work modelling the IRIS observables including travel time is underway via both HYDRAD and RADYN.

Kerr et al., 2016 compared the atmospheric dynamics and radiative output of two RADYN simulations, 1) an electron beam, and 2) a mono-chromatic Alfvén wave. The energy flux of each was set to be 1 × 10¹¹ erg s⁻¹ cm⁻², and the Alfvén wave parameters set to most effectively heat the upper chromosphere. A magnetic field stratification was imposed for the purpose of defining the Alfvén speed and damping lengths; it did not evolve during the simulations. Two spectral lines were compared, the Ca ii 8,542 Å and Mg ii k line, the latter synthesised using RADYN atmospheres with RH. While the atmospheres showed some striking similarities in each model’s ability to heat the chromosphere and drive strong upflows, as was first seen in Reep and Russell (2016), there were intriguing differences in the chromospheric stratification. These differences revealed themselves in the spectral lines also.

The Alfvén waves produced a flatter, more spatially extended energy deposition profile compared to the electron-beam heating profile, resulting in temperature rises at deeper heights than the electron beam simulation. Despite this, it did take time for the electron density in the lower atmosphere to catch up to the electron beam simulation because of the absence of non-thermal collisional ionisation due to the beam itself. Runaway helium ionisation due to the more concentrated electron beam heating removed the 304 Å line as a radiator, resulting in a high temperature bubble forming, flanked by narrow cool high density regions. These flanking regions expanded as a high velocity upflow, and slower downflow (in addition to the initial explosive evaporation). While this did not form in the Alfvén wave simulation, a secondary upflow appeared in the Alfvén wave simulation also, but was more gentle with a shallower spatial gradient. In the electron beam simulation the Mg ii k line had a central reversal that was redshifted during most of the heating phase. The line wings had small optically thin contributions due to the flow patterns. In the Alfvén wave simulation, however, the line formed in a gentle upflowing region of the chromosphere, shifting the absorption profile strongly to the blue. Since the densities in the upflow were relatively weak this did not fully shift the line but instead pushed the core and blue k2v peak closer in formation height until they merged. The upflow produced optically thin contributions through the blue wing. Fewer absorptions due to shifting the absorption profile boosted the red k2r peak in comparison to the heavily suppressed k2v peak, meaning that the whole profile took on a very asymmetric form. The k line could be mistaken as being single peaked with a large blue wing asymmetry. Differences in the shape of the Mg ii k line cores were a direct result of the different flows, that themselves were due to the different stratifications of damping in either the electron beam or Alfvén wave energy transport mechanisms. Kerr et al., 2016 demonstrated that Mg ii can help discriminate between energy transport models, but much more work needs to be done here, particularly studying multiple IRIS spectral lines forming in a wider array of Alfvén wave driven simulations that include the wave travel time. The predictions from each model should also be compared to the k-means classifications of Panos et al., 2018, and we must work to improve the models to include a spectrum of wave frequencies, and to constrain the properties of the waves.

While Alfvén waves are certainly produced during magnetic reconnection they have been detected in situ in the magnetosphere, (e.g., Chaston et al., 2005; Wygant et al., 2002; Gershman et al., 2017), a vital question is how much energy do they carry to the lower atmosphere? Is it enough to compete with electron beams as an important contributor to the flare energy budget, or is it negligible and thus safely ignorable? Thus far, the simulations of Reep and Russell (2016), Reep et al. (2018b) and Kerr et al. (2016) injected a Poynting flux of the level that we know from electron beam driven flare simulations, and bolometric flare observations, is required to significantly heat the chromosphere. An observational constraint on the Poynting flux is required. An upper limit could be placed on this by investigating the width of lines formed at different temperatures (i.e., altitudes). The non-thermal component of the width could result from ion motion in response to an Alfvén wave. To demonstrate what the upcoming EUV observations from the Multi-Slit Solar Explorer (MUSE, scheduled for launch in 2026; De Pontieu et al., 2020) would reveal about solar eruptive events, many flare models synthesised MUSE observables and demonstrated how MUSE might discriminate between model predictions (Cheung et al., 2022). As part of that effort we modelled the broadening that would be induced due to an Alfvén wave propagating down the loops in our RADYN_Arcade model, noting that the line was indeed substantially broadened. This is demonstrated in Figure 6 which shows the RADYN_Arcade model before and after Alfvén wave broadening is included. Coordinated high spatiotemporal resolution observations between MUSE and the High Throughput EUV Solar Telescope (EUVST, also scheduled for a ∼2026 launch) could track the development of non-thermal widths during a flare, placing constraints on the Poynting flux. Knowledge of the coronal magnetic field would also be very advantageous here, to help set the Alfvén speed and damping lengths, and to determine the amplitude of magnetic field perturbations. Finally, it is worth noting that Alfvénic waves have been proposed as the mechanism responsible for the observed elemental fractionation between the photosphere and corona. Low-first ionisation potential (FIP; < 10 eV) species generally have a coronal abundance that is 4 or so times that of the photosphere. The ponderomotive force generated in MHD waves has been suggested as a potential cause of the so-called FIP effect (e.g., Laming, 2015). Observations of the FIP effect in flares (e.g., Doschek et al., 2018) may then shed light on the properties of Alfvénic waves produced during flare reconnection.

The energy flux injected to dynamic flare simulations has typically ranged on the order F = 10^9–11 erg s⁻¹ cm⁻², driven in part, admittedly, because of the computational expense and difficulty of injecting very much stronger values of F into time-dependent models until fairly recently (F > 10¹² erg s⁻¹ cm⁻² fluxes, while computationally demanding, are now possible). This range has been inferred from numerous studies of flares in both the RHESSI era and before, but we are now realising that in some of the strongest flare sources we may be underestimating F, perhaps by an order of magnitude in some cases! The physical rationale and implications behind this, with regard to non-thermal particle production and transport, are beyond the scope of this review, but an important factor in tying down the existence of very high beam fluxes are the modern observations at high spatiotemporal resolution of UV and optical flare sources. IRIS, Hinode/Solar Optical Telescope (SOT), and ground based observatories have revealed flare sources are smaller that typically assumed from older data (particularly so if looking at white light flare data) ¹⁴ . A detailed study comparing source sizes from flare ribbons observed by Hinode/SOT to hard X-ray imaging spectroscopy suggested that the beam flux may very well be F > 10¹² erg s⁻¹ cm⁻² in that flare (Krucker et al., 2011). Further, some groups have started looking at newly activated sources to define the areas into which energy is being injected within some observational window. Newly activated sources might be as small as to be on the order 10¹⁶ cm⁻² or below (Krucker et al., 2011; Sharykin and Kosovichev, 2014; Milligan et al., 2014; Kleint et al., 2016; Kowalski et al., 2017; Graham et al., 2020). IRIS observations can both guide the magnitude to inject based on high resolution observations of source areas, and act as a validation.

Kowalski et al., 2017 injected fluxes of F = [1 × 10¹¹, 5 × 10¹¹] erg s⁻¹ cm⁻² to simulate the two brightest sources in the 2014-March-29th X-class flare, focussing on modelling the NUV continuum response. These fluxes were guided by the hard X-ray analyses of Kleint et al., 2016 and Battaglia et al., 2015, with the range based on arguments of the continuum emitting areas identified by Kowalski et al., 2017. Portions of the NUV continuum in the region λ ∼ [2,814–2,832] Å, observed by IRIS, were first identified by Heinzel and Kleint (2014), who extracted patches of continua free from lines. They determined these line-free regions as likely being part of the Balmer continuum that remained optically thin during the flare and which formed in the mid-upper chromosphere. This means the continuum response would be very sensitive to the electron density throughout the flare chromosphere. Kowalski et al., 2017’s numerical experiments showed that the NUV continuum was too weak in the lower energy flux simulation, and much too weak in a set of experiments in which similar energy flux was instead deposited directly in the corona and allowed to conduct down to the chromosphere (potentially due to the lack of non-thermal collisional ionisations in the conduction-only simulations, though this was not commented on by the authors). In the high energy flux simulation the continuum did reach a sufficient level to match observations by t ∼ 2 s, peaked a few seconds later, before declining thereafter (but still remaining 100–200% above the pre-flare). Thus, a high energy flux was in fact required to produce conditions to raise the continuum intensity to the observed level. An in-depth analysis found that the NUV continuum was formed by hydrogen recombination emission from two distinct layers, both optically thin: a stationary chromospheric layer and a dense condensation that rapidly forms and accrues mass. As time progressed the condensation became responsible for the bulk of the emission, due to the fact that as the density increased an increasing proportion of the non-thermal electrons thermalised in the condensation itself, and consequently the stationary layer cooled somewhat. The conditions inside this condensation were found to be comparable to those of earlier slab model explanations of Balmer continuum enhancements (e.g. Donati-Falchi et al., 1985), suggesting that condensations (and high beam fluxes) are required to explain the brightest continua enhancements.

Various effects should be accounted for when considering very large non-thermal electron flux densities, such as the beam-neutralising return current including the effects of runaways (Zharkova and Gordovskyy, 2005; Holman, 2012; Allred et al., 2020; Alaoui and Holman, 2017; Alaoui et al., 2021), and instabilities that affect the beam propagation (e.g. Hannah et al., 2009; Lee et al., 2008; Li et al., 2014). A discussion of those is included in Kowalski et al., 2017, but are beyond the scope of this review, though I note that the careful model-data analysis of the type performed by Kowalski et al., 2017 is crucial as we explore the impact of these effect in flare loop models.

In the standard flare model there is typically not enough power carried by the highest energy electrons to meaningfully heat the deepest chromospheric/photospheric layers. However, there are some strands of evidence that suggest we do indeed require heating deeper than models currently predict. Listing some examples: excess line widths of chromospheric transitions (e.g., Mg ii h and k) cannot be accounted for; there is evidence of heating at the temperature minimum region (e.g., from inversions of Mg i λλ4571Å and λλ5173Å, Metcalf et al., 1990a; Metcalf et al., 1990b); there is speculation that white light flares (WLFs) may originate from the photosphere via enhanced H⁻ emission following a local temperature increase, or contain significant contributions from the lower atmosphere (see discussions in Neidig, 1989; Machado et al., 1989; Neidig et al., 1993; Martínez Oliveros et al., 2014; Kerr and Fletcher, 2014; Kleint et al., 2016; Jurc̆ák et al., 2018). An alternative explanation to the WLF problem is optically thin bound-free (recombination) radiation resulting from overionisation of the mid-upper chromosphere (e.g., Hudson, 1972). This would produce a Balmer jump at 3,646 Å. Of course, likely both of those mechanisms play a role. If we do need heating to great depth, then we must identify the agent capable of doing so, and constrain how much energy is required.

Given the scarcity of white light continuum observations, the NUV continuum as observed by IRIS is one such means to constrain the need for deep heating (the NUV is thought to be closely related to the optical continuum, albeit we do not yet know if they always originate from the same volume during flares). Heinzel and Kleint (2014) first determined that the Balmer continuum could be observed by IRIS, using observing windows in the NUV near λ = [2,813–2,816, 2,825–2,828, 2,831–2,834] Å. They carefully extracted narrow, line free, portions of the spectrum, finding ∼ 100 − 200 % contrast compared to the pre-flare values, with an impulsive rise and more gradual decay. Comparing to bright, non-flaring features, they note that the continuum rose but lines did not (as in the flare case) suggesting that the continuum patches between the lines are unaffected by the line emission. Using the static flare models of Ricchiazzi and Canfield (1983), processed with the radiation transfer code MALI, they inferred from the similarity in model-to-data intensities and from the formation properties in the model, that the observed NUV spectra did represent the Balmer continuum, and that it was due to optically thin recombination emission in the upper chromosphere. This represented the first detection of the Balmer continuum from space based instruments, and provides a constraint on flare models due its proximity to the Balmer jump.

Following on from this initial detection, Kleint et al., 2016 studied the IRIS Balmer continuum emission alongside other continuum enhancements from both space and ground observatories, spanning the UV through infrared. They performed blackbody fits to data (including modifying the blackbody intensity due to opacity effects) to determine the viability of an upper photospheric origin to the continuum emission, finding that the NUV lay well above the blackbody curve predicted by the optical and IR emission (as expected if the NUV emission was indeed recombination radiation, producing a Balmer jump). Building upon the modelling work started by Heinzel and Kleint (2014), Kleint et al., 2016 selected a few models from Ricchiazzi and Canfield (1983) and calculated the non-LTE hydrogen recombination continuum using the MALI code. Several were consistent with the observed NUV continuum enhancements (with a non-thermal electron flux close to that derived from RHESSI observations for that flare), but those models under-predicted the optical and IR enhancements from SDO/HMI and the Facility Infrared Spectropolarimeter at the Dunn Solar Telescope (DST/FIRS). Instead, a semi-empirical model atmosphere with photospheric temperature rise was required to achieve consistency with the optical/IR observations. Finally, they manually modified atmospheres that were input to RH in an attempt to find a stratification consistent with all three regimes (UV, optical and IR). A model with a modest photospheric temperature increase alongside a strongly heated chromosphere was required, as shown in Figure 7. Thus, IRIS in combination with ground based observations demonstrated that in some events we may indeed need both chromospheric and photospheric heating. To my knowledge no time-dependent flare model has self-consistently produced such an atmosphere (electrons beams typically do not carry enough power to such depths), and this should be a focus of our continuing efforts. However, it should also be noted that since an optically thin source at T ∼ 10 kK produces optical emission with a radiation temperature of 4–6 kK (see e.g. Kowalski and Allred, 2018), some ambiguity remains. At the same time as attempting to model self-consistently the heating throughout the chromosphere and photosphere to determine how to obtain atmospheres similar to the empirical models of Kleint et al., 2016, we should endeavour to obtain IRIS observations in the NUV alongside a broad spectral coverage of the optical (e.g. from DKIST) to determine the spectral shape more accurately, which would help resolve the ambiguity over emission mechanisms.

Since it will likely remain challenging to obtain observations covering the Balmer jump (and thus a guide as to the formation of the optical continuum), Kowalski et al., 2019 have begun to search for alternative metrics that can gauge the extent to which the lower atmosphere is strongly heated. Using the fact that Fe ii lines observed by IRIS form under similar physical conditions as the NUV Balmer continuum (T ∼ 8–18 kK but mostly towards the cooler end, discerned from their earlier modelling work regarding high beam fluxes; Kowalski et al., 2017), they explored the ratio of wavelength-integrated flare excess Fe ii 2,814.45 Å intensity to the average continuum intensity in the region λ = [2,824.5–2,825.9] Å. This ratio was observed to be R _Fe:NUV ∼ 7–8 at the peak of very bright flare sources located in a sunspot umbra during the 2014-Oct-25th X class flare, significantly higher than the prediction from slab models with low-to-moderate densities of ρ < 10^–9 g cm⁻³ which had values R _Fe:NUV ∼ 1. They speculate that this means there is significant heating (to T ∼ 10 kK) at high column depth (log  m ∼ − 2 [g cm⁻²]) where Fe ii can be optically thick. They are currently modelling this ratio in a range of RADYN flares (private communication 2022), but noted that their earlier study of the 2014-March-29th X class flare (Kowalski et al., 2017) only produced an observed ratio of R _Fe:NUV ∼ 1, with a modelled ratio of R _Fe:NUV ∼ 1–1.8. Clearly there is something quite different at work during the 2014-Oct-25th flare. This could be due to the pre-flare atmospheres, since the 2014-October-25th flare sources propagated into the sunspot umbra, allowing a colossal 1,000% NUV contrast, and an excess intensity 20× that of the 2014-March-29th flare. Hopefully further observations from a variety of flares, and modelling of a variety of energy inputs (including varying the pre-flare atmosphere) will lead to a firm diagnostic of deep heating during IRIS flares.

The NUV continuum was forward modelled by Heinzel et al., 2016, using one of the short-pulse experiments of Kašparová et al., 2009. In those FLARIX simulations, a non-thermal energy flux was injected in a trapezoidal form over time, with peak flux F = 4.5 × 10¹⁰ erg s⁻¹ cm⁻². The whole pulse was very short, only lasting 3 s. The resulting Balmer continuum intensity was quite small compared to the observations of Kleint et al., 2016. This was attributed to the magnitude of energy flux deposited in the chromosphere. It was an order of magnitude smaller than that of the flare studied by Kleint et al., 2016. Further, the short duration of this relatively moderate injection meant that evaporation was weak and the pressure in the upper chromosphere did not increase to the level inferred from the best-match Ricchiazzi and Canfield (1983) analysed by Kleint et al., 2016. This could point to the need for either longer electron beam dwell times in large flares, or for a train of short pulses. In those models the hydrogen subordinate continua, in particular the Balmer continuum, were the dominant source of radiative losses throughout the chromosphere, overtaking losses from singly ionised metals such as Ca ii and Mg ii, underscoring the importance of IRIS Balmer continuum observations. Since the Balmer continuum is seemingly optically thin, the radiative losses integrated through the continuum formation heights are directly related to the emergent intensity. Thus, the observed excess intensities impose strict constraints on flare energetics.

The Balmer continuum is also a useful constraint for smaller events where heating is, largely, confined to the upper chromosphere. A ‘mini-flare’ event accompanying a jet was studied by Joshi et al., 2020 who determined that a reconnection event occurred at the base of the jet. A small Balmer continuum excess was present in very localised sources, from both IRIS SG spectra, and the SJI 2832 Å images (Joshi et al., 2021). Comparing to the analysis of Kleint et al., 2016, who had processed the Ricchiazzi and Canfield (1983) atmospheric models using MALI to obtain predictions for the hydrogen recombination continuum, Joshi et al., 2021 found that a few simulations were consistent with their observations. This balanced the continuum intensity, as well as the brightness of the Mg ii line cores. The most well matched models had non-thermal electron energy fluxes F = 1 × 10^9–10 erg s⁻¹ cm⁻², with δ = 5 and E _c = 20 keV (note that the Ricchiazzi and Canfield 1983 did not sample other values of E _c and only a few values of δ for more energetic flares). From the Fermi/GBM (Meegan et al., 2009) hard X-ray observations, the injected non-thermal electron distribution was calculated as having an energy flux F = 6.5 × 10⁹ erg s⁻¹ cm⁻² for E > 20 keV, assuming the area into which the electrons were injected was the same as the continuum enhancement source. While there is some uncertainty in the low-energy cutoff, the parameters derived from the hard X-ray observations are consistent with those models that also provide a well-matched Balmer continuum excess intensity. Joshi et al., 2021 estimate that 82% of the intensity in the IRIS SJI 2832 Å images is continuum emission, contrasting the result from similar analysis of an X-class flare that found significant line emission (Kleint et al., 2017), which makes sense if we consider the weak energy flux involved that was unable to sufficiently excite the many lines within that part of the spectrum. In summary, small-scale reconnection at the base of a jet seemingly was able to accelerate enough electrons to bombard the upper chromosphere, enhancing the Balmer continuum, but was unable to really effect the lower chromosphere.

While there is unambiguous observational evidence for the presence of non-thermal electrons in a great many solar flares, they are not ubiquitously present. This is true both in the global sense, meaning there are some “thermal” flares that do not exhibit strong evidence of hard X-rays or microwaves of non-thermal origin, and the local sense meaning that hard X-rays sources do not appear uniformly along optical/UV flare ribbons (though this latter case may be related to the dynamic range of most hard X-ray observatories that precludes detection of weak sources alongside strong footpoints). In such flares in situ heating of the corona results in a conductive heat flux that transports energy to the lower atmosphere, generating the strong heating and mass flows (e.g. Zarro and Lemen, 1988; Battaglia et al., 2009; Fletcher et al., 2011; Brosius, 2012; Brosius and Holman, 2012; López et al., 2022). This is often referred to as ‘direct heating’, a somewhat nebulous term that refers generally to any heating of the corona (either looptops or along the legs of the loop) following the release of energy during magnetic reconnection ¹⁵ , including the retraction of magnetic loops that produce shocks such as those modelled by PREFT. It is indeed likely that some form of direct coronal heating acts alongside non-thermal electrons even in flares with clear evidence of particle acceleration, but the dominance of each mechanism varies from flare to flare and with spatial location. For example, recent results using GOES soft X-ray observations suggest that in a number of flares the corona is rapidly heated to T ∼ 10–15 MK before the onset of evaporation (so-called “hot onsets” Hudson et al., 2021). Several studies have either modelled flares as being purely driven by thermal conduction (e.g., Cheng et al., 1983; MacNeice, 1986; Gan et al., 1991, to name but a few) or have contrasted predictions between electron beam and conduction driven flares (e.g., Polito et al., 2018b; Kerr et al., 2021; Cheung et al., 2022). This latter exercise should be performed more often, as it is likely that both mechanisms act but the direct heating in the corona is often ignored. As shown in Cheung et al., 2022, the inclusion of direct heating can have impacts on the predicted intensities and Doppler motions of coronal and transition region lines. Some of those effects can only be seen at very high spatial and temporal resolution, such as will be afforded by the MUSE mission (see Section 4).

I summarise in detail here two recent examples of modelling flares driven purely by a conductive heat flux, and what we can learn about the nature of condensations from those simulations and IRIS observations.

The seminal studies of mass flows in flare models by Fisher et al. (1985a), Fisher et al. (1985b), and Fisher (1989) revealed the relationship between flare energy input and the development of both upflows and downflows in the chromosphere. Of those, Fisher (1989) concentrated on chromospheric condensations, developing an analytical model that described the timescales and magnitudes of flare-induced downflows. Fisher (1989) built that equation of motion from generalising various properties that occurred in flare loop models. Notably, Fisher, (1989) discovered that the lifetime of the condensations depend only on the chromospheric conditions, not the energy input, and is τ _life ≈ 2(H/g)^1/2, where H is the chromospheric density scale height and g is gravitational acceleration. For reasonable values of H, τ _life ∼ 60 s. The half-life of the condensation was τ _1/2 ≃ 2(H/g)^1/2/M _peak , where M _peak is the ratio of the peak downflow velocity to the sound speed in the pre-flare chromosphere. Though τ _life does not depend on properties of the energy injection, the peak downflow velocity (and thus τ _1/2) does, varying as u ₀ ∝ F ^1/3, with some dependence on whether energy was transported via non-thermal electrons or was conducted down from a hot corona.

Ashfield and Longcope (2021) used the PREFT gas-dynamic code to study conduction-only driven chromospheric condensations, building upon and complementing the work of Fisher (1989). They found similar relationships, with some differences as described below, but using an alternative approach. They set up a simplified set of physical parameters to explore the dynamics of shocks in the chromosphere using jump conditions, and compared those predictions to numerical experiments. The analytical model informed a subsequent analysis and interpretation of numerical results from PREFT. Unlike Fisher (1989), Ashfield and Longcope (2021) allowed H to be a free parameter in their model, from which they found a similar scaling for the half-life, but with a different factor (2.8 versus 2). Casting in terms of u ₀ this was: τ _1/2 ≈ 2.8H/u ₀. PREFT was set up as a rigid flux tube so that the only energy transport was via energy injected to the loop to mimic direct flare heating, and conduction-driven flares were simulated. Using those simulations with various values of pre-flare H and injected energy flux F, the velocity as a function of time u(t) was fit with a similar functional form as that from their analytical model. The best-fit values suggested that the relationship was instead τ _1/2 ≈ 1.67H/u ₀, with the discrepancy attributed to at least one of the simplifying assumptions in their analytical model not being satisfied. The results of Fisher (1989) and Ashfield and Longcope (2021) demonstrate that both analytically and numerically the pre-flare chromospheric density scale height H ∝ u ₀ τ _1/2, with the numerical results of Ashfield and Longcope (2021) pointing to H ≃ 0.6u ₀ τ _1/2. The slight difference between the results of Ashfield and Longcope (2021) and the earlier work of Fisher (1989) could be due in part to the heating profile assumed. Fisher (1989) assumed short pulses of non-thermal electrons in the models which the analytical expressions used a base, whereas Ashfield and Longcope (2021) assumed direct heating in the corona. Nevertheless, this is a rather powerful diagnostic that suggests variations in the lifetime of condensations could be related to variations in the pre-flare chromospheric densities into which shocks propagate, and that variations along a flare ribbon could reveal corrugation of the pre-flare atmosphere. Note that the condensations referred to here are the relatively strong, transient, downflows that may appear on top of a longer-lived envelope as discussed Paper 1. The peak velocity, u ₀, itself scaled with the input energy flux, without much reliance on H, going as u ₀ ∝ F ^1/2 for weaker energy fluxes (F < 2 × 10¹⁰ erg s⁻¹ cm⁻²), and u ₀ ∝ F ^1/3 for stronger fluxes. The latter was predicted by Fisher (1989), and the former by Longcope (2014) for low-energy fluxes.

Applying their findings to an actual flare, Ashfield et al., 2022 analysed the 2014-October-25th X-class event. During the flare there were persistent redshifts of Si iv, with v _Dopp ∼ 10 km s⁻¹, but on top of which were many transient ( < 1 minute) redshifts of several tens of km s⁻¹. After carefully analysing the IRIS Si iv spectra to extract a candidate condensation event to model, they used SDO/AIA and HMI data to trace a magnetic loop. Imaging the hard X-rays from RHESSI in the range 25–50 keV revealed only a coronal source, with no evidence of non-thermal particles at the footpoints. Thus, they used PREFT to model this event as a ‘thermal’ flare, to confirm that the observationally derived energy flux could drive the observed condensation. The density scale height was estimated from u ₀ and τ _1/2 as H = 369 km, which was used to initialise a rigid flux tube of length 85 Mm. Since no evidence of electron precipitation into that footpoint was present, Ashfield et al., 2022 used the UV Footpoint Calorimeter method (UFC; Qiu et al., 2012; Liu et al., 2013; Zhu et al., 2018; Qiu, 2021) to determine the time-dependent energy flux injected into the loop. In that method, the energy flux is proportional to the intensity in the SDO/AIA 1600 Å passband, with a scaling factor constrained by the SDO/AIA EUV channels and GOES Soft X-ray channels. The Enthalpy-Based Thermal Evolution of Loops (EBTEL; Klimchuk et al., 2008) model was used to synthesise the EUV and soft X-ray intensities for different energy fluxes, from which the best-fit match to the observations is used to define the scaling constant. EBTEL is a 0D model, and thus very quick to run, allowing many thousands of possible solutions to be generated efficiently to obtain the best-fit match. The 1,600 Å peak was fit by a Gaussian function, to obtain an energy flux profile with peak flux F = 6.2 × 10⁹ erg s⁻¹ cm⁻², and duration of a few minutes. The heating profile and observed lightcurves are shown in the lefthand column of Figure 8.

Injecting this derived energy flux profile into the PREFT loop resulted in the formation of a condensation, peaking prior to the peak of the energy input, and reaching a maximum of velocity 45 km s⁻¹, before decreasing to 10 km s⁻¹, again prior to the peak energy input. Only a small fraction of the input energy flux was required to drive the condensation, leading to the conclusion that the timescales of condensations do not necessarily impart knowledge of the duration of energy input, in agreement with the analytical work of Fisher (1989). From the dynamic PREFT simulation, the Si iv spectral lines were synthesised assuming optically thin emission, and summed through the extent of each leg of the loop. Non-equilibrium effects were considered by tracking the change in ionisation state for each Lagrangian grid cell. While qualitatively consistent with the observed Si iv emission from the footpoint, with similar redshifts produced (see the righthand column of Figure 8), the synthetic intensities did not track the condensation evolution and were over an order of magnitude too high. The synthetic profiles were also fully redshifted rather than showing two component behaviour, possibly pointing to the need for many strands along the lines of Reep et al., 2018a’s multi-threaded modelling. Since the synthetic profiles continued to increase in intensity after their peak redshift, it is possible that the Gaussian form assumed for energy input was not correct. Keeping the total energy from the UFC method, but having a more impulsive initial energy release could decrease the energy deposited later in the event, reducing the line’s intensity. Prior experiments with PREFT that included the loop retraction exhibited much more rapid energy release (Longcope and Klimchuk, 2015). Of particular importance here, aside from testing the work of Ashfield and Longcope (2021), and demonstrating that the observed conductions could be produced in a thermal conduction-driven flare, was that this was the first simulation of a flare in which the energy input to the chromosphere was inferred from coronal observations.