Electron Heating by Magnetic Pumping and Whistler-mode Waves

The investigation of mechanisms responsible for the heating of cold solar wind electrons around the Earth’s bow shock is an important problem in heliospheric plasma physics because such heating is vitally required to run the shock drift acceleration at the bow shock. The prospective mechanism for electron heating is magnetic pumping, which considers electron adiabatic (compressional) heating by ultralow-frequency waves and simultaneous scattering by high-frequency fluctuations. Existing models of magnetic pumping have operated with external sources of such fluctuations. In this study, we generalize these models by introducing the self-consistent electron scattering by whistler-mode waves generated due to the anisotropic electron heating process. We consider an electron population captured within a magnetic trap created by ultralow-frequency waves. Periodical adiabatic heating and cooling of this population drives the generation of whistler-mode waves scattering electrons in the pitch-angle space. The combination of adiabatic heating and whistler-driven scattering provides electron acceleration and the formation of a suprathermal electron population that can further participate in the shock drift acceleration.

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.periodically heat and cool the trapped electrons (see examples in Yao et al. 2021;Jiang et al. 2022), whereas whistler-mode waves may scatter electrons in the pitch-angle space and provide redistribution of the electron energy between the crossfield and field-aligned degrees of freedom.Lichko et al. (2017) and Lichko & Egedal (2020) have shown that such a combination of adiabatic (reversible) energy variation and pitch-angle scattering should result in electron heating and the formation of a suprathermal electron population.
In this study, we investigate the magnetic pumping mechanism of electron heating for a system where the electron pitch-angle scattering is provided by self-consistently generated whistler-mode waves.The basic elements of the model are provided in Section 2 and Appendix: Section 2.1 describes the electron adiabatic heating/cooling by ultralow-frequency magnetic field variations, Section 2.2 describes the model of the whistler-mode wave generation by adiabatically heated (anisotropic) electrons, and Section 2.3 describes the model of the electron pitch-angle scattering by whistler-mode waves.The combination of these three models provides a full kinetic description of the evolution of the electron distribution.Numerical solutions of the corresponding kinetic equations are given in Section 3, and the numerical scheme is described in Appendix.Then, in Section 4 we discuss the results obtained and provide general conclusions of this study.

Theoretical Model
The theoretical model combines three effects: an electron adiabatic response to ultralow-frequency magnetic field variations, whistler-mode wave generation by adiabatically heated electrons, and electron scattering by whistler-mode waves.The main assumption of this model is the following hierarchy of timescales: the fastest timescale (after the electron cyclotron  (Angelopoulos 2008) observations of ultralow-frequency magnetic field fluctuations upstream of the Earth's bow shock (x GSE = 13R E ) in the top panel.Middle panels show the electron energy spectrum and the spectrum of high-frequency magnetic field fluctuations (white and black curves show the local electron gyrofrequency and the fraction of this gyrofrequency, respectively).The bottom set of panels zoom-in on the interval with whistlermode bursts embedded into local magnetic dips created by ultralow-frequency magnetic field fluctuations.We use magnetic field measurements from the fluxgate magnetometer (Auster et al. 2008) and the search-coil magnetometer (Le Contel et al. 2008).The electron spectrum is obtained from electrostatic analyzer measurements (McFadden et al. 2008).
period that is implicitly included) is the electron bounce period (∼2π/ω b ), the intermediate timescale is the period of adiabatic (periodic) magnetic field variation (∼T ad ), and the slowest timescale is the electron diffusion time (∼1/D αα ).This hierarchy allows consideration of the bounce-averaged adiabatic variations and the diffusion, and determines the diffusion effect as a weak perturbation on a T ad timescale.
We consider a system consisting of two electron populations: a cold dense component with the density n c and the temperature T ∥c = T ⊥c = T c → 0, and an anisotropic hot component with the density n h = n c and temperature components T ∥h = T ∥ , T ⊥h = T ⊥ determined with respect to the background magnetic field B. Taking into account that T c → 0, the evolution (dynamical) equations should be solved only for the electron distribution function f of the hot component, whereas zeroenergy electrons mostly determine the background plasma density.In the solar wind and around the bow shock an inequality pe ce 2 2  w w is well satisfied (ω pe = 4πn e e 2 / m e ≈ 4πn c e 2 /m e is the electron plasma frequency, ω ce = eB/m e c is the electron cyclotron frequency), and therefore, the effects of the local charge separation can be neglected for processes considered on significantly larger timescales than the plasma period pe 1 w ~-.The density of the hot electron population in such a system is the only parameter regulating the magnitude of electron fluxes resonating with whistler-mode waves, and thus contributing to the wave growth.As the linear wave growth rate is proportional to the hot electron density, this parameter controls temporal scales, but should not change the character of the evolution of the electron distribution in quasilinear theory.A more complicated (and probably more realistic) system should include a set of electron populations with different anisotropies, densities, and temperatures (see, e.g., Vasko et al. 2019, for a discussion of a multicomponent electron population in the solar wind).The interplay of such populations may significantly affect the wave spectrum (e.g., Fu et al. 2014;Li et al. 2016;Sauer et al. 2020;Vasko et al. 2020;Frantsuzov et al. 2022) and alter the evolution of the electron distribution.Therefore, to obtain general results about the efficiency of the magnetic pumping, we use a simplified system with a single hot (resonating with waves) electron component, whereas more realistic multicomponent systems should be considered for specific solar wind (or magnetic foreshock) events.

Adiabatic Evolution of the System
For a magnetized plasma, the inequality ω ce τ ? 1 is satisfied where τ determines the effective time of the Coulomb interactions between particles.This implies that we can consider the collisionless Hamiltonian  of the particles (electrons) moving in a strong magnetic field B. This Hamiltonian  can be written with a pair of conjugate variables (r, p): , 1 where A = A(t, r) is the vector potential, B = rot A. We consider the magnetic field B to be quasi-parallel to the z-axis, and thus, we define the vector potential A xB t z 0, , , 0 , where B t z ,  ( ) determines the magnitude of the magnetic field The trapped electron motion is localized near the magnetic field minima, and thus, a Taylor expansion is valid: , where B 0 (t) defines the adiabatic evolution of the magnetic field amplitude, and L is the typical magnetic length of the system determined by B As the Hamiltonian  does not depend on the coordinate y, the momentum p y has a constant value: p 0 . The value of this invariant is determined by the initial momentum along the y-axis, and in the axisymmetric case, p y = 0.
The electron dynamics described by Hamiltonian 2 can be presented as a combination of three types of motion: a cyclotron rotation, an electron's bounce motion, and an adiabatic energy conversion throughout the evolution of the magnetic field B. The adiabatic invariant, corresponding to the cyclotron motion, is the magnetic moment ).In terms of the energy E and the pitch angle α, it can be written as ) energies (Cary & Brizard 2009): The second invariant, corresponding to the averaging over bounce motion (z-integration), I z can be obtained as


From Equation (6) we can obtain an equation for the bounce motion frequency ω b : . For I z conservation, the inequalities ω ce ?ω b ?ω ad have to be satisfied, where ω ad = 2π/T ad is the frequency of the adiabatic evolution of the background magnetic field B.
Having two invariants, μ and η, we now should write the kinetic equation for the probability density distribution function of electrons f.This equation, being written through μ and η, shall describe the electron adiabatic response to the magnetic field variation.The basic kinetic description of the system is based on the Vlasov equation for the probability density distribution function of electrons f (t, r, v) (Sinitsyn et al. 2011): For the individual particle with a trajectory r = r(t), we get the invariant of the motion: Notice that there are two time-independent functions of the coordinate z, the energy E and the pitch angle α-μ and η.Any other invariant of the system can be represented as a function of these two, otherwise, the system of three invariant functions will be overdetermined.Therefore, f satisfies the equation f (t, z, E, α) = f (μ, η).To describe the evolution of the distribution function f, we introduce a local evolution operator U ˆthat corresponds to the contribution of the background magnetic field B to the dynamics of the system: Let us specify the initial condition f (0, z, E, α) as a Maxwellian distribution with a temperature T 0 : 11 To determine the distribution function f (t, z, E, α) at any time t, the conservation of invariants μ and η has to be considered.We have to substitute μ(0, z, E, α) → μ(t, z, E, α) and η(0, z, E, α) → η(t, z, E, α) in the initial condition to perform a time translation: Using Equation (13), we can obtain an equation for the temperature anisotropy: The anisotropy has a maximum at z = 0 and goes down to zero for z → ∞.Notice that the growth of the magnetic field magnitude results in the compression of electrons near the magnetic field minima, and thus, locally the hot electron density n h is not constant: As the prefactor in Equation (15) has to be independent of time t and coordinate z, we can obtain the functional form of the density n h : Figure 2 shows the dynamics of the electron temperature anisotropy and the density given by Equations ( 14) and ( 16).
Starting with an isotropic state (T ⊥ /T ∥ − 1 = 0), the electron distribution reaches the maximum anisotropy at the moment of the strongest magnetic field compression, T ad /2, going up to T T 1 2 1  -= -^at z = 0 when the amplitude of the magnetic field doubles.The anisotropy maximizes at the magnetic field minimum location, z = 0, where plasma compression also increases the density of hot electrons n h (by a factor of 2 5/4 times relative to the initial level).Thus, we expect that such a dense and transversely anisotropic electron population will drive the whistler-mode wave generation.After one period of magnetic field variation, at t = T ad , both the temperature anisotropy and the density return to the initial (background) level.Therefore, without electron scattering by high-frequency waves, the electron heating is well reversible, as expected for adiabatic electron dynamics.

Generation of Whistler-mode Waves
Figure 2 shows that variations of the magnetic field B result in a periodical variation of the temperature anisotropy T ⊥ /T ∥ − 1, which is one of the drivers of the generation of whistler-mode waves (e.g., Sagdeev & Shafranov 1961;Kennel 1966).To model the expected increase in the intensity of the whistler-mode wave and to describe the basic properties of the generated waves, we consider the general dispersion equation for the electromagnetic waves propagating along the background magnetic field B (Krall & Trivelpiece 1973): where ω = Ω + iγ, the index j ä {c, h} is the number of electron populations, and ω cj = ω ce = eB 0 /m e c for any j and f j0 is an unperturbed electron velocity distribution function.
Considering f j0 to have a bi-Maxwellian form, we can simplify Equation (17): Using a linear approximation γ = Ω, we can separate the real part of Equation (18) from the imaginary one.Notice that I , and therefore, we have to inspect the limit of this function for y → 0 + , where z = x + iy: where we introduce Cauchy principal value of the integral that excludes the singular point from the interval of the integration.The equation for Ω can be determined by setting γ = 0 and ignoring the imaginary part of Equation (18): The imaginary part of Equation ( 18) defines the equation for the growth rate γ.We consider the case of the sufficiently small growth rate, γ = Ω, and thus the expansion of Equation ( 18 The argument (ω ce − Ω)/ω Tj controls the magnitude of the growth rate and has to be sufficiently large to make the approximation of the small growth rate valid.For whistlermode waves, the relation Ω ω ce holds, and as T c → 0, we can estimate the upper bound of the applicability condition: 24 The temperature anisotropy T ⊥ /T ∥ − 1 and (b) the evolution of density n h in time, due to the periodic magnetic field B fluctuation: ).The adiabatic effect reaches its maxima at z = 0 where the compression of hot electrons increases the hot electron density n h by a factor of 2 5/4 times.leading to where n e m 4 ).Within the quasi-linear theory, the maximum growth rate γ m = γ(Ω m ) determines the maximum wave intensity B w 2 .The dependency can be approximated by a power law with constants C and A (Kuzichev et al. 2019): To describe the frequency spectral density function S Ω (Ω), we assume that the wave energy spectrum is determined by the linear regime of the wave generation, and thus S Ω ∝ γ: is the wave magnetic energy density and I Ω is the normalization factor.This approximation is valid only for the restricted frequency interval Ω/ω ce ä (0, 1 − T ∥ /T ⊥ ), where γ > 0. If Ω/ω ce ∉ (0, 1 − T ∥ /T ⊥ ), then S Ω (Ω) = 0.The functional dependency of the frequency spectral density function S Ω (Ω) on the frequency Ω is shown in Figure 3.To fully determine S Ω (Ω), we introduce the normalization The frequency of the maximum growth rate, Ω m = ω ce x m , is determined by two parameters: the temperature anisotropy T ⊥ /T ∥ − 1 and the parameter β.In the case of β ? 1, x m is sufficiently small and can be found from the expansion The applicability of this approximation is determined by the relation T T ) .Note that the temperature anisotropy T ⊥ /T ∥ − 1 limits the upper bound of the frequency interval for the frequency spectral density function S Ω (Ω), and therefore, Ω m /ω ce stays sufficiently small even if the expansion is not valid.
With the growth of the temperature anisotropy T ⊥ /T ∥ − 1 and the parameter β, the extremum point Ω m gets closer to zero.

Diffusive Scattering of the Electrons
The last element of our model is the electron diffusive scattering.The local diffusion process is determined by a Fokker-Planck equation (Lyons & Williams 1984): where D ij is the local diffusion coefficient over i and j variables (the pitch angle α and the momentum p).The equations for local diffusion coefficients can be obtained using the quasilinear theory (Kennel & Engelmann 1966;Summers 2005): There is only a single solution to the resonant equation that determines the resonant frequency Ω 0 .If β ? 1, the effective part of the energy spectrum lies in the boundaries of Ω = ω ce .Thus, the main impact of the scattering process comes from the resonant frequencies of Ω 0 = ω ce .In the linear approximation, the following equation is valid: Ω 0 /ω ce = ν 2 = 1.The error of this approximation can be evaluated by the smallness of the factor Ω 0 /ω ce (see Figure 4).We estimate the factors of the diffusion coefficients D ij using the obtained equation and Equation  diffusion with a coefficient D αα , whereas the momentum (energy) diffusion can be neglected: In the magnetic trap, electrons perform bounce oscillations, and throughout the bounce oscillation, the pitch angle α changes significantly, which affects the impact of the local diffusion.Therefore, the Fokker-Planck equation and diffusion rates should be averaged over bounce oscillations (Lyons et al. 1972).Let us define the local diffusion operator D L ˆwith respect to Equation (35): On the timescale of one bounce period there is no adiabatic change in the energy E, and thus, the pitch angle α = α(z) can be derived as sin const 2 a k = .,where α 0 = α(0) is the pitch angle at the magnetic field minima.Without the scattering, the distribution function f is the invariant of the motion: f (t, z, E, α) = f (t, 0, E, α 0 ) = f 0 (t, E, α 0 ).As the scattering is a perturbation of the adiabatic evolution, its contribution can be evaluated within the linearized Equation (35): where z = z(t) defines the trajectory of the electron bounce.On the timescale of the adiabatic evolution, the equation for f 0 has to be averaged over the bounce oscillations.As f 0 does not depend on the coordinate z nor does it change in the zeroth order of the perturbation, the averaged diffusion operator D ˆcan be defined as As the frequency spectral density function S Ω is nonzero only on the interval Ω/ω ce ä (0, 1 − T ∥ /T ⊥ ), the upper boundary of the integration over l can be set as l b , which corresponds to If the root expression is negative, then the diffusion coefficient D αα is zero on the whole integration interval and D D 0 . Equation (42) defines the equation for D 0 ¶ a aa :

Simulation Results
To evaluate the impact of the bounce-averaged diffusive scattering on the distribution function f, we have to solve the following evolution equation: It is convenient to perform the transformation f f ln  because at large energies the value of the distribution function f is exponentially small, and therefore, the solution can be strongly affected by the numerical fluctuations.Thus, the second pitchangle derivative 2 0 ¶ a has to be rewritten as Equation ( 56) includes a nonlinear term that complicates the numerical integration.To ensure the stability of the numerical solution, a two-step scheme can be used (see Appendix).To analyze the change of the distribution function f due to the scattering, we introduce the energy distribution function f E and the pitch-angle distribution function f We solve Equation (55) with an initially Maxwellian distribution f and system parameters ω pe /ω ce (t = 0) = 100, n h (t = 0)/n e = 10 −3 , and T 0 = T ∥ (t = 0) = 1 keV.Figure 5 shows the evolution of f E (t, E), the electron temperature anisotropy and the frequency spectral density function S Ω (Ω/ω ce ).There is a periodical electron heating due to the magnetic field compression, and each such heating provides the electron anisotropy increase that drives the whistler-mode wave generation.The electron pitch-angle scattering by such waves is much slower than the adiabatic magnetic field variation, and there is only a small change in the energy distribution f E within the 10 periods shown.The irreversible effect of the electron scattering is better seen in the decrease of the anisotropy peak from the first to the tenth adiabatic period.If the adiabatic magnetic field variation is due to ultralow-frequency compressional waves around the bow shock, then T ad ∼ 0.1-10 s, about or below the typical ion gyroperiod.Therefore, we may consider a hundred such periods, corresponding to tens of minutes of the real time.
Figure 6 demonstrates the time evolution of the energy distribution function f E and the pitch-angle distribution function f 0 a .Over the period of 1000 adiabatic oscillations, the diffusive scattering of electrons isotropizes the system, lowering the value of the temperature anisotropy T ⊥ /T ∥ − 1, whereas the density of the highly energetic particles (with energy E > 10T 0 = 10T ∥ (t = 0)) grows significantly.The isotropization is represented in the pitch-angle space as an inflow of electrons in the areas of low density.This process is related to the functional dependency of the averaged diffusion coefficients D 0 0 a a and D 0 a on the energy E and the pitch angle α 0 : the nonzero area is bounded from above in the pitch-angle space for E = const.Due to the decrease in electron anisotropy, the wave growth rate gradually decreases as well, i.e., the system is reaching its stable state.Within the first 200 adiabatic periods T ad , the density of energetic particles with E/T 0 ∼ 50 increases by a factor of 50 times.Therefore, the combination of the adiabatic pumping and the self-consistent scattering by whistler-mode waves provides a quite effective electron acceleration.
To check if the initial electron temperature, T 0 , may change the rate of the wave generation, electron scattering, and acceleration, we compare simulations for three T 0 .Figure 7(a) shows that T 0 does not significantly affect the evolution of the electron distribution, and the normalized distribution f E (1000T ad , E)/f E (0, E) has nearly the same profile versus normalized energy E/T 0 .However, the wave magnetic field magnitude B w grows with the increase in the initial temperature T 0 (with the initial β(t = 0) ∝ T 0 ).Therefore, the initial temperature controls the wave intensity and the scattering rate, but this increase is compensated by the decrease in scattering rate with the electron energy, and the final evolution of the normalized energy distribution function f E (t, E/T 0 ) does not depend on T 0 .This effect shows some universality of the proposed self-consistent pumping mechanism of the electron acceleration in the presence of the adiabatic magnetic field compression.

Conclusion and Discussion
We have investigated a self-consistent electron acceleration process driven by an adiabatic periodic modulation of the electron distribution function.A natural mechanism of such modulation is the dynamics of ultralow-frequency magnetic field structures (e.g., magnetic holes or compressional waves, see Liu et al. 2021;Huang et al. 2022) that magnetically trap hot electron populations and transport them across the inhomogeneous and dynamical plasma medium near the Earth Figure 6.Evolution of the energy distribution function f E over the period of 1000 adiabatic oscillations.The initially Maxwellian distribution f transforms into a new distribution with a highly energetic tail (a).The efficiency of the energetic electron pumping can be estimated from the inserted panel showing the normalized energy distribution f E (t, E)/f E (0, E).The evolution of the electron pitch-angle distribution function f 0 a presents as a transport of the electron density into the area α 0 60°t hat isotropizes the system (b).The system parameters are ω pe /ω ce (t = 0) = 100, n h (t = 0)/n e = 10 −3 , and T 0 = T ∥ (t = 0) = 1 keV.The magnetic field magnitude is described by an equation B t B t T 0 1 sin ) with T ad = 200π/ω ce (t = 0).The integral wave amplitude B w is given by Equation (30) with approximation coefficients C = 1 and A = 0.75.
as a function of E/T 0 for several values of the initial temperature T 0 .Profiles of the distribution function f are nearly identical, within the uncertainties of numerical simulations.(b) Time evolution of the normalized wave amplitude B w /B 0 .The system parameters are: ω pe / ω ce (t = 0) = 100, n h (t = 0)/n e = 10 −3 .The magnetic field magnitude is described by an equation B t B t T 0 1 sin ) with T ad = 200π/ω ce (t = 0).The wave amplitude B w is given by Equation (30) with approximation coefficients C = 1 and A = 0.75.Note small amplitude variations of f E (1000T ad , E)/f E (0, E) are due to the numerical scheme error.
bow shock (e.g., Huang et al. 2017;Yao et al. 2020;Karlsson et al. 2022) and in the solar wind (e.g., Wang et al. 2021;Yu et al. 2021Yu et al. , 2022)).The modulation of the electron distribution periodically makes it unstable for the generation of highfrequency whistler-mode waves (e.g., Huang et al. 2018;Yao et al. 2021;Jiang et al. 2022), which should scatter electrons in the pitch-angle space.Indeed, observations of compressional magnetic structures near the bow shock are often associated with the whistler-mode wave activity (Hull et al. 2012(Hull et al. , 2020)).Actually, a similar type of compressional magnetic structure with trapped hot electrons generating whistler-mode waves is observed within planetary magnetospheres (e.g., Zhang et al. 2017;Shustov et al. 2019;Zhang et al. 2019;Guo et al. 2021), i.e., the proposed mechanism of the coupling of ultralowfrequency and high-frequency waves through the hot electron population is quite universal.The self-consistent electron scattering by whistler-mode waves breaks the reversibility of the electron adiabatic modulation and allows electrons to gain energy from ultralow-frequency magnetic field variations (see a detailed discussion of this mechanism in Lichko et al. 2017;Egedal et al. 2018).Such electron acceleration can be quite effective in increasing the phase space density of the energetic electron population, at least for the initial Maxwellian electron distribution.
Let us discuss the main simplifications of the proposed model of electron acceleration and further possible generalizations of this model: 1.For modeling the wave intensity and the wave spectrum, we use a semiempirical approach based on the utilization of the linear wave growth rate (Tao et al. 2017;Kuzichev et al. 2019).The parameters C and A that are being used have fixed values and do not change with the evolution of the system.In a more general case, a functional form of these coefficients has to be obtained: C = C(β, T ⊥ /T ∥ − 1) and A = A(β, T ⊥ /T ∥ − 1).Although this approach significantly simplifies the evaluation of electron distribution dynamics, it should be further verified or substituted with the full evaluation of the whistler-mode wave spectrum from the self-consistent quasi-linear equations (see Andronov & Trakhtengerts 1964;Kennel & Engelmann 1966;Galeev & Sagdeev 1979), showing if the proposed method is applicable for a certain set of system parameters.There are several successful examples of the application of the full set of quasi-linear diffusion equations (Yoon et al. 2012;Shaaban et al. 2019;Shaaban & Lazar 2020) and full kinetic simulations (Kuzichev et al. 2019;Roberg-Clark et al. 2019;Innocenti et al. 2019Innocenti et al. , 2020) ) for solar wind electrons, and the same equations can be implemented into the proposed model of the electron acceleration by the magnetic pumping.especially effective here could be the approach proposed in Yoon et al. (2017), where the full set of quasi-linear equations and a simplified form of the electron distribution function are combined to model the evaluation of wave spectrum dynamics.2. We have restricted our consideration to the electron generation and the scattering by whistler-mode waves, whereas the transversely anisotropic electron population may also generate electron cyclotron harmonics (e.g., Karpman et al. 1975;Ashour-Abdalla & Thorne 1977), which are quite effective in thermal electron scattering and acceleration (e.g., Horne et al. 2003;Ni et al. 2011).Moreover, electron resonant interaction with intense whistler-mode waves may generate localized structures of the electron phase space density, unstable to the generation of electrostatic waves (e.g., Agapitov et al. 2018;An et al. 2019) that are effective in thermal electron scattering and acceleration.The consideration of these effects requires self-consistent particle-in-cell simulations.3. We are focused on the generation of electron diffusion by low-intensity broadband whistler-mode waves, whereas observations near the Earth bow shock also detect quite intense monochromatic whistler-mode waves (Hull et al. 2012;Oka et al. 2017;Shi et al. 2023b condition, the form of the distribution function f stays the same throughout the adiabatic oscillations with respect to Equation (15).Therefore, the slope of the logarithmic distribution function f ln , d f dE ln , does not depend on the energy E. This applies so that the boundary condition at E = 0.1m e c 2 has to save the slope value, d f dE ln .For the pitch angle α, the boundary conditions have to satisfy the symmetry properties of the distribution function f.The solution is symmetrical with respect to α = π/2 and is invariant to the substitution α → − α.As a result, the following conditions are valid: change their signs during the evolution of the system, and thus, both energy boundary conditions have to be used to obtain a stable difference scheme.We introduce a half-step (τ/2) scheme that separates the energy E and the pitch-angle α contributions for F f ln : where indexes j, n, and m correspond to the time t, the energy E and the pitch-angle α grids, respectively, and τ, ε, and a are the corresponding step sizes; functions with a tilde are the half-step approximations taken at t + τ/2.  ) , where l = i exp wt ( ).The frequency ω is a complex variable, and its imaginary part defines the evolution of the wave amplitude.For a time-dependent operator U ˆ, the stability of the scheme can be achieved if |λ| 1 + Cτ, where C does not explicitly depend on any grid parameters or wave characteristics: where δj determines the time period of the wave amplitude growth, T is the time period of the stable solution, i.e., the last inequality determining the stability of the scheme condition should work on the time period T. For a plane wave, Equation (A2) can be written as where the constant C = C(r) is only dependent on the parameter r.As a result, the scheme is stable for any grid parameters.To solve Equation (A2), the system has to be simplified: . This system has to be solved using the tridiagonal matrix algorithm.

A.2. Diffusive Scattering Scheme
The scattering is characterized by the averaged diffusion operator D ˆ, and in the case of the function F f ln = , includes a nonlinear term.As the diffusion is weak, we consider it as a perturbation on a timescale of one adiabatic period.Thus, we describe the diffusion as a second step of the numerical scheme with the following evolution equation: ( ) The scheme appears to be as unconditionally stable as the one for the adiabatic evolution.

Figure 1 .
Figure 1.An overview of THEMIS (Angelopoulos 2008) observations of ultralow-frequency magnetic field fluctuations upstream of the Earth's bow shock (x GSE = 13R E ) in the top panel.Middle panels show the electron energy spectrum and the spectrum of high-frequency magnetic field fluctuations (white and black curves show the local electron gyrofrequency and the fraction of this gyrofrequency, respectively).The bottom set of panels zoom-in on the interval with whistlermode bursts embedded into local magnetic dips created by ultralow-frequency magnetic field fluctuations.We use magnetic field measurements from the fluxgate magnetometer (Auster et al. 2008) and the search-coil magnetometer (Le Contel et al. 2008).The electron spectrum is obtained from electrostatic analyzer measurements (McFadden et al. 2008).
) is valid.Expanding Equation (18) in a Taylor series over γ and taking the imaginary part of this equation, we obtain an equation for the growth rate γ: Considering the defined applicability interval, the expanded forms of J and J¢ can be substituted into Equation (21): left part of Equation (26) can be set to zero, and we obtain Equation (27) and the expanded forms of J and J¢, Equation (22) can be approximated as
(k)  is the wavenumber spectral density function, ε = B 2 /8π is the magnetic energy density.The wavenumber spectral density function S k (k) can be determined through the relation with the frequency spectral density function S Ω (Ω): 27) determines the relation between the resonant frequency Ω j and the resonant wavenumber k j : This shows that local diffusion coefficients D ij are proportional to different ν powers: D αp /D αα = D pα /D αα ∝ ν 2 and D pp /D αα ∝ ν 4 .As a result, the main contribution of lowfrequency waves, Ω 0 = ω ce , comes from the pitch-angle

Figure 4 .
Figure4.Dependency of the pitch-angle diffusion coefficient D αα ∝ ν 2 F(ν 2 ) and the resonant frequency Ω 0 on the parameter ν 2 .The effective part of the pitch-angle diffusion coefficient D αα gets closer to the zero-point as β increases.The smallness of Ω 0 /ω ce justifies the applicability of the linear approximation described by Equation (41).

Figure 5 .
Figure 5. 10 periods of adiabatic field variations: the evolution of the energy distribution function f E , the frequency spectral density function S Ω (Ω/ω ce ), and the electron anisotropy T ⊥ /T ∥ − 1.
Using the same method to analyze the stability of the scheme, we obtain that it is absolutely stable as 