On the Two Approaches to Incorporate Wave‐Particle Resonant Effects Into Global Test Particle Simulations

Energetic electron dynamics in the Earth's radiation belts and near‐Earth plasma sheet are controlled by multiple processes operating on very different time scales: from storm‐time magnetic field reconfiguration on a timescale of hours to individual resonant wave‐particle interactions on a timescale of milliseconds. The most advanced models for such dynamics either include test particle simulations in electromagnetic fields from global magnetospheric models, or those that solve the Fokker‐Plank equation for long‐term effects of wave‐particle resonant interactions. The most prospective method, however, would be to combine these two classes of models, to allow the inclusion of resonant electron scattering into simulations of electron motion in global magnetospheric fields. However, there are still significant outstanding challenges that remain regarding how to incorporate the long term effects of wave‐particle interactions in test‐particle simulations. In this paper, we describe in details two approaches that incorporate electron scattering in test particle simulations: stochastic differential equation (SDE) approach and the mapping technique. Both approaches assume that wave‐particle interactions can be described as a probabilistic process that changes electron energy, pitch‐angle, and thus modifies the test particle dynamics. To compare these approaches, we model electron resonant interactions with field‐aligned whistler‐mode waves in dipole magnetic fields. This comparison shows advantages of the mapping technique in simulating the nonlinear resonant effects, but also underlines that more significant computational resources are needed for this technique in comparison with the SDE approach. We further discuss applications of both approaches in improving existing models of energetic electron dynamics.

It is actually rather challenging to incorporate wave-particle resonant interactions into test particle simulations of plasma injections, because of their vastly different timescales.Figure 1 shows a schematic of the four main types of electron motions: cyclotron rotation with the timescale of τ g ∼ 10 3 s, bounce motion between magnetic mirrors with the time scale of τ b ∼ 10 1 s (for 100 keV electrons), earthward transport with the time-scale of τ D ∼ 10 3 s, and azimuthal motion around the Earth with the time-scale of τ ϕ ∼ 10 4 s (for 100 keV electrons).The ratio of the timescales of the fastest and the slowest motions can be τ ϕ /τ g ∼ 10 7 or τ D /τ g ∼ 10 6 .Individual waveparticle interaction, for example, via cyclotron resonance, occurs over ∼τ g , and thus such interactions may not be directly incorporated into test particle simulations of the plasma injections, which occur over ∼τ D .
The main limitation of the quasi-linear models is the requirement of a small wave intensity.A significant portion of most intense electromagnetic whistler-mode waves (Tyler et al., 2019;Wilson et al., 2011;Zhang et al., 2018Zhang et al., , 2019)), electrostatic whistler-mode waves (Agapitov et al., 2014(Agapitov et al., , 2015;;Cattell et al., 2008Cattell et al., , 2015;;Cully et al., 2008), and electromagnetic ion cyclotron waves (e.g., Tonoian et al., 2022;Wang et al., 2017) can exceed their threshold amplitudes and likely resonate with electrons nonlinearly.Such nonlinear resonant interactions Four main time-scales: electron gyroperiod τ g = mc/eB ∼ 10 3 s (with a typical magnetic field magnitude of B = 10 nT), electron bounce period τ b ∼ 4L CS /v 0 ∼ 10 1 s (with the current sheet thickness being L CS = 1R E and v 0 being the thermal velocity of 100 keV electrons), τ D ∼ L x /v x ∼ 10 3 s is electron transport time by plasma flows (v x ∼ 100 km/s), τ ϕ = 2πR/ v ϕ ∼ 10 4 s is the electron azimuthal drift period (the radial distance of the injection region is assumed to be R = 6R E and the electron energy is 100 keV).
include effects of phase bunching and phase trapping (see Albert, 1993;Inan & Bell, 1977;Itin et al., 2000;Karpman et al., 1975;Nunn, 1971;Solovev & Shkliar, 1986), which can significantly modify the characteristics of wave-particle interactions (see reviews by Albert et al., 2013;Artemyev, Neishtadt, Vainchtein, et al., 2018;Shklyar & Matsumoto, 2009, and references therein).Nonlinear effects can be incorporated into the kinetic equation (i.e., modified Fokker-Plank equation) for electron distribution functions (e.g., Artemyev, Neishtadt, et al., 2016;Artemyev, Neishtadt, Vasiliev, & Mourenas, 2018;Hsieh & Omura, 2017;Omura et al., 2015), but the corresponding characteristic equations will be different from those in the SDE approach.The main difference is the probability distribution function of pitch-angle/energy jumps due to resonant interactions: nonlinear effects cannot be described by Gaussian probability distributions that are often adopted to model the diffusive scattering within SDE approach.Recently, an alternative to SDE approach was proposed in (Artemyev, Neishtadt, Vainchtein, et al., 2018;Zheng et al., 2019), where a non-Gaussian probability distribution of pitch-angle/energy jumps has been incorporated into equations of motion for resonant electrons.This approach resembles the generalization of the classical mapping technique (e.g., Chirikov, 1979;Khazanov et al., 2014;Zaslavskii et al., 1989) to systems with a finite probability of very large pitch-angle/energy jumps due to phase trapping (see Artemyev et al., 2020).The mapping technique can reproduce many observed effects of nonlinear wave-particle interactions (e.g., Artemyev, Zhang, et al., 2022;Zhang et al., 2022) and in principle can be incorporated into testparticle codes (Artemyev, Neishtadt, & Angelopoulos, 2022).Therefore, both SDE and mapping technique can be used to include wave-particle resonant interactions into models of energetic electron dynamics in global electromagnetic fields provided by MHD/hybrid simulations.The mapping technique should generalize the SDE approach, but it is yet to be investigated whether the mapping equations can describe diffusive scattering and nonlinear resonant effects with the same accuracy level.The mapping technique usually adopts analytical equations to model pitch-angle/energy jumps, which do not include diffusive scattering (Artemyev et al., 2020).Recently Lukin et al. (2021) has generalized the mapping technique for the entire probability distribution function of pitch-angle/energy jumps, but it remains to be verified for electron cyclotron resonances with whistler-mode waves.
In this paper, we combine two approaches from Artemyev et al. (2020) and Lukin et al. (2021) to construct the mapping technique that operates with the probability distribution function of pitch-angle/energy jumps for electron cyclotron resonances with whistler-mode waves.We also compare results from this newly developed mapping technique to those from the SDE approach, which operates by a single characteristic of such probability distribution functions of pitch-angle/energy jumps-distribution variance, which dictates the diffusion rate.We examine two wave intensities: small intensity for the diffusive interaction, when SDE and mapping are expected to provide the same results, and large intensity with nonlinear wave-particle interactions, when the mapping technique is validated by full test particle simulations.Therefore, this paper demonstrates the approaches to incorporate the long-term effects of quasi-linear and nonlinear wave-particle interactions into global scale test-particle models.
The rest of the paper starts with introducing basic characteristics of electron resonant interactions with whistlermode waves in the dipole magnetic field (see Section 2).Then in Section 3, we introduce the main property of electron ensemble dynamics-the probability distribution function of energy/pitch-angle jumps in a single resonance.This distribution is used to introduce SDE approach (in Section 4) and mapping technique (in Section 5).Results obtained from SDE and mapping approaches are compared with test particle simulations in Section 6.Finally, we discuss the applicability of both methods and briefly summarize our conclusions in Section 7.

Resonant Electron Interactions With Whistler-Mode Waves
We examine the interaction of relativistic electrons (rest mass m e , charge e, speed of light c) with field-aligned whistler-mode waves moving in the dipole magnetic field B 0 (λ), where λ is the magnetic latitude.Electron dynamics in such a system can be described by the following Hamiltonian (Albert, 1993;Vainchtein et al., 2018): where (z, p z ) and (ψ, I x ) are two pairs of conjugate variables, the field-aligned coordinate and moment, and electron gyrophase and magnetic moment.The electron cyclotron frequency Ω 0 (λ) is defined as where Ω eq = eB 0 (0)/m e c is the equatorial cyclotron frequency determined by the radial distance from the Earth to the equatorial crossing of the magnetic field line, that is, L-shell.In the dipole field, the magnetic latitude, λ, is related to the field-aligned coordinate z as: where R E is the Earth radius.Note that in Hamiltonian (Equation 1), the magnetic moment, I x , conjugate to the gyrophase, ψ, should be normalized in such a way that I x Ω 0 has a dimension of energy.
The second term in Hamiltonian (Equation 1) describes the wave contribution to the electron dynamics.The wave amplitude B w is modeled as: where ɛ parameter controls B w /B 0 (0), function f(λ) determines the wave latitudinal profile that agrees with the wave generation around the equator, that is, wave amplitude growth within the generation/amplification region Δλ, which is controlled by the value of c λ : Δλ ∼ 1/c λ , followed by saturation (see typical B w (λ) profiles from statistical models in Agapitov et al., 2013Agapitov et al., , 2018)).Note that we assume the waves only exist in one hemisphere, z > 0, because the wave field and wave propagation equations are symmetric in two hemispheres (z → z does not change the system equation).
This assumption allows us to simplify the calculations, because in this case the wave-particle interaction occurs only during a quarter of the bounce period when the particles move northward from the equator.
The wave phase is given by equation φ = k(λ) ż ω, where the wave frequency ω is constant and the wave number k(λ) is determined by the cold plasma dispersion relation (Stix, 1962): (5) The plasma frequency Ω pe is given by the empirical function (Denton et al., 2006) with the equatorial values Ω pe,eq given by the empirical function (Sheeley et al., 2001) Ω pe,eq /Ω 0 (0) ≈ L.
Wave phase ϕ linearly depends on time, ∂ϕ/∂t = ω = const, and thus Hamiltonian (Equation 1) has an integral of motion h (see, e.g., review by Shklyar & Matsumoto, 2009): In the absence of wave perturbation, B w = 0, particle energy m e c 2 (γ 1) and equatorial pitch-angle α eq = arc sin 2I x Ω 0 (0)/ m e c 2 γ 2 1)) are conserved; here is the gamma factor.Wave-particle resonant interactions can change particle's energy and equatorial pitch-angle (change I x ), but due to the conservation of h(γ, I x ) waves move electrons along a specific curve in the energy, pitch-angle space.Therefore, for a fixed value of h (i.e., when all particles have the same initial h), the waveparticle resonant interaction becomes a 1D problem, and we may just examine energy changes, whereas changes of equatorial pitch-angle (or I x ) can be determined from h = const.
To obtain basic characteristics of wave-particle resonant interactions, we integrate electron equations of motion (Hamiltonian equations) with fourth order Runge-Kutta scheme with an adaptive time step (1/20 of the local electron gyroperiod).Throughout the rest of the paper, we use the following dimensionless variables: Thus, the dimensionless Hamiltonian and integral of motion h take the following forms where Ω(λ) → Ω(λ)Ω 0 (0).To demonstrate the result, we use the following parameters throughout: h = 3/2, λ 0 = 5°, and c λ = 180/π.Each integration starts from the equatorial plane (λ = z = 0) and stops there after N res resonant wave-particle interactions.Since waves only exist in the northern hemisphere at a fixed frequency, and propagate along magnetic field lines (i.e., resonate with electrons through the first-order cyclotron resonance only), each electron bounce period corresponds to one wave-particle interaction on a time scale of a quarter bounce period.Effects of multiple resonances within one bounce period (e.g., due to oblique wave propagation or wave frequency variation with time, see Artemyev et al., 2021;Hsieh & Omura, 2023;Shklyar & Matsumoto, 2009) can be incorporated into this approach by including additional wave terms into Equation 9.
Depending on the wave magnitude, there are two possible regimes of wave-particle resonant interactions: diffusive scattering and nonlinear resonant interactions.Figure 2 shows the profiles of electron energy evolution with time for both regimes (all electrons have the same initial energy and pitch-angle, but random initial gyrophases).For small wave amplitudes: (a) the interaction is linear and energy evolution with time is stochastic.After each wave-particle interaction, the electron energy undergoes a small (compared to electron initial energy) positive or negative jump with almost equal probabilities to increase or decrease the energy.This process (with normalized diffusion coefficients) can be approximated by the Wiener stochastic process, that is, the evolution of electron distribution function can be described by the Fokker-Planck equation or, equivalently, by SDEs.For large wave amplitudes: (b) electron dynamics cannot be described with the diffusive approach.Most of the time, particle energy undergoes small negative jumps (but positive jumps are also possible, see Albert et al., 2022) and there is a nonzero probability of large positive jumps, caused by particle phase trapping.

Probability Distributions of Energy Jumps
To quantify variations of electron energy due to wave-particle resonant interactions, we use test particle simulations and evaluate the distributions of energy jumps as a function of initial energy E 0 , ΔE(E 0 ). Figure 3 shows such distributions for several initial electron energies and two different magnitudes of the wave field.For each histogram, we integrated trajectories of N p = 32,768 test particles with random initial gyrophases and the same initial energy and h values.Each integration includes one bounce period (i.e., a single resonant interaction).For small wave amplitudes (panels (a-d)), electron energy jumps are small and randomly (but symmetrically) distributed around zero.For sufficiently high wave amplitudes (panels (e-h)), there appear nonlinear resonant effects: most of the electrons lose their energy due to the phase bunching with ΔE < 0, while a small population of phase trapped particles (panels (g-h)) gain a significant portion of energy with ΔE > 0 comparable to their initial energy, E 0 .Changes of the electron pitch-angle (and I x ) are directly determined by the energy changes due to the conservation of the integral of motion (Equation 7).Therefore, the probability distribution of ΔE(E 0 ) fully characterizes the evolution of electron distribution function and can be used as a basic input for the SDE approach or mapping technique.Note that for intense, but very low-coherent waves, the probability distribution function of ΔE(E 0 ) will be symmetric relative to ΔE(E 0 ) = 0, and thus will largely resemble distributions from panels (a-d) (see An et al., 2022;Frantsuzov et al., 2023;Gan et al., 2022;Zhang et al., 2020).In panels (e, f), we show simulation results for E 0 without trapping: a small population of ΔE > 0 is due to the positive phase bunching effect (see Albert et al., 2022).In panels (g, h), we show simulation results with phase trapping; the inserted panels show the expanded view of the population with large ΔE > 0.

Stochastic Differential Equations
In the limit of small wave field amplitudes, the wave-particle resonant interaction is diffusive (e.g., Allanson et al., 2022;Kennel & Engelmann, 1966;Lyons & Williams, 1984;Schulz & Lanzerotti, 1974, and references therein).Note that in inhomogeneous magnetic field, these interactions are diffusive even for monochromatic waves (see Albert, 2010;Shklyar, 2021).Therefore, for such systems one can use the Fokker-Planck equation to describe the evolution of the electron distribution function f(χ) of the velocity vector χ: where t is time, μ(χ(t), t) is an N-dimension vector of the drift coefficient, D is an N × N matrix of diffusion coefficients (Albert, 2018;Lyons & Williams, 1984;Schulz & Lanzerotti, 1974).Instead of solving the Fokker-Planck equation, we can use the corresponding Ito SDEs to integrate trajectories of quasi-particles (Tao et al., 2008;Zheng et al., 2014): where Δt is the time step over which we calculate the change of χ, σ(χ(t), t) is an N × N-dimension matrix related to the diffusion coefficients written in such a way that D = 1 2 σσ T , W t is an N-dimension standard Wiener process; where N is a vector of standard normal random values, N i ∼ N(0, 1).
The term quasi-particles means that we do not directly integrate the equations of motion, but treat the change of χ as a stochastic process and approximate it by Equation 12.As a result, two quasi-particles having equal initial conditions χ 0 may have different trajectories χ(t) (in the numerical integration of Equation 12, one can fix the seed of the pseudo random generator to preserve the sequence of random numbers and make the results repeatable).We will examine electron distributions in the energy space, and thus the Ito Equation 12 can be rewritten in the following form This equation describes the energy evolution for fixed h given by Equation 10, that is, for a monochromatic wave we may reduce the energy, pitch-angle evolution to the energy-only evolution and calculate the associated pitchangle changes from h conservation.For generic wave spectra, there are three diffusion rates (energy, pitch-angle, and mixed energy-pitch-angle), and thus we would need to solve a system of equations for the energy and pitchangle evolution (see detail in Tao et al., 2008).Note that if we rewrite the Fokker-Planck equation in terms of energy, additional coefficients of variable transformation from velocity (momentum) to energy and pitch-angle (Lamé coefficients) should be added (e.g., Glauert & Horne, 2005), but the Ito equation will still have the same form.
The main limitation of this approach is the requirement on small energy changes ΔE for each wave-particle interaction, that is, it is only applicable for systems without nonlinear effects of phase trapping.Equation 13includes the drift μ E (E 0 ) and bounce averaged diffusion D E (E 0 ) coefficients, which can be estimated using the distribution of energy jumps ΔE(E 0 ): For the drift term, we use which keeps the divergence-free form of the Fokker-Planck equation (Allanson et al., 2022;Lemons, 2012;Lichtenberg & Lieberman, 1983;Sinitsyn et al., 2011;Zheng et al., 2019).Note that using such estimates of drift and diffusion coefficients, we implicitly assume that the ΔE(E 0 ) distributions are symmetric relative to the mean ΔE value.Figures 4a and 4b show the numerically obtained ΔE distributions (black) and their fittings to symmetric Gaussian distributions (this fitting outputs the variance ∼D E (E 0 )).However, this assumption may be violated even in the case of small wave amplitudes (see Figure 3c), when we can over-or underestimate the drift and diffusion coefficients.We will discuss the corresponding uncertainties in Section 6. Figure 4c shows the drift and diffusion coefficients as a function of the electron initial energy E 0 .
Each electron bounce period corresponds to a single wave-particle interaction (as commented above), and thus the time step of integration Δt in Equation 13 should be set equal to τ b (E 0 ).The main advantage of the SDE approach (and also of the mapping technique, which will be discussed in Section 5), in comparison with the test particle approach, is the significant reduction of computational time for long-term simulations (see also discussions in Lukin et al., 2021): we integrate trajectories of quasi-particles with a time step equal to the bounce period, considering only the effects of wave-particle interactions and do not trace particles during adiabatic paths of their motion.So at each step of integration, we recalculate the electron energy as where N i ∼ N(0, 1) is the standard normal random number.
Panels (d, e) in Figure 4 show examples of the electron energy profiles calculated by direct integration of the Hamiltonian equations or using the SDE approach.Particles having equal initial conditions will have different trajectories from these two approaches, due to randomness of resonant interactions, but statistical properties of the evolution of the electron ensemble should be the same.We have verified this property by a set of simulated evolution of the electron ensembles (not shown).

Mapping Technique
For high wave amplitudes, the SDE approach is no longer applicable, because of the nonlinear effects of waveparticle interactions.In this case, the Wiener stochastic process cannot describe the evolution of electron energy as there is a finite probability of large energy jumps caused by the phase trapping.For such systems, the Ito SDE can be generalized by introducing a new stochastic process accounting for the phase trapping and phase bunching.This process can be represented by a series of mapping functions depending on electron energy (the changes of pitch-angle are related to energy changes through the integral of motion (Equation 7)).
Figure 5a shows the distribution of energy jumps for the wave amplitude ɛ = 10 3 and initial electron energy E 0 = 257 keV.For this particular energy, there is no electron phase trapping: the resonant latitude for this energy and h (which determines the equatorial pitch-angle) corresponds to zero probability of electron trapping, because wave intensity increases slower along the resonant trajectory than the background magnetic field inhomogeneity (see detailed equations determining the trapping probability in the Appendix of Vainchtein et al., 2018, and in references therein).Hence the procedure to construct the mapping function is straightforward: using this distribution, one can calculate the cumulative function of energy jumps F(ΔE, E 0 ) ∈ [0, 1] (right vertical axis in Figure 5a) and then use the inverse function ΔE(U) = F 1 (U, E 0 ), U ∈ [0, 1] as the generator of electron energy jumps for each wave-particle interaction.Figure 5b shows the distribution of energy jumps for an initial electron energy of E 0 = 350 keV.There is a small probability of electron phase trapping in this case, thus we can divide the ΔE distribution into bunching (|ΔE|/ E 0 ≪ 1) and trapping (ΔE/E 0 ∼ 1) parts, and then calculate corresponding cumulative functions: F b (ΔE, E 0 ) and F t (ΔE, E 0 ).Note that this separation of ΔE distribution into two parts is not necessary, but it simplifies the calculations by allowing a linear interpolation for the cumulative function.There is a gap of ΔE between two parts of ΔE distribution, that is, for the case shown in Figure 5b the electron energy cannot change by the value between ∼4 and ∼95 keV (prohibited values of ΔE in this case).If we would not separate the ΔE distribution into two parts, this ΔE gap will require a separate treatment.Using trapping ΔE distributions, we may estimate the probability of electron trapping p trap (E 0 ) (shown in Figure 5c) as a function of the initial electron energy: where N p is the total number of particles used to construct the ΔE(E 0 ) distribution and N(ΔE(E 0 )/E 0 ∼ 1) is the number of trapped particles.
We use the subscript b for the bunching part of the ΔE(E 0 ) distribution and also for the entire ΔE(E 0 ) distribution if the probability of trapping is equal to zero.Then the generalization of the Ito SDE is straightforward: we can replace the Wiener process in Equation 13with the constructed mapping functions.At each integration step, we generate two uniform random numbers u 1 , u 2 ∈ U(0, 1) and recalculate electron energy as Note that the gap between ΔE distributions due to phase trapping and bunching (see Figure 5b) depends on system parameters (mostly on the wave field model) and, for example, for the system with small wave-packets, the energy of phase-trapped electrons could not be far away from the initial energy (e.g., Omura et al., 2015).In this case, there is no need to divide the distributions into bunching and trapping parts and thus only one cumulative function, F(ΔE, E 0 ), needs to be constructed.In this case, the computation scheme slightly simplifies: at each iteration one needs to calculate a single, uniform random number u ∈ U(0, 1) and the energy change can be calculated in the same manner: Figures 5d and 5e show several profiles of electron energy as a function of the number of wave-particle resonances, calculated by integration of the test particle trajectories and using the mapping technique, respectively.Both approaches show very similar electron energy dynamics: long-duration electron drift to smaller energy due to the phase bunching and rare large jumps to higher energy due to the phase trapping.Note that for both SDE and mapping techniques, we still need to use test particles to calculate the initial distributions of ΔE(E 0 ), but having these distributions (requiring statistics of single wave-particle interactions), we can use the simplified integration scheme with a time step equal to the bounce period.

Method Validation
For small wave amplitudes, we can apply all three methods (test particles, SDE, and mapping) to simulate the evolution of the electron distribution function, and these three methods are expected to give the same results.The test particles approach is expected to be more precise, because it does not rely on the constructed ΔE distribution and is based on the full set of equations of motion.The main advantage of SDE and mapping techniques is their computational efficiency in long-term simulations, so they can be less accurate, but should still describe the main features of the evolution of electron distributions and their results should statistically repeat those from the test particles approach.Figure 6 (left panels) shows the evolution of the electron distribution function for a small wave amplitude (ɛ = 10 4 ) at four time instants.Hereinafter N p = 32,768 test particles are used to compute the changes of the distribution functions for all mentioned approaches.Without nonlinear resonant effects, both SDE and mapping technique show results that are consistent with the test particle simulation.In this case, the evolution is diffusive (as expected for quasi-linear theory) and shows a spread of the initially localized electron phase space density peak.The difference between SDE and the test particle simulation, most clearly seen around E ∼ 420 keV, is due to the overestimation of the diffusion coefficients.We evaluate diffusion coefficients as a half of the variance of ΔE distributions, and thus we assume that ΔE distributions are symmetric relative to the mean value.However, even in the case of low-amplitude waves (see, e.g., Figure 3c) this assumption may not work, which will result in an overestimation of the diffusion rate.The mapping technique does not require any assumptions about ΔE distributions, and thus it performs better even in the case of low-amplitude waves.
Figure 6 (right panels) shows the evolution of the electron distribution function for a large wave amplitude (ɛ = 10 3 ) at four time instants.For such intense waves, the SDE approach becomes inapplicable, but we can still compare results of test particle simulations and the mapping technique.The mapping technique accounts for nonlinear resonant effects (e.g., phase trapping) and describes well the evolution of electron distributions.After several wave-particle resonant interactions (see top right panel of Figure 6), the main electron population propagates to lower energies due to the phase bunching, while a small population becomes trapped by waves and gains energy.During the drift of the main population to the smaller energies, the probability of particle trapping increases (see Figure 5c) and more particles become trapped and accelerated.Accelerated particles appear in the resonant latitudes where no more phase trapping is possible, and thus these particles start losing their energy due to the phase bunching.Around the time when the main population (at the initial peak of electron phase space Gray color shows the initial distribution, black color shows test particle simulation results, red color shows results obtained with stochastic differential equation, and blue color shows results for the mapping technique.Time (in seconds) is calculated under the assumption of a single resonant interaction per bounce period, τ b (E), where τ b is evaluated at L-shell = 6 and for equatorial pitch-angles derived from Equation 7.
density) reaches the left boundary of the allowed energies, the processes of phase bunching and phase trapping statistically compensate each other.This results in formation of a plateau in the distribution function (see the bottom right panel in Figure 6).Such an evolution of the electron distribution is consistent with theoretical predictions for the system with multiple nonlinear resonances (Artemyev et al., 2019).
The main uncertainty of the mapping technique arises at the edges of the simulation domain due to the finite grid size (discretization) in the ΔE distribution.In this approach, we treat the electron energy change as a probabilistic process and in Equations 13 and 17, electrons can reach the energy E < E min outside of the simulation domain.In the subsequent calculations, therefore, we will use the mapping functions corresponding to the energy E min to calculate electron energy change.This underestimates the probability of positive and overestimate probability of negative energy jumps, so that particles can stay at the edges of the energy grid for a long time.The SDE approach also suffers from this boundary effect.This problem can be eliminated by introducing boundary conditions, for example, reflective boundaries, or by normalizing the ΔE distribution around boundaries (see discussions in the Appendix of Artemyev et al., 2021).

Discussion and Conclusions
In this study we compare two approaches that incorporate the wave-particle resonant effects into test particle simulations: SDE approach and mapping technique.Both approaches simplify the characterization of waveparticle interactions and reduce resonant effects to a ΔE distribution of energy jumps during a single interaction; in particular, SDE further simplifies this description and operates only by the second moment (variance) of this distribution.The comparison shows that these two approaches provide the same results for the system with low-amplitude waves, whereas electron nonlinear resonant interactions with intense waves may be only described well by the mapping technique.Note that the diffusion approximation is not only applicable to low intensity waves, but also to systems with very intense, low-coherent waves where the resonance overlapping results in destruction of nonlinear effects (An et al., 2022;Frantsuzov et al., 2023;Gan et al., 2022;Tao et al., 2013;Zhang et al., 2020).Let us discuss advantages and limitations of this technique.
The mapping technique is based on ΔE distributions, which for fixed system parameters depends only on the initial electron energy E 0 , that is, we deal with 2D distributions of energy (or pitch-angle) jumps, F(ΔE, E 0 ), which are used as described in Section 5.These distributions can be determined with any required accuracy from a set of short-term test particle simulations.However, in realistic space plasma systems, we deal with some ensemble of waves that can be described by the distribution of wave amplitudes and frequencies, P w (B w ,ω) ; for reasonable discretization levels and available statistical wave data sets, this wave distribution usually consists of ∼10 2 -10 3 different pairs of (B w , ω) (see, e.g., Artemyev, Zhang, et al., 2022;Zhang et al., 2022).Therefore, for the mapping technique, we would need to evaluate 10 2 -10 3 realizations of 2D F(ΔE, E 0 ) distributions.Additional system dimensions can be introduced due to the dependence of P w (B w ,ω) on the geomagnetic activity and geophysical coordinates (MLT, L-shell).Therefore, although the mapping technique essentially reduces computations relative to the full test particle simulations, this approach requires significant resources for simulation of electron dynamics during long-term global events (where wave and background characteristics can vary significantly), for example, geomagnetic storms.More suitable application of the mapping technique is simulations of localized (spatially and temporally) events, like plasma sheet injections (Artemyev, Neishtadt, & Angelopoulos, 2022) or strong precipitation bursts (Artemyev, Zhang, et al., 2022;Zhang et al., 2022).For global simulations with SDE, which substitute 2D F(ΔE, E 0 ) distributions by 1D diffusion coefficients D EE (E 0 ), should be much more realistic, although this approach does not account for nonlinear resonant effects.
Both the SDE approach and mapping technique assume the evaluation of F(ΔE, E 0 ) distributions before simulating the electron dynamics, and thus such distributions are usually evaluated for a prescribed background magnetic field.This simplification may work for the inner magnetosphere, where the background dipole field does not vary too much (see Ni et al., 2011;Orlova & Shprits, 2010, for discussions on when this assumption does not work well), but cannot be well justified for plasma injections characterized by rapid, significant variations of the magnetic field configuration (see discussion in Ashour-Abdalla et al., 2013;Birn et al., 2022;Sorathia et al., 2018).The background magnetic field configuration determines the resonant condition and thus should control the efficiency of electron scattering by waves.In principle, the effect of the magnetic field reconfiguration can be included into SDE and mapping approaches, but this would require the evaluation of F(ΔE, E 0 ) distributions for multiple magnetic field configurations, which is usually unlikely with available computational resources.A possible solution will be to analytically evaluate F(ΔE, E 0 ) distributions (e.g., Artemyev et al., 2021;Vainchtein et al., 2018), but this solution requires more developed theoretical models of wave-particle resonant interactions including effects of wave-field deviations from a simple plane wave model (see discussions in Artemyev et al., 2023;Mourenas et al., 2018).So far such effects have been evaluated (e.g., Allanson et al., 2021;An et al., 2022;Gan et al., 2022;Tao et al., 2011;Zhang et al., 2020) and incorporated into F(ΔE, E 0 ) distributions (e.g., Hsieh et al., 2020Hsieh et al., , 2022;;Kubota & Omura, 2018) only via numerical integration of a large test particle ensemble.
Both the SDE approach and mapping technique are based on numerical evaluations of the probability distribution function of energy jumps, F(ΔE, E 0 ).Although we have adopted the plane wave approximation in this study, the procedure of F(ΔE, E 0 ) evaluations is independent of such approximations and F can be calculated for more realistic modes of wave-packets (see observations and simulations in Zhang et al., 2018Zhang et al., , 2021)).This generalization for the finite wave-packet size is quite important for constraining the efficiency of nonlinear resonant effects (see, e.g., An et al., 2022;Kubota & Omura, 2018;Zhang et al., 2020).Therefore, in the next step, one would need to incorporate realistic distributions of whistler-mode wave-packet sizes.

Figure 1 .
Figure 1.Schematic of four main time-scales of electron motion during plasma sheet injections into the inner magnetosphere.Four main time-scales: electron gyroperiod τ g = mc/eB ∼ 10 3 s (with a typical magnetic field magnitude of B = 10 nT), electron bounce period τ b ∼ 4L CS /v 0 ∼ 10 1 s (with the current sheet thickness being L CS = 1R E and v 0 being the thermal velocity of 100 keV electrons), τ D ∼ L x /v x ∼ 10 3 s is electron transport time by plasma flows (v x ∼ 100 km/s), τ ϕ = 2πR/ v ϕ ∼ 10 4 s is the electron azimuthal drift period (the radial distance of the injection region is assumed to be R = 6R E and the electron energy is 100 keV).

Figure 2 .
Figure2.Examples of electron energy dynamics for the system with diffusive scattering (a) and nonlinear resonant effects (b).All electrons have the same initial energy and pitch-angle, but random initial gyrophases.Energy is plotted versus number of resonant interactions (i.e., number of electron bounce periods).

Figure 3 .
Figure3.Examples of ΔE probability distributions for systems without nonlinear resonant effects (a-d) and with nonlinear resonant effects (e-h).In panels (e, f), we show simulation results for E 0 without trapping: a small population of ΔE > 0 is due to the positive phase bunching effect (seeAlbert et al., 2022).In panels (g, h), we show simulation results with phase trapping; the inserted panels show the expanded view of the population with large ΔE > 0.

Figure 4 .
Figure 4. Panels (a, b) show ΔE distributions (in black), obtained after a single wave-particle interaction, and Gaussian distributions with the same variance (in red).Panel (c) shows the energy diffusion rate and drift term versus initial energy.Panels (d, e) show examples of 10 trajectories of long-term electron energy dynamics (∼300 resonances) for a system without nonlinear resonant effects.

Figure 5 .
Figure 5. Panels (a, b) show ΔE distributions and the corresponding cumulative probability distribution function (only on panel (a)).Panel (c) shows the probability of electron phase trapping versus the initial energy.Panels (d, e) show examples of 10 trajectories of long-term electron energy dynamics (∼300 resonances) for a system with nonlinear resonant effects.

Figure 6 .
Figure6.Evolution of the electron energy distribution for a systems without nonlinear resonant effects (left panels) and with nonlinear resonant effects (right panels).Gray color shows the initial distribution, black color shows test particle simulation results, red color shows results obtained with stochastic differential equation, and blue color shows results for the mapping technique.Time (in seconds) is calculated under the assumption of a single resonant interaction per bounce period, τ b (E), where τ b is evaluated at L-shell = 6 and for equatorial pitch-angles derived from Equation7.