Role of Suprathermal Runaway Electrons Returning to the Acceleration Region in Solar Flares

When a sub-Dreicer electric field (${\mathscr{E}}\lt {{\mathscr{E}}}_{D}$) is applied to a thermal plasma, a fraction of the Maxwellian distribution will run away, and the bulk of the Maxwellian distribution will drift under the influence of the electric field. A higher electric field results in a higher runaway rate. Many authors carried out solutions to the spatially homogeneous Fokker–Planck equation where an electric field is applied to a Maxwellian distribution (e.g., Knoepfel & Spong 1979, for a review) and calculated the runaway growth rate.

In this section we benchmark the numerical growth rate from Landreman et al. (2014) and Kulsrud et al. (1973) with the analytical expression of Kruskal & Bernstein (1964, hereafter KB64). We must first adjust these growth rates to consistent definitions of the Dreicer field, collision frequency, and thermal velocity.

KB64 provided an asymptotic solution to the spatially homogeneous Fokker–Planck equation by expanding the distribution function in powers of $\tfrac{{\mathscr{E}}}{{{\mathscr{E}}}_{D}}$ near $\mu =\tfrac{{v}_{\parallel }}{v}=1$ by recognizing that if the electric field is much less than the Dreicer field, then the critical velocity v_cr above which electrons are in the runaway regime have a pitch angle μ ∼ 1. Under these conditions, the critical velocity above which runaway occurs is given by (e.g., Smith & Verwichte 2008)

$$\begin{eqnarray}&&{v}_{{cr}}={v}_{e}\sqrt{{{\mathscr{E}}}_{D}/2{\mathscr{E}}},\end{eqnarray} \tag{ 1 }$$

where ${v}_{e}=\sqrt{2{k}_{B}T/{m}_{e}}$ is the thermal velocity. The Kruskal & Bernstein runaway rate, in units per second, is given in (e.g., KB64 ; Holman 1985),

$$\begin{eqnarray}&&{\gamma }_{\mathrm{run}}=C{\nu }_{e}{\left(\displaystyle \frac{{{\mathscr{E}}}_{D}}{{\mathscr{E}}}\right)}^{3/8}\exp \left(-\sqrt{2\displaystyle \frac{{{\mathscr{E}}}_{D}}{{\mathscr{E}}}}-\displaystyle \frac{1}{4}\displaystyle \frac{{{\mathscr{E}}}_{D}}{{\mathscr{E}}}\right),\end{eqnarray} \tag{ 2 }$$

where C is a constant of order 1 (e.g., Hesslow et al. 2019; Connor & Hastie 1975), and

$$\begin{eqnarray}&&\begin{array}{l}{{\mathscr{E}}}_{D}=4\pi {e}^{3}\mathrm{ln}\,{\rm{\Lambda }}\displaystyle \frac{{n}_{e}}{{k}_{B}T}=\displaystyle \frac{3\sqrt{\pi }}{2}\displaystyle \frac{{{mv}}_{e}{\nu }_{e}}{e}\\ \quad =1.0\times {10}^{-11}\mathrm{ln}\,{\rm{\Lambda }}\displaystyle \frac{{n}_{e}[{\mathrm{cm}}^{-3}]}{T[K]}\end{array}\end{eqnarray} \tag{ 3 }$$

is the Dreicer field in StatV per centimeters (e.g., Somov 2012),

$$\begin{eqnarray}&&{\nu }_{e}=\displaystyle \frac{4}{3}\displaystyle \frac{\sqrt{2\pi }}{\sqrt{{m}_{e}}}\displaystyle \frac{{{ne}}^{4}\mathrm{ln}\,{\rm{\Lambda }}}{{\left({{kT}}_{e}\right)}^{3/2}}\end{eqnarray} \tag{ 4 }$$

is the electron collision frequency given by Huba (2018). In these expressions, m_e is the electron mass, e is the absolute value of the electron charge, $\mathrm{ln}{\rm{\Lambda }}$ is the Coulomb logarithm, n_e is the electron density, k_B is the Boltzmann constant, and T_e is the electron temperature.

Kulsrud et al. (1973) numerically solved the Fokker–Planck equation using a collision operator that is linearized, nonrelativistic, spatially homogeneous, and axially symmetrical in velocity space. The solutions are represented by the blue crosses in Figure 1, and the collision frequency used is ${\nu }_{k73}=\tfrac{3\sqrt{\pi }}{\sqrt{2}}{\nu }_{e}$. Landreman et al. (2014) used the relativistic test-particle collision operator in the limit of nonrelativistic plasma temperature (Papp et al. 2011). The collision operator is valid for arbitrary electron velocities. Their operator neglects self-collisions between runaway electrons, i.e., runaways only collide with the Maxwellian bulk. The code can handle any initial distribution, not only small deviations from a Maxwellian.

The runaway growth rate is given by Equation (17) in Landreman et al. (2014). Their results have been extended to higher fractions of the electric field from 0.1–0.7 in Figure 1 using Collisional Distributions of Electrons (Landreman et al. 2014; Stahl et al. 2015; Embréus et al. 2018).

To sum up, the KB64 solution is equivalent to Landreman's at lower fractions of the electric field: ${{\mathscr{E}}}_{\mathrm{RC}}/{{\mathscr{E}}}_{D}\lesssim 0.1$, if the constant C is equal to 1 instead of 0.35 (e.g., Holman 1985), and the Dreicer field and collision frequency are defined by Equations (3) and (4), respectively.

In this section we derive a set of equations for describing the beam/RC system in which a power-law distribution of nonthermal electrons is injected into a thermal, fully ionized hydrogen plasma modeled by a 1D loop. The assumptions used made by our model are described in the following:

• 1.  
  The accelerated electrons are continuously injected at the apex of a 1D loop model. Their flux density distribution is represented by a single power law with a sharp low-energy cutoff E_c0. That is the flux density distribution, F, is given by

  $$\begin{eqnarray}&&F(E,x=0)=(\delta -1){F}_{e0}{E}_{c0}^{\delta -1}{E}^{-\delta },\end{eqnarray} \tag{ 5 }$$

  where δ is the electron spectral index, F_e0 the total injected electron flux density and E is the electron energy, and x is the geometric distance from the loop apex.
• 2.  
  Steady state: We consider timescales much longer than that determined from the electron–ion collision frequency meaning that RC/beam system has reached a steady state (van den Oord 1990). This is given by the inverse of Equation (4). This is much less than 4 s for usual values of coronal temperature and density combinations in solar flares. For example, for a 30 MK and 10⁹ cm⁻³ plasma the collision time is 2.1 s, and 15 ms for 1 MK and 10⁹ cm⁻³ plasma. Siversky & Zharkova (2009) showed that the time to reach the steady state is between 0.07 and 0.1 s by solving a time-dependent Fokker–Planck equation. This is consistent with the collisional timescales. We are in addition limited by the 4 s time resolution of Ramaty High Energy Solar Spectroscopic Imager (RHESSI; Lin et al. 2002) for a statistically significant spectrum.
• 3.  
  1D loop: All electrons have velocities parallel to the magnetic field. This is useful to provide direct comparisons with previous results of Holman (2012) and AH17, and gain physical insight into the problem. We will assume the temperature and density to be constant in the corona in Section 3 and a model with temperature and density stratification in Section 5.1.
• 4.  
  We focus on modeling the effect of RC and collisional losses on the dynamics of the beam/RC system. Other mechanisms not included in this paper which affect the propagation of the nonthermal beam include magnetic mirroring (Karlicky & Henoux 1993; Zharkova et al. 1995), and wave-particle instabilities (e.g., Hannah et al. 2009; Hannah & Kontar 2011).
• 5.  
  The plasma resistivity, η, is assumed to be the classical Spitzer resistivity, η_S (Spitzer 1962; Spitzer & Härm 1953).
• 6.  
  We assume that beam electrons are thermalized and lost from the beam when they have energy less than E_th = δ kT, where T is the ambient plasma temperature (Kontar et al. 2015; Jeffrey et al. 2019). Allred et al. (2020) showed that energy diffusion dominates over the dynamic friction, i.e., thermalization dominates over energy loss for an electron energy below (δ + 1)kT.
• 7.  
  We require the RC electric field to be weak enough compared to the Dreicer field everywhere along the 1D loop to accelerate less than 0.3 of the local electron density. That is we require the rate of runaways to be $\tfrac{1}{{n}_{e}}\tfrac{{{dn}}_{r}}{{dt}}={\gamma }_{\mathrm{run}}({\nu }_{e})\lt 0.3$.
• 8.  
  We assume that there is not a neutralizing ion beam propagating co-spatially along with the electron beam. Since ions are three orders of magnitude more massive, their velocity, and therefore, current density is likely much smaller than the electron beam's. Therefore, even if considered, their neutralizing effect on the electron beam would be small (see reviews by Miller et al. 1997; Holman 2016).

In a low-density corona, collisions between beam electrons and the background plasma are relatively infrequent, but the RC electric field can produce a potential drop that is a significant fraction of the beam electrons' energies. Therefore, we assume that in the corona the energy lost by beam electrons due to the RC electric field dominates over energy lost from CCs.

We make the assumption that in the steady state the RC locally balances the beam current. That is

$$\begin{eqnarray}&&| {j}_{\mathrm{RC}}(x)| =| {j}_{b}(x)| ,\end{eqnarray} \tag{ 6 }$$

where j_b is the beam current density, and j_RC is the RC density, which can be modeled as the sum of a drifting bulk thermal distribution (j_d ) and the runaway tail of the thermal distribution (j_run):

$$\begin{eqnarray}&&{j}_{\mathrm{RC}}(x)={j}_{d}(x)+{j}_{\mathrm{run}}(x).\end{eqnarray} \tag{ 7 }$$

The electron beam current density, j_b = eF_e ${j}_{b}(x)={{eF}}_{e}(x)\,=e{\int }_{{E}_{c}(x)}^{\infty }F(E,x){dE}$, in the presence of an RC electric field was studied by Holman (2012). He derived an analytic expression describing the beam flux density, F_e (x), as a function of distance along a loop, given the injection of a power-law distribution of electrons. Here we reproduce that result (his Equation (15)) for clarity,

$$\begin{eqnarray}{F}_{e}(x,{E}_{c0})=\left\{\begin{array}{ll}{F}_{e0} & {\rm{;}}x\leqslant {x}_{\mathrm{RC}}\\ {F}_{e0}{\left[\displaystyle \frac{{E}_{\mathrm{th}}+e{\displaystyle \int }_{0}^{x}{{\mathscr{E}}}_{\mathrm{RC}}(x^{\prime} ){dx}^{\prime} }{{E}_{c0}}\right]}^{1-\delta } & {\rm{;}}x\gt {x}_{\mathrm{RC}}\end{array}\right.,\end{eqnarray} \tag{ 8 }$$

where δ is the electron distribution spectral index, F_e0 is the injected electron flux density in square centimeters per second, E_c0 is the low-energy cutoff of the injected electron distribution, x_RC is the distance along the loop where beam electrons are thermalized. Since electrons with an initial energy equal to the low-energy cutoff are thermalized first, x_RC is defined as the distance where the potential drop becomes equal to E_c0 − E_th, i.e., $e{\int }_{0}^{{x}_{\mathrm{RC}}}{{\mathscr{E}}}_{\mathrm{RC}}(x){dx}={E}_{c0}-{E}_{\mathrm{th}}$.

The number density of runaway electrons produced per second at each position x is n_e (x)γ_run(x). These electrons are all accelerated upward by the RC electric field. All runaways produced below position x contribute to the runaway current at x. For our assumed steady-state solution and constant plasma density, the runaway flux and, therefore, current is

$$\begin{eqnarray}&&{j}_{\mathrm{run}}(x)={{en}}_{e}{\int }_{L}^{x}{\gamma }_{\mathrm{run}}(x^{\prime} ){dx}^{\prime} ,\end{eqnarray} \tag{ 9 }$$

where $\tfrac{{\gamma }_{\mathrm{run}}}{{\nu }_{e}}=\tfrac{1}{{\nu }_{e}{n}_{e}}\tfrac{{{dn}}_{\mathrm{run}}}{{dt}}$ is the numerical runaway rate of Landreman et al. (2014) as shown in Figure 1.

For timescales much longer than the collisional timescale, we assume the drifting component of the RC is ohmic, as supported by van den Oord (1990), and is given by

$$\begin{eqnarray}&&{j}_{d}(x)=\displaystyle \frac{{{\mathscr{E}}}_{\mathrm{RC}}(x)}{{\eta }_{s}}.\end{eqnarray} \tag{ 10 }$$

2.3.1. Small Potential Drop: No Thermalization

Let us start with the simple case where the length of the loop is less than x_RC. Substituting the expressions for j_b , j_run, and j_d into Equation (6) gives

$$\begin{eqnarray}&&{{eF}}_{e0}=\displaystyle \frac{{{\mathscr{E}}}_{\mathrm{RC}}(x)}{{\eta }_{s}}-e{\int }_{x}^{L}{n}_{e}(x^{\prime} ){\gamma }_{\mathrm{run}}(x^{\prime} ){dx}^{\prime} .\end{eqnarray} \tag{ 11 }$$

Dividing by ${{\mathscr{E}}}_{D}/{\eta }_{s}$ and setting $z(x)={{\mathscr{E}}}_{\mathrm{RC}}(x)/{{\mathscr{E}}}_{D}(x)$, we obtain the dimensionless expression:

$$\begin{eqnarray}&&\displaystyle \frac{{{eF}}_{e0}{\eta }_{s}}{{{\mathscr{E}}}_{D}}=z(x)+\displaystyle \frac{e{\eta }_{s}}{{{\mathscr{E}}}_{D}}{\int }_{L}^{x}{n}_{e}(x^{\prime} ){\gamma }_{\mathrm{run}}(x^{\prime} ){dx}^{\prime} .\end{eqnarray} \tag{ 12 }$$

We solve for z in Equation (12) using the boundary condition at the footpoints that j_run = 0, so

$$\begin{eqnarray}&&z(x=L)=\displaystyle \frac{{{eF}}_{e0}{\eta }_{s}(L)}{{{\mathscr{E}}}_{D}(L)}.\end{eqnarray} \tag{ 13 }$$

The RC heats the ambient plasma via Joule heating. That is the resistive dissipation of energy from the drifting component. We assume that the runaway component does not significantly heat the ambient plasma.

2.3.2. Large Potential Drop: Electron Thermalization

Now let us consider the case where the loop length is greater than x_RC. As in the previous section, substituting the expressions for j_b , j_run, and j_d into Equation (6), multiplying by ${\eta }_{s}/{{\mathscr{E}}}_{D}$, and using $z(x)={{\mathscr{E}}}_{\mathrm{RC}}(x)/{{\mathscr{E}}}_{D}(x)$, we obtain the dimensionless equation:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{{eF}}_{e0}{\eta }_{s}(x)}{{{\mathscr{E}}}_{D}(x)}{\left[\displaystyle \frac{{E}_{\mathrm{th}}(x)}{{E}_{c0}}+\displaystyle \frac{e}{{E}_{c0}}{\displaystyle \int }_{0}^{x}{{\mathscr{E}}}_{D}(x^{\prime} )z(x^{\prime} ){dx}^{\prime} \right]}^{1-\delta }\\ \quad =\,z(x)-\displaystyle \frac{e{\eta }_{s}(x)}{{{\mathscr{E}}}_{D}(x)}{\displaystyle \int }_{x}^{L}{n}_{e}(x^{\prime} ){\gamma }_{\mathrm{run}}(x^{\prime} ){dx}^{\prime} .\end{array}\end{eqnarray} \tag{ 14 }$$

In this case, the boundary condition and thermalization length are not known a priori, because the solution at any position along the loop is sensitive to the behavior of the runaway current which is integrated from the footpoints upward, whereas the beam current is integrated from the loop top downward. The dynamics of the runaways inform the dynamics of the beam, and vice versa.

We solve this equation iteratively. The solution without thermalization is used to determine the critical values of the low-energy cutoff E_cc and injected electron flux density F_e0c for which the thermalization distance is x_RC = L. For all values of the low-energy cutoff E_c0 < E_cc and F_e0 > F_e0c , Equation (14) should be used. Using Equation (12) as initial guess, and noting that the drift current is proportional to z, this value is updated at each iteration as to minimize ∣j_b − j_RC∣ everywhere along the loop until a satisfactory tolerance is achieved. As in the previous section, the RC heating rate in this case is simply Joule heating from the drifting component. For simplicity we will refer to this model as RA.RC which stands for the RC model with runaway acceleration. The RC model without runaways presented in Holman (2012) is labeled H12.

Since energy lost by beam electrons due to CCs is sensitive to the electrons' energies and the ambient plasma density, for low-energy electrons in a relatively dense loop, energy lost by CC can be comparable to the RC losses and affect the dynamics of the beam/RC system. We derive here the equations which include CC self-consistently with RC losses.

The energy loss equation of a nonthermal beam electron undergoing RC and collisional losses in a hydrogen atmosphere with degree of ionization α is

$$\begin{eqnarray}&&\displaystyle \frac{{dE}}{{dx}}=-e{{\mathscr{E}}}_{\mathrm{RC}}(x)-2\pi {e}^{4}\displaystyle \frac{{n}_{e}(x)}{E}(\alpha {{\rm{\Lambda }}}_{{ee}}+(1-\alpha ){{\rm{\Lambda }}}_{{eH}}),\end{eqnarray} \tag{ 15 }$$

where Λ_ee and Λ_eH are the electron–electron Coulomb logarithm for a fully ionized plasma and the effective Coulomb logarithm for electron-hydrogen atom collisions, $\alpha =\tfrac{{n}_{i}}{{n}_{i}+{n}_{n}}$ is the degree of ionization, where n_i and n_n are the ionized and neutral hydrogen densities, respectively. The second term in Equation (15) is the energy loss by CCs using Equation (23) in Emslie (1978). The beam current density is affected by both RC and CC losses:

$$\begin{eqnarray}{J}_{b}=\left\{\begin{array}{ll}{{eF}}_{e0} & {\rm{;}}x\leqslant {x}_{C}\\ {{eF}}_{e0}{\left(\displaystyle \frac{{E}_{\min }(x)}{{E}_{c0}}\right)}^{1-\delta } & {\rm{;}}x\gt {x}_{C}\end{array}\right.,\end{eqnarray} \tag{ 16 }$$

where x_C is the thermalization distance which is equivalent to x_RC if CC losses are negligible. If RC losses are negligible, x_C is defined as the distance where the energy lost by an electron with initial energy E_c0 becomes equal to E_th, i.e.,

$$\begin{eqnarray*}&&{E}_{\mathrm{th}}^{2}-{E}_{c0}^{2}=-2\pi {e}^{4}{\int }_{0}^{{x}_{C}}n(x^{\prime} )(\alpha {{\rm{\Lambda }}}_{{ee}}+(1-\alpha ){{\rm{\Lambda }}}_{{eH}}){dx}^{\prime} .\end{eqnarray*}$$

${E}_{\min }(x)$ is calculated numerically using Equation (15) and the flux conservation equation F(E, x)dE = F₀(E₀, x = 0)dE₀. When Coulomb collisional losses are negligible,

$$\begin{eqnarray*}&&{E}_{\min }(x)={E}_{\mathrm{th}}(x)+e{\int }_{0}^{x}{{\mathscr{E}}}_{\mathrm{RC}}(x^{\prime} ){dx}^{\prime} .\end{eqnarray*}$$

Similarly, if RC losses are negligible,

$$\begin{eqnarray*}&&\begin{array}{l}{E}_{\min }(x)=\left(\Space{0ex}{3.00ex}{0ex}{E}_{\mathrm{th}}^{2}(x)\right.\\ \quad {\left.+\,4\pi {e}^{4}{\displaystyle \int }_{0}^{x}{n}_{e}(x^{\prime} )(\alpha {{\rm{\Lambda }}}_{{ee}}+(1-\alpha ){{\rm{\Lambda }}}_{{eH}}){dx}^{\prime} \right)}^{1/2}.\end{array}\end{eqnarray*}$$

The dimensionless current balance equation, which includes RC and collisional losses, acceleration of runaway electrons, and thermalization of electrons is the following:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{{eF}}_{e0}{\eta }_{s}(x)}{{{\mathscr{E}}}_{D}(x)}{\left(\displaystyle \frac{{E}_{\min }(x)}{{E}_{c0}}\right)}^{1-\delta }\\ \quad =\,z(x)-\displaystyle \frac{e{\eta }_{s}(x)}{{{\mathscr{E}}}_{D}(x)}{\displaystyle \int }_{x}^{L}{n}_{e}(x^{\prime} ){\gamma }_{\mathrm{run}}(x^{\prime} ){dx}^{\prime} .\end{array}\end{eqnarray} \tag{ 17 }$$

This equation generalizes Equation (14) to include Coulomb losses. It is solved using the same method described in Section 2.3.2, with the difference being that the initial guess for the beam current density is calculated numerically using Equation (16). We will refer to this model as RA.RC.CC, and the case without runaways as RC.CC where RC and CC losses are considered but J_b = − J_d .

We assume that the energy lost by the beam (Equation (15)) is partially gained as heat by the ambient plasma and accelerates the remaining fraction. This is calculated self-consistently. So in this case, the volumetric energy flux lost by the nonthermal beam is

$$\begin{eqnarray}&&q(x)=-\displaystyle \frac{d}{{dx}}\int {EF}(E,x){dE},\end{eqnarray} \tag{ 18 }$$

whereas the Joule heating rate is

$$\begin{eqnarray}&&Q(x)={{\mathscr{E}}}_{\mathrm{RC}}^{2}/{\eta }_{s}.\end{eqnarray} \tag{ 19 }$$

Table 1 summarizes the limiting cases and nomenclature used. The fourth column shows other papers with similar assumptions about the mechanisms taken into account.

As a first example we choose a test demonstrating the case with energy loss less than the thermalization energy (Section 2.3.1). We model the transport of an electron beam with high flux density propagating in a cool/less dense loop , in which CC losses are negligible. We inject a high flux density of beamed electrons to produce a large RC density (since j_RC must balance j_b ). We use a cool loop to increase the resistivity and therefore the electric field. We find that these conditions also show that the effects of runaway electrons are critical.

Figure 2 shows examples of the predicted spatial evolution of RC densities, electric fields, and potential drops, which include and exclude the effects of runaways as defined in Section 2.3.1. In each model, the injected electron beam is characterized by a power law with F_e0 = 2.66 × 10¹⁸ cm ⁻² s⁻¹, E_c = 17 keV, and δ = 4 propagating in a 20 Mm half-loop with apex temperature, T = 3.5 MK, and apex density, n_e = 4.5 × 10⁹ cm⁻³. These beam and background parameters will be referred to as case 1 in the remainder of the paper. In each panel the red and black lines correspond to the RA.RC.CC and RC.CC models, respectively.

Zoom In Zoom Out Reset image size

The sudden decreases in current density (panel (a)) and electric field (panel b) in the RC.CC model are due to the thermalization of electrons out of the beam occurring at about 17 Mm from the loop top. Only the drifting component of the RC density is ohmic, so the inclusion of runaways results in an electric field that is less than that predicted by RC.CC in most of the upper part of the loop. In fact, in this example the electric field produces a total potential drop that is smaller than the thermalization energy, so in the RA model no electrons are lost from the beam throughout their transport. The runaway current increases from the footpoints to the loop top as the electric field accelerates more background electrons, and the drifting current is reduced as a consequence of the current balance assumption. Note that the runaway current density dominates the drift density over much of the loop. At the loop top the runaway flux density is 33% of the total RC flux density.

Finally notice the RC electric field is lower than $0.11{{\mathscr{E}}}_{D}$ everywhere along the loop.

This example shows that for a beam with injected flux density of 2.66 × 10¹⁸ cm⁻² s⁻¹ propagating in a cool, relatively low-density coronal loop, the runaway electrons significantly change the dynamics of the direct beam/RC system by reducing the local electric field.

For the next test case, we have varied the half-length of our model loop. Of course, the longer the distance the beam has to propagate, the higher the resulting potential drop. For this case, we use the same injected electron distribution into loops with the same temperature and density as was used in Figure 2. The results of this experiment are shown in Figure 3, which plots the potential drops for the RA.RC.CC and RC.CC models as a function of loop half-length.

Zoom In Zoom Out Reset image size

In the RC.CC case, the change in potential drop versus loop length slows for loops longer than about 22 Mm. At that height the energy lost due to the potential drop is equal to the thermalization energy, so the beam current and hence RC densities lessen. Note that even though we have considered a maximum loop half-length of 30 Mm, the loop half-length as determined from RHESSI observations is up to 37 Mm (AH17).

The fraction of the RC flux carried by runaway electrons saturates for loop lengths above 25 Mm. This is because no beam electrons are thermalized for loop lengths less than 25 Mm and as beam electrons have to propagate along longer loops, the potential drop increases resulting in thermalization of the lower energy electrons. The longer the loop, the higher (closer to the loop top) the thermalization distance in both models with and without runaways. However, in the runaway case, the potential drop gained by traveling along a longer loop is compensated for by the loss of electrons from the beam.

As a final experiment in this section, we compare the RA.RC.CC and RC.CC models. We have constructed this experiment to demonstrate a case in which beamed electron thermalization begins to occur in the coronal portion of the loop (i.e., applicable to the theories in Section 2.3.2 and Section 2.4) rather than the footpoint (as was the case in Figure 2). We inject an electron beam characterized by flux density, power-law index, and cutoff energy given by, F_e0 = 2 × 10¹⁹ cm⁻² s⁻¹, δ = 4, and E_c0 = 20 keV, respectively. The corresponding injected energy flux density is 1 × 10¹² erg cm⁻² s⁻¹. The beam is injected into a 20 Mm half-length loop with apex temperature and density of T = 9 MK n_e = 2 × 10¹⁰ cm⁻³, respectively. These beam and background parameters will be referred to as case 2 in the remainder of the paper. The results of this experiment are shown in Figure 4. In the RA.RC.CC case, the electric field reaches a maximum at the position where electron thermalization begins to occur, x = 9.4 Mm. As electrons are thermalized, the beam and RC densities and therefore electric field begin to drop.

Zoom In Zoom Out Reset image size

The local contribution to the runaway electrons, γ_run, (most easily inferred from the analytical solution of Equation (2)), depends on the normalized electric field and the collision frequency. In some cases the local contribution of runaways is a small fraction of the local drifting Maxwellian, but the cumulative buildup of runaways along the loop means that as they reach the loop top, their total contribution may be substantial. This is the case for the experiment plotted in Figure 2. In other circumstances, the local runaway contribution may become a significant fraction of local Maxwellian distribution (i.e., where the square of the ratio of the thermal speed to the critical speed is not much less than 1, resulting in ${{\mathscr{E}}}_{\mathrm{RC}}/{{\mathscr{E}}}_{D}$ close to 1). In these cases, the assumption to treat the drifting and runaway components as separate distributions breaks down. The RC density, J_RC, should be calculated as a moment of a combined distribution. Since the runaway growth rate also depends on the collision frequency (Equation (4)), for lower temperature and/or higher density, more runaway electrons are accelerated locally for the same normalized electric field strength. An example demonstrating how the electron distribution differs from a Maxwellian as a function of normalized electric field strength is shown Figure 3 in Stahl et al. (2017).

It is important to recognize the difference between the local and cumulative fraction of runaway electrons. The local accelerated value depends on the local density, temperature, and magnitude of the electric field, whereas the fraction of runaway electrons at the loop top is the sum of all accelerated runaways. By the same token, the minimum energy of electrons in the runaway population at position x that were accelerated out of thermal distribution at position L is

$$\begin{eqnarray}&&{E}_{* }(x)=\displaystyle \frac{m}{2}{v}_{{cr}}^{2}(x)+{\int }_{L}^{x}e{{\mathscr{E}}}_{\mathrm{RC}}(x^{\prime} ){dx}^{\prime} ,\end{eqnarray} \tag{ 20 }$$

where v_cr is the critical velocity in Equation (1). E_*(x) depends on how the electric field and the Dreicer field vary along the loop. The critical energy is the lowest energy of runaway electrons. In a system with finite number of electrons and spatially homogeneous, such as a laboratory plasma, the lowest energy of runaways is given by the first term on the right side of Equation (20). This is the case if the system reaches a steady state, i.e., for ${\mathscr{E}}\ll {{\mathscr{E}}}_{D}$. However, since the beam/RC system in our case evolves along the loop and continuously draws electrons from the chromosphere to form the RC, the lowest energy of runaway electrons is given by Equation (20): The second term on the right side gives the energy gain by runaway electrons due to the electric field acceleration at position x. For values of the RC electric field strength ${{\mathscr{E}}}_{\mathrm{RC}}\lesssim {{\mathscr{E}}}_{D}$, the lowest energy of runaway electrons is determined by the local distribution function of the background plasma. This is because this distribution is not simply the sum of a drifting and a runaway components, as is the case in this paper and Equation (20).

For the examples shown in Figures 2 and 4, the minimum runaway energies at the loop top are given by E_*(0) = 13 and E_*(0) = 21 keV, respectively. These energies are of the same order of magnitude as the accelerated nonthermal beam electrons. For comparison, the thermal energy in the first and second cases are 0.3 and 0.8 keV.

For the experiment in Figure 4, note that even though the total potential drops from the loop top to the footpoints are similar in the RA.RC.CC and RC.CC models (20 and 21 kV, respectively), the nature of the electrons at the loop top (and assumed to be returning to the acceleration region) is significantly different. In the RA.RC.CC model, 22.7% of the electrons returning to the acceleration region are suprathermal with a minimum energy of 21 keV. Note that the assumption that electrons return to the acceleration region is reasonable if the time to reach the steady state is lower than the acceleration time. The time to reach the steady state in atmospheres 1 and 2 are 21 and 20 ms, respectively.

The energy gained by the runaway electrons at the loop top is equal to the total potential drop and is just over 20 keV (second term in Equation (20)). The minimum energy in the runaway current distribution depend on the values of E_cr and the actual shape of the drift current electron distribution. Given that the electrons start with an energy of at least 0.8 keV, 21.8 keV is the energy of the bulk of the runaway electron distribution at the top of the loop.

Here we include some more detailed comparisons of results from our analytical (CC.only and RC.only) and numerical (RC.CC, RA.RC, and RA.RC.CC) models. In Figures 5(a) and (b), we show the volumetric heating rates and low-energy cutoffs as functions of position for our five models listed in Table 1. The loop half-length, apex temperature, and density, and injected beam distribution are identical to those used in the experiment presented in Figure 2 (case 1). Figures 5(c) and (d) show similar quantities for loop conditions and injected beam distribution identical to the case presented in Figure 4 (case 2). The case in the left panels is chosen to illustrate conditions where energy losses from the RC are significant but CC losses are negligible. In the right panels, RC losses, CC, and runaways, all significantly affect the solution.

Zoom In Zoom Out Reset image size

Runaway electrons reduce Joule heating overall, panel (a) in Figure 5 shows that, in the models that exclude runaways, the heating in the upper corona is overestimated by a factor of 2.4, while in the lower corona the heating rate is underestimated by as much as a factor of 1.7 when compared to models that include runaways. In panel (b), for the RA.RC.CC and RA.RC models, the total potential drop is reduced enough by the presence of runaways for all electrons to reach the thick target without being thermalized. In contrast, in the RC.CC model, the beam flux loses ∼32% of its electrons to thermalization before it reaches the footpoints (panel (a) of Figure 2). This results in a factor of 1.7 decrease in the heating rate in the lower corona (panel (a) of Figure 5). As already mentioned, Coulomb collisional losses are negligible in this example, as evidenced by little change in the low-energy cutoff and nearly zero heating for the CC.only case. But they do have create a small increase of the heating rate of RC.CC compared to RC.only.

For the case presented in panels (c) and (d), the heating rate is lower by only a factor 1.8 in the upper corona in the runaway case (RA.RC.CC) compared to the model without runaways (RC.CC). As expected, the thermalization distance is lowest when RCs and CCs are considered without runaways. The thermalization distance in the RA.RC.CC model is lower down the loop because RC losses dominate over CCs even when runaways are accounted for. For the models excluding CC losses (RA.RC and RC.only), the thermalization distance is lower in the case with highest energy losses, i.e., without runaways. Finally, note that since the density is higher in loop used in the right panels compared that used in the left panels, the heating by CC losses is more significant. This is demonstrated in panel (d), for the CC.only case an electron with initial energy of 20 keV loses 7 keV before it reaches the footpoint. These examples show that runaway electrons decrease the heating rate compared to collisional RC models (RC.CC and RC.only), but the heating is significantly higher than models neglecting RC losses altogether. We chose relatively high but reasonable injected flux densities (as inferred by spectral fitting of X-rays observed during solar flares) to demonstrate the effect of runaway electrons. In the next section, we show how these cases compare to other injected beam flux densities and coronal atmospheres.

Finally, Figure 6 shows the total volumetric power lost by the nonthermal beam as calculated by Equation (18) and the volumetric power gained as heat by the ambient plasma as calculated by Equation (19). The remaining energy lost by the beam served to accelerate ambient electrons into runaways. Note that the total energy lost by the nonthermal beam in the runaway model is lower than that lost by the beam in the model without runaways (RC.CC, dashed red curves in panels (a) and (c) of Figure 5).

Zoom In Zoom Out Reset image size

Interestingly, the influence of runaway electrons on the RC electric field is self-reducing. A higher injected flux produces a higher electric field, resulting in a higher rate of runaways. But a higher contribution of runaways decreases the RC electric field, thereby decreasing the runaway rate. Thus, the solution is found where the self-regulating effect is balanced. This can be seen in the upper panel of Figure 7. Increasing the injected electron flux density results in the eventual decrease in electric field between the loop top and footpoints if no beam electrons are thermalized. However, if thermalization occurs, the potential drop between the loop top and footpoints is reduced but the local electric field induced by the lower injected flux density beam is higher below the thermalization distance. For the atmosphere 1 case, this result is seen for beam fluxes above ∼2 × 10¹⁸ electrons cm⁻² s⁻¹, and in the atmosphere 2 case above ∼1.3 × 10¹⁹ electrons cm⁻² s⁻¹ (dashed blue curve).

Another self-reducing effect is caused by the thermalization of beamed electrons. A larger beamed electron flux produces a stronger electric field, which thermalizes more electrons out of the beam, thereby reducing the electric field. This can be inferred by differences in the potential drops between the E_c0 = 20 and E_c0 = 50 keV cases that are plotted in the top panel of Figure 7. For the 20 keV case electrons are thermalized for lower injected flux densities, whereas in the 50 keV case many beam electrons reach the footpoints without thermalization, even for the highest injected fluxes considered. Increasing the injected flux density, moves the thermalization distance closer to the injection site, so the beam flux density decreases producing a weaker electric field.

In this section, we discuss the relative effect of our 1D assumption (i.e., assuming all electrons move parallel to the magnetic field) compared to a 1.5D assumption in which electrons may have nonzero pitch angles. We aim to compare the fraction of runaways at the loop top (using RA.RC.CC) to the reversed nonthermal beam electrons (using FP 1.5D RC&CC).

Because of the RC electric field, initially beamed electrons systematically lose the component of their energy parallel to the magnetic field. They are decelerated and turned backward if the potential drop is greater than their parallel energy. Additionally, CCs scatter these electrons in all directions, producing a broader pitch-angle distribution. Scattering, the reversal due to deceleration, and the reacceleration by the RC electric field are significant. These reversed electrons are effectively part of the RC, reducing the needed balancing drifting component, the RC electric field, and RC losses.

To compare the electron beam transport in 1.5D with our 1D model, we use the Fokker–Planck kinetic transport model, FP, developed by Allred et al. (2020). FP models the distribution of nonthermal electrons injected at the loop top as they transport to the footpoints, under the influence of CCs and pitch-angle scattering and the RC electric field force, but without considering the runaway component. Figure 8 shows an example of this comparison. The loop temperature and density stratification are shown in the upper right panel (Allred et al. 2005, 2015). We label FP configured to include or exclude forces changing the pitch angle as FP 1.5D and FP 1D, respectively. We inject a beam with 2.6 × 10¹⁸ electrons cm⁻² s⁻¹, a low-energy cutoff of 25 keV and a spectral index δ = 3. First, note the similarities in current density, electric field, and heating rate of the 1D models excluding runaways (RC.only and FP 1D). Since, FP 1D models the CC force and RC.only does not, thermalization is predicted to occur higher in the corona in FP 1D than in RC.only. The heating rate rises slowly over the first ∼18 Mm because the decreasing temperature results in increased resistivity. The decrease in heating rate beginning between 18 and 21 Mm is due to thermalization of the beam electrons. The decrease is smoother in the FP 1D case because energy diffusion is considered and treated self-consistently.

Zoom In Zoom Out Reset image size

Since FP 1.5D models the evolution of the pitch angle, μ, it is necessary to specify an injected pitch-angle distribution. To best compare with our 1D models, in FP 1.5D we use a distribution with all injected electrons having zero pitch angle (i.e., μ₀ = 0). The heating rate predicted by FP 1.5D is dominated by RC losses but is reduced by a factor 1.2 at the loop top and up to 3.3 at 19 Mm compared to the 1D case because the pitch-angle scattering and subsequent potential drop reduce the mean downward velocity of the beam electrons. The initial sudden drop in current density and heating rate between 0 and 2 Mm is due to collisional scattering and the reduction of the heating rate for distances just greater than 2 Mm is due to the deceleration by the RC electric field.

The upper panel Figure 9 shows the FP heating rates of the 1.5D and 1D models, in blue and black, respectively. The solid and dashed curves include both RC and CC forces, and only CC frictional forces, respectively. Note the decrease in heating rate within 2 Mm of the loop top in both 1.5D solutions, with and without RCs, as evidence for the scattering being responsible for the initial drop in current density and heating rate from Figure 8. The lower panel of Figure 9, shows the downward and upward propagating electrons. The latter are 19% of total RC flux returning to the loop top.

Zoom In Zoom Out Reset image size

The current density in FP 1.5D is reduced compared to the FP 1D case (upper left panel of Figure 8). The initial scattering of electrons reduces the downward nonthermal beam flux, which results in a lower RC electric field compared to the 1D case. For the same reason the thermalization of beam electrons starts farther from the loop top. The difference in thermalization distance is also affected by the fact that electrons with nonzero pitch angles lose less energy by the RC deceleration.

Finally, we find that since the RC electric field predicted by the RA.RC.CC model is smaller throughout the top half of the loop compared to all other models, in this case runaways more effectively reduce the electric field than backscattering or thermalization. The heating rate predicted by RA.RC.CC is lower than FP 1D and RC.only by almost a factor 3 in the upper corona. Additionally, the fraction of the RC in the runaway regime at the loop top is 39% (as inferred from the upper left panel of Figure 8) compared to 19% of the injected nonthermal electrons becoming reversed as predicted by FP 1.5D (lower panel Figure 9).

The fraction of reversed electrons depends on the injected pitch-angle distribution as well, and has been studied by Zharkova et al. (1995, 2010), Siversky & Zharkova (2009), and Zharkova & Dobranskis (2016). The authors find that the fraction of upward propagating beam electrons is higher closer to the chromosphere. Of course the collisional scattering is stronger lower in the atmosphere and the beam electrons are decelerated by the electric field (e.g., lower panel of Figure 9 and Tables 1–3 in Zharkova & Dobranskis 2016). Conversely, the fraction of runaway electrons accelerated out of the RC plasma is highest at the loop top as more runaways are accumulated from the chromosphere toward the primary acceleration region in the corona. Refer for example to the upper left panel of Figure 8. The runaway electrons accelerated at the footpoint gain an energy equal to the total potential drop.

The total potential drop in the corona is 35 kV in RA.RC.CC versus 38 kV in FP 1.5D. The current density is constant from the loop top to the thermalization distance in all 1D models (RA.RC.CC, RC.only, and FP 1D).

The reversal of nonthermal beam electrons by the RC electric field and collisions was also studied by Karlicky (1993), who recognized that both reversed electrons and runaways should be studied self-consistently but focused on the reversed component. He found that these electrons are higher by an order of magnitude compared to the runaway electrons accelerated out of the background plasma distribution. However, his estimate of the runaways density is simply the fractional part of the total density in the runaway tail (Norman & Smith 1978). As shown by our model however, the runaways are accelerated and accumulate along the full length of the loop, and the RC electric field evolves as a function of the total runaway density.

Although runaway electrons can dominate the RC, their number will be affected by the reversed electrons as they both contribute to reducing the effect of RCs and the heating rates. While both contribute suprathermal electrons returning to the loop top, the number returning and the reduction of the RC electric field is larger in our runaway model. These are nonlinear effects however, so they should be studied self-consistently.

In addition, two problems arise for ${{\mathscr{E}}}_{\mathrm{RC}}\lesssim {{\mathscr{E}}}_{D}$ as opposed to ${{\mathscr{E}}}_{\mathrm{RC}}\ll {{\mathscr{E}}}_{D}$ when considering 1D versus 1.5D treatments: (1) the velocity diffusion becomes non-negligible (Smith & Verwichte 2008) because the critical velocity of electrons approaches the thermal velocity, so that a 1.5D model is more appropriate. In fact the critical velocity is no longer a well-defined value for higher electric fields. (2) The runaway electrons can move into the slide-away regime (Stahl et al. 2016). This is a consequence of locally depleting the bulk distribution which effectively reduces the density, and therefore increases the effective Dreicer field. The result is that the entire background distribution is freely accelerated even though ${{\mathscr{E}}}_{\mathrm{RC}}\lt {{\mathscr{E}}}_{D}$. We will extend the model to 1.5 D (2D in velocity space) in an upcoming paper.

Langmuir wave-beam instabilities can flatten the electron distribution at lower energies (Hannah & Kontar 2011) on timescales orders of magnitude lower than the collisional timescales, which are relevant for the beam/RC dynamics of interest in this paper (e.g., Hannah et al. 2009; Siversky & Zharkova 2009; Holman et al. 2011; Zharkova & Siversky 2011; Karlický & Kontar 2012). Other beam instabilities can also affect RC losses such as the Weibel instability (e.g., Karlický & Kašparová 2009) and pitch-angle instability (e.g., Holman et al. 1982). A more comprehensive simulation is needed to fully understand the effect of beam instabilities and runaway electrons on the beam/RC system.

Rather, this section addresses the question of whether the drift speed is sufficiently low so that the beam/RC system is stable to parallel current-driven instabilities. The goal is to check whether these instabilities may arise under the conditions of our runaway model and discuss the considerations involved in analyzing the nonthermal beam/RC system.

Classical thresholds for current-driven instabilities are calculated for a current propagating in a plasma that is a simpler scenario compared with the coevolution of the heating and the beam/RC system. (e.g., Figure 1 in Holman 1985). Another complication is that the ion temperature is also unknown and some instabilities are sensitive to the electron-to-ion temperature ratio, like the electrostatic ion acoustic instability (Duijveman et al. 1981; Cromwell et al. 1988; Hoyng et al. 1978), while the Buneman and electrostatic ion cyclotron are insensitive to the temperature ratio (e.g., Papadopoulos 1977; Norman & Smith 1978; Benz 2002).

For simplicity, we consider the lowest threshold for a current-driven instability to occur if the electron and ion temperatures are comparable, T_e ∼ T_i , which is about the electron thermal speed v_e , and note that the threshold is lower as the temperature ratio T_e /T_i increases. It is lower by an order of magnitude if T_e ∼ 5 T_i (e.g., Holman 1985; Benz 2002). Note that the thermal speed in Holman is lower than the thermal speed in this paper by a factor $\sqrt{2}$.

We compare the drift velocities to the electron thermal speed calculated in the cases discussed in previous sections. Remember that all cases in this paper satisfy the condition that the injected electron beam density, n_b , is at most 10% of the background density in the upper corona, n_e .

Figure 10 shows the drift speed for the beam/atmosphere cases in Figures 2 and 4, and the corresponding electron thermal speeds in atmospheres 1 and 2 (see Table 2). This shows that the drift speed is lower than the thermal speed in both cases everywhere along the loop.

Zoom In Zoom Out Reset image size

Although not shown in the figure, the drift speeds for all beam cases injected in all atmospheres are lower than the electron thermal speed.

The drift velocity of the RC in the runaway model is reduced compared to purely collisional RC models. However, this does not guarantee that the drift speed is low enough to keep the beam/RC system stable to RC driven instabilities.

If the ratio T_e /T_i is initially greater than ∼5, a current-driven instability is likely to arise. Interestingly, if current-driven instabilities play a role in the beam/RC dynamics, the runaway electrons will be reduced because the collisions are enhanced (e.g., Büchner & Elkina 2005), thereby increasing the Dreicer field. However, for most cases, the drift speed is lower than the thermal speed if T_e ∼ T_i .

A time-dependent treatment might be necessary for higher injected electron flux densities. Consider the following dynamic scenarios. The increase in the electron heating rate due to the RC, decreases the Dreicer field and increases the fraction of runaways. In higher injected flux densities, the thermalized electrons can become a significant fraction of the background electrons. This increases the Dreicer field, and results in fewer runaways. Finally, a current-driven instability may arise, thereby heating the ions and then electrons to higher temperatures.

It is important to understand the underlying physics for higher injected flux densities because if more runaway electrons are accelerated, the heating is reduced in the corona compared to collisional RC models, whereas the heating is increased in the unstable regime. Another important distinction is the nature of electrons returning to the acceleration region, i.e., partially or fully suprathermal in the runaway case, and likely fully thermal in the unstable case.