Butterfly Distributions of Energetic Electrons Driven by Ducted and Nonducted Chorus Waves

Bursts of electron butterfly distributions at 10s keV correlated with chorus waves are frequently observed in the Earth's magnetosphere. Strictly ducted (parallel) upper‐band chorus waves are proposed to cause them by nonlinear cyclotron trapping. However, chorus waves in these events are probably nonducted or not strictly ducted. In this study, test‐particle simulations are conducted to investigate electron scattering driven by ducted (quasi‐parallel) and nonducted upper‐band chorus waves. Simulation results show butterfly distributions of 10s keV electrons can be created by both ducted and nonducted upper‐band chorus waves in seconds. Ducted upper‐band chorus waves cause these butterfly distributions mainly by accelerating electrons due to cyclotron phase trapping. However, nonducted waves tend to decelerate electrons to form these butterfly distributions via cyclotron phase bunching. Our study provides new insights into the formation mechanisms of electron butterfly distributions and demonstrates the importance of nonlinear interactions in the Earth's magnetosphere.

Drift shell splitting and magnetopause shadowing are theoretically energy independent, which can cause the butterfly PADs of tens of keV to several MeV electrons at larger L-shells (e.g., Fritz et al., 2003;Sibeck et al., 1987).Wave-particle interactions that are energy dependent are also an important contributor to the formation of electron butterfly PADs especially at lower L-shells (e.g., Horne et al., 2007;J. Li et al., 2016;Xiao et al., 2015).Many previous studies have shown that the butterfly PADs of hundreds of keV to several MeV electrons are usually induced by magnetosonic (MS) waves, plasmaspheric hiss, and chorus waves (e.g., Albert et al., 2016;Hua et al., 2019;J. Li et al., 2014J. Li et al., , 2016;;Ma et al., 2015;Ni et al., 2017Ni et al., , 2018Ni et al., , 2020;;Xiao et al., 2015).In recent years, the bursts (within ∼30 s) of electron butterfly PADs at tens of keV, with high correlations with chorus waves, are frequently observed in the Earth's magnetosphere (Fennell et al., 2014;Kurita et al., 2018;Peng et al., 2022).Observational statistics reveal that more than 80% of these events are upper-band chorus dominated events and these events tend to occur over ∼21-05 MLT at L = ∼4.5-6.5 (Peng et al., 2022).By assuming these chorus waves are strictly ducted (parallel), previous test-particle simulations suggest nonlinear cyclotron trapping is mainly responsible for the bursts of these butterfly PADs (Gan, Li, Ma, Artemyev, & Albert, 2020;Saito et al., 2021).However, chorus waves in these events are probably nonducted or not strictly ducted, since many of these events are observed without a density structure (see Figure S1).Chorus waves can be ducted (quasi-parallel) in a density structure and propagate nearly along the magnetic field lines (e.g., R. Chen et al., 2021;Ke et al., 2021;Liu et al., 2021;Streltsov et al., 2006).However, nonducted chorus waves tend to become obliquepropagating gradually (e.g., Breuillard et al., 2012;Gao, Lu, et al., 2022;Lu et al., 2019) after leaving their equatorial sources (e.g., LeDocq et al., 1998;W. Li et al., 2009;Santolík et al., 2005).The relations between nonducted chorus waves and electron butterfly PADs at tens of keV is less well-known.In this study, test-particle simulations combined with electron magnetohydrodynamics (EMHD) simulations have been used to investigate electron scattering driven by ducted and nonducted upper-band chorus waves.Simulation results show that ducted and nonducted upper-band chorus waves lead to the butterfly PADs of tens of keV electrons through different nonlinear processes.

Simulation Model
To model chorus waves in the inner magnetosphere, the geomagnetic field is set as a dipole magnetic field , where B 0eq ≈ 3 × 10 4 nT/L 3 and λ is the magnetic latitude.The electron density (in cm 3 ) in the plasma trough is set based on an empirical model (Carpenter & Anderson, 1992;Denton et al., 2002Denton et al., , 2006)), (1) here a is set to 1 and n e0 = 5,800 + 300 MLT (0 ≤ MLT < 6).In normal (nonducted) cases, n d = 0.In the presence of a density duct, n d is given as where L 0 and D L are respectively the central location and the radial width of this duct.We perform three simulation cases: a ducted case with L 0 = 5.5, D L = 0.035, and δn = 0.1n eq (L 0 ), and two nonducted cases with different wave source scales.Our simulations are carried out around L = 5.5 at MLT = 0, which are typical conditions for bursts of electron butterfly PADs based on observational statistics (Peng et al., 2022).
A two-dimensional (2-D) EMHD model (Ke, Gao, et al., 2022) is used to simulate the propagation of upper-band chorus waves.2-D EMHD simulations are widely used to study the propagation of chorus waves (e.g., Hanzelka & Santolík, 2022;Hosseini et al., 2021;Katoh, 2014).Chorus waves are excited by the energetic electrons injected from the magnetotail and drifting eastward around the Earth, which usually consist of quasi-periodical discrete elements with frequency chirping and a majority of them exhibit rising-tone elements (Burtis & Helliwell, 1969;W. Li et al., 2012;Tsurutani & Smith, 1974).The repetition period of chorus elements is found to be correlated with the drift velocity of energetic electrons (Gao, Chen, et al., 2022;Lu et al., 2021), which has a typical value of ∼0.5 s, and each element lasts from ∼0.1 to 1 s (Shue et al., 2015(Shue et al., , 2019;;Teng et al., 2017).The transverse source scale of chorus elements is estimated in the range ∼100-800 km (∼0.016-0.126R E , R E is the earth radius) based on multiple-satellite measurements (e.g., Agapitov et al., 2017Agapitov et al., , 2018;;Santolík & Gurnett, 2003;Shen et al., 2019).In our simulations, parallel upper-band chorus waves with a repetition period

Geophysical Research Letters
10.1029/2024GL108307 T RP = 0.5 s and a typical element duration T D = 0.3 s are launched from an equatorial source region.For each chorus element, its frequency f rises evenly from 0.5f ce to 0.65f ce and its amplitude remains constant B w0 = 0.0011B 0eq,0 = 0.2 nT in the middle phase, but rises (drops) in the initial (end) phase T ini (T end ) = 0.1T D (B 0eq,0 is B 0eq (L 0 ) and f ce is the equatorial electron gyrofrequency at L 0 ).The wave spectra are shown in Figure 1a.
The chorus source scale is set as ΔL = 0.042 (L = 5.479-5.521) in Ducted case and Nonducted Case 1, and 3ΔL = 0.126 (L = 5.479-5.605) in Nonducted Case 2.Moreover, the wave amplitudes decrease from B w0 to 0 in the inner and outer edges (with a width δ L = 0.0042) of each source region.The simulation domain within λ ≈ 15°-15°, has L = 5.4-5.55 in Ducted case and Nonducted Case 1, and L = 5.4-5.63 in Nonducted Case 2. The L-shell width in Case 2 is 1.533 times that in Case 1. Absorbing boundary conditions are applied for waves.
Assuming the wave fields are confined to |λ| < 15°is reasonable for upper-band chorus waves according to satellite statistics (e.g., Meredith et al., 2009;Teng et al., 2019).The simulation grid numbers in parallel and perpendicular directions are N ∥ = 16,000 and N ⊥ = 1,000 (or 1,533) for Ducted case and Nonducted Case 1 (or Nonducted Case 2).
Test-particle simulations are used in combination with the 2-D EMHD simulations.Initially, test electrons are located at the equator at L 0 = 5.5.They are uniformly distributed in the energy E from 2 to 90 keV with ΔE = 2 keV, the equatorial pitch angle α eq from 5°to 89°with Δα eq = 2°(the loss cone ∼3°), and the gyroangle φ from 0°to 354°with Δφ = 6°.A particle weight is given to each test electron according to the initial electron flux distribution assumed as where j 0 is the differential flux at E 0 = 10 keV and α eq = 90°.The index p is set as 1 and 2 at energies less and greater than E 0 , respectively.The initial electron flux distribution is shown in Figure 1b.The similar power law distributions of energetic electrons are widely used in previous works (e.g., Bortnik et al., 2011;L. Chen et al., 2012;Saito & Miyoshi, 2022).The simulation time step is 1.5771 × 10 6 s and the total simulation time is 3 s.

Simulation Results
Figure 1 shows an overview of ducted and nonducted propagations of the chorus waves in our EMHD simulations.Figures 1c-1e show the spatial profiles of magnetic field amplitudes of upper-band chorus waves at t = 0.2 s (marked by the dashed line in Figure 1a) in three simulation cases.Figure 1c illustrates that upper-band chorus waves propagate nearly along the magnetic field lines in a density duct of 10% density reduction.Obviously, these upper-band chorus waves are ducted by the depleted duct, consistent with previous studies (e.g., Liu et al., 2021;Smith et al., 1960).Without a density duct, upper-band chorus waves gradually deviate from the magnetic field lines across the wave source region (simply called source field lines) during their propagation.Most of these nonducted waves completely deviate from source field lines at |λ| < 5°in Nonducted Case 1 and at |λ| < 10°in Nonducted Case 2 (Figures 1d and 1e).The averaged amplitude B w and wave normal angle θ of these chorus waves at f/f ce = 0.55-0.6along L 0 = 5.5 are estimated (in Ducted case and Nonducted Case 2) and shown in Figures 1f and 1g, respectively.In Ducted case, B w and θ of these chorus waves fluctuate around 0.0011B 0eq,0 and 8°at |λ| ≤ 11.5°, respectively.Both B w and θ are not shown at |λ| > 11.5°, since these waves decay sharply due to absorbing boundary conditions.In Nonducted Case 2, B w / B 0eq,0 of these chorus waves decreases from 0.0011 to 0.0006 within |λ| ≤ 7.5°, and then drops quickly to 0.0002 within |λ| = 7.5°-9°, since more nonducted waves deviate from source field lines at higher latitudes.Besides, θ of these nonducted waves increases rapidly with the latitude, reaching up to 36°at |λ| = 9°.
Figure 2 shows the simulation results of the flux distributions of energetic electrons projected on the equatorial plane.Figures 2a-2c present the electron flux distributions at t = 3 s in three simulation cases.Compared to the initial flux distributions (Figure 1b), the flux distributions at t = 3 s are greatly modified by the chorus waves, which are clearly different in ducted and nonducted cases.In Ducted case, chorus waves cause significant flux decreases and increases at α eq ∼ 40°-60°and α eq ∼ 60°-80°, respectively.In Nonducted cases, chorus waves cause both significant flux increases and decreases at α eq ∼ 60°-80°.Besides, the significant flux increases in Nonducted cases occur at higher energies than those in Ducted case.The greatest relative flux increase j t=3s / j t=0s  time.These PADs at E = 18 keV become butterfly PADs only in Ducted case.We identify a butterfly PAD based on two criteria.The first is similar to that of Ni et al. (2016): j(90°) < β × j avg (α: 90°), where j(90°) is the electron flux at α eq = 90°, j avg (α: 90°) is the averaged electron flux over the pitch angle range between α and 90°, α is selected from 45°to 85°to determine the maximum value of j avg (α: 90°), and the index β is set as 0.95.This criterion is to identify an electron PAD with a minimum flux at 90°.The second is j max ≥ η × j(90°), where j max is the maximum flux at α eq = 0-90°and the threshold η is set as 1.5.This criterion is to exclude those PADs with slight flux variations, like the PADs at E = 18 keV in Nonducted Case 1 at t ≥ 2 s (Figure 2e).Based on the criteria, the energy ranges (widths) of the butterfly PADs are estimated as 16-44 keV (28 keV), 36-58 keV (22 keV), and 36-54 keV (18 keV) in Ducted case, Nonducted Case 1, and Nonducted Case 2, respectively.
The ratios of the electron fluxes at t = 3 s and t = 0 s in three simulation cases are shown in Figures 3a-3c, which indicate where the fluxes increase or decrease.In Ducted case, the large flux increases appear at E ∼ 16-46 keV and α eq ∼ 60°-80°, while the evident flux decreases appear at lower energies and pitch angles (Figure 3a).In Nonducted Case 1 and Case 2, the large flux increases occur at E ∼ 38-58 keV and α eq ∼ 60°-75°, while the main flux decreases occur at higher energies and pitch angles (Figures 3b and 3c).A box with ΔE = 8 keV and Δα eq = 8°is located in a region with the largest averaged flux ratio in each panel of Figure 3.The central location (E, α eq ) of this box is (39 keV, 74°) in Ducted case, (47 keV, 71°) in Nonducted Case 1, and (46 keV, 66°) in Nonducted Case 2. Figures 3d-3f show the initial flux distributions of the electrons scattered into the box (at t = 3 s) in three cases.In Ducted case, 91% of the scattered electrons are accelerated, and their α eq also increase.
Most of them are initially distributed at E ∼ 15-35 keV and α eq ∼ 15°-70°.However, 77% (or 84%) of the scattered electrons are decelerated in Nonducted Case 1 (or Case 2), and their α eq also decrease.Most of them are initially distributed at energies up to ∼63 keV and pitch angles up to ∼82°.Obviously, these electrons scattered into the box are accelerated or decelerated by upper-band chorus waves through the cyclotron resonance rather than the Landau resonance, which causes the opposite variations of E and α eq .In a word, ducted upper-band chorus waves cause these butterfly PADs of 10s keV electrons mainly by accelerating electrons, while nonducted waves do the opposite.
The typical trajectories of an accelerated electron in Ducted case and a decelerated electron in Nonducted Case 2 are presented in Figure 4, which shows the energy E, pitch angle α eq , and parallel velocity v ∥ as functions of the

Geophysical Research Letters
10.1029/2024GL108307 magnetic latitude.This electron in Ducted case, with an initial energy of E ∼ 30 keV, is first accelerated to ∼50 keV from λ = 8°to λ = 2°, and then undergoes several rapid decelerations at |λ| ≤ 5°, finally reaching 40 keV (Figure 4a).The variation trend of α eq is very similar to that of E for this electron (Figure 4b).The parallel velocity v ∥ decreases with fluctuating around the cyclotron resonant velocities v R in the acceleration phase, and increases sharply with intersecting v R in the rapid deceleration phases (Figure 4c).Obviously, this electron first undergoes nonlinear phase trapping and then potentially experiences phase bunching.This electron in Nonducted Case 2, with an initial energy of E = 52 keV, only undergoes several rapid decelerations, like the electron decelerations in Ducted case.It probably experiences several times of phase bunching at |λ| ≤ 7°, and finally drops to 42 keV (Figures 4d-4f).
If an electron is phased trapped, the phase ζ between its perpendicular velocity v ⊥ and the wave perpendicular magnetic field B w⊥ will change periodically and satisfy 0 < ζ < 2π.In simulation results, phase trapping of an electron is identified by the criterion: the phase ζ satisfies 0 < ζ < 2π for more than five periods.A group of electrons with the same initial E and α eq but uniformly distributed gyrophases are scattered roughly symmetrically by quasilinear scattering, but most of them are scattered to lower E and α eq by phase bunching (e.g., Bortnik et al., 2008).Generally, phase bunching leads to significantly larger energy and pitch angle changes than those caused by quasilinear scattering.We assume such a group of electrons are in a quasilinear scattering period until , where ΔE is the energy change of each electron within a half bounce period T hb (Gan, Li, Ma, Albert, et al., 2020).An electron is assumed to be phase bunched if ΔE < |ΔE| max and Δα eq < |Δα eq | max within T hb , where |ΔE| max and |Δα eq | max are the maximum energy and pitch angle changes within T hb of all the electrons in their quasilinear scattering periods.In three simulation cases, |ΔE| max and |Δα eq | max are found to be about 1 keV and 1°.For the electrons scattered into the box in Ducted case, 78% of them have experienced phase trapping, and 81% of them have experienced phase bunching.Finally, most of them get accelerated as a result of competition between phase trapping and phase bunching.Besides, 64% of the trapped electrons begin their trapping at |λ| > 9°.For the electrons scattered into the box in Nonducted Case 1 (or Case 2), 72% (or 86%) of them have experienced phase bunching, but only 22% (or 25%) of them have experienced phase trapping.Finally, most of them get decelerated due to multiple phase bunchings, which mainly occur at |λ| < 5°.

Discussion and Conclusions
In this study, test-particle simulations in combination with 2-D EMHD simulations have been performed to study the scattering of 10s keV electrons driven by ducted and nonducted upper-band chorus waves.Our simulation results suggest ducted waves cause the butterfly PADs of 10s keV electrons mainly by accelerating electrons via cyclotron phase trapping, while nonducted waves cause these butterfly PADs mainly by decelerating electrons via cyclotron phase bunching.We explain the difference based on nonlinear resonant conditions (Bell, 1984(Bell, , 1986;;Tao & Bortnik, 2010).Cyclotron phase trapping is possible under the condition: which are described in Equations 2-4 of Bell (1986).When the chorus waves resonate with an electron of 10s keV at α eq < ∼60°at |λ| < 10°, there are , where μ 0 is the permeability of vacuum.The blue and red dotted lines mark the cyclotron resonant velocities v R at f = 0.5f ce and f = 0.7f ce .The black and blue asterisks mark the start and end of the trajectories.

10.1029/2024GL108307
where q e and m e are the electron charge and mass, k ∥ is the parallel wave number, B w⊥1 and B w⊥2 are two components of the perpendicular wave magnetic field.The inhomogeneity factors h(z, t) at |λ| < 10°for ducted and nonducted waves are slightly different, which are mainly related to the inhomogeneity of the background magnetic field ∂Ω e ∂z and the wave frequency chirping rate ∂ω ∂t .Obviously, the wave magnetic field amplitudes B w dominate the difference in ω 2 t / h(z,t) of ducted and nonducted waves.Phase bunching is easy to occur when ω 2 t / h(z,t) is much greater than 1.Thus, it easily occurs for both ducted and nonducted waves at lower latitudes (| λ| ≤ 5°) because of the smaller h(z,t) due to the smaller ∂Ω e ∂z ( ∂Ω e ∂z = 0 at λ = 0°) .The amplitudes B w of ducted waves decrease slightly along the magnetic field lines due to the ducted propagation, and easily satisfy ω 2 t / h(z, t) > 1 even at higher latitudes (|λ| > 5°).Thus, ducted waves trap and accelerate electrons to produce the butterfly PADs.Although phase bunching is also involved in the electron dynamics, the acceleration effect due to phase trapping is dominant.The amplitudes B w of nonducted waves decrease rapidly along the magnetic field lines due to the nonducted propagation, and hardly satisfy ω 2 t / h(z, t) > 1 at higher latitudes.Thus, nonducted waves mainly decelerate electrons to causes the butterfly PADs, since phase bunching is dominant in the electron dynamics.Besides, we have also performed two simulation cases with the smaller amplitude B w0 = 0.05 nT (the other parameters are same to those of Ducted case and Nonducted Case 1) and another two simulation cases with L-shell around L 0 = 6.5, and the simulation results give the same conclusions.
Nonducted chorus waves potentially accelerate electrons by Landau trapping, which contributes less to the formation of electron butterfly PADs in our simulations.In Ke, Lu, et al. (2022), a lower-band chorus with a constant frequency f/f ce = 0.4 and the amplitude B w0 = 0.005B 0eq,0 can form the butterfly PADs of tens of keV electrons by Landau trapping.There are two main reasons for the difference.One reason is that the lower-band chorus almost propagates along the magnetic field line at |λ| < 15°.Another reason is that the amplitude of the lower-band chorus is much larger.Therefore, the lower-band chorus is easier to satisfy the condition of Landau trapping for tens of keV electrons at large pitch angles.
In the realistic magnetosphere, the magnetic amplitude of nonducted upper-band chorus waves may decrease faster due to Landau damping (Bortnik et al., 2007), which is not included in our EMHD simulations.Thus, nonducted upper-band chorus waves are more likely to form the electron butterfly PADs through phase bunching.However, the observational evidence is still lacking.Kurita et al. (2018) provided an observation event detected by Arase, which indicates the electron butterfly PADs result from acceleration of lower energy electrons.Unfortunately, the plasma density data are absent in the same period.It is not easy to determine whether the electron butterfly PADs observed by Van Allen Probes in previous studies originate from lower or higher energy electrons due to low resolution (Fennell et al., 2014;Peng et al., 2022).Thus, the observational evidence is also required as a future work.The primary conclusions are summarized as follows.
1. Test-particle simulations in combination with 2-D EMHD simulations demonstrate both ducted and nonducted upper-band chorus waves can form significant butterfly distributions of tens of keV electrons within seconds.2. Ducted upper-band chorus waves tend to cause the butterfly distributions of tens of keV electrons by accelerating electrons via cyclotron phase trapping.3. Nonducted upper-band chorus waves tend to decelerate electrons to cause the butterfly distributions of tens of keV electrons via cyclotron phase bunching.

Figure 1 .
Figure 1.(a) The frequency-time spectrogram of upper-band chorus waves launched from the magnetic equator and (b) the initial electron flux distribution at the magnetic equator in our simulations.(c)-(e) The spatial profiles of magnetic field amplitudes of chorus waves at t = 0.2 s in three simulation cases.(f), (g) The averaged amplitude B w and wave normal angle θ of these chorus waves at f/f ce = 0.55-0.6along L 0 = 5.5 in Ducted case and Nonducted Case 2.

Figure 2 .
Figure 2. (a-c) The electron flux distributions as a function of the equatorial pitch angle and the energy at t = 3 s and (d-f) the pitch angle distributions of electrons at two energy channels at t = 0, 1, 2, and 3 s in three simulation cases.

Figure 3 .
Figure 3. (a-c) The ratios of the electron fluxes at t = 3 s and t = 0 s in three simulation cases.A box with ΔE = 8 keV and Δα eq = 8°is located in a region with the largest averaged flux ratio in each panel.(d-f) The initial flux distributions of the electrons scattered into the box (at t = 3 s) in three simulation cases.

Figure 4 .
Figure 4. (a-c) The energy E, equatorial pitch angle α eq , and parallel velocity v ∥ of a typical accelerated electron in Ducted case as functions of the magnetic latitude λ. (d-f) The E, α eq , and v ∥ of a typical decelerated electron in Nonducted Case 2 as functions of λ.The unit ofv ∥ is V Ae = B 0eq,0 / ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ μ 0 m e n eq (L 0 ) √, where μ 0 is the permeability of vacuum.The blue and red dotted lines mark the cyclotron resonant velocities v R at f = 0.5f ce and f = 0.7f ce .The black and blue asterisks mark the start and end of the trajectories.