Solar flares are the most powerful energy-release phenomena on the Sun (e.g., Fletcher et al., 2011; Benz, 2017). A large number of charged particles are accelerated to high energies, including electrons, protons, and heavy ions, which can further excite nonthermal emissions from radio to gamma-rays via different radiation mechanisms. Although the primary acceleration mechanism remains unclear (Miller et al., 1997; Zharkova et al., 2011; Kong et al., 2019; Li et al., 2021), it is generally believed that electrons are accelerated in the coronal region above flare loops. After being injected at the top of flare loops, accelerated electrons travel to the loop footpoints and deposit energy in the high-density chromosphere, resulting in chromospheric evaporation. In some strong flares, energetic particles can deliver energy to the deeper atmosphere and have impact on the photosphere, suggested as the driver of white-light flares and sunquakes (e.g., Wu et al., 2023). Therefore, the acceleration and transport of energetic electrons plays a central role in the solar flare dynamics.

Energetic electrons are not free-streaming and subject to various effects during their transport from the looptop to the footpoints. The transport effects include magnetic mirroring due to the convergence in magnetic field, pitch-angle scattering by magnetic turbulence, energy loss and pitch-angle scattering via Coulomb collisions, return current, etc (e.g., Fletcher and Martens, 1998; Minoshima et al., 2011; Jeffrey et al., 2014; Kontar et al., 2014; Bian et al., 2017; Effenberger and Petrosian, 2018; Musset et al., 2018; Allred et al., 2020; Tang et al., 2020; Kong et al., 2022). Therefore, the pitch-angle distribution of energetic electrons is time dependent and should not be isotropic. Anisotropy in the velocity space is of critical importance to nonthermal emissions. It plays a fundamental role in plasma wave excitation in coherent emission mechanisms of solar radio bursts (Melrose, 2017). For example, electron-cyclotron maser emission requires a positive gradient of perpendicular direction, such as loss cone and horseshoe distributions (e.g., Melrose andWheatland, 2016; Zhao G. Q. et al., 2016; Ning et al., 2021a; Ning et al., 2021b; Tang et al., 2024). Anisotropic distribution can also affect the intensity, spectrum, and polarization of incoherent emissions, e.g., in microwave (e.g., Kuznetsov and Fleishman, 2021) and X-rays (e.g., Kuznetsov and Fleishman, 2021) and X-rays (e.g., Charikov et al., 2012; Melnikov et al., 2013).

Magnetic turbulence is an essential element both in particle acceleration (e.g., stochastic or shock acceleration) and transport processes in solar flares. Recent observations from nonthermal broadening of spectral lines by Hinode/EIS (e.g., Stores et al., 2021) revealed the presence of turbulence throughout the flare loop, although the strongest is at the looptop. In MHD simulations of magnetic reconnection in solar flares, the impact of reconnection outflows on the flare loop can trigger various instabilities and cause a highly turbulent plasma environment (e.g., Ruan et al., 2023; Wang et al., 2023; Ye et al., 2023). Recently, Effenberger and Petrosian (2018) studied the particle escape time for different initial pitch-angle distributions by solving the Fokker-Planck transport equation and assuming isotropic pitch-angle scattering by magnetic turbulence. Melnikov and Filatov (2020) investigated the conditions for the generation of whistler turbulence in the flare loop, which can resonate with energetic electrons and significantly affect their spectral and pitch-angle distributions (Melnikov and Filatov, 2021).

In this work, we numerically model the propagation of energetic electrons after being injected into the flare loop, and focus on the effects of magnetic mirror and turbulent scattering on the transport and anisotropic distribution of electrons. The paper is organized as follows. Section 2 describes our numerical model and Section 3 presents the simulation results. Summary and discussion are given in Section 4.

For the flare loop, we use an analytical two-dimensional magnetic field model in the x − y plane (Lin et al., 1995; Minoshima et al., 2010), A z x , y = − y + d x 2 + y + d 2 − y a + d 2 , B x x , y = ∂ A z ∂ y = − x 2 − y + d 2 x 2 + y + d 2 2 + 1 a + d 2 , B y x , y = − ∂ A z ∂ x = − 2 x y + d x 2 + y + d 2 2 , where A z is the flux function, B x and B y are two components of the magnetic field, a is the height of X-type neutral line at the top of flare loops, and d is the depth of the dipole below the photosphere. In Figures 1A, 2, the thin curves are contours of A z and illustrate the magnetic field lines of the flare loop model. Here we assume a = 80 Mm, d = 50 Mm.

Following our previous work Kong et al. (2022), we model the transport of energetic electrons in the flare loop by numerically solving the focused transport equation (Roelof, 1969; Skilling, 1971; van den Berg et al., 2020). The equation includes various transport effects, such as streaming along the magnetic field, advection with the solar wind, pitch-angle scattering, magnetic focusing/mirroring, and adiabatic cooling. Therefore, it has been widely applied to study the acceleration and transport of solar energetic particles (e.g., Qin et al., 2006; Zhang et al., 2009; Dröge et al., 2010; Zuo et al., 2011; Wang et al., 2012; Zhao L. et al., 2016; Hu et al., 2017; Zhang and Zhao, 2017; Wei et al., 2019; Wijsen et al., 2019). A similar Fokker-Planck transport equation has also been used in modeling energetic electrons in solar flares, in which the effects of magnetic mirroring, Coulomb collisions, and pitch-angle scattering are often included (e.g., Hamilton and Petrosian, 1990; Fletcher, 1995; Kontar et al., 2014; Effenberger and Petrosian, 2018; Melnikov and Filatov, 2021).

In this work, we focus on the effect of pitch-angle scattering on electron trapping/precipitation and the anisotropic distribution of energetic electrons. We neglect the advection term and the energy change due to Coulomb collisions, compression and shear in plasma flow (Kong et al., 2022). Test-particle simulations in synthetic turbulence suggested that the perpendicular diffusion coefficient is a few percent of the parallel diffusion coefficient (Giacalone and Jokipii, 1999). Cross-field diffusion may affect both the size and energy dependence of nonthermal emissions (Kontar et al., 2011), and the escape of electrons to the open field line. Here it is neglected for simplicity. The reduced particle transport equation can be written as (Roelof, 1969; Effenberger and Petrosian, 2018), ∂ f ∂ t = − v μ b ̂ ⋅ ∇ f − v 1 − μ 2 2 L B ∂ f ∂ μ + ∂ ∂ μ D μ μ ∂ f ∂ μ , where f is the distribution function of charged particles, v is the particle speed, μ is the pitch-angle cosine, and t is the time. The terms on the right-hand side describe the electron streaming along the direction of magnetic field b ̂ , the magnetic mirroring effect with the focusing length L B = ( b ̂ ⋅ ∇ l n B ) − 1 , and the pitch-angle diffusion with a coefficient D μ μ .

The pitch-angle diffusion coefficient D μ μ describes the resonant interaction between the particle and the turbulent magnetic field. In the quasi-linear theory, it is given by (Jokipii, 1971), D μ μ = π 4 Ω 0 1 − μ 2 k r P k r B 0 2 , where Ω 0 = q B 0 / m is the particle gyrofrequency with the mass m and the charge q , P ( k ) is the turbulence power spectrum, and k r = Ω 0 / ( v | μ | ) is the resonant wavenumber. We assume the form of Kolmogorov turbulence spectrum with the spectral index Γ = 5/3. In the non-relativistic limit, the pitch-angle diffusion coefficient can be expressed as (Beeck and Wibberenz, 1986), D μ μ = D μ μ 0 p p 0 Γ − 1 1 − μ 2 | μ | Γ − 1 + h 0 . D μ μ 0 is a constant describing the scattering rate and depends on the level of magnetic field fluctuation. p 0 is the particle momentum at the energy E 0 = 10 keV. The parameter h 0 is added to describe the finite scattering through μ = 0 and here we set h 0 = 0.05.

Because the transport equation is essentially a Fokker-Planck equation, it can be recast into a set of stochastic differential equations (SDEs) (e.g., Zhang, 1999; Strauss and Effenberger, 2017). Here we use the following time-forward SDEs to trace the particle’s position and pitch-angle (Kong et al., 2022). d X = v μ b ̂ d t , d μ = − v 1 − μ 2 2 L B + ∂ D μ μ ∂ μ d t + 2 D μ μ d W μ t , where d W μ is a Wiener process.

In the simulations, we assume that electrons have been accelerated near the top of the flare loop and only consider the transport process in the loop. Energetic electrons with a power-law energy spectrum, f ( E ) ∼ E − δ , are impulsively injected in the looptop region, given by x = [-2, 2] Mm and y = [48, 52] Mm. Here we set the electron energy spectral index δ = 3, and the energy range is between 0.7 and 153 keV (electron velocity between 0.05 c and 0.64 c , c is the speed of light). The initial electron pitch-angle distribution is assumed to be isotropic. In each simulation, a total of 6 million pseudo-particles are injected. We note that the injection of accelerated electrons is not necessarily at the top of the loop in a realistic solar flare. For example, when the reconnection takes place between a closed loop with other loops or an open field line. This may give rise to asymmetric distribution in space and anisotropy of energetic electrons.

To study the effect of turbulent scattering on electron transport and anisotropic distribution, we conduct five simulation runs with different levels of magnetic fluctuations by changing the value of D μ μ 0 . We take D μ μ 0 = 0 ( R u n A ) , 0.0272 s − 1 ( R u n B ) , 0.272 s − 1 ( R u n C ) , 2.72 s − 1 ( R u n D ) , and 27.2 s − 1 ( R u n E ) . Then, the time scale of turbulent scattering is approximately τ d = 1/ D μ μ 0 , varying between 36.8 s and 0.0368 s from R u n B to R u n E . Note that for electrons with the energy E 0 = 10 keV, v e 0 = 0.195 c = 5.85 × 1 0 7 m s − 1 , and a loop length of L 0 = 100 Mm, the crossing time scale in the loop where they are injected, τ c = L 0 / v e 0 = 1.71 s.

Three regimes of turbulent pitch-angle diffusion was defined in Bespalov et al. (1987), weak ( τ d > σ τ c ) , moderate ( τ c < τ d < σ τ c ) , and strong ( τ d < τ c ) , where σ is the mirror ratio of the flare loop. For the field lines where electrons are injected, the magnitudes of magnetic field in the looptop and at the footpoint are 40.8 G and 241 G, respectively. Then, the mirror ratio is σ = B F P / B L T = 5.9, and the critical pitch angle is θ c = a r c s i n 1 / σ = 24.3 ° . Therefore, R u n B corresponds to the weak diffusion regime, R u n C in the moderate diffusion regime, and R u n D and R u n E in the strong diffusion regime. We can also calculate the particle mean free path for 10 keV electrons, λ ‖ 0 = 3 κ ‖ 0 / v e 0 , where κ ‖ 0 is the spatial diffusion coefficient along the direction of the magnetic field and related to the pitch-angle diffusion coefficient D μ μ (Kong et al., 2022). Then, we can get λ ‖ 0 = 4,100 Mm in R u n B , much larger than the loop length; λ ‖ 0 = 410 Mm in R u n C , comparable to the loop length; λ ‖ 0 = 41 Mm and 4.1 Mm in R u n D and R u n E , smaller than the loop length.

To test the validity of the simulation, particularly the pitch-angle scattering through μ = 0, we first examine the trajectory of a single electron. Figure 1 shows the simulation results for three representative electrons, trajectories plotted over the magnetic field lines, variations of y position and μ as a function of time. The three electrons are injected in the same position at the looptop with the same energy of 10 keV and initial pitch-angle of 45 ° . For the electron without turbulent scattering (as in R u n A ), the curves are plotted in black. It is reflected at y R = 35.7 Mm, where the magnetic field strength is B R = 81.7 G. We then can get the critical pitch-angle at the reflection point is a r c s i n B L T / B R = 45 ° , consistent with the initial condition. For the two electrons with weak turbulent scattering (as in R u n B ), the curves are plotted in red and blue, respectively. Due to the pitch-angle diffusion, electrons can be scattered into the loss cone. Therefore, the two electrons can go deeper than the expected reflection position y R . For the electron plotted in blue, it is not reflected while moving to the left footpoint and finally precipitates into the solar surface. As seen from the evolution of μ in panel (c), the electrons can be scattered smoothly through μ = 0.

Figure 2 displays the spatial distributions of energetic electrons at three energy ranges, 5–10 keV, 20–30 keV, and > 50 keV, in R u n A . The simulation times in panels (a)-(c) are τ c = 1.71 s, 2 τ c = 3.42 s, and 3 τ c = 5.13 s, respectively. Due to the trapping effect of magnetic mirror, most electrons are concentrated around the top of the flare loop. Since the initial pitch-angle distribution of injected electrons is isotropic, electrons with larger pitch-angles take much more time as they move from the loop top to lower altitudes. Therefore, we can see multiple streams of electrons bouncing back and forth in the loop, as denoted by the black arrows (only the right side is marked). The number density of streaming electrons is smaller than that trapped at the looptop and the pattern varies with energy. As shown below, the streams of electrons are the reason for the presence of stripes in the velocity space distribution. For different simulation runs, the spatial distribution is generally similar. With increasing scattering rate, the distribution gets smoother and streaming electrons are harder to be distinguished.

We now analyze the effect of different scattering rates (as described by D μ μ 0 ) on the trapping and precipitation of electrons in the flare loop. Figure 3 shows the fractions of electrons trapped in the looptop and precipitating to the solar surface for various energies, 5 keV, 10 keV, 50 keV, and 100 keV, respectively, at the end of the simulation (3 τ c ). For each energy, the trapped fraction is defined as F t = N y > 45 / N inject , where N y > 45 is the number of electrons that remain trapped at y > 45 Mm in the simulation domain and N inject is the injected population. For the precipitating fraction, it is defined as F p = 1 − N loop / N inject , where N loop is the number of electrons that remain bouncing in the loop and have not reached the bottom boundary. As noted above, for the field line where electrons are initially injected, the mirror ratio σ = 5.9 and the critical pitch angle θ c = 24.3 ° . The electrons with pitch angle smaller than θ c fall into the loss cone and can escape. Therefore, it results in an expected precipitating fraction F p c = 27%, as denoted by the dashed line in Figure 3B. In the simulation of R u n A without turbulent scattering, the precipitating fractions at different energies agree well with the theoretical predication.

As shown in Figure 3, with increasing scattering rate, the variations of both the trapped fraction ( F t ) and precipitation fraction ( F p ) display a non-monotonic pattern. From the non-scattering case in R u n A to the weak and intermediate scattering cases in R u n B and R u n C , the trapped fraction decreases and the precipitation fraction increases. Due to pitch-angle scattering, more and more electrons with initial pitch-angle larger than the critical value θ c = 24.3 ° are scattered into the loss cone and escape. However, for the low-energy electrons of 5 keV, although the trapped fraction decreases, there is no rise in the precipitation fraction. It suggests that the electrons have not reached the solar surface while they have left the looptop, possibly due to their low speed. For the strong scattering cases in R u n D and R u n E , the pitch-angle scattering is so frequent that electrons should stay at the looptop for much longer time before moving to lower altitudes. Thus, from moderate to strong scattering, the trapped fraction increases and the precipitation fraction decreases. This indicates that the precipitation fraction is highest in the moderate diffusion regime, therefore, in favor of high intensity of nonthermal HXR and microwave emissions in the footpoints. In contrast, to reproduce a bright nonthermal source in the looptop, either weak or strong scattering is required. We also find that the magnitude of variation is energy dependent. The trapped fraction decreases at higher energies, while the precipitation fraction increases with energy.

Due to the effects of magnetic mirror and turbulent scattering, the particle distribution in the velocity space varies along the flare loop and with time. Figures 4, 5 show the velocity space distribution in the looptop and loopleg regions, respectively. v ‖ and v ⊥ are velocity components in the parallel and perpendicular directions. Here the looptop is integrated over x = [-5, 5] Mm and y = [45, 55] Mm, and the loopleg on the right side is integrated over x = [20, 40] Mm and y = [20, 30] Mm. The simulation results in R u n A , R u n B , and R u n E are displayed in panels (a), (b), and (c), respectively. The left and right columns are at two different times, τ c and 3 τ c , respectively. The dashed line in each panel illustrates the critical pitch angle θ c = 24.3 ° for the magnetic field lines where electrons are injected.

In the non-scattering ( R u n A ) and weak scattering ( R u n B ) runs, the velocity space distributions are obviously anisotropic, but exhibit different patterns in the looptop and loopleg regions. As shown in panel (a) in Figure 4, at the looptop, multiple narrow bands (which resemble branches or fishbone) stretch out from the vertical axis and present positive gradients in the perpendicular direction, i.e., ∂ f / ∂ v ⊥ > 0 . The number of bands increases with time and they gather towards the origin of the coordinate system. If we continue to run the simulation, the gap between stripes gets smaller and the distribution will evolve into a double-sided loss cone. Those electrons at the looptop are mainly reflected and trapped electrons, therefore most electrons are distributed in the perpendicular direction. Since electrons with larger v ‖ can leave the looptop faster or be reflected faster, multiple streams of electrons can be observed as shown in Figure 2, leading to fishbone-like multiple bands as time goes on. As shown in panel (b) in Figure 4, with weak scattering, similar stripes can be seen at the early time, which also exhibits ∂ f / ∂ v ⊥ > 0 . Due to turbulent scattering, the width of stripes increases and some electrons fill in the gaps between stripes. At later time as shown on the right, the gaps between stripes are nearly smoothed out and the distribution resemble a double-sided loss cone. As shown in panel (c), for the simulation with strong scattering ( R u n E ) , the distribution has already become nearly isotropic at the early time. In addition, compared with the non-scattering case in panel (a), electrons are scattered into the loss cone (below the dashed line) at various energies as a result of pitch-angle scattering, consistent with the results as shown in Figure 3.

As shown in Figure 5, in the loopleg region, the velocity space distributions look different from that in the looptop. In panels (a) and (b), at the early time, the first stripe (close to the origin of the coordinate system) is circular-shaped and represents the contribution from the beam-like electrons before getting reflected. The distribution resembles the so-called horseshoe distribution as observed in the source of auroral kilometric radiation (AKR, see, e.g., Ergun et al., 2000; Treumann, 2006). It contains positive gradients in both parallel and perpendicular directions, i.e., ∂ f / ∂ v ‖ > 0 and ∂ f / ∂ v ⊥ > 0 . At later time, after being reflected, v ‖ evolves from positive to negative values and the distribution presents mainly ∂ f / ∂ v ⊥ > 0 . Similarly, the other stripes appearing later also have ∂ f / ∂ v ⊥ > 0 . As in the looptop, if we continue to run the simulation, the distribution will eventually evolve into a double-sided loss cone, but it is asymmetric. For the case with strong scattering, as shown in panel (c), the distribution is nearly isotropic.

In this work, we numerically model the transport of energetic electrons in the flare loop after being injected around the top of the loop. We examine the effect of turbulent pitch-angle scattering on the trap/precipitation fraction and velocity space distribution by introducing different levels of scattering rates. We find that both the fractions of electrons trapped in the looptop and precipitating into the solar surface vary in a non-monotonic way with increasing scattering rate. From non-/weak to intermediate scattering, the trapped fraction decreases and the precipitation fraction increases, while from intermediate to strong scattering, the trapped fraction increases and the precipitation fraction decreases. Therefore, in the moderate diffusion regime, we expect intense nonthermal HXR and microwave emissions in the footpoints due to the highest precipitation fraction. In addition, the trap/precipitation fraction apparently shows energy dependence, which will affect the electron energy spectra both in the looptop and loopleg regions. We also find that the velocity space distribution varies both along the loop and with time. With non-/weak turbulent scattering, it presents different patterns of stripes and is highly anisotropic both in the looptop and loopleg, and gradually evolves into a double-sided loss-cone as the simulation continues. In the case of enhanced turbulent scattering, the distribution becomes nearly isotropic because a large number of electrons can be scattered into the loss-cone.

In flare regions with strong magnetic fields, the plasma frequency can be smaller compared to the electron gyro-frequency, i.e., ω pe / Ω ce < 1 . Different anisotropic features in the velocity space can drive the electron cyclotron maser instability (ECMI) in different manners. In the looptop region, most electrons are distributed in the perpendicular direction, with horizontal branches stretching outward. Positive gradients along the loss-cone boundaries mainly generate fundamental X-mode emissions via ECMI, propagating along the parallel and oblique directions (see, e.g., Yoon and Ziebell, 1995; Ning et al., 2021b). According to the plasma kinetic theory, the linear growth rates of ECMI can be approximated with the integral of the velocity distribution function gradient ( ∂ f / ∂ v ⊥ ) along the resonance curve in the phase space (Wu and Lee, 1979; Wu, 1985). We note that in the non-scattering case, the distribution presents branch features with sharp gradients where the resonance curve could pass through. This could drive ECMI with higher efficiency, compared to the classical loss-cone distribution. In the loopleg region, the distribution resembles the horseshoe distribution in the source of planetary AKR. Recently, the horseshoe-driven ECMI has been applied to explain the solar spikes (e.g., Melrose and Wheatland, 2016; Ning et al., 2021a). Multi-stripe distribution has been demonstrated in earlier studies (e.g., White et al., 1983). Yousefzadeh et al. (2021) carried out kinetic simulations and found that such electrons mainly generate second harmonic X-mode emissions (X2), which could solve the escaping difficulty of fundamental emission in solar corona. For the strong scattering case, the distributions in both regions are nearly isotropic, making it hard to drive the ECMI.

We used the reduced transport equation that includes magnetic mirroring and turbulent pitch-angle scattering in this study. Other effects such as Coulomb collisions and cross-filed diffusion have been neglected and will be discussed in future work. We considered different regimes of turbulent pitch-angle diffusion as defined in Bespalov et al. (1987). However, the level of magnetic turbulence in realistic flares remains unclear. Recently, some studies (Kontar et al., 2017; Stores et al., 2021) investigated the spatial and temporal distributions of turbulence in one solar flare from the observations of nonthermal broadening of spectral lines. They calculated the turbulent kinetic energy density from the nonthermal broadening velocity ( v nth ) , which approximates the energy density associated with the magnetic field fluctuations. Then, one can estimate the level of turbulent magnetic fluctuation δ B / B ∼ v nth / v A (Kontar et al., 2017). Taking the Alfven speed v A ∼ 2000 km s − 1 , while v nth ranging between ∼ 10–100 km s − 1 , we can get δ B / B is about 0.05%–5%. This indicates that the turbulence is relatively weak on average, at least for this flare event, and may provide the required condition for anisotropic distribution.

A similar electron trap and precipitation process occurs in the Earth’s magnetosphere, where the resonant interaction between energetic electrons and plasma waves such as chorus waves has been applied to explain the characteristics of aurora (e.g., Thorne et al., 2010; Zhang et al., 2022). Our simulation results may provide helpful insights to the dynamics of energetic particles in the radiation belts of magnetosphere.