Simulations of Electron Energization and Injection by BBFs Using High‐Resolution LFM MHD Fields

We study electron injection and energization by bursty bulk flows (BBFs), by tracing electron trajectories using magnetohydrodynamic (MHD) field output from the Lyon‐Fedder‐Mobarry (LFM) code. The LFM MHD simulations were performed using idealized solar wind conditions to produce BBFs. We show that BBFs can inject energetic electrons of few to 100 keV from the magnetotatail beyond −24 RE to inward of geosynchronous, while accelerating them in the process. We also show the dependence of energization and injection on the initial relative position of the electrons to the magnetic field structure of the BBF, the initial pitch angle, and the initial energy. In addition, we have shown that the process can be nonadiabatic with violation of the first adiabatic invariant (μ). Further, we discuss the mechanism of energization and injection in order to give generalized insight into the process.


Introduction
The plasma flow in the magnetosphere is highly structured in space and time, and sometimes it can be bursty. Angelopoulos et al. (1992) used the term bursty bulk flows (BBFs) to characterized the flows in the plasma sheet. They found that BBFs are enhanced (>100 km/s) flow for a 10-min time scale with a peak flow of >400 km/s for 1-min time scale. BBFs have strong induced electric field, and they are associated with dipolarization (Nakamura et al., 2002;Ohtani et al., 2004;Runov et al., 2009). The enhanced B z is usually preceded with a decrease in B z and weakened magnetic field strength (Ohtani et al., 2004;Yao et al., 2015;Zhou et al., 2014). BBFs are believed to be created by localized magnetic reconnection (Chen & Wolf, 1993). They are also believed to represent flux tubes of low entropy, which convect faster than the bulk of the plasma (Pontius & Wolf, 1990). MHD simulations using the Lyon-Fedder-Mobarry (LFM) code under idealized solar wind conditions have reproduced the observed qualitative features of the flow and the magnetic field pattern of the BBFs (Wiltberger et al., 2015).
BBFs are associated with important magnetospheric phenomena including injection of ring current ions Ukhorskiy et al. (2018) and electrons . Injection of tens to hundreds of kiloelectron volt electrons is a signature of BBFs (Gabrielse et al., 2014;Sergeev et al., 2012). The role of the fields of BBFs in accelerating and injecting plasma sheet electrons is a subject of many papers (Ashour-Abdalla et al., 2011;Birn et al., 1998Birn et al., , 2013Birn et al., , 2014Birn et al., , 1998Birn et al., , 2004Gabrielse et al., 2016Gabrielse et al., , 2017Li et al., 2003;Sarris et al., 2002). Li et al. (2003) have done test particle simulations by adding a propagating electromagnetic magnetic field pulse with a constant speed of 100 km/s to the dipole field of the Earth. Sarris et al. (2002) have performed a similar study with a radially dependent velocity of the pulse. Gabrielse et al. (2017Gabrielse et al. ( , 2016 have traced electrons using modeled dipolarizing flux bundles that self-consistently define the velocity and the electromagnetic field. They have shown that dipolarizing flux bundles can explain short (10-min) time scale electron energization and injection. Gabrielse et al. (2016Gabrielse et al. ( , 2017, Li et al. (2003), and Sarris et al. (2002) traced 90 • pitch angle electrons assuming that electrons conserve their first adiabatic invariant ( ). Birn et al. (2004) have studied the role of diplorization in acceleration and injection of electrons by tracing electrons with fields of a three-dimensional Magnetohydrodynamic (MHD) simulation of neutral line formation. They traced all pitch angle electrons and did not assume the conservation of . They concluded that dipolarization accelerates electrons via betatron and Fermi mecanisms for high and low pitch angle electrons respectively. Ashour-Abdalla et al. (2011) traced electrons with MHD fields including electron kinetic effects. They concluded that electrons are energized by an earthward propagating dipolarization front consistent with the betatron mechanism.

LFM MHD Simulation
The LFM MHD simulation solves the ideal MHD equations on an irregular grid optimized for magnetosphere studies and subject to specified solar wind input (Lyon et al., 2004). Simulations were performed to reproduce the observed qualitative features of BBFs using idealized solar wind conditions (Wiltberger et al., 2015). The simulation was done with a high-resolution grid of (212, 192, and256) in radial, polar, and azimuthal directions, respectively. The LFM-Magnetosphere-Ionosphere Coupler/Solver (MIX) version of the code, which treats the inner boundary at 2R E as a spherical surface with known conductivity, is used (Merkin & Lyon, 2010). The dipole field of the Earth was aligned along the solar magnetic z axis of the LFM coordinate. The LFM x axis points in the direction of the Sun, and the y axis completes the right-hand coordinate system. There is no corotational electric field in this simulation. The solar wind density was set at 5 cm −3 , the x component of velocity at −400 km/s. The y and z components of the velocity were set to 0. The density and the velocity were kept fixed throughout the simulations. The x and y components of the interplanetary magnetic field (IMF) were also set to 0. The z component of the IMF is −5 nT for the first 2 hr and flipped to 5 nT for the next 2 hr. Then, after the first 4 hr, the IMF was flipped back to −10 nT for the remaining time of the simulations. The LFM fields were dumped at every 1 s of the simulations. This time was chosen so that it is smaller than the ratio of the grid width to the propagation speed of the structure. Using a propagation speed of 400 km/s gives a time of ≈3 s. Using a time less than this value will allow the interpolated field to preserve the sharp gradients of the field structures.
The results of the MHD simulation show earthward propagating BBFs emerging from the magnetotail after about 5:00 hr of simulation time (ST), which proliferate later. For more detailed characterization of these BBFs, see the study by Wiltberger et al. (2015). Kress et al. (2007) have developed a code rbelt3d, to follow energetic particle trajectories in prescribed fields. The code integrates either the Lorentz equations or the guiding center equations based on a parameter,

Particle Tracing Code
, characterizing spatial gradients using either a prescribed field or output fields from MHD simulations. The guiding center equations (equations 6.22 to 6.26) of Cary and Brizard (2009) are used. The guiding center equations of Cary and Brizard (2009) are integrated numerically by the fourth-order Runge Kutta algorithm.
We added some modification to the code, highlighted as follows. In Kress et al. (2007), the Lorentz equation was integrated using various orders of Runge-Kutta methods. Here we integrate the Lorentz equation using a volume preserving method that can be applied to systems with time-dependent electromagnetic fields with good accuracy and conservative properties over long times of simulation (He et al., 2016). The method has been shown to have a long-term fidelity and higher accuracy compared to the traditional Runge-Kutta methods and Boris algorithms.  Littlejohn (1983) algorithm, which approximates the transformation up to an order of L , and the red line is zeroth-order approximations. Left is over several switching, and the right is zoomed in figure around a switch.
The switching between the Lorentz and the guiding center quantities in Kress et al. (2007) uses zeroth-order approximation, where the guiding center quantities are approximated by the particle variables. Here we use the relativistic version of the Littlejohn (1983) algorithm; see Boghosian (2003). The transformation between guiding center and particle variables is approximated up to order of L . The error of the algorithm is order of ( L ) 2 . The method allows us to track of an electron with a greater accuracy than the zeroth-order switch used by Kress et al. (2007). A comparison of a full Lorentz calculation with a switch based on the Littlejohn (1983) algorithm and switch of zeroth order is shown Figure 1. The result shows that the zeroth-order switch introduces error due to the gyrooscillation every time the switching is applied, whereas the switch based on the Littlejohn (1983) algorithm agrees quite well with the full Lorentz equation for several switches.
The main difference between the switch using the Littlejohn (1983) algorithm and the full Lorentz calculations is seen in the part where the Lorentz equations are used in both cases. This is because we do not follow the gyrophase. When we switch from the guiding center to the Lorentz trajectories, a phase difference will be introduced. This phase difference is manifested as different oscillations between the Lorentz trajectories of the full Lorentz calculations and the one using the switch. Even though the two oscillate with different phases, their phase-averaged values are about the same. In fact, the difference in their guiding center values is only of the order of ( L ) 2 . When switching back to guiding center, the difference becomes smaller than what appears in the Lorentz trajectories. The gyrophase is not necessary for higher-order approximations of and also the momentum. Thus, after many switches back and forth, the guiding center parallel velocity, pitch angle, and can be approximated well without the gyrophase information.
The parameter used was chosen in order to take into account the magnetic rotational shear, curvature, and scale length effect in one parameter in an ad hoc manner. However, there is no theoretical grounds for the choice. A more rigorous analysis of a parameter is given by Pfefferle et al. (2015). However, the evaluation of their parameter is too computationally intensive for the general electromagnetic field. Here we use in terms of the magnetic curvature, gradient, and shear scale lengths. Even though this parameter is not exactly the same as that of Pfefferle et al. (2015), it bounds all the possible scale lengths. The code switches to integrating guiding center equations when < 0.005 and to the Lorentz equation if ≥ 0.01.
The LFM MHD output fields have a localized nonzero E || due to numerical resistivity. In addition, the grid-to-particle interpolation, which is quadrilinear in the three spatial and time coordinates, introduces E || even if the fields at the grid position do not have E || . This E || , even though small compared to E ⟂ , has nonphysical origin. We set E || = 0 at the particle position to eliminate numerical and nonphysical E || . The numerical resistivity introduces electric field in the perpendicular direction comparable to the neglected E || , about a factor of 2 larger, but it is much smaller than the inductive electric field. The red circle is the geosynchronous orbit. The snapshots are taken to show some of the main features described in the text; the full movie can be seen on Movie S1 of the supporting information. ST = simulation time.
When the particle is following a guiding center trajectory, we can separate the energization terms based on the following equation (Northrop, 1963) useful in describing physical processes: Since we do not have E || in the MHD fields, the term containing ⃗ v || can be dropped. The equation can be used if the guiding center position, and hence the drift velocity, (⃗ v D ), is well defined. However, the simulation under BBF fields shows that there is a spatially confined region where the electron motion can be chaotic due to = √ R c being close to 1. This happens near the magnetic equator and is localized in the x-y plane adjacent to the enhanced z component magnetic field perturbationB z . In chaotic motion the conservation of is violated (Buchner & Zelenyi, 1989), and the guiding center position and the drift velocity are not well defined. We integrate the terms of equation (2) separately when the particle is in its guiding center trajectory. However, in order to integrate the two terms of equation (2), we need the guiding center position and momentum, which are found from the guiding center equation of Cary and Brizard (2009). The change in energy calculated from equation (2) is not exactly the same as that of the guiding center equation of Cary and Brizard (2009). Even though the guiding center equation of Cary and Brizard (2009) has several advantages over equation (2), it does not separate the different energy terms like that of equation (2), which is useful for describing the acceleration mechanism. To complete the energy terms, we integrate q ⃗ E · ⃗ v with time during the Lorentz trajectory, where ⃗ v is particle velocity.

Simulations and Results
We traced energetic electrons within BBFs for different initial energy and position relative to one selected BBF structure chosen as a typical example from the LFM MHD simulation results. One BBF is chosen because the qualitative features of the interaction of energetic electrons with BBFs are similar. In addition, the BBFs were produced by MHD simulations using idealized solar wind conditions so that they do not represent any particular event. The simulations were done for initial energies of 3, 10, and 100 keV chosen to cover the wide range of superthermal energetic electrons seen in the magnetotail. The distribution in pitch angle was isotropic at the magnetic equator. For each energy, the electrons were distributed randomly on a chosen x-y location of 0.4R E × 0.4R E . The spatial width of 0.4R E is chosen based on the grid resolution of the LFM in the midtail region, which is ≈ 0.2 R E . Three spatial positions were chosen for each initial energy to demonstrate how the interaction of the BBF with the electrons depends on the spatial position relative to the BBF. The regions are the center, the leading edge of the enhancedB z , and a magnetic field minimum position on the duskside of the BBF. For each energy and position, 3,600 electrons were traced. These electrons are different in initial pitch angle and location within spatial width of 0.4R E × 0.4R E . Since the spatial width is within a ± of the grid resolution of LFM, the number is mainly to cover the initial pitch angle distribution.  time cadence is shown. If an electron does not cross the equator within the time cadence, it is not shown on the snapshot. The snapshot shows how the bulk of the simulated electrons behave. Supporting movies for all cases of the simulations are presented in the supporting information. In general, the result shows that electrons are injected and accelerated with the BBF structure. The detail of the energization and injection depends on the initial position and initial energy of the electrons, which will be discussed in the following sections. The mechanism of the injection and the acceleration will also be presented in the following sections.

Mechanism of Injection
We use the drift approximation to explain the mechanism of the injection even though is not always conserved for injected electrons. The mechanism for the violation of conservation of is described later.
The assumption here is that the violation of the conservation of happens in a localized region close to the magnetic equator and a localized region on the x-y plane, while for most electrons the guiding center drift trajectories can explain the overall drift motions. The main drifts are ⃗ E × ⃗ B drift, which is the same as the perpendicular component of the plasma flow velocity, and the ∇B and the curvature drifts, which depend on the position relative to the BBF. Figures 2 -7 show time sequences of electrons injected from three radial locations relative to a BBF located in the dusk-midnight sector for three different initial energies: 10, 100, and 3 keV. The simulations for initial energy of 3 keV are shown only in the supporting movies included with the supporting information.
B dt and the displacement of the electron shows good agreement for this region; see Figure 8. Earthward, the ⃗ E × ⃗ B drift, which is the same as the perpendicular component of the flow weakened, so that the ΔB drift starts to dominate. The injection inward of r ≈ 10R E is a combination of ΔB drift and ⃗ E × ⃗ B drift directed toward the Earth. As the electrons move inward, the dipole magnetic field starts to dominate such that electrons drift along the ΔB drift of the dipole. This drift, which is directed eastward, causes electrons to move away from the BBF when they are on the east side of the BBF. We refer to this as detrapping. The position where the electron starts to detrap depends on the electron energy, with higher energies being detrapped farther away from the Earth as the ΔB drift increases with energy. Figure 3 shows the results of the injection for initial energy of 10-keV electrons, initially placed at the minimum B on the duskside of the BBF shown in Figure 2. For this position, the ΔB drift is tailward, which is mainly opposite to the flow ( ⃗ E × ⃗ B drift). This result shows the competition between the ΔB and ⃗ E × ⃗ B drifts, where electrons with larger ΔB drift move tailward and those with larger ⃗ E × ⃗ B move earthward with the BBF. This result is consistent to Gabrielse et al. (2016) of tailward drifting electron when they are on the duskside of the BBF. However, their result is only for 90 • pitch angle electrons. The main difference between tailward-and earthward-moving electrons is their initial equatorial pitch angle. Electrons with larger equatorial pitch angle stay near the equator where the ΔB drift is large; thus, their motion is dominated by the ΔB drift, and they move tailward relative to the BBF. Even though electrons with small pitch angle have large curvature drift velocity, the integrated drift is not as important as the curvature drift. Further discussion of the dependence on pitch angle and the pitch angle evolution is presented in the section 7.
Electrons that move tailward encircle the enhancedB z , but they do not travel with the BBF structure. Most of these electrons are not injected with the BBF, whereas those electrons that move earthward drift to the enhancedB z . Subsequently, they are injected with the BBF and energized like the case of starting in the enhancedB z region of Figure 2.
The third initial position investigated is the leading edge of the enhancedB z region. For these electrons, initially, the ΔB drift dominates over the ⃗ E × ⃗ B drift. As the snapshots show, electrons drift in the ΔB drift direction encircling the enhancedB z of the BBF. Later, some of the electrons move to a position where the ⃗ E × ⃗ B drift dominates so that they move with the plasma flow. The overall efficiency of the BBF in trapping and injecting is lower than for the initial location within the BBF shown in Figure 2.
We repeated the above simulations with initial energy of 100 and 3 keV. Figure 5 shows the result when the 100-keV electrons are initiated on the enhancedB z structure. The results show electrons initially following mainly the ⃗ E × ⃗ B drift like the 10-keV case. However, in this case, the ΔB drift is larger as a result of higher energy, so for some electrons, the ΔB drift dominates the ⃗ E × ⃗ B drift, which is independent of energy. The other difference is that at 100-keV initial energy electrons, which move inward with the BBF, start to detrap farther away from the Earth than in the 10-keV case, because of higher ΔB drift due to their larger energy. The result is that the injection does not go as deep into the inner magnetosphere for 100-keV initial energy electrons as in the 10-keV case. This finding is similar to that of Gabrielse et al. (2016), who showed that the electron motion is controlled by competition of ⃗ E × ⃗ B and ΔB drift with ⃗ E × ⃗ B drift dominating at low energy (<10 keV) and ΔB drift at high energy. Figure 6 shows the result when the electrons are initiated in the minimum B structure. The result is very similar to the 10-keV case except that the tailward drifting electron number increases. This is because the increased ΔB drift causes more electrons to have higher ΔB drift than ⃗ E× ⃗ B drift compared to the 10-keV case. Figure 7 shows the result when the electrons are started at the leading edge of the BBF. This result is similar to the 10-keV case for electrons that drift around the enhancedB z , but the ΔB drift is faster due to larger energy. Like the 10-keV case, the injection is not efficient.
Results for 3-keV initial energy electrons are very similar to the 10-keV electrons except that ΔB effects are weaker and more electrons move earthward due to ⃗ E × ⃗ B drift. The other difference is that when they are initiated at the enhancedB z , their energy by the time they reach geosynchronous is so small that they are not detrapped. Instead, the electrons evolve with the complex flow structure of the BBF, which moves duskward and then back toward the tail. When electrons move tailward they loose energy as they move to a weaker magnetic field strength region. The result can be seen in the Movies S7-S9 in the supporting information.
Electrons are injected as a result of ⃗ E × ⃗ B and ΔB (including curvature) drifts. The injection is efficient when ⃗ E × ⃗ B drift dominates over the ΔB drift, as it is the same as the perpendicular component of the flow and keeps the electrons within the BBF. Where the ΔB drift dominates, the direction depends on the position relative to the BBF and energy. In general, ΔB drift does not keep the electrons with the BBF. The ΔB drift becomes eastward closer to the Earth as the dipole field dominates closer to the Earth. This drift is away

Mechanism of Energization
To explain the energization, we integrate with respect to time the two terms of equation (2) when the electrons satisfy the guiding center approximation, and when electrons follow Lorentz trajectories, we integrate q ⃗ E · ⃗ v, where ⃗ v is the particle velocity. The first term of equation (2) is the energization due to the electric force acting along the gyromotion as a result of nonzero Δ × ⃗ E. It is also referred to as gyrobetatron acceleration (Fillius and McIlwain, 1967). The second term is the work done by the electric force in displacing the guiding center along with it. The ⃗ E · ⃗ v D term causes the electron to gain (lose) energy as a result of the work done by the electric field causing the electron to move to a stronger (weaker) magnetic field region. This can also be interpreted as the electric force acting along the drift path due to nonzero Δ × ⃗ E, and it is referred to as drift betatron acceleration (Fillius and McIlwain, 1967). The nature of the energization during the chaotic motion is a combination of the ⃗ E acting along the gyration and the displacement of the guiding center, but the two terms in this case are not separated.
To see the mechanism of energization, first, we choose the electrons that are injected and accelerated from the simulations. By injected and accelerated, we mean electrons that show a change of energy by at least a factor of 3 and that move inward of 7.5R E . This position is chosen because the flow brakes around this region. The criterion of selection is somewhat arbitrary, but the qualitative result does not change significantly with the criterion. The ensemble average of the change of energy and the energization terms are calculated. This was done for 3-, 10-, and 100-keV electrons initiated at the enhancedB z and minimum ofB z . We did not show the case when electrons are initiated at the leading edge of the BBF, as most of those electrons are not injected and it does not give further insight into how BBFs contribute to populate the radiation belt.
The results of the ensemble average energy terms and the average for all cases are shown in Figures 9 -12. First, we note that the change in energy calculated using equations of Cary and Brizard (2009) and the sum of the three energization terms in equation (2) are approximately equal even though the two methods are different. This indicates that splitting the energy terms using equation (2) does not have much error compared to the equations of Cary and Brizard (2009) derived from Hamiltonian theory. In general, the energization when the electrons are following Lorentz trajectories is the smallest among the three terms. This term is relatively larger when the electrons are initiated at the minimum B locations and when the initial energy is larger. This is because the parameter is small under these conditions that the electrons undergo more chaotic motion. The details of the chaotic motion will be summarized in the following section. The drift betatron term is always positive and largest. The result of the gyrobetatron term depends on the initial energy and initial position. This is because this term depends on the location of the electrons relative Figure 9. The ensemble average energy terms (left) and (right) of the injected electrons, for 10 keV initial energy started on the enhancedB z position of the BBF. E tot is the total energy, E b t is energization due to gyro-betatron, E drif t is energization due to drift betatron, and E chaos energization during chaotic motion.
to the peak ofB z . In general the gyrobetatron term is positive when the electrons are ahead of the peak iñ B z , and negative when they are on the tail of the peak. The reason is that B t is positive ahead of the peak and negative on the tail. The overall contribution of this term is in general less than that of the drift betatron term, and it can even be negative when the electrons are trailing the peak of the magnetic field. Figure 9 shows the ensemble average and energy of the injected electrons where electrons are injected at the enhancedB z . < > is constant showing that the injection is mainly adiabatic for this case. The result shows that the energization comes mainly from the drift betatron term (red), whereas the gyrobetatron term (yellow) is negative. Since the drift betatron term is larger and positive, there is a net gain in energy of electrons as they are injected. The gyrobetatron term is negative because during most of the injection, these electrons were behind the peak inB z and the B t is negative. At the same time there is a positive q ⃗ E · ⃗ v D as the electron drift direction is along the electric force. The later term dominates so that there is a net increase in energy. This can also be interpreted as that the ΔB seen by electrons is positive, as a result of the spatial gradient of B, even though the temporal ( B t ) change is negative. Figure 10 shows the results when the 10-keV electrons are initiated at the minimum B location; < > is not conserved for this case. For the first ≈ 100 s, the electrons were moving in the chaotic region so that is not conserved. At the same time, electrons were ⃗ E × ⃗ B drifting to the enhancedB z region. This drift is directed to the enhancedB z region because of ⃗ E × ⃗ B drift (the flow velocity perpendicular to ⃗ B) at the minimum B relative to the flow at the enhancedB z , which is in the direction of the enhancedB z . After the electrons move out of the minimum B region to the enhancedB z region, becomes constant. This figure also shows that the energization is mainly from drift betatron acceleration (red), even though is not conserved. The gyrobetatron (yellow) term is negative during this time as electrons are behind the peak ofB z so that B t is negative. The energization for this case is higher than the case when the electrons were initiated at the enhancedB z . This result is as expected from adiabatic energization, even though is not conserved.   Figure 11 shows the result of the case with initial energy of 100 keV initiated at the enhancedB z . The result shows that the energization is coming mainly from drift betatron acceleration (red). Unlike the 10-keV case, the gyrobetatron term (yellow) is positive. This is because these electrons stayed mainly ahead ofB z so that B t is positive; see the supporting information Movie S4; is conserved for most of the injection farther out, but at later stage, the conservation of is violated. The reason for the violation of conservation of at the later stage is because electrons move through the edge of the enhancedB z . In this region, the field is highly curved compared to the central region of the enhancedB z . Thus, electrons acquire chaotic motion in these regions; is conserved for the 10-keV initial energy at a similar stage in the injection as their gyroradii ( ) are smaller compared to the 100-keV initial energy case. Figure 12 shows electrons initiated with 100 keV in the minimum B region. Like the 10-keV case, electrons drift toward enhancedB z , while undergoing chaotic motion. This behavior is followed by adiabatic transport due to ⃗ E × ⃗ B drift until the electrons reach megaelectron volt energy range.
The results for the 3-keV electrons shown in the supporting information Movies S7-S9 are similar to the 10-keV case. The only difference occurs when the electrons were initiated at the minimum B region where the fractional change in is smaller. This is because 3-keV electrons have smaller gyroradii so that they are more adiabatic.
To summarize, electrons are mainly injected by ⃗ E × ⃗ B drift. Electrons are mainly energized by drift betatron acceleration as they drift in the direction of the electric force (the second term of the right-hand side of equation (2). The gyrobetatron energization (the first term of the right-hand side of equation (2) is also significant, but depending on the relative position of the electrons from the peakB z , it too can be negative. The overall importance of the gyrobetatron term is smaller than that of the drift betatron acceleration, contrasting the first and the second terms in equation (2). Electrons can also undergo chaotic motion that violates the conservation of . Chaotic motion occurs in a region of high curvature and low B value. Electron motion  becomes more chaotic when the energy increases. The chaotic motion, which happens when the electrons are initiated from a low B value, causes the energization to be nonadiabatic. For example, 10-keV electrons initiated from the minimum B region and injected to geosynchronous orbit have shown a decrease in < > by a factor of ≈ 0.6 as seen in Figure 10. This implies that these electrons have less energy in the perpendicular direction than what is expected from adiabatic motion. Thus, they move to the stronger magnetic field region without the constraint of the energy requirement by adiabatic motion. This effect is stronger for higher-energy electrons as can be seen from a decrease of < > by a factor of >5 for 100-keV initial-energy electrons in the first minute. The total energy is also smaller than what is expected from adiabatic energization, as it is smaller than the expected energy in the perpendicular direction. Thus, the energization by BBF cannot be completely described by the adiabatic theory.
MHD does not resolve oscillations with a frequency higher than the ion gyrofrequency. As a result high frequency waves, which could resonantly interact with the gyromotion of the electrons are absent in these simulations. Since these waves interact with the gyromotion like the gyrobetatron energization, inclusion of such wave effects in the simulations could change the importance of local acceleration. By local acceleration, we mean the acceleration due to time change of the fields or waves as compared to the drift betatron, which needs particles to drift to a stronger magnetic field region.

Chaotic Electron Motion
The conservation of can be violated when the fields vary on the time scale of the gyroperiod or when the gyroradius ( ) is comparable to either of L, R c , or L s , the magnetic gradient, curvature, and shear scale lengths. The nonadiabatic motion of charged particles in highly curved magnetic fields has been previously studied (Birmingham, 1984;Buchner & Zelenyi, 1989;Chen, 1992;Chirkov, 1978;Hudson et al., 1997;Speiser, 1991;Young et al., 2002;Engel et al., 2016). Nonadiabatic motion of particles can also exist in sheared magnetic fields (Smets, 2000;Pfefferle et al., 2015) or for short scale length Northrop (1963). Most of the study of chaotic motion was for the highly stretched field of the magnetotail. However, highly curved fields can also exist in a localized region in the vicinity of the BBFs. We demonstrate that these regions can cause nonadiabatic motion of energetic electrons for energy from as low as a few kiloelectron volts; see Eshetu et al. (2018). As an example, we traced 3-keV electrons near the BBF coming from the duskside of the magnetotail by solving the Lorentz equation. The electrons were started with a pitch angle of 45 • with randomly chosen gyrophase; and the corresponding versus time are shown in Figure 13. We plot because we found that the curvature effect is the main reason for chaotic motion in this case. The result shows that changes appreciably every time the electrons move through the magnetic equator, where is close to 1. We refer to this as a jump in . The jump of (Δ ) is sensitive to initial condition, and it even depends on the gyrophase. It is also random in which the value of Δ every jump is uncorrelated. All are qualitative features of the chaotic trajectories described by Buchner and Zelenyi (1989).

Pitch Angle Dependence and Evolution
The evolution of the pitch angle distribution depends on whether the energization is mainly in the parallel or perpendicular direction to the magnetic field. The gyrobetatron term energizes only in the perpendicular direction, whereas the drift betatron acceleration energizes both in parallel and perpendicular directions Northrop (1963), due to electric force acting along the ΔB drift direction resulting in perpendicular energization and along the curvature drift direction in to parallel energization (Birn et al., 2013). In addition, the perpendicular energy and the parallel energy are mediated by the mirror force. However, by looking at the pitch angle distribution at the magnetic equator, this effect can be neglected. Thus, the distribution in pitch angle at the equator is a result of perpendicular and parallel accelerations by betatron and Fermi acceleration of type B (Northrop, 1963), hereafter referred to as just Fermi acceleration.
A movie of the evolution of the pitch angle of electrons for the case of 10 keV initiated at the enhancedB z is shown in the Movie S10 in the supporting information, corresponding to Movie S1 for the energization. We use cos( eq ) for the color bar. The reason is that the number of particles per range in the color bar is uniform when the distribution is isotropic. The initial distribution in pitch angle at the magnetic equator is isotropic. From the movie S10, it can be seen that the injection of electrons is clustered by pitch angle, as a result of different ΔB and curvature drifts with pitch angle. In addition, the time electrons spent in the BBF region depends on pitch angle as the BBF is localized in the z direction. From the movie, it can be seen that high pitch angle electrons mostly stay behind the peak of the enhancedB z trailing the low pitch angle electrons until the BBF reaches x ≈ − 8R E . Farther inward, the flow becomes more complicated; however, the different pitch angle electrons are clustered in different regions. Looking at the movie, at the end of the simulation (05:38 ST), it can be seen that electrons that drift farthest eastward are high pitch angle electrons and energized to ≈ 60-70 keV. There are also electrons that are energized as much but injected later. These electrons are seen around y ≈ 5 inside of geosynchronous orbit, and they have low pitch angle. The high pitch angle electrons that are seen farthest eastward are energized by the betatron mechanism, whereas the low pitch angle electrons are energized mainly by the Fermi mechanism. In order to confirm this, we included a movie showing the ratio of the change of energy in the perpendicular to parallel direction ( Movie S11 in the supporting information. Note that if ΔW ⟂ or ΔW || are negative value, we use saturated red and blue, respectively. Looking at Movie S11 at 05:38 ST, the high pitch angle electrons are mainly energized in the perpendicular direction, whereas the low pitch angle electrons are energized in the parallel direction consistent with the betatron and Fermi accelerations.
A movie of the evolution of the pitch angle of electrons for the case of 10-keV initial-energy electrons initiated at the minimum B locations can be seen in Movie S12 in the supporting information. From the movie, it can be seen that electrons are clustered with pitch angle. The electrons with the highest pitch angle drift tailward because of higher ΔB drift; see, for example, at 05:23:47 ST. The electrons trailing the peakB z have the lowest pitch angle. Even though the curvature drift is larger for small pitch angle electrons, they did not drift tailward like the larger equatorial pitch angle electrons. The reason is that if the pitch angle is too small, the electrons stay in the BBFs for shorter time as the BBFs are localized in latitude. Thus, the integrated curvature drift is not large. For these electrons, the integrated ⃗ E × ⃗ B drift is also less than the propagation speed of the BBFs so that they trail the peak in the enhancedB z .
The electrons that drift toward the enhancedB z have higher pitch angle than those trailing the peak. For these electrons, the ⃗ E × ⃗ B drift is larger than those with lowest pitch angle. In addition, their ⃗ E × ⃗ B dominates over ΔB and curvature drifts. For x < −10R E , ahead (behind) the enhancedB z the electrons are mainly high (low) pitch angle. We did not use ΔW ⟂ ΔW || to see if the energization is betatron or Fermi, because the electrons undergo chaotic motion where the equatorial pitch angle can change without energization (Eshetu et al., 2018). The relative position of the high and the low pitch angle electrons relative to the enhancedB z is contrary to the case of electrons initiated at the enhancedB z . This is because the complicated structure of the flow also changes direction, and a study of this convection pattern is outside the scope of this work.
In general, the Fermi and betatron accelerations are location dependent relative to the magnetic structure of the BBF in addition to pitch angle. A more detailed analysis of the importance of these two terms and the evolution in pitch angle for a realistic initial distribution in space and energy will be deferred for the future.

Discussion and Conclusions
We have shown that when the electrons are initiated at the location of BBFs where the magnetic field is localized in dipolarization flux bundles (Gabrielse et al., 2016), they mainly drift inward due to ⃗ E × ⃗ B. These electrons conserve if their initial energy is 10 keV or lower. The energization in this case is mainly from drift betatron acceleration. The gyrobetatron term, the first term in equation (2), is also significant, but this contribution depends on where the particle is located relative to the peakB z while injected. When the electrons trail the peak ofB z , the gyrobetatron term is negative, and when they are ahead of the peak, it is positive. Thus, the overall contribution of the gyrobetatron acceleration is small. The total energy gain in this case depends on the ratio of the magnetic field at the injection location to the initial value. Since electrons were initiated in the stronger magnetic field of the BBF, the energization for this case is limited.
When electrons are initiated at the minimum B location on the duskside of the BBF, some of the electrons move tailward because of ΔB drift. The rest move earthward while also moving into the enhancedB z . Electrons that moved into the enhancedB z are eventually injected as when they are initiated there. The motion of the electrons while bouncing through the minimum B region does not conserve . Even though the electrons are mainly accelerated by the drift and gyrobetatron terms, their chaotic motion causes to change significantly.
When the electrons are started ahead of the BBF structure, electrons can easily ΔB drift around the enhanced B z region. There is also ⃗ E × ⃗ B drift, which pushes the electrons ahead of it. The curvature and the ΔB drift dominate for all initial energies simulated, and its effect is larger for higher energy. Electrons can move around the enhancedB z , and some of them can also move into the enhancedB z region. The electrons that move into the enhancedB z are injected by ⃗ E × ⃗ B drift like the other cases. But the efficiency of the BBF in injecting the electrons is weaker than in the other cases.
In general, the electron motion can be described qualitatively with drift motion, even though is not strictly conserved. The injection is efficient when the electrons are located where the ⃗ E × ⃗ B drift dominates. ⃗ E × ⃗ B dominates when electrons are in the enhancedB z region and when electrons have lower energy. However, the energy gain is smaller for lower initial energy and when they are initiated in the enhancedB z region. The ΔB drift is important not only for energization but also for detraping electrons from the flow and injecting them into the inner magnetosphere. Electrons that are not energized enough to be detrapped after being injected will evolve with the complicated flow structure of the BBF. These electrons move inward up to ≈ 7R E , but they move back tailward after the flow brakes. Electron motion can be chaotic, and the conservation of can be violated. This happens in the highly curved, low B region, where the electron energy is high such that is close to 1. The violation of the conservation of is stronger farther away from the Earth and for higher energy.
Overall, Figures 2 -7 show that BBFs play a significant role in earthward transport of plasma sheet electrons, while Figures 9 -12 show how the electrons are energized from drift betatron and gyrobetatron accelerations. The energization in chaotic motion of electrons in the BBF sampled in Figures 2 -7 is small compared to the the gyro and drift betatron acceleration terms of equation (2). The pitch angle distribution of the injection and effect of chaotic motion on pitch angle scattering and loss will be examined separately.