A Field-Particle Correlation Analysis of a Perpendicular Magnetized Collisionless Shock

Using the field-particle correlation technique, we examine the particle energization in a 1D-2V continuum Vlasov--Maxwell simulation of a perpendicular magnetized collisionless shock. The combination of the field-particle correlation technique with the high fidelity representation of the particle distribution function provided by a direct discretization of the Vlasov equation allows us to ascertain the details of the exchange of energy between the electromagnetic fields and the particles in phase space. We identify the velocity-space signatures of shock-drift acceleration of the ions and adiabatic heating of the electrons due to the perpendicular collisionless shock by constructing a simplified model with the minimum ingredients necessary to produce the observed energization signatures in the self-consistent Vlasov-Maxwell simulation. We are thus able to completely characterize the energy transfer in the perpendicular collisionless shock considered here and provide predictions for the application of the field-particle correlation technique to spacecraft measurements of collisionless shocks.


Introduction
Shock-waves, disturbances propagating faster than the largest local wave speed, are ubiquitous in space and astrophysical plasmas. From supernova remnants to the Earth's bow shock, we observe a variety of plasma environments in which these shock-waves efficiently convert the bulk kinetic energy of their supersonic flows into other forms of energy, e.g., plasma heat, accelerated particles, and electromagnetic radiation. Critically, a commonality amongst this variety of plasma environments is that they are weakly collisional, i.e., the collisional mean free path is much larger than the relevant plasma length scales, such as the gyroradius. Thus, this energy conversion must be mediated by collisionless interactions. Here, collisionless energy transfer refers to the myriad of mechanisms for transferring energy between the particles in the plasma and the electromagnetic fields, or vice versa, from wave-particle resonances to instabilities.
Diagnosing collisionless energy transfer is a grand challenge in plasma physics, as the plasma has many routes at its disposal for converting energy from one form to another. For collisionless shocks, a number of processes have been identified as potential energy transfer mechanisms, and the efficiency of each of these mechanisms is strongly dependent upon factors such as the shock geometry and the fast magnetosonic Mach number M f = U shock /v f , where U shock is the shock velocity and v f is the fast magnetosonic wave velocity. In lower Mach number shocks, dispersive radiation and wave-particle interactions provide an effective resistivity through the shock ramp and the requisite dissipation routes for the shock's energy conversion (e.g., Kennel et al. 1985;Balogh & Treumann 2013, and references therein).
To ascertain the details of the energy transfer in collisionless shocks, we perform a first principles, continuum kinetic simulation of a perpendicular collisionless shock and utilize the field-particle correlation technique (Klein & Howes 2016;Klein 2017;Klein et al. 2017;Howes et al. 2017;Howes et al. 2018;Li et al. 2019;Chen et al. 2019;Klein et al. 2020;Horvath et al. 2020) to characterize this energy exchange directly in phase space. We consider a reduced dimensionality and simplified geometry to isolate the available energization mechanisms available to the plasma, focusing on the energization mechanisms of shock-drift acceleration, for the ions, and adiabatic heating, for the electrons. We emphasize that, while these processes have been studied previously using kinetic simulations and the particle-in-cell numerical method in higher dimensionality and greater generality (e.g., Park et al. 2013;Guo et al. 2014a,b;Park et al. 2015;Xu et al. 2020), this study is the first direct diagnosis of the energy transfer in a collisionless shock in phase space and identification of the velocity-space signatures of shock-drift acceleration and adiabatic heating.
This Eulerian perspective-focusing on individual regions of phase space for determin-ing the details of the energy exchange, in contrast to the more commonly used Lagrangian perspective of integrating particle trajectories to identify how individual particles are energized-is of high utility for interpreting spacecraft data. For example, using Magnetospheric Multiscale mission measurements of the electron distribution function in the Earth's turbulent magnetosheath and the field-particle correlation technique, Chen et al. (2019) found both the velocity-space signature of electron Landau damping and determined that the observed turbulence was principally dissipating via electron Landau damping. In this regard, the work presented here is the beginning of a broader program of study to identify the velocity-space signatures of energization mechanisms in collisionless shocks and deploy the field-particle correlation technique for the analysis of energy exchange using in situ measurements of collisionless shocks.
To understand the observed velocity-space signatures and connect the resulting signatures to known mechanisms for plasma energization, we construct simplified analytical models for ions and electrons being energized by similar processes absent the complications of a fully self-consistent shock, computing the field-particle correlation on the particle distribution functions predicted by these idealized models. These simplified models allow us to proceed pedagogically and connect the two distinct pictures, the Eulerian point-of-view for identifying where in phase space the particles are being energized and the Lagrangian point-of-view for analyzing how individual particles gain and lose energy. And while significant intuition is gained from the Lagrangian perspective, this novel Eulerian perspective provided by the field-particle correlation technique has some advantages, chief among them the ability to easily distinguish how different regions of phase space are being energized. We will show how the field-particle correlation technique allows us to easily separate the energy exchange occurring between the electromagnetic fields and the multi-component distribution functions (e.g., incoming beam vs. reflected ions in the shock foot and ramp) which frequently characterize collisionless shock dynamics. We may thus distinguish the different effects the same electromagnetic fields are having on different parts of the distribution function, from how the cross-shock electric field decelerates the incoming bulk flow and accelerates reflected ions, to how the motional electric field supporting the upstream E × B motion energizes both the reflected ions and bulk electrons.
The rest of the paper is organized as follows. In Section 2 we provide details of the simulations performed, followed by a broad overview of the results of the simulations examining the shock structure in the electromagnetic fields and electron and ion distribution functions. We then present an overview of the central analysis tool of this paper, the field-particle correlation technique, in Section 3. We apply the field-particle correlation technique to obtain the key results of the paper in Sections 4 and 5: the velocity-space signatures of (i) shock-drift acceleration of the ions and (ii) adiabatic heating of the electrons. We conclude in Section 6 with a discussion of the implications of the results presented for spacecraft observations and future avenues of research applying field-particle correlations to a larger range of shock parameters.

Computational Model and Overview of Results
To perform a self-consistent simulation of a perpendicular collisionless shock, we employ the continuum Vlasov-Maxwell solver in the Gkeyll framework . We emphasize that, unlike traditional particle based approaches such as the particle-in-cell method, Gkeyll directly discretizes the Vlasov-Maxwell system of equations on a phase-space grid to obtain a high fidelity representation of the distribution function, free of the shot noise introduced by finite sized particles. In other words, we solve the following system of equations with a grid-based method for every species s in the plasma, where f s = f s (x, v) is the particle distribution function for species s, q s and m s are the charge and mass of species s respectively, E = E(x) and B = B(x) are the electric and magnetic fields respectively, and the coupling between the electromagnetic fields and particles is given by velocity moments of the particle distribution function, i.e., the charge and current density. This approach has been previously leveraged in the study of electrostatic collisionless shocks (Pusztai et al. 2018;Sundström et al. 2019), allowing for a detailed study of the phase-space dynamics which result from the evolution of the shock. Comparisons to the particle-in-cell method for the study of kinetic instabilities have clearly demonstrated the advantages of using a continuum representation to eliminate discrete particle noise in the particle velocity distributions (Skoutnev et al. 2019;Juno et al. 2020).
Here, a perpendicular shock refers to the orientation of the magnetic field with respect to the shock normal. Since the magnetic field is perpendicular to the shock normal, in one spatial dimension, we require only the two velocity components perpendicular to the magnetic field to fully describe the dynamics of the system, i.e., 1D-2V. The particular 1D geometry we choose is the one spatial coordinate along the shock normal in the x direction, with the initial magnetic field in the z direction, B(t = 0) = B 0ẑ . For completeness, in this dimensionality and field geometry, the Vlasov-Maxwell system of equations is where we have added a collision operator C[f s ] on the right-hand side of of the Vlasov equation †. The electrons and ions are initialized with the same supersonic flow directed in the negative x direction towards a reflecting wall, which leads to a shock wave that propagates in the positive x direction in our simulation. Note that the particles reflect from the wall, † The addition of a collision operator on the right-hand side of the Vlasov equation introduces some semantic ambiguity of the name of this equation, which, with the inclusion of a collision operator, is now formally the Boltzmann equation and the system of equations the Boltzmann-Maxwell system of equations-see Hénon (1982) for a discussion of this linguistic history. Because the collision operator is principally employed for numerical reasons and to provide velocity-space regularization, we will continue to refer to the equation system of interest as the Vlasov-Maxwell system of equations to emphasize our focus on collisionless physics. but the "reflecting wall" boundary condition for the electromagnetic fields is a conducting wall boundary condition in the traditional sense, with zero normal magnetic field and zero tangential electric field. This method of initialization is often called the "reflectingwall" setup ‡ and has been previously employed in numerous particle-in-cell studies of collisionless shocks (e.g., Papadopoulos et al. 1971;Spitkovsky 2005Spitkovsky , 2007. Detailed parameters are as follows: the reflecting wall for the particles and conducting wall for the electromagnetic fields are at x = 0, and plasma is injected with a copy boundary condition † at x = 25d i , where d i = c/ω pi is the ion collisionless skin depth. Here, c is the speed of light, and ω pi = e 2 n 0 /ǫ 0 m i is the ion plasma frequency. Note that the subscript 0 denotes the upstream value, e.g., n 0 is the upstream density and B 0 is the upstream magnetic field magnitude. We use a reduced mass ratio between the ions and electrons, m i /m e = 100. The total plasma beta, β = 2µ 0 n 0 (T e0 + T i0 )/B 2 0 = 2, with the ion beta, β i = 1.3, and electron beta, β e = 0.7. Both the ions and electrons are non-relativistic, with v te /c = 1/16, where v ts = 2T s0 /m s . With this choice of electron beta and v te /c, the ratio of the electron plasma frequency, ω pe = e 2 n 0 /ǫ 0 m e , to the electron cyclotron frequency, Ω ce = −eB 0 /m e , is ω pe /Ω ce ∼ 13.4. The in-flow velocity in the simulation frame to initialize the perpendicular, electromagnetic shock is Note that the in-flow velocity is negative because the plasma initially flows in the negative x direction. Since the plasma is initialized with a flow in a background magnetic field, we initialize the corresponding electric field to support this flow, For the grid resolution in configuration space, we use N x = 1536 grid cells, corresponding to ∆x ∼ d e /6 ∼ 3.7λ D , where d e = c/ω pe and λ D = v te /( √ 2ω pe ) are the electron inertial length and electron Debye length respectively, and we employ piecewise quadratic Serendipity elements for the discontinuous Galerkin basis expansion (Arnold & Awanou 2011). In velocity space, the electron extents are ±8v te , and the ion extents are ±16v ti , with zero-flux boundary conditions at the velocity-space limits and N vx = N vy = 64 grid cells for both species, corresponding to ∆v = v te /4 for the electrons and ∆v = v ti /2 for the ions. The basis expansion in velocity space is also piecewise quadratic Serendipity elements. For further details about the algorithm and the choice of basis expansion, we refer the reader to Juno et al. (2018) and .
We have run the simulation with a small amount of collisions to regularize velocity space. In this case, we choose an electron-electron collision frequency, ν ee = 0.01Ω ci , much less than the ion cyclotron frequency, Ω ci = eB 0 /m i , with the ion-ion collision frequency correspondingly smaller based on the square root of the mass ratio, ν ii = 0.001Ω ci . Note that because the ions are hotter than the electrons, they should formally be even more collisionless than the electrons; however, these collisionalities are larger than typical solar wind collisionalities  and are not chosen to be realistically small, but instead chosen to be just large enough to provide regularization of velocity-space structure given finite velocity resolution. Details on the implementation of the collision operator and its conservation properties can be found in Hakim et al. (2019). ‡ In contrast to the "injection" set-up where an injection boundary condition would be employed on each side of the domain and the plasma blocks would collide and form a shock in the middle of the domain.
† For Gkeyll, this means that the value in the layer of cells beyond the x = 25di edge (the "ghost" or "halo" layer) is exactly equal to the value in the layer of cells at the x = 25di edge, for all the quantities being evolved: the distribution functions for the electrons and ions, and the electromagnetic fields. Because the plasma is initialized with a flow propagating in the negative x direction, this boundary condition leads to a continuous injection of plasma from the right wall with the corresponding electric field and magnetic field to support the E × B flow. Figure 1. The x-electric field (a), y-electric field (b), z-magnetic field (c), ion distribution function integrated in v ′ y (d), and electron distribution function integrated in v ′ y (e) after the perpendicular shock has formed and propagated through the simulation domain, t = 11Ω −1 cp . Note that the distribution functions are plotted in the simulation frame fs(x, v ′ x , v ′ y ) for each species s. We have marked an approximate transition from upstream of the shock to the shocked plasma (dashed-dotted line), and likewise an approximate transition from the shock to the dowstream region (dashed line). To mark the oscillation of the electromagnetic fields and the sloshing of energy between the fields and particles in the downstream region, we have used a solid black line to mark the approximate compression of the magnetic field, along with E = 0. We expect the y-electric field to roughly oscillate about zero in the frame of the simulation, as the "reflecting-wall" set-up is performed in the frame of the downstream plasma, where the E × B velocity is zero.
In Figure 1, we show the electromagnetic fields and particle distribution functions for the electrons and ions in x− v x phase space after the perpendicular shock has formed and propagated through the simulation domain, t end = 11Ω −1 ci . Although the downstream of the shock is fairly oscillatory as the energy injected into the plasma by the shock sloshes back and forth between the electromagnetic fields and particles, we can estimate the compression ratio of this low Mach number shock from the average downstream magnetic field to be roughly r ∼ 2.5 (solid black line in B z plot in Figure 1). With this estimate for the compression ratio, we calculate the shock velocity in the simulation frame to be U shock = U x /(r − 1) = 2v A . Note that in this reflecting wall set-up, the simulation frame is equivalent to the frame in which the downstream plasma is at rest.
Thus, combining the velocity of the incoming flow with the velocity of the shock in the downstream frame, this self-consistently produced perpendicular shock is a M A = 5 shock, where M A is the Alfvén Mach number. Equivalently, using the definitions for the sound speed and magnetosonic speed given by where γ i = γ e = 1 + 2/VDIM = 2 since the simulation domain has two velocity dimensions, we find this shock has fast magnetosonic Mach number, M f ≈ 2.89. With these plasma parameters and this magnetosonic Mach number, we note that this shock is supercritical M f > M fcrit ≃ 2 (Wilson III 2016), similar to the Earth's bow shock, and thus bodes well for the ultimate goal of predicting velocity-space signatures of energization mechanisms in spacecraft observations of heliospheric shocks. In this regard, we focus our attention now on the particle distribution functions and the phase-space structure generated through the shock. While the particle distribution functions in x−v x phase space shown in Figure 1 are illustrative of the dynamics through the shock, showing a reflected population of ions (d) and a clear compression of the electrons (e), we can gain further insight into the dynamics of this shock by looking at the distribution function in v x − v y at fixed points in configuration space through the shock. In Figure 2, we plot the ion ((a)-(f)) and electron ((g)-(i)) distribution functions in velocity space at several points through the shock, from upstream through the ramp to downstream.
As an example of the wealth of data contained in the distribution function, we draw special attention to the ion distribution function in the shock ramp. As we move from upstream into the shock, at the beginning of the ramp at x = 22.5d i , we begin to see a small population of reflected ions, forming a small "crescent" distribution in the lower right quadrant of v x − v y space. Further up the ramp at x = 21.5d i , we observe that the incoming ion beam begins to be deflected by the fields in the shock transition, generating a "boomerang" distribution that smoothly connects the decelerated incoming ion beam with the reflected ion population. It is this reflected population, in agreement with previous studies of super-critical shocks (e.g., Ball & Melrose 2001;Balogh & Treumann 2013), which dominates the ion energization and provides a segue to our key result: diagnosing the velocity-space signatures of particle energization in this perpendicular electromagnetic shock. To obtain these velocity-space signatures, we now describe our tool of choice for our analysis of the high quality distribution function data provided by the continuum kinetic simulation: the field-particle correlation technique. Figure 2. The ion (panels (a)-(f)) and electron (panels (g)-(i)) distribution functions in the simulation frame (downstream frame), fs(x, v ′ x , v ′ y ) for each species s, plotted at various points through the shock at t = 11Ω −1 cp . As we move from upstream, x = 24.5di, through the shock ramp, x = 21.5di, we can identify the reflected ion population as well as a broadening of the electron distribution function.

The Field-Particle Correlation Technique
From combining the Vlasov equation and Maxwell's equations, we can obtain a conservation equation for the total energy of the kinetic plasma (e.g., Klein et al. 2017), The first integral represents the energy of electromagnetic fields in the plasma and the second accounts for the combined microscopic kinetic energy † of all plasma species s. In the absence of particle collisions, the net microscopic kinetic energy of a given plasma † Note that the microscopic kinetic energy of a plasma species s includes contributions from the kinetic energy of bulk plasma flows (associated with the first-moment of the distribution) as well as thermal and non-thermal energy contained in the second moment of the particle velocity distribution fs(v). One cannot extract energy from the thermal component, of course, and the irreversible, entropy increasing conversion of free energy in the non-thermal component species may only be changed through collisionless interactions between the particles of that species and the electromagnetic fields.
To explore the energy transfer between fields and particles, we define the phase-space energy density for a particle species s by w s (x, v, t) ≡ m s v 2 f s (x, v, t)/2 in the nonrelativistic limit. Multiplying the Vlasov equation by m s v 2 /2, we obtain an expression for the rate of change of this phase-space energy density, This equation describes the mechanisms that govern how the energy density in the 3D-3V phase space (x, v) evolves, where each term has a clear physical interpretation. The first term on the right-hand side of (3.2) describes how w s (x, v, t) changes due to particle advection from other spatial regions, giving rise in fluid theory to the energy change through pressure forces and heat fluxes †. Because this term describes the advection of particle kinetic energy as particles move from one spatial position to another, when integrated over the full plasma volume, this term yields zero net change of the total kinetic energy of particle species s, W s = dx dv 1 2 m s v 2 f s . The third term on the right-hand side of (3.2) describes the magnetic forces on the particles. Although this term can move kinetic energy from one location in velocity space to another, when integrated over all velocity space, this term does zero net work on the particles, as expected for the magnetic force.
The second term on the right-hand side of (3.2) describes the work done on the plasma species s by the electric field. When (3.2) is integrated over all velocity space and all physical space to obtain the rate of change of the total kinetic energy W s of a particle species s, the first and third terms have zero net contribution (Klein & Howes 2016;Howes et al. 2017), yielding This expression makes clear that the change in species energy W s is due to work done on that species by the electric field, j s · E. In our exploration of particle energization at collisionless shocks, we choose to focus on the second term in (3.2) to investigate the energization of the particles by the electric field. The form of that term demonstrates that the rate of particle energization can be computed at a single-point in physical space x 0 by measuring the electric field at that position E(x 0 ) and the particle velocity distribution at the same position f s (x 0 , v). This fundamental fact underlies the field-particle correlation (FPC) technique (Klein & Howes 2016;Howes et al. 2017;Klein et al. 2017), where the unnormalized correlation (essentially a time-average) of the product of the electric field E(x 0 ) and a term that depends on the particle velocity distribution f s (x 0 , v) over some correlation interval τ is computed by to thermal energy is dictated by the physics of nonequilibrium thermodynamics in this kinetic system. † In the context of fluid theory, it has been shown that these pressure forces can mediate the conversion of bulk flow kinetic energy to random kinetic energy in the velocity distribution (Yang et al. 2017). This distinction between energy converstion and energization is further discussed in Appendix C.
The resulting correlation C E (x 0 , v, t, τ ) directly measures the rate of change of phasespace energy density at position x 0 as a function of 3V particle velocity space v, producing a velocity-space signature that is characteristic of the mechanism of energization and can be used to identify a particular, locally-occurring energization process, e.g., Landau damping Klein et al. 2017) and cyclotron damping (Klein et al. 2020). We note that as part of this identification, further analysis may be required to ascertain certain details; for example, if one obtains velocity-space signatures corresponding to the presence of Landau damping, the resonant velocity the velocity-space signature is concentrated around is necessary to determine what wave modes are Landau damping in the plasma, as one can find similar structure whether a Langmuir wave  or kinetic Alfvén wave Horvath et al. 2020) is undergoing Landau damping. However, even with this caveat, a key advantage of the FPC method to diagnose particle energization is that it requires only measurements at a single spatial point x 0 to determine the energization by the electric field. An appropriately instrumented single spacecraft mission can provide the requisite full 3V particle velocity distribution f s (x 0 , v, t) and electric field E(x 0 , t) at the spacecraft position x 0 as a function of time. Thus, the velocity-space signatures determined here using kinetic numerical simulations, our key results, may be directly sought using spacecraft observations. In the case of particle energization as a consequence of the dissipation of weakly collisional plasma turbulence, the rate of particle energization represented by the second term in (3.2) generally includes two distinct contributions: (i) an often large-amplitude oscillatory component that leads to zero net energization which is associated with undamped wave motion, and (ii) a typically smaller amplitude secular component that corresponds to the net collisionless transfer of energy from the fields to the particles (Klein & Howes 2016;Howes et al. 2017). By an appropriate choice of the correlation interval τ , the oscillatory energy transfer is largely eliminated by the time-average, exposing the secular energy transfer associated with the collisionless damping of the turbulent fluctuations. For the perpendicular collisionless shock in this study, the shock is quasi-stationary in the shock-rest frame of reference, with smooth electromagnetic fields through the shock as seen in Figure 1, and thus we need not time-average the correlation, but instead take the instantaneous correlation (the limit τ → 0). We will thus suppress the dependence of the correlation on τ henceforth. We note that the FPC with τ = 0 is simply the instantaneous rate of change of the phase-space energy density, ∂w s /∂t, due to work done on the particles by the electric field. If kinetic instabilities were to arise upstream or within the shock transition region, or if the shock itself were to become nonstationary, then it is likely that taking a correlation interval τ longer than either the unstable wave period or the shock reformation time would be necessary to recover a meaningful velocity-space signature of the net particle energization.
In addition, we adopt two final modifications of the FPC analysis that are well suited for the study of collisionless shocks: (i) we separate the contributions to the rate of energization by the different components of the electric field, E x and E y ; and (ii) we replace v 2 in Eq. (3.4) by the component associated with the electric field, e.g., using v 2 x for the correlation using E x . We refer the reader to Appendix A for a discussion of the validity and usefulness of this transformation. Therefore, the form of the FPCs implemented here for a position x = x 0 in our 1D-2V Gkeyll simulation is given by An issue which cannot be overemphasized in performing the FPC analysis of a collisionless shock is making a judicious choice of the frame of reference in which to calculate Eq. (3.5) and Eq. (3.6) (Goodrich & Scudder 1984). We choose to evaluate the correlations in the frame of reference in which the shock is at rest (the shock-rest frame, unprimed variables), as opposed to the frame of reference of the simulation, in which the plasma is at rest downstream of the shock (the simulation frame or downstream frame, primed variables) †. For clarity, the shock velocity in the simulation frame is given by U shock = U shockx = 2v Ax . It is critical not only that the velocity coordinates are transformed to the shock-rest frame, v = v ′ − U shock , but also that the electromagnetic fields are appropriately Galilean transformed to the shock-rest frame, and B = B ′ . We note that our discussion and application of the FPC to the collisionless shock is principally concerned with how the plasma is energized via the electromagnetic fields and focuses on the phase-space dynamics governed by the electric field term in the Vlasov equation. As mentioned previously, plasmas additionally convert bulk kinetic to thermal energy, and vice versa, via other terms in (3.2) such as the v · ∇ term, which gives rise to pressure forces and heat fluxes. To explore how these other physical mechanisms impact the flow of energy through 3D-3V phase space, one can perform complementary correlations with these other terms. Correlating with the magnetic term in the Lorentz force allows the determination of how the magnetic field leads to changes in w s (x, v, t) as a function of velocity v-e.g., energy can be moved between different degrees of freedom by the magnetic field, even though the net energy change (integrated over velocity space) must always be zero. Similarly, if spatial gradients of f s (x, v, t) are available, the velocity-space signatures of the work done on the particles by the pressure tensor can be determined. Of course, computing the total rate of change of the phase-space energy density w s (x, v, t) at a particular point in configuration and velocity space requires all terms of (3.2).
Our focus here, however, is on the term in the Vlasov equation which produces net energization of a plasma species s, the electric field term in (3.2). In fact, as shown in Appendix B, these additional terms such as the v ×B term, can have a cancellation effect on the evolution of the phase-space energy density, so that the net energization due to, for example, an E × B drift is identically zero, as it should be. As such, we are well justified in formulating the FPC to focus only on the net energization and avoid obfuscating the signatures of energization with the additional motion of phase-space energy density due to these other terms in the Vlasov equation. For a further discussion of energization versus energy conversion, we refer the reader to Appendix C.  Figure 3. The ion distribution function fi(vx, vy) (a), and the CE x (b) and CE y (c) components of the FPC, Eq. (3.5) and Eq. (3.6), computed at x = 22.9di from the self-consistent Gkeyll simulation. Note that the FPC is computed in the shock-rest frame. While the bulk incoming ions are slowed down by the cross-shock electric field, Ex, we see the distribution of reflected ions gain energy due the motional electric field, Ey, which supports the incoming supersonic E × B flow.

Velocity-Space Signature of Ion Energization
In Figure 3, we plot the ion distribution function f i (v x , v y ) (a), and the C Ex (b) and C Ey (c) components of the FPC, Eq. (3.5) and Eq. (3.6), at x = 22.9d i . The ion distribution function includes both a component from the incoming beam of ions upstream, as well as the aforementioned "crescent" population of reflected ions. The E x contribution to the , shows that the incoming ion beam is being acted upon strongly by the cross-shock electric field E x , but that E x has little effect on the reflected ion population at this position. On the other hand, the E y contribution to the FPC, C Ey (v x , v y ), at position x = 22.9d i , panels (c), shows that the reflected ions principally interact with this component of the electric field, i.e., the motional electric field which supports the incoming E × B flow.
To understand this visual representation of the rate of ion energization over velocity space, recall that the FPC determines the rate of change of the phase-space energy density of a particular plasma species, w s (x, v, t) = m s |v| 2 f s (x, v, t)/2 due to the electric field. The phase-space energy density of the ions, w i , can only change if the number of ions in that volume of phase space changes. Therefore, nested blue and red crescents in Figure 3(c) indicate that ions are accelerated by E y from the blue region to the red region. Conservation of particle number requires that the number of ions lost from the blue region is the same as the number gained in the red region, but because that red region is at higher velocity v y , the net effect, obtained by integrating C Ey over velocity space (v x , v y ), is an increase in the ion phase-space energy density w i . We also note that the observed C Ey signature is a larger amplitude than the observed C Ex , such that E y dominates the energy exchange at this particular point in space. Furthermore, the FPC method computes the rate of change of energy density, so the rate of energization per ion in the low density population of reflected ions is much higher in amplitude than the loss of energy per ion by the much more dense incoming beam. As a first attempt to understand this signature, consider that the gradient length scale of the collisionless shock in our simulation is L shock ∼ ρ i , where ρ i = v ti /Ω ci is the ion Larmor, or gyro-, radius. Therefore, ions encountering this gradient in the magnetic field will not necessarily have closed orbits and smoothly transition downstream. Depending on an ion's gyrophase when it encounters this magnetic field gradient, the ion's new Larmor orbit may cause the ion to move back upstream, where the magnetic field magnitude is smaller. The increased Larmor radius of this reflected ion in the upstream region then allows the ion to gain energy along the motional electric field supporting the incoming E × B motion. This energization of the reflected ion population via E y is consistent with the well-known energization mechanism shock-drift acceleration (Paschmann et al. 1982;Sckopke et al. 1983;Anagnostopoulos & Kaliabetsos 1994;Anagnostopoulos et al. 1998;Ball & Melrose 2001;Anagnostopoulos et al. 2009;Park et al. 2013). But to understand why shock-drift acceleration would produce the particular velocity-space signature observed in panel (c) of Figure 3, we turn to a simplified analytic model to connect the well-known Lagrangian picture for shock-drift acceleration with the new Eulerian perspective granted by the FPC.

Shock-Drift Acceleration in an Idealized Perpendicular Shock
We consider now a simplified reduction of the electromagnetic fields observed in our self-consistent simulation to a step function in the magnetic field, with amplitude jump B d /B u = 4. We will also continue to work exclusively in the shock-rest frame, where to a good approximation the motional electric field, E y , is a constant through the entire shock. The value of the constant E y , as well as the ion and electron plasma betas, are chosen so that the shock velocity is similar to the selfconsistent simulation, M A = 4.9 and M f = 3.0. This reduced model corresponds to the limit L shock /ρ i ≪ 1 and allows us to decompose the ion motion more easily between upstream and downstream gyro-and E × B motion. To mimic the geometry of the selfconsistent simulation, we take E y < 0 and B z > 0 so that the inflow E × B is in the negative x direction.
In Figure 4, we plot (a) the trajectory of an ion in the (x, y) plane and (b) its corresponding trajectory in (v x , v y ) velocity space in the shock-rest frame, where the colors indicate the corresponding segments of the trajectory. In the upstream region at x > 0 (black), the black circle centered about the upstream E × B velocity (black star) corre- sponds to the Larmor orbit of the ion about the upstream inflow velocity in the (v x , v y ) plane. Upon first crossing the magnetic discontinuity to x < 0, the ion changes to a Larmor gyration in the (v x , v y ) plane (blue) about the downstream E × B velocity (green star). In the larger amplitude downstream perpendicular magnetic field, the radius of the Larmor motion in the (x, y) plane in the shock-rest frame is reduced (blue). Depending on the ion's gyrophase when the ion crosses the magnetic discontinuity, the ion passes back upstream to x > 0, and once again undergoes a Larmor orbit in the (v x , v y ) plane (red) about the upstream E×B velocity (black star). In this segment of the trajectory (red), the ion gains perpendicular energy in the shock-rest frame, graphically represented by the distance in velocity space of the ion from the origin of the (v x , v y ) plane. Finally, the ion will eventually cross back into the downstream region to x < 0 (green), resuming its Larmor orbit in the (v x , v y ) plane (green) about the downstream E× B velocity (green star). Without any additional crossings of the magnetic discontinuity, the ion will simply E × B drift downstream.
In the segment of the trajectory where the ion can gain energy, it is the motional electric field, E y , that is doing positive work on the ion, exactly like in our self-consistent simulation. We note that this ion's dynamics-the reflection due to the magnetic gradient and energy gain from its traversal upstream and alignment with the motional electric field-is the well-known single-particle picture of shock-drift acceleration. In fact, this picture in velocity space of where a single ion gains energy via this reflection by a magnetic gradient has been previously noted (Gedalin 1996a). We wish now to connect this Lagrangian perspective on how a single ion gains energy from this reflection off a magnetic gradient to the Eulerian point-of-view we have from the FPC.

Velocity-Space Signature of Shock Drift Acceleration
To connect the single-particle picture of shock-drift acceleration with how a distribution of ions is energized, we employ the Vlasov-mapping technique (Scudder et al. 1986;Kletzing 1994;Hull et al. 1998;Hull & Scudder 2000;Hull et al. 2001;Mitchell & Schwartz 2013, described in Appendix D, to determine the velocity distribution function in our simplified model for the electromagnetic fields through the shock. We show in Fig (c) (d) Figure 5. Comparison of the reconstructed ion distribution function (a) and CE y component of the FPC (b) computed from this reconstruction to the self-consistently produced ion distribution function (c) and CE y component of the FPC (d) from the Gkeyll simulation. Using the Vlasov-mapping technique we can connect the single-particle orbits (over-plotted white (a) and black (b) lines) to the distribution function dynamics. CE y integrated over velocity space is net positive, meaning the observed velocity-space signature corresponds to an energization process. We identify this particular velocity-space signature as the signature of shock-drift acceleration, energization of the reflected ions via the motional electric field in the upstream, via the connection between where in velocity space a single ion is energized and the specific region of velocity space where the strongest energy exchange is occurring.
computed from the motional electric field, E y , and gradients of this reconstructed distribution function. In addition, we repeat Figure 3, panels (a) and (c), for reference in comparing the distribution function and generated velocity-space signature between the simplified model and self-consistent simulation.
In the reconstructed distribution function from the idealized model, we identify, in addition to the incoming upstream population centered at the upstream E × B velocity, a component of reflected particles that have returned upstream, exactly like in the self-consistent simulation. Overplotted on the ion distribution function and computed FPC from the Vlasov-mapping technique is the trajectory in (v x , v y ) for the ion analyzed in Figure 4, showing that this reflected population and velocity-space signature are coincident with the red segment of the trajectory in Figure 4. Integrating this field-particle correlation over velocity space simply yields the net rate of work done by E y , C Ey (v x , v y )dv x dv y = j y E y , and we find the integration to be positive. We thus identify the whole population of reflected ions as experiencing net energization, with the velocityspace signature of this energization process, shock-drift acceleration, given by Figure 5 (b) and (d).
We have now connected the Lagrangian picture of shock-drift acceleration with the Eulerian picture provided by the FPC technique, and we conclude this section noting that while shock-drift acceleration has been studied extensive theoretically and numerically (e.g., Gedalin 1996aGedalin ,b, 1997Gedalin et al. 2000;Park et al. 2013;Guo et al. 2014a,b;Park et al. 2015;Gedalin et al. 2018;Xu et al. 2020), the velocity-space signature of shock-drift acceleration provides a new perspective on the energization of the ions in phase space via this process. In both cases, we understand that a portion of the distribution of ions are reflected via the magnetic field gradient and return upstream, where they can gain energy via the motional electric field. Although the single-particle trajectory in phase space guides our understanding of where we expect the ions to be gaining energy, using the FPC technique enables us to see clearly the exact region of phase space in which ions are being energized via shock-drift acceleration.
This confirmation of the velocity-space signature of shock-drift acceleration in a selfconsistent simulation is a vitally important step for comparison to measured velocityspace signatures of energy exchange using in situ spacecraft measurements; however, it is also interesting that the velocity-space signature of shock-drift acceleration is unchanged between the idealized model and a self-consistent simulation given the additional physics of the self-consistent simulation: the finite shock width and cross-shock electric field. We explore the reasons for the excellent agreement despite these two key differences between the simulation and the idealized model in Appendix E, where we find the cross-shock electric field assists in reflecting ions, allowing the ions to traverse further back upstream and gain additional energy via shock-drift acceleration. Thus, while the combination of the finite shock width and cross-shock electric field quantitatively changes the population of ions that are reflected, the qualitative signature of energization in velocity space via shock-drift acceleration remains unchanged.

Velocity-Space Signature of Electron Energization
We now examine the energization of the electrons by the simulated perpendicular collisionless shock. Similar to Figure 3 for the ions, we show in Figure 6(a) the electron distribution function f e (v x , v y ) and (b) the C Ex and (c) the C Ey components of the FPC, Eq. (3.5) and Eq. (3.6). As shown in Figure 1(e), the electron distribution broadens through the entire shock ramp, so we plot in Figure 6 the results of the FPC analysis at x B = 21.8d i , where the cross-shock electric field peaks. In panel (a), the center of the distribution in the shock-rest frame is displaced away from the origin to v x < 0 and v y < 0 due to the particle drifts in the varying electric and magnetic fields through the shock ramp. Because the thermal width of the electron velocity distribution is much larger than the net drift of the distribution in the shock-rest frame, computing the C Ex and C Ey correlations using Eq. (3.5) and Eq. (3.6) leads to the qualitative "two-lobed" velocity-space signatures observed in Figure 6(b) and (c). The small drifts-i.e., |v x /v te | ≪ 1-lead to a slight asymmetry of the two-lobed structure, so we over-plot contours of constant C Ex,y to make these slight asymmetries more visually apparent. Although the gain (red) (c) Figure 6. The electron distribution function fe(vx, vy) (a), and the CE x (b) and CE y (c) components of the FPC, Eq. (3.5) and Eq. (3.6), computed at xB = 21.8di from the self-consistent Gkeyll simulation. Note that the FPC is computed in the shock-rest frame. The contours on the FPC plots, which are the same in both CE x and CE y , make clear that Ey leads to a net loss of electron energy, whereas Ex yields a net increase of electron energy.
and loss (blue) of electron energy largely cancels out upon integration over velocity space (v x , v y ), asymmetries in the two lobes lead to a non-zero net energization. In Figure 6(b), the asymmetry of the velocity-space signature leads to a net energization of the electrons by the cross-shock component of the electric field E x . And, in contrast to the shock-drift acceleration of the ions by the motional electric field E y seen in Figure 3(c), we see in Figure 6(c) that the electrons experience a net loss of energy due to the E y component of the electric field. As with the ion analysis, we now turn to an idealized model for the electron dynamics through the shock layer to understand these velocity-space signatures for the electron energization.

Adiabatic Heating in an Idealized Perpendicular Shock
Unlike the ions for which L shock ρ i , the electron gyroradius is much smaller than the gradient length scale of the shock, L shock ≫ ρ e , and we thus expect the electrons to stay well magnetized through the shock. Therefore, we adopt a simplified model for the electron dynamics through the shock by taking a shock ramp with a linearly increasing magnetic field over a length L = 2d i and mass ratio m i /m e = 1836, satisfying the limit L shock ≫ ρ e . The other parameters are the same as the simple shock model used for the ion analysis in Section 4.2: a magnetic field increase of B d /B u = 4, a constant and uniform motional electric field E y < 0 in the shock-rest frame, with the same ion and electron plasma betas such that the shock velocity is comparable to the self-consistent simulation M A = 4.9 and M f = 3.0. Note that this idealized model for the electrons has no cross-shock component of the electric field, E x = 0; the implications of this choice are discussed at length in the upcoming subsections. In this model, the increase in the magnetic field magnitude through the ramp leads to a steady decrease in the E × B velocity as the plasma flows through the shock transition. In addition, the gradient of the magnetic field magnitude in the shock-normal direction induces a ∇B drift in the +y direction. In Figure 7, we plot (a) the profile of the perpendicular magnetic field B z (x) (blue) and the motional electric field E y (x) (red) along the shock normal direction, as well as (b) the trajectory of an electron in the (x, y) plane as it flows through the shock ramp over 0 x/d i 2. The trajectory plot shows clearly the ∇B drift in the +y direction. We note that, as the electrons flow through the shock ramp over 0 x/d i 2 and undergo a ∇B drift in the +y direction, (c) the net energization for a distribution of electrons j y E y is positive.
In the region where the perpendicular magnetic field changes magnitude, 0 x/d i 2, the ∇B drift is anti-aligned with the motional electric field, thus leading to net energization of electrons. In fact, the rate of energization of the electrons by the ∇B drift in the motional electric field is precisely that needed to conserve the first adiabatic invariant of the electron, i.e., the magnetic moment µ = m e v 2 ⊥ /2B z . This energization via conservation of the electron's adiabatic invariant is thus often referred to as adiabatic heating.
We can show the relationship between the energization via the ∇B drift and the conservation of the electron's magnetic moment by considering the change in the perpendicular kinetic energy of the electrons, where the magnitude of the ∇B drift in the +y direction is given by For the static fields in this model, the total time derivative is dominated by the E × B velocity, proving that the electron's magnetic moment µ is conserved. As before, we now wish to connect this single-particle, Lagrangian picture of adiabatic heating with the Eulerian picture provided by the FPC technique.

Cross-Shock Electric Field Impact on Velocity-Space Signature of Adiabatic Heating
We again employ the Vlasov-mapping technique, as in Section 4.3, to reconstruct the electron distribution function through the idealized shock model and the resulting FPC velocity-space signatures. In Figure 8, we present (a) the resulting electron distribution function f e (v x , v y ) and (b) the C Ey correlation at position x A = 1.8d i in the model, where the velocity-integrated energy transfer rate j y E y is positive, as seen in Figure 7(c). Note that the small shift (relative to v te ) of the distribution to v y > 0 due to the ∇B drift leads to an asymmetry in the velocity-space signature because of the v 2 y weighting in Eq. (3.6) for C Ey †. Although a large part of the energy transfer rate represented by this two-lobed velocity-space signature cancels out upon integration over v y , the slight asymmetry leads to a net positive energization of the electrons, yielding j y E y > 0, as plotted in Figure 7 (c).
We compare the predicted velocity-space signature of electron adiabatic heating from the idealized shock model in Figure 8(b) to the Gkeyll simulation results, where we plot in Figure 8(c) the f e (v x , v y ) and (d) the C Ey from the simulation at position x B = 21.8d i . While the correlation C Ey from the simulation has the same qualitative, two-lobed structure indicative of drift energization of the electrons, the net drift of the electron velocity distribution has v y < 0, and so therefore the resulting asymmetry in C Ey has the opposite overall sign, leading to a net loss of energy for the electrons due to E y . How can we reconcile these apparently contradictory results for the electron energization by the motional electric field E y , particularly given that we have shown in Eq. (5.3) that the ∇B drift in the +y direction leads to adiabatic heating of the electrons?  Figure 5, the model CE y displays the opposite asymmetry from the CE y computed from the self-consistent simulation. Even though the signatures are qualitatively similar, this first computation of CE y suggests that Ey is responsible for a net loss of energy through the shock.
We can check if the adiabatic invariant of a distribution of electrons, is constant through the shock, and whether our intuition about how the electron temperature should increase through a magnetic field gradient has merit. Note that in the 1D-2V geometry of the simulated perpendicular shock where both velocity coordinates are perpendicular to the magnetic field, the adiabatic invariant of the distribution of electrons can be simplified to as shown in Figure 9, the adiabatic invariant for the distribution of electrons, µ e , is well conserved through the shock transition. This result suggests that the cross-shock electric field is complicating the electron dynamics and energy exchange through the shock. Here we can leverage the Vlasovmapping model to explore the effect of the cross-shock electric field on the energetics of the electrons. In Figure 10, we present a comparison of (b) the electron trajectories, (c) the rate of energization of the electrons by the electric field j e · E, and (d) the cumulative electron energization integrated from upstream x xup dx j e · E for two models: (i) a "full model" (solid) which integrates electron trajectories in the self-consistently produced electromagnetic fields in Figure 10(a); and (ii) a "zero E x model" (dashed), in which we artificially set the cross-shock electric field to zero.
The trajectories in Figure 10(b) show a clear qualitative difference between the zero E x model and the full model: in the zero E x model (blue), the transverse drift through the shock ramp is relatively weak and in the +y direction, whereas the full model (black) yields a much larger amplitude drift in the opposite direction. This qualitative difference in the transverse drift direction explains the stark differences in C Ey between the idealized model in Figure 8 (b) and the simulation in (d).
Looking at the rate of electron energization by the electric field in Figure 10(c), we indeed see that j ye E y (red dashed) is positive for the zero E x model (and is the only means of energization of the electrons since E x = 0), but j ye E y is negative for the full model (red solid). However, when the energization by j xe E x (blue solid) is combined with j ye E y (red solid), the net rate of energization j e · E of the two models is exactly the same (black solid and red dashed overlap). Therefore, although the dynamics of the electrons differ qualitatively in the presence or absence of the cross-shock electric field, their energization is the same in either case. To explain this puzzling finding, we exploit the limit L shock ≫ ρ e to execute a guiding-center drift analysis of the electron energization.

Guiding-Center Drift Analysis of Electron Energization
In the idealized model presented in Section 5.2, there are only two drifts: an E × B drift in the x direction due a constant E y through the magnetic field ramp and a ∇B drift in the y direction due to the linearly increasing magnetic field. The introduction of a cross-shock electric field, along with the transition from a single-particle picture to a distribution of particles, adds several new drifts to the full list of potentially dynamically important drifts. We now not only have an E × B drift that has a component in the x direction due to the motional electric field E y supporting the incoming supersonic flow, but also a new y component due to the cross-shock electric field E x . For a distribution of electrons the ∇B drift in the y direction is modified from its single particle form, where p ⊥,e is the electron perpendicular pressure, We now must also consider a polarization drift in the x direction, and we note that the total time derivative as the electrons flow through the shock is dominated by the convective contribution d/dt = ∂/∂t + U · ∇ ≃ U x ∂/∂x, giving Finally, we have the magnetization drift, ∇ × M, a bulk drift in the y direction due to the increasing density through the shock ramp, (5.11) Although there can, in principle, be a polarization drift in the y direction due to the variation in E y through the shock, in the shock rest frame E y changes very little, so this drift is negligible. We note that the ∇B drift and magnetization drift can be combined to form the diamagnetic drift, This generalization from the single-particle picture to a distribution of particles is somewhat subtle and often dubbed Spitzer's paradox from the early work on plasma equilibria in a fusion context (Spitzer 1952(Spitzer , 1962Qin et al. 2000). While the concept of pressure is ill-defined for a single particle, magnetic field gradients must inevitably be balanced by pressure gradients in equilibrium because the ∇B drift depends on the particles' velocity, and thus different parts of the distribution of particles will experience different ∇B drifts. As a consequence, the pressure of the plasma must change through the magnetic field gradient, presuming of course that the plasma is magnetized. To understand this generalization from the single-particle picture to a distribution of particles,  (c) (d) Figure 11. (a) A comparison of the major drifts through the shock evaluated in the shock-rest frame of reference: (i) E×B drift in x (black) and y (green dashed), (ii) the ∇B drift in y (blue), (iii) the magnetization drift in y (red dashed), and (iv) the polarization drift in x (magenta dashed-dotted). We check that these drifts sum to the total first moment computed from the electron distribution function (b) as well as determine how each of these drifts contributes to the overall energy exchange, je · E (c) and compare the je · E computed from these drifts to the total je · E computed from moments of the electron distribution function. Note that in comparing how each of these drifts contributes to je · E, we plot the polarization drift multiplied by Ex (green dashed-dotted), the ∇B and magnetization drifts multiplied by Ey in the shock-rest frame (blue and red dashed respectively), and the total je · E arising from both components of the E × B flow (black), as we expect the total energization due to the E × B flow to be zero.
a detailed derivation of these drifts from the first moment of the Vlasov equation can be found in Appendix F. We plot these drifts in the shock-rest frame in Figure 11(a) as a function of the distance through the shock, showing that there is a clear ordering of the magnitude of these drifts: both components of the E × B drift are dominant, the ∇B and magnetization ∇ × M drifts provide equal contributions at a smaller amplitude, and the polarization drift is yet smaller. To demonstrate that the sum of these drifts indeed completely describes the bulk flow of the plasma in the shock-rest frame, we show in (b) that the first moment of the electron distribution agrees well with the sum of the drifts for both x and y components.
The rate of energization of the electrons in the drift approximation is equal to the rate of work done by the electric field on the drifting electrons, given by j d,e · E = q e n e U d · E (5.13) where U d is the total drift motion, which can be decomposed into the contributions by the individual drifts identified above in, e.g., Eqns.

2019
). However, we emphasize again that we have generalized from a single-particle picture to a distribution function picture as part of our goal to understand the energization of the electrons in phase space, and thus our equivalent guiding center energization equation has additional terms such as the energization of the plasma via the magnetization drift Eq. (5.10), which do not appear in the guiding center energy equation of a single particle, e.g., Eq. (1) in Dahlin et al. (2014) and Dahlin et al. (2016) †. We reiterate that the perspective provided by the Lagrangian point-of-view, wherein one considers the energization of individual particles, has significant merit, but that the Eulerian perspective has its own advantages, and to obtain the Eulerian point-of-view we must generalize from single particles to distributions of particles.
In Figure 11(c), we plot the rate of electron energization by each of these four drifts in the shock-rest frame. The first key takeaway from this figure is that, although the two components of the E × B drift dominate the total drift motion, there is no net work done by the E × B drift as we expect because (E × B) · E = 0. So, while the energization due to one component of the E × B drift may appear large, summed over all components the energization must be zero, as shown in Figure 11(c). In this regard, separately plotting the contributions to the rate of electron energization j e · E by the different components of the E × B drift can be somewhat misleading, as both components of the E × B flow are much larger than the other drifts.
That the full E × B drift leads to zero net energization of the particles also explains the puzzling finding in our single-particle modeling of the electron energization shown in Figure 10(c). Although the j ex E x and j ey E y from the full model were larger than and significantly different from the zero E x model, when summed they yielded the same net rate of energization of the electrons as the zero E x model. This cancellation is exactly the result of the two components of "energization" from the E × B flow cancelling, as we know they must. The remaining net energization is then solely from the other drifts and their alignment with the motional electric field E y .
In Figure 11(d), we check that the total rate of energization of the electrons by the electric field in the shock-rest frame, j e · E, agrees with the sum of the energization by the ∇B, magnetization, and polarization drifts, finding good agreement. Note that for this comparison, we are computing the electron current in the shock-rest frame from the electron distribution function, i.e., (5.14) After the comparison between the full model and zero-E x model in Figure 10(c) and (d) revealed that the total j e · E was roughly equivalent between the two models, when the zero-E x model only had energy exchange due to the alignment of the ∇B drift with the motional electric field, we might have anticipated that the only energy gain was due to this same adiabatic heating process from the idealized model in Section 5.2. Importantly though, we see from the drift analysis in Figure 11(c) that the model in Section 5.2 must be generalized to the case when a distribution of particles is drifting. While the energy gain by a single electron ∇B drifting in the model fields is exactly † In fact, this generalization from a single-particle picture to a distribution function picture is also pointed out, at least implicitly, by Dahlin et al. (2014) (and subsequent studies) when summing over particles, thus permitting the definition of the pressure of the plasma, and measuring the total heating of the plasma. While the curvature drift term of interest in these studies is proportional to v 2 for a single particle, the evolution of the energy density, or integrated energy density as in Eq. (5) in Dahlin et al. (2014), transforms this term to be proportional to p because particles of different parallel velocities experience different curvature drifts. the energy gain required for that single electron's magnetic moment µ to be conserved, in the self-consistent simulation it is not only the ∇B drift that ensures the electron distribution's adiabatic invariant µ e is well conserved in the shock. We also have an equal contribution to the energization from the magnetization drift. Together, as shown by Eq. (5.12), the ∇B and magnetization drift are equivalent to the diamagnetic drift, u diamag = u ∇B + u ∇×M , so another perspective on the electron energization via adiabatic heating is the adiabatic invariant of a distribution of electrons, µ e in Eq. (5.4), is conserved due to the alignment of the diamagnetic drift and the motional electric field, E y .
Two important questions remain: (i) what is responsible for the slight disagreement between the energy gain due to the ∇B and magnetization drifts and j e · E from velocity moments of the electron distribution function in Figure 11(d)?; and (ii) even if the physics of adiabatic heating is the same with the electrons gaining energy via the alignment of drifts with the motional electric field, is the velocity-space signature of adiabatic heating the same as that predicted by the Vlasov-mapping model in Figure 10(b)? To answer the first question, we have performed a more realistic mass ratio simulation, m i /m e = 400, in Appendix G where we find better agreement between the energy gain due to the ∇B and magnetization drifts and j e ·E from velocity moments of the electron distribution function. Thus, the small disagreement between these methods of measuring the electron's energy gain in the m i /m e = 100 simulation is simply due to the fact that the electron gyroradius is not asymptotically smaller than the shock's extent. To answer the second question, we seek a means of eliminating the E × B component of the energy exchange in the FPC.

Velocity-Space Signature of Adiabatic Electron Heating
If the large E × B flows are polluting the analysis of the overall exchange of the energy, when fundamentally the electron heating is principally due to the alignment of the ∇B and magnetization drifts with the motional electric field, we return to the velocity-space signature plotted in Figure 6(c) and Figure 8(d) to determine how we might remove the contribution from the large E × B flows to the FPC signal. Guided by the fact that the transverse electric field component E y governs the adiabatic heating of the electrons, we seek to eliminate the large contribution to the rate of energization associated with the y-component of the E × B drift (which is ultimately canceled by energization associated with its x-component, as we show in Appendix B). Therefore, we transform to a frame of reference moving in the transverse direction at the same velocity as the y-component of the E × B drift, U td = −E x /B zŷ . We define the transverse drift frame of reference: (i) in the x direction, or shock-normal direction, the shock is at rest; (ii) in the y direction, or transverse direction, the frame moves at a velocity equal to the y-component of the local E × B drift in the shock-rest frame.
Critically, the electric field transforms from the shock-rest frame (sf , unprimed) to the transverse drift frame (td, double primed), as where we have assumed E z = 0 in this 1D-2V perpendicular shock. Therefore, in the transverse drift frame, the cross-shock electric field is zero, E ′′ x = 0, and the motional electric field is unchanged from the shock-rest frame, E ′′ y = E y . The y-component of the velocity coordinate in the transverse drift frame is to our previous computation of the FPC from Ey using only a frame transformation in vx to the shock-rest frame, CE y (vx, vy) (b). Note that panel (b) is a repeat of panel (c) of Figure 6 and panel (d) of Figure 8. While the previous correlation, CE y , suggested the electrons were losing energy in this degree of freedom, the newly transformed correlation, C E ′′ y , demonstrates that the electrons are in fact gaining energy once the E × B motion in this degree of freedom is subtracted. This energization mechanism, caused by an alignment of the ∇B and magnetization drifts with the motional electric field, Ey, is the same mechanism responsible for energizing the electrons in the idealized model, and the velocity-space signature for this adiabatic heating now exactly matches the results of the idealized model presented in Sections 5.2 and 5.3.
For completeness, the x velocity coordinate and magnetic field are unchanged from the shock-rest frame, v ′′ x = v x and B ′′ = B. Although the transverse drift frame changes with position through the shock as the cross-shock electric field changes along the normal direction, one can determine this frame of reference at any position from a local, single-point measurement of the electric field. In contrast, other drifts generally depend on gradients, and therefore cannot be uniquely specified using only single-point measurements. The benefit of the transverse drift frame of reference is not only that the rate of energization associated with the total E × B drift is equal to zero, which is true in any frame of reference, but also that the rates of energization associated with each component of the E × B drift are separately zero. † Since E ′′ x = 0 in the transverse drift frame, we need only examine the y correlation, to explore the energization of the electrons. We plot in Figure 12(a) the correlation C E ′′ y (v ′′ x , v ′′ y ) in the transverse drift frame at position x B = 21.8d i compared to (b) the corresponding correlation C Ey (v x , v y ) in the shock-rest frame (a repeat of Figure 6(c)).
We note two things immediately from Figure 12. First, transformation to the transverse drift frame has switched the sign of the net rate of electron energization (from integrating the correlation over velocity space) relative to the case in the shock-rest frame, made clear † Note that the transverse drift frame is not the only frame of reference in which the contributions to the energization due to each component of the E × B drift are zero. One could also define a normal drift frame moving at the normal component of the local E × B drift, which would yield E ′′ y = 0. This could be useful in determining the energization associated with the polarization drift in the x direction, but since this is a subdominant contribution, we do not pursue that line of investigation here. by counting the over-plotted equally space contours in the blue-red, two-lobe structure of each case. Second, the correlation C E ′′ y (v ′′ x , v ′′ y ) in the transverse drift frame is strikingly similar to the correlation found in Figure 8(b) computed from the reconstructed distribution function from the idealized model. Both points are perhaps unsurprising, as this transformation has removed the energy exchange in the v y degree of freedom due to the E × B flow, and we are left with a similar energization mechanism found in the idealized model: alignment of the ∇B and magnetization drifts with the motional electric field, E y . The alignment of these two combined drifts (constituting the diamagnetic drift) with the motional electric field is exactly what is required for the electron distribution's adiabatic invariant µ e to be well-conserved and for the electrons to gain energy through the increasing magnetic field of the perpendicular shock. Figure 12(a) contains the second key result of this paper, the velocity-space signature of adiabatic electron heating.
It is worth emphasizing that transformation to the transverse drift frame not only revealed the same velocity-space signature for adiabatic heating as in the idealized model analyzed in Sections 5.2 and 5.3, but also allowed the extraction of an energization signature which was buried in the large background energy exchange from the E × B flow. Ultimately, the adiabatic heating via alignment of the ∇B and magnetization drifts with the motional electric field was masked by the large oscillation of energy between the v x and v y degrees of freedom due to the two components of the E×B motion. The need to carefully consider the frame of reference in which the energization analysis is performed, especially due to the impact the frame of reference choice has on the cross-shock electric field, has been pointed out in previous studies (Goodrich & Scudder 1984).
One must therefore exercise extreme caution in attributing energization to a particular drift when separately considering the work done by the different components of the electric field. In Figure 10(c), for example, it would be easy to attribute erroneously the energization of the electrons to the cross-shock electric field E x . Rather, the energization due to the x-component of the E × B drift and E x is exactly canceled by a loss of energy due to the y-component of the E × B drift and E y . The net result is a much smaller positive rate of energization due to the motional electric field E y and the remaining drifts in the transverse direction.
That the addition of the magnetization drift did not qualitatively alter the velocityspace signature of adiabatic electron heating can be understood by considering the general qualitative features of drift energization. In the drift approximation for a warm electron distribution relevant to most heliospheric plasmas, it is often true that the magnitude of the drift (not including the dominant E × B drift, which cannot energize particles) is much less than the thermal velocity of the electrons, U d,e ≪ v te . In this case, the center of the drifting velocity distribution is offset from the origin of velocity space, but this offset will be much smaller than the thermal width of the distribution. When the correlation is computed by taking the velocity derivatives of f e and multiplying by the appropriate component of the electric field weighted by |v| 2 , the result will generally produce a twolobed velocity-space signature, qualitatively similar to that shown in Figure 12(a).
A careful consideration of the drifts is essential to interpret properly the mechanism responsible for the adiabatic heating. But much of the interest in collisionless shocks is focused not on adiabatic heating and instead on identifying mechanisms of non-adiabatic heating. In fact, because the observed temperature increase in the electrons is entirely in the perpendicular temperature to conserve the electron distribution's adiabatic invariant, even if the electron response is initially adiabatic, this anisotropy will itself be a source of instabilities, such as the whistler anisotropy instability (Gary 1993;Gary et al. 2011;Wilson III et al. 2013b;Wilson et al. 2020), that will further complicate the energy exchange.
Likewise, the transverse drift frame relies on the electrons being magnetized so that their E × B motion is well-defined. While the transformation is sensible here because L shock ≫ ρ e , there are many observations of collisionless shocks where the shock ramp is not asymptotically larger than the electron gyroradius (Hobara et al. 2010) and thus would warrant caution in the application of the transverse drift frame to reveal any masked energization signatures such as the velocity-space signature of adiabatic heating in Figure 12(a). In addition, the transverse-drift frame is a non-inertial frame because the E × B velocity is changing through the shock. While the transformation to the instantaneous transverse-drift frame at a single point in configuration space in a simulation is a perfectly reasonable frame transformation, care will be required applying this technique to analysis of spacecraft data, which inevitably average over a small volume of configuration space.
Nevertheless, the velocity-space signatures of the mechanisms that govern non-adiabatic heating are likely to be entirely distinct from the simple two-lobed appearance of adiabatic heating, and so will be easy to distinguish using a field-particle correlation analysis. In addition, it is useful to first characterize the "background" signature of adiabatic heating, as represented by the typical velocity-space signature shown in Figure 12(a). In cases where the adiabatic response produces a temperature anisotropy that drives instabilities, the energetic response, and thus the electron's velocity-space signature through the shock, is likely to be characterized by a combination of both adiabatic and non-adiabatic signatures.

Summary and Future Outlook
This paper presents the first attempt to characterize the energy exchange in a collisionless shock using the field-particle correlation technique (Klein & Howes 2016;Howes et al. 2017). We have examined a self-consistent perpendicular collisionless shock using the continuum Vlasov-Maxwell solver in the Gkeyll simulation framework and our results can be summarized as follows: (i) We have identified the velocity-space signatures of shock-drift acceleration of ions in Figure 5(d) and adiabatic heating of electrons in Figure 12(a).
(ii) Using simplified models of single-particle motion through idealized models of the electromagnetic fields through the shock transition, we identified the conditions under which we expect to observe these velocity-space signatures for these energization processes.
(a) We determined that the velocity-space signature of shock-drift acceleration can be seen clearly in a reflected ion population and is robust to the presence of the finite shock width and cross-shock electric field that arise in the self-consistent simulation.
(b) For the electrons, we determined it was critical to eliminate the energization associated with the separate components of the E × B drift by transforming to the transverse drift frame of reference, as the large E × B drifts in both x and y due to the incoming supersonic flow and cross-shock electric field significantly obscured the effect of the energization of the electrons via other drifts. (iii) Finally, we observed a general strength of this method of analysis using the fieldparticle correlation technique, going beyond identification of where j · E is positive and which components are positive, and furthermore captured the subtleties of a generalization from a single-particle picture to a distribution of particles.
(a) In the case of the ions, a cursory glance at the x component of j · E, which is overall negative due to the slowing down of the bulk distribution, obscures the role the cross-shock electric field plays in increasing the reflected fraction of ions-see Appendix E for further details.
(b) In the case of electrons, just using the components of j · E would completely miss the actual source of energization, as one might expect that the positive x component of j · E corresponds to net energization via the cross-shock electric field. In fact, it is not the cross-shock component of the electric field E x which leads to the observed energization of the electrons through the shock, but rather the motional electric field E y . And while in the single-particle picture the energization is simply the alignment of said motional electric field with the ∇B drift, the generalization to a distribution of electrons introduces an additional drift, the magnetization drift, which when combined with the ∇B drift, forms the diamagnetic drift. When aligned with the motional electric field, these two drifts provide the necessary energization for the adiabatic invariant of the electron distribution to be well conserved through the shock. The work presented here is only the beginning of a program of study to determine, in general, how we may be able to leverage the full information contained in the particle velocity distribution function to ascertain the details of the energy exchange in a collisionless shock. While historically the Lagrangian perspective of examining how individual particles gain and lose energy has led to enormous improvements in our understanding of the dynamics and energetics of collisionless shocks, this complementary Eulerian approach, directly analyzing the energy exchange in phase space using the field-particle correlation technique, has significant value for interpreting both simulation and spacecraft data. Especially when advances in spacecraft instrumentation provide ever higher resolution and higher cadence particle velocity distribution measurements of collisionless shocks (Chen et al. 2018;Goodrich et al. 2018), the time is ripe for fully exploiting the information contained in phase space to provide a deeper of understanding of the mechanisms of particle energization at a collisionless shock.
Further studies of higher dimensional, higher magnetosonic Mach number, and more general geometry collisionless shocks are of the utmost importance. As reviewed in the introduction, there is a large variety of processes not considered in this study that have been studied previously as potential energization mechanisms, from shock surfing acceleration to diffusive shock acceleration. As with the body of work utilizing the field-particle correlation technique for analyzing dissipation via resonant processes, we will require a systematic study of all of these processes if we have any hopes of distinguishing their velocity-space signatures. We may expect certain energization processes such as diffusive shock acceleration and the Fast Fermi process, which also rely on particle reflection, may have qualitatively similar velocity-space signatures to shock-drift acceleration but still contain all the requisite information to identify the particular process locally occurring, just as we can utilize information such as the velocity around which a resonant wave process is identified to characterize the particular waves which are resonantly energizing the plasma.
We expect as we move to higher dimensionality, higher magnetosonic Mach number, and more general shock geometry, our analysis will also be further complicated by upstream kinetic instabilities such as those observed in the Earth's bow shock. The velocityspace signatures of these instabilities are of equal importance to characterize and study using the field-particle correlation technique. While the instantaneous field-particle correlation technique employed in this study was well-suited to the impulsive ion energization via shock-drift acceleration and the steady electron energization via adiabatic heating, we may require finite time correlations to characterize the energy exchange within the upstream fluctuations and turbulence of the shocks in exact analogy with previous field- Ex (vx, vy) using the full v 2 weighting, both computed from the ion distribution function of the self-consistent Gkeyll simulation.
particle correlation studies (Klein et al. , 2020Horvath et al. 2020). Nevertheless, an exciting frontier awaits in applying the field-particle correlation technique to the distribution functions produced by more realistic collisionless shock simulations and classifying the observed velocity-space signatures of particle energization.

Appendix A. Component-wise Separation of the Field-Particle Correlation Technique
In Eq. (3.4) derived in Section 3, the field-particle correlation was expressed as a dot product between the electric field and the velocity-space gradient of the particle distribution function. While this gives the total energy exchange, it is often useful to identify components of the energy exchange by decomposing the field-particle correlation technique like so, similar to Eq. (3.5) and Eq. (3.6), but without the substitution of the components of v 2 . We now justify this additional substitution to obtain the form of the field-particle correlation technique employed throughout this manuscript. Although this substitution alters the rate of change of phase-space energy density as a function of velocity space (v x , v y ), the difference in these two forms vanishes upon integration of the field-particle correlation over velocity space. In other words, the change does not alter the net rate of particle energization at a given spatial position x 0 . That this replacement does not alter the net rate of energization is easily seen by examining the x contribution to the dot product in the second term of Eq. (3.2) when integrated where the v x integral of the factor with v 2 y is a perfect differential, and thus contributes nothing assuming appropriate velocity-space boundary conditions, lim vx→±∞ f (v x , v y ) = 0.
The primary motivation for this substitution of the components of v 2 is to highlight the regions in velocity space that contribute to the net energy transfer between the particles and fields. Let us compare the velocity-space signature obtained with Eq. (A 2) to Eq. (3.6) for the case considered in Section 4. In Figure 13, we plot the two forms of the field-particle correlation (a) Eq. (3.6), C Ey , and (b) Eq. (A 2), C (v 2 ) Ey for the same case shown in Figure 3, panel(e). Using the alternative form in Figure 13 Ey given by Eq. (A 2), we see that there is a large feature in the velocity-space signature of the ion energization associated with the incoming ion flow, but that significant feature leads to zero net ion energization. In fact, apparent energy transfer associated with E y in this form is actually canceled exactly by the magnetic field term (v × B) x ∂f /∂v x in the Lorentz force, as discussed in Appendix B.
Using the preferred form in Figure 13(a) C Ey , this net zero energy transfer associated with the incoming ion beam does not appear. Only the energy transfer associated with the reflected ions appears when using the form in Eq. (3.6). Therefore, although using only the v 2 y contribution does not capture the full energy flow in velocity space, it does capture the energy transfer associated with the net rate of energization, and so this form is preferable for the study of particle energization.

Appendix B. Calculation of Field-Particle Correlation for the E × B Drift
Here we calculate the field-particle correlation for the rate of change of phase-space energy density of a plasma undergoing uniform E×B motion. Consider the case, relevant to the particular transverse magnetized shock problem addressed here, of a constant transverse magnetic field B = B z0ẑ and a constant electric field E = −E y0ŷ where E y0 > 0, giving an upstream E×B velocity of u E×B = −(E y0 /B z0 )x. The 2V Maxwellian distribution drifting with this u E×B velocity is given by where we have assumed no spatial variation, in analogy with the upstream region of the perpendicular shock studied here. With no spatial variation, the rate of change of phase-space energy density w s ( Substituting in the fields and the velocity-space derivatives we obtain the following result where we see that the term from (v × B) y (∂f s /∂v y ) cancels with the contribution from where the change of phase-space energy density due to the electric field in the first term is canceled by the change of phase-space energy density due to the magnetic field acting on the E × B flow. Importantly, Eq. (B 6) demonstrates that the instantaneous rate of change of the phasespace energy density, at every point in velocity space, is zero for a Maxwellian plasma simply undergoing uniform E × B motion. Of course, we expect that an E × B flow produces no net energization. But, in combination with our intuition that E × B flows produce no net energization, the result presented here more strongly motivates the form of the field-particle correlation technique presented in Appendix A, and other transformations employed throughout this study to eliminate the contribution of E × B flows to individual components of the energization, such as the transformation to the transversedrift frame in Section 5. Using these transformations, we can then gain further insight into the energy exchange in phase space without having to sum over components, as is required in Eq. (B 6) to completely cancel the E × B contribution to the rate of change of the phase-space energy density. x/di −2 in the idealized shock model yields (a) the distribution function fi(vx, vy) and (b) the field-particle correlation CE y of the averaged distribution. The averaged field-particle correlation is approximately symmetric, corresponding to zero net energization, in agreement with a spatial average of jyEy in Figure 14 (c), which suggests that the ions experience no further energization once they have crossed downstream.
we have marked with x A the point that the field-particle correlation was calculated in Section 4.
We now plot in Figure 15 that spatially averaged downstream distribution function and field-particle correlation, averaging over −4 x/d i −2. We observe a broadened ring distribution and an approximately anti-symmetric velocity-space signature, corresponding to zero net energization over this spatial interval. Importantly, if we transform to the downstream frame of reference and compare the distance to the origin of the upstream distribution (white circle) to the distance to the origin of the ring distribution, we find that both distributions, upstream and spatially-averaged downstream, are roughly equidistant to the origin. In other words, the energy of these two distributions is roughly equivalent. In the upstream, the energy is dominantly bulk kinetic, while in the downstream the energy is mostly nonthermal, an increase in the effective perpendicular temperature of the distribution. But, this energy conversion is conservative, not changing the total net energy of the ions. The only process in the simplified model which changes the total microscopic kinetic energy of the ions is the energization via shock-drift acceleration.
It is worth expanding upon this subtlety of energy conversion versus energization by considering a more generic idealized model in which we vary the amplitude of the magnetic field increase at the magnetic discontinuity. In general, as an ion E × B drifts through a magnetic discontinuity, the perpendicular velocity in the local bulk-flow frame of reference increases at the expense of the diminished bulk-flow E × B velocity. The downstream perpendicular velocity relative to the upstream bulk-flow velocity v ⊥d /U u is determined by three dimensionless parameters for this idealized problem: (i) the ratio of the downstream to the upstream magnetic field magnitude B d /B u ; (ii) the ratio of the upstream perpendicular velocity to the upstream bulk-flow velocity v ⊥u /U u ; and (iii) the gyrophase θ of the ion † when it first reaches the magnetic discontinuity.
In Figure 14(b), the specific ion trajectory plotted returns upstream (red segment) due to the increased magnetic field downstream of the discontinuity. If the ion does not return upstream, then one can compute the downstream perpendicular velocity v ⊥d /U u as the difference between the velocity upon crossing the discontinuity and the downstream E×B velocity, yielding Note that, although the perpendicular energy relative to the local (upstream or downstream) bulk-flow velocity generally increases, this increase comes at the expense of the kinetic energy of the incoming bulk flow, and the total kinetic energy of each ion is conserved in this process. This statement can be proven for a ring of ions with perpendicular velocity v ⊥u about upstream velocity U u by squaring Eq. (C 1), substituting B d /B u = U u /U d , integrating the gyrophase θ over 2π, and multiplying by m i /2, to obtain the expression The conservation of energy is obvious when evaluated in the downstream frame (U d = 0), where Eq. (C 2) proves that the downstream perpendicular energy of the ring of ions is simply the sum of the upstream perpendicular energy plus the "bulk" kinetic energy of the ring distribution moving at U u . Of course, if an ion does return upstream, it can gain energy by the process of shockdrift acceleration via the alignment of transverse (to the shock normal) component of the Larmor velocity and motional electric field that supports the inflow at the E × B velocity. Only if an ion returns upstream is there any net energy exchange between the electromagnetic fields and the particles. We demonstrate this energy gain for the idealized model in Figure 16, which plots the gain of perpendicular energy (v ⊥d /v ⊥d,th ) 2 as a function of (v x , v y ) for B d /B u = 2, 3, 4 in panels (a,b,c).
In this figure, the black contours separate regions with different numbers of crossings of the magnetic discontinuity (x = 0 in Figure 14), where every ion must cross the magnetic discontinuity an odd number of times to eventually cross downstream, given by the large numbers on the plot. The increase of the perpendicular energy due to shock drift acceleration is given by the colorbar. In the self-consistent simulation, this energy gain comes at the expense of the field energy, while for the idealized model, this energy gain via shock-drift acceleration is not conservative. Nevertheless, the idealized model helpfully illustrates where in phase space we observe merely energy conversion versus where we expect to see actual energization due to the electromagnetic fields. All ions that cross the magnetic discontinuity only once conserve their energy, leading to the increase of the perpendicular energy (relative to the downstream frame) predicted by Eq. (C 1) (green color).
We conclude this appendix by making a few final notes on this distinction between energy conversion and energization. The movement of particles from one position in velocity space to another requires acceleration, so this energy conversion is still mediated by forces in the plasma: the v × B force in the single-particle picture, which combines with the v · ∇ x streaming term in the distribution function picture. The latter term corresponds to the traditional picture of pressure work, as once we have a distribution of particles, the pressure can participate in this energy conversion.
Therefore, θ = 0 corresponds to a perpendicular velocity that increases the magnitude of the inflow velocity, and θ = 180 • decreases the magnitude of the inflow velocity. It is important though to distinguish pressure work, which simply converts one form of energy into another, and a pressure supported electric field, which requires gradients in the pressure and can energize the plasma. For example, the cross-shock electric field which arises in the self-consistent simulation is a result of the electron pressure gradient. This pressure-supported cross-shock electric field both increases the reflection of ions-see Appendix E for further details-and is a critical component to the increase in T ⊥ of the electrons via adiabatic heating, where the pressure gradient provides the relevant drifts for the electron distribution's adiabatic invariant to be conserved through the magnetic field gradient.
Fundamentally, we seek to be as precise as possible in what we are diagnosing with the field-particle correlation technique by focusing exclusively on the electric field component of the evolution of the phase-space energy density. While the energy conversion that occurs in collisionless shocks is a component of the overall increase in the temperature of, e.g., the ions, this process of energy conversion is distinct from the energization processes that occur, such as shock-drift acceleration. And it is these energization processes that we seek the velocity-space signatures of, as we may then be able to leverage this same

Appendix E. Effect of Finite Ramp Width and Cross-Shock Electric Field on Ion Energization
To understand the effects of the finite shock width and the cross-shock electric field on the ion energization, we first use the Vlasov mapping model to separate out the effects of different components of the electric field on the ion trajectories, similar to the analysis of the electrons presented in Section 5.3. In Figure 17(b), we compare the ion trajectories between two Vlasov-mapping models: (i) the "full model" (solid) which computes the ion trajectories and evolution of the ion velocity distribution using the full electromagnetic fields from the Gkeyll simulation; and (ii) the "zero E x model" (dashed), in which we artificially set the cross-shock electric field to zero. In Figure 17(c), we plot the rates of energization of the ion distribution by the electric field, j x E x (blue) and j y E y (red), along with the total energization, j·E (black), for the two models. We also show the cumulative total energization of the ion distribution by integrating from upstream x xup j · Edx, along with the separate contributions from each component of the electric field for both models in Figure 17(d).
First, we adopt a Lagrangian approach to examine the ion energization along its trajectory. The comparison of the example reflected ion trajectories in Figure 17 (c) Rate of work done by the components of the electric field, jyEy (red) and jxEx (blue) for the full model (solid) and zero Ex model (dashed), along with total j · E (black). (d) Cumulative work done integrated from the upstream x j · E. Inclusion of the cross-shock electric field enhances ion reflection, thereby achieving a larger energy gain due to the motional electric field Ey through shock-drift acceleration.
how the cross-shock electric field alters the ion trajectory. Because E x opposes the flow of ions into the shock, the ion penetrates less deeply into the shock before turning back upstream (solid blue segment of trajectory) than for the zero E x model (blue dashed). For the full model including E x , the ion returns further upstream (red solid), where the lower magnitude of the magnetic field leads to a larger Larmor radius of its orbit. This return further upstream in the full model is particularly important when the shock ramp has a finite width. The enhancement of the ion reflection by E x significantly affects energization of these reflected ions by E y through the shock drift acceleration mechanism, where Figure 17(c) shows that the rate of ion energization j y E y (red) in the foot and ramp region, 22 x/d i 24, is much larger for the full model (solid) than for the zero E x model (dashed). This increased energization is a direct result of the larger distance the full model ion traverses in y upstream as its gyroradius is increased by the combination of acceleration by E x and the decreased magnetic field amplitude upstream.
Another way to view the effect of the cross-shock electric field in increasing the efficiency of shock drift acceleration is to employ a complementary Eulerian point of view to examine the energization as a function of velocity space (v x , v y ) at a single point in configuration space. Following this approach, we explore the enhanced reflection due to the cross-shock electric field by examining C Ex to understand how E x accelerates or decelerates ions in different regions of velocity space. In Figure 18, we plot (a) the ion distribution function f i (v x , v y ) and (b) the correlation with the cross-shock electric field C Ex (v x , v y ) from the simulation at the position x B = 21.8d i (vertical red line in Figure 17 and the same point where the electron analysis in Section 5 was performed) where the cross-shock electric field peaks.
The ion distribution at this position is dominated by the incoming beam, with a small fraction of reflected ions forming a "boomerang" shaped distribution. The dominant effect is that E x decelerates the incoming ion beam. But, the population of reflected ions with v x > 0 at x B -which corresponds to the upper crossing of the red segment of the trajectory with the vertical line at x B in Figure 17(b)-is being accelerated by E x , causing these ions both to return further upstream and to increase their perpendicular velocity, thereby leading to a larger Larmor radius. These two effects E x has on the reflected ions with v x > 0 reinforce the enhanced reflection and increased energization of these ions by the shock drift acceleration mechanism.
In this regard, we reiterate a powerful feature of the FPC: the velocity-space signatures produced by the FPC reveal how electric fields energize different components of the ion distribution in qualitatively different ways. The cross-shock electric field decelerates the incoming beam while accelerating the reflected population with v x > 0, as shown by the blue and red signatures respectively in these regions of phase space. The separation of the energization of different populations of the ion velocity distribution from an Eulerian perspective, provided by the FPC method, enables a deeper understanding of the underlying mechanisms of ion acceleration at the shock. We emphasize that by looking only at the velocity-integrated rate of energization by E x -given by j x E x in Figure 17(c) at x B -one sees just the net loss of ion energy due to E x , masking the important effect that the cross-shock electric field plays to enhance the ion reflection. While the role of the cross-shock electric field in enhancing the reflection of the ions has been previously theorized to be an important component of energizing the reflected ion population (Cohen et al. 2019), the Eulerian perspective provided by the FPC makes the physics of the cross-shock electric field especially clear by illustrating where the ions are gaining and losing energy in phase space.  Figure 18. The ion distribution function from the Gkeyll simulation (left), and CE x computed from the Gkeyll simulation (right) plotted at xB = 21.8di, near the peak of the cross-shock electric field. We note two features in the velocity-space signature found from computing CE x : the strong negative correlation coincident with the incoming beam, denoting the deceleration of the incoming flow and transfer of energy from the bulk upstream kinetic energy to electromagnetic energy, and the modest positive correlation at vy < 0, vx > 0 where particles can now be accelerated by the cross-shock electric field and pushed back upstream. This acceleration of ions of particular velocities is the principal reason for the increased efficiency of shock-drift acceleration despite the finite shock width, as the cross-shock electric field assists in increasing the phase-space density of reflected ions that can gain energy along the motional electric field upstream.
If the plasma is magnetized, or at least the electrons are as in Section 5, it is natural to split the pressure tensor as where P C s = (I − bb)p s,⊥ + bbp s, = Ip s,⊥ + bb(p s, − p s,⊥ ), (F 7) is the Chew-Goldberger-Low (CGL) pressure tensor (Chew et al. 1956), b = B/|B| is the direction of the magnetic field, and Π s is the agyrotropic part of the pressure tensor. Note that Tr (P C s ) = 2p s,⊥ + p s, = 3p s , where p s is the scalar pressure. From this definition, we can also see that Tr (Π s ) = 0. In the 1D-2V simulation of interest in this study, where B = B z (x)ẑ, the trace of the pressure tensor is instead Tr (P C s ) = 2p s,⊥ = 2p s because we are not evolving the degree of freedom parallel to the magnetic field. Thus, p s,⊥ = p s in this geometry, but for generality we will retain the subscript ⊥ for the remainder of the derivation.
The divergence of the CGL pressure tensor is ∇ · P C s = ∇p s,⊥ + (p s, − p s,⊥ ) ∇ · (bb) (∇·b)b+∇ b +b∇ (p s, − p s,⊥ ), (F 8) where ∇ ≡ b · ∇. Hence, we can calculate the contribution of the CGL pressure tensor to the bulk drift − ∇ · P C s × B q s n s |B| 2 = − ∇p s,⊥ × B q s n s |B| 2 + (p s,⊥ − p s, ) ∇ b × B q s n s |B| 2 , (F 9) where the terms in the direction of the magnetic field in Eq. (F 8) are eliminated upon taking the cross product with B. Putting everything together, we obtain u s⊥ = E × B |B| 2 − ∇p s,⊥ × B q s n s |B| 2 + (p s,⊥ − p s, ) ∇ b × B q s n s |B| 2 − ∇ · Π s × B q s n s |B| 2 − m s q s |B| 2 du s dt × B.

(F 10)
We now define the magnetization vector (Hazeltine & Waelbroeck 1998), which is a generalization of the definition in Eq. (5.11). We note that The first three terms, the E × B drift, the magnetization drift, and the ∇B drift † are the dominant three drifts in the 1D-2V perpendicular shock of interest in this study. The † We can see this is identical to the ∇B drift definition employed in Eq. (5.7) for a magnetic field only in the z direction with a bit of vector calculus, ∇ × (ẑ/Bz) = −∇Bz ×ẑ/B 2 z .
other terms: the curvature drift, agyrotropy drift, and polarization drift are all either identically zero in this geometry or small. For example, we demonstrated in Section 5 that the polarization drift is small through the shock, and because the electron's adiabatic invariant is well conserved the agyrotropic component of the drift must be small. We note again that the combination of the bulk ∇B drift and the magnetization drift produce the familiar diamagnetic drift, and that an alternative interpretation of the results of Section 5 is that the electron distribution's adiabatic invariant is conserved via the alignment of the diamagnetic drift with the motional electric field. In other words, whereas for a single particle only the ∇B drift was important for that single particle to heat adiabatically, the generalization to a distribution of particles leads to a bulk drift, the diamagnetic drift, which is a combination of the ∇B drift and the magnetization drift, being the principally important drift for the distribution of particles to heat adiabatically.
Appendix G. Checking µ e Conservation with a m i /m e = 400 Simulation Here, we repeat some of the analysis of Section 5 for a more realistic mass ratio simulation, m i /m e = 400, to determine a possible source for the slight disagreement between the energization due to the ∇B and the magnetization drifts and the energization, j e · E computed from moments of the electron distribution function. All other parameters are the same, e.g., box size, L x = 25d i , plasma betas, β i = 1.3, β e = 0.7, and electronelectron collisionality, ν ee = 0.01Ω ci . Note that with the increased mass ratio, the ion-ion collisionality is commensurately reduced. In addition, because the ions are more massive and there is more scale separation between the ions and electrons, we have doubled the configuration space resolution to N x = 3072 to keep ∆x ∼ d e /6. For computational convenience, we perform our analysis just after the shock is formed, t = 4.3Ω −1 ci . We plot in Figure 19 an identical figure to Figure 11 for the m i /m e = 400 simulation to compare strengths of the same drifts of interest in Section 5 in panel(a): E × B in x and y, the ∇B drift in y, the magnetization drift, ∇ × M, in y, and the polarization drift in x. We also repeat the comparison of these drifts to the computed first moment from the electron distribution function (b), alongside a comparison of the amount of bulk energization arising from these drifts (c), and how it compares with the bulk energization, j e · E computed from moments of the electron distribution function (d). We note that the agreement between the energization of the electrons arising solely from the alignment of the ∇B and the magnetization drifts with the motional electric field and the total j e · E computed from moments of the electron distribution function is better than what was observed in Figure 11 for the m i /m e = 100 simulation. The more realistic mass ratio increases the scale separation between the shock-width, which remains L shock ∼ d i , and the electron gyroradius, and thus we expect the electron adiabatic invariant to be more strongly conserved through the shock. This stronger conservation is indeed the case, as we show in Figure 20 comparing µ e computed from both the m i /m e = 100 and m i /m e = 400 at the same time t = 4.3Ω −1 ci . (c) (d) Figure 19. A comparison of the strength of the major single-particle drifts through the shock (a), E × B in x (black) and y (green), the ∇B drift in y (blue), magnetization drift in y (red dashed), and the polarization drift in x (magenta. dashed-dotted) for a mi/me = 400. We check that these drifts sum to the total first moment computed from the electron distribution function (b) as well as determine how each of these drifts contributes to the overall energy exchange, je ·E (c), and compare the je ·E computed from these drifts to the total je ·E computed from moments of the electron distribution function. As before in Figure 11, we sum the energy exchange due to E × B flows to demonstrate that this total energization is zero, as it should be. We note that the energization due to the combination of the ∇B and magnetization drifts more closely agrees with the energy exchange, je · E computed from moments of the electron distribution function in comparison to Figure 11. . The electron adiabatic invariant, µ = T ⊥ /Bz, for the mi/me = 100 simulation (blue) and mi/me = 400 simulation (red dashed). The conservation of the adiabatic invariant is within ∼ 1 percent in the mi/me = 100 simulation, while the conservation is even better for the mi/me = 400 simulation, suggesting that the mi/me = 400 is even more strongly in the asymptotic limit of ρe ≪ L shock .