Particle energization due to magnetic reconnection is an important unsolved problem for myriad space and astrophysical plasmas. Electron energization in magnetic reconnection has traditionally been examined from a particle, or Lagrangian, perspective using particle-in-cell (PIC) simulations. Guiding-center analyses of ensembles of PIC particles have suggested that Fermi (curvature drift) acceleration and direct acceleration via the reconnection electric field are the primary electron energization mechanisms. However, both PIC guiding-center ensemble analyses and spacecraft observations are performed in an Eulerian perspective. For this work, we employ the continuum Vlasov–Maxwell solver within the Gkeyll simulation framework to reexamine electron energization from a kinetic continuum, Eulerian, perspective. We separately examine the contribution of each drift energization component to determine the dominant electron energization mechanisms in a moderate guide-field Gkeyll reconnection simulation. In the Eulerian perspective, we find that the diamagnetic and agyrotropic drifts are the primary electron energization mechanisms away from the reconnection x-point, where direct acceleration dominates. We compare the Eulerian (Vlasov Gkeyll) results with the wisdom gained from Lagrangian (PIC) analyses.

## I. INTRODUCTION

Magnetic reconnection is a ubiquitous process in space and astrophysical plasmas, and it plays a fundamental role in transforming stored magnetic energy into particle kinetic and thermal energies. Reconnection is thought to produce high energy, non-thermal particles in a variety of systems, including gamma ray bursts,^{1} solar and stellar flares,^{2} and pulsar magnetospheres.^{3} Non-thermal electrons are responsible for producing a significant fraction of the emitted light from high-energy astrophysical objects; therefore, understanding how these non-thermal electrons are produced is fundamentally important. Solar observational results suggest that reconnection very efficiently produces non-thermal electrons.^{4,5} The process often associated with the acceleration of the electrons is Fermi acceleration due to the strongly curved, contracting magnetic island structures that form in the outflow of magnetic reconnection.^{6,7} These outwardly flowing magnetic structures are believed to lead to local enhancements of curvature drift acceleration. Similar structures are observed in relativistic turbulence simulations^{8,9} and are also sites of enhanced thermal and non-thermal electron production.

In addition to non-thermal electron production, reconnection also leads to significant thermal energization of both ions and electrons. Heating due to reconnection likely plays an important role in the heating of the solar corona and acceleration of the solar wind.^{10–12} Electron heating due to small-scale reconnecting current sheets in turbulence plays a fundamental role in dissipating energy injected at large scales in myriad space and astrophysical plasmas.^{13–16} Recent spacecraft missions like Cluster,^{17} THEMIS,^{18} and the magnetospheric multiscale (MMS) mission^{19} have also identified significant electron heating (and particle acceleration) at reconnection sites in the magnetosheath^{20,21} and in the magnetotail.^{22–25}

Whether thermal or non-thermal, electron energization has traditionally been examined from a particle, or Lagrangian, perspective using particle-in-cell (PIC) simulations. Guiding-center based analyses of the PIC simulations have suggested that Fermi (curvature drift) acceleration and direct acceleration via the reconnection electric field are the primary electron energization mechanisms.^{26–30} Due to its spatially limited extent in the vicinity of the reconnection x-point or x-line, direct acceleration is generally believed to play a sub-dominant role in reconnection, except when the guide field is much larger than the reconnection field.^{31,32}

The guiding-center analyses are traditionally based on the work of Northrop,^{33} who semi-rigorously derives the energy evolution equation for a single electron in the non-relativistic limit, averaged over a gyration period. However, the derivation includes only the $ \u2207 B$ and curvature drifts, while excluding other single-particle drifts. The single-particle energization equation is then summed over local ensembles of guiding-centers to produce a fluid-like energization equation for the particle moments, i.e., $ d E e / d t = j e \xb7 E$, where $ j e$ is the electron current density, **E** is the electric field, and *E _{e}* is the total electron kinetic energy. Li

*et al.*

^{29}effectively extend the number of drifts included in the guiding-center model by adding drifts to the electron current following derivations by Parker

^{34}and Blandford

*et al.*

^{35}for the electron drift contributions to the current in a pressure gyrotropic plasma. Li

*et al.*

^{30}add additional drift contributions due to pressure agyrotropy. In all of these studies, guiding-center ensemble averaged analyses of the PIC results indicate curvature drift is generally the dominant guiding-center energization mechanism. Interestingly, Li

*et al.*

^{36}perform a different separation of the drift energization mechanisms into compressional and shear (curvature drift accounting for pressure anisotropy) and find that compressional energization dominates for weak-to-moderate guide fields ( $ B g \u2264 0.2 B 0$), while for larger guide fields, the compressional and shear terms contribute equally.

However, PIC simulations follow individual particles rather than guiding centers, so particle energization terms due to the motion of the particles relative to their guiding-center rest frames are absent from existing guiding-center analyses. Furthermore, by performing local ensemble averages and then examining particle moments, one has moved from a Lagrangian to an Eulerian perspective. Therefore, existing kinetic studies to identify the dominant electron energization mechanism in reconnection miss energization terms due to being performed in the particle rather than guiding-center limit, and the studies mix Lagrangian and Eulerian perspectives. Formally, in the Eulerian perspective, one is no longer directly quantifying particle energization, rather the Eulerian perspective provides information about the particle distribution function, as well as bulk flow and thermal energization. To date, the only fully Eulerian reconnection simulations focused on identifying the primary energization mechanism have been performed using magnetohydrodynamics (MHD). The MHD simulations examining the importance of magnetic curvature energization on bulk flow energization arrive at conflicting conclusions, finding that curvature dominates the energization^{37} and that curvature and magnetic pressure expansion contribute approximately equally.^{38}

Recent phase-space diagnostics like the field-particle correlation^{39,40} (FPC) have been developed from an Eulerian perspective to identify the signatures of particle energization due to mechanisms such as Landau^{41,42} and cyclotron^{43} damping, direct acceleration in reconnection,^{32} and shock drift acceleration.^{44,45} However, the initial identification of a particular FPC signature requires knowledge of the active energization mechanism at a specific configuration point. Therefore, it is essential to derive an Eulerian picture of the energization to interpret velocity-space signatures provided by the FPC. Development and tabulation of these signatures will enable single-spacecraft observations to explore the nature of the electron energization in reconnection, turbulence, and shocks.

To address this dichotomy between Lagrangian and Eulerian perspectives of electron energization in reconnection, we use the continuum (Eulerian) Vlasov–Maxwell solver within Gkeyll to perform 2D-3V kinetic simulations of magnetic reconnection with a moderate guide field. We identify the dominant energization mechanisms locally and globally in an Eulerian simulation by employing an Eulerian-based analysis, including a rigorously complete description of all fluid drifts: diamagnetic, curvature, polarization, and agyrotropic. The results herein are straightforwardly generalizable to 3D for use in other simulation or spacecraft data analysis.

In Sec. II, we derive and discuss the physical significance of the guiding-center and Eulerian energization equations, including the presentation of alternative formulations of the Eulerian energization. We describe the continuum kinetic simulation code, Gkeyll, and initial conditions in Sec. III. In Sec. IV, we present the results of the Gkeyll simulations. Finally, in Sec. V, we discuss the implications and significance of the results.

## II. ENERGIZATION MODELS

### A. Guiding-center energization

*et al.*

^{26–28}and Northrop

^{33}to construct the energy evolution equation for a single electron in the non-relativistic, guiding-center limit averaged over a gyration period,

*q*, and

_{e}*m*are the electron charge and mass, $ \mu = m e v \u22a5 2 / 2 B$ is the magnetic moment, $ v c$ and $ v g$ are the curvature and grad-B drifts,

_{e}*ε*correspond to the guiding-center kinetic energy, while

*μB*is related to the average energy of rotation about the guiding center.

*et al.*

^{6}

Before describing the drifts in an Eulerian perspective, we note that this guiding-center approach is somewhat ambiguous because only the curvature and $ \u2207 B$ drifts have been included. Li *et al.*^{29,30} ameliorate this shortcoming by including additional drifts, some of which they find to be particularly important in the weak or zero guide field limits. However, regardless of the number of drifts included, Eq. (4) only accounts for the kinetic energy of the guiding centers plus the particle energy averaged over the gyroperiod.

### B. Eulerian energization

*s*to find the bulk drifts

^{46}and $ \Pi s a$ is the remaining, agyrotropic portion of the pressure tensor. Finally, to compute the drifts, one must cross Eq. (5) with

**B**. After manipulation (see Appendix C of Juno

*et al.*

^{44}for details),

*et al.*,

^{44}the RHS terms of Eq. (10) are referred to as the $ E \xd7 B$, diamagnetic, curvature, agyrotropic, and polarization drifts, respectively. We note that while there is some ambiguity in this nomenclature, Eq. (10) is a complete description of all perpendicular bulk motion in the plasma.

^{47}for a full derivation). Despite this physical difference in reference frames, we can re-organize the terms in Eq. (11) to more closely resemble the guiding-center derivation:

*W*now contains all of the $ p \u22a5$ dependence, including curvature-related components. To make the comparison more direct, we can also further manipulate the guiding-center drift energization equation [see Eq. (A3) in Appendix].

_{beta}To aid in following the three different groupings of the curvature, betatron, and diamagnetic drift terms, Fig. 1 provides a visual illustration of the groupings. The color-coding corresponds to that used in Fig. 4.

*W*contributes little to the electron energization. Since

_{pol}*W*is small, this term is itself either small or mostly canceled by other terms appearing in

_{pol}*W*. Therefore, we do not separate it when discussing Eulerian energization.

_{pol}## III. SIMULATION SETUP

To study electron energization in magnetic reconnection, we employ the fully kinetic continuum Vlasov–Maxwell solver portion of the Gkeyll framework.^{48,49} Gkeyll directly evolves the full Vlasov–Maxwell system in up to 3D-3V, three configuration and three velocity space dimensions, for an arbitrary number of species using the discontinuous Galerkin discretization scheme in all phase space dimensions.

^{50}Therefore, we initialize a Harris sheet with the following magnetic field and density profiles:

*T*and

_{i}*T*, satisfying $ T i / T e = 5$. To maintain electron magnetization throughout the layer,

_{e}^{51}we add a uniform, moderate guide field, $ B g = 0.1 B 0 z \u0302$, to the above magnetic field profile. Pressure balance requires that $ n 0 ( T i + T e ) = B 0 2 / 2 \mu 0$, which implies $ \beta i = 2 \mu 0 n 0 T i / B 0 2 = 5 / 6$. We choose $ n 0 / n \u221e = 5$, a reduced mass ratio, $ m i / m e = 25$, and speed of light, $ v A e / c = 0.25$, where $ v A e = B 0 / \mu 0 m e n e$ is the in-plane, upstream, electron Alfvén speed. The layer thickness,

*λ*, is taken to be $ \lambda = 0.5 d i$, where $ d i = v A i / \Omega c i$ is the ion inertial length. The simulation domain is chosen to be $ \u2212 L x / 2 \u2264 x \u2264 L x / 2 , \u2212 L y / 2 \u2264 y \u2264 L y / 2 , \u2212 6 v t h s \u2264 v s \u2264 6 v t h s$, where $ L x = 2 L y = 8 \pi d i , \u2009 v t h s = 2 T s / m s$, and the same velocity space extents are used in each species' three velocity dimensions, except we extend the velocity space extents of the electrons in the guide field direction to $ \u2212 8 v t h e \u2264 v z , e \u2264 8 v t h e$ to resolve additional particle acceleration. Periodic boundary conditions are employed in the

*x*direction and conducting (reflecting) walls for the fields (particles) are used in the

*y*direction, while we employ zero-flux boundary conditions in velocity space. $ ( n x , n y ) = ( 112 , 56 )$ configuration space cells and $ n v 3 = 24 3$ velocity space cells are used for the protons and electrons with piecewise quadratic Serendipity polynomial basis functions

^{52}to discretize the fields and particle distribution functions. We also employ a conservative Lenard–Bernstein collision operator

^{53,54}with electron–electron collision frequency $ \nu e = 0.0004 \Omega c e$ and ion–ion collision frequency $ \nu i = \nu e / m i / m e$. The collision frequencies are chosen to satisfy the weakly collisional limit while being sufficient to regularize velocity space and maintain numerical stability.

*x*and

*y*magnetic field modes are randomly perturbed with a root mean square amplitude $ B pert RMS = 0.01 B 0$.

## IV. SIMULATION RESULTS

### A. Overall evolution

In Fig. 2, we present the evolution of the out-of-plane current density, *j _{z}*, in color along with contours of the out-of-plane vector potential,

*A*, at $ t \Omega c i = 4 , 9 , 14 ,$ and 19 in panels (a)–(d), respectively. In Fig. 3, we present the reconnection rate (a), integrated and normalized change in energy components (b), and the $ j \xb7 E$ work (c). The components of the spatially integrated and normalized change in energy in panel (b) are defined as: $ E B = B 2 / 2 \mu 0$ is the magnetic energy, $ E E = \epsilon 0 E 2 / 2$ is the electric field energy, $ E u = m n u 2 / 2$ is the flow energy, $ E T = Tr ( P ) / 2$ is the thermal energy, $ E tot = E B + E E + \u2211 s ( E T s + E u s ) , \u2009 E 0 = E tot ( t = 0 )$, and $ \Delta E * ( t ) = E * ( t ) \u2212 E * ( t = 0 )$. In panels (b) and (c), red (blue) lines correspond to ion (electron) quantities. From these panels, we can see that the electron work peaks at $ t \Omega C i \u2243 19$, shortly after the peak reconnection rate. We will focus on this time for our local analyses of electron energization.

_{z}### B. Electron energization

Following Refs. 26, 29, and 30, we integrate over the entire domain to examine the drift energization terms and their sums in Eqs. (4), (11), and (12), which are plotted in Fig. 4. The guiding-center energization terms [Eq. (4)] plotted in panel (a) display the expected results based on PIC simulations: the curvature drift (red) is large, positive, and the primary contributor to the sum (dashed black); the $ W beta \u2212 gradB$ (blue) is negative due to the conversion of magnetic energy into particle energy; and *W _{par}* (green) is positive but small. However, the sum and the directly computed $ j e \xb7 E$ (magenta) track each other well but display quantitative differences of up to 20%, suggesting that non-negligible energization terms may be missing.

Moving to the Eulerian energization formulations in Eqs. (11) and (12) in panels (b) and (c) respectively, we first note that the sum (black dashed) and $ j e \xb7 E$ (purple) agree very well at all times, indicating that all significant energization is included, as expected for a rigorous accounting of all bulk flow motion. Next, we note that the diamagnetic drift (lime), *W _{diam}*, is the primary contributor to the domain integrated work in panel (b). For the base Eulerian energization formulation, Eq. (11), the curvature drift (orange), $ W curv 1$, is at all times negative in panel (b). The agyrotropic (cyan) and polarization (yellow) drifts as well as

*W*contribute little to the total. Mirroring the guiding-center formulation using Eq. (12) and splitting the curvature drift, $ W curv 1$, into parallel curvature drift, $ W curv 2$ (magenta) and betatron acceleration,

_{par}*W*, (gray) in panel (c), we now see that this form of the curvature drift is large and positive, dominating all other drifts, except at late times. In panel (d), we repeat the decomposition found in panel (b), except we separate the diamagnetic ( $ W diam = u E \xb7 \u2207 p s , \u22a5$) and curvature-like betatron (pink), $ W beta \u2212 curv = \u2212 p s , \u22a5 u E \xb7 \kappa $, terms (see Fig. 1 for a visual illustration of the groupings of curvature, betatron, and diamagnetic drifts). Here, we see clearly that the curvature-like betatron term (pink), $ W beta \u2212 curv$, and $ W curv 2$ (magenta) have similar magnitudes and opposite signs, with $ W beta \u2212 curv$ being larger for most of the simulation, leading to $ W curv 1 = W curv 2 + W beta \u2212 curv$ being small and negative in the panel (b).

_{beta}To make the comparison between drift models and their terms more quantitative, Table I presents the same domain integrated values as Fig. 4 but focusing on $ t \Omega c i = 19$. Each term of the different energization formulations is presented separately, normalized to $ j e \xb7 E$, and the sum column entries refer to the sum of all contributing terms in that row. Here, we can see more clearly that the guiding-center formulation overestimates the energization by $ \u223c 17$%, while the fluid energization approaches are equivalent and accurate to within 1%. We can also see that the diamagnetic drift contributes $ \u223c 90$% of the overall energization. Finally, using Eq. (13), it is clear that $ W curv 2$ and $ W beta \u2212 curv$ almost cancel each other, leaving a remainder of $ W curv 1 \u2243 \u2212 0.11$, an overall loss of energy due to the fluid curvature drift.

Equation . | $ j e \xb7 E$ . | W
. _{par} | $ W beta \u2212 gradB$ . | $ W curv 0$ . | W
. _{diam} | $ W curv 1$ . | W
. _{agyro} | W
. _{pol} | $ W curv 2$ . | W
. _{beta} | $ W beta \u2212 curv$ . | Sum . |
---|---|---|---|---|---|---|---|---|---|---|---|---|

W_{0} [Eq. (4)] | 1 | 0.0484 | −0.4866 | 1.6062 | ⋯ | ⋯ | ⋯ | ⋯ | ⋯ | ⋯ | ⋯ | 1.1681 |

W_{1} [Eq. (11)] | 1 | 0.0484 | ⋯ | ⋯ | 0.8959 | −0.1144 | 0.1136 | 0.0644 | ⋯ | ⋯ | ⋯ | 1.0079 |

W_{2} [Eq. (12)] | 1 | 0.0484 | ⋯ | ⋯ | ⋯ | ⋯ | 0.1136 | 0.0644 | 1.1971 | −0.4157 | ⋯ | 1.0079 |

W_{3} [Eq. (13)] | 1 | 0.0484 | ⋯ | ⋯ | 0.8959 | ⋯ | 0.1136 | 0.0644 | 1.1971 | ⋯ | −1.3116 | 1.0079 |

Equation . | $ j e \xb7 E$ . | W
. _{par} | $ W beta \u2212 gradB$ . | $ W curv 0$ . | W
. _{diam} | $ W curv 1$ . | W
. _{agyro} | W
. _{pol} | $ W curv 2$ . | W
. _{beta} | $ W beta \u2212 curv$ . | Sum . |
---|---|---|---|---|---|---|---|---|---|---|---|---|

W_{0} [Eq. (4)] | 1 | 0.0484 | −0.4866 | 1.6062 | ⋯ | ⋯ | ⋯ | ⋯ | ⋯ | ⋯ | ⋯ | 1.1681 |

W_{1} [Eq. (11)] | 1 | 0.0484 | ⋯ | ⋯ | 0.8959 | −0.1144 | 0.1136 | 0.0644 | ⋯ | ⋯ | ⋯ | 1.0079 |

W_{2} [Eq. (12)] | 1 | 0.0484 | ⋯ | ⋯ | ⋯ | ⋯ | 0.1136 | 0.0644 | 1.1971 | −0.4157 | ⋯ | 1.0079 |

W_{3} [Eq. (13)] | 1 | 0.0484 | ⋯ | ⋯ | 0.8959 | ⋯ | 0.1136 | 0.0644 | 1.1971 | ⋯ | −1.3116 | 1.0079 |

Next, we look at cuts of the energization terms along the midplane at *y* = 0 at the time of maximum energization ( $ t \Omega c i = 19$). In Fig. 5, we present these cuts using the same four drift formulations presented in Fig. 4. First, at the x-point, $ x \u2243 1 d i$, all of the energization is due to the parallel acceleration term (green), *W _{par}*. Using the guiding center formulation Eq. (4) (a), we see that away from the x-point, the curvature drift term is the largest contributor; however, we also see that the sum (dashed black) does not well agree with the directly computed $ j e \xb7 E$ at any point along the cut. Using the Eulerian description Eq. (11) (b), we see that the largest positive contributor is the agyrotropic drift (cyan), and the curvature drift (orange) is generally large and negative at the same position in

*x*. Using the alternative Eulerian description Eq. (12) (c), again the curvature drift (magenta) is generally large and positive but is mostly balanced by an even larger, negative contribution from the betatron term (gray). In panel (d), we again repeat the decomposition found in panel (b), except we separate the diamagnetic (lime) and curvature-like betatron terms (pink). In all cases, we note that at most points away from the x-point, several drifts contribute simultaneously with comparable amplitudes, with the agyrotropic drift generally dominating and the curvature drift ( $ W curv 1$) being generally negative.

Finally, we examine the full spatial distribution of *W* and its components at $ t \Omega c i = 19$ in Fig. 6 with black contours of the out-of-plane vector potential, *A _{z}*, overplotted for reference. Superscript values indicate the drift energization equation in which the term appears. Due to the wide range of values, we employ a mixed log-linear color bar which is logarithmic for values $ | W | \u2265 1 \xd7 10 \u2212 4$ and linear for values $ | W | < 1 \xd7 10 \u2212 4$. From these panels, we can see: (i) $ W \u2225$ is large and positive in the vicinity of the x-point; (ii)

*W*is mostly positive and is the most space-filling term; (iii)

_{diam}*W*is large and positive in the entire diffusion region; (iv) $ W curv 1$ is mostly negative; (v) $ W curv 2$ is the largest positive contributor in the downstream region and is generally negated by; and (vi)

_{agyro}*W*, which is the largest negative contributor, essentially canceling $ W curv 2$ everywhere. In Fig. 7, we present the sums of the guiding center drifts,

_{beta}*W*

_{0}, and the Eulerian drifts,

*W*

_{1}in the top and middle panels. For comparison, we again plot $ j e \xb7 E$.

*W*

_{0}is entirely dominated by $ W curv 0$ and does not follow the actual $ j e \xb7 E$ energization, whereas

*W*

_{1}is a nearly perfect match to the electron work.

## V. DISCUSSION AND CONCLUSIONS

The guiding-center description of electron energization, Eq. (4), clearly shows that curvature drift, $ W curv 0$, is the dominant energization mechanism in the reconnection simulation, both globally and locally, away from the x-point. However, neither globally nor locally does the sum of the energization mechanisms agree well with $ j e \xb7 E$, suggesting that the Lagrangian description, unsurprisingly, does not work well in an Eulerian simulation. The disagreement can be partially ameliorated by including additional guiding-center drifts in Eq. (1), e.g., polarization, and properly transforming from guiding-centers to the fluid frame.

On the other hand, the Eulerian energization perspective, Eq. (11), is a rigorously complete description of all bulk-drift motion. This perspective suggests that the diamagnetic drift, *W _{diam}*, dominates globally with large local spikes of agyrotropic drift,

*W*, energization. We note that both of these drifts are unique to the Eulerian perspective and have no direct equivalent in a guiding-center model, even when summed over ensembles of guiding centers. Despite being rigorously complete, the nomenclature and grouping of drift terms in Eq. (11) is somewhat ambiguous, as demonstrated by the additional formulations in Eqs. (12) and (13). We also note that the diamagnetic drift can be separated into $ \u2207 B$ and magnetization drifts, c.f., Appendix F of Juno

_{agyro}*et al.*

^{44}

In the Eulerian perspective, the curvature drift, $ W curv 1 = W curv 2 + W beta \u2212 curv$, manifestly depends on the pressure anisotropy, differing significantly from the guiding-center description, which depends on only $ p \u2225$. The dependence on the pressure anisotropy leads to significant local and global cancellation, with $ W curv 2 = p \u2225 u E \xb7 \kappa $ and $ W beta \u2212 curv = \u2212 p \u22a5 u E \xb7 \kappa $ nearly canceling each other everywhere in the domain. The strong negative correlation between the two components suggests that they should be combined as $ W curv 1$ rather than separated. Although $ W curv 1$ is generally negative in the simulation, clearly indicating that the curvature drift leads to a net loss of energy, that does not mean that high energy particles are not being produced in the Eulerian picture. Since the Eulerian perspective provides information about bulk energization, it is possible for $ W curv 2$ to produce an electron power-law tail leading to a bulk drift while $ W beta \u2212 curv$ produces a greater loss of bulk flow energy from the core of the distribution, resulting in an overall loss of energy due to the curvature drift. Furthermore, we note that in the pressure isotropic limit, the fluid curvature drift mostly vanishes—the component of the curvature drift proportional to $ u \u2225 2$ that formally appears in the polarization drift remains in the pressure isotropic limit. Indeed, this cancellation of curvature-drift terms in a pressure isotropic fluid is even noted in Goldston and Rutherford's textbook.^{47}

The guiding-center formulation can be manipulated to more closely resemble the Eulerian description, including the appearance of the $ W beta \u2212 curv = \u2212 p \u22a5 u E \xb7 \kappa $ term in Eq. (A3). However, this manipulation is not well motivated physically and obfuscates the relevant guiding-center physics by re-expressing betatron acceleration as a series of other terms. The re-expression of the betatron term does, however, serve to highlight the problematic nature of attempting to apply a guiding-center description to an Eulerian system, whether that Eulerian system is a simulation or the solar wind. Because Eqs. (4) and (A3) are derived by summing over ensembles of guiding centers, fluid moments like pressure and bulk flow arise, which provides a false impression that these guiding-center descriptions can be connected directly to fluid-like bulk evolution of the plasma [Eqs. (4) and (A3)], but that is not simply the case. Terms like the diamagnetic drift and agyrotropic drift will always be absent from any formal derivation based on guiding-center motion since guiding-center (or even particle) descriptions cannot account for bulk effects. Furthermore, the guiding-center energization equation, Eq. (1), includes only the motion (drift) of the guiding center plus an induction term ( $ \mu \u2202 B / \u2202 t$) that is due to the electric field acting on the gyroperiod integrated motion of an electron, which is not the complete motion of the particle. The motion of the particle in the reference frame in which the guiding-center is at rest is not included in Eq. (1), but to connect to an Eulerian perspective, this motion must be included. When this motion is included,^{47} one arrives at the Eulerian description described herein for terms like the curvature drift.

So, where does this leave us? Are the prior PIC results indicating that curvature drift leads to significant particle energization incorrect? The short answer is no, because the particle guiding centers are indeed energized (accelerated) by $ W curv 0$, meaning that where that term is positive and large, particle guiding centers are gaining significant energy. However, that is not the full story, because the analysis typically leading to that conclusion is incomplete and typically performed in an Eulerian framework by integrating over regions of particles—this fact applies to both PIC simulations and spacecraft data. In the Eulerian perspective, one must also include additional terms arising from the full particle motion, not just the guiding-center drift plus gyroperiod integrated motion as included in the guiding-center energization equation. When doing so, the full Eulerian curvature drift, $ W curv 1$, may indicate a very different scenario for the electrons. In the case of the simulation examined herein, $ W curv 1$ is relatively small and negative, whereas $ W curv 0$ is large, positive, and the dominant contribution to electron energization. Fundamentally, this contradiction arises because the Eulerian picture describes the generation of bulk flows rather than guiding-center acceleration described by Eq. (1). The confusion follows from attempting to directly relate guiding-center acceleration described by Eq. (1) to bulk quantities described by the guiding-center ensemble integrated Eq. (4). The guiding-center ensemble equations do not include the full particle motion, nor do they account for bulk plasma properties that lead to bulk drifts, e.g., diamagnetic and agyrotropic effects. Therefore, the guiding-center ensemble energization Eq. (4) can be useful for highlighting regions of the plasma in which electrons are being accelerated by the curvature drift, potentially leading to highly energetic electrons via Fermi acceleration;^{6,27,28} however, the guiding-center ensemble energization does not necessarily indicate regions in which electron bulk flows are large nor where electron heating (thermal broadening of the electron distribution) may be occurring. On the other hand, the Eulerian framework [Eq. (11)] only indicates where electrons are being energized in a bulk sense (bulk flows and heating), since in the Eulerian picture, there is no conception of single-particle or guiding-center motion.

In summary, the guiding-center approach can provide insights into high-energy particles, while the Eulerian approach gives a more complete perspective on the bulk energization and heating of the plasma. The Eulerian picture does describe the production of high energy particles, but only in a bulk sense as power-law tails of the distribution function, which in terms of moments, also lead to bulk flows. The guiding-center description of heating is incomplete, while the Eulerian approach provides no information about individual particle energization to high energies. We, thus, encourage a greater awareness of the desired goals when applying these diagnostics to the analysis of spacecraft data and simulations so that the proper diagnostic and perspective is applied. Finally, we note that field-particle correlations, by highlighting the different regions of velocity space in which particles are energized, may help to bridge this gap between the Lagrangian and Eulerian pictures.^{44}

## ACKNOWLEDGMENTS

The authors thank Hui Li for helpful discussions. This research is part of the Frontera computing project at the Texas Advanced Computing Center. Frontera is made possible by National Science Foundation (NSF) Award OAC-1818253. J.M.T. was supported by the NSF under Grant No. AGS-1842638. G.G.H. was supported by the NSF Grant AGS-1842561. J. Juno was supported by the U.S. Department of Energy under Contract No. DE-AC02-09CH1146 via an LDRD grant. The development of Gkeyll was partly funded by the NSF-CSSI program, Award No. 2209471.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Jason M. TenBarge:** Formal analysis (lead); Writing – original draft (lead); Writing – review & editing (equal). **James Juno:** Writing – review & editing (equal). **Gregory G. Howes:** Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are openly available in Zenodo at https://doi.org/10.5281/zenodo.10524773 (Ref. 55).

### APPENDIX: ALTERNATIVE GUIDING CENTER FORMULATION

*et al.*

^{36}However, we further note that $ \u2207 \xb7 ( p \u22a5 u E ) = p \u22a5 \u2207 \xb7 u E + u E \xb7 \u2207 p \u22a5$, which integrated over all space implies $ \u2212 \u222b R p \u22a5 \u2207 \xb7 u E d r = \u222b R u E \xb7 \u2207 p \u22a5 d r = \u222b R W diam d r$, assuming there is no net flux of thermal energy into or out of the domain. Therefore, in a spatially integrated sense, as often presented in guiding-center analyses, this compressional betatron-related term is related to the Eulerian diamagnetic drift energization. Thus, the electron guiding-center energization can be expressed approximately as

## REFERENCES

*The Adiabatic Motion of Charged Particles*

*Introduction to Plasma Physics*