© ESO 2018

1 Introduction

Solar flares are spectacular events of enhanced emission in a wide range of the electromagnetic spectrum including hard X-rays (HXR; Kontar et al. 2011) and γ-rays (Vilmer et al. 2011), soft X-rays (SXR; Culhane et al. 1992; Sylwester & Sylwester 1999; Khan et al. 2004; Tomczak & Ciborski 2007) and ultraviolet (UV) radiation (Fletcher et al. 2011; Li et al. 2015; Battaglia et al. 2015; Brosius et al. 2016), optical (Fletcher et al. 2011) and radio (Smith et al. 2011) emission, as well as a seismic response to flaring events from the solar interior (sunquakes) (Kosovichev & Zharkova 1998; Zharkova 2008; Donea 2011). While high energy HXR, SXR, and UV emission have been intensively investigated in recent decades, observations of optical emission were recently intensified with the launch of the Interface Region Imaging Spectrograph (IRIS; De Pontieu et al. 2014). The discoveries made with IRIS have renewed motivation for ground-based observations of hydrogen emission of flares.

Observations of hydrogen spectral lines provide a diagnostic of the conditions and dynamics of the lower solar chromosphere at various heights. During the observations of impulsive phases of solar flares, taken using narrow wavelength windows, there are frequently reported delays of the order of 30 s in the maximum of the Balmer alpha line (Hα) intensity, when compared to the HXR signal (Kaempfer & Magun 1983; Veronig et al. 2002; Radziszewski et al. 2011). In observations with larger spectral windows, the Hα line often shows large red-shifts around the time of peaks in HXR spectra (Ichimoto & Kurokawa 1984; Wuelser & Marti 1989; Canfield et al. 1990; Wulser et al. 1992; Wuelser et al. 1994; Druett et al. 2017). Such sizeable Doppler shifts indicate large velocities of the local plasma motion towards the photosphere. These line profiles can be used to analyse the proxies of energy deposition in a flaring atmosphere.

Wuelser & Marti (1989) attributed the red-shifts to energy release in the higher atmosphere and a complex structure of the flaring region. Kaempfer & Magun (1983) suggested that delays between the HXR peaks and Hα maximum can be used to establish different transport mechanisms of energy. Modern instruments such as the Big Bear Solar Observatory(BBSO) and Swedish 1 m Solar Telescope (SST) are capable of high spatial and temporal resolution imaging that can be used to carry out observations in the early stages of flares and to distinguish the sources of observed Hα profiles withstrong red-shifts. In this paper the observations and interpretations made by Wuelser & Marti (1989) and Kaempfer & Magun (1983) are re-evaluated in light of new observational and modelling advances.

Early observations of white light (WL) flares by Uchida & Hudson (1972; see also Hiei 1987; Kurokawa et al. 1988) show that WL sources often occur simultaneously with HXR bursts, gamma rays, and optical emission in the hydrogen lines. Balmer continuum enhancement has been observed to be co-temporal and co-spatial with HXR sources (Kleint et al. 2016; Kotrč et al. 2016). Hudson (1972) suggested that an energetic electron beam, accelerated in the corona and injected into the chromosphere, could ionise the flaring plasma causing the temporal correlation of HXR and WL. However, observations of co-temporal HXR and WL on the solar limb (Martínez Oliveros et al. 2012) indicate that both can have sources located close to the height of the quiet sun photosphere, below the stopping depths for beam electrons derived from the quiet sun atmospheres of Fontenla et al. (2009).

Authors have also sought alternative mechanisms that could increase the lower chromospheric and photospheric ionisation degree co-temporally with HXR sources, suggesting energy transport from the corona by Alfven waves (Fletcher & Hudson 2008; Kerr et al. 2016). The IRIS and the Image Selector offer greater scope for the analysis of high-cadence changes in continuous emission in WL and near-UV, as reported in Heinzel & Kleint (2014) and Kotrč et al. (2016). Their observations have provided an ideal opportunity to examine the predictions of the electron beam delivery model using updated simulations with a dynamic hydrodynamic response to beam heating combined with fully non-local thermodynamic equilibrium (NLTE) radiative modelling (Varady et al. 2010).

The first radiative models were developed in the 1980s (Ricchiazzi & Canfield 1983; Aboudarham & Henoux 1986, 1987) for the interpretation of hydrogen emission, suggesting that hydrogen ionisation and excitation in a semi-empirical flaring atmosphere is caused by non-thermal electrons precipitating from the corona to chromosphere. These models used electron beam kinetics for a flux conservation approach (FCA) developed by Brown (1971); Emslie (1978), which was later shown (Mauas & Gómez 1997) to overestimate the plasma heating of upper atmosphere and to lead to infinite singularity at the stopping depths of electrons with lower energies. To avoid the latter, beam electrons in FCA were restricted from reaching the middle and lower chromosphere (Kontar et al. 2011). Using this model and approximate hydrogen excitation and ionisation rates derived from the beam heating functions, Ricchiazzi & Canfield (1983) and later Aboudarham & Henoux (1986, 1987) showed that the ionisation rates for collisions between beam electrons and ambient hydrogen atoms are orders of magnitude higher than thermal collisional ionisation rates in the middle and low chromosphere. These authors indicated that the ionisation from the non-thermal beam would significantly modify the intensities of the hydrogen continua and chromospheric line wings. Moreover, Aboudarham & Henoux (1986, 1987) suggested that this ionisation could increase the H⁻ population and enhance white light emission at photospheric and chromospheric depths, producing a 5% contrast at disk centre at a wavelength of 600 nm.

Using the same beam kinetics as above, Canfield & Gayley (1987) computed time-dependent Hα profiles using the hydrodynamic model atmospheres of Fisher et al. (1985) for an impulsive injection of power-law beam electrons. In their models plasma macro velocities generated Hα line profiles that show Doppler shifted contributions in the wings. Line centres were shifted to the red wing during the impulsive phase, but the relative intensities always peaked in the blue wing, in contrast to the observations of Ichimoto & Kurokawa (1984); Wuelser & Marti (1989); Canfield et al. (1990); Wulser et al. (1992). Heinzel et al. (1994) simulated Hα line profiles with slight blue asymmetry, which they attributed to the downflow of chromospheric plasma absorbing emission in the red wing. In particular these authors compared their simulations with profiles from a flare on 4 October, 1991. The timing of the observed profiles coincided with a peak in HXR at 09:37 and had been preceded by a HXR burst at 09:35. However, none of the synthetic line profiles could reproduce the observed Hα line profiles with large red-shifts.

Alternatively, Zharkova & Kobylinskii (1991, 1993) presented a NLTE model of hydrogen emission caused by precipitation of beam electrons with kinetics following the continuity equation approach (CEA; Syrovatskii & Shmeleva 1972; Dobranskis & Zharkova 2015). This approach revealed that beam electrons can easily reach lower chromosphere and photosphere providing a smooth atmospheric heating throughout the whole flaring atmosphere. The analytical CEA solutions for electron beam precipitation and their ability to reach photospheric levels were later validated by numeric simulations with the Fokker-Planck-Landau approach (see, for example Zharkova et al. 1995; Zharkova & Gordovskyy 2005, 2006; Siversky & Zharkova 2009, and references therein). Moreover, Zharkova & Kobylinskii (1989) calculated analytically the accurate non-thermal hydrogen excitation and ionisation rates by beam electrons based on electron power law energy distributions and semi-empirical cross sections used for thermal electrons on the Sun.

Using the accurate non-thermal hydrogen excitation and ionisation rates and considering hydrogen ionisation to be governed by a radiative transfer in Lyman continuum (see for example Morozhenko & Zharkova 1982; Zharkova 1989, and references therein), Zharkova & Kobylinskii (1993) applied the NLTE simulations to the hydrodynamic atmosphere heated by electron beams following CEA (Somov et al. 1981) and demonstrated a dominance of non-thermal processes for hydrogen excitation and ionisation below the transition region. These authors directly compared the hydrogen ionisation degrees for thermal and non-thermalelectrons and for the first time presented the evidence that non-thermal beam electrons significantly increase the ionisation degree at the lower chromosphere and photosphere. The additional free electrons generated by non-thermal ionisation were suggested as a plausible source of WL flares (Aboudarham & Henoux 1986; Zharkova & Kobylinskii 1993; Zharkova 2008). However, the authors did not simulate line profiles with strong red-shifts produced in the dynamic atmospheres with strong downflow motions.

In the modern era, Allred et al. (2005) produced the new radiative hydrodynamic code RADYN, which incorporates the multi-element NLTE radiative code MULTI and accurately synthesises spectral lines of hydrogen, helium, and metals using the method of linearised non-local statistical equilibrium equations for population numbers. The hydrogen atom was treated in RADYN by solving continuity equations combined with steady-state equations for continuous transitions in each levels from 1 to 6 without considering radiative transfer in continua and neglecting non-thermal ionisation from the excited (>1) states of hydrogen atoms. The initial flaring atmosphere was assumed to be formed from the semi-empirical flaring atmosphere VAL F (Vernazza et al. 1981), to which the quiet sun corona was attached. This initial flaring atmosphere was heated by beam electrons considered in FCA discussed above with the heating function truncated at the upper chromosphere to avoid the infinite discontinuity (Mauas & Gómez 1997) and in the corona to avoid overheating. Then hydrodynamic equations of energy, flux, and momentum were solved in linear coordinates to ease the radiative transfer treatment. This approach has been successful in interpreting metal lines in flares caused by weaker beams with lower initial fluxes about 10¹⁰ erg cm⁻² s⁻¹ (Kašparová & Heinzel 2002; Allred et al. 2005).

Using RADYN Kuridze et al. (2015) have also been successful in describing the Hα line profile during a thermal heating phase of a flare (without red-shifts), while Rubio da Costa et al. (2016) have proposed the idea of multiple “threads”, or multiple beam injection sites, for possible interpretation of the observed Hα spectral line profiles with large red-shifts. Recently Druett et al. (2017) have reported a very close agreement between the simulated Hα line profiles for the F10 model and the Hα line profiles with large red-shifts observed in a C1.5 flare on June 30, 2013 with the CRisp Imaging Spectro-Polarimeter (CRISP) at the SST (Scharmer et al. 2003). The flaring atmosphere was assumed to be generated by a hydrodynamic response to heating by beam electrons in CEA and these authors used analytical non-thermal excitation and ionisation rates (Zharkova & Kobylinskii 1993).

At the same time Kowalski et al. (2017) have successfully simulated the timing and scale of the red-shifts observed in Fe II chromospheric lines at the time of a secondary HXR peak (their “BFP2”), away from the brightest kernel at the time, which was likely to mark the main source of the injection. Simulations at another footpoint “BFP1” were shown from 0–8 s, i.e. during the first beam injection, although the observations were taken 28 s after the initial HXR impulse (see their Fig. 2), by which time their simulations would have shown negligible red-shifts (see their Fig. 9). As yet, their models have not managed to account for Hα line profiles with large Doppler red-shifts and highly enhanced wing emission observed during an impulsive phase of solar flares (Ichimoto & Kurokawa 1984; Wuelser & Marti 1989; Canfield et al. 1990; Wulser et al. 1992; Wuelser et al. 1994), or co-temporal HXR and WL sources located close to the height of the quiet sun photosphere (Martínez Oliveros et al. 2012; Krucker et al. 2015).

Inspired by our successful interpretation (Druett et al. 2017) of Hα line profiles from a C-class solar flare observed with SST (Sweden), in this paper we proceed with the attempt to explain the other observations with larger red-shifts reported above. For this purpose we follow the model described by Druett et al. (2017) by employing the heating function by electron beams developed from CEA (Syrovatskii & Shmeleva 1972). Electron beams are assumed to be injected into a quiet Sun chromosphere and their precipitation causes strong hydrodynamic responses of the ambient plasma from the corona to the photosphere (Somov et al. 1981; Zharkova & Zharkov 2007). This modelling technique allows us to investigate hydrodynamic models heated by stronger and harder proton and mixed beams (Zharkova & Zharkov 2015) and to consider greater beam electron densities at given atmospheric depths, resulting in greater non-thermal excitation and ionisation rates compared to the approach used for FCA (Ricchiazzi & Canfield 1983; Aboudarham & Henoux 1987).

In this paper we present the updated NLTE radiative models for flaring atmospheres, which are hydrodynamic responses to electron beam injections. We focus on the beam injection phase, during which the quiet Sun atmosphere is converted into a flaring atmosphere. We describe accurately radiative mechanisms and emission that characterise the responses in continua and lines of the Balmer and Paschen series. The radiative and hydrodynamic models and solution methods are described in Sect. 2, simulated line and continuum profiles are discussed in Sect. 3, and their comparison with observations are discussed in Sect. 4. Conclusions are drawn in Sect. 5.

2 Hydrogen in a hydrodynamic atmosphere

2.1 Physical conditions

For the physical conditions of flaring atmospheres we use models of hydrodynamic responses of ambient plasma to energy deposition by energetic electron beams precipitating from the corona to the lower atmosphere (Zharkova & Zharkov 2007). These models have the initial conditions of a quiet Sun chromosphere from a column depth just below the transition region of the quiet Sun (ξ = 10¹⁷ cm⁻²) to the beginning of the upper photosphere (ξ = 10²² cm⁻²). Details of these initial conditions are given in Somov et al. (1981) and include (a) a constant temperature of 6,700K derived from semi-empirical calculations; (b) hydrostatic equilibrium v(0, ξ) = 0; and (c) a density distribution as defined by the straight line in the logarithmic plot of Fig. 1b. The beams are assumed to be power laws in energy with initial fluxes of F₀, a lower energy cut-off of E_low = 10 keV, and a spectral index γ = 3 for the electron spectrum. During the beam injection of 10 s, their initial fluxes are modelled to linearly increase for the first 5 s and then linearly decrease until the beam is switched off at 10 s. The simulations are carried out for fluxes of 10¹⁰, 10¹¹, and 10¹² erg cm⁻² s⁻¹, hereafter referred to using the shorthand names of F10, F11, and F12, respectively. The models provide the physical conditions of flaring atmospheres including hydrogen number density, n(ξ), electron kinetic temperature, T(ξ), and plasma macro velocities, v_m(ξ), as functions of the column depth, ξ.

The column depth, ξ, is defined in terms of a total number of particles per cm squared above the depth in question in the flux tube, ξ = ∫ n₀ds. Here n₀ (s) is the total number density of the plasma and s is a linear depth; the integration is carried out from a given point to the top of a flaring atmosphere.

The simulated hydrodynamic responses are shown in Fig. 1, where blue lines represent the F12 flare model and green and red represent the F11 and F10 models, respectively. The rows display the response from 1 to 5 s after the onset of the beam injection from top to bottom, respectively. The left panels show the electron kinetic temperature response (in Kelvin), the central panels represent the number density of the plasma (cm⁻³), and the right panels represent the macro velocity responses of the plasma (km s⁻¹), each plotted against column depth, ξ on the x-axis (cm⁻²).

It can be seen from Fig. 1 that within seconds an electron beam converts the quiet Sun chromosphere into a flaring atmosphere with its own flaring corona, transition region, and chromosphere. After five seconds of beam injection these regions can be discerned most clearly from the temperature responses with the coronal temperatures (T >10⁶ K) dropping sharply in the flaring transition regions to the temperatures of flaring chromospheres (T ~ O(10⁴) K). Figure 1 (left panels) shows the temperature responses against column depth from times from 1 to 5 s. Magnitudes of the temperature increases are positively correlated with the initial flux of a beam and anti-correlated with its spectral index.

The heating caused by beam electrons sweeps the coronal plasma to deeper atmospheric levels while creating a hydrodynamic shock travelling towards the photosphere and beneath (Zharkova & Zharkov 2007, 2015). The column depth of the flare transition region is shifted towards deeper atmospheric levels with the increase of the initial flux of the beam from 3 × 10¹⁸ cm⁻² for the F10 model to 1 × 10²⁰ cm⁻² for the F12 model. The ambient plasma density also changes following a hydrodynamic response to the beam heating (Fig. 1, central panel column) showing, first, a decrease of chromospheric densities to the coronal densities due to the formation of a flaring corona above the flare transition region and, second, a high density and low temperaturecondensation below the flare transition layer due to the formation of shock moving downward to the photosphere.

The plasma heated by beam electrons is swept to deeper atmospheric levels causing upflows of chromospheric plasma into the corona and a formation of hydrodynamic shock of condensed plasma moving downward to the photosphere, which are seen in the macro velocity responses (Fig. 1, right panels). Since the quiet Sun plasma densities are lower at smaller column depths, the beams of smaller incident flux sweep less plasma, leading to a lower plasma density below the transition region of the flare than for more powerful beams (see Fig. 1, central column, t = 4 and 5 s). At 5 s we see that the high density peaks that result from downward travelling condensations have plasma densities of about 10¹³ cm⁻³ in each model (Fig. 1n). It is evident that for the F10 model the density increase in the condensation appears more noticeable compared to the surrounding plasma than in the F11 and F12 models.

Similar to the temperature increase, the velocities of generated hydrodynamic shocks are also scaled rather closely by the initial flux of an electron beam. For example, at 5 s after a beam onset the maximum downflow speed varies from 57 km s⁻¹ for the F10 model to 382 km s⁻¹ for the F12 model. The hydrodynamic shocks move downward to the photosphere and penetrate to larger column depths in the photosphere and beneath (Zharkova & Zharkov 2015). At the same time there are also strong upflow motions reflecting either gentle (for weaker beam) or explosive (for stronger beams) evaporation of the swept chromospheric plasma back into the corona with the velocities proportional to the initial flux of the beam (Somov et al. 1981; Fisher et al. 1985; Zharkova & Zharkov 2007). This evaporation starts 100–300 s after the beam onset and lasts for a really long time of few thousand seconds or longer until the initial quiet sun atmosphere is restored (Somov et al. 1981; Fisher et al. 1985).

Fig. 1 Hydrodynamic responses of the ambient plasma to electron beam heating for the F10 (red), F11 (green), and F12 (blue) models for times t = 1–5 s. These are shown from the top row to the bottom row, respectively, as follows: the left panels indicate electron kinetic temperature; central panels indicate plasma density; right panels indicate plasma macro velocity plotted vs. plasma column depth (cm−2 ).

2.2 Non-local thermodynamic equilibrium model for hydrogen

2.2.1 Statistical equilibrium

We considered a 5-level plus continuum hydrogen model. It is assumed that statistical equilibrium is established for each level of the hydrogen atom modelled at each instant of the hydrodynamic responses described in Sect. 2.1, similar to Zharkova & Kobylinskii (1993, see Appendix A.1). This assumption is based on the fact that the timescale for each radiative process is much shorter than the hydrodynamic characteristic time of tens of seconds (Shmeleva & Syrovatskii 1973; Somov et al. 1981; McClymont & Canfield 1983). The exception to this is the recombination rates, which are not so important for the period when beam electrons are injected because their ionisation rates are two orders of magnitude higher than the recombination rates. We did not consider the period after the beam injection stops here. However, after a beam stops, the slow recombination rate and optically thick Lyman continuum radiative transfer keeps the hydrogen excitation and ionisation level at high levels (Zharkova & Kobylinskii 1993). In Druett et al. (2017) we considered a slow recombination which sustains the high hydrogen ionisation for tens of seconds.

Statistical equilibrium of the hydrogen atomic level populations and ionisation during the electron beam injection phase is formed by considering electron transitions due to the following processes: excitation and de-excitation rates for collisions between hydrogen atoms and thermal electrons (Johnson 1972). Three-body recombination, in which electrons transfer energy and momentum to another electron in the vicinity of a hydrogen ion, result in recombination to a bound level (Johnson 1972). Excitation and ionisation rates due to collisions between hydrogen atoms and beam electrons are calculated with the rate coefficients found for semi-empirical cross sections using analytical formulae derived by Zharkova & Kobylinskii (1993). The depth distribution of the beam electron number densities used with these rates were modelled with the work of Zharkova & Kobylinskii (1993). We also considered hydrogen atomic transitions due to spontaneous radiative rates (Allen 1977) including recombination rates. Stimulated excitation or de-excitation and ionisation are included from external radiation (sometimes called backwarming) entering from the levels above and below the flux tube model. The rate coefficients for these processes were taken from Zharkova & Kobylinskii (1991, see Appendix A.1). We also accounted for transitions and ionisation due to internal diffusive radiation. These parameters were incorporated by solving radiative transfer equations and converting the steady-state algebraic equations into Fredholm integral radiative transfer equations of the second kind. The Fredholm integral equations were solved via the L2 approximation method described in Sect. 2.2.3.

Zharkova & Kobylinskii (1993) compared the overall collisional rates in a simulated flare, s⁻¹, between the model with only thermal collisional processes and those including non-thermal electron beams using a variety of single power-law energy spectrum parameters. These comparisons were made in a differential study, that is in which the temperature and density profiles were taken from a single hydrodynamic response model, rather than varied to match the response to the particular beam used in the radiative code. Figures 3a to d in Zharkova & Kobylinskii (1993) compared the collisional rates for hydrogen in the Lyman alpha transition, ionisation from ground state, the Hα transition, and ionisation from the second state, respectively.

We find the following relationships based on these figures: (a) Non-thermal collisional processes dominate thermal collisional processes below the flare transition layer. This is particularly evident for ionisation from the ground state of hydrogen in which the collisional rates from non-thermal processes are more than 5 orders of magnitude greater than those for thermal processes for all the models used. (b) The magnitude of the dominance noted in point (a) is increased for a beam with a higher initial flux if both have the same spectral index. This is due to the greater beam electron density at each depth. (c) Beams of a lower spectral index (harder beams) deposit their energy deeper in the atmosphere than softer beams, as they contain more high energy particles that penetrate to greater column depths. This means that the ionisation and excitation rates are higher in the deep atmosphere for harder beams and higher in the upper atmosphere for softer beams.

The equations of statistical equilibrium, (A.1) and (A.2), are used alongside the particle conservation number Eq. (A.11).

2.2.2 Radiative transfer

We considered radiative transfer for all optically thick lines and continua (see Appendix A.3) and expressed in integral form. Assuming complete redistribution of frequencies, the mean intensities give the following expressions, as stated in Eqs. (11) and (12) of Zharkova & Kobylinskii (1991), for the simulated rates of excitation, de-excitation, and ionisation due to internal diffusive radiation: for excitation and de-excitation in the lines, $$B_{ik}{j}_{ik}^{\partial }(\tau )=\frac{A_{ki}}{2}{\displaystyle \underset{0}{\overset{{\tau }_{ik}^0}{\int }}\frac{n_k}{n_i}}(t)K_1(|\tau -t|)\text{d}t,$$(1)

and for ionisation in the optically thick Lyman continuum, $$B_{1c}{j}_{1c}^{\partial }(\tau )=\frac{A_{c1}}{2}{\displaystyle \underset{0}{\overset{{\tau }_{1c}^0}{\int }}\frac{n_e{n}_{+}}{n_1}}(t)K_1(|\tau -t|)\text{d}t,$$(2)

with first kernel functions in the lines, $$K_1(|\tau |)=A{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\alpha }^2}(\nu )E_1(\alpha (\nu )|\tau |)\text{d}\nu ,$$(3)

and in the continuum $$K_1(|\tau |)=F(T){\displaystyle \underset{{\nu }_{ic}}{\overset{\infty }{\int }}{f}_i}{\nu }^2\exp \left(-\frac{h(\nu -{\nu }_{ic})}{k_{B}T(\tau )}\right)E_1(f_i|\tau |)\text{d}\nu .$$(4)

The absolute value of the input to the kernel function K₁(|τ − t|) is used in Eqs. (1) and (2) so that contributions from upward and downward travelling radiation are both included in one integral. These formulae are substituted into the steady-state equations (Eq. (A.1)) via the diffusive rates of excitation D_ij, de-excitation D_ji, and ionisation D_ic, for all the diffusive terms relating to radiation that is not optically thin, i.e. non-zero D values. This substitution leads to the integral form of the radiative transfer equations as below $$S(\tau )=\frac{P_{\text{s}}(\tau )}{2}{\displaystyle \underset{0}{\overset{{\tau }_0}{\int }}{K}_1}(|\tau -t|)S(t)\text{d}t+S^{*}(\tau ).$$(5)

The variable S(τ) is a source function, including diffusive radiation. The primary source function is S^* (τ), which defines all other sources besides the diffusive sources. The limits of integration are taken from the top and bottom of the model, where optical depths for the line centre or continuum head are zero and τ₀, respectively. The value P_s (ξ) is the survival probability of a scattered photon at optical depth τ, i.e. the chance that an absorbed photon of this wavelength is re-emitted with the same frequency in the next atomic transition, $$P_{\text{s}}(\tau )=\frac{A_{ci}(\tau )}{{\displaystyle {\sum }_{k=1}^5{R}_{ck}(\tau )}}\cdot$$(6)

Equation (5) is a Fredholm integral equation of the second kind, α(ν) is the line absorption coefficient, and f_i(ν) are the relevant continuum absorption coefficients described above. The value F(T) is the normalisation function in the continuum and A is normalisation coefficient in lines defined as $$A{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}\alpha }(\nu )\text{d}\nu =1.$$(7)

In the continuum the source functions are related to the relative emission measure $\frac{n_e{n}_{+}}{n_i}(\xi )$ (REM) via the relationship $$\begin{array}{ll}\frac{n_e{n}_{+}}{n_i}(\xi )\hfill & =\frac{P_{\text{s}}(\tau )}{2}{\displaystyle \underset{0}{\overset{{\tau }_0}{\int }}\frac{n_e{n}_{+}}{n_i}}(t)K_1(|\tau -t|)\text{d}t\hfill \\ \hfill & +\frac{P_{\text{s}}(\xi )}{A_{ci}(\xi )}{\displaystyle \sum _{k=1,k\ne i}^5\frac{n_k}{n_i}(\xi )R_{kc}^{*}(\xi )}\hfill \end{array}$$(8)

where $R_{kc}^{*}$ are the rate coefficients for ionisation from level k (see Eq. (A.8)) without considering internal diffusive radiation.

2.2.3 Solution method

The system of 10 equations from Sect. 2.2.2, including Fredholm integrals (for optically thick transitions) and algebraic steady-state equations (for optically thin transitions), are solved using the iterative method first proposed by Zharkova & Kobylinskii (1993). For all the integral radiative transfer equations, the solutions are found in terms of the source functions S using the approximate method of L2 solutions developed by Ivanov & Serbin (1984) as follows: $$\begin{array}{ll}S(\tau )\hfill & \simeq S^{*}(\tau )\times {(1-P_{\text{s}}(\tau )+P_{\text{s}}(\tau )K_2(\tau ))}^{-1/2}\hfill \\ \hfill & \times {(1-P_{\text{s}}(\tau )+P_{\text{s}}(\tau )K_2({\tau }_{\text{max}}-\tau ))}^{-1/2}.\hfill \end{array}$$(9)

The use of L2 approximation significantly improves the accuracy of the solutions (Ivanov & Serbin 1984) in comparison with L1 solutions used by Zharkova & Kobylinskii (1993). The L2 approximation is related to the second order escape probability methods described by Rybicki (1984), but does not suffer from the same shortcomings of methods described therein. This is because it is based on a more rigorous analysis of the asymptotic behaviour of the solutions of the equations of radiative transfer, and is accurate with maximum error 20%, which is only present in the few top layers of the model. The second kernel functions employed in this technique are related to the first kernel functions by the relationship $$K_2(\tau )={\displaystyle \underset{\tau }{\overset{\infty }{\int }}{K}_1}(t)\text{d}t,$$(10)

and are found numerically via Chebyshev-Laguerre quadrature. The ratios of atoms with electrons in levels i and j, $\frac{n_j}{n_i}(\xi )$ are initially set via the Boltzmann distribution. After each iteration of the solution method the new ratios are calculated and used in the next iteration. The source functions are related to the relative populations of atoms with electrons in each level $\frac{n_j}{n_i}(\xi )$ in the transition lines via Eq. (A.14), and the relative emission measure $\frac{n_e{n}_{+}}{n_i}(\xi )$ (REM) in the continuum (Eq. (8)). The solution method is repeated until the source functions at each of the depth points become constant with the relative accuracy of 10⁻⁵.

2.3 Intensities of emission in the lines and continua

2.3.1 Line emission

The converged source functions are then used to calculate the overall intensities of emission, integrating over contributions from all depth points and angles. The formula for transition lines is $$\begin{array}{ll}{I}_{ij}(\lambda )\hfill & =2\pi \frac{c}{{\lambda }^2}\frac{2hc}{{\lambda }^3}\hfill \\ \hfill & \times {\displaystyle \underset{0}{\overset{{\tau }_{\text{max}}}{\int }}{\left(\frac{n_i{g}_j}{n_j{g}_i}(\tau )-1\right)}^{-1}}{e}^{\left(-\frac{{\alpha }_{ij}(\lambda )}{\mu }\tau \right)}\frac{{\alpha }_{ij}(\lambda )}{\mu }\text{d}\tau .\hfill \end{array}$$(11)

In this equation, the coefficient 2π accounts for axial symmetry of the generated emission from different azimuthal angles, and the relationship $\text{d}\nu =-\frac{c}{{\lambda }^2}\text{d}\lambda$ generates a factor of $\frac{c}{{\lambda }^2}$ to convert between the intensity formulae given per unit frequency (in Eqs. (A.18) and (A.19)) and per unit wavelength used. Additionally the formulas for the source functions in the lines (Eq. (A.14)) were converted into the functions of a wavelength, where g_i is the degeneracy of the quantum state i.

2.3.2 Continuum emission

The formula for the intensity of continuum emission escaping from the top of the model is derived in Appendix B and stated below in wavelength units (Eq. (12)) with the explicit τ dependence removed for compactness. The relationship $\text{d}\nu =-\frac{c}{{\lambda }^2}\text{d}\lambda$ is used to convert the intensity formulae for use with the results in wavelength units, and thus we have $$\begin{array}{l}\text{for}\hspace{0.17em}i=1,2,3\\ I_{ic}(\lambda )=\frac{h}{2}\frac{c}{{\lambda }^2}\\ \times {\displaystyle \underset{0}{\overset{{\tau }_{\text{max}}}{\int }}\frac{A_{ci}}{{\kappa }_{ic}}}\frac{n_e{n}_{+}}{n_i}\frac{{\lambda }_{ic}^3}{{\lambda }^3}\frac{\exp \left(-\frac{hc}{k_{B}T\lambda }\right)\frac{f_i}{\mu }\exp \left(-\frac{f_i}{\mu }{\tau }^{\prime }\right)}{E_1\left(\frac{hc}{k_{B}T{\lambda }_{ic}}\right)}\text{d}{\tau }^{\prime }\\ \text{for}\hspace{0.17em}i=4,5\end{array}$$ $$\begin{array}{l}\\ I_{ic}(\lambda )=\frac{h}{2}\frac{c}{{\lambda }^2}\\ \times {\displaystyle \underset{0}{\overset{{\tau }_{\text{max}}}{\int }}\frac{A_{ci}}{{\kappa }_{ic}}}\frac{n_e{n}_{+}}{n_i}\frac{{\lambda }_{ic}^2}{{\lambda }^2}\frac{\exp \left(-\frac{hc\left(\frac{1}{\lambda }-\frac{1}{{\lambda }_{ic}}\right)}{k_{B}T}\right)\frac{f_i}{\mu }\exp \left(-\frac{f_i}{\mu }{\tau }^{\prime }\right)}{\frac{k_{B}T\lambda }{hc}}\text{d}{\tau }^{\prime }.\end{array}$$(12)

This formulation uses radiative transfer in the continuum, governed by the optically thick Lyman continuum, to define the ambient ionisation and is affected by collisions, external, and internal diffusive radiation. This approach is different from the RADYN code using combined continuity and steady-state equations for six levels to define the ambient ionisation while neglecting the non-thermal ionisation by beam electrons from the excited states (2–6) of hydrogen (Allred et al. 2005). Moreover our code uses ionisation and excitation rates for inelastic collisions by beam electrons following the analytical formulae derived by Zharkova & Kobylinskii (1989, 1993). These rates were carefully compared with those by thermal electrons in Zharkova & Kobylinskii (1993), for one-to one and volume rates. Our code also differs from RADYN as a result of the number of beam electrons at a given depth and the depth where these electrons can reach (see Sect. 5).

2.3.3 Contribution functions

We investigated the emission of hydrogen in Balmer and Paschen series formed in model atmospheres with 150 points of column depth, converted into an optical depth. The emission from each layer contributes to the overall intensity at each wavelength, and this emission from the layer is a contribution function, which gives the contribution from a specific depth point to the total emission measured from the top of the flaring atmosphere at a given wavelength. The contribution functions are functions of the depth point number and wavelength.

Thus a contribution function for the depth point m is derivedby taking the relevant intensity integral (Eqs. (11) and (12)) with the integration limits between the optical depths of the layer with the depth point m, i.e. τ(m− 1) and τ(m). These contribution functions are utilised to investigate the origin of emergent radiation observed from the top of a flaring model.

3 Results of simulations

3.1 Analysis of optical depths in hydrogen lines

The optical depths of hydrogen lines are investigated using two approaches. In the first instance, using a differential study, the radiative responses to thermal collisions and to collisions with electron beams with the initial energy fluxes F₀, ranging 10⁹ −10¹² erg cm⁻² s⁻¹, are calculated for the same flaring atmosphere with the physical conditions defined by the hydrodynamic response to a beam with F₀ = 10¹⁰ erg cm⁻² s⁻¹. In the second instance, we use the matched hydrodynamic and radiative responses to complete the investigation.

There are several effects that can be derived from the expressions for optical depth (see Appendix C): (1) A higher ionisation degree implies a lower optical depth, all other things being equal, as this results in fewer hydrogen atoms being present that can absorb photons of the given wavelength. (2) The greater the proportion of neutral hydrogen atoms with electrons in the lower excited state of the transition, n_i, the higher the optical depth, as there is more material present that can destroy photons of the relevant wavelength by stimulated excitation or ionisation. (3) The only difference in the formulae for optical depths of lines in the same series is the absorption coefficient in the central wavelength. Therefore, aside from a multiplicative constant, the profiles of optical depths in the central wavelength (plotted against column depth) are the same for the lines and continuum of a series.

In regards to the first part of this investigation, Fig. 2 shows the optical depths calculated in the differential study for (a) the Lyman alpha line, (b) the Hα line, (c) the Pβ line, and 5 s after the onset of the beam. The optical depths are plotted against column depth, ξ (cm⁻²). In each subfigure the optical depth in the central wavelength of the line profile is shown for the thermal response (black line) and for models including excitation and ionisation rates from electron beams with fluxes in the range 10⁹ −10¹² erg cm⁻² s⁻¹.

The major radiative processes occurring directly as a result of inelastic collisions with beam electrons, i.e. non-thermal collisional excitation and ionisation, both cause a decrease in optical depth in the Lyman series. Both act to reduce the number density of hydrogen atoms with electrons in the ground state, which could otherwise absorb Lyman radiation. Thus, the stronger the beam the lower the optical depth for all Lyman lines. This is shown in Fig. 2a, where for all column depths modelled the optical depth that results from a model using a beam with a higher initial flux is less than or equal to the optical depth from a model using a lower initial flux. However, since the Lyman line cores become very optically thick just below the transition region, the effect this has on the column depths of their formation regions is minimal.

Whereas both non-thermal ionisation and excitation decrease optical depth for Lyman lines, the processes work against each other for the Balmer series and continuum (Fig. 2b). Ionisation continues to decrease the amount of absorbing material, but excitation from the ground state into the second level creates more absorbing material for Balmer radiation and, thus, increases the optical depth. At the depth where the plasma drops below the temperature of full ionisation, optical depths rise steeply (Fig. 2b). However, because a higher proportion of neutral hydrogen (and thus absorbing material with electrons in level 2 of a hydrogen atom) is ionised by stronger beams, this increase of optical depth is reduced for the simulations with higher initial fluxes of electron beams (see Fig. 2b, at column depths around 10¹⁹ cm⁻²).

In a deeper flaring atmosphere the plasma density becomes higher and the hydrogen ionisation degree is lower (Fig. 1, central panels, and Fig. 6). As a result, the excitation of electrons to level 2 of hydrogen bycollisions with beam electrons greatly increases. The number of hydrogen atoms with electrons in level 2 defines the optical thickness of the Balmer transitions (see Appendix C). Thus, stronger beams result in Balmer emission with higher optical thickness (see Fig. 2b, at column depths around 10²¹ −10²¹ cm⁻²).

The excitation by non-thermal beam electrons is an important mechanism for increasing the optical depth in Paschen lines. Non-thermal excitation dominates higher in the atmosphere and is more significant for the Paschen series than for the Balmer one (Fig. 2c). Without the additional excitation supplied by non-thermal collisions there are many fewer electrons in level 3 of hydrogen atoms that can absorb Paschen radiation. In the events excited by stronger electron beams Paschen Pα and Pβ line cores become optically thick (Fig. 3b). Also the formation regions of these lines can expanded closer to the flare’s transition region for more intense events.

In the second part of this investigation, the hydrodynamic and radiative models are calculated for the same beam parameters. The optical depths 5 s after beam onset are shown in Fig. 3 for (a) the Lyman alpha line and (b) the Pβ line, plotted against plasma column depth (cm⁻²). Results are shown for beams with fluxes of 10¹⁰ (red line), 10¹¹ (green line), and 10¹² (blue line) erg cm⁻² s⁻¹.

When the matched radiative and hydrodynamic simulations are used, we can still observe the effects derived from the differential study earlier in this section; i.e. the reduction in optical depths of Lyman lines due to ionisation and excitation by collisions with non-thermal electrons (Fig. 3a) and the increase in Paschen line optical depths due to excitation by beam electrons (see Fig. 3b at column depths deeper than 2 × 10²⁰ cm⁻²). However, understanding of the radiative effects derived from the differential study must be combined with an appreciation of differences in hydrodynamic atmospheres generated by appropriate beams, such as the different column depths of the flare transition layer for each flux (see Fig. 3). This helps to complete our understanding of the effects of the beam on optical depths in the hydrogen lines in real dynamic atmospheres.

Fig. 2 Optical depths calculated in the differential study for (a) the Lyman alpha line; (b) the Hα line; and (c) the Pβ line 5 s afterthe onset of the beam, plotted against column depth (cm−2). Each subfigure shows the optical depth in the central wavelength of the line profile calculated for the thermal response (black line)and for models including excitation and ionisation rates from electron beams with fluxes in the range 109 − 1012 erg cm−2 s−1.

3.2 Line intensity profiles: core, wings, and macro velocity

To analyse emission from a flaring atmosphere in the hydrogen Balmer and Paschen lines and continua we use the full NLTE radiative code for the relevant hydrodynamic atmospheres with matching beam parameters. At first, to discern direct effects of energetic particle beams on the line profiles, we consider emission without including any shifts due to plasma macro velocities as presented inSect. 3.2.1. Subsequently, in Sect. 3.2.2, we consider line profiles with the Doppler shifts included to produce a more complete picture of the resultant line profiles and to compare these line profiles with observations.

3.2.1 Balmer and Paschen lines: core and wing responses

Figure 4 shows the simulated enhancements of the Hα line profiles for 1–5 s after a beam onset (panels a–e, respectively) with intensity (I) presented on the y-axis in erg cm⁻² s⁻¹ sr⁻¹ Å⁻¹ and distancefrom a line centre (λ − λ₀) in Å on thex-axis. Hereafter, the lines show the profiles for the F12 model (blue), F11 model (green), and F10 model (red). When a beam of energetic electrons is injected down through the solar atmosphere from the corona, inelastic collisions between beam electrons and neutral hydrogen causes electrons in hydrogen atom to move into excited states of the upper atomic levels and to become quickly ionised.

This increases the source functions of the hydrogen line and continua so that lines, such as Hα, which are inabsorption in the quiet sun are converted into emission lines in flaring atmospheres. When in emission, these line profiles havea Doppler-broadened core around the central wavelength (see Fig. 4a, at wavelengths between ±0.5 Å). If the optical thickness is high in this transition, there is significant self-absorption in the central wavelength of the line with horns of increased intensity appearing in the near wings of the profile (see Figs. 4b and e at ±0.7 Å and ±0.5 Å in the F12 model). There are also extended wings in these profiles due to the Stark effect. These extended wings can be seen between 1 and 2.4 Å from the central wavelength in Fig. 4.

We examined these profiles in a greater detail, as shown in Table 1, which presents the number density (n_tot), ionisation degree (χ), and ratio of hydrogen atom numbers with electrons in level 3 to those with electrons in 2 ($\frac{n_3}{n_2}$), one second after beam onset. This information is presented for each of the models at the point where the optical depth (τ) of the Hα line is one. Simulations showed that one second after a beam onset all the Hα line cores become optically thick at the depths where the kinetic temperature approaches 18 000 K.

Zharkova & Zharkov (2007) show that electron beams with higher initial fluxes of electron beams result in the flare transition region forming at a greater column depths. This effect can lead to an increase of the column depth and, therefore, densities of the formation regions for hydrogen lines (see Table 1). This can be an essential addition to the influence of radiative processes on resulting intensities in lower transitions (Lyman series) or in opposition to radiative processes acting to form the core higher in the atmosphere in upper transitions (Paschen Series) (see Sect. 3.1). In the Balmer line formation regions there are greater non-thermal excitation and ionisation rates for beams with higher initial fluxes. Therefore, there is more hydrogen atoms with electrons in the upper excited levels (Table 1). The greater ratio of atoms with electrons in the upper levels of Balmer lines to those in the lower state (n = 2) explains an increase of the emission in Balmer lines for stronger beams (Fig. 4).

The line broadening results from collisional broadening caused by ambient and beam electrons and the Stark effect caused by local electric fields produced by these electrons. These are both dependent on a number of free electrons. In stronger flares the hydrodynamic response to a beam injection shifts the line formation regions into deeper layers with higher densities. Combining this with the increased ionisation degree that results from a stronger beam (evident in Table 1 and Fig. 6), we have ideal conditions for the broadened lines with stronger wings as the initial energy flux of a beam causing it is increased. Figure 5 presents the simulated enhancements in line profiles in the Hα line (top row), Hβ line (middle row), and Hγ line (bottom row) at 4 and 5 s after the onset of a beam. There is an evident pattern in the Balmer series showing the beams with larger initialfluxes producing broader, less pronounced cores of emission lines and more extended wings (Figs. 4 and 5) throughout the beam injection phase.

Figure 6 shows the ionisation degree plotted against column depth for the thermal (black line), F10, F11, and F12 models.Ionisation rates below the flare transition layer are dominated by collisions between hydrogen atoms and non-thermal electrons (Zharkova & Kobylinskii 1993). The ionisation of hydrogen increases owing to non-thermal inelastic collisions with beam electrons compared to the pure thermal case. Additionally, excitation caused by non-thermal electrons increases the thermal ionisation rates from higher atomic excited states. This leads to a wing intensity increase that is strongly dependent on the initial energy flux of the beam, which is evident across all the Balmer lines at all times (Figs. 4 and 5).

The pattern of wing intensity dependence on non-thermal ionisation and excitation by beam electrons is enhanced in the differential studies, where the beam parameters are varied in the radiative code while using the same hydrodynamic response. When we use the joint solutions of radiative and hydrodynamic problems, the wing pattern is still observed. The enhanced Hα wings without strong red-shifts are well known and observed as “moustaches” (Rust & Keil 1992; Zharkova & Kashapova 2005). However, after 4 and 5 s the Hα line profile produced bythe F12 flare becomes broadened so much by the Stark effect and shifted to the red wing by Doppler effect to wavelengths that extend well beyond the 2.4 Å displacement from the line centre plotted in simulated figures (Fig 4).

At each column depth where the beam causes increased emission, the optical depths linked to that column depth becomes lower for the Balmer lines through the series from Hα to Hβ to Hγ. As a result the emission profile of the Hα line shows greater core self-absorption than the Hβ line and both show greater core self-absorption than the Hγ line profile (Figs. 5a–c).

Fig. 3 Optical depths in coupled hydrodynamic and radiative models, 5 s after beam onset in a) the Lyman alpha line and b) the Pβ line, plotted against plasma column depth (cm−2). Simulations are made with beams with fluxes of 1010 (red line), 1011 (green line), and 1012 (blue line) erg cm−2 s−1.

Fig. 4 Profiles of Hα line intensityenhancements, in erg cm−2 s−1 sr−1 Å−1, plotted vs. distance in Å from the line central wavelength for times from 1 to 5 s in panels (a) to (e), respectively.Results are shown for the F12 (blue line), F11 (green line), and F10 (red line) beam models. The colour scheme for the lines representing these models is also used in Figs. 5 to 10.

Fig. 5 Intensity enhancements in erg cm−2 s−1 sr−1 Å−1 against distance from line centre in Å for the Hα (top panels, a and b), Balmer beta (Hβ) (middle panels, c and d), and Hγ (bottom panels, e and f) lines 4 s (left panels) and 5 s (right panels) after the beam injection begins.

3.2.2 Macro velocity: effects on line profiles

We generated the hydrogen emission profiles for lines with a central wavelength λ₀ using the modelled depth points in the flaring atmosphere. The emission from each layer has its own profile and associated Doppler shift $\Delta {\lambda }_i=\frac{V_m(i)}{c}\mu {\lambda }_0$ from the layer i, due to its macro-velocity V _m (i). To smooth the contributions of the various layers considered, cubic splines were applied to the line profile merged from these layers. Figure 7 shows the Hα line profile intensity enhancements and Fig. 9 the Pα for the F10 (red), F11 (green), and F12 (blue) models at 2, 3, 4, and 5 s from panel (a) to (d), respectively.

Large errors are known to occur in Doppler velocity calculations using the bisector methods when working with the line profiles that containlarge deformations (Deng et al. 2013). To avoid these errors in our simulations, it is practical to use the maxima in the emissionintensity profiles. Red-shifts occur when the hydrodynamic shock moving downward enters the line formation region. In the Hα line there isa slight skewing of the core to the red wing for the F10 flare model at 3 s and 4 s (Figs. 7b and c, red line). At 5 s the red-shifted core peak is visible at a wavelength of 1 Å from the line centre wavelength and the blue wing forms a slight horn at − 0.5 Å. The peak shift suggests a Doppler velocity of 4.5 × 10⁶ cm s⁻¹ (45 km s⁻¹). At the column depths that contribute most to the emission, we find macro velocities in the range 3.91−5.36 × 10⁶ cm s⁻¹.

By comparing profiles of the F10 flare at 5 s and the F11 flare at 3 s (Fig. 7d, red line and Fig. 7b, green line), we see that a + 0.8 Å shift in the maximum of the line profile occurs earlier in the beam injection phase in the F11 case. Likewise, comparing the F12 and F11 cases at3 and 5 s, respectively (Fig. 7d, green line and Fig. 7b, blue line), we see that the delay between beginning of the impulsive phase and the time when the hydrodynamic shock enters the Hα formation region decreases with the increasing initial energy flux of a beam. This is due to larger macro velocities generated by hydrodynamic responses of flaring atmospheres to the injection of beam electrons with larger initial fluxes.

In the F12 flare model the core is shifted from the central wavelength of the spectral line out of the wavelength window used in our simulation 2–3 s after the beam injection phase begins, as the hydrodynamic shock hits the Hα line core formation region. Thus, one would expect the macro velocities in the line formation region could exceed the maximum value captured by our wavelength window (±2 Å), and reach 1.05 × 10⁷ cm s⁻¹, or 3.5 Å, for this model, which is beyond the 2 Å window in the figure. The main formation region (identified via the contribution functions) has the temperatures around 18 000 K and associated macro velocities of around 2.1 × 10⁷ cm s⁻¹.

Hence, summarising the above we can conclude that the Hα red wing enhancements in these simulations are clearly associated with the plasma downflows, which can be found entering the formation region between 2 and 3 s. The scale of the line deformation remains less noticeable for weaker flares, which generate lower maximum macro velocities in their downflows. Maximum values of downflow velocities in the Hα line formation region are found at the peak of the initial flux, after 5 s. For a comparison of the above simulations of Hα line emission with observations using SST we refer to Sect. 4.1 and to Druett et al. (2017).

If the core emission is red-shifted thus having a smaller optical thickness appropriate to the wing wavelength, this leaves the wing intensity contribution from this region overlapping with the core emission coming from greater optical depths in the central wavelength. This results in a self-absorbed profile (Figs. 7c,d, and 9c,d, and the F11 and F12 models). However, if there is a very strong intensity in the line wing generated in the core formation region, for example due to Stark’s wings and collisional effects, as was the case in the F12 Hα line, then this can compensate for the red-shifted core intensity being obscured by these effects (Figs. 7c,d).

Figure 8a shows the enhancements to the Hγ line profiles for the F11 model from t = 2 s (light green line) to t = 5 s (black line). During the beam injection the Hβ and Hγ line profiles show similar temporal evolution to Hα line profiles. The Hγ line profiles in the F11 model show Doppler red-shifts with the emission enhancements peaking for Doppler velocities around 150 km s⁻¹ at 4 s and 5 s.

Figure 8b shows the enhancements to the profiles of the first three lines in the Balmer series (Hα in orange, Hβ in red, and Hγ in black) calculated for the F10 model at t = 5. At 5 s The core of the Hα line is red-shifted because the hydrodynamic shock travelling downward has entered the Hα core formation region (Fig. 8, orange line), whereas only around half of the core emission produced in the Hβ line has been red-shifted (Fig. 8, red line). This is because at 5 s the core formation region of Hβ extends to depths below the hydrodynamic shock. Only a small proportion of the core emission in the Hγ line is red-shifted (Fig. 8, black line) as the majority of core intensity is produced from depths below the large downward macro velocities. This figure highlights the greater importance of contributions to emission from greater column depths in the higher lines of the Balmer series. In the F10 model at 5 s the ratio of the maximum intensity in the red wing compared to that in the centre or blue horn is 2.8 for Hα, 1.1 for Hβ, and 0.62 for Hγ. Hence, using simultaneous, high cadence observations of Hα, Hβ, and Hγ line profiles during a beam injection it is plausible to track the propagation of a hydrodynamic shock through the chromosphere by comparing the line profiles.

Additional information about the beam can be provided by Paschen lines, which are formed over a broader range of the chromospheric depths if affected by beams. We examine the results for the Pα line (Fig. 9).

By examining the contribution functions for the emission it can be seen that the higher part of the Pα formation region overlaps with the formation region of Balmer lines. This causes a co-temporal onset of the red-shifts appearing in the line profiles for each flare model. The fact that the formation region of Paschen line cores extends to greater column depths results in a smaller proportion of the emission from the line core being red-shifted at these times. This can be seen by comparing the line profiles from 5 s for the F10 profile of the Hα (Fig. 7d, red line) and Pα (Fig. 9d, red line) lines. Indeed, in the F10 flare model the hydrodynamic shock caused by the beam does not penetrate into a significant proportion of the Pα formation region at 5 s. The effect of the hydrodynamic shock entering the formation region can be seen in the Doppler shift of the bowing of the red wing and the shape of the line core at this time (Fig. 9d, red line).

However, using the same comparison for the F11 and F12 models, one can observe that both the Pα and Hα line have a large proportion of their line cores red-shifted away from the central wavelength. This is because, as shown in Sect. 3.1 and Figs. 2 and 3, the Paschen line cores become highly optically thick in the case of powerful beams. Therefore, in strong flares (with powerful beams) the Paschen and Balmer line profiles reflect the conditions in the upper flaring chromosphere, where a similar proportion of the line core is Doppler shifted from the central wavelength. However, for weaker beams with optically thinner Paschen line cores one observes a smaller fraction of the Paschen line core being Doppler shifted away from the central wavelength than in the Balmer lines at the same time.

Fig. 6 Ionisation degree at 4 (a) and 5 (b) seconds plotted against plasma column depth (cm−2) for the thermal (black line), F10 (red line), F11 (green line), and F12 (blue line) models.

Fig. 7 Hα line intensity enhancements with macro-velocity effects included for the F10 (red), F11 (green), and F12 (blue) models. Profiles are shown for a) t = 2 s; b) t = 3 s; c) t = 4 s; and d) t = 5 s in erg cm−2 s−1 sr−1 Å−1 against distance from line centre in Å.

Fig. 8 Panel a: spectral profiles of the Hγ line for t = 2 s (light green) to t = 5 s (black). Panel b: Hα, Hβ, and Hγ line profiles at t = 5 for the F10 model. Intensities (I) are simulated in erg cm−2 s−1 sr−1 Å−1 against distance from line centre (λ − λ0) in Å.

3.3 Simulations of Balmer and Paschen Continua

The simulated results for the Balmer and Paschen continua (Fig. 10) reveal the profiles with intensities that are highly sensitive to the initial energy flux of the beam. This is caused by the fact that hydrogen continuum radiation is generated by recombination of free electrons with an ionised hydrogen atom (proton); the intensity of emission produced by this process depends strongly on the ionisation degree of the plasma. In our models the ionisation rates are dominated by collisions of hydrogen atoms with non-thermal electrons (see Sect. 2.2.1), increasing within a very short timescale the ambient ionisation by orders of magnitude and the optical thickness of continuum emission in Lyman and, sometimes, Balmer continua. This enhances the role of radiative transfer in continuum, which governs hydrogen ionisation throughout a flaring atmosphere (Zharkova & Kobylinskii 1993) for rather long times while slowly emitting continuous radiation from a thin atmospheric layer where τ_1c ≈ 1.

For all beams the simulated Balmer continuum emission was the most effectively enhanced and contained a greater intensity in the continuum head than any other hydrogen continuum. The Lyman continuum head is optically thick. Therefore, contributions do not emerge from a large fraction of the depths, at which the beam causes increased ionisation, thus the emerging intensity of Lyman continuum emission is less enhanced than the Balmer continuum.

To investigate the relationships between the intensities of emission in the other continua we used the ratios of the intensity in the continuum heads. Tables 2 to 4 list the ratios of the emission intensity in the Balmer continuum head to those in the Paschen, Bracket, and Pfund continuum heads, respectively. The top rows lists results for the F12 model, the middle rows contain results for the F11 model, and the bottom rows show the F10 model data. The columns indicate the time since the beam injection began, from 1 to 5 s. Over the first 5 s of the injection the ratio of Balmer continuum head intensity to other continua head intensities is reduced by between 18 and 20% for the F12 beam model. For the F11 model the ratio is reduced by between 3 and 4%. For the F10 model there is a negligible change in the ratio over time (Tables 2, 3, and 4). This decrease, or its absence, results from the amount of increase in the ionisation rates by collisions with beam electrons at varying depths in the various beam models.

In all the simulated models, the ambient hydrogen is efficiently ionised for a large proportion of the depths that make up the Balmer continuum formation region, but extend to a smaller proportion of the formation regions for the other continua. In models with stronger beams, this increase extends downward into the depths where subordinate continua are more effectively contributed to, reducingthe ratio of the subordinate continuum head intensities compared to the Balmer continuum head intensity.

With similar reasoning one expects that the same pattern of ratios is present between Paschen head intensity and the subordinate Brackett head intensities. Table 5 lists the results in the same format as Tables 2 to 4 and although there is a similar relationship, the magnitude of the effect is negligible even in the F12 flare, as the flaring chromosphere does not span a large portion of the formation regions of either continuum (Table 5).

The Balmer jump is the difference in intensities on either side of the Balmer continuum head wavelength (3646 Å) in the continuum spectrum of a star. In order to assess how the initial flux of an electron beam affects the Balmer jump, we examined our results at the Balmer continuum head wavelength for both Balmer and Paschen continuum radiation. The ratio of intensity of the Balmer continuum to the intensity in the Paschen continuum at the wavelength of the Balmer continuum head (BP ratio) does not show a clear pattern across all flare models (Table 6). This results from the fact that the profiles of the continua are affected somewhat by the detailed hydrodynamic response to beam deposition and by the ionisation degree. However, we see that the ratio is lowered by a greater amount for the F12 model than the F10 and F11 models. This suggests that one would expect to see a smaller Balmer jump developing during beam injection, in the case of beams with high initial fluxes because of the strong ionisation increase in the photosphere.

Fig. 9 Pα line profiles vs. distance in Å from the line centre wavelength calculated at a) t = 2; b) t = 3; c) t = 4; and d) t = 5, in erg cm−2 s−1 sr−1 Å−1.

Fig. 10 Intensities of Balmer continuum enhancement at a) t = 1; c) t = 3; e) t = 5, and Paschen continuum enhancement at b) t = 1; d) t = 3; f) t = 5, in erg cm−2 s−1 sr−1 Å−1 against wavelength in Å.

3.4 Formation regions of Balmer and Paschen continua

Insight into the origin of continuous emission for Balmer and Paschen continua can be gained by examining the contribution functions of the hydrogen continua. We only considered Paschen continuum as the source of WL flares. Figure 11 shows the Paschen continuum contribution functions at t = 5 s for the F10, F11, and F12 model simulations in panels a and b, c and d, and e and f, respectively. The intensity contributions are plotted on the z-axis, with wavelength λ, and the column depth in logarithmic scale log₁₀ξ on the x- and y-axes, respectively.The left panels show the contribution functions for the models with radiative rates including non-thermal excitation and ionisation by a beam of electrons. Those on the right show the results for “thermal flares” using the same hydrodynamic model as shown in the left panel, but without non-thermal beam excitation and ionisation included.

In the F10 beam model, non-thermal ionisation causes a great increase in continuum contributions higher in the flaring chromosphere for column depths in a broad range from around ξ = 10¹⁹ to 10²² cm⁻² (compare the beam and thermal model contribution functions in Figs. 11a and b, respectively, for column depths from ξ = 10¹⁹ to 10²² cm⁻²). However, the contributions from the deeper atmosphere are barely enhanced compared to the thermal model (see Figs. 11a,b, ξ = 10²² cm⁻²). This is because, for the beam with a low initial flux, the beam electron densities are reduced to a point wherethey are not ionising a large amount of the neutral hydrogen at a column depth of ξ = 10²² cm⁻². Additionallythere is minimal backwarming of the photosphere from the radiative response in the chromosphere.

In the F11 model, contributions from the mid-chromosphere are even more greatly enhanced as a direct resultof the increased non-thermal ionisation rates at these depths. Additionally, we observe an increase in contributions from deep in the atmosphere (Figs. 11c,d). This results from both radiative transfer and a more significant number of beam electrons penetrating to depths greater than ξ = 10²² cm⁻².

The F12 model generates a very strong Paschen emission, or WL flare, coinciding with large increase of Balmer continuum radiation produced in the flaring chromosphere. A combination of the effects of radiative transfer with strong non-thermal ionisation in the lower atmosphere produces the extremely large increase of upper and lower atmospheric contributions (Figs. 11e,f). The F12 flare model produces a non-negligible increase in a temperature in the deep atmosphere (photosphere) (Fig. 1, left panels), but this is not the source of WL enhancement suggested by the contribution functions presented here. The thermal flare model (only lacking non-thermal excitation and ionisation) produces a much smaller increase in the Paschen continuum, or WL, compared to the model considering the effect of a powerful beam of non-thermal electrons. Hence only the presence of beam electrons can explain the occurrence of WL flares simultaneously with HXR emission reported by many authors (Hudson 1972; Martínez Oliveros et al. 2012).

The Paschen continuum head intensity ratios for the model with the beam to the thermal flare are of factor 3.53 for the F10 model, 14.4 for the F11 model, and 15.5 for the F12 model. A similar study was conducted for the Balmer continuum, which produced similar results for the ratios of head intensity in models with beams to those for the thermal model, 3.55, 14.0, and 13.0 for the F10, F11, and F12 flares, respectively. The variation of Balmer continuum contribution functions from different column depths also showed similar patterns to those for the Paschen continuum in each simulation (Fig. 11). The differences were (1) that the contribution functions had higher magnitudes for the Balmer continuum and (2) that there were some differences in the variation of the profiles with wavelength. The results for the Paschen continuum were selected for display in Fig. 11 to aid the comparison with WL observations carried out in Sect. 4.3.

Thus, one can observe that electron beams are highly effective agents for the production of WL emission co-temporally with the beam onset in a flaring atmosphere at higher atmospheric depths, compared with the pure thermal heating. Comparison with observations of Balmer continuum and WL in flares is carried out below, in Sects. 4.2 and 4.3.

Fig. 11 Paschen continuum contribution functions (see Sect. 2.3.3) calculated at t = 5 s. The z-axis shows the contribution function value, the x- and y-axes are in terms of wavelength, λ and the logarithm of column depth log10ξ. The top row shows results for the F10 flare model, the central row shows data for the F11 flare model, and the bottom row shows data for the F12 model. The left panels show the models including excitation and ionisation of hydrogen atoms by non-thermal beam electrons and the right panels represents the “thermal flare” simulations without excitation and ionisation by beam electrons.

4 Comparison with observations

4.1 Comparison with Hα line observations

4.1.1 Line profiles

We present a comparison of our simulations with the observations of Hα line profiles taken at a flare onset. Wuelser & Marti (1989) reported the observations of flare kernels in Hα wavelengths for the M1 solar flare on May 24, 1987 taken with the Specola Solare Ticinese at Locarno-Monti (see Figs. 12c,d). In particular, the profile presented for the A1 kernel at the peak of the HXR burst shows a little sign of any pre-flare heating. The normalised Hα profile extracted from data in the A1 kernel at the time of the peak in the HXR signal (15:24:55 UT) shown in their Fig. 3 was digitally recorded (Fig. 12c, green line). The quiet Sun profile presented intheir Fig. 3 (Fig. 12c, grey line) was subtracted to find the enhancement in Hα emission (Fig. 12c, pink line). The resulting profiles are shown with normalised intensity on the y-axis and the wavelength relative to the Hα line central wavelength is shown on the x-axis.

As shown in a number of studies (Jess et al. 2008; Druett et al. 2017; Kowalski et al. 2017) the small (<1^″) impulsive flaring kernels associated with beam electrons are also particularly associated with chromospheric line profiles with asymmetric red wing emission. The pixel size of the observations by Wuelser & Marti (1989) was 3^″ × 3^″; this implied thatthe authors were unable to resolve the locations of individual beam injections, which covered much wider areas instead produced by the background solar atmosphere. This fact explains the large enhancement in the Hα emission in a central wavelength shown in Fig. 12c because only a small part of the region used to generate the emission profile undergoes a large red-shift due to the hydrodynamic shocks caused by non-thermal beam injections.

Therefore, to identify the emission from regions with beam injections and avoid accounting for the background, we had to calculate the excess emission in the red wing at the wavelength λ₀ + Δλ as shown in Fig. 12d by subtracting the normalised intensity enhancement of the profile at a wavelength of ${\lambda }_0-\left|\Delta \lambda \right|$ from the intensity enhancement at λ₀ + Δλ. Of course, such comparisons should be approached with caution, but they are preferable to those simulations directly compared from 1D flux tube models with the complete emission profiles or background subtracted profiles that have been integrated over a large area of the whole active region. This is because during the subtraction we reduced the background emission and picked out the feature of line profiles specifically associated with the beam injection rather than looking at the kernel emission aggregated with emission from bright ribbons and neighbouring thermal contributions in addition to that from a beam. Selected simulated Hα line profiles that show red wing excesses inside the modelled wavelength window are shown in Fig. 12a for the F10 model at 5 s (red line), the F11 model at 3 and 5 s (light and dark green lines), and the F12 model at 3 s (blue line). The red wing excesses for these profiles is shown in Fig. 12b. Both the graphs are plotted with enhancement to emission intensity on the y-axis, against wavelength relative to the central frequency of the Hα line, λ₀, in Å on the x-axis.

The A1 kernel (Wuelser & Marti 1989) shows a great red wing enhancement in the Hα line at 15:24:55 UT, co-temporal with the HXR peak intensity. The red-shifted emission had an intensity peak at around 2 Å from the line centre wavelength at these times, corresponding to a downflow of around 90 km s⁻¹ that represents velocities between the simulated maximum downflow velocities in the F10 and F11 hydrodynamic models. Additionally, there was enhancement in emission extending up to around 4 Å, suggesting that strong Stark wings had been generated. Both of these features are evident in the profile of the excess emission in the red wing (Fig. 12d). This description fits the Hα line profile enhancements and red wing excesses for beams with fluxes between 10¹⁰ erg cm⁻² s⁻¹ and 10¹¹ erg cm⁻² s⁻¹.

Additionally, a peak in the red excess at around 1–3 Å also indicates the presence of a shock caused by a beam with a lower initial flux in the kernel area. We used the profile produced by the F10 model at 5 s as a template for the shock resulting from an electron beam injection with the initial flux slightly lower than 10¹⁰ erg cm⁻² s⁻¹ to fit this peak. Then, we combined the red excess from this profile linearly with the red excess modelled for a flare with the incident flux somewhat lower than the F11 model. The resulting red excess is shown by the black line in Fig. 12d and fits the observed red wing emission excess very closely. This suggests that the large kernel area used in the observation could have injections from at least two beams with fluxes around 10¹⁰ −10¹¹ erg cm⁻² s⁻¹.

Recently Druett et al. (2017) reported a very close agreement between the simulated Hα line profiles and the Hα line profiles observed 7 s after in a C1.5 flare onset on June 30, 2013 using the SST (Scharmer et al. 2003). The observations exhibit a large Doppler shift of ~1 Å, which is in agreement with the timing and extent of the enhancement simulated in the red wing of the Hα line for the beam injection with an initial flux 10¹⁰ erg cm⁻² s⁻¹. This successful verification of the Hα line profile highlights a need of high resolution Hα observations with a wider spectral window, covering at least ±2 Å around the line central wavelength.

Fig. 12 Panel a: simulated Hα profiles showing red wing enhancement for the F10 (red line), F11 (green line), and F12 (blue line) models. Normalised intensity enhancement is on the y-axis, and wavelength relative to the Hα profile line centre on the x-axis. Panel b: excess of emission in the red wing of the profiles in panel “a” relative to the emission in the blue wing. Panel c: normalised Hα profiles from the M1 solar flare on May 24, 1987 reported in Wuelser & Marti (1989). The normalised Hα profile (green line), quiet Sun background level (grey line), and enhancement in emission (pink line) are shown for the A1 kernel at the time of the maximum of the HXR (15:24:55 UT). The y-axis shows normalised intensity and the wavelength relative to the Hα line central wavelength is shown on the x-axis. Panel d: excess emission in the red wing compared to the blue wing for the A1 kernel at 15:24:55 UT (green line). The red line shows a simulated fit of this excess based on a linear combination of contributions based on the F10 and F11 models.

4.1.2 Temporal variations

The observations of Kaempfer & Magun (1983) report three kernels for the flare on June 17, 1982. One of the kernels (L) has co-temporal peaks in HXR and Hα intensities in the 0.25 Å window. The other two kernels (M and R) exhibit a 30 s delay in the maximum of emission intensity of the Hα line. From this, the authors suggested that different energy transport mechanisms are at play in the different flare kernels. Although, this cause cannot be ruled out, imagining the observations by Wuelser & Marti (1989) conducted using a narrow wavelength band, such as the 0.25 Å window used in Kaempfer & Magun (1983), provides a cautionary tale.

The maximum of the HXR spectra occurred at around 15:25:00 UT, according to Fig. 3 in Wuelser & Marti (1989). In the same figure it is shown that a maximum of the emission intensities in the red wing at 1.5 and 3 Å occurred at around 15:25:09 UT in the A1 kernel. For the same kernel there was a second set of intensity maxima in the line centre and the red and blue wings of the Hα line profile at around 15:25:30 UT, which would be the only peak recorded using a spectral window of 0.25 Å around the line centre. Therefore, a 30 s delay between the peak of the response in the Hα line core and the peak of the HXR would have been reported.

This time delay occurs because the emission profile has a maximum value that is Doppler shifted into the red wing, well beyond 0.25 Å (see, for example, the simulated profiles in Fig. 7). The B1 kernel with its pre-heated enhancements exhibited a very slight delay between the HXR peak intensity and that in Hα line profile. It is perfectly possible, rather it is a prediction of our models, for delays such those reported in Kaempfer & Magun (1983) to be evident at the onset of a beam injection. This delay results from the red-shifts caused by plasma downflows resulting from a hydrodynamic response of flaring atmosphere to an non-thermal beam injection. The Hα line emission is shifted to the red wing and the filter with the given size observes only the blue wing of this shifted profile. The line profile only returns to the profile centred on the central wavelength, which can be observed by the given filter, at later times after the beam is switched off and the plasma stops moving downward.

We suggest that for capturing the full dynamics of Hα line profiles during the impulsive phase of solar flares larger wavelength windows are essential and recommended in new instruments.

4.2 Comparison with the observations of Balmer continuum enhancement

Heinzel & Kleint (2014) identified Balmer continuum enhancement during the X1 class flare on March 29, 2014 using IRIS, at a disk position of μ = 0.83 and a wavelength of 2826 Å. They report an enhancement of 4.1 × 10⁵ erg cm⁻² s⁻¹ sr⁻¹ Å⁻¹ over the background quiet Sun contribution of 3.7 × 10⁵ erg cm⁻² s⁻¹ sr⁻¹ Å⁻¹. Heinzel & Kleint (2014) state that the RHESSI data suggests, to an order of magnitude, the beam flux for electrons with energies greater than 20 keV was 10¹¹erg cm⁻² s⁻¹ with spectral index = 5.

To model this observation Heinzel et al. (2016) simulated Balmer continuum emission using the FLARIX code (Varady et al. 2010). At the peak of an initial flux, F₀ = 4.5 × 10¹⁰ erg cm⁻² s⁻¹ sr⁻¹ Å⁻¹ FLARIX produces Balmer continuum enhancement of around 6 × 10⁴ erg cm⁻² s⁻¹ sr⁻¹ Å⁻¹, which is in reasonable agreement with the observation. To evaluate our model, we examined the simulated Balmer continuum at a wavelength of 2826.5 Å and at the peak of the beam injection (t = 5 s), when a significant amount of the ionisation occurs owing to non-thermal inelastic collisions. With μ = 0.83, the F11 model produces intensities of Balmer continuum that are even closer fits to the observation than those produced using FLARIX 8.7 × 10⁴ and 3.13 × 10⁵ erg cm⁻² s⁻¹ sr⁻¹ Å⁻¹ for γ = 5 and γ = 3, respectively. The Balmer continuum intensity produced using FLARIX was very close to the value produced with our F10 model.

Sources of errors between the simulated emission intensities and those observed include: (1) The duration of the beam injection appears to be significantly longer for the observations discussed than the 10 s used in our model. The HXR spectra in Fig. 4 of Heinzel & Kleint (2014) show HXR bursts lasting for times in the order of minutes. (2) The loop footpoint areas are unresolved in RHESSI data, giving lower bounds of the initial energy flux for a beam. (3) The calibration process used to convert the observed data number (DN) units into the CGS units. (4) There is a difference in timing of the observed and simulated continua, with our model representing the first seconds of the flare onset and with an 8 s exposure time and a 75 s cadence in the observations used by Heinzel & Kleint (2014).

4.3 Comparison with the observations of white light enhancement

Martínez Oliveros et al. (2012) reported observations of an M3.5 class flare, which occurred on the solar limb, on February 24, 2011 at 07:35 UT. Data from the 6173 Å channel of the Solar Dynamics Observatory (SDO) were used to analyse the WL enhancement. HXR spectra were generated from RHESSI data in the 30–80 KeV channel over a 45 s exposure that was co-temporal with enhancements in WL. Simultaneous imaging from the Solar-TErrestrial RElations Observatory (STEREO) provided heliographic coordinates of the flare footpoints. A height scale relative to the photosphere was defined using τ = 1 at 5000 Å as the surface of the photosphere.

The HXR sources had centroids located 420 km and 210 km above the quiet Sun photosphere and the continuum sources had centroids with heights of 230 km and 160 km in the northern and southern footpoints of the flare, respectively. The uncertainty in the measurements of the heights of the HXR sources is ± 240 km and for the heights of the WL sources the uncertainty is ± 100 km. Thus the sources of the HXR and WL were co-spatial to within instrumental resolution. The authors state that their result, “strongly associates the WL continuum enhancement with the collisional losses of the non-thermal electrons observed via bremsstrahlung HXRs in the impulsive phase of the flare” (Martínez Oliveros et al. 2012). They also note that the observed heights of the sources occurred well below the 800 km height that represents the stopping distance of 50 keV electrons found using the quiet Sun atmosphere of Fontenla et al. (2009). We compared these observations with the strongest beam model presented in this paper (F12), which has the greatest sweeping of chromospheric plasma towards the photosphere (see Figs. 1k,l, blue lines), in order to assess whether a non-thermal electron beam is a feasible agent for the production of the observed sources of deep, co-temporal HXR and WL.

The column depths of the simulated WL signal, from the contribution functions due to Paschen continuum recombination in the F12 model at the peak of HXR, are shown in Fig. 11e. At 6173 Å the secondary source, produced in the flaring chromosphere, is seen to extend from the column depth of 1.28× 10²⁰ cm⁻² to 2.52× 10²² cm⁻² with a peak at 3.44× 10²¹ cm⁻² (Fig. 11e). The heights in the F12 hydrodynamic model (Fig. 1, blue lines) are used to identify the vertical heights of these positions above the quiet Sun photospheric level. The top and bottom of this WL signal are 340 km and 180 km above the quiet sun, respectively, (to the nearest 10 km) with the centroid of the source at a height of 260 km. Thus, for a strong, hard electron beam such as in the F12 beam model, the contribution functions of the simulated Paschen continuum radiation is consistent with the observed height of WL in Martínez Oliveros et al. (2012).

The stopping depths of beam electrons with energies of 30 and 80 keV due to collisional losses are calculated using the method presented in Zharkova & Gordovskyy (2005). The column depths calculated are 1.67 × 10²⁰ cm⁻² and 1.19 × 10²¹ cm⁻², respectively.The heights of these column depths in the F12 model are 320 km and 280 km, respectively, at the peak of the beam flux. Hence the stopping depths of the electrons due to collisional losses in Zharkova & Gordovskyy (2005) are also consistent with the simulated column depths at the height of the signal in the 30–80 keV RHESSI channel used in Martínez Oliveros et al. (2012).

Krucker et al. (2015) report the heights of WL (λ = 6173 Å) and HXR sources measured for three limb flares at the peak of the HXR signal in the 30–100 keV range. The M1.7 class flare presented in their work has the HXR profile that most closely matches the form of the beams used in our models: little pre-heating and a sharp HXR peak followed by a swift decline without repeated injections. The WL source was observed at a height of 799 ± 70 km above the photosphere and had a radial extent of ~652 km for a full width half maximum of intensity. The F11 model has a centroid of WL at height 771 km with an extent of around 984 km for the full width half maximum (ignoring the deeper, photospheric signal, Fig. 11c) at the time of the HXR maximum. Therefore, this model presents good agreement with the observations of WL source height for the M1.7 flare and a reasonable agreement with the observations of the other flares presented in Krucker et al. (2015). It is also possible to alter the extent and height of the WL sources in our models by varying the flux and spectral index of the beam.

5 Discussion and conclusions

Our simulations illustrate that during the impulsive phase of solar flares the Balmer and Paschen line profiles are strongly affected by a dual action of an injected electron beam on a flaring atmosphere via hydrodynamic and radiative responses.

Hydrodynamics: We show that within seconds the injected electron beam converts the quiet Sun chromosphere into a flaring atmosphere with its own flaring corona, transition region, and chromosphere (see Fig. 1). The column depth of the flare transition region that forms through the hydrodynamic response to a non-thermal electron beam is greater by 2 orders of magnitude for the initial flux of the beam 10¹² erg cm⁻² s⁻¹ than it is for a beam with flux 10¹⁰ erg cm⁻² s⁻¹ (Fig. 1, left panels). The differences in column depth of the flare transition region for beams with like parameters are comparable to other models with beam heating functions derived using a continuity equation (Fisher et al. 1985).

The heating caused by beam electrons sweeps the coronal plasma to deeper atmospheric levels, creating a hydrodynamic shock travelling towards the photosphere and beneath. This condensation moves downward as a shock acting like a hammer or piston with the velocities scaled by the initial flux of a beam (Fig. 1, right panels). The sweeping of the plasma and beam in our code also results in the upward motion of the chromospheric plasma to the flaring corona as either gentle (F10 model) or explosive (F12 model) evaporation (see Fig. 1, right panels). This evaporation continues for hundreds of seconds after the beam is switched off, reaching flow speeds that are comparable with those in Fisher et al. (1985).

In orderto put our results in context with the previous hydrodynamic and radiative model (Fisher et al. 1985; Allred et al. 2005), we discuss below some basic discrepancies in these models and their implications on the simulation results. Since the 1970s there were two analytical approaches for electron kinetics. The first approach was revealed by Brown (1971) and Emslie (1978); the second was found by Syrovatskii & Shmeleva (1972), considering pure collisional losses, and was updated recently by Dobranskis & Zharkova (2015); Zharkova & Dobranskis (2016), which considered collisional and Ohmic losses. The Brown and Emslie solutions are restricted to the upper chromosphere because of the FCA limitations discussed in the Introduction. Models using the flux conservation equation approach (FCA) to beam heating by Coulomb collisions, such as Allred et al. (2005), have lower velocity of downward moving shocks below the transition region.

In this paper, we use the heating functions by Syrovatskii & Shmeleva (1972) found from CEA, which show that beam electrons can easily reach the lower chromosphere and photosphere because the thermalisation of beam electrons happens in the chromosphere or even the photosphere at the collisional or Ohmic stopping depths for electrons with given energies; see Table 1 of Zharkova & Gordovskyy (2006) or tables in Zharkova & Dobranskis (2016) for a comparison of stopping depths for collisional and Ohmic losses. The macro velocities of the downward moving hydrodynamic shocks in the F10, F11, and F12 HYDRO2GEN models are similar to those in the models presented in Fisher et al. (1985) for beams of electrons with equivalent parameters and can reach the values derived from observations by Lee et al. (2017). This is essentially different from the work of Allred et al. (2005), which uses the heating function by electrons beams derived by Emslie (1978) and Nagai & Emslie (1984) truncated in the top flaring atmosphere to keep a flaring corona not overheated (see Fisher et al. 1985)and at the upper chromospheric column depths just before the stopping depth of the electrons with lower cut-off energy to avoid the infinity problem (Mauas & Gómez 1997; Kontar et al. 2011, p6, footnote). Nevertheless, the beam heating functions based on CEA (Somov et al. 1981; Zharkova & Zharkov 2007) avoid such problems.

With these CEA heating functions the hydrodynamic modelling shows the downward chromospheric velocities of 50 km s⁻¹ easily achievable at the very first 5 s after a beam onset for a moderate F10 beam (our Fig. 4b). However, the simulation for the F10 model by Allred et al. (2005) does not show any significant shocks moving downward in the chromosphere (see their Fig. 8). This is because in the FCA hydrodynamic model of Allred et al. (2005), the maximal downward velocity of the shocks generated in the chromosphere is 40 km s⁻¹ even for the electron flux 10 times higher (F11 model, shown in their Fig. 9) and this velocity is reached much later in time than with the HYDRO2GEN approach.

Hα radiative response: Hydrogen line emission originates in the low-temperature condensation of a flaring chromosphere below the transition region of the flare. Therefore, in the HYDRO2GEN models the beams with greater initial fluxes push the chromospheric plasma of the quiet Sun, or the hydrogen line formation regions, to greater column depths and larger densities (see Table 1). However, our simulations show that these beams still do not reach the regions with very high densities (10¹⁵−10¹⁶ cm⁻³), in which the negative hydrogen ions appear (Aboudarham & Henoux 1987). Electrons beams with larger initial fluxes cause higher non-thermal excitation to the upper atomic states of the transitions in the Balmer and Paschen series, raising the ratio of the electron abundances in the upper state compared to those in the lower state (see Table 1). This increases the intensity of hydrogen emission produced in the cores and wings of spectral lines (Figs. 4 and 5). Since the ambient densities and temperatures in the flaring chromosphere are also increased by the beam heating, hydrogen spectral lines become broadened.

Besides heating, beam electrons also raise the ionisation degree of hydrogen atoms below the flare transition region because the ionisation rates from inelastic collisions between beam electrons and hydrogen atoms dominate over the ionisation rates caused by thermal electrons (Fig. 6). A growth of ionisation degree caused by non-thermal electrons extends the spectral line wings, owing to the increased local electric fields (the Stark effect). Therefore, the beams with higher initial fluxes produce a larger intensity in wings and more broadened and flattened cores of Balmer and Pashen lines (Figs. 4 and 5) that confirm the previous conclusions by (Aboudarham & Henoux 1986; Zharkova & Kobylinskii 1993). Additionally, in the current study we are able to identify the hydrodynamic effects of beam electrons in the line cores (see Sect. 3.2.1).

The magnitudes of Doppler shifts in the spectral line profiles indicate the macro velocities of the shocks when they pass through the formation regions of a given spectral line. As a result, the Hα line undergoes largest red wing enhancements among the hydrogen lines presented in our simulations because this line is formed at upper chromospheric depths where the macro velocities of the shock are the largest. The lines formed at deeper atmospheric depths, such as the Hγ or Pα lines, have a smaller proportion of their core emission being red-shifted (Fig. 8b). By comparing the Hα, Hβ, and Hγ line profiles observed with a high cadence during the first minute or two after a flare (beam injection) onset, the observers can track propagation of the hydrodynamic shocks through the whole flaring chromosphere (Fig. 8). The time delay between a beam onset and the instance of hydrodynamic shocks entering the formation regions of hydrogen lines is much shorter for the beams with greater initial fluxes. This occurs because the beams with higher initial fluxes produce the shocks with larger downward macro velocities (Sect. 3.2.2).

In our F10 model the downward chromospheric shocks produce red wing enhancement in the Hα profiles with a 1 Å Doppler shift at the very first 5 s after a beam onset (our Fig. 4b). The scale of the Doppler shifts in these profiles agrees with those of Canfield & Gayley (1987) for a similar flux. However, in Canfield & Gayley (1987) the peak in emission intensity remains firmly in the blue horn of the simulated profiles. The first radiative hydrodynamic simulations by Heinzel et al. (1994) showed blue horn asymmetries due to absorption in the red wing of Hα at wavelengths around 1 Å. A more recent radiative hydrodynamic model by Allred et al. (2005) simulated for the same beam parameters as in our model F10, does not show any significant shocks moving downward in the chromosphere and, thus, no red-shifted Hα line profiles (see their Fig. 8). Moreover, the Allred et al. (2005) F11 model yields a downward velocity only of 40 km s⁻¹, while our F11 model produces downward velocities up to 200 km s⁻¹, providing a red-shift up to 3 Å in Hα line profiles,similar to those reported from observations.

In addition, apart from the chromospheric evaporation simulated for coronal temperatures detectable in soft X-rays and extra ultraviolet emission, our hydrodynamic simulations do not show any significant upward motions in the flaring chromosphere, a formation region of the Balmer, or Paschen series during a beam injection phase (Fig. 1). Thus, our models do not produce not observable blue-shifted Balmer or Paschen line emission, in contrast to the other radiative hydrodynamic simulations (Allred et al. 2005; Rubio da Costa et al. 2016).

Hα observations: During the impulsive phase of solar flares observed Hα line profiles often reveal strong red-shifts of 1–4 Å (Ichimoto & Kurokawa 1984; Wuelser & Marti 1989; Wuelser et al. 1994). These profiles cannot be explained bychanges in the maximum opacity (Kuridze et al. 2015). The HYDRO2GEN models presented in this paper explain these profiles by the Doppler effect of the plasma macro motion due to hydrodynamic shocks produced by a non-thermal electron beam (Figs. 7–9).

The HYDRO2GEN F10 model produces enhancements in Hα line profiles that are in excellent agreement with those observed with SST for the C1.5 class flare on June 30, 2013 (Druett et al. 2017); this model addresses the concerns raised by Rubio da Costa et al. (2016), while interpreting Hα emission from an X1.0-class flare. Rubio da Costa et al. (2016) suggested red asymmetries peaking around 1 Å are to be likely caused by the downflows missing from their simulations. The enhancements of the Hα emission intensity in the core and red wing observed by Rubio da Costa et al. (2016) at these times are within a factor of 2–4 of those generated via the HYDRO2GEN F10 model. Moreover, the wing enhancements of the Hα profiles presented in Rubio da Costa et al. (2016) are naturally explained by our simulated profiles, owing to the Stark effect, caused by the additional ionisation of a flaring chromosphere by beam electrons (see our Figs. 4, 7).

The maximal red-shift, somewhat under 1 Å, in their Hα line profiles is achieved for the F11 model by Allred et al. (2005) occurring 5 s after the beam onset and disappearing 6 s later, being replaced with a large blue-shift (see Fig. 9 in Allred et al. 2005). Our simulations with higher initial fluxes, however, can match the maximum Doppler shifts of 3–4 Å reported by Ichimoto & Kurokawa (1984); Wuelser & Marti (1989) at the peak of the HXR emission and lasting for minutes (Druett et al. 2017). For stronger beams HYDRO2GEN also generates stronger Hα line wings, similar to those often observed and explains that these line wings are caused by an increased number of the ambient electrons (the Stark effect) that appear owing to strong ionisation by electron beams (Sect. 3.2.1 and Fig. 4).

By extracting the relevant line intensities from our full Hα line profiles using a narrow wavelength window of about 0.25 Å and examining the observations by Kaempfer & Magun (1983); Wuelser & Marti (1989), we can derive that with such the narrow wavelength window it is not possible to capture any large red-shifted emission caused by a strong chromospheric downflow. Instead, the instrument with such the narrow spectral window would observe blue wings of the shifted Hα line profiles with very low intensity until the red-shifted Hα line returns back to its normal spectral position as reported by these observations. Since downflow velocities of the shocks can be rather high, they can produce large red-shifts of up to 4 Å. This offers a very physical alternative explanation of the delays between HXR and Hα emission reported by Kaempfer & Magun (1983). Hence, larger wavelength windows are essential to capturing the full behaviour of Hα line profiles during the impulsive phase of solar flares.

PaschenPα radiative response: The Pashen emission in a flaring chromosphere is also strongly affected by electron beams, which cause the Paschen Pα and Pβ line cores to become highly optically thick (Figs. 2c and 3b). For the beams with larger initial fluxes (F11 and F12 models) the Paschen Pα and Pβ line cores are formed in the upper flaring chromosphere, whereas for the F10 model they are formed throughout the lower chromosphere. As a result, in our simulations there is a great red wing excess occurring at the peak of beam fluxes as shown in the Pα line for the F11 and F12 models, whereas only a small proportion of the Paschen Pα core emission was red-shifted in the F10 model (Fig. 9), where the hydrodynamic shock had not yet reached the lower chromosphere.

Radiative responses of the Balmer and Paschen continua: We show that the effects of electron beams (Fig. 10) caused by non-thermal collisions between beam electrons and hydrogen atoms define the hydrogen ionisation rates. As result, the Balmer continuum contributions in the flaring chromosphere are most effectively enhanced by the energy delivered by electron beams. The ratio of the intensity enhancement in the Balmer continuum head to that in the heads of other continua is found to decrease with the increasing initial flux of a beam. This occurs because stronger beams, which produce stronger flares, deliver greater energies to deeper atmospheric depths where the other continua are formed (Tables 2 to 4).

Moreover, our simulations demonstrate that a power-law electron beam is the much more effective agent for production of Pashen continuum or WL emission than thermal electrons caused by a temperature increase during flares (Fig. 11). The non-thermal collisions with beam electrons immediately generate WL emission that begins co-temporally with HXR emission as indicated by Aboudarham & Henoux (1986). In weaker flares, the WL enhancement due to non-thermal beams is principally generated at upper chromospheric densities. The WL emission for stronger F11 and F12 beam models (derived from the Pashen continuum contribution functions, see Fig. 11) is also more effectively enhanced at upper chromospheric column depths, in addition to the normal WL emission occurring at the photopheric depths in thermal flares. This is caused by the combined radiative transfer and non-thermal ionisation effects in the Lyman continuum governing hydrogen ionisation affected by electron beams penetrating into the deeper atmospheric layers.

We do not replace the source functions for any transitions with those for black bodies. Instead we consider radiative transfer in the Lyman continuum because this radiation controls the ionization of hydrogen in a flaring atmosphere. This is an essential difference from many other approaches, such as those used in Ricchiazzi & Canfield (1983); Allred et al. (2005).

Observations of Balmer and Pashen continua and HXR: The chromospheric Balmer continuum enhancement in our models results in the close agreement with the enhancement observed using IRIS at a wavelength of 2830Å reported in Heinzel & Kleint (2014). The column depths at which the WL emission for the F11 and F12 models largely originate are at the locations of stopping depths of lower cut-off electrons (10–20 keV). This is not the case for the thermal flare models in which WL sources principally occur at photospheric column depths.

Since HXR emission is mainly produced at the chromospheric footpoints by electrons with such energies of 10–20 keV, this explains the close correlation between the heights of HXR and WL emission observed in the limb flares (Martínez Oliveros et al. 2012; Krucker et al. 2015). The heights of these chromospheric sources of WL emission in our models, resulting from the sweeping of the QS plasma downward to the photosphere, due to heating by non-thermal electron beams, are found to be close to the photospheric heights in the quiet Sun seen on the limb. This confirms that electrons beams are the most plausible agents, which can explain solar flare observations on the limb by Martínez Oliveros et al. (2012) and Krucker et al. (2015).

Summary: Therefore, non-thermal particle beams are undoubtedly shown for the first time to be the very likely agents producing very large Doppler red-shifts (1–4 Å in Hα) observed in the hydrogen spectral lines during the impulsive phase of solar flares. These large red-shifts are the signature of a high-velocity macro motion of the hydrodynamic shocks generated by such beams and moving to the QS photosphere. These beams can also account for the intensities and depths of the formation of WL emission being rather close to the depths of maximum HXR emission, as reported from many observations of limb flares.

Acknowledgements

M.D. and V.Z. wish to acknowledge the Ph.D. studentship support from Northumbria University via Research Councils UK (RCUK). The authors express our deepest gratitude for the constructive discussion and comments by Philip Judge, from which the paper strongly benefitted.

Appendix A: Non-local thermodynamic equilibrium model for hydrogen

A.1: Statistical equilibrium

Statistical equilibrium for each bound level of the model hydrogen atom, i, at each plasma column depth, ξ (cm⁻²), is constrained as follows: $$\begin{array}{lll}\frac{\partial n_i(\xi )}{\partial t}\hfill & =-n_i(\xi )\left({\displaystyle \sum _{k=1,k\ne i}^N{R}_{ik}}(\xi )+R_{ic}\right)\hfill & +{\displaystyle \sum _{k=1,k\ne i}^N{n}_k}(\xi )R_{ki}(\xi )+n_e{n}^{+}(\xi )R_{ci}(\xi )=0,\hfill \end{array}$$(A.1)

and in the continuum, $$\frac{\partial n_{+}(\xi )}{\partial t}=-n_e{n}_{+}(\xi ){\displaystyle \sum _{i=1}^N{R}_{ci}}(\xi )+{\displaystyle \sum _{i=1}^N{n}_i}(\xi )R_{ic}(\xi ),$$(A.2)

where N = 5 for a 5 level plus continuum model atom and R_ik is the net probability of an atom moving from state i to state k. The recombination terms R_ci are multiplied by the product of the free electron and ionised hydrogen number densities n_e n⁺(ξ). This is because both constituents are required for the recombination to occur and therefore the recombination rate coefficients R_ci are defined below in units cm⁻⁶ s⁻¹, whereas the other rate coefficients R_ik and R_ic are defined in units cm⁻³ s⁻¹. In the case of bound levels (Eq. (A.1)) the first summation term represents deactivation, out of level i, to other levels. The second summation term represents activation, into level i, from other levels. In the continuum case (Eq. (A.2)) these terms represent recombination and ionisation, respectively. The rate coefficients are calculated using the formulas below: $$\begin{array}{l}k<i,\\ R_{ik}(\xi )=n_e(\xi )C_{ik}(\xi )+A_{ik}\\ k>i,\\ R_{ik}(\xi )=n_e(\xi )C_{ik}(\xi )+n_e^b(\xi )C_{ik}^b(\xi )+Q_{ik}(\xi )+D_{ik}(\xi )\end{array}$$(A.3) $$\begin{array}{l}k>i,\\ R_{ki}(\xi )=n_e(\xi )C_{ki}(\xi )+A_{ki}\\ k<i,\\ R_{ki}(\xi )=n_e(\xi )C_{ki}(\xi )+n_e^b(\xi )C_{ki}^b(\xi )+Q_{ki}(\xi )+D_{ki}(\xi )\end{array}$$(A.4) $$\begin{array}{l}\text{ionisation},\\ R_{ic}(\xi )=n_e(\xi )C_{ic}(\xi )+n_e^b(\xi )C_{ic}^b(\xi )+Q_{ic}(\xi )+D_{ic}(\xi )\end{array}$$ $$\begin{array}{l}\text{recombination},\\ R_{ci}(\xi )=A_{ci}(\xi )+n_e(\xi )C_{ci}(\xi ).\end{array}$$

Excitation and de-excitation rates for collisions between hydrogen atoms and thermal electrons, C_ik, C_ki, are taken from Johnson (1972). For excitation, de-excitation, and ionisation the rates per second, n_i C_ik and n_i C_ic are formed by multiplying these coefficients by the abundances in the initial state, i, of each process. In three-body recombination electrons transfer energy and momentum to another electron in the vicinity of a hydrogen ion, resulting in recombination to a bound level, i. Coefficients for this process are taken from Johnson (1972) and multiplied by the product of abundances of hydrogen ions and free electrons, n_e n⁺, and again by the abundance of free electrons, n_e. This process produces lower rates per second than others considered in all conditions used in our simulations. For excitation and ionisation due to collisions between hydrogen atoms and the beam electrons $C_{ik}^b$, $C_{ic}^b$, we use the analytical formulae devised by Zharkova & Kobylinskii (1993) with hydrogen photoionisation cross section formulae from Beynon (1965); Burgess (1965); Ludde & Dreizler (1982). The depth distribution of the beam electron number densities, $n_e^b$, is modelled using the work of Zharkova & Kobylinskii (1993). To generate rates per second, the coefficients are multiplied by the abundance of hydrogen in the initial state for a transition to state k or ionisation to continuum, c, $n_e^b{C}_{ik}^b$, $n_e^b{C}_{ic}^b$.

The hydrogen spontaneous radiative rates (Allen 1977) are calculated by multiplying the coefficient of single interaction, A_ji, by the number density of hydrogen atoms in the upper state n_j of the transition from j to i, n_j A_ji. Recombination rates are derived by multiplying the coefficient by the product of the abundances of the hydrogen ions and free thermal electrons required for recombination to occur, n_en⁺, to level i, A_ci n_en⁺. Stimulated excitation or de-excitation, Q_ik, and ionisation, Q_ic, coefficients in the model due to the external radiation (sometimes called “backwarming”) entering from the levels above and below are taken from Zharkova & Kobylinskii (1991) as follows: $$\begin{array}{ll}{Q}_{ik}^{\odot }(\xi )=B_{ik}{j}_{ik}^{\odot }(\xi )\hfill & =\frac{g_k}{g_i}\frac{c^2{A}_{ki}}{2h{\nu }_{ik}^3}\frac{W_{ik}}{2\sqrt{\pi }}\hfill \\ \hfill & \times {\displaystyle {\int }_{\infty }^{\infty }{I}_{ik}^{\odot }}\alpha (x)E_2(\alpha (x){\tau }_{ik}(\xi ))\text{d}x\hfill \end{array}$$(A.9) $$\begin{array}{ll}{Q}_{ic}^{\odot }(\xi )=B_{ic}{j}_{ic}^{\odot }(\xi )\hfill & =4\pi {\kappa }_{ic}^0\frac{W_{ic}}{2h}\hfill \\ \hfill & \times {\displaystyle {\int }_{{\nu }_{ic}}^{\infty }{I}_{ic}^{\odot }}{f}_i(\nu )E_2(f_i(\nu ){\tau }_{ic}(\xi ))\frac{\text{d}\nu }{\nu }\cdot \hfill \end{array}$$(A.10)

They take the forms $Q_{ik}=B_{ik}{j}_{ik}^{\odot }$, $Q_{ki}=B_{ki}{j}_{ik}^{\odot }$, and $Q_{ic}=B_{ic}{j}_{ic}^{\odot }$ where the coefficients B_nm are the Einstein coefficients for stimulated emission and absorption, $j_{nm}^{\odot }$ is the average intensity of radiation at depth ξ from sources entering from the top and bottom of the model for the transition n → m, and W_ik and W_ic are the dilution factors of the lines and continua, respectively. These are multiplied by the abundance in the initial states of the processes to produce excitation and ionisation rates, n_iQ_ik and n_iQ_ic.

The values D_ik and D_ic are the coefficients for transitions and ionisation due to internal diffusive radiation. They take the forms $D_{ik}=B_{ik}{j}_{ik}^{\partial }$, $D_{ki}=B_{ki}{j}_{ik}^{\partial }$, and $D_{ic}=B_{ic}{j}_{ic}^{\partial }$, where $j_{nm}^{\partial }$ is the average emission intensity of diffusive radiation within the model for the transition n → m. These rates are incorporated by solving radiative transfer equations and converting the steady-state algebraic equations into integral radiative transfer equations of Fredholm of the second kind. The Fredholm integral equations are solved using the L2 approximation method described in Sect. 2.2.3.

A.2: Particle conservation

Since the flaring plasma remains neutral, the sum of the number densities of hydrogen ions and atoms, n(ξ), at a column depth ξ is conserved during the radiative processes leading to the following equation to hold: $$n(\xi )={\displaystyle \sum _{i=1}^N{n}_i}(\xi )+n^{+}(\xi ).$$(A.11)

A.3: Radiative transfer

The radiative transfer equation for an intensity of light with frequency ν, travelling towards the corona is $$\mu \frac{\partial I(\nu ,\tau )}{\partial \tau }=\alpha (\nu )I(\nu ,\tau )-\alpha (\nu )S(\tau ),$$(A.12)

where τ is the optical depth (see Appendix C), α(ν) is the absorption coefficient at frequency ν, and the source function S is the coefficient of emissivity divided by that for absorption. In the hydrogen transition lines this becomes $$\mu \frac{\partial I_{ij}(\nu ,\tau )}{\partial \tau }={\alpha }_{ij}(\nu )I_{ij}(\nu ,\tau )-{\alpha }_{ij}(\nu )S_{ij}(\tau ).$$(A.13)

Here α_ij is the absorption profile of radiation from a spectral line for the transition from j → i, assumed to be Voigtian (Eq. (D.1) in Appendix D), μ is the cosine of a propagation angle of radiation with respect to the local vertical in the flaring atmosphere, I_ij is the intensity in the spectral line for the transition j → i. The value τ is the optical depth in the relevant line centre, ν = c∕λ is the frequency, and λ a wavelength of line radiation. The source functions S_ij(τ) are related to the relative populations of the atomic levels for the bound-bound transitions (Eq. (A.14)) $$S_{ij}(\nu )=\frac{2h{\nu }_{ij}^3}{c^2}{\left(\frac{n_i{g}_j}{n_j{g}_i}-1\right)}^{-1}.$$(A.14)

In the continua the source functions are related to the relative emission measure derived from Eq. (8). The intensity of radiation travelling towards the corona in the ith continuum of hydrogen I_ic is derived from Sect. 7.5 in Mihalas (1978) $$\begin{array}{ll}\mu \frac{\partial I_{ic}(\nu ,\tau )}{\partial \tau }\hfill & =f_i(\nu )I_{ic}(\nu ,\tau )\hfill \\ \hfill & -\frac{{ϵ}_{ic}({\nu }_{ic},\tau )}{n_i(\tau ){\kappa }_{ic}(\tau )}\exp \left(-\frac{h(\nu -{\nu }_{ic})}{k_{B}T(\tau )}\right)f_i(\nu ){\nu }^3.\hfill \end{array}$$(A.15)

With emissivity ϵ_ic(ν_ic, τ^′), $${ϵ}_{ic}({\nu }_{ic},{\tau }^{\prime })=\frac{n_e{n}_{+}({\tau }^{\prime })A_{ci}({\tau }^{\prime })h}{4\pi {\displaystyle \underset{{\nu }_{ic}}{\overset{\infty }{\int }}{f}_i}{\nu }^2\exp \left(-\frac{h(\nu -{\nu }_{ic})}{k_{B}T({\tau }^{\prime })}\right)\text{d}\nu },$$(A.16)

where τ is the optical depth in the continuum head, c is speed of light, and ν_ic is the frequency in the head of the ith continuum. The value κ_ic (τ) is the absorption coefficient in the ith continuum head. The value h is Planck’s constant, k_B is Boltzmann’s constant, and T(τ) is the kinetic temperature of the ambient plasma. The value f_i(ν) is the absorption profile of the ith continuum, the form of which is taken from the work of Morozhenko & Zharkova (1980, 1982); Morozhenko (1983, 1984), and was discussed by Canfield & Athay (1974); Canfield & Puetter (1981); Neidig & Wiborg (1984) as follows: $$f_i=\{\begin{array}{ll}\frac{{\lambda }^3}{{\lambda }_{ic}^3}=\frac{{\nu }_{ic}^3}{{\nu }^3}\hfill & i=1,2,3(Ly,Ba,Pa)\hfill \\ \frac{{\lambda }^2}{{\lambda }_{ic}^2}=\frac{{\nu }_{ic}^2}{{\nu }^2}\hfill & i=4,5(Br,Pf).\hfill \end{array}$$(A.17)

The lines of Lyman series and Lyman continuum are optically thick in all of the models presented. The α lines of the Balmer, Paschen, and for some models even the Bracket series become optically thick. Moreover, the Paschen beta (Pβ) spectral line becomes optically thick if affected by powerful beams (model F12) in strong flares. The formal solutions of radiative transfer equations for upward travelling intensity are $$I(\nu ,{\tau }_{\nu })={\displaystyle \underset{{\tau }_{\nu }}{\overset{\infty }{\int }}\frac{\alpha }{\mu }}{e}^{-\frac{\alpha }{\mu }(t_{\nu }-{\tau }_{\nu })}{S}_{\nu }\text{d}{t}_{\nu }$$(A.18)

and in the downward travelling intensity $$I(\nu ,{\tau }_{\nu })={\displaystyle \underset{0}{\overset{{\tau }_{\nu }}{\int }}\frac{\alpha }{\mu }}{e}^{-\frac{\alpha }{\mu }({\tau }_{\nu }-t_{\nu })}{S}_{\nu }\text{d}{t}_{\nu }.$$(A.19)

Assuming complete redistribution of frequencies, the mean intensities give the expressions for the simulated rates of excitation, de-excitation, and ionisation due to internal diffusive radiation, as stated in Sect. 2.2.2 and Eqs. (11) and (12) of Zharkova & Kobylinskii (1991).

Appendix B: Intensities of emission in the lines and continua

The formulafor the intensity of continuum emission escaping from the top of the model is derived below. We begin from the radiative transfer equation in the continuum for radiation travelling towards the corona (see Eq. (A.15)). This is solved usingthe integrator factor method and we arrive at the formulation, considering contributions from all depth points and angles, as follows: $$\begin{array}{ll}{I}_{ic}(0,\nu )\hfill & =\frac{h}{2}{\displaystyle \underset{0}{\overset{{\tau }_{\text{max}}}{\int }}[\frac{A_{ci}({\tau }^{\prime })}{{\kappa }_{ic}}\frac{n_e{n}_{+}}{n_i}({\tau }^{\prime })}\hfill \\ \hfill & \times \frac{\exp \left(-\frac{h(\nu -{\nu }_{ic})}{k_{B}T({\tau }^{\prime })}\right)\frac{f_i}{\mu }{\nu }^3\exp \left(-\frac{f_i}{\mu }{\tau }^{\prime }\right)}{{\displaystyle \underset{{\nu }_{ic}}{\overset{\infty }{\int }}{f}_i}{\nu }^2\exp \left(-\frac{h(\nu -{\nu }_{ic})}{k_{B}T({\tau }^{\prime })}\right)\text{d}\nu }]\text{d}{\tau }^{\prime }.\hfill \end{array}$$(B.1)

The normalisation term $I_N={\left[{\displaystyle \underset{{\nu }_{ic}}{\overset{\infty }{\int }}{f}_i}{\nu }^2\exp \left(-\frac{h(\nu -{\nu }_{ic})}{k_{B}T({\tau }^{\prime })}\right)\text{d}\nu \right]}^{-1}$ depends on the form of the absorption coefficient for the continuum f_i (Eq. (A.17)), and takes the form $$\begin{array}{l}\text{for}\hspace{0.17em}i=1,2,3\\ {\left[{\nu }_{ic}^3\exp \left(\frac{h{\nu }_{ic}}{k_{B}T({\tau }^{\prime })}\right)E_1\left(\frac{h{\nu }_{ic}}{k_{B}T({\tau }^{\prime })}\right)\right]}^{-1}\\ \text{for}\hspace{0.17em}i=4,5\\ {\left[\frac{k_{B}T({\tau }^{\prime })}{h}{\nu }_{ic}^2\right]}^{-1}\cdot \end{array}$$(B.2)

The resulting formula is given in wavelength units in Sect. 2.3.2 and included below (Eq. (B.3)) for reference, with the explicit τ dependence is removed, for compactness. The relationship $\text{d}\nu =-\frac{c}{{\lambda }^2}\text{d}\lambda$ is used to convert these formula for use with the results in wavelength units written as $$\begin{array}{l}\text{for}\hspace{0.17em}i=1,2,3\\ I_{ic}(\lambda )=\frac{h}{2}\frac{c}{{\lambda }^2}\\ \times {\displaystyle \underset{0}{\overset{{\tau }_{\text{max}}}{\int }}\frac{A_{ci}}{{\kappa }_{ic}}}\frac{n_e{n}_{+}}{n_i}\frac{{\lambda }_{ic}^3}{{\lambda }^3}\frac{\exp \left(-\frac{hc}{k_{B}T\lambda }\right)\frac{f_i}{\mu }\exp \left(-\frac{f_i}{\mu }{\tau }^{\prime }\right)}{E_1\left(\frac{hc}{k_{B}T{\lambda }_{ic}}\right)}\text{d}{\tau }^{\prime }\\ \text{for}\hspace{0.17em}i=4,5\\ I_{ic}(\lambda )=\frac{h}{2}\frac{c}{{\lambda }^2}\\ \times {\displaystyle \underset{0}{\overset{{\tau }_{\text{max}}}{\int }}\frac{A_{ci}}{{\kappa }_{ic}}}\frac{n_e{n}_{+}}{n_i}\frac{{\lambda }_{ic}^2}{{\lambda }^2}\frac{\exp \left(-\frac{hc\left(\frac{1}{\lambda }-\frac{1}{{\lambda }_{ic}}\right)}{k_{B}T}\right)\frac{f_i}{\mu }\exp \left(-\frac{f_i}{\mu }{\tau }^{\prime }\right)}{\frac{k_{B}T\lambda }{hc}}\text{d}{\tau }^{\prime }.\end{array}$$(B.3)

Appendix C: Optical depth

The optical depth τ of a particular wavelength of light, λ, is the product of the absorption coefficient at that wavelength, κ and the number density of the absorbing medium through which it is passing, n_λ multiplied by the distance through which it passes, ΔS. This leads to the well-known formula for optical depth ${\tau }_{\lambda }(s)={\displaystyle \underset{s_0}{\overset{s}{\int }}{\kappa }_{\lambda }}{n}_{\lambda }\text{d}{s}^{\prime }$. The value s₀ is the height at the top of the model and s that at the point under inspection

Hydrogen in excitation state i is able to absorb light produced by the atomic transition, j → i with lower state i. Modelling hydrogen in the lower state of a line or continuum transition as the only absorbing material for photons produced by that transition, then n_λ = n_i. Thus we have, for a 5-level plus continuum hydrogen atom in a stationary medium, in the line central wavelength λ₀, $${\tau }_{ij}(\xi )={\displaystyle \underset{{\xi }_{\text{min}}}{\overset{\xi }{\int }}\frac{{\kappa }_{ij}({\lambda }_0)(1-\chi )\frac{n}{n_0}}{{\displaystyle \sum _{j=1}^5\frac{n_j}{n_i}}}}\text{d}{\xi }^{\prime },$$(C.1)

where τ_ij is the line centre optical depth or continuum head optical depth in the case j = c, ξ_min is the column depth at the top of the model and ξ that at the point in question, n is the total number density of hydrogen atoms and protons, and χ is the ionisation degree of hydrogen atoms.

Appendix D: Voigt profile

The Voigt profile, used to describe the shape of a broadened spectral line is given in Ivanov & Serbin (1984) as $$U(a,x)=\frac{a}{{\pi }^{3/2}}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}\frac{e^{-y^2}}{{(x-y)}^2+a^2}}\text{d}y.$$(D.1)

In our context x is a dimensionless wavelength measurement that is measured in Doppler half widths. The value a is a broadening term for Doppler and Lorentzian broadening, the ratio of natural line width to Doppler width in frequency, in our code, is written as $$a=\frac{\Delta {\nu }_{\text{natural}}}{\Delta {\nu }_{\text{Doppler}}}.$$(D.2)

An un-normalised form of this function was analysed by Harris (1948) using a Taylor Series $$H(a,x)=\frac{a}{\pi }{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}\frac{e^{-y^2}}{{(x-y)}^2+a^2}}\text{d}y.$$(D.3)

Harris gave the form $$\begin{array}{ll}H(a,x)\hfill & =H_0(x)+a{H}_1(x)+a^2{H}_2(x)+a^3{H}_3(x)\hfill \\ \hfill & +\mathrm{...}+a^n{H}_n(x)+\mathrm{...},\hfill \end{array}$$

with $$\begin{array}{l}{H}_0(x)=e^{-x^2}\hfill \\ H_1(x)=-\frac{2}{\sqrt{\pi }}(1-2xF(x))\hfill \\ H_2(x)=(1-2x^2)e^{-x^2}\hfill \\ H_3(x)=-\frac{2}{\sqrt{\pi }}[\frac{2}{3}(1-x^2)-2x(1-\frac{2}{3}{x}^2)F(x)]\hfill \\ H_4(x)=\left(\frac{1}{2}-2x^2+\frac{2}{3}{x}^4\right)e^{-x^2},\hfill \end{array}$$(D.6)

where F(x) is the Dawson function $$F(x)=e^{-x^2}{\displaystyle \underset{0}{\overset{x}{\int }}{e}^{-t^2}}\text{d}t.$$(D.7)

D.1: Dawson function

In our code the modified Dawson function, MD(x) = 1 − 2xF(x), is calculated using procedures for the Dawson function outlined in Cody et al. (1970). For the domain |x| ≤ 2.5 we use the n = 5 form. For the 2.5 < |x|≤ 3.5 scheme we use the n = 4 form. We use n = 3 for cases 3.5 < |x|≤ 5 and n = 1 for 5 < |x|.

D.2: Taylor expansion coefficients

We now calculate of the Taylor expansion terms. Mihalas (1978) (p. 280, 9–40) gives that $$H_n(x)\equiv \frac{{(-1)}^n}{\sqrt{\pi }(n!)}{\displaystyle \underset{0}{\overset{\infty }{\int }}{y}^n}{e}^{-\frac{y^2}{4}}\cos (xy)\text{d}y.$$(D.8)

Let us introduce the notation $$\begin{array}{l}{C}_n(x)\equiv {\displaystyle \underset{0}{\overset{\infty }{\int }}{y}^n}{e}^{-\frac{y^2}{4}}\cos (xy)\text{d}y\hfill \\ S_n(x)\equiv {\displaystyle \underset{0}{\overset{\infty }{\int }}{y}^n}{e}^{-\frac{y^2}{4}}\sin (xy)\text{d}y.\hfill \end{array}$$(D.9)

This gives $$\begin{array}{ll}{C}_0=\sqrt{\pi }{e}^{-x^2}\hfill & S_0=2F(x).\hfill \end{array}$$(D.10)

Integrating by parts gives $$\begin{array}{l}{C}_1=2(1-2xF(x))\\ S_1=2x\sqrt{\pi }{e}^{-x^2}\\ C_n=2(n-1)C_{n-2}-2x{S}_{n-1}\\ S_n=2(n-1)S_{n-2}+2x{C}_{n-1}.\end{array}$$(D.11)

It would be preferable to create equivalent recurrence formulae for the Taylor expansion coefficients themselves, where $$H_n(x)\equiv \frac{{(-1)}^n}{\sqrt{\pi }(n!)}{\displaystyle \underset{0}{\overset{\infty }{\int }}{y}^n}{e}^{-\frac{y^2}{4}}\cos (xy)\text{d}y\equiv \frac{{(-1)}^n}{\sqrt{\pi }(n!)}{C}_n(x).$$(D.12)

Let us also define $$I_n(x)\equiv \frac{{(-1)}^n}{\sqrt{\pi }(n!)}{\displaystyle \underset{0}{\overset{\infty }{\int }}{y}^n}{e}^{-\frac{y^2}{4}}\sin (xy)\text{d}y\equiv \frac{{(-1)}^n}{\sqrt{\pi }(n!)}{S}_n(x).$$(D.13)

Then we also see that $$\begin{array}{ll}{S}_n={(-1)}^n\sqrt{\pi }(n!)I_n\hfill & C_n={(-1)}^n\sqrt{\pi }(n!)H_n.\hfill \end{array}$$(D.14)

And, by substitution $$\begin{array}{l}{H}_0=\frac{{(-1)}^0}{\sqrt{\pi }(0!)}\sqrt{\pi }{e}^{-x^2}=e^{-x^2}\\ I_0=\frac{{(-1)}^0}{\sqrt{\pi }(0!)}2F(x)=\frac{2}{\sqrt{\pi }}F(x)\\ H_1=\frac{{(-1)}^1}{\sqrt{\pi }(1!)}2(1-2xF(x))=\frac{2}{\sqrt{\pi }}(2xF(x)-1)\\ I_1=-2x{e}^{-x^2}\\ H_n=\frac{2H_{n-2}+2x{I}_{n-1}}{n}\\ I_n=\frac{2I_{n-2}-2x{H}_{n-1}}{n}\cdot \end{array}$$(D.15)

The coefficients in the Taylor expansion of the line broadening function can be calculated using these recurrence relations. The area under the curve is normalised to $\sqrt{\pi }$ so when the normalisation required is 1, division of the result by $\sqrt{\pi }$ is necessary. The advantages of this scheme over the differentiation based recurrence method discussed in Mihalas (1978, exercise 9.2) are that the calculations do not require retention of the explicit, ever-expanding, polynomial forms of factors in the coefficients. The additional calculation cost of terms is linear beyond the second term. Values from the process are compared against tables in Finn & Mugglestone (1965) and are seen to agree to the five significant figures used in saving values from our Fortran process: for a values 0 to 1 steps of 0.1 and for x values from 0 to 6 in steps of 0.5 and to a slightly lower, but tolerable accuracy for x from 6.5 to 15.

References

All Tables

All Figures

Fig. 1 Hydrodynamic responses of the ambient plasma to electron beam heating for the F10 (red), F11 (green), and F12 (blue) models for times t = 1–5 s. These are shown from the top row to the bottom row, respectively, as follows: the left panels indicate electron kinetic temperature; central panels indicate plasma density; right panels indicate plasma macro velocity plotted vs. plasma column depth (cm−2 ).

In the text

Fig. 2 Optical depths calculated in the differential study for (a) the Lyman alpha line; (b) the Hα line; and (c) the Pβ line 5 s afterthe onset of the beam, plotted against column depth (cm−2). Each subfigure shows the optical depth in the central wavelength of the line profile calculated for the thermal response (black line)and for models including excitation and ionisation rates from electron beams with fluxes in the range 109 − 1012 erg cm−2 s−1.

In the text

	Fig. 3 Optical depths in coupled hydrodynamic and radiative models, 5 s after beam onset in a) the Lyman alpha line and b) the Pβ line, plotted against plasma column depth (cm−2). Simulations are made with beams with fluxes of 1010 (red line), 1011 (green line), and 1012 (blue line) erg cm−2 s−1.
In the text

	Fig. 4 Profiles of Hα line intensityenhancements, in erg cm−2 s−1 sr−1 Å−1, plotted vs. distance in Å from the line central wavelength for times from 1 to 5 s in panels (a) to (e), respectively.Results are shown for the F12 (blue line), F11 (green line), and F10 (red line) beam models. The colour scheme for the lines representing these models is also used in Figs. 5 to 10.
In the text

	Fig. 5 Intensity enhancements in erg cm−2 s−1 sr−1 Å−1 against distance from line centre in Å for the Hα (top panels, a and b), Balmer beta (Hβ) (middle panels, c and d), and Hγ (bottom panels, e and f) lines 4 s (left panels) and 5 s (right panels) after the beam injection begins.
In the text

	Fig. 6 Ionisation degree at 4 (a) and 5 (b) seconds plotted against plasma column depth (cm−2) for the thermal (black line), F10 (red line), F11 (green line), and F12 (blue line) models.
In the text

	Fig. 7 Hα line intensity enhancements with macro-velocity effects included for the F10 (red), F11 (green), and F12 (blue) models. Profiles are shown for a) t = 2 s; b) t = 3 s; c) t = 4 s; and d) t = 5 s in erg cm−2 s−1 sr−1 Å−1 against distance from line centre in Å.
In the text

	Fig. 8 Panel a: spectral profiles of the Hγ line for t = 2 s (light green) to t = 5 s (black). Panel b: Hα, Hβ, and Hγ line profiles at t = 5 for the F10 model. Intensities (I) are simulated in erg cm−2 s−1 sr−1 Å−1 against distance from line centre (λ − λ0) in Å.
In the text

	Fig. 9 Pα line profiles vs. distance in Å from the line centre wavelength calculated at a) t = 2; b) t = 3; c) t = 4; and d) t = 5, in erg cm−2 s−1 sr−1 Å−1.
In the text

	Fig. 10 Intensities of Balmer continuum enhancement at a) t = 1; c) t = 3; e) t = 5, and Paschen continuum enhancement at b) t = 1; d) t = 3; f) t = 5, in erg cm−2 s−1 sr−1 Å−1 against wavelength in Å.
In the text

Fig. 11 Paschen continuum contribution functions (see Sect. 2.3.3) calculated at t = 5 s. The z-axis shows the contribution function value, the x- and y-axes are in terms of wavelength, λ and the logarithm of column depth log10ξ. The top row shows results for the F10 flare model, the central row shows data for the F11 flare model, and the bottom row shows data for the F12 model. The left panels show the models including excitation and ionisation of hydrogen atoms by non-thermal beam electrons and the right panels represents the “thermal flare” simulations without excitation and ionisation by beam electrons.

In the text

Fig. 12 Panel a: simulated Hα profiles showing red wing enhancement for the F10 (red line), F11 (green line), and F12 (blue line) models. Normalised intensity enhancement is on the y-axis, and wavelength relative to the Hα profile line centre on the x-axis. Panel b: excess of emission in the red wing of the profiles in panel “a” relative to the emission in the blue wing. Panel c: normalised Hα profiles from the M1 solar flare on May 24, 1987 reported in Wuelser & Marti (1989). The normalised Hα profile (green line), quiet Sun background level (grey line), and enhancement in emission (pink line) are shown for the A1 kernel at the time of the maximum of the HXR (15:24:55 UT). The y-axis shows normalised intensity and the wavelength relative to the Hα line central wavelength is shown on the x-axis. Panel d: excess emission in the red wing compared to the blue wing for the A1 kernel at 15:24:55 UT (green line). The red line shows a simulated fit of this excess based on a linear combination of contributions based on the F10 and F11 models.

In the text