Magnetic pumping model for energizing superthermal particles applied to observations of the Earth's bow shock

Energetic particle generation is an important component of a variety of astrophysical systems, from seed particle generation in shocks to the heating of the solar wind. It has been shown that magnetic pumping is an efficient mechanism for heating thermal particles, using the largest-scale magnetic fluctuations. Here we show that when magnetic pumping is extended to a spatially-varying magnetic flux tube, magnetic trapping of superthermal particles renders pumping an effective energization method for particles moving faster than the speed of the waves and naturally generates power-law distributions. We validated the theory by spacecraft observations of the strong, compressional magnetic fluctuations near the Earth’s bow shock from the Magnetospheric Multiscale mission. Given the ubiquity of magnetic fluctuations in different astrophysical systems, this mechanism has the potential to be transformative to our understanding of how the most energetic particles in the universe are generated.

S uperthermal populations of ions and electrons are abundant in a wide variety of astrophysical systems throughout the universe, often with distributions characterized by energetic power-law tails 1,2 . However, the consensus on the physical mechanisms that energize these particles is far from settled. While it is well known that plasma can be energized by waves, most theories of wave-particle energization (excluding shocks and magnetic reconnection) are only effective for energizing particles moving at velocities close to the phase velocity of the waves [3][4][5] . We here present an analysis of particle energization by magnetic pumping. While previous work suggests that this mechanism is only effective up to the phase velocity of the wave 6 , we show that the addition of magnetic trapping of particle orbits renders pumping effective for heating particles moving far faster than the wave speed, an example of which is shown in Fig. 1. This mathematical treatment reveals the underlying Fermi heating mechanism consistent with the formation of energetic power-law distributions.
Diffusive shock acceleration (DSA) 7 is a classic example of such a Fermi heating mechanism. Given an initial seed population, DSA yields power-law distributions of often relativistic electrons with Larmor radii larger or comparable to the width of the shock front. In this Letter we find that magnetic pumping in the fluctuations generated in the vicinity of Earth's bow shock can provide a seed population of pre-energized magnetized electrons. Thus, along with mechanisms such as stochastic shock drift acceleration (SSDA) [8][9][10] , magnetic pumping may help address the injection problem of DSA 11,12 .
Magnetic pumping is directly related to the pressure anisotropy that naturally forms when a plasma and its magnetic field are compressed 13 . This anisotropy can be moderated by an effective scattering frequency, ν, caused by processes such as pitch-anglemixing by Whistler waves 14 , or a limited confinement time of the electrons in the magnetic wells, yielding a phase delay between the perpendicular pressure, p ⊥ , and the flow perpendicular to the magnetic field, u ⊥ . Given this phase delay, when averaged over a pump cycle, mechanical work by p ⊥ ∇ ⊥ ⋅ u ⊥ then becomes finite and is the source of the energy for the magnetic pumping process. In this process energy is transferred from magnetic fluctuations to directly heat the plasma, bypassing the turbulent cascade [15][16][17] .
Extensive prior work, however, suggests that magnetic pumping is not effective for energizing superthermal particles. For example magnetic pumping by compressional waves, related to transit-time damping 18 , is shown to be a Landau damping process, where heating is limited to particles moving at the magneticfield-aligned phase velocity of the wave considered, v p = ω/k ∥ , where ω is the angular frequency of the wave and k ∥ is the wavenumber in the direction parallel to the magnetic field 19 . In the framework of quasilinear theory Landau damping causes velocity diffusion limited to particles near the resonance velocity, v~ω/k ∥ , and is derived based on the standard procedure of integrating the plasma kinetic equations along unperturbed particle trajectories. In contrast, here we consider a standing wave geometry and apply the fast transit-time limit 20 to retain anisotropic effects related to the full electron orbit motion.
Here we show that a quasilinear analysis with trapped electron effects and v ≫ ω/k ∥ yields a velocity diffusion equation similar to that obtained by Lichko 2017 6 for the opposite limit, v ∥ ≪ ω/k ∥ . More specifically, the slowly-varying background distribution f 0 is governed by a diffusion equation of the form where G is independent of v, but is a function of ν/ω, and the magnitude of the magnetic perturbations relative to the background magnetic field, δB/B 0 . The result that D ∝ ωv 2 is evidence of a Fermi heating process with a diffusive step size proportional to the velocity, Δv ∝ v.

Results
Model of electron trapping. Below we will outline how Eq. (1) is obtained and provide an evaluation of G, a metric of the effectiveness of the pumping process. Meanwhile, we will compare these predictions to observations by the Magnetospheric Multiscale (MMS) mission in the region of the Earth's bow shock 21 .
Within the electron foreshock there are ripples, variations in the magnetic field itself, that have been shown to be a source of electron acceleration, as well as other large-amplitude magneticfield fluctuations [22][23][24][25][26] , an example of which was recorded on October 7th, 2015 27 . From the time histories of |B|, n, T obs , and T obs /T adiabatic , as shown in Fig. 2a-d, the temperature increase is greater than would be expected from compressional heating alone. The evolving pitch-angle-averaged distribution functions in the foreshock are shown in Fig. 2e for the times marked in Figs. 2a-d. These distribution functions demonstrate energy transfer consistent with a Fermi heating mechanism, where Δv ∝ v, which in log-log format can be seen in the shift of the energetic power-law tails at a constant slope.
We first aim to demonstrate that the observed anisotropy is consistent with electron trapping in a standing wave perturbation and is representative of the perturbed distribution function driven during each pumping cycle. We consider the limit where the bounce time, τ b , is much smaller than the time scales associated with the waves. The instantaneous particle orbits are then described by the magnetic moment, μ ¼ mv 2 ? =ð2BÞ and the total energy U ¼ E À eΦ, where E ¼ 1 2 mv 2 , e is the positive electron charge, and Φ is the electrostatic potential. We assume electron energies larger than the electrostatic potential associated with the perturbations 28 such that the v × B part of the Lorentz force dominates the orbit motion. Using the additional assumption that electrons are well-magnetized, with the Larmor radius smaller than the perpendicular wave length of the perturbations, ρ L < λ ⊥ , for the parameters in this bow shock crossing the model is appropriate for energies 100 eV ≲ E ≲ 100 keV. At each time point during the magnetic pump cycle an instantaneous electron orbit is fully characterized by μ and E. In turn, from the multiple time scale method 20,29 the distribution must be constant along these instantaneous orbits, reducing the dimensionality of the problem.
Starting with the drift kinetic equation 30 with pitch-angle mixing, i.e., df dt we change variables from f ðt; x; v ? ; v k Þ ¼ f ðt; E; χÞ, where χ = Λ/ (j 2 + Λ), Λ ¼ μB 0 =E, and j ¼ J=ð4vLÞ ¼ 1=ðvLÞ R L b 0 v k dx with L b denoting the bounce point and L the length of the flux-tube element. Here E and χ are both constant of motion variables where χ is representative of v 2 ? =v 2 along the instantaneous orbits.
Using these new variables and following the approach of Montag et al. 31 and Egedal et al. 32 , we obtain an orbit-averaged form of Eq. (2) ∂f ∂t À Hðt; χÞE ∂f ∂E where H ¼ ð2=jÞð∂j=∂tÞj χ , and the orbit average ð:::Þ h i x is defined in the Methods subsection Orbit averaging of the Lorentz operator.
Electron trapping in spacecraft data. We solve Eq. (3) numerically, assuming an initial isotropic distribution and a standing wave magnetic field, Bðx; tÞ Bðx; tÞ=B 0 ¼ 1 À ðδB=B 0 Þ sinðωtÞ cosðk k xÞ.  Fig. 3d. A more detailed description of this matching process is detailed in the Methods subsection Obtaining an estimate for the effective scattering frequency, as well as ref. 33 .
The agreement demonstrated above suggests that the model is capturing the anisotropic features of the observed distribution functions. While the inferred amount of scattering is too low for SSDA 10 to be effective, it is near optimal for magnetic pumping and the analysis can now be extended to address the energization of the electrons over many cycles. Following the blueprint of the quasilinear method, we then separate the distribution function into the slowly-varying, isotropic background distribution, f 0 , and the anisotropic portion of the distribution function, f 1 , where the pitch-angle averaging ð:::Þ h i χ is defined in the Method subsection Particle conservation in ð:::Þ h i χ . To make Eq. (3) more analytically tractable, the Lorentz operator is approximated with the Krook operator,  found, Inserting f 1 back into Eq. (3), an evolution equation is obtained for the slowly-varying background distribution More explicitly we recover Eq. (1), with where we recall G is a measure of the energization from a single pump cycle. An approximate form of G that is easy to evaluate is given in the Methods subsection Fitting the results of G.
We validate the analytical model of Eqs. (1) and (6) 7), where T = 2π/ω and the numerator and denominator are found to be nearly linearly dependent functions of v before averaging. In Comparison model to spacecraft data. Using the earliest spectrum (yellow) selected in Fig. 2e with a small amount of smoothing as the initial condition, we apply the model in Eqs. (1) and (6) to predict the evolution of the pitch-angle-averaged distributions recorded by MMS. Compressional heating is modeled by including À _ n=ð3nÞvð∂f =∂vÞ on the right-hand-side of Eq. (1) 6 . The result of this calculation is shown in Fig. 4b and is in good agreement with the observation from Fig. 2e, repeated for convenience in Fig. 4c. While there are some differences at lower  energies, likely due to the effects of electric fields, which are neglected in our approach, the model provides a good account for the energization of electrons at large energies. The green dashed lines in Figs. 2e and 4c are the final spectra when considering compression alone, clearly underestimating the level of energization. Downstream of the shock front (for t > 11:44:44UT) the energization rate by pumping declines with the frequency of the fluctuations, and in combination with parallel streaming losses is consistent with a drop in the fluxes of superthermal electrons.

Discussion
We have here presented an energization mechanism, magnetic pumping, that becomes applicable to superthermal particles, v ≫ ω/k, when the effects of trapping are retained. The MMS observations provide evidence that magnetic pumping has a significant role in electron energization in the region of the Earth's bow shock. Given the potential universal applicability of the model, this could have a far-reaching impact on our understanding of electron and superthermal ion energization in many other plasma environments where particles with v ≫ ω/k are observed, such as the solar corona, cosmic ray generation pumped by magnetic turbulence in the interstellar medium, or possibly shocks driven by supernova explosions.

Methods
Finding j as a function of χ. The kinetic description applied in this work is limited to the superthermal particles characterized by speeds, v, sufficiently large that the Lorentz force is dominated by the magnetic term, i.e., vB ≫ E. To be more specific, in our drift kinetic analysis we are concerned with the parallel motion along the magnetic field, in general governed by forces due to the parallel electric field and the magnetic mirror force. For plasma variations of scale length L, the magnitude of these forces can be estimated as |e∇ ∥ Φ| ≃ T e /L and |μ∇ ∥ B| ≃ mv 2 δB/(2B 0 L), and it follows that the superthermal limit requires v 2 ) v 2 t B 0 =ðδBÞ, where v t is the electron thermal speed and δB/B 0 is the normalized magnetic fluctuation amplitude.
In the superthermal limit the description of the orbit motion is significantly simplified as the value of v k =v ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 1 À ΛB p along an orbit is only dependent on Λ ¼ μB 0 =E (and independent of the electron energy E with the assumption E ∥ = 0). This strongly simplifies the calculation of the second adiabatic invariant J(v, Λ) because j = J/(4vL) is then a function of only Λ, readily evaluated numerically for the slowly evolving magnetic perturbation considered, as illustrated in Fig. 5, where jðΛ; tÞ ¼ 1 4L Given this calculation of j(Λ, t) over the course of a full pump cycle we can determine Λ(χ, t) and in turn obtain j(χ, t) as shown in Fig. 5c. This function is fundamental to our analysis as it is related to the instantaneous rate of particle energization. Because dJ/dt = d(vj)/dt = 0 it follows that j dv/dt + v dj/dt = 0. Furthermore, as j = j(χ, t) and dχ/dt = 0, we have dj=dt ¼ ∂j=∂tj χ , such that v À1 dv=dt ¼ Àj À1 ∂j=∂tj χ . We then obtain the result implicit in Eq. (3) that Particle conservation in ð:::Þ h i χ . Without loss of generality, in our analysis magnetic-field lines in the (x, z)-plane are characterized by the flux function Ψ, and we consider a flux-tube with constant width Δy in the y direction. The width in the z direction varies as 1/B and is parameterized by ΔΨ. The total number of particles within our flux tube of length L must be conserved, i.e., For the point x = L/2, B = B 0 we can evaluate ΔΨ = B 0 Δz, where Δz is the width of the flux tube at that x location. Furthermore, becauseτ b ¼ vτ b =ð4LÞ and dEdμ ¼ m 2 v 3 dvdΛ=B 0 we get Rewriting this in terms of our preferred variables, (v, χ) the expression becomes This motivates the averaging operator ð::: such that the pitch-angle averaged distribution becomes In turn, the total number of particles can then be written as From the form of Eq. (13) it is clear that Fv 2 dv is proportional to the number of particles in the flux tube within a differential speed interval, dv.
Obtaining an estimate for the effective scattering frequency. We can estimate how much scattering is needed to match the spacecraft observations by comparing the MMS distributions to theoretical distributions generated by integrating Eq. (3) for a range of scattering frequencies, ν. By comparing the anisotropic features in the MMS distribution to the features in this set of theoretical distribution functions, we can estimate the effective scattering frequency. Based on a least-squared-fit analysis, the scattering frequency that best fit this data varies as a function of velocity, where ν/(ω/2π) ∈ [0.25, 1.5]. A scattering frequency within this range, ν/(ω/2π) = 0.75, is chosen to generate a comparison with the MMS data, where the resulting scattered distribution is shown in Fig. 3(k). A more thorough investigation into how the estimated scattering frequency varies as a function of velocity can be found in ref. 33 , while examples of theoretical distributions scattered at a set of different rates are shown in Fig. 6.
Observations of distributions throughout the encounter. In Fig. 3 of the main text we demonstrated anisotropic distribution with features matching those of the numerical model. To further demonstrate that these type of anisotropic distributions are representative for the full event, we here include in Fig. 7 electron distribution functions measured by MMS along the entire bow shock crossing. Compared to the peak and valley distributions in Fig. 3 these follow the expectations from the theory.
Fitting the results of G. For an easy-to-evaluate approximation of the theory, the results of Eq. (6) can be fit as where C K ¼ 1:15=ðδB=B 0 Þ 1:13 , as discussed in the text and the last part of this methods section. This fitting, plotted alongside the results from Fig. 4a are plotted in Fig. 8. We note that this approximation for G is also valid in the limit of δB/B 0 → 0.
Orbit averaging of the Lorentz operator. As outlined above, our model is averaged over the fast electron orbit motion. Locally we assume that the scattering process is described by the familiar Lorentz pitch-angle scattering operator, L, which is then subject to orbit averaging 20,29 . For this orbit averaging, it is convenient to express L in terms of the constant of motion variable Λ. Starting with the Lorentz pitch-angle scattering operator in terms of ξ = v ∥ /v, we rewrite the scattering operator in terms of the new variables, Λ and E which then takes the form After averaging along the spatial dimension of the flux tube we then obtain the orbit-averaged operator, where ð::: is the integral over the orbits, whereτ b ¼ vτ b =ð4LÞ is the speed-normalized bounce time.
Calibrating the Krook operator as a function of δB/B 0 . As is clear from the form of Eq. (14), the Lorentz operator describes diffusion of anisotropic features of the distribution function, f. Accordingly, the diffusion time scales as τ D ∝ (δξ) 2 , where δξ is the typical scale for the anisotropic features of f. As the magnetic field increases, the trapped portion of the distribution function will increase commensurately, yielding larger values of δξ.
The Krook operator does not have this same dependence on δξ, as the rate of isotropization is in fact independent of δξ. When computing the Krook collision operator the new, scattered distribution function is formed from a linear combination of the original distribution function, and a fully isotropized version of  the distribution function, where α ¼ expðÀνC K ΔtÞ determines the rate of isotropization.
To approximate the δξ dependency of L we have introduced the coefficient C K . From an analysis of the kinetic equation it follows that δξ of f is similar to δξ of g = j 2 + Λ, which we used to provide a calibration for the efficiency of the Krook operator: By repeating this computation for a multiple δB/B 0 we find the result scales as which is the result we used to compute the Krook distribution curves earlier in the paper.

Data availability
All relevant data are available from the corresponding author upon reasonable request. Additionally, the MMS spacecraft observations used in the production of this work can be found in the MMS Science Data Center (https://lasp.colorado.edu/mms/sdc/public/ datasets/fpi/ and https://www.lasp.colorado.edu/mms/sdc/public/datasets/fields/).

Code availability
All relevant code and Matlab routines are available from the corresponding author upon reasonable request.