© The Authors 2024

1. Introduction

Accelerated particles are one of the least understood outcomes of the magnetic field reconnecting in the solar corona. Magnetic reconnection can bring the particle distribution out of thermal equilibrium and produce beams of non-thermal particles travelling along the magnetic field. The beams can, in some cases, contain a significant portion of the released magnetic energy (Lin & Hudson 1971; Emslie et al. 2004, 2012). The effects of accelerated particles are routinely observed in large flares, but their role in smaller energy-release events is unclear.

The interactions between the non-thermal particle beams and their environment have been modelled with ever-increasing levels of sophistication. Early models calculated the mean rate of change in the velocities of the beam particles due to Coulomb collisions with free electrons and protons (Brown 1972; Syrovatskii & Shmeleva 1972). This was later generalised to include collisions with neutral hydrogen, initially for when the fraction of ionised hydrogen is uniform (Emslie 1978), and subsequently for a non-uniform ionisation fraction (Hawley & Fisher 1994). Other studies incorporated the effects of a neutralising return current (Knight & Sturrock 1977; Emslie 1980) and deflection of particle velocities due to variations in magnetic flux density (Leach & Petrosian 1981; Chandrashekar & Emslie 1986). Instead of calculating mean rates of change, Leach & Petrosian (1981) numerically solved the Fokker–Planck equation governing the evolution of a distribution of non-thermal particles. By doing this, they could account for pitch angle diffusion (the randomisation of directions resulting from the stochastic nature of collisions) in their beam transport simulations. Further improvements to their model include the incorporation of relativistic effects and energy losses caused by the emission of synchrotron radiation (Petrosian 1985; McTiernan & Petrosian 1990). While early treatments of collisions typically ignored the thermal motion of the ambient plasma, corrected expressions accounting for the ambient temperature were later employed for the mean rate of change in velocity and the rates of energy and pitch angle diffusion (Emslie 2003; Galloway et al. 2005; Jeffrey et al. 2014).

While the models for propagation of non-thermal particle beams have become highly developed, they have mainly been applied to individual beams in isolated one-dimensional (1D) atmospheres, typically combined with hydrodynamics and radiative transfer for the purpose of simulating flares (Hawley & Fisher 1994; Abbett & Hawley 1999; Allred et al. 2005, 2015, 2020; Liu et al. 2009; Jeffrey et al. 2019) or nanoflares (Testa et al. 2014; Polito et al. 2018; Bakke et al. 2022). In recent times, three-dimensional (3D) atmospheric simulation codes have been developed that reproduce many of the observed features of both the solar corona and lower atmosphere in a small patch of the Sun outside of active regions, by solving the equations of magnetohydrodynamics (MHD) combined with radiative transfer and thermal conduction (Gudiksen et al. 2011; Rempel 2017). However, these simulation codes did not originally account for accelerated particles.

In Frogner et al. (2020, hereafter Paper I), we presented a first approach for integrating the modelling of accelerated electrons into a 3D atmospheric simulation. To let the simulated atmosphere drive the spatial distribution and energetics of the non-thermal electron beams, we detected magnetic reconnection sites and applied a simple parametric model for particle acceleration based on local conditions. We then modelled the collective energy transport by millions of beams by tracing beam trajectories from the reconnection sites and computing the energy deposited into the ambient plasma through collisions between non-thermal electrons and ambient particles. To compute the deposited energy, we employed the relatively basic analytical transport model of Emslie (1978), with the extension by Hawley & Fisher (1994) to support a non-uniform ionisation fraction.

Here we present a more realistic model for the propagation of non-thermal electron beams that is still computationally efficient enough to be applied to millions of beams. The model is based on solving the continuity equation for electron flux by transforming it into a set of ordinary differential equations using the method of characteristics. This method has been applied by Craig et al. (1985), Dobranskis & Zharkova (2014, 2015), and Zharkova & Dobranskis (2016) to find analytical solutions to the continuity equation (although see the correction to Dobranskis & Zharkova 2014 by Emslie et al. 2014). We solve the characteristic equations numerically. By deriving the continuity equation from the Fokker–Planck equation, we show how to incorporate most physical effects accounted for in state-of-the-art models, including magnetic gradient forces, the return current, corrected collisional rates due to the ambient temperature, and contributions to collisions by elements heavier than hydrogen. We additionally consider acceleration along the magnetic field by the ambient electric field, which, to our knowledge, has never been accounted for in models applied to 1D atmospheres. The most notable effects we leave out are collisional randomisation and relativistic effects. We investigate the impact of each new effect on the energy transported by the diverse beams in a 3D atmosphere identical to that used in Paper I.

2. Methods

2.1. Atmospheric simulation

We used a snapshot from a 3D simulation of the solar atmosphere to provide a realistic environment for non-thermal particle acceleration and transport. The simulation was performed using the Bifrost code (Gudiksen et al. 2011), which solves the non-ideal equations of MHD, accounting for energy transport through field-aligned thermal conduction and radiative transfer (Hayek et al. 2010; Carlsson & Leenaarts 2012). The region of the atmosphere simulated begins at the top of the convection zone, 2.5 Mm below the photosphere, and ends in the corona, 14.3 Mm above the photosphere. It extends 24 Mm in each horizontal direction and uses horizontally periodic boundary conditions. The spatial resolution is 31 km horizontally and between 12 and 80 km vertically, with the finest vertical resolution near the transition region and the coarsest in the upper corona.

The simulation reproduces the basic structure of the quiet solar atmosphere, with a chromosphere and coronal loops heated by magnetic reconnection and acoustic shocks resulting from convective motions below the photosphere. In the snapshot employed for this paper, the atmospheric structure has been further shaped by a magnetic flux emergence event where, 137 min of solar time earlier, a 2000 G horizontal magnetic field was injected through the bottom boundary. As it rose to the photosphere, the injected magnetic sheet was broken up by convective motions, enhancing the magnetic field in the network regions between convection cells. The current snapshot, which is the same as that used in Paper I, is characterised by a central bubble of cool plasma (carried upwards by the emerging magnetic field), with increased magnetic reconnection at the upper boundaries of the bubble where the emerging and pre-existing coronal fields meet at a large angle. Further details on this particular simulation can be found in Hansteen et al. (2019). However, we note that the reconnection events in this simulation at the most release energies of order 10²⁵ erg, while most events are much less energetic (see Paper I). Thus, we are operating in the energy regime of nanoflares (Parker 1988).

2.2. Accelerated particles

Magnetic reconnection, the diffusion of magnetic tension occurring at the interface of domains with oppositely directed magnetic field, tends to produce environments favourable for the acceleration of ambient charged particles to speeds far exceeding their initial thermal speeds. The strong rotation of the magnetic field induces an electric field in the diffusion region, which, in addition to driving currents in the plasma, may directly accelerate electrons and ions to non-thermal speeds (Speiser 1965; Litvinenko & Somov 1993). Trapping of charged particles in shrinking magnetic islands (Drake et al. 2006) and scattering within magnetic turbulence (Dmitruk et al. 2003) may respectively lead to first- and second-order Fermi acceleration (Fermi 1954, 1949). Modelling typically indicates that the resulting energy spectra of accelerated particles are shaped as a power-law.

The particles, confined by the Lorentz force to gyrating trajectories around magnetic field lines, may leave the reconnection sites in coherent beams. In this work, we restrict our attention to accelerated electrons, as these are more readily accelerated to non-thermal speeds than ions due to their low mass. Consequently, they are more likely to transport their energy a significant distance. We can represent a beam of non-thermal electrons in terms of its phase-space distribution function f(r, v, t), defined such that f(r, v, t) d³r d³v is the number of beam electrons within the volume d³r around the position r with velocities within d³v of the velocity v at time t. The CGS unit of f is electrons cm⁻³/(cm s⁻¹)³. We assume that the beam reaches a steady state before the ambient plasma can respond to its presence, so we omit the explicit time dependence. By assuming that the distribution is independent of the azimuthal angle and radial distance of the electrons in their gyrating trajectory around the magnetic field line, we can reduce the phase-space coordinates (r, v) down to the three independent coordinates (s, E, μ). Here s is the distance of the electron from its starting position along the field line. E = m_ev²/2 is the kinetic energy, where m_e is the electron mass and v = ∥v∥ is the total speed. μ = cos(β) is the cosine of the pitch angle β between the velocity and the magnetic field, such that ds/dt = μv. With this choice of coordinates, the volume element is d³r = ds dA, where A is the cross-sectional area of the beam, and the velocity volume element is d³v = (v/m_e) dϕ_v dμ dE, where ϕ_v is the azimuthal angle of v. When using E and μ as the independent velocity coordinates, it is convenient to use the field-aligned electron flux spectrum, defined by

$$\begin{array}{c}\hfill F(s,E,\mu )\mathrm{d}\mu \mathrm{d}E={\displaystyle {\int }_{{\varphi }_v}\mu vf(s,\mathit{v}){\mathrm{d}}^{3}v,}\end{array}$$(1)

to represent the beam rather than the phase-space distribution. The quantity in Eq. (1) is the rate of electrons with energies within dE of E and pitch angle cosines within dμ of μ flowing through a unit cross-sectional area in the positive magnetic field direction. The CGS unit of F is electrons s⁻¹ cm⁻² erg⁻¹.

2.2.1. Injected distributions

To specify the field-aligned flux spectrum F₀(E₀, μ₀)≡F(s = 0, E, μ) of non-thermal electrons injected into each beam in the simulated atmosphere, we apply the same simple parametric acceleration model as in Paper I. We classify a grid cell in the simulation box as a reconnection site if the quantity

$$\begin{array}{c}\hfill K=\left||\mathit{B}\times (\mathrm{\nabla }\times \left(\left(\frac{\mathit{E}·\mathit{B}}{\mathit{B}·\mathit{B}}\right)\mathit{B}\right))\right||\end{array}$$(2)

exceeds a small threshold. Here E and B are the local macroscopic electric and magnetic fields, respectively. In MHD theory, K = 0 is a criterion for the magnetic topology to be conserved (Biskamp 2005), so we detect reconnection by finding where this criterion is sufficiently violated, with K exceeding some small threshold (of the order of 10⁻¹⁵ G statV cm⁻²). We discuss how we chose the threshold value in Paper I.

At each reconnection site, we assume that some unspecified time-independent acceleration mechanism maintains a non-thermal power-law tail on top of the local thermal electron distribution. We express this non-thermal component in terms of an injected electron flux spectrum in the form

$$\begin{array}{c}\hfill F_0(E_0\ge E_{\mathrm{c}},{\mu }_0)={\mathcal{F}}_{\mathrm{beam},0}\frac{\delta -2}{{E_{\mathrm{c}}}^2}{\left(\frac{E_0}{E_{\mathrm{c}}}\right)}^{-\delta }\delta ({\mu }_0-\overline{{\mu }_0}),\end{array}$$(3)

where ℱ_beam, 0 is the field-aligned energy flux injected into the beam, E_c is the lower cut-off energy below which the distribution is empty, and δ is the power-law index describing how strongly the distribution diminishes with increasing energy. We assume that every injected non-thermal electron has the same initial pitch angle cosine ${\mu }_0=\overline{{\mu }_0}$. This is expressed through the Dirac delta function $\delta ({\mu }_0-\overline{{\mu }_0})$.

To determine the injected energy flux ℱ_beam, 0, we assume that 20% of the power released from the magnetic field during reconnection goes into electron acceleration (we discuss the choice of this percentage in Paper I). We partition this power between two beams pointed in opposite directions along the magnetic field based on the local alignment of the electric field with the magnetic field. The opposite directions are encoded into the sign that each of the two beams receives for $\overline{{\mu }_0}$ (and ℱ_beam, 0). The power assigned to each beam is then converted to an energy flux ℱ_beam, 0 using the cross-sectional area A of the beam, which is computed based on the extent of the grid cell at the reconnection site. We calculate the lower cut-off energy E_c as the energy where the resulting non-thermal power-law distribution would intersect the local thermal distribution. The resulting value of E_c is approximately proportional to the local temperature T, with E_c between 1 and 3 keV at T = 1 MK. For the power-law index δ, we assume the same fixed value for every beam. Finally, we estimate the magnitude of the initial pitch angle cosine, $|\overline{{\mu }_0}|$, by assuming that the acceleration process does not alter the initial perpendicular velocity component due to randomised thermal motion. In practice, this results in $|\overline{{\mu }_0}|$ taking a value slightly below unity for most beams.

2.2.2. Transport equations

The evolution of the steady-state field-aligned electron flux spectrum F(s, E, μ) with distance s is governed by the Fokker–Planck equation:

$$\begin{array}{cc}\hfill \frac{\partial F}{\partial s}& +{\left(\frac{\mathrm{d}E}{\mathrm{d}s}\right)}_{!\mathrm{C}}\frac{\partial F}{\partial E}+{\left(\frac{\mathrm{d}\mu }{\mathrm{d}s}\right)}_{!\mathrm{C}}\frac{\partial F}{\partial \mu }=(\frac{1}{E}{\left(\frac{\mathrm{d}E}{\mathrm{d}s}\right)}_{!\mathrm{C}}+\frac{1}{\mu }{\left(\frac{\mathrm{d}\mu }{\mathrm{d}s}\right)}_{!\mathrm{C}})F\hfill \\ \hfill & +\frac{\partial }{\partial E}\left(C_{E}F\right)+\frac{\partial }{\partial \mu }\left(C_{\mu }F\right)+\frac{{\partial }^2}{\partial E^2}\left(C_{E^2}F\right)+\frac{{\partial }^2}{\partial {\mu }^2}\left(C_{{\mu }^2}F\right)\hfill \end{array}$$(4)

(see Appendix A for the derivation of this equation from the more general Fokker–Planck equation expressed in terms of f(r, v)). In Eq. (4), (dE/ds)_!C and (dμ/ds)_!C are the rates of change in E and μ with distance s due to non-collisional forces. The main such force is the Lorentz force F_L on an electron due to the macroscopic electric and magnetic field. In CGS units, it is given as

$$\begin{array}{c}\hfill {\mathit{F}}_{\mathrm{L}}=-e(\mathit{E}+\frac{\mathit{v}}{c}\times \mathit{B}),\end{array}$$(5)

where e is the elementary charge and c is the speed of light. The electric field E can accelerate electrons in any direction. We only consider the influence of the component ℰ of the electric field parallel with the magnetic field direction

$$\begin{array}{c}\hfill \mathcal{E}=\frac{\mathit{E}·\mathit{B}}{B},\end{array}$$(6)

where B = ∥B∥. We ignore the transverse component of the electric field. A non-zero transverse component causes the electron orbit to drift with velocity v_d = cE × B/B², either breaking azimuthal symmetry or (if the electric field is azimuthally symmetric) gradually changing the gyroradius. In either case, the drift speed is tiny compared to the parallel speed of the non-thermal electrons. The highest drift speeds in our atmospheric simulation are of order 10³ cm s⁻¹, far below even the average thermal speed.

We consider two potential contributions to the electric field E. The first is the MHD electric field E_fluid, which may have a parallel component ℰ_fluid when the resistivity is significant (as during magnetic reconnection). The second contribution is the electric field produced in response to the charge displacement that occurs when the accelerated electrons leave the acceleration site. This electric field, which we denote E_beam, produces a counter-flowing return current in the background plasma that cancels out the current associated with the beam (Knight & Sturrock 1977; Emslie 1980). It is always aligned with the magnetic field, so E_beam = ℰ_beamB/B. As discussed in Emslie (1980), we can compute ℰ_beam from the beam electron flux. The field-aligned current density associated with the beam at distance s is

$$\begin{array}{c}\hfill J_{\mathrm{beam}}(s)=-e{F}_{\mathrm{beam}}(s),\end{array}$$(7)

where F_beam(s) is the total field-aligned flux of the beam electrons at distance s (given by Eq. (43)). We note that F_beam is positive if the beam travels along the positive magnetic field direction. Assuming that the return current J_r exactly neutralises the beam current, we have J_r(s) =  − J_beam(s). From Ohm’s law, the electric field is thus

$$\begin{array}{c}\hfill {\mathcal{E}}_{\mathrm{beam}}(s)=\eta (s)J_{\mathrm{r}}(s)=\eta (s)e{F}_{\mathrm{beam}}(s),\end{array}$$(8)

where η is the electrical resistivity of the ambient plasma parallel to the magnetic field. We compute η accounting for electron collisions with protons and neutral hydrogen (see, for example, Chae & Litvinenko 2021).

The work the total parallel electric field ℰ = ℰ_fluid + ℰ_beam does on a beam electron moving with velocity ds/dt = μv along the magnetic field line is

$$\begin{array}{c}\hfill {\left(\frac{\mathrm{d}E}{\mathrm{d}t}\right)}_{\mathrm{E}}=-\mu ve\mathcal{E},\end{array}$$(9)

where the subfix E stands for electric. The acceleration along the magnetic field line is

$$\begin{array}{c}\hfill {\left(\frac{\mathrm{d}(\mu v)}{\mathrm{d}t}\right)}_{\mathrm{E}}=-\frac{e\mathcal{E}}{m_{\mathrm{e}}}·\end{array}$$(10)

The magnetic component of the Lorentz force does no work on the electron since v × B is always perpendicular to the velocity v, so we have

$$\begin{array}{c}\hfill {\left(\frac{\mathrm{d}E}{\mathrm{d}t}\right)}_{\mathrm{M}}=0,\end{array}$$(11)

where the subfix M stands for magnetic. The only direct effect of the magnetic force is to produce the angular velocity dϕ_v/dt = eB/(m_ec) responsible for the gyrating motion. If the magnetic field is uniform, it does not affect the parallel velocity of the electron. However, if there is a weak gradient in the magnetic field, analysis of the linearised equation of motion shows that the electron experiences a net parallel acceleration. This acceleration, due to the magnetic force averaged over one gyration period, is given by (e.g. Bittencourt 2004)

$$\begin{array}{c}\hfill {\left(\frac{\mathrm{d}(\mu v)}{\mathrm{d}t}\right)}_{\mathrm{M}}=-\frac{M}{m_{\mathrm{e}}}\frac{\mathrm{d}B}{\mathrm{d}s}=-(1-{\mu }^2)\frac{E}{m_{\mathrm{e}}B}\frac{\mathrm{d}B}{\mathrm{d}s},\end{array}$$(12)

where M is the gyromagnetic moment of the electron. This assumes that the electron experiences no other forces than the magnetic force, which is not strictly true in our model, considering that the electron can be affected by a parallel electric field and collisions. However, any parallel electric field is weak compared to the magnetic field, and a comparison of the gyrofrequency and collisional frequency for typical atmospheric conditions shows that the electrons gyrate many times for every collision. Equation (12) thus remains a reasonable approximation even in the presence of other forces.

In addition to the Lorentz force, we could also include the reaction force due to the emission of gyromagnetic radiation. This reaction force is given by the Abraham–Lorentz formula (e.g. Barut 1980). However, as discussed by Petrosian (1985), who derived expressions for the corresponding rates of change in energy and pitch angle, the influence of the reaction force is significant only for highly relativistic electrons. Hence, there is nothing to gain by accounting for gyromagnetic radiation until the transport model is extended to include relativistic effects. We note that in the simulations performed for this paper, the most relativistic beams have δ = 4 and E_c ≈ 5 keV. Relativistic electrons with a Lorentz factor exceeding 1.1 account for only 1% of the energy flux injected into these beams, so disregarding relativistic effects is justified.

Adding Eqs. (9) and (11), the full rate of change in energy due to non-collisional forces becomes

$$\begin{array}{c}\hfill {\left(\frac{\mathrm{d}E}{\mathrm{d}t}\right)}_{!\mathrm{C}}=-\mu ve\mathcal{E},\end{array}$$(13)

while from Eqs. (10) and (12) the total parallel acceleration becomes

$$\begin{array}{c}\hfill {\left(\frac{\mathrm{d}(\mu v)}{\mathrm{d}t}\right)}_{!\mathrm{C}}=-\frac{e\mathcal{E}}{m_{\mathrm{e}}}-(1-{\mu }^2)\frac{E}{m_{\mathrm{e}}B}\frac{\mathrm{d}B}{\mathrm{d}s}·\end{array}$$(14)

Finally, using ds/dt = μv, expanding the derivative in Eq. (14), applying Eq. (13), we find

$$\begin{array}{cc}\hfill {\left(\frac{\mathrm{d}E}{\mathrm{d}s}\right)}_{!\mathrm{C}}& =-e\mathcal{E},\hfill \end{array}$$(15)

$$\begin{array}{cc}\hfill {\left(\frac{\mathrm{d}\mu }{\mathrm{d}s}\right)}_{!\mathrm{C}}& =-\frac{1-{\mu }^2}{2\mu }(\frac{e\mathcal{E}}{E}+\frac{\mathrm{d}lnB}{\mathrm{d}s})·\hfill \end{array}$$(16)

Returning to the Fokker–Planck equation, the terms involving the coefficients C_E, C_μ, C_E², and C_μ² in Eq. (4) express the evolution of the flux spectrum due to Coulomb collisions between the beam electrons and ambient plasma particles. The coefficients C_E and C_μ govern the advection of the flux spectrum in energy and pitch angle space, while C_E² and C_μ² govern the diffusion. General expressions for all the coefficients are derived in Appendix A. The results for C_E and C_μ are

$$\begin{array}{cc}& C_E=\frac{2\pi e^4}{\mu E}{\displaystyle (\sum _{\mathrm{c}}(\frac{m_{\mathrm{e}}}{m_{\mathrm{c}}}\mathrm{erf}(u_{\mathrm{c}})-(1+\frac{m_{\mathrm{e}}}{m_{\mathrm{c}}})u_{\mathrm{c}}{\mathrm{erf}}^{\prime }(u_{\mathrm{c}})){z_{\mathrm{c}}}^2{n}_{\mathrm{c}}ln{\mathrm{\Lambda }}_{\mathrm{c}}}\hfill \\ \hfill & +{\displaystyle \sum _{\mathrm{N}}{Z}_{\mathrm{N}}{n}_{\mathrm{N}}ln{\mathrm{\Lambda }}_{\mathrm{N}}^{\prime }),}\hfill \end{array}$$(17)

$$\begin{array}{cc}& C_{\mu }=\frac{\pi e^4\mu }{\mu E^2}{\displaystyle (\sum _{\mathrm{c}}(\mathrm{erf}(u_{\mathrm{c}})-G(u_{\mathrm{c}})){z_{\mathrm{c}}}^2{n}_{\mathrm{c}}ln{\mathrm{\Lambda }}_{\mathrm{c}}+\sum _{\mathrm{N}}{Z_{\mathrm{N}}}^2{n}_{\mathrm{N}}ln{\mathrm{\Lambda }}_{\mathrm{N}}^{″}).}\hfill \end{array}$$(18)

Each expression has a sum of contributions from collisions with ambient particles of different species. The subscript c denotes a particular species of charged particles, while N denotes a particular species of neutral particles. For each charged particle species c, m_c is the particle mass, z_c is the charge number, n_c is the number density, and lnΛ_c is the Coulomb logarithm, for which we follow Emslie (1978) and use the expression

$$\begin{array}{c}\hfill ln{\mathrm{\Lambda }}_{\mathrm{c}}=ln\left(\frac{m_{\mathrm{e}}{m}_{\mathrm{c}}}{m_{\mathrm{e}}+m_{\mathrm{c}}}\frac{v^2\eta }{z_{\mathrm{c}}{e}^2}\right),\end{array}$$(19)

where η = v/ν is the electron mean free path and ν is the plasma frequency. The quantity $u_{\mathrm{c}}=v/(\sqrt{2}{v}_{\mathrm{tc}})$ is the beam electron speed normalised by the thermal speed $v_{\mathrm{tc}}=\sqrt{k_{\mathrm{B}}{T}_{\mathrm{c}}/m_{\mathrm{c}}}$, where k_B is the Boltzmann constant, and T_c is the temperature of the species c particles. The functions erf(u) and erf′(u) are, respectively, the error function and its derivative, while

$$\begin{array}{c}\hfill G(u)=\frac{\mathrm{erf}(u)-u{\mathrm{erf}}^{\prime }(u)}{2u^2}·\end{array}$$(20)

For each neutral particle species N, Z_N is the atomic number, and n_N is the number density. The effective Coulomb logarithms $ln{\mathrm{\Lambda }}_{\mathrm{N}}^{\prime }$ and $ln{\mathrm{\Lambda }}_{\mathrm{N}}^{″}$ are associated respectively with friction and velocity diffusion due to collisions with the neutral particles. They are given by (Evans 1955; Snyder & Scott 1949)

$$\begin{array}{cc}& ln{\mathrm{\Lambda }}_{\mathrm{N}}^{\prime }=ln\left(\frac{m_{\mathrm{e}}{v}^2}{\sqrt{2}{I}_{\mathrm{N}}}\right),\hfill \end{array}$$(21)

$$\begin{array}{cc}& ln{\mathrm{\Lambda }}_{\mathrm{N}}^{″}=ln\left(\frac{v}{\sqrt{2}{Z_{\mathrm{N}}}^{1/3}c\alpha }\right),\hfill \end{array}$$(22)

where I_N is the ionisation potential of a species N particle, c is the speed of light, and α is the fine structure constant. We note that while the Coulomb logarithms lnΛ_c, $ln{\mathrm{\Lambda }}_{\mathrm{N}}^{\prime }$, and $ln{\mathrm{\Lambda }}_{\mathrm{N}}^{″}$ technically vary with beam electron speed v, this dependence is relatively slight, so we ignore it and use a representative value corresponding to the mean speed of the initial distribution. Furthermore, the variations due to the specifics of the particular particle species (m_c, z_c, I_N, and Z_N) are even more minor than the variations due to the electron speed. We, therefore, use the electron’s value of lnΛ_c for all charged particle species and the hydrogen atom’s value of $ln{\mathrm{\Lambda }}_{\mathrm{N}}^{\prime }$ and $ln{\mathrm{\Lambda }}_{\mathrm{N}}^{″}$ for all neutral particle species and omit the subscripts.

In this work, we account for collisions with free electrons (c = e), free protons (c = p), singly and doubly ionised helium atoms (c = HeI and c = HeII), neutral hydrogen atoms (N = NH), and neutral helium atoms (N = NHe). Because the ambient hydrogen and helium atoms have minimal velocities compared to beam electrons, we have u_c ≫ 1 for all c except c = e. We can thus use the asymptotic values erf(u)→1, erf′(u)→0 and G(u)→0 for u → ∞ for collisions with all charged atoms. Since only u_e remains relevant, we use the symbol u as a shorthand for u_e. We can neglect the terms in Eq. (17) that include the factor m_e/m_c for collisions with the charged atoms since the hydrogen mass m_H and helium mass m_He greatly exceed the electron mass. The number densities of the charged atoms can be expressed as n_p = x_Hn_H, n_HeI = x_HeIn_He, and n_HeII = x_HeIIn_He, where x_H, x_HeI, and x_HeII are ionisation fractions, and n_H and n_He are the total number densities of respectively hydrogen and helium (including both charged and neutral atoms). For the neutral atoms, we can then write the number densities as n_NH = (1 − x_H)n_H and n_NHe = (1 − x_HeI − x_HeII)n_He. To express all number densities in terms of the hydrogen density n_H, we further define the electron-to-hydrogen ratio r_e and helium-to-hydrogen ratio r_He such that the electron density becomes n_e = r_en_H and the helium density becomes n_He = r_Hen_H. Equations (17) and (18) then become

$$\begin{array}{cc}& C_E=\frac{2\pi e^4}{\mu E}{n}_{\mathrm{H}}{\gamma }_E(u),\hfill \end{array}$$(23)

$$\begin{array}{cc}& C_{\mu }=\frac{\pi e^4}{E^2}{n}_{\mathrm{H}}{\gamma }_{\mu }(u),\hfill \end{array}$$(24)

where

$$\begin{array}{cc}\hfill {\gamma }_E(u)& =W_E(u)r_{\mathrm{e}}ln\mathrm{\Lambda }\hfill \\ \hfill & +(1-x_{\mathrm{H}}+2(1-x_{\mathrm{HeI}}-x_{\mathrm{HeII}})r_{\mathrm{He}})ln{\mathrm{\Lambda }}^{\prime }\hfill \end{array}$$(25)

and

$$\begin{array}{cc}\hfill {\gamma }_{\mu }(u)=& (W_{\mu }(u)r_{\mathrm{e}}+x_{\mathrm{H}}+(x_{\mathrm{HeI}}+4x_{\mathrm{HeII}})r_{\mathrm{He}})ln\mathrm{\Lambda }\hfill \\ \hfill & +(1-x_{\mathrm{H}}+4(1-x_{\mathrm{HeI}}-x_{\mathrm{HeII}})r_{\mathrm{He}})ln{\mathrm{\Lambda }}^{″}\hfill \end{array}$$(26)

can be considered effective Coulomb logarithms. The functions

$$\begin{array}{cc}& W_E(u)=\mathrm{erf}(u)-2u{\mathrm{erf}}^{\prime }(u),\hfill \end{array}$$(27)

$$\begin{array}{cc}& W_{\mu }(u)=\mathrm{erf}(u)-G(u)\hfill \end{array}$$(28)

encapsulate the dependence of the collision coefficients on the ambient temperature.

Equation (4) is a linear second-order partial differential equation. It can be solved numerically in many ways, including finite difference methods (e.g. Allred et al. 2020) and stochastic methods (e.g. Jeffrey et al. 2020). Unfortunately, these methods are too computationally expensive for simulating huge numbers of beams. To simplify the problem, we ignore the collisional diffusion of the distribution in velocity space by setting the diffusion coefficients C_E² and C_μ² to zero (we consider the consequences of this approximation in our discussion in Sect. 4). Equation (4) can then be written as

$$\begin{array}{c}\hfill \frac{\partial F}{\partial s}+\frac{\partial }{\partial E}\left(({\left(\frac{\mathrm{d}E}{\mathrm{d}s}\right)}_{!\mathrm{C}}-C_E)F\right)+\frac{\partial }{\partial \mu }\left(({\left(\frac{\mathrm{d}\mu }{\mathrm{d}s}\right)}_{!\mathrm{C}}-C_{\mu })F\right)=S,\end{array}$$(29)

where

$$\begin{array}{c}\hfill S=((\frac{1}{E}+\frac{\partial }{\partial E}){\left(\frac{\mathrm{d}E}{\mathrm{d}s}\right)}_{!\mathrm{C}}+(\frac{1}{\mu }+\frac{\partial }{\partial \mu }){\left(\frac{\mathrm{d}\mu }{\mathrm{d}s}\right)}_{!\mathrm{C}})F.\end{array}$$(30)

Equation (29) is the continuity equation for the flux spectrum in phase space. The right-hand side, given by Eq. (30), is the source term, which describes the addition or removal of electron flux due to the macroscopic electric field. Without a parallel electric field, the source term evaluates to zero, meaning that the flux spectrum is conserved as it evolves in phase space.

To solve the continuity equation, we follow Craig et al. (1985) and transform it into a set of ordinary differential equations using the method of characteristics. The characteristics are curves in the (s, E, μ) coordinate space that the flux spectrum F follows, given a set of initial conditions. It is convenient to parametrise the curves by the distance s. Each curve then describes the evolution of energy E(s, E₀, μ₀), pitch angle cosine μ(s, E₀, μ₀), and flux F(s, E₀, μ₀) with distance s for a group of beam electrons given their initial values E₀, μ₀, and F₀ at s = 0. The set of ordinary differential equations governing these characteristic curves for Eq. (29) are given by

$$\begin{array}{cc}& \frac{\mathrm{d}E}{\mathrm{d}s}={\left(\frac{\mathrm{d}E}{\mathrm{d}s}\right)}_{!\mathrm{C}}-C_E,\hfill \end{array}$$(31)

$$\begin{array}{cc}& \frac{\mathrm{d}\mu }{\mathrm{d}s}={\left(\frac{\mathrm{d}\mu }{\mathrm{d}s}\right)}_{!\mathrm{C}}-C_{\mu },\hfill \end{array}$$(32)

$$\begin{array}{cc}& \frac{\mathrm{d}F}{\mathrm{d}s}=(\frac{1}{E}{\left(\frac{\mathrm{d}E}{\mathrm{d}s}\right)}_{!\mathrm{C}}+\frac{1}{\mu }{\left(\frac{\mathrm{d}\mu }{\mathrm{d}s}\right)}_{!\mathrm{C}}+\frac{\partial C_E}{\partial E}+\frac{\partial C_{\mu }}{\partial \mu })F.\hfill \end{array}$$(33)

Inserting Eqs. (15), (16), (23), and (24), this becomes

$$\begin{array}{cc}& \frac{\mathrm{d}E}{\mathrm{d}s}=-\frac{2\pi e^4}{\mu E}{n}_{\mathrm{H}}{\gamma }_E(u)-e\mathcal{E},\hfill \end{array}$$(34)

$$\begin{array}{cc}& \frac{\mathrm{d}\mu }{\mathrm{d}s}=-\frac{\pi e^4}{E^2}{n}_{\mathrm{H}}{\gamma }_{\mu }(u)-\frac{1-{\mu }^2}{2\mu }(\frac{e\mathcal{E}}{E}+\frac{\mathrm{d}lnB}{\mathrm{d}s}),\hfill \end{array}$$(35)

$$\begin{array}{cc}& \frac{\mathrm{d}F}{\mathrm{d}s}=-(\frac{2\pi e^4}{\mu E^2}{n}_{\mathrm{H}}{\gamma }_F(u)+\frac{1+{\mu }^2}{2{\mu }^2}\frac{e\mathcal{E}}{E}+\frac{1-{\mu }^2}{2{\mu }^2}\frac{\mathrm{d}lnB}{\mathrm{d}s})F,\hfill \end{array}$$(36)

where

$$\begin{array}{cc}\hfill {\gamma }_F(u)& =W_F(u)r_{\mathrm{e}}ln\mathrm{\Lambda }\hfill \\ \hfill & +(1-x_{\mathrm{H}}+2(1-x_{\mathrm{HeI}}-x_{\mathrm{HeII}})r_{\mathrm{He}})ln{\mathrm{\Lambda }}^{\prime }\hfill \end{array}$$(37)

and

$$\begin{array}{c}\hfill W_F(u)=\mathrm{erf}(u)-(2u^2+\frac{3}{2})u{\mathrm{erf}}^{\prime }(u).\end{array}$$(38)

The warm-target contribution functions W_E, W_μ, and W_F are plotted in Fig. 1. In the cold-target limit, when the beam electron speeds are very high compared to the thermal electron speeds (u ≫ 1), they all approach unity. When u decreases below ∼4 (this corresponds to a beam electron energy of 2 to 3 keV for typical coronal temperatures in our simulation) the functions deviate from unity. As the energy approaches k_BT, W_E and W_μ decrease, suppressing the collisional rate of energy loss and pitch angle increase. At E ≈ k_BT, the beam electrons experience no average loss in energy from collisions with ambient electrons, and for lower energies, they start to gain energy. As evident from the minimum in W_F at E ≈ k_BT (giving a maximally positive contribution to dF/ds), the influence of the thermal distribution on collisions causes more beam electrons to find themselves with energies close to the mean thermal energy.

Fig. 1. Plots of the warm-target contribution functions WE (orange), Wμ (purple), and WF (green) as functions of the normalised speed $u=v/(\sqrt{2}{v}_{\mathrm{t}})=\sqrt{E/k_{\mathrm{B}}T}$. The dotted grey curves show the mapping between u and energy E for different temperatures T. The dashed blue line highlights the location of u = 1, corresponding to E = kBT.

2.2.3. Computing macroscopic quantities

The ultimate objective of determining the electron flux spectrum F(s, E, μ) is to compute macroscopic quantities associated with the non-thermal electron beam. We can evaluate a macroscopic quantity at a specific distance s by appropriately weighting F(s, E, μ) and integrating it over E and μ. In general, we can compute the sum Σ(χ) of any individual electron property χ(s, E, μ) over all the non-thermal electrons at the distance s by evaluating the following integral:

$$\begin{array}{c}\hfill \mathrm{\Sigma }(\chi )={\displaystyle {\int }_E{\int }_{\mu }\chi (s,E,\mu )\frac{F(s,E,\mu )}{\mu v(E)}\mathrm{d}\mu \mathrm{d}E.}\end{array}$$(39)

From Eq. (1), we see that this is equivalent to weighting the phase-space distribution function f by χ and integrating over velocity space. Performing a change of variables from μ to μ₀, we can write the integral as

$$\begin{array}{c}\hfill \mathrm{\Sigma }(\chi )={\displaystyle {\int }_E{\int }_{{\mu }_0}\chi (s,E,{\mu }_0)\frac{F(s,E,{\mu }_0)}{\mu (s,E,{\mu }_0)v(E)}|\frac{\partial \mu (s,E,{\mu }_0)}{\partial {\mu }_0}|\mathrm{d}{\mu }_0\mathrm{d}E.}\end{array}$$(40)

By using (E, μ₀) as primary independent variables instead of (E, μ) and performing the integral over E and μ₀, we can take advantage of our assumption that every electron in the distribution starts with the same pitch angle cosine ${\mu }_0=\overline{{\mu }_0}$. The Dirac delta function $\delta ({\mu }_0-\overline{{\mu }_0})$ in Eq. (3) is present in the flux spectrum at all distances s. Factoring it out, we can write

$$\begin{array}{c}\hfill F(s,E,{\mu }_0)=\overline{F}(s,E,{\mu }_0)\delta ({\mu }_0-\overline{{\mu }_0}).\end{array}$$(41)

Inserting this into Eq. (40) and performing the integral over μ₀, we get

$$\begin{array}{c}\hfill \mathrm{\Sigma }(\chi )={\displaystyle {\int }_E\chi (s,E,\overline{{\mu }_0})\frac{\overline{F}(s,E,\overline{{\mu }_0})}{\mu (s,E,\overline{{\mu }_0})v(E)}|\frac{\partial \mu (s,E,\overline{{\mu }_0})}{\partial \overline{{\mu }_0}}|\mathrm{d}E.}\end{array}$$(42)

The direction-integrated flux spectrum $\overline{F}(s,E,\overline{{\mu }_0})$ corresponds to the flux spectrum F(s, E, μ₀) integrated over μ₀. We have thus reduced the problem of solving for F, which is two-dimensional in velocity space, to solving for $\overline{F}$, which is one-dimensional in velocity space, with the additional complication of having to determine the Jacobian determinant $|\partial \mu (s,E,\overline{{\mu }_0})/\partial \overline{{\mu }_0}|$.

In terms of the generic integral in Eq. (42), the total field-aligned electron flux at the distance s can be written

$$\begin{array}{c}\hfill F_{\mathrm{beam}}(s)=\mathrm{\Sigma }(\mu v),\end{array}$$(43)

while the field-aligned energy flux becomes

$$\begin{array}{c}\hfill {\mathcal{F}}_{\mathrm{beam}}(s)=\mathrm{\Sigma }(E\mu v).\end{array}$$(44)

The non-thermal electron number density n_beam(s) at the distance s is simply

$$\begin{array}{c}\hfill n_{\mathrm{beam}}(s)=\mathrm{\Sigma }(1).\end{array}$$(45)

The corresponding integrand in Eq. (42) is

$$\begin{array}{c}\hfill \frac{\mathrm{d}n}{\mathrm{d}E}=\frac{\overline{F}(s,E,\overline{{\mu }_0})}{\mu (s,E,\overline{{\mu }_0})v(E)}|\frac{\partial \mu (s,E,\overline{{\mu }_0})}{\partial \overline{{\mu }_0}}|·\end{array}$$(46)

This quantity, the number density of non-thermal electrons per energy interval, provides an informative way of visualising the distribution.

The total power density Q_beam(s) deposited by the non-thermal electrons at the distance s is found by integrating up the collisional rate of energy loss for the individual electrons:

$$\begin{array}{c}\hfill Q_{\mathrm{beam}}(s)=\mathrm{\Sigma }(-{\left(\frac{\mathrm{d}E}{\mathrm{d}t}\right)}_{\mathrm{coll}}).\end{array}$$(47)

The collisional rate of change in energy is found by taking the collisional part of Eq. (34) and multiplying with ds/dt = μv:

$$\begin{array}{c}\hfill {\left(\frac{\mathrm{d}E}{\mathrm{d}t}\right)}_{\mathrm{coll}}=-\frac{2\pi e^{4}v}{E}{n}_{\mathrm{H}}{\gamma }_E(u).\end{array}$$(48)

We find the total heating power per volume Q(s) resulting from the beam by adding the contribution Q_r(s) from the resistive heating due to the return current J_r,

$$\begin{array}{c}\hfill Q(s)=Q_{\mathrm{beam}}(s)+Q_{\mathrm{r}}(s),\end{array}$$(49)

where

$$\begin{array}{c}\hfill Q_{\mathrm{r}}(s)=\eta (s){J_{\mathrm{r}}(s)}^2=e^2\eta (s){F_{\mathrm{beam}}(s)}^2.\end{array}$$(50)

2.2.4. Numerical solution

For each non-thermal electron beam, we obtain its trajectory by tracing the magnetic field line from the reconnection site in the appropriate direction (as determined by the sign of $\overline{{\mu }_0}$) in the same manner as in Paper I. In Frogner & Gudiksen (2022), we covered the procedure for accurately tracing the magnetic field lines in detail. In tandem with tracing the trajectory, we integrate Eqs. (34)–(36) simultaneously to obtain the evolution of energy $E(s,E_0,\overline{{\mu }_0})$, pitch angle cosine $\mu (s,E_0,\overline{{\mu }_0})$ and flux $\overline{F}(s,E_0,\overline{{\mu }_0})$ with distance s.

We start with a set of initial energies $E_0^{(i)}$ and pitch angle cosines ${\mu }_0^{(i)}$ for i = 0…(N − 1) at s = 0. The energies are evenly distributed in log space with a spacing Δlog₁₀E, and the pitch angle cosines all have the same value $\overline{{\mu }_0}$. We then evaluate ${\overline{F}}_0(E_0,\overline{{\mu }_0})$ (Eq. (3) without the Dirac delta function) to obtain a corresponding set of injected flux values, ${\overline{F}}_0^{(i)}={\overline{F}}_0(E_0^{(i)},\overline{{\mu }_0})$. Using Eqs. (34)–(36), we advance $E_0^{(i)}$, ${\mu }_0^{(i)}$, and ${\overline{F}}_0^{(i)}$ a small distance Δs (with the same sign as $\overline{{\mu }_0}$) using a third-order Runge–Kutta method, to obtain their values ${\stackrel{\sim }{E_1}}^{(i)}$, ${\stackrel{\sim }{{\mu }_1}}^{(i)}$, and ${\stackrel{\sim }{F_1}}^{(i)}$ at s₁ = Δs. After the step, some of the lowest energies may have reached zero, meaning that the corresponding particles have been thermalised. If we kept advancing the remaining particles in the same manner, we would soon be left with too few non-zero energies to represent the distribution properly. Moreover, because the derivatives are larger in magnitude for smaller E and μ, the spacing between the lower energies would increase rapidly and amplify the undersampling problem.

To avoid this, we apply a re-meshing procedure after each step. We calculate a new set of energies $E_1^{(i)}={10}^{m\mathrm{\Delta }{log}_{10}E}{E}_0^{(i)}$ for i = 0…(N − 1), where m is an integer that may be positive or negative. The new energies are equal to the previous energies $E_0^{(i)}$ except shifted up or down in log space by a whole number m of the interval Δlog₁₀E. We set m so that the first re-meshed energy $E_1^{(0)}$ is as small as possible while still exceeding the smallest non-zero advanced energy ${\stackrel{\sim }{E_1}}^{(i)}$. In this way, we ensure good sampling coverage even if the distribution shifts significantly in energy. We then calculate the piecewise linear functions $\stackrel{\sim }{\mu }(E)$ and $\stackrel{\sim }{F}(E)$ that yield $\stackrel{\sim }{\mu }({\stackrel{\sim }{E_1}}^{(i)})={\stackrel{\sim }{{\mu }_1}}^{(i)}$ and $\stackrel{\sim }{F}({\stackrel{\sim }{E_1}}^{(i)})={\stackrel{\sim }{F_1}}^{(i)}$ for all non-zero ${\stackrel{\sim }{E_1}}^{(i)}$. For all re-meshed energies $E_1^{(i)}$ not exceeding the highest advanced energy ${\stackrel{\sim }{E_1}}^{(n-1)}$, we use these interpolating functions to calculate ${\mu }_1^{(i)}=\stackrel{\sim }{\mu }(E_1^{(i)})$ and ${\overline{F}}_1^{(i)}=\stackrel{\sim }{F}(E_1^{(i)})$.

For the re-meshed energies exceeding ${\stackrel{\sim }{E_1}}^{(N-1)}$, we have no samples of the flux spectrum and thus can no longer rely on interpolation. Instead, by selecting a sufficiently high upper limit for the initial energies, we can ensure that the highest advanced energy ${\stackrel{\sim }{E_1}}^{(N-1)}$ is large enough that the following high-energy limit versions of Eqs. (34)–(36) are valid:

$$\begin{array}{cc}& {\left(\frac{\mathrm{d}E}{\mathrm{d}s}\right)}_{\mathrm{high}}=-e\mathcal{E},\hfill \end{array}$$(51)

$$\begin{array}{cc}& {\left(\frac{\mathrm{d}\mu }{\mathrm{d}s}\right)}_{\mathrm{high}}=-\frac{1-{\mu }^2}{2\mu }\frac{\mathrm{d}lnB}{\mathrm{d}s},\hfill \end{array}$$(52)

$$\begin{array}{cc}& {\left(\frac{\mathrm{d}F}{\mathrm{d}s}\right)}_{\mathrm{high}}=-\frac{1-{\mu }^2}{2{\mu }^2}\frac{\mathrm{d}lnB}{\mathrm{d}s}F.\hfill \end{array}$$(53)

In the high-energy limit, the influence of collisions is negligible. Only the field-aligned electric field affects the energy of a particle, and only the magnetic gradient force affects the pitch angle and flux. From Eq. (51), the initial energy of a high-energy electron with energy E at the distance s is given by

$$\begin{array}{c}\hfill E_{\mathrm{high},0}(s,E)=E+{\displaystyle {\int }_0^{s}e\mathcal{E}(s^{\prime })\mathrm{d}{s}^{\prime }.}\end{array}$$(54)

Integrating Eq. (52), the pitch angle cosine at the distance s for a high-energy electron with initial pitch angle cosine μ₀ becomes

$$\begin{array}{c}\hfill {\mu }_{\mathrm{high}}(s,{\mu }_0)=\mathrm{sgn}\left({\mu }_0\right)\sqrt{1-(1-{{\mu }_0}^2)exp\left({\int }_0^s\frac{\mathrm{d}lnB}{\mathrm{d}s}(s^{\prime })\mathrm{d}{s}^{\prime }\right)}.\end{array}$$(55)

Combining Eqs. (52) and (53) and integrating, we obtain the expression for the flux of a high-energy electron at the distance s given an initial flux F₀:

$$\begin{array}{c}\hfill F_{\mathrm{high}}(s,{\mu }_0,F_0)=\frac{{\mu }_{\mathrm{high}}(s,{\mu }_0)}{{\mu }_0}{F}_0.\end{array}$$(56)

For all $E_1^{(i)}$ above the highest advanced energy ${\stackrel{\sim }{E_1}}^{(N-1)}$, we can obtain the value of ${\mu }_1^{(i)}$ from Eq. (55). To calculate ${\overline{F}}_1^{(i)}$, we evaluate Eq. (54) to find the initial energy of the electron, sample the initial flux spectrum ${\overline{F}}_0$ at this energy and insert the resulting initial flux into Eq. (56):

$$\begin{array}{cc}& {\mu }_1^{(i)}={\mu }_{\mathrm{high}}(s_1,\overline{{\mu }_0}),\hfill \end{array}$$(57)

$$\begin{array}{cc}& {\overline{F}}_1^{(i)}=\frac{{\mu }_1^{(i)}}{\overline{{\mu }_0}}{\overline{F}}_0(E_{\mathrm{high},0}(s_1,E_1^{(i)}),\overline{{\mu }_0}).\hfill \end{array}$$(58)

When the re-meshing is complete, we advance $E_1^{(i)}$, ${\mu }_1^{(i)}$, and ${\overline{F}}_1^{(i)}$ to obtain ${\stackrel{\sim }{E_2}}^{(i)}$, ${\stackrel{\sim }{{\mu }_2}}^{(i)}$, and ${\stackrel{\sim }{F_2}}^{(i)}$ at s₂ = 2Δs, re-mesh these to $E_2^{(i)}$, ${\mu }_2^{(i)}$, and ${\overline{F}}_2^{(i)}$, and repeat the procedure. We update the values of the integrals in Eqs. (54) and (55) continually during propagation. With the sets of values $E_k^{(i)}$ and ${\overline{F}}_k^{(i)}$, we then have a discrete version of the flux spectrum $\overline{F}(s,E,\overline{{\mu }_0})$ at s_k = kΔs.

In addition to the flux spectrum $\overline{F}(s,E,\overline{{\mu }_0})$, we need to determine the Jacobian determinant $|\partial \mu (s,E,\overline{{\mu }_0})/\partial \overline{{\mu }_0}|$. It can be calculated analytically for simplified versions of Eqs. (34) and (35) (ignoring ionisation fraction variations, magnetic gradient forces, electric fields and the ambient temperature) by solving for E(s, E₀, μ₀) and μ(s, E₀, μ₀), combining these to obtain μ(s, E, μ₀) and computing its partial derivative with respect to μ₀. Unfortunately, this is not possible for the more general case. Instead, we compute the Jacobian numerically by evolving a second set of electrons with the same initial energies but with the initial pitch angle cosine ${\overline{{\mu }_0}}^{\prime }=(1-\epsilon )\overline{{\mu }_0}$ reduced by a small fraction ε¹. So at any given depth s_k, we have the pitch angle cosines ${\mu }_k^{\prime (i)}$ for the “perturbed” electrons with ${\mu }_0^{(i)}={\overline{{\mu }_0}}^{\prime }$ in addition to the original solution values $E_k^{(i)}$, ${\mu }_k^{(i)}$, and $E_{0,k}^{(i)}$ for ${\mu }_0^{(i)}=\overline{{\mu }_0}$. We note that we advance the perturbed electrons in tandem with those in the original distribution and use the same value for m when re-meshing them. Consequently, after re-meshing, the perturbed electrons always have the same energies $E_k^{(i)}$ as the original electrons. We then compute the Jacobian determinant as follows:

$$\begin{array}{c}\hfill |\frac{\partial \mu }{\partial \overline{{\mu }_0}}{|}_k^{(i)}=|\frac{{\mu }_k^{(i)}-{\mu }_k^{\prime (i)}}{\overline{{\mu }_0}-{\overline{{\mu }_0}}^{\prime }}|.\end{array}$$(59)

At each distance s_k, we use the evolved electron properties $E_k^{(i)}$, ${\mu }_k^{(i)}$, and $E_{0,k}^{(i)}$ and flux spectrum ${\overline{F}}_k^{(i)}$ along with the Jacobian determinant ${|\partial \mu /\partial \overline{{\mu }_0}|}_k^{(i)}$ to numerically evaluate Eq. (43) for the total field-aligned electron flux F_beam(s_k) and Eq. (49) for the heating power density Q(s_k). We use F_beam(s_k) to calculate the return current resistive heating Q_r(s_k) (Eq. (50)) and to estimate the return current electric field ℰ_beam(s_k + 1) (Eq. (8)) for the following distance s_k + 1. We also use the total flux to decide when the beam has lost enough energy that we can terminate propagation.

3. Results

The transport model presented in this paper, hereafter referred to as the continuity equation characteristics (CEC) model, can account for a variety of additional physical effects that are ignored in the analytical transport model employed in Paper I (originally from Emslie 1978; Hawley & Fisher 1994). While significantly more computationally intensive than the analytical model, the CEC model remains sufficiently lightweight to be applied on the scale of millions of beams. We can thus compare it directly with the analytical model in the same atmospheric simulation. Our main objective with the results presented in this paper is to highlight the changes in the spatial distribution of beam energy deposition Q resulting from accounting for the additional physical effects.

To verify the CEC model, we ran it under the same assumptions as the analytical model² and compared the resulting energy deposition. To match the CEC model with the analytical model, we neglected collisions with helium by setting the helium-to-hydrogen ratio r_He to zero, we ignored the ambient temperature by setting W_E, W_μ, and W_F to unity, we omitted electric fields and magnetic gradient forces by setting ℰ and dlnB/ds to zero, and we left out ambient electrons from other elements than hydrogen by setting the electron-to-hydrogen ratio r_e equal to the hydrogen ionisation fraction x_H. The results showed good agreement between the analytical and CEC models, indicating that the latter is sound.

In each of the following sections, we present the results from running the CEC model with one of the additional physical effects included and the others ignored. All runs use N = 300 energies to represent the electron flux spectrum $\overline{F}$ and a power-law index of δ = 4 for the injected flux spectrum. To avoid potential confusion, we do not distinguish whether beams travel along the positive or negative magnetic field direction in the results. Hence, μ and s can always be considered positive.

3.1. With magnetic gradient forces

The run that included magnetic gradient forces, where the dlnB/ds factor in Eq. (35) was allowed to be non-zero, revealed some significant changes to the energy deposition in the atmospheric simulation compared to the results from the analytical model. Figure 2 shows the net beam heating power accumulated over the y-axis of the simulation box for the run with magnetic gradient forces. This figure is analogous to Fig. 10 in Paper I, which shows the corresponding result produced with the analytical model. In both cases, regions of particle acceleration, indicated by negative power (blue colour) due to the absorption of reconnection energy by new beams, appear along the interfaces of misaligned coronal loops and near loop footpoints directly above the transition region (TR). No difference in the distribution of acceleration regions is to be expected between the models, as the treatment of particle acceleration is identical. The distribution of deposited energy (orange colour) is broadly similar for both models, with the most intense beam heating occurring close to the boundaries of the acceleration regions and in the TR. However, for the model with magnetic gradient forces, the energy deposition in the corona is somewhat higher. At the same time, the beams do not penetrate as deep into the chromosphere, stopping abruptly around 500 km above the photosphere.

Fig. 2. Net change in heating power due to non-thermal electrons in the atmospheric simulation accumulated over the y-axis of the simulation box. Approximately 1.25 million beams were simulated using the Fokker–Planck-based model, with magnetic gradient forces accounted for in the electron transport equations. Blue regions have a net reduction in heating power resulting from released reconnection energy that would otherwise be released as local heating instead of going into particle acceleration. Orange regions indicate where the non-thermal energy is deposited as heat.

To enable a closer analysis of the differences made by magnetic gradient forces, we extracted three separate sets of beams demonstrating different energy transport scenarios. Figure 3 displays the three sets of beams. The figure is akin to Fig. 2, but contains only the net beam heating power due to the selected sets of beams. In set 1, the electrons are accelerated at the top of a coronal loop and follow one of the loop legs down into the TR. In set 2, acceleration occurs at various locations in the corona, but the resulting bundles of beams converge at the same location in the TR. In set 3, the acceleration region is a strong current sheet situated just above the TR, and the non-thermal electrons are ejected downwards in a coherent bundle. These sets correspond to those selected for the results in Paper I. In the following, we analyse the energy deposition for representative beams in the three sets.

Fig. 3. Three selected sets of electron beams, plotted in the same manner as in Fig. 2. Sets 1, 2, and 3 consist of approximately 4700, 1300, and 2800 beams, respectively.

3.1.1. Coronal loop (set 1)

Due to the high coherency of the beams in set 1, they all have a very similar heating profile, so it suffices to look at one of them. A comparison of Q(s) between the runs with and without magnetic gradient forces for one of the beams is shown in Fig. 4. In the basic model, Q has a peak in the corona around s = 6.5 Mm and another, much stronger peak at s = 20.15 Mm in the TR. The location of the TR and chromosphere can be inferred from the mass density profile in the figure. When magnetic gradient forces are included, a second coronal peak appears around s = 15.5 Mm, followed by a steep drop in Q and then a peak in the TR and chromosphere that is smaller but approximately 100 km deeper than the peak in the basic model.

Fig. 4. Deposited power density Q as a function of propagation distance for a representative beam in set 1 in Fig. 3 for two separate simulations, one with magnetic gradient forces (points with colour varying between blue and red) and one without (green points). The colours for the former simulation indicate the local relative change in magnetic flux density B with distance s. Blue colours imply that the magnetic field decreases in strength in the direction of propagation, while red colours imply that the magnetic field increases in strength. The local plasma mass density is shown as a solid grey curve. The right panel with a yellow background gives a magnified view of the plot in the left panel around the lower atmosphere, over the distance range indicated by the yellow band in the left panel. The labelled vertical lines mark locations of interest along the trajectory: (a) is near the site of injection, (b) is at the first coronal peak in Q, (c) is at the second coronal peak present only in the simulation with magnetic gradient forces, and (d) is near the peaks in Q below the transition region. The beam has ℱbeam, 0 ≈ 1.4 × 104 erg s−1 cm−2 and Ec ≈ 5 keV.

It is informative to consider the microscopic physical mechanisms leading to these features in the heating profile. To aid the discussion, Fig. 5 displays the non-thermal electron number distribution dn/dE as a function of energy E for the two models at some selected distances. The injected distribution has the most electrons at the lower cut-off energy E_c ≈ 5 keV, about 100 times more than at 15 keV (see panel a in Fig. 5, where the distributions are very close to the injected ones). As the electrons begin to propagate, the least energetic electrons lose energy and parallel velocity rapidly, their rate of loss with distance increasing as their energy decreases (due to the 1/E factor in Eq. (34) and 1/E² factor in Eq. (35)) in a runaway deceleration. This produces a steady rise in Q with distance, peaking at the point where the abundant least energetic electrons have lost all their excess energy and joined the thermal distribution (panel b in Fig. 5).

Fig. 5. Distribution dn/dE of beam electrons over energy E at the four separate distances s (with a separate panel for each distance) indicated in Fig. 4, both for a simulation that ignores (solid green curve) and includes (dashed purple curve) magnetic gradient forces. Also included is the collisional rate of energy loss for a single electron, −(dE/dt)coll, from Eq. (48) (solid grey curve). Summing a distribution curve and energy loss rate curve would yield a log-space curve of the integrand for the deposited power density Qbeam (Eq. (47)).

The initially more energetic electrons experience a slower rate of energy and parallel velocity loss at first (contributing a small part to Q) but eventually slow down enough to go through the same rapid deceleration that leads to thermalisation. Because there are fewer of these more energetic electrons (the peak number densities decrease from panel b to c in Fig. 5), the total power deposited as they become thermal is smaller than that deposited when the least energetic electrons were thermalised, so Q now decreases with distance.

When the electrons enter a region of denser plasma, as in panel d in Fig. 5, the collision rate increases, leading to more rapid energy loss and hence an instant proportional increase in Q (due to the n_H factor in Eq. (48)). The increased collision rate is evident when comparing the rate of energy loss (dE/dt)_coll in panel d to the other panels. But with the higher rate of collisions, thermalisation depletes the distribution from the lowest-energy end more rapidly, causing a steeper decrease in Q with distance. This is the reason for the falloff of Q following the TR peak for the basic model in Fig. 4. Entering a region of lower plasma density would have the opposite effect, with less energy deposition and slower depletion of the distribution.

When the electrons experience a rising magnetic flux density as they propagate, the magnetic gradient force deflects the electron velocities, transforming some parallel velocity into transverse velocity and hence reducing their pitch angle cosine μ (a positive dlnB/ds gives a proportionally negative contribution to dμ/ds in Eq. (35)). A reduction in μ means less of the electron’s velocity is directed forwards, so it must endure more collisions to advance along the field line. This leads to an increased rate of energy loss with distance (hence the 1/μ factor in Eq. (34)).

Less energetic electrons tend to be more susceptible to magnetic deflection because their velocities typically are partially deflected already due to their more frequent collisions. As Eq. (35) shows, the magnetic contribution to dμ/ds increases with decreasing μ. Still, because the deflection rate is not directly dependent on electron energy, arbitrarily energetic electrons in the distribution can be affected if their velocities are not entirely parallel to the magnetic field. And once magnetic deflection begins, it proceeds equally rapidly regardless of energy. As a result, the magnetic gradient force causes even the more energetic electrons in the distribution to lose energy faster with distance, making them join the less energetic part of the distribution sooner. This positive contribution to the less energetic part of the distribution (evident from comparing the two models in panel c in Fig. 5), where most of the energy transfer to the ambient plasma takes place, leads to an increase in Q. This is the reason for the second coronal peak in Q for the model with magnetic gradient forces in Fig. 4.

The increase in Q is only temporary, as it is counteracted by the loss of electrons whose velocities become entirely transverse, μ = 0, and thus are left behind. This initially only happens for electrons whose energies are nearly lost already (as in panel c in Fig. 5, where the distribution affected by magnetic gradient forces vanishes below E ≈ 1 keV), but subsequently for electrons with more and more energy remaining (the distribution vanishes at increasingly large energies between panels c and d). The eventual result of the loss of increasingly energetic electrons is a steep drop in Q with distance, as seen after the second coronal peak for the model with magnetic gradient forces in Fig. 4.

After leaving behind the forward-propagating beam, the lost electrons keep accelerating in the opposite direction due to the magnetic gradient force and begin travelling back along the field line in the direction they came from. This phenomenon, known as magnetic mirroring, is not accounted for in the CEC model because the transport equations (Eqs. (34)–(36)), being parametrised by distance s, become singular for μ(s) = 0. In the model, the energy contained in reflected electrons is lost. For the beam we are considering, the fraction of the injected flux lost to reflected electrons is only a few percent.

An implication of the gradual reflection of increasingly energetic electrons caused by magnetic gradient forces can be seen when the distribution enters a region dense enough to absorb its remaining energy. Because the electrons remaining are relatively energetic (evident from the very high average energy in the distribution affected by magnetic gradient forces compared to the unaffected one in panel d in Fig. 5), they penetrate deeper into the region before becoming thermal. Consequently, the peak of deposited power occurs deeper than it would if the less energetic electrons filtered out by magnetic gradient forces remained in the distribution. This explains the deeper peak location in Q for the model with magnetic gradient forces in the right panel of Fig. 4.

3.1.2. Converging bundles (set 2)

Most of the beams in set 2 (including the majority of beams accelerated near the height of 5 Mm to the far right in Fig. 3) follow a field line downwards in a converging magnetic field and thus have a heating profile qualitatively similar to that in Fig. 4. Some beams, however, experience a diverging magnetic field (that is, a reduction in the magnetic flux density) in the first part of their trajectory. This applies to the beams accelerated between the heights of 10.5 and 12.5 Mm in Fig. 3, as the coronal loop leg they travel along expands with depth for the first 10 Mm, before contracting again along the last stretch to the lower atmosphere. Figure 6 shows the heating profile for one of the beams following such a trajectory.

Fig. 6. As in Fig. 4, but for an electron beam in set 2 starting at x = 7.1 Mm and a height of 12.1 Mm in Fig. 3. The beam has ℱbeam, 0 ≈ 1.1 × 104 erg s−1 cm−2 and Ec ≈ 5.4 keV.

For this beam, Q remains lower throughout the corona in the model with magnetic gradient forces than in the model without. This is caused by the initial reduction in magnetic flux density with distance. From Eq. (35), the negative value of dlnB/ds counteracts the decrease in μ due to collisions. All electrons, including the more energetic ones, thus retain a higher value of μ, causing them to lose energy more slowly with distance. Since the more energetic electrons join the lower-energy part of the distribution at a reduced pace, the population at low energies decreases, reducing the energy deposition rate Q. This is precisely the opposite of what happens in a converging magnetic field.

Although Q does exhibit a second coronal peak in the model with magnetic gradient forces due to the eventual convergence of the magnetic field, similarly to the beam in Fig. 4, the energy deposition is nearly always lower than in the basic model. The only exception is at the end of the peak in the TR, as the energetic electrons that reach this depth have been subjected to slightly fewer collisions than in the basic model due to the earlier diverging magnetic field. Despite this, a significant proportion of the injected energy is clearly never deposited for the beam affected by magnetic gradient forces; upon inspection, around 50%. The culprit is heavy magnetic mirroring in the last 7 Mm of the trajectory, where the magnetic convergence is strong. This is clear from Fig. 7, which shows the minimum energy of electrons remaining in the distribution as a function of distance. When the magnetic convergence is negative or slightly positive, the minimum energy is close to zero, meaning that electrons are being thermalised. After s = 16 Mm, when the magnetic convergence increases significantly, the minimum energy increases as the electrons reach μ = 0 and get reflected before they lose their energy to collisions. The energy of reflected electrons increases more and more rapidly with distance. The figure also shows with colour the number density dn/dE of electrons with the minimum energy, indicating the number of electrons reflected at each distance. It peaks shortly after the onset of mirroring when the abundant electrons with energies around 4 keV are reflected. The number of reflected electrons then decreases rapidly with distance as the distribution gets depleted of increasingly energetic electrons.

Fig. 7. Minimum energy E of unreflected electrons (μ > 0) remaining in the distribution as a function of distance s for the beam in Fig. 6. The colours of the points indicate the number density dn/dE of electrons with that energy. Included as a dark grey curve is the magnetic convergence dlnB/ds.

The main reason this beam is particularly strongly affected by magnetic mirroring is the initial diverging magnetic field. Collisions, being more frequent at lower energies, remove less energetic electrons more efficiently and thus flatten the distribution as it propagates. By increasing μ and thus decreasing the number of collisions per unit of distance, the diverging magnetic field enables the distribution to retain not only more of its total flux but also more of its original steepness. This means that more beam flux is concentrated in the least energetic electrons when magnetic mirroring kicks in, and these are the electrons that get reflected first. The result is that more of the energy injected into the beam is lost to reflected electrons compared to what would have been lost in the absence of a diverging magnetic field.

The other notable group of beams subjected to a diverging magnetic field in set 2 are those accelerated near the height of 2.5 Mm around x = 14 to 17 Mm in Fig. 3 and ejected upwards along magnetic field lines that diverge with height. Some travel the length of the periodic simulation domain along the near horizontal magnetic field at the top of the domain before re-entering a coronal loop in the same downward trajectory as the beam in Fig. 6. They produce tiny Q values in their long journey through the upper corona, particularly for the model accounting for magnetic gradient forces. Otherwise, their heating profiles exhibit the same features as seen in Fig. 6, and they typically lose around 25% of their injected energy to magnetic mirroring.

3.1.3. Current sheet leg (set 3)

The beams in set 3 are very different from those in the other two sets in that they originate in a relatively dense environment and are ejected almost directly into the TR from the acceleration region. Figure 8 shows a typical heating profile from set 3. In both models, the vast majority of energy deposition occurs near the beginning of the trajectory owing to the immediately high density. In the basic model, Q decreases monotonically with distance after the first 50 km. However, in the model with magnetic gradient forces, the decrease in Q is strongly counteracted once the magnetic convergence becomes sufficiently high. The mechanism is the same as that responsible for the second coronal peak in Q for the beams starting higher in the corona: a faster decrease in μ with distance for all electrons, leading to a slight steepening of the distribution and thus more efficient energy deposition. The beam is effectively depleted after about 1 Mm of travel with enhanced energy deposition. In the basic model, the beam continues with a steady decrease in Q.

Fig. 8. As in Fig. 4, but for an electron beam in set 3 in Fig. 3 that is ejected straight down into the transition region. The beam has ℱbeam, 0 ≈ 6.4 × 105 erg s−1 cm−2 and Ec ≈ 1.8 keV.

Some of the beams in set 3 travel upwards along a diverging magnetic field for a short distance before the field line turns sharply downwards, leading the beam down into the TR. The heating profile for such a case is shown in Fig. 9. In the basic model, the evolution of Q with distance is effectively the same as in Fig. 8. In the model with magnetic gradient forces, the initial strong decrease in magnetic flux density prevents much of the early collisional deflection of parallel velocity so that μ remains higher than in the basic model throughout the distribution. This reduces the subsequent efficiency of energy loss by collisions. The intense magnetic field convergence ensuing as soon as the beam trajectory turns downwards quickly reflects the least energetic electrons, which are relatively abundant as collisions have not significantly depleted them. Hence, the beam loses most of its electrons and about 50% of its total energy to magnetic mirroring over the first hundreds of kilometres following the switch from magnetic field line divergence to convergence. Near s = 400 km, the magnetic convergence starts to diminish. This reduces the rate of magnetic reflection and allows the least energetic electrons to lose more energy to collisions before being reflected. Consequently, the deposited power Q increases. After about 300 km, magnetic convergence intensifies again, and magnetic reflection rapidly depletes the remaining beam electrons.

Fig. 9. As in Fig. 4, but for an electron beam in set 3 in Fig. 3 whose trajectory initially is directed upwards. The beam has ℱbeam, 0 ≈ 1.7 × 107 erg s−1 cm−2 and Ec ≈ 2.7 keV.

3.2. With non-hydrogen electrons

We performed another run with the basic version of the CEC model configured to use an electron-to-hydrogen ratio r_e computed from the actual electron and hydrogen densities in the atmospheric simulation rather than using r_e = x_H. Since the electron density in the simulation includes electrons contributed from elements other than hydrogen, r_e exceeds the hydrogen ionisation fraction x_H in the parts of the atmosphere hot enough to ionise the heavier elements. In the hot corona, we have x_H = 1 and r_e ≈ 1.1. As a result, the beams in this run typically showed a slightly elevated and earlier peak in energy deposition in the corona compared to the beams in the analytical model (where r_e = x_H). For example, the coronal peak in Q for the beam in Fig. 4 was elevated by 5% and occurred 1 Mm earlier.

3.3. With helium collisions

We performed a run to investigate the impact of including collisions with ambient helium. The only difference between the parameters of this run and the run with the most basic CEC model was that we calculated the helium-to-hydrogen ratio r_He ≈ 0.085 from the hydrogen and helium mass fractions rather than setting it to zero. The low value of r_He suggests that the influence of helium collisions is minor compared to the influence of collisions with electrons and hydrogen. Our results confirm this, showing only a minor increase in peak coronal energy deposition due to the extra collisions with fully ionised helium and a tiny decrease in chromospheric penetration depth caused by the extra collisions with singly ionised and neutral helium.

3.4. With non-zero ambient temperature

To explore the effect of accounting for the non-zero temperature of the ambient plasma, we performed another run with the basic CEC model configured to employ the warm-target transport equations (evaluating W_E, W_μ, and W_F from Eqs. (27), (28), and (38)) instead of the cold-target equations. The most affected beams were those with a high temperature in the coronal part of their trajectory, such as the beam in Fig. 4. They experienced a modest decrease in Q throughout the corona (no more than 10%, typically less). This is caused by the suppressed energy loss rate and slower decrease in μ for beam electrons with energies approaching the mean thermal energy, as explained by the curves in Fig. 1.

3.5. With the ambient electric field

To investigate the influence of the parallel component of the ambient electric field, we performed a run with ℰ = ℰ_fluid, computed using the electric field E_fluid supplied by the MHD simulation. Otherwise, the model was configured in the same way as the most basic version of the CEC model. The results showed that most beams have two sections of notable parallel electric field along their trajectory; at the very beginning and at the end.

The parallel electric field at the beginning is caused by the same reconnection event that produced the beam. Depending on the details of the reconnection event and the direction of the beam, the resulting value of ℰ_fluid can be either positive or negative, with a magnitude typically between 10⁻⁹ and 10⁻⁸ statV cm⁻¹. The distance after which ℰ_fluid vanishes corresponds roughly to the extent of the corresponding blue acceleration region in Fig. 2, but naturally depends on how deep within the region the beam originates. For some beams, the parallel electric field near their origin is sufficiently strong to accelerate or decelerate the non-thermal electrons slightly. The result is a minor shift in energy deposition towards smaller (if decelerated) or larger (if accelerated) distances. This shift is typically no longer than 1 Mm.

The parallel electric field often present at the end of a beam trajectory is also associated with reconnection, this time in the TR and chromosphere. Here the resistivity is generally higher than in the corona, so reconnection is more pervasive. However, the resulting electric fields are typically not strong enough to significantly influence the beam evolution, which is heavily dominated by the local high collision rates.

3.6. With the return current

We also performed a run including only the effects of the return current on top of the most basic model. We thus used ℰ = ℰ_beam, with ℰ_beam computed from Eq. (8), and included the return current heating term Q_r from Eq. (50) in the total heating power Q. The results were practically identical to those obtained without the return current. In the corona, ℰ_beam is of order 10⁻¹⁴ statV cm⁻¹ for most beams. As the beams enter the TR, where the resistivity increases enormously from the extremely low coronal values, ℰ_beam increases by roughly three orders of magnitude. The value decreases rapidly with depth as the beams lose flux to collisions, except for the beams accelerated just above the TR (such as the beams in set 3), which retain their flux for longer. Even for the most affected beams, the return current electric field is too weak to alter the flux spectrum in any way, and the return current resistive heating Q_r is lower than the collisional heating Q_beam by at least five orders of magnitude.

4. Discussion and conclusions

Magnetic gradient forces are the most consequential phenomenon ignored in the analytical model of Paper I that can be accounted for in the CEC model. For beams travelling along coronal loops, the converging magnetic field in the lower part of the loop legs deflects the electrons, causing a peak in energy deposition followed by a substantial dip as electrons of increasingly high energy get reflected and leave the beam. When the beam enters the lower atmosphere, the resulting peak in energy deposition occurs slightly deeper since magnetic mirroring has filtered out all but the relatively energetic electrons. After this peak, energy deposition drops immediately to zero as the last electrons are reflected. For beams ejected directly into the TR from the acceleration region, the heating peaks immediately due to the high density and then decreases rapidly, but magnetic convergence can hamper this decrease for a distance until all electrons have been reflected. A diverging magnetic field in the initial part of the trajectory focuses the beam. This makes it less susceptible to collisions, reducing the collisional flattening of the distribution and immediate energy deposition. When the beam enters a converging magnetic field with a steeper distribution, more electrons and energy are lost to magnetic mirroring.

An increase in coronal energy deposition combined with reduced penetration depth in the lower atmosphere, as we see in our heating profiles with magnetic convergence, was reported by Chandrashekar & Emslie (1986; and later Emslie et al. 1992), who modified the analytical expression for Q derived in Emslie (1978) to account for a converging magnetic field with dlnB/ds ∝ n_H. In their case, however, the increase in energy deposition was more evenly distributed. It did not include the second coronal peak and associated dip evident in many of our heating profiles. This discrepancy is because dlnB/ds typically increases significantly more steeply with depth than n_H in the lower parts of the coronal loops in our atmospheric simulation. Magnetic gradient forces can thus assert their influence on the beam quite abruptly before it becomes heavily depleted by collisions, which is the cause of the second coronal peak.

Numerical simulations of electron beam transport in isolated coronal loops with varying magnetic field strength (e.g. Leach & Petrosian 1981; Emslie et al. 1992; Allred et al. 2015) typically assume a purely converging magnetic field along the beam trajectory. In our atmospheric simulation, however, beams commonly find themselves in a diverging magnetic field. Our results (specifically Sect. 3.1.2) indicate that even a weakly diverging magnetic field in the initial parts of the trajectory can focus the distribution enough to resist collisional flattening and make it more susceptible to magnetic mirroring when the field subsequently converges. This suggests that magnetic mirroring may be of greater importance for collective non-thermal energy transport in realistic 3D atmospheres than one would infer from most existing simulations based on isolated loops.

Although the CEC model does not account for the reflected electrons, we can reason how they would affect the energy deposition. For the beam in Fig. 6 (from set 2), which is particularly strongly affected by magnetic mirroring, about 45% of the total injected energy is deposited in the corona, 5% is deposited in the lower atmosphere, and 50% is lost to magnetic mirroring. Figure 7 suggests that most of the mirrored energy is contained in electrons with relatively low energies reflected several megametres above the TR. Most of the energy lost to magnetic mirroring would thus be deposited in the corona. Consequently, this beam’s total coronal energy deposition would roughly double if reflected electrons were accounted for. A similar analysis for the beam in Fig. 4 (from set 1), where magnetic mirroring is less prominent, suggests a more modest increase of about 5% in coronal energy deposition due to reflected electrons. For beams ejected directly into the TR from the acceleration region (such as the beams in set 3), the electrons are reflected in a relatively dense environment, and most of the reflected energy would likely end up close to the acceleration region. Overall, the total fraction of injected energy lost to magnetic mirroring in our simulations is not very high: about 15% in beam set 1, 5% in set 2, 10% in set 3, and 5% in the atmosphere at large.

Likely, the most prominent changes the inclusion of magnetic gradient forces would make to synthetic observables from the electron beam simulations (such as those presented in Bakke et al. 2023) would be in the thermal and non-thermal emission produced at the locations of peak energy deposition in the TR and chromosphere. While the peak deposited energy is higher without magnetic convergence, the peak typically occurs a couple hundred kilometres deeper when magnetic convergence influences the beam. This extra depth could easily mean an order of magnitude higher density and a significantly lower ambient temperature, which would affect both the observed intensity and spectral signature.

The other physical effects supported by the CEC model that were neglected in the analytical model all had a much lower impact on the energy deposition than magnetic gradient forces. Coronal energy deposition was slightly increased by including electrons from ionised elements heavier than hydrogen or collisions with ambient helium and modestly reduced by accounting for a non-zero ambient temperature. Incorporating the parallel ambient electric field shifted energy deposition slightly towards smaller or larger distances for some beams. Modelling resistive heating from the return current and energy loss to the electric field driving it had no visible effect on the results.

The insignificance of the return current in our results is not surprising. As we argued in Paper I, the fluxes of the non-thermal electron beams accelerated in our atmospheric simulation are too small to produce significant charge separation and the associated return current. Nevertheless, being able to account for this effect opens up for applying the model to much more energetic beams where the return current plays an important role. To assess when the return current becomes significant, we artificially increased the energy flux ℱ_beam, 0 injected into the beam in Fig. 4. For this beam, ℱ_beam, 0 is originally of order 10¹² erg s⁻¹ cm⁻². After increasing this flux by five orders of magnitude, the influence of the return current started to become visible as an increase in Q along the first stretch of the trajectory and a decrease along the remainder.

From the limited difference using the warm-target transport equations made to our results, it is clear that the relatively high coronal temperatures of up to 2 million kelvin present in our atmospheric simulation are not sufficient to considerably influence the mean rates of change in energy and pitch angle for the non-thermal electrons. Even when we artificially increased the temperature along the coronal part of the trajectory of the beam in Fig. 4 from 2 to 10 million kelvin, the resulting reduction in coronal energy deposition remained minor. While the potential importance of the ambient temperature for non-thermal electron transport has been clearly demonstrated (Galloway et al. 2005; Jeffrey 2014; Kontar et al. 2015), our results attest that this mainly comes into play in more accurate transport models where collisional diffusion of energy and pitch angle is included.

The primary weakness of the non-thermal electron transport model presented here is that it neglects velocity randomisation by setting the collisional diffusion coefficients C_E² and C_μ² in Eq. (4) to zero. While the randomisation of energy only affects the evolution of the non-thermal distribution for energies approaching the mean thermal energy (see Eq. (A.60)), pitch angle diffusion becomes significant as soon as |μ| deviates appreciably from unity (see Eq. (A.61), as well as Bian et al. 2016). Emslie et al. (2018) show that the most prominent change to the energy deposition caused by accounting for pitch angle diffusion is the disappearance of the peak associated with the almost simultaneous thermalisation of electrons with energies near the lower cut-off energy E_c (visible, for example, at distance (b) in Fig. 4). Instead, these beam electrons scatter to different pitch angles, many being prevented from propagating far due to low values of μ. Hence, the highest energy deposition becomes concentrated near the point of injection. In addition, it seems plausible that by spreading out the electrons with a given energy across a range of different pitch angles, pitch angle diffusion would somewhat dampen the features associated with magnetic gradient forces apparent in the heating profiles presented here.

For the results presented in this paper, we used a power-law index of δ = 4. This gives a relatively flat injected distribution that is probably not representative of the weak beam type present in our simulation (Lin et al. 2001; Krucker et al. 2002). We did this to make the beams somewhat more resistant to collisions, making it easier to isolate and present the non-collisional influences on beam transport emphasised in this paper. This also makes the presented results more relevant for understanding more energetic beams. That being said, we also performed runs with higher values of δ. With the exception that a higher proportion of the injected energy is deposited in the corona for higher δ (Paper I covers the effect of varying δ in the analytical model), the results were qualitatively similar to those obtained with δ = 4.

Unlike the analytical transport model used in Paper I, we did not embed the CEC model directly into the Bifrost code. Instead, we applied the model on snapshots of the simulated atmosphere outputted from Bifrost, using the standalone Rust-based tool Backstaff³. As discussed in Frogner & Gudiksen (2022), there are significant challenges in efficiently implementing global energy transport in the domain decomposition-based parallelisation scheme used in Bifrost. These difficulties are magnified with the CEC model due to the considerable amount of data associated with each beam that must be communicated between processes. While the data usage of the model could likely be optimised, this is outside the scope of this paper.

Because the electron transport simulations are run after the fact, we are currently unable to model the magnetohydrodynamic response of the ambient plasma to the non-thermal energy transport computed with the CEC model. However, our results in this paper show that the analytical transport model running inside Bifrost does a decent job of matching the energy deposition computed with the CEC model, with the only significant deviations coming from the absence of magnetic gradient forces in the analytical model.

The fact that the CEC model represents the electron flux spectrum explicitly opens up applications that the analytical model does not support. For example, the injected flux spectrum does not have to be a simple power-law but can take on an arbitrary distribution over energy. Hence, the model could be combined with a more sophisticated treatment of particle acceleration than the simple parametric acceleration model employed here. Another example is the computation of non-thermal bremsstrahlung spectra for comparison with observations.

The CEC model is a major step towards bridging the gap between the basic analytical beam propagation model used in Paper I and the state-of-the-art models based on the direct numerical solution of the Fokker–Planck equation. Our starting point for the CEC model is a fairly comprehensive version of the Fokker–Planck equation, including all relevant contributions from the Lorentz force as well as a detailed treatment of collisions. By neglecting diffusion in velocity space, we can convert the Fokker–Planck equation into a set of ordinary differential equations for the mean evolution of beam electrons. We can solve these efficiently enough to simulate millions of beams. This enables us to model non-thermal electrons in a realistic 3D atmosphere with their spatial distribution, energetics, and trajectories emerging from the atmospheric simulation rather than being prescribed ad hoc. Thus, we sacrifice some realism in the physical modelling of electron transport (no velocity randomisation) for increased realism in the context in which this modelling is applied. The trajectories that beams follow in our simulation and the conditions along them are highly diverse. As we have shown, the diversity in conditions leads to equal diversity in the behaviour of the electron beams. This sometimes exposes interesting phenomena (such as the amplification of magnetic mirroring by a preceding magnetic field line divergence) that are easily overlooked in electron beam modelling with a more manually prescribed atmospheric structure along the trajectory.

Acknowledgments

We thank Dr. Gordon Emslie for providing helpful clarifications on mean scattering theory. This research was supported by the Research Council of Norway through its Centres of Excellence scheme, project number 262622, and through grants of computing time from the Programme for Supercomputing.

References

Appendix A: Specialising the Fokker–Planck equation

Let the electron distribution f(r, v, t) be defined such that f(r, v, t) d³r d³v is the number of electrons within the volume element d³r centred on position r with velocities within d³v of velocity v at time t. The evolution of f(r, v, t) is governed by the Fokker–Planck equation,

$$\begin{array}{c}\hfill \frac{\partial f(\mathbf{r},\mathbf{v},t)}{\partial t}+\mathbf{v}·\mathrm{\nabla }f(\mathbf{r},\mathbf{v},t)+{\left(\frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t}\right)}_{!\mathrm{C}}·{\mathrm{\nabla }}_{v}f(\mathbf{r},\mathbf{v},t)={\left(\frac{\partial f(\mathbf{r},\mathbf{v},t)}{\partial t}\right)}_{\mathrm{coll}},\end{array}$$(A.1)

where (dv/dt)_!C is the acceleration of the electrons due to non-collisional forces, and (∂f/∂t)_coll is the rate of change in the electron distribution due to collisions with ambient particles. We assume that the distribution reaches a steady state quickly compared to the rate of change in the background plasma so that we can ignore the explicit time dependence of the distribution:

$$\begin{array}{c}\hfill \frac{\partial f(\mathbf{r},\mathbf{v},t)}{\partial t}=0.\end{array}$$(A.2)

Since the electrons follow gyrating trajectories along a magnetic field line, it is helpful to express the position r in terms of the coordinates (s, r_⊥, ϕ), where (r_⊥, ϕ) are the polar coordinates of the electron within the plane perpendicular to the magnetic field vector B at a distance s along the field line. We then have d³r = r_⊥ dϕ dr_⊥ ds =  dA ds. In the case of gyromotion, it is reasonable to assume azimuthal symmetry, implying that no quantities in Eq. (A.1) depend on ϕ. If we furthermore assume that the gyrating trajectories are relatively tight so that r_⊥ remains very small compared to the scale of distances s, we can ignore the dependence of the distribution on r_⊥. We can then consider f(s, v) constant over its local cross-sectional area A(s) and write the second term in Eq. (A.1) as

$$\begin{array}{c}\hfill \mathbf{v}·\mathrm{\nabla }f(\mathbf{r},\mathbf{v})=\frac{\mathrm{d}s}{\mathrm{d}t}\frac{\partial f(s,\mathbf{v})}{\partial s}.\end{array}$$(A.3)

The right-hand side in the Fokker–Planck equation is given by (omitting for brevity the dependency of f on s)

$$\begin{array}{c}\hfill {\left(\frac{\partial f(\mathbf{v})}{\partial t}\right)}_{\mathrm{coll}}=-{\displaystyle \sum _i\frac{\partial }{\partial v_i}\left(\langle \mathrm{\Delta }{v}_i\rangle f(\mathbf{v})\right)+\frac{1}{2}\sum _{\mathit{ij}}\frac{{\partial }^2}{\partial v_i\partial v_j}\left(\langle \mathrm{\Delta }{v}_i\mathrm{\Delta }{v}_j\rangle f(\mathbf{v})\right),}\end{array}$$(A.4)

where

$$\begin{array}{c}\hfill \langle \mathrm{\Delta }{v}_i\rangle ={\displaystyle \sum _{\alpha }{\langle \mathrm{\Delta }{v}_i\rangle }_{\alpha }\mathrm{and}\langle \mathrm{\Delta }{v}_i\mathrm{\Delta }{v}_j\rangle =\sum _{\alpha }{\langle \mathrm{\Delta }{v}_i\mathrm{\Delta }{v}_j\rangle }_{\alpha }.}\end{array}$$(A.5)

Here ⟨χ⟩_α is the mean rate of a quantity χ due to collisions with ambient particles of species α, and is calculated as

$$\begin{array}{c}\hfill {\langle \chi \rangle }_{\alpha }(\mathbf{v})={\displaystyle {\int }_{v_{\alpha }}{\int }_{\mathrm{\Omega }}\chi \left||\mathbf{v}-{\mathbf{v}}_{\alpha }\right||f_{\alpha }({\mathbf{v}}_{\alpha }){\left(\frac{\mathrm{d}\sigma }{\mathrm{d}\mathrm{\Omega }}\right)}_{\alpha }(\mathbf{v},{\mathbf{v}}_{\alpha },\mathrm{\Omega })\mathrm{d}\mathrm{\Omega }{\mathrm{d}}^3{v}_{\alpha },}\end{array}$$(A.6)

where (dσ/dΩ)_α is the differential scattering cross-section for a collision between a beam particle and an ambient particle of species α with velocity v_α, and f_α is the distribution of the species α ambient particles. In Eq. (A.4), the quantity Δv_i for which the mean rate is computed is the i’th component of the change in velocity v due to a collision.

We can write Eq. (A.5) as

$$\begin{array}{c}\hfill \langle \mathrm{\Delta }{v}_i\rangle ={\displaystyle \sum _{\mathrm{c}}{\langle \mathrm{\Delta }{v}_i\rangle }_{\mathrm{c}}+\sum _{\mathrm{N}}{\langle \mathrm{\Delta }{v}_i\rangle }_{\mathrm{N}}\mathrm{and}\langle \mathrm{\Delta }{v}_i\mathrm{\Delta }{v}_j\rangle =\sum _{\mathrm{c}}{\langle \mathrm{\Delta }{v}_i\mathrm{\Delta }{v}_j\rangle }_{\mathrm{c}}+\sum _{\mathrm{N}}{\langle \mathrm{\Delta }{v}_i\mathrm{\Delta }{v}_j\rangle }_{\mathrm{N}},}\end{array}$$(A.7)

where c represents charged particle species, and N represents neutral particle species. By inserting the Rutherford differential scattering cross-section into Eq. (A.6), Rosenbluth et al. (1957) derive expressions for ⟨Δv_i⟩_c and ⟨Δv_iΔv_j⟩_c:

$$\begin{array}{cc}\hfill {\langle \mathrm{\Delta }{v}_i\rangle }_{\mathrm{c}}& =-4\pi (1+\frac{m_{\mathrm{e}}}{m_{\mathrm{c}}})\frac{{\mathrm{\Gamma }}_{\mathrm{c}}}{{m_{\mathrm{e}}}^2}\frac{\partial {\varphi }_{\mathrm{c}}(\mathbf{v})}{\partial v_i},\hfill \end{array}$$(A.8)

$$\begin{array}{cc}\hfill {\langle \mathrm{\Delta }{v}_i\mathrm{\Delta }{v}_j\rangle }_{\mathrm{c}}& =-8\pi \frac{{\mathrm{\Gamma }}_{\mathrm{c}}}{{m_{\mathrm{e}}}^2}\frac{{\partial }^2{\psi }_{\mathrm{c}}(\mathbf{v})}{\partial v_i\partial v_j},\hfill \end{array}$$(A.9)

where m_e is the electron mass, m_c is the mass of a species c particle, and

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{c}}=4\pi e^4{z_{\mathrm{c}}}^{2}ln{\mathrm{\Lambda }}_{\mathrm{c}}.\end{array}$$(A.10)

Here e is the elementary charge, z_c is the charge number, and lnΛ_c is the Coulomb logarithm. The functions ϕ_c and ψ_c, defined as

$$\begin{array}{cc}\hfill {\varphi }_{\mathrm{c}}(\mathbf{v})& =-\frac{1}{4\pi }{\displaystyle {\int }_{v_{\mathrm{c}}}\frac{f_{\mathrm{c}}({\mathbf{v}}_{\mathrm{c}})}{\left||\mathbf{v}-{\mathbf{v}}_{\mathrm{c}}\right||}{\mathrm{d}}^3{v}_{\mathrm{c}},}\hfill \end{array}$$(A.11)

$$\begin{array}{cc}\hfill {\psi }_{\mathrm{c}}(\mathbf{v})& =-\frac{1}{8\pi }{\displaystyle {\int }_{v_{\mathrm{c}}}\left||\mathbf{v}-{\mathbf{v}}_{\mathrm{c}}\right||f_{\mathrm{c}}({\mathbf{v}}_{\mathrm{c}}){\mathrm{d}}^3{v}_{\mathrm{c}},}\hfill \end{array}$$(A.12)

are potentials satisfying Poisson’s equations:

$$\begin{array}{cc}\hfill {{\mathrm{\nabla }}_v}^2{\varphi }_{\mathrm{c}}(\mathbf{v})& =f_{\mathrm{c}}(\mathbf{v}),\hfill \end{array}$$(A.13)

$$\begin{array}{cc}\hfill {{\mathrm{\nabla }}_v}^2{\psi }_{\mathrm{c}}(\mathbf{v})& ={\varphi }_{\mathrm{c}}(\mathbf{v}).\hfill \end{array}$$(A.14)

Disregarding, for now, the contribution of collisions with neutral particles, we can apply Eqs. (A.8) and (A.9) so that Eq. (A.4) becomes

$$\begin{array}{c}\hfill {\left(\frac{\partial f(\mathbf{v})}{\partial t}\right)}_{\mathrm{coll}}={\mathrm{\nabla }}_v·{\displaystyle (\left(\sum _{\mathrm{c}}\frac{4\pi {\mathrm{\Gamma }}_{\mathrm{c}}}{{m_{\mathrm{e}}}^2}(1+\frac{m_{\mathrm{e}}}{m_{\mathrm{c}}}){\mathrm{\nabla }}_v{\varphi }_{\mathrm{c}}(\mathbf{v})\right)f(\mathbf{v})-{\mathrm{\nabla }}_v·\left(f(\mathbf{v})\sum _{\mathrm{c}}\frac{4\pi {\mathrm{\Gamma }}_{\mathrm{c}}}{{m_{\mathrm{e}}}^2}{\mathrm{\nabla }}_v{\mathrm{\nabla }}_v{\psi }_{\mathrm{c}}(\mathbf{v})\right)),}\end{array}$$(A.15)