Nonlinear resonant interactions of radiation belt electrons with intense whistler-mode waves.

A. V. Artemyev<fnref n="1" /> · D. Mourenas<fnref n="2" />,<fnref n="3" /> ·

X.-J. Zhang<fnref n="4" />,<fnref n="1" /> · O. Agapitov<fnref n="5" /> · A. I.

Neishtadt<fnref n="6" /> · D. L. Vainchtein<fnref n="7" /> · A. A.

Vasiliev<fnref n="8" /> · X. Zhang<fnref n="1" /> · Q. Ma<fnref n="1" /> · J.

Bortnik<fnref n="1" />, V. V. Krasnoselskikh<fnref n="9" />

Received: date / Accepted: date

**Abstract** The dynamics of the Earth's outer radiation belt, filled by energetic electron fluxes, is largely controlled by electron resonant interactions with electromagnetic whistler-mode waves. The most coherent and intense waves resonantly interact with electrons nonlinearly, and the observable effects of such nonlinear interactions cannot be described within the frame of classical quasi-linear models. This paper provides an overview of the current stage of the theory of nonlinear resonant interactions and discusses different possible approaches for incorporating these nonlinear interactions into global radiation belt simulations. We focused on observational properties of whistler-mode waves and theoretical aspects of electron nonlinear resonant interactions between such waves and energetic electrons.

**Keywords** Relativistic electron precipitation · Radiation Belts · Magnetosphere · Electromagnetic Ion Cyclotron Waves · Whistler-mode chorus · Plasma waves

**Mathematics Subject Classification (2020)** MSC code1 · MSC code2 · more

# 1 Introduction

The outer radiation belt is a near-Earth magnetospheric region filled with energetic electrons reaching relativistic and even ultra-relativistic energies [Allison and Shprits, 2020; Baker et al., 2014, 2016; Horne, 2007; Horne et al., 2005].

<footnote n="1">1 University of California, Los Angeles, Los Angeles, CA 90095, USA</footnote>

<footnote>E-mail: aartemyev@igpp.ucla.edu,</footnote>

<footnote n="2">2 CEA, DAM, DIF, Arpajon, France,</footnote>

<footnote n="3">3 Laboratoire Matière en Conditions Extrêmes, Université Paris-Saclay, CEA, Bruyères-le-Châtel, France</footnote>

<footnote n="4">4 University of Texas at Dallas, Richardson, TX 75080,</footnote>

<footnote n="5">5 Space Sciences Laboratory, University of California, 94720, Berkeley, USA,</footnote>

<footnote n="6">6 Department of Mathematical Sciences, Loughborough University, Loughborough, LE11 3TU, United Kingdom,</footnote>

<footnote n="7">7 Nyheim Plasma Institute, Drexel University, Camden, New Jersey 08103, USA,</footnote>

<footnote n="8">8 Harbour University, Carrer de Rosa Sensat 9-11, 08005 Barcelona Spain</footnote>

<footnote n="9">9 LPC2E, CNRS-University of Orléans-CNES, 45071, Orléans, France</footnote>
<figure><img src="image_1.png" /><figcaption>**Fig. 1** Schematic view of electron motion in dipole magnetic field. Three types of periodical motions are shown (from left to right): azimuthal drift around the Earth with the invariant $\Phi = \oint P_D dr_D$ where $P_D = eA_\Phi/c + mv_D \approx eA_\Phi/c$ is the azimuthal momentum and $dr_D = Rd\Phi$ is the element of azimuthal trajectory; bounce motion along magnetic field lines between mirror points with the invariant $J_\| = \oint p_\| ds$ where $s$ is the field aligned coordinate; gyrorotation around the magnetic field with the invariant $\mu = \oint p_\perp dr_\perp$ where $(p_\|, p_\perp)$ are parallel and perpendicular momentum components. The direction of the magnetic field is shown in blue. In dipole field $R = R_E L$ ($R_E$ is the Earth radius and $L$ is the L-shell parameter), the magnetic field inhomogeneity scale along magnetic field lines is the curvature radius $R_C = R/3$, the electron gyroradius $\rho \ll R$.</figcaption></figure>

The potentially significant damaging effects of relativistic electron fluxes for the many satellites on orbit continuously drives investigation, modelling, and forecasting of the radiation belt dynamics [Baker et al., 2018; Horne, Glauert, Meredith, Boscher, Maget, Heynderickx and Pitchford, 2013]. Although the magnetic field in the outer radiation belt is dominated by the quite stable dipole field of the Earth, energetic electron fluxes in this region may vary by several orders of magnitude. Wave-particle resonant interaction is the main driver of such flux variations: the radial drift of energetic electrons is provided by drift resonance with ultra-low-frequency (ULF) waves, whereas bounce and cyclotron resonances with extremely and very-low-frequency (ELF/VLF) waves are responsible for electron pitch-angle scattering and energization [Lyons and Williams, 1984; Schulz and Lanzerotti, 1974; Trakhtengerts and Rycroft, 2008; Tverskoy, 1969]. Without wave-particle resonant interactions, since electrons are magnetized by the strong dipolar geomagnetic field, three adiabatic invariants are conserved during electron motion, corresponding to three types of electron periodical motions: the magnetic moment, $\mu$, corresponds to the gyrorotation, the bounce invariant, $J_\|$, corresponds to bounce oscillations along field lines, and the third invariant, $\Phi$, corresponds to the azimuthal drift motion around the Earth (see schematic 1 and [Schulz and Lanzerotti, 1974; Tverskoy, 1969]). The conservation of these invariants can fix the electron energy $\gamma = E/m_e c^2$ and equatorial pitch-angle $\alpha_{eq}$, and the L-shell, the normalized distance between the Earth center and the farthest (equatorial) point of the magnetic field line along which electrons are bouncing at every point along the electron drift orbit. Therefore, any change of electron phase space density should be attributed to destruction of one or more of these invariants.

Adiabatic invariants are conserved exponentially well [Cohen et al., 1978; Dykhne, 1960; Kulsrud, 1957; Lenard, 1959; Neishtadt, 2000; Slutskin, 1964], and their de-
<figure><img src="image_2.png" /><figcaption>**Fig. 2** Example of electron resonant interactions with a field-aligned whistler-mode wave: effect of electron scattering and magnetic moment violation is shown. Plots show the electron magnetic moment $I_x = m_e c^2 (\gamma^2 - 1) \sin^2 \alpha / \Omega_0$ along the particle trajectory during half of an electron bounce period in a dipole magnetic field. System parameters are: whistler-mode wave frequency equal to 0.35 of equatorial electron gyrofrequency, wave magnetic field amplitude is 10 pT, the background magnetic field is dipolar [Bell, 1984] with L-shell equal to 6, equatorial plasma density is given by an empirical model [Sheeley et al., 2001] and it is constant along magnetic field lines. The magnetic moment is $\mu = eI_x/m_e c$, and time is normalized to $c/R$ with $R = LR_E$, with $R_E$ the Earth radius. The total time of this simulation is $\sim 41$s.</figcaption></figure>

struction requires the action of a force varying in space or time faster than the
spatial/temporal scale of the periodic motion corresponding to the specific invari-
ant [e.g., Roberts, 1969]. For example, destruction of $\Phi$ requires some external
force with a temporal scale comparable to the electron azimuthal drift period,
whereas violation of $\mu$ requires an external force with a temporal scale comparable
to the electron gyroperiod. Such external forces can be Lorentz force of wave elec-
tromagnetic fields varying with the corresponding temporal scales. Figure 2 shows
an example of such invariant violation for the electron resonant interaction with
circularly polarized electromagnetic waves. Each resonant interaction lasts for a
short interval during which the adiabatic invariant experiences a random jump,
whereas the longer time intervals between resonant interactions are characterized
by conservation of the adiabatic invariant, when the particle is bouncing or drifting
far from the wave region.

Therefore, the primary theoretical problem for evaluation of the radiation belt
dynamics is to develop an approach describing the evolution of the electron phase
space density *f* due to multiple resonances with realistic electromagnetic waves.
Let us start our introduction to this problem with a brief description of a well devel-
oped, and most frequently used, theoretical concept of such evolution – quasi-linear
theory [Drummond and Pines, 1962; Kennel and Engelmann, 1966; Kennel and
Petschek, 1966; MacDonald and Walt, 1961; Roberts, 1969; Trakhtengerts, 1963;
Vedenov et al., 1962; Wentworth, 1963]. This theory describes the self-consistent
dynamics of the charged particle distribution and of the spectrum of electromag-
netic waves: waves are generated by unstable particle populations and scatter these
populations, moving them in parameter space toward the equilibrium state. The
two main equations of the quasi-linear theory are the Fokker-Planck (diffusion)
equation

$$ \frac{\partial f_0}{\partial t} = \frac{\partial}{\partial p_{\perp}} \left( D_{\perp} \frac{\partial f_0}{\partial p_{\perp}} + D_{\parallel} \frac{\partial f_0}{\partial p_{\parallel}} \right) + \frac{\partial}{\partial p_{\parallel}} \left( D_{\parallel} \frac{\partial f_0}{\partial p_{\parallel}} + D_{\perp} \frac{\partial f_0}{\partial p_{\perp}} \right) \\ + \frac{D_{\perp \parallel}}{p_{\perp}} \frac{\partial f_0}{\partial p_{\parallel}} + \frac{D_{\perp \perp}}{p_{\perp}} \frac{\partial f_0}{\partial p_{\perp}} = \frac{\partial}{\partial \mathbf{p}} \left( \hat{D}_{pp} \frac{\partial f}{\partial \mathbf{p}} \right) $$

and the wave spectrum equation

$$ \frac{\partial \mathcal{U}_w}{\partial t} + (\nabla \cdot \mathbf{v}_g \mathcal{U}_w) = 2\gamma \mathcal{U}_w $$

where $\hat{D}_{pp}$ is the 2×2 tensor consisting of diffusion coefficients in momentum space ($p_{\parallel}$, $p_{\perp}$); and $\parallel, \perp$ are relative to the direction of the background magnetic field; $\mathcal{U}_w$ is the wave energy density, that can be expressed through the wave magnetic field energy $B_k^2$ and the Hermitian part of the dielectric tensor $\hat{\varepsilon}(\omega)$ [Shklyar and Matsumoto, 2009], as

$$ \mathcal{U}_w = \frac{1}{16\pi\omega} \left( \mathbf{E} \cdot \frac{d\omega^2 \hat{\varepsilon}(\omega)}{d\omega} \mathbf{E} \right) = \frac{1}{16\pi} \left( B_k^2 \mathbf{a} \cdot \frac{d\omega \hat{\varepsilon}(\omega)}{d\omega} \mathbf{a} + \mathbf{B}^2 \right) = B_k^2 Y, $$

where the vector $\mathbf{a}$ is given by the relationship between wave electric field vector and wave magnetic field magnitude $\mathbf{E}_w = \mathbf{a}B_k$ and is determined by the wave dispersion relation, as well as coefficient $Y$. Note that the mostly used quasi-linear equations are written for cyclotron resonance between charged particles and electromagnetic waves and, thus, describe the system averaged over electron gyrophase, i.e., these equations reduce the initially 3D momentum space to 2D ($p_{\parallel}, p_{\perp}$) space. Diffusion coefficients $\hat{D}_{pp}$ can be derived from the linear perturbation theory, the basic assumption of any quasi-linear model [Drummond and Pines, 1962; Kennel and Engelmann, 1966; Trakhtengerts, 1963; Vedenov et al., 1962]. Although the initial formulation of quasi-linear theory assumes that $\hat{D}_{pp}$ should be derived for a self-consistent wave spectrum $B_k^2$, the diffusion equation is often solved only for the most energetic (or relativistic) particle population, which usually does not contribute significantly to the variation of $B_k^2$. Therefore, statistical models of the waves, measured by spacecraft, can then be used instead of numerically evaluating the self-consistent evolution of $B_k^2$ [see examples in Glauert and Horne, 2005; Ni et al., 2008; Summers, 2005; Summers et al., 2007].

The quasi-linear equations describe only diffusion $\hat{D}_{pp} \propto B_w^2/B_0^2$ and drift $\hat{V}_p \nabla_p \hat{D}_{pp} \propto B_w^2/B_0^2$ with $B_w^2 = \langle B_k^2 \rangle$ the wave spectrum intensity. These two processes, diffusion and drift, are the results of integration of the wave Lorentz force along unperturbed particle trajectories and, thus, do not include any nonlinear effects. Nonlinearity consists in a significant role played by the wave Lorentz force in charged particle dynamics, and should manifest itself in a nonlinear dependence of diffusion on wave intensity, $\hat{D}_{pp} \propto (B_w/B_0)^\kappa$ with $\kappa \neq 2$, strong drifts $\hat{V}_p \nabla_p \hat{D}_{pp} \propto (B_w/B_0)^\kappa$ with $\kappa < 2$, and nondiffusive/drift terms in the full kinetic equation. The latter terms are the most difficult to include in the basic numerical models of radiation belt dynamics [see discussion in Artemyev, Vasiliev, Mourenas, Agapitov, Krasnoselskikh, Boscher and Rolland, 2014; Furuya et al., 2008; Omura et al., 2015]. To explain a possible generalization of the Fokker-Planck equation
suitable for including nonlinear wave-particle interactions, let us start with the general form of Smoluchowski coagulation equation [Sinitsyn et al., 2011]

$$ \frac{\partial f_0}{\partial t} = \int_{-\infty}^{\infty} \left( \hat{K}(\mathbf{p}|\mathbf{p}') f_0(\mathbf{p}') - \hat{K}(\mathbf{p}'|\mathbf{p}) f_0(\mathbf{p}) \right) d\mathbf{p}' \quad (1) $$

describing the evolution of $f(\mathbf{p})$. Here, $\hat{K}(\mathbf{p}|\mathbf{p}')$ is the coagulation kernel that describes the rate at which particle positions change in 2D ($p_\parallel, p_\perp$) space from $\mathbf{p}$ to $\mathbf{p}'$. Thus, the first term in Eq. (1) describes the particle flux toward $\mathbf{p}$ and the second term describes the particle fluxes away from $\mathbf{p}$.

Let us assume that each resonant interaction slightly changes particle momentum, $\Delta p_\parallel, \perp/p \ll 1$, so that we can expand $f(\mathbf{p}')$ as

$$ f_0(\mathbf{p}') = f_0(\mathbf{p} + \Delta\mathbf{p}) = f_0(\mathbf{p}) + \frac{\partial f_0}{\partial \mathbf{p}} \cdot \Delta\mathbf{p} + \frac{1}{2} \frac{\partial^2 f_0}{\partial \mathbf{p}^2} (\Delta\mathbf{p})^2 $$

Substituting this expression into Eq. (1), we obtain

$$ \frac{\partial f_0}{\partial t} = f_0(\mathbf{p}) \int_{-\infty}^{\infty} \left( \hat{K}(\mathbf{p}|\mathbf{p}') - \hat{K}(\mathbf{p}'|\mathbf{p}) \right) d\mathbf{p}' + \frac{\partial f_0}{\partial \mathbf{p}} \int_{-\infty}^{\infty} \hat{K}(\mathbf{p}|\mathbf{p}') \Delta\mathbf{p} d\mathbf{p}' \\ + \frac{1}{2} \frac{\partial^2 f_0}{\partial \mathbf{p}^2} \int_{-\infty}^{\infty} \hat{K}(\mathbf{p}|\mathbf{p}') (\Delta\mathbf{p})^2 d\mathbf{p}' = \mathbf{V} \frac{\partial f_0}{\partial \mathbf{p}} + \hat{D}_{pp} \frac{\partial^2 f_0}{\partial \mathbf{p}^2} = \frac{\partial}{\partial \mathbf{p}} \left( \hat{D}_{pp} \frac{\partial f_0}{\partial \mathbf{p}} \right) $$

where the last equality is provided by the divergence free condition $\mathbf{V} = \partial \hat{D}_{pp}/\partial \mathbf{p}$ [Lichtenberg and Lieberman, 1983; Sagdeev et al., 1988; Sinitsyn et al., 2011]. In this case, the Smoluchowski coagulation equation can be reduced to the Fokker-Planck diffusion equation:

$$ \frac{\partial f_0}{\partial t} = \frac{\partial}{\partial \mathbf{p}} \left( \hat{D}_{pp} \frac{\partial f_0}{\partial \mathbf{p}} \right) $$

This is the limit of the quasi-linear theory. If the resonant interaction is nonlinear, but can still be described by $\Delta p_\parallel, \perp/p \ll 1$, the diffusion equation is still usable, although diffusion coefficients should then be derived from test particle models including nonlinear wave field effects [see examples in Allanson et al., 2022; Frantsuzov et al., 2023; Inan, 1987; Karpman and Shklyar, 1977; Shklyar, 2021].

The most sophisticated case is when $\hat{K}$ describe large momentum changes $\Delta p_\parallel, \perp/p \sim 1$, and the integral operator $\sim \int \hat{K}$ cannot be reduced to the differential one. Such situation is common for electron nonlinear resonant interaction with intense whistler-mode waves [e.g., Agapitov et al., 2015b; Bortnik et al., 2008; Demekhov et al., 2006; Omura et al., 2007] and EMIC [e.g., Albert and Bortnik, 2009; Grach and Demekhov, 2018a, 2020; Omura and Zhao, 2012]. Figure 3 shows examples of electron trajectories with the large resonant changes of momentum due to the so-called phase trapping effect. A comparison of black (scattered electrons) and red (phase trapped electron) trajectories illustrates the main problem for the description of nonlinear resonant interactions with the differential operators of the
<figure><img src="image_3.png" /><figcaption>**Fig. 3** Examples of electron resonant interactions with a field-aligned whistler-mode wave: the effect of phase trapping is shown by red color, whereas all black trajectories show the phase bunching effect (small momentum decrease). Plot shows electron momentum $p = m_e c \sqrt{\gamma^2 - 1}$ for a half of an electron bounce period in a dipole magnetic field. System parameters are the same as in Fig. 2, but wave amplitude is 300pT.</figcaption></figure>

Fokker-Planck equation: the momentum change for a single resonant interaction (i.e., during the interval between electron trapping into the resonance and escape from this trapping) is comparable to the initial momentum amplitude. Therefore, the inclusion of such large momentum jumps into the Fokker-Planck equation would either require a significant decrease of the typical time-step of the simulation [such that the time-step of electron distribution evolution is much smaller than the electron bounce period, and phase trapping is modelled as a combination of small consecutive energy changes, see examples in Foster et al., 2017; Shklyar, 1981], or the development of a non-differential (integral) operator describing the large energy change on the smallest system time-scale, during an electron bounce period. Therefore, the main challenge for radiation belt models is to construct $\hat{K}$ or find an approach for taking into account the effects of the integral operator $\sim \int \hat{K}$ into the diffusive Fokker-Planck equation. This review is devoted to possible solutions to this challenge.

The most direct, and quite effective, approach consists in numerically evaluating the $\hat{K}$ function. This approach has been applied for test particle simulations of electron interactions with EMIC and electrostatic waves [Artemyev, Vasiliev and Neishtadt, 2019; Zheng et al., 2019], but the most elaborate variant of this approach, called the *Green function* approach, has been proposed for electron resonant interactions with chorus waves [Furuya et al., 2008; Omura et al., 2015]. This is the most developed and advanced approach accounting for multiple properties of resonances: wave frequency drift [Hsieh and Omura, 2017b], wave propagation in the form of a train of short wave packets [Hiraga and Omura, 2020; Kubota and Omura, 2018], wave oblique propagation [Hsieh et al., 2020; Hsieh and Omura, 2017a], multiple resonances due to wave obliqueness [Hsieh and Omura, 2023; Hsieh et al., 2022]. The main advance of this approach is that the numerical evaluation of $\hat{K}$ can be performed for an arbitrary and very realistic wave field
model. The main disadvantage is that such numerical evaluation requires a discretization of the momentum space with a sufficient statistics (sufficiently large number) of resonant interactions inside each bin, whereas the lowest bin size is determined by weak diffusive scattering $\Delta p \propto B_w/B_0$ and the range of $\Delta p$ is determined by the phase trapping property with $\Delta p \propto O(B_w/B_0) \gg B_w/B_0$. Thus, a purely numerical evaluation of $\hat{K}$ may sometimes miss some weak diffusive effects, and this approach should be mostly effective for modeling brief events with not-widely-varying system characteristics (to avoid a recalculation of $\hat{K}$ for multiple realizations of system parameters).

In this review, we examine the theoretical properties of the $\hat{K}$ function, and explore different approaches for its analytical evaluation (Sections 2, 3, and 4). We provide a detailed investigation of $\hat{K}$ for nonlinear electron interactions with monochromatic intense whistler-mode waves, and provide asymptotic solutions for the kinetic equation including such $\hat{K}$ (Sections 4). Then, we generalize $\hat{K}$ for systems with a large wave ensemble, and perform such a generalization via the mapping technique for nonlinear resonant interactions (Section 4.4). In Appendix E, we provide several examples of application of this technique for simulations of the observed dynamics of the electron flux. The next natural generalization of the theoretical approach for nonlinear wave-particle interactions consists of the inclusion of the effects of short wave-packets. We discuss the main aspects of this generalization (Section 5) and of the theoretical approaches for inclusion of short wave-packets into the mapping technique (Section 6). Next, we consider an approach allowing the incorporation of nonlinear resonant interactions into existing global numerical models of the radiation belts (Section 7). Finally we discuss several aspects of nonlinear resonant interactions that are not included in this review, but can be important for specific plasma systems (Section 8). The review also contains Appendix A, with the main equations of the Hamiltonian approach for wave-particle resonant interactions, Appendix B considering a special situation of nonlinear interactions for field-aligned particles, Appendix C describing analytical estimates for electron resonant interaction with short wave-packets, Appendix D describing the problem of the electron phase gain between two resonances, and Appendix E with several examples of observations of nonlinear resonant effects.

## 2 Basic properties of electron resonant interactions

We start with general information about whistler-mode waves observed in the inner magnetosphere, and specifically within the outer radiation belt, outside the plasmasphere. These right-hand circularly polarized electromagnetic waves are mainly generated in the frequency range from 0.1 to 0.7 times the equatorial electron gyrofrequency $\Omega_0$ under the form of repetitive rising tones, and have been called chorus waves [Burton and Holzer, 1974; Helliwell, 1965; Storey, 1953; Tsurutani and Smith, 1974]. There are two main modes of these waves: electromagnetic mode with nearly parallel propagation relative to the background magnetic field and quasi-electrostatic mode with strongly oblique propagation. Figure 4 shows examples of both wave modes. More detailed information about statistics of these two wave modes and their relative occurrence rates can be found in [Agapitov et al., 2013, 2018; Artemyev, Agapitov, Mourenas, Krasnoselskikh, Shastun and
Mozer, 2016; Li, Bortnik, Thorne and Angelopoulos, 2011; Li, Santolik, Bortnik, Thorne, Kletzing, Kurth and Hospodarsky, 2016; Meredith et al., 2012].

Field-aligned (i.e., propagating parallel to the background magnetic field) whistler-mode waves are generated at the equator by transversely anisotropic electron populations [Kennel, 1966; Sagdeev and Shafranov, 1961], which are either injected from the plasma sheet [Fu et al., 2014; Tao, Thorne, Li, Ni, Meredith and Horne, 2011] or generated by dayside magnetosphere compression [Li et al., 2015]. After an initial linear wave growth [Kennel, 1966], nonlinear wave growth takes over once the generated wave reaches a threshold amplitude for electron trapping in the inhomogeneous magnetic field, leading to the formation of characteristic rising or falling tone elements [see reviews in Demekhov, Taubenschuss and Santolík, 2017; Helliwell and Crystal, 1973; Nunn, 1974; Omura, 2021; Omura et al., 2008, 2013; Tao et al., 2021, 2020]. The electron azimuthal drift from the injection region to the day side and such day side compression determine the domain of presence of near-equatorial field-aligned whistler-mode waves [Agapitov et al., 2013; Li, Bortnik, Thorne and Angelopoulos, 2011; Li, Bortnik, Thorne, Cully, Chen, Angelopoulos, Nishimura, Tao, Bonnell and Lecontel, 2013; Meredith et al., 2012]. Propagating away from their equatorial source region, these waves become oblique [Alekhin and Shklyar, 1980; Bell et al., 2002; Bortnik et al., 2006; Shklyar et al., 2004] and experience Landau damping by suprathermal electrons [Bortnik et al., 2007; Chen et al., 2013; Watt et al., 2013]. Such damping is stronger on the night side due to larger magnetic field line curvature, leading to a confinement of these waves near the equator (< 15°), whereas on the day side field-aligned and weakly oblique waves may propagate up to middle latitudes of ~ 30° [Agapitov et al., 2013, 2018]. A potentially important sub-population of field-aligned waves consists of ducted waves, which are trapped within plasma density perturbations and can propagate without damping to high latitudes [Karpman and Kaufman, 1982; Laird and Nunn, 1975]. Such waves have been observed in-situ [Chen, Gao, Lu, Chen, Tsurutani, Li, Ni and Wang, 2021; Chen, Gao, Lu, Tsurutani and Wang, 2021; Streltsov and Bengtson, 2020], reproduced in numerical simulations [Hanzelka and Santolík, 2019; Ke et al., 2021; Streltsov and Goyal, 2021], and detected by ground-based stations [Collier et al., 2011; Demekhov, Manninen, Santolík and Titova, 2017; Martinez-Calderon et al., 2020, 2015; Titova et al., 2017, 2015]. However, the population of ducted whistler-mode waves has not yet been precisely quantified and their occurrence rate in each region is not known [see discussion in Artemyev, Demekhov, Zhang, Angelopoulos, Mourenas, Fedorenko, Maninnen, Tsai, Wilkins, Kasahara, Miyoshi, Matsuoka, Kasahara, Mitani, Yokota, Keika, Hori, Matsuda, Nakamura, Kitahara, Takashima and Shinohara, 2021; Artemyev et al., 2024; Zhang, Meng, Artemyev, Zou and Mourenas, 2023].

Very oblique waves observed at high latitudes likely result from the diffraction of initially field-aligned waves during their propagation along the inhomogeneous magnetic field [Agapitov et al., 2013; Breuillard et al., 2014, 2012; Chen et al., 2013]. However, additionally to this high-latitude population, there are also near-equatorial very oblique waves [Agapitov et al., 2013; Cattell et al., 2008; Cully, Bonnell and Ergun, 2008], which are likely generated by transversely anisotropic electrons in the presence of field-aligned electron streams that reduce the Landau damping [Artemyev, Agapitov, Mourenas, Krasnoselskikh, Shastun and Mozer, 2016; Gao et al., 2016; Li, Mourenas, Artemyev, Bortnik, Thorne, Kletzing, Kurth, Hospodarsky, Reeves, Funsten and Spence, 2016; Mourenas et al., 2015]. Genera-
tion of very oblique waves, propagating around the resonance cone angle [Sazhin, 1993], require specific distributions of suprathermal (100 eV to a few keV) electrons with a plateau in parallel velocity space [Artemyev, Agapitov, Mourenas, Krasnoselskikh, Shastun and Mozer, 2016; Chen et al., 2019; Ke et al., 2022; Kong et al., 2021; Mourenas et al., 2015]. Such field-aligned electron streams can be formed either by electrostatic parallel fields often observed around plasma sheet injections [see discussion in Artemyev and Mourenas, 2020] or by ionosphere outflow [see discussion in Artemyev, Zhang, Angelopoulos, Mourenas, Vainchtein, Shen, Vasko and Runov, 2020]. Although both scenarios assume specific conditions for very oblique wave generation, this wave population is quite widespread in observations [Agapitov et al., 2013; Li, Santolik, Bortnik, Thorne, Kletzing, Kurth and Hospodarsky, 2016] and important for energetic electron flux dynamics [Agapitov et al., 2015b; Artemyev, Agapitov, Mourenas, Krasnoselskikh and Mozer, 2015; Artemyev, Agapitov, Mourenas, Krasnoselskikh and Zelenyi, 2013; Hsieh et al., 2020; Mourenas, Artemyev, Agapitov and Krasnoselskikh, 2014]. Nevertheless, intense very oblique waves are rarely observed simultaneously with intense field-aligned waves, probably due to Landau damping and nonlinear effects [Agapitov et al., 2016].

Electron resonant interactions with these two wave modes are quite different [Albert, 2017; Artemyev, Agapitov, Mourenas, Krasnoselskikh and Mozer, 2015; Bell, 1984, 1986; Shklyar and Matsumoto, 2009; Solovev and Shkliar, 1986]. Therefore, nonlinear effects will be considered separately for each wave mode. Both wave modes share an important property – they are coherent and quite narrow band waves, with a high intensity. The importance of this property will become clear if we consider the applicability criteria for quasi-linear theory. This theory is based on the concept of resonance overlap for a wide spectrum of waves [Shapiro and Sagdeev, 1997]. In momentum space, the resonance width is $\Delta p_R \approx \sqrt{eU_w m_e}$ where $U_w = A_w \sqrt{2I_x \Omega_0 / m_e c^2}$ for the cyclotron resonance with electromagnetic waves [Karimabadi et al., 1990] and $U_w = \phi$ for Landau resonance with electrostatic waves [Palmadesso, 1972]. The width of the Landau resonance for an electrostatic wave does not depend on electron characteristics and is entirely determined by the wave electrostatic potential, $\phi$, whereas the width of the cyclotron resonance depends not only on the wave vector potential, $A_w = B_w/k$, but also on electron magnetic moment $I_x$ (i.e., on pitch-angle and energy). Note, however, that Landau resonance with electromagnetic waves is also characterized by a resonant width depending on wave characteristics [Shklyar and Matsumoto, 2009]. The distance between resonances with two nearby waves in the spectrum would be $\Delta p_\omega = m_e (\Delta\omega/k) \cdot (1 - p_R/m_e v_g)$, where $\Delta\omega$ is the distance between waves (e.g., spectral width of wave-packet), and $v_g = \partial\omega/\partial k$ is the wave group velocity [Karpman, 1974]. The overlap condition requires that there are many $\Delta p_R$ within $\Delta p_\omega$, i.e., that the wave spectrum be wide enough ($\Delta\omega$ is large) or that the wave amplitude be weak enough ($U_w$ is small). This condition is generally satisfied for low amplitude whistler-mode waves observed in the near-Earth plasma sheet [Gao et al., 2022; Waheed et al., 2023], but it is often not satisfied for narrow band intense whistler-mode chorus waves in the plasma injection regions and Earth's outer radiation belt, see Fig. 5. Besides quite small wave spectrum width, $\Delta\omega/\omega$, whistler-mode wave packets have peak amplitudes of about $\in [10^{-2}, 10^{-3}] B_0$ [e.g., Agapitov et al., 2014; Tyler et al., 2019; Wilson et al., 2011; Zhang, Mourenas, Artemyev, Angelopoulos, Bortnik, Thorne, Kurth, Kletz-
<figure><img src="image_4.png" /><figcaption>**Fig. 4** Two examples of whistler-mode waves (from top to bottom): (a) wave electric field spectrum, (b) wave magnetic field spectrum, (c) wave normal angle, (d) examples of individual wave-packets (electric field components are shown), (e) examples of electric field polarization (two transverse components are plotted). Left panels show the example of field-aligned intense waves, and right panels show the example of very oblique waves. In panels (a-c) $f_{ee} = \Omega_0/2\pi$ is the electron gyrofrequency. Data are collected by Van Allen Radiation Belt Probes (RBSP) [Mauk et al., 2013] in the outer radiation belt (see coordinates in L-shell, magnetic local time, and latitudes below panels (c)). We use wave field measurements by the Electric and Magnetic Field Instrument Suite and Integrated Science (EMFISIS) on RBSP [Kletzing et al., 2013]. Details of wave data processing can be found in [Agapitov et al., 2015b, 2016].</figcaption></figure>

ing and Hospodarsky, 2019], that is, a factor ×100 higher than mean wave amplitudes derived from the averaged wave spectra (e.g., compare wave packet statistics in [Zhang, Thorne, Artemyev, Mourenas, Angelopoulos, Bortnik, Kletzing, Kurth and Hospodarsky, 2018] or [Zhang, Mourenas, Artemyev, Angelopoulos, Bortnik, Thorne, Kurth, Kletzing and Hospodarsky, 2019] with time-averaged wave statistics in [Agapitov et al., 2018]). Such high intensity wave packets may nonlinearly interact with electrons through cyclotron or Landau resonances.
<figure><img src="image_5.png" /><figcaption>Fig. 5 Distribution of statistics of field-aligned intense whistler-mode waves in $(\Delta\omega/\omega, \sigma)$ parametric space where $\Delta\omega$ is the spectrum width, $\sigma = \sqrt{B_w/B_0\beta_e^{1/4}}(\Omega_0/\omega - 1)^{-1/4}$ and $\beta_e$ is the electron plasma beta parameter (ratio of electron thermal pressure and background magnetic field pressure). Curves show the $\sigma = \Delta\omega/\omega$ threshold. Panel (a) shows statistics of field-aligned whistler-mode wave packets in the Earth's outer radiation belt [see Zhang, Mourenas, Artemyev, Angelopoulos, Bortnik, Thorne, Kurth, Kletzing and Hospodarsky, 2019]; panel (b) shows statistics of whistler-mode waves observed in the plasma sheet, around injection regions [Zhang, Angelopoulos, Artemyev and Liu, 2018]; panel (c) shows statistics of very oblique whistler-mode waves observed in the outer radiation belt (we use THEMIS [Angelopoulos, 2008] electric field and magnetic field [Bonnell et al., 2008; Le Contel et al., 2008] wave measurements and the same criteria for wave identification as [Li, Santolik, Bortnik, Thorne, Kletzing, Kurth and Hospodarsky, 2016]).</figcaption></figure>

Therefore, wave-particle resonant interactions should be considered under the assumption that electrons interact resonantly with individual intense waves. The corresponding Hamiltonian for wave particle interactions in dipole field can be written as (see Appendix A and [Albert, 1993; Albert et al., 2013; Vainchtein et al., 2018]):

$$H = m_e c^2 \gamma - eU_w (I_x, s) \cos (\phi - n_r \psi), \quad \gamma = \sqrt{1 + \left(\frac{p_{\parallel}}{m_e c}\right)^2 + \frac{2I_x \Omega_0 (s)}{m_e c^2}} \quad (2)$$

where $\Omega_0(s) = eB_0(s)/m_e c$ is the electron gyrofrequency, $(s, p_{\parallel})$ are conjugated field-aligned coordinate and momentum, $(\psi, I_x)$ are conjugated gyrophase and magnetic moment $(I_x = (\gamma^2 - 1)m_e c^2 \sin^2 \alpha/\Omega_0$ and $\alpha$ is a local pitch-angle), $n_r = 0, \pm 1, \pm 2, ...$ is the resonance number. The resonance condition $\phi = n_r \psi$, the wave phase definition $\dot{\phi} = k_{\parallel}(s)\dot{s} - \omega$, and Hamiltonian equations $\dot{\psi} = \partial H/\partial I_x = \Omega_0/\gamma$, $\dot{s} = \partial H/\partial p_{\parallel} = p_{\parallel}/m_e \gamma$ determine the resonant momentum $p_{\parallel} = p_R$:

$$p_R = m_e \frac{\gamma\omega + n_r\Omega_0}{k_{\parallel}}$$

The function $U_w$ is the generalized wave amplitude, including effects of whistler-mode dispersion (see Appendix A and [Albert, 1993; Artemyev, Neishtadt, Vasiliev
and Mourenas, 2018)]:

$$U_w = A_w \sum_{n_r=-\infty}^{+\infty} h^{(n_r)}$$

and

$$h^{(n_r)} = \frac{\rho \Omega_0}{2c\gamma} \left( (C_1 - \cos \theta) J_{n_r-1} - (C_1 + \cos \theta) J_{n_r+1} \right) - \left( \frac{p_{\parallel}}{\gamma m_e c} + C_2 \right) \sin \theta J_{n_r}$$

where $J_{n_r}$ are Bessel functions with argument $k_{\perp} \rho$, $\rho = \sqrt{2I_x/m_e \Omega_0}$, $\theta$ is the wave normal angle (wave number has two components: field-aligned $k_{\parallel} = k \cos \theta$ and transverse $k_{\perp} = k \sin \theta$). This equation for $u_w^{(n_r)}$ has been derived using the relations between wave magnetic and electric field components in a cold plasma [Tao and Bortnik, 2010; Williams and Lyons, 1974] and, thus, the coefficients $C_{1,2}$, given in Appendix A, depend only on wave characteristics. An important property of $C_1$ is that for $\theta = 0$ (field-aligned propagation) it is equal to one, whereas $U_w^{(n_r)} = 0$ for all $n_r$ except $n_r = -1$, for which we have $h^{(-1)} = \rho \Omega_0/c\gamma = \sqrt{2I_x \Omega_0 / m_e c^2}$.

To describe the wave dispersion, $\omega = \omega(k, \theta)$, we again use the cold plasma approximation [Stix, 1962]:

$$N^2 = \frac{1}{2} \frac{(RL - PS) \sin^2 \theta + 2PS}{S \sin^2 \theta + P \cos^2 \theta} + \frac{1}{2} \frac{\sqrt{(RL - PS)^2 \sin^4 \theta + 4P^2 D^2 \cos^2 \theta}}{S \sin^2 \theta + P \cos^2 \theta}$$

where $N = kc/\omega$ is the refractive index, and Stix coefficients are

$$R = 1 - \sum_j \frac{\Omega_{pj}^2}{\omega \omega + \Omega_{0j}}, \quad L = 1 - \sum_j \frac{\Omega_{pj}^2}{\omega \omega - \Omega_{0j}}$$

$$S = \frac{1}{2}(R+L), \quad D = \frac{1}{2}(R-L), \quad P = 1 - \sum_j \frac{\Omega_{pj}^2}{\omega^2}$$

Here $\Omega_{0j} = e_j B/m_j c$ is the gyrofrequency for particles with a charge $j = e$ for electrons and $j = i$ or $j = p$ for ions or protons, $\Omega_{pj} = \sqrt{4\pi n_0 e_j^2 / m_j}$, $n_0$ is the background density ($n_0 = n_e = n_i$). Note $e_j = -e$ for electrons and $e_j = e$ for protons.

A simplified form of this dispersion relation for a dense plasma is

$$\omega^2 = \frac{\Omega_0^2 \cos^2 \theta}{(1 + (\Omega_{pe}/kc)^2)^2} + \frac{\nu \Omega_0^2}{1 + (\Omega_{pe}/kc)^2}$$

where $\Omega_{pe}(s) \gg \Omega_0(s)$ is the electron plasma frequency, $\Omega_0 = \Omega_{0e}$, and $\nu = m_e/m_i$, $m_i$ is the mass of the ion mixture (i.e., this is the proton mass for purely proton-electron plasma). For quasi-parallel wave propagation, $\theta \sim 0$, the last term in the simplified dispersion relation can be omitted:

$$\omega = \frac{\Omega_0 \cos \theta}{1 + (\Omega_{pe}/kc)^2}$$

Let us consider the Hamiltonian (2) under the assumption of constant wave frequency, $\omega = \text{const.}$ Although this assumption is not applicable for the most
intense whistler-mode chorus waves [Omura, 2021; Omura et al., 2008; Tao et al., 2021], it is quite useful to describe the basic properties of wave-particle interactions. In this case, the time dependence is included only as a linear term $\omega t$ in the wave phase $\phi$ and, therefore, we may remove this dependence (to obtain a conservative Hamiltonian) by changing the variables: from $(\psi, I_x)$ to $(\zeta, I)$ where $\zeta = n_r\psi + \phi$. The corresponding generating function is

$$W = (\phi - n_r\psi) I + sp + \psi \tilde{I}_x$$

and the new Hamiltonian is

$$\mathcal{H}_I = -\omega I + m_e c^2 \gamma - eU_w (\tilde{I}_x - n_r I, s) \cos \zeta \quad (3)$$

$$\gamma = \sqrt{1 + \left(\frac{p + k_{\parallel} I}{m_e c}\right)^2 + \frac{2(\tilde{I}_x - n_r I) \Omega_0(s)}{m_e c^2}}$$

where $p = p_{\parallel} - k_{\parallel} I$ is the new momentum conjugated to the new coordinate $\tilde{s} = s$ (keeping the $s$ notation), $I_x = \tilde{I}_x - n_r I$, and $I$ is conjugated to the new phase $\zeta$ [see, e.g., Artemyev, Neishtadt, Vasiliev and Mourenas, 2018].

The Hamiltonian (3) does not depend on time, and thus $\mathcal{H}_I=\text{const}$, with the integral of electron motion

$$h = m_e c^2 \gamma - \omega I \quad (4)$$

Note that $\tilde{I}_x$ is a constant, because $\mathcal{H}_I$ does not depend on $\psi$. For the case of Landau resonance, $n_r = 0$, $\tilde{I}_x = I_x$, i.e., $I_x$ is conserved. Using this conservation law in Eq. (4) gives $h = \Omega_{eq}I_x = \text{const}$. For the first cyclotron resonance, $n_r = -1$, $\tilde{I}_x = I_x - I$ and can be set equal to zero, i.e. $I_x = I$. In this case, Eq. (4) can be written as:

$$h = m_e c^2 \gamma - \omega I_x$$

The integral given by Eq. (4) describes trajectories, in the momentum space or in energy, pitch-angle space $(E, \alpha_{eq})$, along which wave-particle interactions are transporting electrons. Note that without such interactions, the energy $E = m_e c^2(\gamma-1)$ and equatorial pitch-angle $\alpha_{eq} = \arcsin(2I_x\Omega_{eq}/(\gamma^2 - 1))$ would be constants of motion. Figure 6 shows these trajectories (also called resonance curves; e.g., [Walker, 1993]) for the Landau ($n_r = 0$) and first cyclotron ($n_r = -1$) resonances of electrons with whistler-mode waves. The main feature of the cyclotron resonance is that the resonance curves are almost parallel to the pitch-angle axis for small pitch-angles, and show energy increase with pitch-angle increase for larger pitch-angles. Thus, electrons transported to smaller pitch-angles (and finally into the loss-cone, thereby precipitating into the atmosphere) lose their energy, but around the loss-cone (at small pitch-angles) this transport occurs almost without energy change. The main feature of the Landau resonance is the large energy increase of electrons transported along the resonant curves toward smaller pitch-angles, i.e., electron precipitation due to the Landau resonance is accompanied by electron acceleration. In contrast to the cyclotron resonance, electron transport to higher pitch-angles via the Landau resonance is associated with electron energy loss. For the Landau resonance, the shape of resonance curves is dictated by magnetic moment conservation and does not depend on the wave frequency, whereas
<figure><img src="image_6.png" /><figcaption>Fig. 6 Resonance curves given by Eq. (4) for Landau ($n_r = 0$; dashed black line) and first cyclotron ($n_r = -1$; solid red and blue lines) resonances. Two wave frequencies are shown: $\omega/\Omega_{eq} = 0.15$ (red) and $\omega/\Omega_{eq} = 0.35$ (blue).</figcaption></figure>

for the cyclotron resonance a smaller wave frequency means less energy change for the same pitch-angle change. At high energies for cyclotron resonance, the resonance curves show a change of direction: there, energy increase corresponds to pitch-angle decrease. This is the so-called *turning acceleration* effect [Omura et al., 2007], which occurs when $\gamma\omega/\Omega_0$ is larger than one and $d\gamma/d\alpha$ changes sign for large $\alpha$:

$$ \frac{d\gamma}{d\alpha} = \frac{(\omega/\Omega_0)(\gamma^2 - 1)\sin\alpha \cos\alpha}{1 - (\omega\gamma/\Omega_0)\sin^2\alpha} $$

This equation is obtained by differentiation of Eq. (4) over $\alpha$. Note that $\gamma\omega/\Omega_0 = 1$ also separates cyclotron resonances with negative resonant momentum (waves and resonant particles are moving in opposite directions) and with positive resonant momentum (waves and resonant particles are moving in the same direction), $p_R = m_e\gamma\omega(1 - \Omega_0/\gamma\omega)/k_\parallel$.

The resonance curves in Fig. 6 show that for the same wave characteristics, electrons with different energies and pitch-angles (i.e., with different $h$ integrals of motion) can interact resonantly with a same monochromatic ($\omega = \text{const}$) wave. This is due to background magnetic field inhomogeneity: electrons bounce along magnetic field lines and have a different local pitch-angle $\alpha(s) = \arcsin(\sin(\alpha_{eq})\sqrt{\Omega_0(s)/\Omega_{eq}})$ at different magnetic latitudes (coordinate $s$). Accordingly, electrons with a fixed energy will be able to find a specific latitude where their pitch-angles will satisfy the resonance condition with a whistler-mode wave of fixed frequency. There are two important consequences of such magnetic field inhomogeneity: (1) a monochromatic wave may interact resonantly with electrons in a wide energy, pitch-angle range (in contrast to a homogeneous plasma, where only a wide wave spectrum may provide resonances over a wide energy, pitch-angle range), (2) electron bounce motion moves particles into the resonances and moves them out of the resonances, which are therefore of limited size in latitude, and of limited duration. This second effect actually replaces the limitation of the duration of the resonant interaction
of electrons with whistler-mode waves that is due to the existence of a broad wave spectrum in classical quasi-linear models [see discussion in Albert, 2001, 2010; Allanson et al., 2022; Karpman, 1974; Shklyar, 1981, 2011].

## 2.1 Resonant interaction with monochromatic wave

Let us now consider the effects of a high wave amplitude on resonant particle dynamics for an arbitrary wave propagation angle. We start with the Hamiltonian (3) and follow the standard procedure for analysis of resonant Hamiltonian systems [Albert et al., 2013; Artemyev, Neishtadt, Vasiliev and Mourenas, 2018; Neishtadt, 2014; Neishtadt and Vasiliev, 2005]. First, let us write the equation for the resonant condition, $\dot{\zeta} = 0$

$$ \dot{\zeta} = \frac{\partial H_I}{\partial I} \approx -\omega + m_e \frac{n_r \Omega_0 + k_{\parallel} (p + k_{\parallel} I)}{m_e \gamma} = 0 $$

where we omit the term

$$ \sim \frac{\partial eU_w}{\partial I} \cos \zeta $$

that is proportional to the wave amplitude and should not be included into the definition of the resonance (see discussion of exceptions in [Li et al., 2022] and discussion of the importance of this term in Appendix B).

Equation $\dot{\zeta} = 0$ has a solution $I = I_R$, where

$$ \frac{k_{\parallel} I_R}{m_e c} = -\frac{p}{m_e c} - \frac{n_r \Omega_0}{k_{\parallel} c} + \frac{1}{\sqrt{N_{\parallel}^2 - 1}} \sqrt{1 + \frac{2 \tilde{I}_x \Omega_0}{m_e c^2} - \left(\frac{n_r \Omega_0}{k_{\parallel} c}\right)^2 - 2 \frac{n_r \Omega_0}{k_{\parallel} c} \frac{p}{m_e c}} $$

where $N_{\parallel} = k_{\parallel}c/\omega$.

Expanding Hamiltonian (3) around $I = I_R$, we obtain

$$ \mathcal{H}_I \approx \Lambda + \frac{1}{2M} (I - I_R)^2 + eU_w (I_R, s) \cos \zeta $$

$$ \Lambda = -\omega I_R + m_e c^2 \gamma_R, \quad \gamma_R = \frac{N_{\parallel}}{\sqrt{N_{\parallel}^2 - 1}} \sqrt{1 + \frac{2 \tilde{I}_x \Omega_0}{m_e c^2} - \left(\frac{n_r \Omega_0}{k_{\parallel} c}\right)^2 - 2 \frac{n_r \Omega_0}{k_{\parallel} c} \frac{p}{m_e c}} $$

In Appendix A, we describe the variable change and the corresponding generating function: $(\zeta, I) \rightarrow (\zeta, P_{\zeta})$ where $P_{\zeta} = I - I_R$. After this change and an expansion relative to the resonance, the Hamiltonian $\mathcal{H}_I$ is separated into two parts. The first part

$$ \Lambda = -\omega I_R(s, p) + m_e c^2 \gamma_R(s, p) \quad (5) $$

describes the dynamics of slow variables, $(s, p)$. The second part

$$ \mathcal{H}_{\zeta} = \frac{1}{2M} P_{\zeta}^2 + A\zeta + B \cos \zeta \quad (6) $$

describes the dynamics of fast variables $(\zeta, P_{\zeta})$, with coefficients depending on slow variables $(p, s)$. Coefficients $M$, $A$, and $B$ are derived in Appendix A for the general
case of arbitrary wave propagation direction. For the first cyclotron resonance with field-aligned waves, these coefficients are

$$A = \{\Lambda, I_R\}_{s,p} \approx -\frac{m_e c^2 D}{\gamma_R} \frac{N_{\parallel}^2}{N_{\parallel}^2 - 1} \left( \left( \frac{p_{\parallel,R}}{m_e c} \right)^2 \frac{\partial \ln N_{\parallel}}{\partial \ln \Omega_0} + \frac{\Omega_0}{\omega N_{\parallel}} \frac{p_{\parallel,R}}{m_e c} - \frac{\Omega_0 I_{x,R}}{m_e c^2} \right) \\ B = eU_w(I_R, s_R), \quad M = \left( \frac{\partial^2 \mathcal{H}_I}{\partial I^2} \right)_{I_R}^{-1} = \frac{m_e c^2}{\omega^2} \frac{\gamma_R}{N_{\parallel}^2 - 1}$$

where $N_{\parallel} = k_{\parallel}c/\omega$, $p_{\parallel,R} = (\gamma_R - n_r\Omega_0/\omega)/N_{\parallel}$, $D = c(\partial \ln \Omega_0/\partial s)/N_{\parallel}\omega \ll 1$ is a dimensionless factor of system inhomogeneity, and $s_R$ and $I_{x,R}$ are coordinate and momentum in the resonance.

These coefficients determine the character of wave-particle resonant interactions. Figure 7 shows the phase portrait for two regimes: $a = |B/A| < 1$ and $a = |B/A| > 1$, where

$$\frac{B}{A} = -\frac{eU_w(I_R, s_R) \gamma_R N_{\parallel}^2 - 1}{m_e c^2 D} \left( \left( \frac{p_{\parallel,R}}{m_e c} \right)^2 \frac{\partial \ln N_{\parallel}}{\partial \ln \Omega_0} + \frac{\Omega_0}{\omega N_{\parallel}} \frac{p_{\parallel,R}}{m_e c} - \frac{\Omega_0 I_{x,R}}{m_e c^2} \right)^{-1}$$

For $a < 1$ the phase portrait is filled by trajectories crossing the resonance, $\dot{\zeta} = P_c/M = 0$, only once. These are so-called transient trajectories, and as particles moving along these trajectories spend a quite limited time around the resonance, the wave field cannot significantly change their orbits. Such regime of wave-particle resonant interactions is principally similar to linear scattering, which is evaluated for unperturbed particle trajectories [Albert, 2001; Karpman and Shklyar, 1977].

For $a > 1$ the phase portrait is divided into two domains: the internal domain is filled by closed trajectories with multiple resonance $\dot{\zeta} = 0$ crossings, whereas the external domain is filled by open trajectories with a single resonance crossing, but particles on these trajectories need to move around the interval domain and thus spend much more time around the resonance in comparison with particles moving along transient trajectories. Particles inside the internal domain are called phase trapped and may spend a very long time in the resonance (oscillating around the resonance). Particles in the external domain will be scattered, but this scattering cannot be evaluated under the approximation of unperturbed trajectories. This nonlinear scattering is very different from the linear one, appropriate for transient trajectories. The entire problematic of nonlinear wave-particle interaction consists in developing an accurate description of such nonlinear scattering and phase trapping effects for a large ensemble of electrons.

The factor $1/a$ is essentially the same parameter as the S-parameter used in previous studies of wave-particle resonant interactions by Helliwell [1967]; Omura et al. [2009, 2008] and as the $\tau/\alpha$-parameter used by Karpman et al. [1975, 1974]; Shklyar [2011]; Shklyar and Matsumoto [2009]. We can rewrite the expression for
<figure><img src="image_7.png" /><figcaption>Fig. 7 Phase portraits for Hamiltonian (6) for $a = |B/A| < 1/2$ (a) and $a = |B/A| = 2$ (b).</figcaption></figure>

a as

$$ \begin{aligned} \frac{1}{a} &= \frac{m_e c^2}{eB_w \gamma_R} \sqrt{\frac{m_e c^2}{2\Omega_0 I_x}} \frac{k}{\omega} \frac{N_\parallel}{N_\parallel^2 - 1} \left( \frac{\Omega_0 I_x}{m_e c^2} - \left( \frac{p_\parallel, R}{m_e c} \right)^2 \frac{\partial \ln N_a}{\partial \ln \Omega_0} - \frac{\Omega_0}{\omega N_a} \frac{p_\parallel, R}{m_e c} \right) c \frac{\partial \ln \Omega_0}{\partial s} \\ &= \frac{c}{s'_0 \Omega_w \omega} \frac{\partial \Omega_0}{\partial s} \frac{1}{2 \xi' \delta'} \left( \frac{\gamma_R \omega}{\Omega_0} \left( \frac{V'_{\perp 0}}{c} \right)^2 - 2 \gamma_R \frac{\omega}{\Omega_0} \left( \frac{V'_R}{c} \right)^2 \frac{\partial \ln N_\parallel}{\partial \ln \Omega_0} - 2 \frac{V'_R V'_p}{c^2} \right) \\ &= \frac{c}{s'_0 \Omega_w \omega} \frac{\partial \Omega_0}{\partial s} \frac{1}{2 \xi' \delta'} \left( \frac{\gamma_R \omega}{\Omega_0} \left( \frac{V'_{\perp 0}}{c} \right)^2 - \left( 2 + \frac{\omega}{\Omega_0} \frac{\Omega_0 - \gamma_R \omega}{\Omega_0 - \omega} \right) \frac{V'_R V'_p}{c^2} \right) \end{aligned} $$

where we used notation from [Omura et al., 2009]:

$$ \begin{aligned} s'_0 &= \frac{N_\parallel^2 - 1}{N_\parallel} \frac{V_{\perp 0}}{c}, \quad \xi' = \sqrt{N_\parallel^2 - 1}, \quad \delta' = \sqrt{1 - N_\parallel^{-2}} \\ V'_R &= p_\parallel, R/m_e \gamma_R, \quad V'_{\perp 0} = c \sqrt{2I_x \Omega_0 / m_e c^2 / \gamma_R}, \quad V'_p = \omega / k = c / N_\parallel \end{aligned} $$

and a simplified dispersion relation with constant plasma density (see Appendix A):

$$ \frac{\partial \ln N_\parallel}{\partial \ln \Omega_0} = -\frac{1}{2} \frac{\omega}{\Omega_0 - \omega} $$

In this form, $1/a$ fully coincides with $S$ given for $\partial\omega/\partial t = 0$ by Eqs. (10)-(13) from Omura et al. [2009].

This parameter determines whether the wave field $\sim B$ is sufficiently strong to overcome the inhomogeneity of the background magnetic field $\sim A$. For $a > 1$, the wave-particle resonant interaction is nonlinear, and this nonlinearity both changes the electron diffusion and introduces new effects of phase bunching and trapping. The parameter $a$ depends on background plasma and wave characteristics, but also on electron energy and pitch-angle. Figure 8 shows the probability distribution of $a > 1$ for field-aligned whistler-mode waves in the outer radiation belt (this probability is displayed in $(E, \alpha_{eq})$ space). Medium to high pitch-angle electrons with
<figure><img src="image_8.png" /><figcaption>**Fig. 8** Probability distribution of $a > 1$ (i.e., of nonlinear wave-particle interaction) for electrons and intense field-aligned whistler-mode waves in the outer radiation belt, for L-shell $\epsilon \in [5, 7]$. Dataset is from [Zhang, Mourenas, Artemyev, Angelopoulos, Bortnik, Thorne, Kurth, Kletzing and Hospodarsky, 2019].</figcaption></figure>

not-too-high energies interact resonantly with waves close to the equator (where $\partial\Omega_0/\partial s$ is small), and this increases the probability of nonlinear interaction. Small pitch-angles and/or very high energy electrons interact resonantly with waves at high latitudes (farther from the equator), where the large $\partial\Omega_0/\partial s$ reduces the probability of nonlinear interaction. Note that the equation for $a$ has been derived from Hamiltonian (6), whereas this Hamiltonian would not work for very small pitch-angle particles having $I_x\Omega_0 \sim eU_w$ (see details in Section 3.1 and Appendix B). In this limit of small pitch-angles, another criterion for nonlinear interactions can be derived [Gan et al., 2024], which shows the possibility for such interactions.

## 2.2 Diffusion by monochromatic waves

To demonstrate how the inhomogeneity can affect wave-particle interactions, let us derive the electron diffusion rates in the simple case of a field-aligned whistler mode wave. Equations of motion for Hamiltonian (2) have the form:

$$ \mathcal{H} = mc^2\gamma - \sqrt{\frac{2I_x\Omega_0}{m_e c^2}} \frac{eA_w}{\gamma} \cos(\phi + \psi), \quad \gamma = \sqrt{1 + \left(\frac{p_\|}{m_e c}\right)^2 + \frac{2I_x\Omega_0(z)}{m_e c^2}} $$

and we can use $U_w = \sqrt{2I_x\Omega_0/m_e c^2}A_w/\gamma$.

Let us consider the approximation of a small wave amplitude $|A| \gg |B|$ at the resonance. Then, we can consider unperturbed ($U_w = 0$) particle trajectories to determine the phase $\zeta$ variation around the resonance $t = t_R$, $\zeta_R = \zeta(t_R)$. Expanding $\zeta$ around the resonance, we obtain

$$ \zeta \approx \zeta_R + \frac{1}{2}(t - t_R)^2 \ddot{\zeta}|_{t=t_R} $$

because $\dot{\zeta}|_{t=t_R} = 0$. For Hamiltonian $\mathcal{H}_\zeta$ from Eqs. (6) $\ddot{\zeta}|_{t=t_R}$ can be expressed as

$$ \ddot{\zeta} = \frac{\dot{P}_\zeta}{M} = \frac{\{A, I_R\}}{M} + \frac{U_w}{M} \cos \zeta $$
Thus, for $|A| \gg |B|$ (i.e., $|A, I_R| \gg |eU_w|$) we have

$$ \ddot{\zeta}|_{t=t_R} = \frac{A}{M} = \frac{\{A, I_R\}}{M} $$

Therefore, the $\Delta I$ change in the resonance is given by the time integral of $-\partial H_I / \partial \zeta$ and can be written as

$$ \begin{aligned} \Delta I &= eU_w \int_{-\infty}^{+\infty} \sin \zeta dt = eU_w \int_{-\infty}^{+\infty} \sin \left( \zeta_R + \frac{1}{2} \ddot{\zeta} \bigg|_{t=t_R} (t-t_R)^2 \right) dt \\ &= 2eU_w \left( \frac{1}{2} \ddot{\zeta} \bigg|_{t=t_R} \right)^{-1/2} \left( \sin \zeta_R \int_0^{+\infty} \cos(q^2) dq + \cos \zeta_R \int_0^{+\infty} \sin(q^2) dq \right) \\ &= \sqrt{2\pi} e U_w \left( \frac{1}{2} \ddot{\zeta} \bigg|_{t=t_R} \right)^{-1/2} \sin \left( \zeta_R + \frac{\pi}{4} \right) \end{aligned} \quad (7) $$

where the $\pm\infty$ limits of integration correspond to the large $\Delta t\omega \gg 1$ approximation for the time $\Delta t$ of motion around the resonance, and we use the Fresnel integral equation

$$ \int_{0}^{+\infty} \cos(q^2) dq = \int_{0}^{+\infty} \sin(q^2) dq = \sqrt{\frac{\pi}{8}} $$

The averaged momentum change is zero, $\langle \Delta I \rangle_{\zeta_R} = 0$, and we have

$$ \langle (\Delta I)^2 \rangle_{\zeta_R} = \frac{\pi (eU_w)^2}{\ddot{\zeta}|_{t=t_R}} = \pi \left( \frac{m_e c^2}{\omega} \right)^2 \frac{B_w^2 2I_x \Omega_0 \Omega_0^2}{B_0^2 m_e c^2 \omega^2 N_\parallel^4} \frac{c}{\delta v_R} \quad (8) $$

where $\delta v_R$ is the effective resonance width:

$$ \delta v_R = c \frac{\partial \ln \Omega_0}{k_\parallel \partial s} \left( \left( \frac{p_{\parallel,R}}{m_e c} \right)^2 \left( \frac{\partial \ln N_\parallel}{\partial \ln \Omega_0} + \frac{1}{2} \right) + \frac{\Omega_0}{\omega N_\parallel} \frac{p_{\parallel,R}}{m_e c} + \frac{1}{2} (1 - \gamma^2) \right). $$

In Eq. (8), all variables depending on $s$ should be evaluated at the resonance location.

Taking into account Eq. (4), the equation for $\langle (\Delta \gamma)^2 \rangle$ can be rewritten as

$$ \langle (\Delta \gamma)^2 \rangle_{\zeta_R} = \left\langle \left( \frac{\omega \Delta I}{mc^2} \right)^2 \right\rangle_{\zeta_R} = \pi \frac{B_w^2 2I_x \Omega_0 \Omega_0^2 c}{B_0^2 m_e c^2 \omega^2 N_\parallel^4} $$

The conservation of magnetic moment $I_x = (\gamma^2 - 1) \sin^2 \alpha_{eq} / 2\Omega_0(0)$ provides an additional relation between energy $mc^2\gamma$ and pitch-angle $\alpha_{eq}$ changes:

$$ \left. \frac{\partial \alpha_{eq}}{\partial \gamma} \right|_{I_x = \text{const}} = \frac{\gamma \tan \alpha_{eq}}{\gamma^2 - 1} $$
This relation allows the recalculation of $\langle (\Delta\gamma)^2 \rangle_{\zeta_R}$ to $\langle (\Delta\alpha_{eq})^2 \rangle_{\zeta_R}$. Therefore, we may characterize wave-particle interactions by diffusion coefficients

$$D_{\gamma\gamma} = \frac{\langle (\Delta\gamma)^2 \rangle_{\zeta_R}}{2\tau}, \quad D_{\alpha_{eq}\alpha_{eq}} = \frac{\langle (\Delta\alpha_{eq})^2 \rangle_{\zeta_R}}{2\tau} = \left( \frac{\gamma \tan \alpha_{eq}}{\gamma^2 - 1} \right)^2 D_{\gamma\gamma} \quad (9)$$

where $\tau = \tau(\gamma, \alpha_{eq})$ is the time between two resonant interactions for a single resonance within half of a bounce period, $\tau = \tau_{bounce}/2$, and Eqs. (9) provide the bounce averaged energy and pitch-angle diffusion rates [see also Albert, 2010].

Note that Eq. (7) provides the $\Delta I$ change for unperturbed particle trajectories (i.e., when particle coordinate and velocity do not depend on $eU_w$) and, thus, $\Delta I \sim \int eU_w \sin\zeta dt$ with $\zeta = \zeta_0(t)$ does not depend on $eU_w$. The mean value of such a change is $\langle \Delta I \rangle = 0$, whereas its variance is $\langle (\Delta I)^2 \rangle \propto (eU_w)^2$. The next order of particle trajectory perturbations should include terms linearly proportional to $eU_w$, and thus $\zeta = \zeta_0(t) + eU_w \cdot C \sin(\zeta_0(t))$ where $C$ is a function of slow coordinates. At this next order, $\Delta I \sim \int eU_w \sin\zeta dt$ will contain two terms $\int eU_w \sin\zeta_0 dt$ and $\sim \int eU_w^2 \sin^2\zeta dt$, with a mean value $\langle \Delta I \rangle \propto (eU_w)^2$ provided by the second term. Thus, to estimate $\langle \Delta I \rangle$ of the same order as $\langle (\Delta I)^2 \rangle \propto (eU_w)^2$, one needs to consider a second order perturbation of particle trajectories. The absence of a mean value $\langle \Delta I \rangle$ for the unperturbed trajectories does not mean that $\langle \Delta I \rangle = 0$, but only implies that one has not $\langle \Delta I \rangle \propto eU_w$, but instead $\langle \Delta I \rangle \propto (eU_w)^2$.

### 3 Nonlinear resonant characteristics

Let us describe the main characteristics of nonlinear resonant interactions. We start with the resonant energy change experienced by the particles crossing the resonance $\dot{\zeta} = 0$, the so-called transient particles. Using Eq. (4), we can write $m_e c^2 \Delta\gamma = \omega \Delta I$ where $\Delta I$ is given by integration of the Hamiltonian equation for Hamiltonian (3):

$$\Delta I = 2eU_w \int_{-\infty}^{t_R} \sin \zeta dt$$

where $t_R$ is the time of resonance crossing ($\dot{\zeta} = 0$). Using Hamiltonian (6) for $dt = d\zeta/\dot{\zeta} = Md\zeta/P_\zeta$, we obtain

$$\Delta I = 2eU_w M \int_{-\infty}^{\zeta_R} \frac{\sin \zeta}{P_\zeta} d\zeta = \int_{-\infty}^{\zeta_R} \frac{eU_w \sqrt{2M} \sin \zeta d\zeta}{\sqrt{2\pi\xi A - A\zeta - B \cos \zeta}}$$

where we introduce $\xi = \mathcal{H}_\zeta/2\pi A$, the energy of particles in the phase portrait from Fig. 7. Equation (10) can be written as $\Delta I = \Delta I_0 F(\xi, a)$, where

$$\Delta I_0 = \sqrt{2M |A|}$$

and

$$F = \int_{-\infty}^{\zeta_R} \frac{a \sin \zeta d\zeta}{\sqrt{2\pi\xi - \zeta - a \cos \zeta}} \quad (10)$$
The function F is shown in Fig. 9(a). It is a periodic function of ξ with the period 1:

$$F(a, \xi + 1) = \int_{-\infty}^{\zeta_R} d\zeta \frac{\sqrt{a \sin \zeta}}{\sqrt{2\pi\xi + 2\pi - \zeta - a \cos \zeta}} = \int_{-\infty}^{\zeta_R - 2\pi} d\tilde{\zeta} \frac{\sqrt{a \sin \tilde{\zeta}}}{\sqrt{2\pi\tilde{\xi} - \tilde{\zeta} - a \cos \tilde{\zeta}}} = \\ = \int_{-\infty}^{\tilde{\zeta}_R} d\tilde{\zeta} \frac{\sqrt{a \sin \tilde{\zeta}}}{\sqrt{2\pi\tilde{\xi} - \tilde{\zeta} - a \cos \tilde{\zeta}}} = F(a, \xi)$$

For $a < 1$, this periodic function is such that its average $\langle F(a, \xi) \rangle_{\xi \in [0,1]}$ is zero, whereas for $a > 1$ the profile of F becomes asymmetric relative to zero and $\langle f(a, \xi) \rangle_{\xi \in [0,1]}$ is finite (see Appendix A and [Neishtadt, 1975]). As function F takes both positive and negative values, resonant electrons can increase and decrease their moment I, with $\Delta I \sim F$. The $\xi$-averaged F is always non positive, and this determines the conventional description of electron resonant scattering (for $a > 1$ and the first cyclotron resonance this scattering is called phase bunching, see [Matsumoto and Omura, 1981; Omura and Matsumoto, 1982; Winglee, 1985]) as a process with electron momentum (and energy) decrease [e.g., Albert, 2001]. However, for individual $\xi$ values (i.e., for specific values of electron resonant energy $H_\zeta$) $\Delta I$ is positive, and this effect is called positive phase bunching and has been considered in [Albert et al., 2022; Vargas et al., 2023]. Although such positive phase bunching may be important for transient electron scattering by very intense whistler-mode wave packets [see discussion in Inan et al., 1978; Lundin and Shkliar, 1977], the overall electron ensemble dynamics can be described by $\xi$-averaged system characteristics that do not include positive phase bunching [Vargas et al., 2023].

The energy $\xi$ depends on initial particle gyrophase and can be considered as a random variable with uniform distribution within $[0, 1]$ (see Fig. 8 in [Itin et al., 2000] and Fig. 5 in [Frantsuzov et al., 2023]). Therefore, to estimate the actual energy variation of transient particles we shall average the function F over $\xi$ [see also Artemyev, Neishtadt, Vainchtein, Vasiliev, Vasko and Zelenyi, 2018; Frantsuzov et al., 2023; Neishtadt, 2014]:

$$\langle F \rangle_{\xi} = -\frac{1}{\pi} \int_{\zeta_{-}}^{\zeta_{+}} \sqrt{\zeta_{+} - \zeta + a(\cos \zeta_{+} - \cos \zeta)} d\zeta \quad (11)$$

This equation has been derived in [Neishtadt, 1975], and we repeat the derivations in Appendix A.

Equation (11) shows that $\langle \Delta I \rangle_{\xi} = \Delta I_0 \langle F \rangle_{\xi}$ can be written as

$$\langle \Delta I \rangle_{\xi} = \sqrt{2M|A|} \langle F \rangle_{\xi} \quad (12)$$

and is equal to $-S/2\pi$, where S is the area surrounded by the separatrix in Fig. 7(b):

$$S = \sqrt{8M|A|} \int_{\zeta_{-}}^{\zeta_{+}} \sqrt{\zeta_{+} - \zeta + a(\cos \zeta_{+} - \cos \zeta)} d\zeta, \quad (13)$$
<figure><img src="image_9.png" /><figcaption>**Fig. 9** Panel (a) shows function $F(\xi, a)$ given by Eq. (10). Panel (b) shows functions $\langle F(\xi, a) \rangle_{\xi}$, $\langle F^2(\xi, a) \rangle_{\xi}^{1/2}$, and $\mathrm{Var}(F) = \sqrt{\langle F(\xi, a)^2 \rangle_{\xi} - \langle F(\xi, a) \rangle_{\xi}^2}$. Blue dashed lines present fittings of the displayed functions.</figcaption></figure>

Here the two $\zeta_{\pm}$ values are the coordinate of the saddle point ($\sin \zeta_{-} = 1/a$) and a solution of $\zeta + a \cos \zeta = \zeta_{-} + a \cos \zeta_{-}$ equation, different from $\zeta_{-}$. From the definition of $S$ it is clear that $S = 0$ for $a \le 1$. The equality $\langle \Delta I \rangle_{\xi} = -S/2\pi$ is an important property of the Hamiltonian system (6) that determines a balance between phase trapping and phase bunching processes [see Artemyev, Neishtadt, Vasiliev and Mourenas, 2016; Itin et al., 2000; Solovev and Shkliar, 1986]. A useful asymptotic expression of $S$ is

$$S_{a \gg 1} \approx \sqrt{8M|B|} \int_{0}^{2\pi} \sqrt{1 - \cos \zeta} d\zeta = 2^{5/2} \sqrt{8M|B|} = 16\sqrt{M|B|}$$

where $\zeta_{-} \to 0$ and $\zeta_{+} \to 2\pi$. This asymptote shows that $S$ scales with wave amplitude as $S \propto \sqrt{|B|} \propto \sqrt{U_w} \propto \sqrt{B_w/B_0}$.

Figure 9(b) shows functions $\langle F(\xi, a) \rangle_{\xi}$, $\langle F^2(\xi, a) \rangle_{\xi}^{1/2}$, and $\sqrt{\langle F(\xi, a)^2 \rangle_{\xi} - \langle F(\xi, a) \rangle_{\xi}^2}$ that describe the mean energy change, variance, and energy dispersion. Comparing Figs. 9(a) and (b), we see that although $F$ can be both positive and negative, the averaged $\langle F(\xi, a) \rangle_{\xi}$ is always positive (or zero, for $a < 1$). Thus, individual particles may experience $\Delta I$ increase due to the phase bunching, but the $\xi$-averaged effect of such bunching corresponds to $\Delta I < 0$. Such positive phase bunching ($\Delta I > 0$) has been explained and investigated in Albert et al. [2022], whereas Vargas et al. [2023] demonstrated that positive bunching does change the overall dynamics of a charged particle ensemble, which can be fully described by the averaged $\langle F(\xi, a) \rangle_{\xi}$ function.

Let us now consider particles experiencing phase trapping. In contrast with transient particles merely crossing the resonance $\dot{\zeta} = 0$, trapped particles cross the separatrix (the curve demarcating open and closed trajectories in Fig. 7(b)), and their characteristic motion is qualitatively changed, i.e., they start oscillating around the resonance $\dot{\zeta} = 0$. To estimate the amount of such particles, we shall
compare the variation of the area surrounded by the separatrix, $S$, and the total phase space flux $\mathcal{F}$ crossing the resonance:

$$ \mathcal{F} = \int_{0}^{2\pi} \dot{P}_{\zeta} d\zeta = 2\pi A $$

This comparison provides the so-called probability of trapping:

$$ \Pi \approx \frac{\{S, A\}}{|\mathcal{F}| + \{S, A\}/2} = \frac{2\{S, A\}}{4\pi|\{I_R, A\}| + \{S, A\}} \quad (14) $$

where $\dot{S} = \{S, A\}$ and $\Pi = 0$ for area decrease, $\dot{S} < 0$. If $S$ changes slowly enough ($\dot{S} < 4\pi|\{I_R, A\}|$), which is the case in almost all systems under consideration, we can use the approximation

$$ \Pi \approx \frac{\{S, A\}}{2\pi|\{I_R, A\}|} \quad (15) $$

We note that $I_R, A = \dot{I}_R$, and thus the probability of trapping can be rewritten as

$$ \Pi \approx \frac{\dot{S}}{2\pi \dot{I}_R} = \frac{1}{2\pi} \frac{dS}{dI} \Big|_{I=I_R} $$

Taking into account that $\langle \Delta I \rangle_{\xi} = -S/2\pi$, we obtain

$$ \Pi \approx -\frac{d \langle \Delta I \rangle_{\xi}}{dI} $$

This equation provides a direct relationship between the probability of trapping and the average variation $\langle \Delta I \rangle_{\xi}$ due to nonlinear scattering.

The probability of trapping can be understood as the ratio of particles experiencing trapping for a single resonant interaction to the total number of particles crossing the resonance. For fixed system characteristics and $h$ constant, this probability depends only on the initial electron energy. Figure 10 shows examples of verification of Eq. (15) for the first cyclotron resonance with field-aligned waves and for the Landau resonance with very oblique waves [see more examples in Artemyev et al., 2012; Artemyev, Vasiliev, Mourenas, Agapitov, Krasnoselskikh, Boscher and Rolland, 2014; Artemyev, Vasiliev, Mourenas, Neishtadt, Agapitov and Krasnoselskikh, 2015; Vainchtein et al., 2018]. The usual scheme for an evaluation of the probability of trapping in numerical test particle simulations includes an integration of trajectories for a large particle ensemble with the same initial energy and pitch-angles, but random initial gyrophases. This ensemble passes through the resonance once and we can count the number of particles trapped into the resonance. The ratio of such trapped particles to the initial total number of particles in the ensemble will provide an estimate of the probability of trapping. Each colored circle in Fig. 10 has been obtained via such a scheme, applied to different particle (energy, pitch-angle) and system (L-shell, wave amplitude) parameters. This numerically evaluated probability of phase trapping is compared with the analytical formula of $\Pi$, confirming that Eq. (15) accurately describes the trapping probability. Note that we can resort to probabilities for the description of trapping,
<figure><img src="image_10.png" /><figcaption>**Fig. 10** Four sets of examples of numerical verification of the equation for the probability of trapping: curves show theoretical results from Eq. (15) and circles show results of test-particle simulations (a relative number of particles experiencing the phase trapping for a single resonant interaction). Panel (a) shows results for Landau resonance from [Artemyev, Vasiliev, Mourenas, Agapitov and Krasnoselskikh, 2013] (wave amplitudes are the same for black and blue curves, but a factor of 2 larger for the red curve), panel (b) shows results for Landau resonance from [Artemyev, Vasiliev, Mourenas, Agapitov, Krasnoselskikh, Boscher and Rolland, 2014], panel (c) shows results for the cyclotron resonance from [Artemyev, Vasiliev, Mourenas, Neishtadt, Agapitov and Krasnoselskikh, 2015], and panel (d) shows results for the cyclotron resonance from [Vainchtein et al., 2018]. Details of wave models and background magnetic field conditions can be found in the corresponding studies.</figcaption></figure>

because in any realistic system the initial electron gyrophase is an unknown parameter. Therefore, we can average the system over this gyrophase to reduce its dimensionality, which leads to an inherently probabilistic description of trapping [Neishtadt, 1975; Shklyar, 1981].

To evaluate the energy variation of an electron due to phase trapping, let us consider the motion of a trapped electrons. After crossing the separatrix in the phase portrait from Fig. 7(b), electrons start rotating around the resonance along closed trajectories. This rotation occurs with a frequency $\Omega_{tr} \sim \sqrt{B/M} \propto (B_w/B_0)^{1/2} \Omega_0$, whereas the phase portrait evolves with the rate of (s, p) change, that is $\sim 2\pi/\tau_b \gg \Omega_0$ (the bounce period being the longest time scale in this system). Note that the resonance condition for whistler-mode waves assumes that $\Omega_0$
<figure><img src="image_11.png" /><figcaption>**Fig. 11** Schematic view of electron trapping/de-trapping (escape) for the Landau (a) and cyclotron (b) resonances. Top panels show typical profiles of area $S(\lambda)$, and bottom panels show unperturbed electron trajectory $p_\parallel(\lambda)$ and resonant condition $p_\parallel = p_R(\lambda)$. Trapping and escape positions are shown by vertical red lines.</figcaption></figure>

is of the same order as $k_\parallel \dot{s} \sim k_\parallel R/\tau_b$, where R is the spatial scale of the background magnetic field inhomogeneity. On the other hand, the condition for nonlinear resonant interaction, $a > 1$, assumes that $k_\parallel^{-1} \partial \ln \Omega_0 / \partial s \sim 1/k_\parallel R$ is of the same order as the wave strength $\sim B_w/B_0$, i.e. $k_\parallel R \sim B_0/B_w$ and $\Omega_0 \tau_b \sim B_0/B_w \gg 1$. Thus, trapped electron oscillations around the resonance are much faster than the phase portrait evolution, $\Omega_{tr} \tau_b \propto (B_w/B_0)^{1/2} \Omega_0 \tau_b \propto (B_w/B_0)^{-1/2} \gg 1$. Such a fast periodical motion should introduce an adiabatic invariant $I_\zeta = (2\pi)^{-1} \oint P_\zeta d\zeta$ [Landau and Lifshitz, 1960], which is equal to $S_{tr}/2\pi$ evaluated at the time of trapping. Thus, trapped electrons move within the region surrounded by the separatrix, and at the time of trapping their invariant (area surrounded by their trajectories) is $I_\zeta = S_{tr}/2\pi$ with an increasing $S_{tr}$. Electrons will stay trapped until $S/2\pi$ becomes larger than $I_\zeta$. Accordingly, electrons will escape from the trapping regime when $S = S_{tr}$ and $S$ decreases. Illustration of this trapping/de-trapping dynamics is shown in Fig. 11.

Figures 9-11 show that all nonlinear resonant effects are described by the $S(I)$ (or $S(\gamma)$) function: the energy variation of phase bunched particles is $\Delta\gamma = -\omega S/2\pi m_e c^2$, the probability of phase trapping is $\Pi \approx (2\pi)^{-1} dS/dI = (2\pi)^{-1}(dS/d\gamma)$. $(\omega/m_e c^2)$, and the energy variation due to trapping is determined by the equation $S(\gamma) = S(\gamma + \Delta\gamma_{tr})$ [see also Artemyev, Neishtadt, Vasiliev and Mourenas, 2018]. Conversely, electron diffusion requires knowing $\langle F^2 \rangle_\xi$, which cannot be expressed through the $S = -2\pi\sqrt{2|A|\bar{M}\langle F \rangle_\xi}$ function, and should be evaluated separately. In the next section, we will use a quite universal property of $S$ to construct a kinetic equation including the effects of nonlinear resonant interactions.
## 3.1 Small pitch-angle limit

There is one important limitation of system (6, 5): the factor ~ $B_w \sqrt{2I\Omega_0}$ in the wave term is implicitly assumed to be constant within a typical time-scale of resonant interaction ~ $\Omega_{tr}^{-1}$. This assumption is valid as long as the phase variation $\dot{\zeta}$ is controlled by the $O(B_w/B_0)$ term, while the term ~ $B_w \sqrt{2\Omega_0/I} \sin \zeta$ remains a small correction. This assumption is naturally violated for sufficiently small $I$, when the term ~ $B_w \sqrt{2\Omega_0/I}$ becomes important and the variation rate of $I$ becomes comparable to the $\zeta$ variation rate [Lundin and Shkliar, 1977]. A detailed description of this case can be found in [Albert et al., 2021; Artemyev, Neishtadt, Albert, Gan, Li and Ma, 2021] and in Appendix B, whereas multiple important effects of small $I$ (small $I_x$) on resonant electron motion are described in [Gan et al., 2022; Grach and Demekhov, 2020; Kitahara and Katoh, 2019]. Here, we only briefly discuss these effects.

Let us consider resonant nonlinear scattering of electrons with small $I_x = I$ and $n_r = -1$, and expand the Hamiltonian (3) for small $I$:

$$ \tilde{\mathcal{H}}_I \approx \Lambda_0 + \frac{1}{2M_0} (I - I_R)^2 + \sqrt{\frac{2I\Omega_0}{m_ec^2}} \frac{eB_w}{k} \cos \zeta $$

where $M_0 = \partial^2 \mathcal{H}_I / \partial I^2|_{I=0}$, $\Lambda_0 = \mathcal{H}_I|_{I=0}$, and $I_R$ is given in Appendix B. The main difference from the Hamiltonian expanded around the resonant $I_R$ is that wave amplitude in $\tilde{\mathcal{H}}_I$ depends on $I$. The Hamiltonian $\tilde{\mathcal{H}}_I$ describes fast $\zeta$ motion and slow $(s, p, I)$. An important property of this Hamiltonian is that for $I \sim (B_w/B_0)^{2/3}$ the dynamics of $I$ becomes as fast as the dynamics of $\zeta$ and, thus, there is no time separation between $\zeta$ and $I$, which become fast variables. In Appendix B we show that the Hamiltonian of these fast variables is

$$ \mathcal{L} = \frac{1}{2} \left( \frac{1}{2} P^2 + \frac{1}{2} q^2 - Y_R \right)^2 + uq \quad (16) $$

where $P = (B_w/B_0)^{-1/3} \sqrt{2I\Omega_0/m_ec^2} \sin \zeta$, $q = -(B_w/B_0)^{-1/3} \sqrt{2I\Omega_0/m_ec^2} \cos \zeta$, $Y_R = (B_w/B_0)^{-2/3} I_R \Omega_0 / m_e c^2$, and $u$ is a normalization constant. The phase portrait of Hamiltonian (16) is shown in Fig. 12(c). Let us compare this phase portrait with phase portraits of the Hamiltonian expanded around $I_R$,

$$ \mathcal{H}_I = -\omega I_R + \gamma_R + \frac{1}{2M} (I - I_R)^2 + \sqrt{\frac{2I_R \Omega_0}{m_e c^2}} \frac{eB_w}{k} \sin \zeta \quad (17) $$

The main difference between this Hamiltonian and the Hamiltonian from Eq. (16) is that the effective wave amplitude does not depend on $I$. The phase portrait of Hamiltonian (17) with frozen slow variables is shown in Fig. 12(a).

Instead of introducing $(q, P)$ coordinates, it is more convenient to introduce $P_\zeta = I - I_R$ and rewrite Hamiltonian (17) into Eq. (6). The phase portrait of Hamiltonian (6) with frozen slow variables is shown in Fig. 12(b). This is the classical portrait of the pendulum with torque with three main phase space regions: before resonance $P_\zeta = 0$ crossing, particles are in $G_1$, and resonance crossing can result in trapping (particles appear in $G_0$) or scattering (particles appear in $G_2$). Therefore, there is a direct relation between three regions $G_{0,1,2}$ of Hamiltonian of (16) and Hamiltonian (6).
<figure><img src="image_12.png" /><figcaption>Fig. 12 Comparison of phase portraits of (a) Hamiltonian (3), (b) Hamiltonian (6), and (c) Hamiltonian (16).</figcaption></figure>

For the initial system given by Eq. (2) with $n_r = -1$ the phase bunching (transition from $G_1$ to $G_2$) always appears with $I_x$ decrease, and this effect is well seen in the phase portrait (c) of Fig. 12: the area $S_1$ is always larger than the area $S_2$. But when the initial invariant $I_x$ is sufficiently small, particles become trapped within region $G_0$ as soon as this region appears during particle motion along their trajectories. In phase portrait (c) of Fig. 12 this trapping means that the area surrounded by the particle trajectory $2\pi \oint Pdq \sim 2\pi\sqrt{2I\Omega_0/m_ec^2}$ is smaller than $S_0$ at the moment when $G_0$ appears. The threshold $I_x$ value is $2I_x\Omega_0/m_ec^2 \sim (B_w/B_0)^{2/3}$ (see Appendix B and [Albert et al., 2021; Artemyev, Neishtadt, Albert, Gan, Li and Ma, 2021]).

This is the so-called autoresonance phenomena of 100% probability of trapping in resonant systems [see, e.g., Fajans and Frièdland, 2001; Friedland, 2009; Neishtadt, 1975; Neishtadt et al., 2013; Sinclair, 1972, and references therein]. Being trapped into $G_0$, particles can both increase or decrease their $I_x$ during the transition from $G_0$ to $G_2$: the $I_x$ change depends on the ratio of $S_0/S_2$ at the moment when $S_0$ becomes equal to $2\pi \oint Pdq$. For sufficiently small $2\pi \oint Pdq$ (small $I_x$), the ratio $S_0/S_2 = 2\pi \oint Pdq/S_2$ will be below one, and particles will increase their $I_x$ due to the resonant interaction. Formally, this interaction cannot be called bunching, because particles are trapped into $G_0$ from the beginning.

## 4 Kinetic equation with nonlinear resonant effects

In this section, we use the function of $S(\gamma)$ to derive a kinetic equation for the distribution function of electrons nonlinearly interacting resonantly with whistler-mode waves for a constant $h$ given by Eq. (4). Let us start with the simplified situation of a constant bounce period $\partial\tau_b/\partial\gamma = 0$. We separate the $\hat{K}$ kernel into two parts: one part describing diffusion and the phase bunching effect $\hat{K}_{bc}$ and another part describing the trapping effect $\hat{K}_{tr}$. Using the approximation of small energy change due to bunching, we expand $\hat{K}_{bc}$ as

$$\int_0^\infty \hat{K}_{bc}(\gamma, \gamma') f(\gamma') d\gamma' = \frac{\partial}{\partial\gamma} \left( D_{\gamma\gamma} \frac{\partial f}{\partial\gamma} \right) - \frac{\partial}{\partial\gamma} (V_\gamma f)$$
Note that there are two main contributions to the particle drift term $V_γ$: the gradient of the diffusion rate, $V_γ \propto \partial D_{γγ}/\partial \gamma \propto (B_w/B_0)^2$, and nonlinear scattering (phase bunching effect), $V_γ \propto S \propto (B_w/B_0)^{1/2}$. The first term determines the divergence-free condition for the Fokker-Planck diffusion equation and is taken into account in the form of a $\propto D_{γγ}$ operator [Lieberman and Lichtenberg, 1973; Sagdeev et al., 1988]. Thus, in the equation for $\hat{K}_{bc}$ we may include only the second term $V_γ \propto (B_w/B_0)^{1/2}$, which is much larger than the first one, because $B_w/B_0 \ll 1$.

For the trapping kernel, we use the following formulation:

$$ \int_{0}^{\infty} \hat{K}_{tr+} (\gamma | \gamma') f(\gamma') d\gamma' - \int_{0}^{\infty} \hat{K}_{tr-} (\gamma' | \gamma) f(\gamma) d\gamma' $$

where the second term describes particle transport from the energy of trapping, $\gamma$, and the first term describes particle transport toward the energy of escape from the resonance, $\gamma$, from the energy of trapping, $\gamma'$, given by $S(\gamma') = S(\gamma)$. Substituting $K_{bn}$ and $K_{tr}$ into the Smoluchowski equation (1), we obtain

$$ \frac{\partial f}{\partial t} = \frac{\partial}{\partial \gamma} \left( D_{\gamma\gamma} \frac{\partial f}{\partial \gamma} \right) - \frac{\partial}{\partial \gamma} (V_{\gamma} f) + \int_{0}^{\infty} \hat{K}_{tr+} (\gamma | \gamma') f(\gamma') d\gamma' - \int_{0}^{\infty} \hat{K}_{tr-} (\gamma' | \gamma) f(\gamma) d\gamma' $$

The term responsible for particle transport from the energy of trapping, $\gamma$, can be written as

$$ \int_{0}^{\infty} \hat{K}_{tr-} (\gamma' | \gamma) f(\gamma) d\gamma' = \int_{0}^{\infty} \frac{\Pi(\gamma)}{\tau_{b}} \delta(\gamma - \gamma') f(\gamma) d\gamma' = \frac{\Pi(\gamma)}{\tau_{b}} f(\gamma) (1 - \Theta) $$

where

$$ \Theta(\gamma) = \begin{cases} 0, & dS/d\gamma < 0 \\ 1, & dS/d\gamma > 0 \end{cases} \qquad (18) $$

The term responsible for particle transport from the energy $\gamma'$ into the energy of escape, $\gamma$, can be written as

$$ \int_{0}^{\infty} \hat{K}_{tr+} (\gamma | \gamma') f(\gamma') d\gamma' = \int_{0}^{\infty} \Pi(\gamma') \delta(\gamma - \gamma') f(\gamma') d\gamma' \\ = \frac{\Pi(\gamma^*)}{\tau_b} \left| \frac{dT(\gamma^*)}{d\gamma^*} \right|^{-1} f(\gamma^*) $$

where $T(\gamma^*) = \gamma$ is the solution of the equation $S(\gamma^*) = S(\gamma)$. Note that Eq. (4) provides a linear relation between electron energy $\gamma$ and momentum $I$: $d\gamma/dI = \omega/m_e c^2 = \text{const}$. This relation allows us to use the equation $d\gamma^*/d\gamma = dI^*/dI$ to write:

$$ \frac{dT(\gamma^*)}{d\gamma^*} = \frac{d\gamma}{d\gamma^*} = \frac{dI}{dI^*} = \frac{dS(I^*)/dI^*}{dS(I)/dI} = -\frac{\Pi(I^*)/\tau_b}{dS(I)/dI} = -\frac{\Pi(\gamma^*)/\tau_b}{dV_\gamma(\gamma)/d\gamma} $$
Thus, for $\hat{K}_{tr+}$ we obtain for $\Theta = 1$

$$ \int_{0}^{\infty} K_{tr+} (\gamma | \gamma') f(\gamma') d\gamma' = \frac{dS(\gamma)}{d\gamma} f(\gamma^*) $$

Therefore, combining all these terms in the equation for $\partial f / \partial t$, we obtain:

$$ \frac{\partial f}{\partial t} = \frac{\partial}{\partial \gamma} \left( D_{\gamma\gamma} \frac{\partial f}{\partial \gamma} \right) - \frac{\partial}{\partial \gamma} (V_{\gamma} f) + \frac{dV_{\gamma}(\gamma)}{d\gamma} f(\gamma^*) \Theta - \frac{\Pi(\gamma)}{\tau_b} f(\gamma) (1 - \Theta) \quad (19) $$

Taking into account that $\Pi = (2\pi)^{-1}(dS/dI) = -dV_{\gamma}/d\gamma$, we can rewrite Eq. (19) as:

$$ \frac{\partial f}{\partial t} = \frac{\partial}{\partial \gamma} \left( D_{\gamma\gamma} \frac{\partial f}{\partial \gamma} \right) - V_{\gamma} \frac{\partial f}{\partial \gamma} + \frac{dV(\gamma)}{d\gamma} \left( f(\gamma^*) - f(\gamma) \right) \Theta $$

This is a basic formulation of a kinetic equation describing electron distribution dynamics due to multiple nonlinear resonant interactions. The numerical solution of this equation has been verified by comparisons with test particle simulations in [Artemyev, Neishtadt, Vasiliev and Mourenas, 2016; Artemyev et al., 2017b; Leoncini et al., 2018]. Figure 13 shows two examples of such verifications. We use two different initial conditions for the $f(\gamma)$ function and solve Eq. (19) within the energy range of $S \neq 0$, i.e., we consider electron dynamics only within the energy range of nonlinear resonant interactions, $\gamma \in [\gamma_{min}, \gamma_{max}]$. The diffusion rate $D_{\gamma\gamma}$ does not necessarily vanish at $\gamma_{min, max}$ and particles can diffusively move in/out of the energy range of nonlinear resonant interactions. Thus, to verify Eq. (19) within the range $\gamma \in [\gamma_{min}, \gamma_{max}]$ we numerically modify $D_{\gamma\gamma} \rightarrow D_{\gamma\gamma} \cdot 4(\gamma - \gamma_{min})^2 (\gamma_{max} - \gamma)^2 / (\gamma_{max} - \gamma_{min})^4$, to guarantee that electrons will not leave the range $\gamma \in [\gamma_{min}, \gamma_{max}]$ of nonlinear resonant interactions. We also numerically integrate $10^6$ test particle trajectories described by the Hamiltonian (2).

The solutions shown in Fig. 13(a) start with a power law distribution $f \sim \gamma^{-3}$. Phase trapping forms an accelerated population of electrons at large energies, and with time this population drifts to smaller energies due to the phase bunching. The fine balance of bunching speed and trapping probability prevents an accumulation of electrons in the region where electrons are released from trapping acceleration: the new accelerated population only replaces a previously accelerated population moved to smaller energy, and the electron population at large energy does not increase in magnitude but occupies a larger energy range as time goes on. At the time $50R/c$, corresponding to about fifty bounce periods with a single resonant interaction during each period, we compare the solution of Eq. (19) with test particle simulation results: the corresponding red and black curves are quite close, demonstrating that the kinetic equation describes well the dynamics of the electron energy distribution.

The solutions shown in Fig. 13(b) start with a distribution having a small maximum at intermediate energies and a plateau. The resonant dynamics of the electron energy distribution includes the same components shown in Fig. 13(a): formation of an accelerated population due to phase trapping and the following drift of this population to smaller energies due to phase bunching. This dynamics is supplemented by the evolution of the maximum initially present at intermediate energies: due to an absence of phase trapping at these energies, this maximum
<figure><img src="image_13.png" /><figcaption>**Fig. 13** Solutions of Eq. (19) for two initial conditions (shown in blue). The function $f$ is normalized to its initial value at the left boundary, $f_b$. Results are shown at three different times, and the solution for $tc/R = 50$ is compared with results of test particle simulations (red) with $10^6$ trajectories. Solutions are obtained for $h/m_ec^2 = 1.454$ (this value of $h$ corresponds to, e.g., an equatorial pitch-angle $\alpha_{eq} = 45^\circ$ for 300 keV electron energy). We use a curvature-free dipole magnetic field [Bell, 1984] at $L$-shell= 4.5. The whistler-mode wave frequency is 0.35 times the electron cyclotron frequency at the equator, and the plasma frequency is 4.5 times the electron cyclotron frequency at the equator. Wave amplitude is 300 pT. The distribution of the wave amplitude along magnetic field lines is modeled by the empirical function $\tanh((\lambda/\delta\lambda_1)^2)\exp(-(\lambda/\delta\lambda_2)^2)$ with $\delta\lambda_1 = 2^\circ$, $\delta\lambda_2 = 20^\circ$. This empirical function fits the typically observed whistler-mode wave intensity distribution [Agapitov et al., 2013].</figcaption></figure>

merely drifts toward smaller energies via phase bunching. The enhanced electron population at smaller energies (due to the initial plateau) provides more particles for trapping acceleration, and the accelerated population in Fig. 13(b) has a larger magnitude in comparison with the population shown in Fig. 13(a). The comparison of solution of Eq. (19) with results of test particle simulations (at time $50R/c$) confirms that the kinetic equation describes correctly the dynamics of the electron energy distribution.

We note that in the case of test particle simulations, we cannot perform the modification of the diffusion coefficient that prevents particle escape from the energy range of nonlinear wave-particle interactions in the numerical solution of Eq. (19). Thus, to perform fair comparisons between test particle simulations and solutions of Eq. (19), we rerun each particle escaping $\gamma \in [\gamma_{min}, \gamma_{max}]$ and save the total number of particles within this energy range. The diffusion $D_{\gamma\gamma}$ changes the electron distribution much slower than nonlinear resonant phase trapping and bunching, and such a weak effect of diffusion may help the comparison of test particle simulation results with solutions of Eq. (19), despite the different descriptions of diffusive scattering in the frame of these two approaches (due to $D_{\gamma\gamma}$ modification in the kinetic equation).
## 4.1 Boundary conditions

The profile $S(\gamma)$ determines all the main terms of the kinetic equation (19) and, therefore, the properties of $S(\gamma)$ are very important for understanding the solution of this equation. The important characteristic of $S(\gamma)$ is its asymptotic behavior near zero values. Depending on the specific distribution of wave electromagnetic field (the term $\sim B_w$ in Eq. (2)0, there are different numbers of zeros of $S(\gamma)$, but the simplest case corresponds to the cyclotron resonance with field-aligned whistler-mode waves when the term $\sim B_w$ varies monotonically along the magnetic latitude. Note that for a conserved invariant $h$ from Eq. (4), the resonant energy $\gamma$ is a linear function of $I = I_x$ and a monotonic function of magnetic latitude $\lambda$ (or coordinate $s$). Figure 14(a) shows examples of latitudinal profiles of resonant energy for different $h$ values.

Near the equatorial plane, $S(\gamma)$ should drop to zero because $B_w \to 0$ in the wave source region, and thus $a \propto B_w$ is less than one [see Shklyar and Luzhkovskiy, 2023; Shklyar, 2017, 2021, for descriptions of alternative $B_w$ profiles for whistler-mode waves generated by lightning in the ionosphere and propagating from high latitudes to the equator]. At high latitudes, the gradient $\partial\Omega_0/\partial s$ becomes quite large and thus $a \propto (\partial\Omega_0/\partial s)^{-1}$ becomes less than one. Thus, for sufficiently large $B_w$ we find that $S$ is above zero between the equatorial plane and some high-latitude location. Let us consider the behavior of $S$ around these two zeros, where $a \to 1$. First, we expand Hamiltonian (6) as [see Artemyev, Neishtadt and Vasiliev, 2019; Artemyev, Neishtadt, Vasiliev and Mourenas, 2021]:

$$ \mathcal{H}_{\zeta} = \frac{1}{2M} P_{\zeta}^2 + B \cdot \left( \frac{1}{a}\zeta + \cos\zeta \right) \approx \frac{1}{2M} P_{\zeta}^2 - B \cdot \left( \left(1 - \frac{1}{a}\right)\zeta - \frac{1}{6}\zeta^3 \right) \quad (20) $$

where we expand $\cos\zeta$ around $\zeta = \pi/2$ and shift $\zeta \to \zeta - \pi/2$. The profile of potential energy $U/B = -(1 - a^{-1})\zeta + \zeta^3/6$ for this Hamiltonian is shown in Fig. 14(b): potential energy has a local minimum in the interval $(\zeta_-, \zeta_+)$. The corresponding area surrounded by the separatrix takes the form

$$ \begin{aligned} S_{a \approx 1} &= \sqrt{8M|B|} \int_{\zeta_{-}}^{\zeta_{+}} \sqrt{(1 - a^{-1})(\zeta - \zeta_{-}) - \frac{1}{6}(\zeta^3 - \zeta_{-}^3)} d\zeta \\ &= \sqrt{8M|B|} (1 - a^{-1})^{5/4} \int_{\zeta_{-}}^{\zeta_{+}} \sqrt{(\tilde{\zeta} - \tilde{\zeta}_{-}) - \frac{1}{6}(\tilde{\zeta}^3 - \tilde{\zeta}_{-}^3)} d\tilde{\zeta} \\ &= \sqrt{8M|B|} (1 - a^{-1})^{5/4} \frac{12 \cdot 3^{3/4}}{5} = \frac{12 \cdot 3^{3/4}}{5} \sqrt{8M \frac{|A|^{5/2}}{|B|^{3/2}} (a - 1)^{5/4}} \end{aligned} $$

where $\tilde{\zeta} = \zeta / (a^{-1} - 1)^{1/2}$. This equation provides the asymptotic variation of $S$ for $a \to 1$, and we can rewrite this asymptote as a function of the resonant energy. Indeed, $a$ is a function of $s$, and for a conserved $h$ given by Eq. (4) there is a profile of resonant energy $\gamma(s)$, see Fig. 14(a). Therefore, we can write $a = a(\gamma)$ and then expand $a$ around the energy $\gamma_b$ corresponding to $a = 1$ and $S = 0$:
<figure><img src="image_14.png" /><figcaption>Fig. 14 (a) Resonant energy ($\gamma - 1)m_e c^2$ as a function of magnetic latitude for first cyclotron resonance with field-aligned waves ($\omega/\Omega_{eq} = 0.35$, $L = 6$) for several $h$ (equatorial pitch-angle $60^\circ$, $45^\circ$, and $30^\circ$ for 200keV electron) values. (b) Profile of $\omega/\Omega_{eq} = 0.35$ for Hamiltonian (20).</figcaption></figure>

$a \approx 1 + C \cdot (\gamma - \gamma_b)$, where $C = \text{const.}$ Substituting this expansion into equation $S_{a \approx 1}$, we obtain

$$S_{a \approx 1} \propto |\gamma - \gamma_b|^{5/4}$$

Note that $S \neq 0$ within some $\gamma \in [\gamma_{min}, \gamma_{max}]$ domain, and we should have an asymptotic form $S \propto |\gamma - \gamma_{min,max}|^{5/4}$ at both boundaries. This asymptotic form of $S$ around $a \approx 1$ shows that $\Pi \propto dS/d\gamma \propto |\gamma - \gamma_b|^{1/4} \to 0$ at the energy boundary, i.e., both drift $V \propto S$ and probability of trapping $\Pi$ drop to zero at the boundary.

The equation for $S_{a \approx 1}$ determines the $S$ profile around the zeros, and can be adopted to describe systems with intrinsically small $S$. The most natural example of such systems is a system with wave-packets. The corresponding schematic is provided in Fig. 15(a). The plane wave $\sim \sin(\zeta)$ is a very simplified approximation of much more complicated and realistic situations, where the wave field often consists of a series of wave-packets [Zhang, Agapitov, Artemyev, Mourenas, Angelopoulos, Kurth, Bonnell and Hospodarsky, 2020; Zhang, Mourenas, Artemyev, Angelopoulos, Bortnik, Thorne, Kurth, Kletzing and Hospodarsky, 2019]. Each packet has a finite length, and small-scale packets will correspond to a small latitudinal range of $S \neq 0$, i.e., the $a \approx 1$ regime will be applicable for the entire range of energies with $S \neq 0$. Then, $S$ can be approximated as

$$S_{small} = S_0 \cdot (\gamma_{max} - \gamma)^{5/4} (\gamma - \gamma_{min})^{5/4}$$

where $S_0$ is the magnitude of $S$ and $\gamma_{max,min}$ are boundary values of energies where $a = 1$. Using changes of variables, $x = 2\gamma/(\gamma_{max} - \gamma_{min}) - (\gamma_{max} + \gamma_{min})/(\gamma_{max} - \gamma_{min})$ and $S_0 \to S_0 \cdot ((\gamma_{max} - \gamma_{min})/2)^{5/2}$, we can write

$$S_{small} = S_0(1 - x^2)^{5/4} \qquad (21)$$
<figure><img src="image_15.png" /><figcaption>Fig. 15 Panel (a) shows a schematic view of S profiles for plane wave (red) and wave-packet (blue). Panel (b) shows the difference of S(x)/2π = −Δx<sub>bn</sub>(x) profiles for a long wave-packet having Δx<sub>tr</sub> ≫ |Δx<sub>bn</sub>| and for short wave-packets having Δx<sub>tr</sub> ∼ |Δx<sub>bn</sub>|</figcaption></figure>

with the equation for the probability of trapping Π = dS<sub>small</sub>/dx = −(5 / 2) · xS<sub>0</sub>(1 − x<sup>2</sup>)<sup>1 / 4</sup>. Equation (21) provides a very convenient and useful model for investigating the properties of systems with nonlinear wave particle resonant inter- actions. This equation works as long as the area is sufficiently large that Δx<sub>tr</sub> ≫ Δx<sub>bn</sub> ∼ S<sub>0</sub> (see scheme in Fig. 15(b)), and with the decrease of (γ<sub>max</sub> − γ<sub>min</sub>) (with increase of S<sub>0</sub> this equation becomes inapplicable (see [Artemyev, Neish- tadt, Vasiliev and Mourenas, 2021] and Appendix C)).

4.2 Asymptotic solutions

Let us consider the asymptote of the solution of Eq. (19) for t → ∞ [see details in Artemyev, Neishtadt and Vasiliev, 2019]. We use the normalized variable x = 2 γ/ ( γ<sub>max</sub> − γ<sub>min</sub> ) − ( γ<sub>max</sub> + γ<sub>min</sub> ) / ( γ<sub>max</sub> − γ<sub>min</sub> ) and rewrite this equation as

$$
\frac{\partial f}{\partial t} = \begin{cases} -V(x) \frac{\partial f}{\partial x} + \frac{1}{2} \frac{\partial}{\partial x} \left( D(x) \frac{\partial f}{\partial x} \right), & x \le 0 \\ -V(x) \frac{\partial f}{\partial x} - \frac{dV(x)}{dx} (f-f') + \frac{1}{2} \frac{\partial}{\partial x} \left( D(x) \frac{\partial f}{\partial x} \right), & x \ge 0 \end{cases} \quad (22)
$$

where V(x) = V<sub>0</sub>(1 − x<sup>2</sup>)<sup>5 / 4</sup> for x ∈ [ − 1 , 1], f* = f(x*) and x* = −x. This equation, as well as Eq. (19), conserves the total number of particles:

$$
\begin{align*}
\frac{d}{dt} \int_{-1}^{1} f(x,t) dx &= \int_{-1}^{1} \frac{\partial f(x,t)}{\partial t} dx = \int_{-1}^{1} \frac{1}{2} \frac{\partial}{\partial x} \left( D \frac{\partial f}{\partial x} \right) dx + \int_{-1}^{0} \left( -\frac{\partial (Vf)}{\partial x} + f \frac{\partial V}{\partial x} \right) dx \\
&\quad + \int_{0}^{1} \left( -\frac{\partial (Vf)}{\partial x} + f^* \frac{\partial V}{\partial x} \right) dx = \int_{-1}^{1} f(x,t) \frac{\partial V(x)}{\partial x} dx + \int_{0}^{1} f^*(y,t) \frac{\partial V(y)}{\partial y} dy
\end{align*}
$$
Changing the integration variable in the second integral $y \to x$ and using $V(x) = V(y)$, $f(x, t) = f'(y, t)$ we get

$$ \frac{d}{dt} \int_{-1}^{1} f(x, t) dx = \int_{-1}^{0} f(x, t) \frac{dV(x)}{dx} dx + \int_{0}^{1} f(x, t) \frac{dV(x)}{dx} dx = 0 $$

To find an asymptotic solution to Eq. (22) for $t \to \infty$, we restrict our consideration to the case with $D = 0$. This is a reasonable approximation, because diffusion leads to much slower changes in the electron distribution than nonlinear resonant phase bunching and trapping. This approximation also requires that bunching and trapping do not compensate each other, leaving diffusion the main process (see discussion in Appendix C). Thus, we consider

$$ \frac{\partial f}{\partial t} = \begin{cases} -V(x) \frac{\partial f}{\partial x}, & x \le 0 \\ -V(x) \frac{\partial f}{\partial x} - \frac{dV(x)}{dt} (f - f^*), & x \ge 0 \end{cases} $$

We start with the construction of the general solution to this equation. The characteristic curves of this equation are

$$ \begin{aligned} \frac{dt}{1} &= \frac{dx}{V(x)} = \frac{df}{0}, & x \le 0 \\ \frac{dt}{1} &= \frac{dx}{V(x)} = \frac{df}{-\frac{dV}{dx}(f-f^*)}, & x \ge 0 \end{aligned} $$

The first equation gives

$$ t - \int_{0}^{x} \frac{dy}{V(y)} = C_{1} = \text{const}, $$

and thus for $x \le 0$ we have the general solution

$$ f(x, t) = W \left( t - \int_{0}^{x} \frac{dy}{V(y)} \right), \quad x \le 0 $$

where $W$ is an arbitrary smooth function.

The characteristics for $x \ge 0$ give

$$ t - \int_{0}^{x} \frac{dy}{V(y)} = \text{const} = C_{2}, \quad \frac{df}{dx} = - \frac{1}{V(x)} \frac{dV}{dx} (f - f^{*}), $$

and thus for $x \ge 0$ we have the general solution

$$ f(x,t) = \frac{1}{V(x)} \int_{0}^{x} \frac{dV(y)}{dy} f^{*} \left( y, t + \int_{x}^{y} \frac{dV(z)}{dz} \right) dy + \frac{1}{V(x)} \mathcal{Y} \left( t - \int_{0}^{x} \frac{dV(z)}{dz} \right), \quad x \ge 0 $$

where $\mathcal{Y}$ is an arbitrary smooth function. This solution differs from one for $x \le 0$ by the term $\sim (f - f^*)$ where $f^*$ can be considered as a source term. The requirement that the solution to be continuous at $x = 0$ gives $W = V_0 \mathcal{Y}$ where $V_0 = V(0)$.
We consider a simple symmetric $V(x)$ function, such that if $x' < 0$ is the trapping value then $x = -x^*$ is the value of release from the trapping. Thus, we can write (see the more general case with arbitrary $V(x)$ function in [Artemyev, Neishtadt and Vasiliev, 2019]):

$$f^*(x,t) = f(-x,t) = W \left( t - \int_0^{-x} \frac{dy}{V(y)} \right)$$

Using this equation and $W = V_0 \mathcal{Y}$, we obtain for the general solution at $x > 0$:

$$f(x,t) = \frac{1}{V(x)} \int_{0}^{x} \frac{dV(y)}{dy} W \left( t - \int_{0}^{x} \frac{dz}{V(z)} - \int_{y}^{-y} \frac{dz}{V(z)} \right) + \frac{V_0}{V(x)} W \left( t - \int_{0}^{x} \frac{dy}{V(y)} \right)$$

Let $f_0(x)$ be an initial distribution function, $f(x, 0) = f_0(x)$. For negative arguments of $W$ we can define $W$ through $f_0(x)$:

$$f_0(x) = W \left( - \int_0^x \frac{dy}{V(y)} \right)$$

and thus $W$ is a bounded function for negative values of its argument. We define the function $W$ for positive argument through the initial distribution as a solution of the equation

$$f_0(x) = \frac{1}{V(x)} \int_{0}^{x} \frac{dy}{V(y)} W \left( - \int_{0}^{x} \frac{dz}{V(z)} - \int_{y}^{-y} \frac{dz}{V(z)} \right) dz + \frac{V_0}{V(x)} W \left( - \int_{0}^{x} \frac{dz}{V(z)} \right),$$

One can show that this equation determines a unique bounded function $W$.

Let us define

$$\bar{W} = \lim_{x \to 1} \int_{0}^{x} \frac{dV(y)}{dy} W \left( - \int_{0}^{x} \frac{dz}{V(z)} - \int_{y}^{-y} \frac{dz}{V(z)} \right) dy$$

This is a finite limit, because $dV/dy = O(|y \mp 1|^{1/4})$ if $y \sim \pm 1$ and $dV/dy = O(y)$ if $y \sim 0$.

One can show that $W(arg) \to -\bar{W}/V_0$ for $arg \to +\infty$.

Let us consider the limit $\lim_{t \to \infty} f(x, t)$ at a fixed $x \neq \pm 1$. For $x \le 0$ we can write

$$f(x,t) = W \left( t - \int_{0}^{x} \frac{dy}{V(y)} \right) \to -\frac{\bar{W}}{V_0} = \text{const}, \quad \text{as } t \to \infty$$

Thus, for $x \ge 0$ we can write

$$f^*(x,t) \to -\frac{\bar{W}}{V_0} = \text{const}, \quad \text{as } t \to \infty$$
and

$$f(x,t) \rightarrow \frac{1}{V(x)} \left( -\frac{\bar{W}}{V_0} \int_0^x \frac{dV(y)}{dy} dy - \bar{W} \right) = -\frac{\bar{W}}{V_0}, \quad \text{as } t \rightarrow \infty$$

Therefore, $f(x,t) \rightarrow \text{const}$ as $t \rightarrow \infty$, and the asymptotic solution to Eq. (19) is a constant function $f = \text{const}$. This is an important result, because it theoretically proves that the final stage of the evolution of the electron distribution in the presence of multiple nonlinear resonant interactions is identical to the final stage of a diffusive evolution, that is, a plateau with a null gradient along the resonance curve. Figure 16 shows a numerical verification of this result. For the numerical solution, we use $S(x)$ given by Eq. (21), but we do not make any assumption of smallness $\gamma_{\max} - \gamma_{\min}$, i.e., the function $S(x)$ does not necessarily corresponds to the short wave-packet approximation. Numerical results show that the initially localized maximum of $f(x)|_{t=0}$ quickly (over a time scale $\sim (B_w/B_0)^{-1/2}$) evolves toward $f = \text{const}$. Therefore, there are two main differences between nonlinear and quasi-linear evolution: (1) the formation of new phase space gradients (like beam structures) due to nonlinear resonant interactions in the transient initial phase, (2) the much shorter time-scale of evolution to the final stage in the case of nonlinear interactions. The first difference is important mainly when we compare very quick phenomena associated with resonant wave-particle interactions which include only a few such interactions; a good example is the phenomenon of microburst precipitation [see details in the Section 8 and in Chen et al., 2020, 2022; Kang et al., 2022; Shumko et al., 2021]. The second difference is more important for long-term radiation belt dynamics [see details in the Section 7 and in Artemyev, Mourenas, Zhang and Vainchtein, 2022].

## 4.3 Effect of a non-constant bounce period

Equation (19) can be generalized for the case when the period between resonant interactions depends on energy. For electron resonance with a monochromatic whistler-mode wave, this period $\tau$ is equal to the half of the bounce period along a magnetic field line (in case of waves propagating away from the equatorial source region). Thus, in dipole magnetic field we have

$$\tau = \frac{\tau_b}{2} = 2 \int_0^{s_{\max}} \frac{m_e \gamma ds}{p_{\parallel}} = \frac{2 m_e}{c \sqrt{1 - \gamma^{-2}}} \int_0^{s_{\max}} \left( 1 - \sin^2 \alpha_{eq} \frac{\Omega_0(s)}{\Omega_0(0)} \right)^{-1/2} ds$$

where $s_{\max}$ is determined by $\sin^2 \alpha_{eq} \Omega_0(s_{\max}) = \Omega(0)$. Thus, for $h = \text{const}$ from Eq. (4) we have $\tau = \tau(\gamma; h)$:

$$\tau = \frac{2m_e}{c\sqrt{1-\gamma^{-2}}} \int_{0}^{s_{\max}} \left( 1 - \frac{2(\gamma - h/m_e c^2) \Omega_0(s)}{\gamma^2 - 1} \frac{\Omega_0(s)}{\omega} \right)^{-1/2} ds$$
The generalized form of Eq. (19) can be written as [Artemyev, Neishtadt and Vasiliev, 2021; Artemyev et al., 2017a]:

$$ \begin{aligned} \frac{\partial f}{\partial t} &= \frac{1}{2} \frac{\partial}{\partial \gamma} \left( D_{\gamma\gamma} \frac{\partial f}{\partial \gamma} \right) - \frac{\partial}{\partial \gamma} \left( \frac{\langle \Delta\gamma \rangle}{\tau} f \right) + \frac{f}{\tau} \frac{d \langle \Delta\gamma \rangle}{d\gamma} (1 - \Theta) + \frac{d \langle \Delta\gamma \rangle}{d\gamma} \frac{f^*}{\tau^*} \Theta \\ &= \frac{1}{2} \frac{\partial}{\partial \gamma} \left( D_{\gamma\gamma} \frac{\partial f}{\partial \gamma} \right) - V \frac{\partial f}{\partial \gamma} - \frac{dV}{d\gamma} \left( f - f^* \frac{\tau}{\tau^*} \right) \Theta + V \frac{d \ln \tau}{d\gamma} \left( \frac{\tau}{\tau^*} f^* \Theta - f (1 - \Theta) \right) \end{aligned} \quad (23) $$

where $V = \langle \Delta\gamma \rangle/\tau$, $f^* = f(\gamma^*)$ and $\tau^* = \tau(\gamma^*)$ with $S(\gamma) = S(\gamma^*)$, and $\Theta$ is defined by Eq. (18). Note $\langle \Delta\gamma \rangle = \omega\langle \Delta I \rangle_\xi/m_e c^2$ and $\langle \Delta I \rangle_\xi$ is given by Eq. (12). We may introduce a new variable $J$ as

$$ \frac{d\gamma}{dJ} = \frac{2\pi}{\tau(\gamma)} $$

Using this new variable, we introduce $\tilde{f}(J, t)$, $D_{JJ}$, $V_J$:

$$ f = \frac{\tilde{f}\tau}{2\pi}, \quad V = \frac{2\pi V_J}{\tau}, \quad D_{\gamma\gamma} = \frac{4\pi^2 D_{JJ}}{\tau^2} $$

Then, we can rewrite Eq. (23) as

$$ \frac{\partial \tilde{f}}{\partial t} = -V_J \frac{\partial \tilde{f}}{\partial J} + \frac{1}{2} \frac{\partial}{\partial J} \left( D_{JJ} \frac{\partial \tilde{f}}{\partial J} \right) + \frac{\partial V_J}{\partial J} (f^* - f) \Theta \quad (24) $$

This equation coincides with Eq. (19), but instead of $\gamma$ (or $I = m_e c^2 \gamma / \omega - \text{const}$) we should use the variable $J$. The asymptotic solutions of Eqs. (19) and (24) are the same. Therefore, the system equations are not changed in the more general case with $\tau_b = \tau_b(\gamma)$, but the energy space is modified: instead of having $\gamma$ linearly proportional to $I_x$ through Eq. (4), we now have $J$ given by $dJ/d\gamma = \tau_b(\gamma)/2\pi$.

## 4.4 Simulation Techniques

In this section, we briefly review several possible schemes for the numerical simulation of electron dynamics in a system with multiple nonlinear resonances. We focus on the mapping technique, which is first introduced in the case of a single monochromatic wave, and then generalized for a wave ensemble. But we also compare this technique with the well-developed and quite powerful Green function approach [Furuya et al., 2008; Omura et al., 2015] and with an analytically derived version of the generalized Fokker-Planck equation, which relies on a Probabilistic approach. All these techniques are based on the same equations and physical concepts, and the main (if not only) difference concerns the numerical implementation of these techniques. Note that, although we do not provide results of simulations of realistic (observed by spacecraft) events in this section, Appendix E includes a detailed analysis of two observational events modeled with the mapping technique.
## 4.5 Mapping for a single wave

Kinetic equations (19) and (23) describe the dynamics of the electron distribution function in a system with multiple nonlinear resonances. Since the evolution of the distribution consists of information about multiple electron trajectories, instead of solving the kinetic equation we may solve a large set of equations describing each individual electron trajectory. The most detailed Hamiltonian equations (2) trace all electron coordinates, fast and slow, but kinetic equations (19,23) describe only the electron energy evolution (or I) for h = const. Therefore, the corresponding equation for electron trajectories should also involve only equations for energy, and should not describe electron motion between resonant interactions. The closest analog of such equations describing electron energy change at resonances is the Chirikov map [Chirikov, 1987] that should give $\gamma_{n+1} = \gamma_n + \Delta\gamma(\gamma_n)$, where n is the number of resonant interactions [see such maps for electron scattering by whistler-mode waves in Benkadda et al., 1996; Khazanov et al., 2014]. For kinetic equations (19,23) this map can be written as [Artemyev, Neishtadt and Vasiliev, 2020a,b]:

$$ \gamma_{n+1} = \gamma_n + \begin{cases} \Delta\gamma_{trap}(\gamma_n), & \xi_n \in [0, \Pi(\gamma_n)) \\ \Delta\gamma_{bunch}(\gamma_n), & \xi_n \in (\Pi(\gamma_n), 1] \end{cases} \qquad (25) \\ \Pi(\gamma_n) = -d\Delta\gamma_{bunch}(\gamma_n)/d\gamma $$

where $\xi_n$ is a random variable with a uniform distribution within [0, 1]. The value of $\xi_n$ is indeed determined by the phase gain between two resonant interactions, and this gain can be considered as a random variable, which is a nontrivial result. A rigorous proof of $\xi_n$ properties can be found in [Gao et al., 2023]. A simplified version of this proof is provided in Appendix D.

For the simplified model (21) this map can be rewritten as

$$ x_{n+1} = x_n + \begin{cases} -2x_n, & \xi_n \in [0, \Pi(x_n)) \\ -\Delta x_0 \cdot (1 - x_n^2)^{5/4}, & \xi_n \in (\Pi(x_n), 1] \end{cases} \qquad (26) \\ \Pi(x_n) = -\Delta x_0 \frac{5}{2} x_n (1 - x_n^2)^{1/4} $$

with $\Delta\gamma_{bunch} = -S/2\pi$ and $\Delta x_0 = (S_0/2\pi) \propto \sqrt{B_w/B_0}$. Figure 16(a) shows several examples of $x_n$ trajectories obtained with the map (26). The dynamics of $x$ consists of rare and large jumps (due to phase trapping) with $x$ increase and regular $x$ drift to smaller values (due to phase bunching). This dynamics consists of basic elements, phase trapping and phase bunching, resembling well electron energy changes due to a single resonant interaction, see Fig. 3. The fine balance between trapping and bunching (the probability of trapping $\Pi = (2\pi)^{-1}dS/dx$, and magnitude of bunching $\Delta x = -S/2\pi$) results in a quasi-periodical $x$ motion between small ($x \sim -1$) and large ($x \sim 1$) values.

For realistic systems, $S(\gamma)$ should be derived based on actual wave field and background magnetic field latitudinal profiles. Using such realistic $S(\gamma)$, we plot $\gamma_n$ trajectories in Fig. 16(b) and compare them with trajectories obtained by direct numerical integration of Hamiltonian equations (2), which are plotted in Fig. 16(c). This comparison demonstrates that the map (25) describes well the main constitutive elements of the dynamics of $\gamma$. Note that the map (25) includes a significant randomization factor (the random $\xi$) and, therefore, we cannot expect a one-to-one similarity (at each time) between the $\gamma_n$ profiles obtained from the map
<figure><img src="image_16.png" /><figcaption>**Fig. 16** Panel (a) shows three trajectories $x_n$ for map (26) with $\Delta x_0 = 10^{-3/2}$. Panel (b) shows three trajectories $m_e c^2 (\gamma_n - 1)$ for map (25). Panel (c) shows three trajectories obtained by numerical integration of Hamiltonian equations (2). Results in panels (b,c) are obtained for initial energy 100 keV, initial pitch-angle $\alpha_{eq} = 60^\circ$, wave amplitude 300 pT, and L-shell= 6. Other systems parameters are the same as in Fig. 2.</figcaption></figure>

and from numerical integration of Hamiltonian equations, even if both trajectories start with the same initial conditions. However, the most important point is that the statistical properties of the $\gamma_n$ dynamics are the same for trajectories obtained by the mapping technique and by direct numerical integration of Hamiltonian equations.

## 4.6 Nonlinear resonances with multiple waves

Kinetic equations (19) and (23) have been derived for a 1D system with $h = \text{const}$, i.e., for a system including only one monochromatic wave. If we consider a more realistic situation where electrons interact resonantly with various waves
<figure><img src="image_17.png" /><figcaption>**Fig. 17** Top panels show resonance curves (black) and electron trajectories in the energy/pitch-angle space for Hamiltonian (2), when the electron interacts only with the first whistler-mode wave (a), only with the second whistler-mode wave (b), or with both whistler-mode waves (c). Wave frequencies are $\omega_1 = 0.4\Omega_{eq}$, $\omega_2 = 0.2\Omega_{eq}$, both wave amplitudes are 300 pT, $L = 6$, and other system parameters are from Fig. 2.</figcaption></figure>

having different frequencies, we cannot use this 1D approximation, because for each wave frequency we will correspond to a different $h$ given by Eq. (4). This means that the frequency value determines the shape of resonance curves in the velocity, energy space, and for different frequencies the wave-particle resonance move electrons along different curves. Moreover, $I$ also depends on wave frequency (or wave number $k$), and for each wave frequency we have a corresponding $I$, i.e., in the system with two wave frequencies we have $I_1$ and $I_2$. During the resonant interaction with the first wave, $I_1$ changes but $I_2$ is conserved, and vice versa. Therefore, for each resonance electrons move along the corresponding resonance curve. The conservation of $h$ and of one of the momenta ($I_1$ or $I_2$) leads to a one-dimensional electron dynamics in the energy/pitch-angle space, and if there is only one wave in the system electrons never leave the corresponding single resonance curve. However, electron dynamics becomes 2D in the presence of two waves, when both $I_1$ and $I_2$ change. Figure 17 illustrates this effect by showing electron resonant interactions with a single wave and with two waves. The electron moves along resonance curves and jumps between these curves due to $I_{1,2}$ jumps. Accumulation of such jumps between resonance curves ultimately leads a single electron trajectory to cover the entire energy/pitch-angle space [see more details in Artemyev, Neishtadt, Vasiliev, Zhang, Mourenas and Vainchtein, 2021].

As emphasized earlier, to generalize kinetic equations (19) and (23) for systems with multiple waves, we need to use a probabilistic approach. We introduce a probability distribution function $P(B_w, \omega)$ that determines the wave characteristics $(\omega, B_w)$ during the next resonant interaction. And then we average all operators in these equations over $P$. Since such averaging significantly complicates the kinetic equations, we need an alternative approach. Let us discuss three possible methods for modeling the evolution of the electron distribution due to multiple resonances with a wave ensemble. Although we are speaking about a wave ensemble,
ble, it is assumed that electrons interact resonantly with only one monochromatic wave at a time, without the wave resonance overlap [see discussion in An, Wu and Tao, 2022; Gan et al., 2022; Tao et al., 2013], while the resonant waves may have different properties during different bounce period.

**Green function approach** The Green function approach [Furuya et al., 2008; Omura et al., 2015] assumes that systems with multiple different waves (or with a single wave with evolving characteristics) can be described using Eq. (1) with a kernel derived from test-particle simulations. Figure 18(a) shows the basic scheme of this approach: test particle simulations provide the probability distribution function of energy and pitch-angle changes, and this function is used to construct the kernel and rewrite the Smoluchowski equation as

$$f_{n+1}(E, \alpha_{eq}) = \int_{0}^{\infty} dE' \int_{0}^{\pi/2} f_n(E', \alpha'_{eq}) G(E, \alpha_{eq} | E', \alpha'_{eq}) d\alpha'_{eq} \\ G(E, \alpha_{eq} | E', \alpha'_{eq}) = C[\delta(E - E', \alpha_{eq} - \alpha'_{eq})]$$

where $\partial f / \partial t$ is replaced by a discrete difference $f_{n+1} - f_n$ during each time step of one bounce period, and $f_n$ is included into the Green function

$$\int_{0}^{\infty} dE' \int_{0}^{\pi/2} G(E, \alpha_{eq} | E', \alpha'_{eq}) d\alpha'_{eq} = 1$$

while $C$ is the particle scattering operator describing energy and pitch-angle change during a single bounce period. This operator can be obtained from test particle simulations: a large particle ensemble can be traced across the resonance and the variations of their energy and pitch-angles can be combined to determine the probabilities of all possible transports $(E', \alpha'_{eq}) \rightarrow (E, \alpha_{eq})$ [Furuya et al., 2008]. This is quite a powerful approach for quantitatively describing multiple nonlinear resonant interactions affecting the dynamics of energetic electron fluxes. Several important results have been obtained this way, like that description of the formation of relativistic/ultra-relativistic electron populations due to turning acceleration by field-aligned waves [Omura et al., 2015] or due to the Landau and high-order resonances with oblique waves [Hsieh et al., 2020; Hsieh and Omura, 2017a], and investigations of rapid scattering and losses of energetic electrons due to a combination of the Landau and cyclotron nonlinear resonances [Hsieh and Omura, 2023; Hsieh et al., 2022]. These simulations demonstrated that the Green function approach is really promising and, combined with the observed wave distributions, it can potentially replace and supersede standard simulations of radiation belt dynamics based on the Fokker-Planck diffusion equation. The main technical difficulty of this approach is the need to predefine the resolution in energy and pitch-angle for the $C$ operator that will be determined using test particle simulations. For instance, typical nonlinear wave-particle interactions with intense waves often cover at least three order of magnitude of energies, [1, 1000] keV, whereas a simultaneous inclusion of electron scattering by weak waves requires a minimum energy bin size about 10–100eV [see typical diffusion rate magnitudes in Glauert and Horne, 2005; Horne, Kersten, Glauert, Meredith, Boscher, Sicard-Piet, Thorne and Li, 2013]. Therefore, to accurately incorporate the effects from both intense and weak waves, one would need $10^4–10^5$ energy bins and about the same number
of pitch-angle bins. Thus filling the corresponding 2D matrix for $\mathcal{C}$ is computationally very expensive. Consequently, the Green function approach is useful mostly for describing electron flux dynamics in a system with intense waves (during active geomagnetic conditions), when the energy and pitch-angle bin sizes can remain sufficiently large. The two alternative methods described below, the Probabilistic approach and the Mapping technique, would require similarly huge numbers of small energy and pitch-angle bins for an accurate description of electron dynamics in the presence of both intense and weak waves, but as we will see, with potentially different intrinsic accuracy and total CPU time. The main differences of the (Probabilistic approach and Mapping technique) from the Green function approach consists in the analytical evaluation of the basic properties (characteristics) of wave-particle resonant interactions. This improves the accuracy of the evaluation of such characteristics, but reduces the flexibility for including comprehensive details of wave-particle interactions (like wave frequency drift and wave-packet structure). Roughly speaking, the Green function approach is ideal for a detailed modeling of short-term dynamics of electron fluxes, when peculiarities of wave-field can play the most important role, the Probabilistic approach is optimal for the inclusion of nonlinear resonant effects into existing numerical schemes of radiation belt dynamics (an even simpler technique for such an inclusion is discussed in Section 7), and the Mapping technique is the most suitable for observation-based modeling of meso-scale events, when wave characteristics are derived from spacecraft observations, and for incorporation of wave-particle resonant interactions into global test-particle simulations (see discussion in [Artemyev, Neishtadt and Angelopoulos, 2022] and in Section 6.3).

**Probabilistic approach** An alternative to the Green function approach and a different way of rewriting the Smoluchowski equation was proposed in Vainchtein et al. [2018]. This approach is based on the idea of separating phase trapping and bunching processes and their analytical evaluations. In this case, the discretization of the electron energy, pitch-angle space allows constructing a matrix $\hat{R}$ of energy, pitch-angle probability jumps. Within this approach, Eq. (1) takes the form

$$ \frac{\partial f_{ij}}{\partial t} = -\frac{N_{ij}}{\tau_{ij}} \int_{W} \mathcal{P}(W) dW + \left( \int_{W} \sum_{k,l} \hat{M}_{ij}^{kl}(W) \mathcal{P}(W) dW \right) f_{ij} \quad (27) $$

where $f_{i,j} = f(E_i, \alpha_{eq,j})$, $N_{ij} = N(E_i, \alpha_{eq,j})$ is the number of resonant interactions that particles undergo during a single bounce period $\tau_{ij} = \tau(E_i, \alpha_{eq,j})$, $W = (B_w, \omega)$ are wave characteristics and $\mathcal{P}(W)$ is the probability distribution function of wave characteristics normalized in such a way that $\int_W \mathcal{P}(W)dW$ is the ratio of the total time interval of spacecraft wave measurements to the cumulative time interval of observations of intense waves resonating with electrons nonlinearly. The operator $\hat{M}_{ij}^{kl}(W)$ is a 4D matrix describing the transformation of the 2D matrix $f_{ij}$ of initial energy, pitch-angle to a $f_{kl}$ matrix of energy, pitch-angle after one bounce period.

A schematic view of this probabilistic approach is provided in Fig. 18(b): for each energy, pitch-angle bin and each pair of wave characteristics $(B_w, \omega)$, we determine the probabilities of particle coming and leaving due to the phase trapping and phase bunching. These probabilities are determined using the analytical model
<figure><img src="image_18.png" /><figcaption>**Fig. 18** (a) Schematic view of the Green function approach. **Numerical simulations** are performed to obtain energy, pitch-angle distribution of electrons after resonant interactions (bottom panel) for electrons having the same initial energy, pitch-angle (top panel). Such distributions, obtained for different wave characteristics and for all energy, pitch-angle bins are combined into the Green function $G(E, \alpha_{eq} | E', \alpha'_{eq})$ that describes the probability of $(E', \alpha'_{eq}) \rightarrow (E, \alpha_{eq})$ transition [see details in Omura et al., 2015]. (b) Schematic view of the probabilistic approach, where for each energy, pitch-angle bin the analytical model (for a given distribution of wave characteristics) provides the probability $\Pi$ and energy, pitch-angle changes due to phase trapping (green) and due to phase bunching (red) [see details in Vainchtein et al., 2018].</figcaption></figure>

of wave-particle nonlinear interaction and, therefore, do not require large numerical test-particle simulations. Next, the matrices $\hat{R}(W)$ are averaged over the probability distributions of wave characteristics. The averaged matrix is used in the equation for the electron distribution function. The main advantage of this approach is that it is based on an analytical model of energy changes and trapping probabilities that can be rapidly evaluated over a wide parametric domain with arbitrarily high accuracy. The main disadvantage is that such analytical equations do not allow a simple generalization for more sophisticated and realistic waveforms (see discussion in Section 5). However, even such an idealized approach can still provide important information about systems with nonlinear wave-particle interactions. Specifically, numerical simulations of electron distribution dynamics with realistic and wave distributions, $P(B_w, \omega)$, have shown the importance of appropriately choosing the initial electron distributions [Vainchtein et al., 2018]. The evolution of the electron distribution function due to nonlinear wave-particle interactions leads to a rapid change of the total electron energy, but the energy of a whistler-mode wave is not sufficient in itself to provide such an deceleration/acceleration [see discussion and estimates in Shklyar, 2011]. Therefore, the large population of scattered (decelerated) electrons provides the only available energy source for trapped electron acceleration. The corresponding energy balance imposes several constraints on the initial electron distribution function for models describing nonlinear interactions. This initial distribution should be consistent with the conservation of the total energy for solutions of Eq. (27). Vainchtein et al. [2018] checked the total particle energy variation for different initial pitch-angle distributions: an initial distribu-
tion with transverse anisotropy ($F \propto \sin \alpha_{eq}$) immediately loses energy, whereas a field-aligned anisotropic distribution ($f \propto \cos \alpha_{eq}$) first gains energy before losing it. More specific distributions (butterfly $f \propto \sin 2\alpha_{eq}$ sometimes observed in the radiation belts [see statistics in Åsnes et al., 2005] and a distribution with field-aligned beams $f \propto \sin \alpha_{eq} \cos^4 \alpha_{eq}$ observed during particle injections [e.g., Mozer et al., 2016]) show energy increase during the first stage of the evolution. This suggests that for a $f$ to result in total energy conservation it should be more isotropic (i.e., an intermediate state between transverse anisotropy $\sim \sin \alpha_{eq}$ and field-aligned anisotropy $\sim \cos \alpha_{eq}$). Moreover, under realistic conditions, particle anisotropy is energy dependent, and can be different for $< 100$ keV particles and $\sim 1$ MeV particles. Therefore, further investigations of particle distributions consistent with the assumption of total energy conservation are needed to accurately model nonlinear wave-particle interactions [see, also, Shklyar and Luzhkovskiy, 2023; Shklyar, 2017, for details on the relation between trapped and scattered electron populations in a system with small variations of wave intensity].

**Mapping technique** The mapping technique can easily be generalized to multi-wave systems: instead of Eq. (25) for a single wave, we should use

$$ \gamma_{n+1} = \gamma_n + \begin{cases} \Delta\gamma_{trap}(\eta_n, \gamma_n), & \xi_n \in [0, \Pi(\eta_n, \gamma_n)) \\ \Delta\gamma_{bunch}(\eta_n, \gamma_n), & \xi_n \in (\Pi(\eta_n, \gamma_n), 1] \end{cases} \quad (28) $$

where the index $\eta_n$ is a random number with a probability distribution such that $\eta$ determines the probability of resonant interaction with different waves, which characteristics are given by the distribution $P(B_w, \omega)$. Thus, before each resonant interaction (each $\gamma$ change) with a particular particle, we determine from $P(B_w, \omega)$ the corresponding wave, and use the corresponding $S$ to determine energy, pitch-angle change for this particle. Particle pitch-angles should be recalculated with new $\gamma_{n+1}$ and corresponding wave frequency $\omega$ accordingly to Eq. (4). Figure 19(a) shows several examples of electron trajectories evaluated with the mapping (28) for the wave probability distribution $P(B_w, \omega)$ from panel (c). This distribution contains waves of different frequencies, and for each frequency the resonant interaction will move electrons along the $h(\omega) = \text{const}$ curve (see Eq. (4)). Phase bunching moves electrons along $h = \text{const}$ curves with energy and pitch-angle decrease (see the prolonged intervals with energy decrease in panel (a)), whereas phase trapping moves electrons along the same curves, but with an energy increase (see the rare positive energy jumps in panel (b)). For a single frequency system, there would be only one curve $h = \text{const}$ for each initial electron energy, pitch-angle, and thus electrons would move along this curve in the energy, pitch-angle space. But the system in Fig. 19 contains multiple frequencies and, thus, resonant interactions will provide complex electron dynamics in the 2D energy, pitch-angle space. Figure 19(b) shows two trajectories from panel (a) in the energy, pitch-angle space: drifting from one resonant curve to another one, electrons may increase their energy with a net pitch-angle decrease (see also Fig. 17(c)), which is purely due to resonant interactions with different waves at multiple frequencies. Note that the distribution shown in Fig. 19(c) contains many waves with an amplitude insufficiently large to resonate with a particular electron nonlinearly. Such waves should provide electron diffusive scattering, which is not included into the mapping given by Eq. (28). Such diffusive scattering can be incorporated into the mapping scheme
<figure><img src="image_19.png" /><figcaption>**Fig. 19** Panel (a) shows five examples of electron trajectories (energy vs. time) for cyclotron resonance with waves having the distribution shown in panel (c). All electrons have the same initial energy and pitch-angle. Panel (b) shows two trajectories from panel (a) in the (energy, pitch-angle) space. Background magnetic field and plasma density characteristics are the same as in Fig. 2.</figcaption></figure>

using standard (Chirikov-type, see [Chirikov, 1987]) maps. Examples of such maps can be found in [Khazanov et al., 2013, 2014; Vasilev et al., 1988; Zaslavskii et al., 1989].

The map (28) can be used with a 3D wave characteristic distribution $\mathcal{P} = \mathcal{P}(B_w, \omega, \theta)$, but then, we need to use a different Eq. (4) for a different resonance number. Generally two main resonances, Landau resonance with very oblique waves ($n_r = 0$) and the first cyclotron resonance with quasi-parallel waves ($n_r = 1$), are included [Artemyev, Neishtadt, Vasiliev, Zhang, Mourenas and Vainchtein, 2021]. The difference is in relation of $I$ from Eq. (4) and the electron magnetic moment $I_x$: for $n_r = 0$ Eq. (4) becomes $I_x = \text{const}$, whereas for $n_r = 1$ Eq. (4) becomes $m_e c^2 \gamma - \omega I_x = \text{const}$.

Using a similar approach with a probability distribution for the waves, we may generalize the map (28) for a system with both types of whistler-mode waves: quasi-parallel and very oblique. In this case, before each resonant interaction, we determine the wave type, which in turn determines the type of resonance – Landau or cyclotron. Figure 20(a) shows examples of electron trajectories in the case
<figure><img src="image_20.png" /><figcaption>**Fig. 20** Panel (a) shows five examples of electron trajectories (energy vs. time) for Landau and cyclotron resonances with waves having the distribution shown in panel (c). All electrons have the same initial energy and pitch-angle. Panel (b) shows two trajectories from panel (a) in the (energy, pitch-angle) space. Background magnetic field and plasma density characteristics are the same as in Fig. 2.</figcaption></figure>

where the probability for electrons to meet very oblique and quasi-parallel waves is the same. Electron energy dynamics resembles results obtained for a system including only cyclotron resonance (see Fig.19(a)): prolonged periods of energy decrease due to the phase bunching (due to both Landau and cyclotron resonances) are intermittently interrupted by rare positive jumps of energy due to phase trapping (related to both Landau and cyclotron resonances). The main difference with the cyclotron-resonance-only system (including only parallel waves) from Fig. 19 is that energy increases/decreases in Landau resonance correspond to pitch-angle decreases/increases, such that Landau phase trapping, in particular, results in simultaneous acceleration and pitch-angle decrease. Together, nonlinear wave-particle interactions through cyclotron and Landau resonances with parallel and oblique waves provide both significant electron acceleration and efficient transport toward the loss-cone (see Fig.20(b) and compare with Fig.19(b)) and, therefore, lead to an enhanced precipitation of energetic (accelerated) electrons compared with a system including only cyclotron resonance with parallel waves [see more details of this effect in Hsieh et al., 2022; Mourenas et al., 2016].
Such mapping simulations of individual electron trajectories (energy as a function of the resonant interaction number or time) can be combined to comprehensively describe the full dynamics of the electron distribution function [Artemyev, Neishtadt, Vasiliev, Zhang, Mourenas and Vainchtein, 2021; Artemyev, Zhang, Zou, Mourenas, Angelopoulos, Vainchtein, Tsai and Wilkins, 2022; Zhang, Artemyev, Angelopoulos, Tsai, Wilkins, Kasahara, Mourenas, Yokota, Keika, Hori, Miyoshi, Shinohara and Matsuoka, 2022]. Therefore, the mapping technique may provide the same level of system description as the probabilistic approach and Green function approach.

## 5 Effects of short wave-packets

The general theoretical results described by Eqs. (19) and (28) were obtained in the simplified case of a plane wave (such that the term ~ $B_w$ contains only $\sin \zeta$ in Hamiltonian (2)). However, in realistic plasma systems, such very long wave-packets (such that the wave packet envelope can be neglected) are rarely observed, and most of the whistler-mode waves present in the outer radiation belt are propagating in the form of short wave-packets (see examples in Fig. 21(a,b) and [Zhang, Mourenas, Artemyev, Angelopoulos, Bortnik, Thorne, Kurth, Kletzing and Hospodarsky, 2019; Zhang, Thorne, Artemyev, Mourenas, Angelopoulos, Bortnik, Kletzing, Kurth and Hospodarsky, 2018]). Typical wave-packets (or subpackets) comprise several (up to ten) wave periods and reach moderate peak amplitudes ($B_w \approx 100-300$ pT), but the most intense ($B_w > 500$ T) wave-packets can be quite long and may include more than hundred wave periods. The formation of such long wave-packets is associated with the classical nonlinear mechanism of whistler-mode chorus wave generation [see reviews in Demekhov, Taubenschuss and Santolík, 2017; Omura, 2021; Omura et al., 2013; Tao et al., 2021, 2020], and it is also the case for moderately long but intense subpackets [Chen et al., 2024]. On the other hand, most of the moderate intensity short wave-packets are likely generated due to wave superposition, also called wave beating [see details in An, Wu and Tao, 2022; Mourenas et al., 2022; Nunn et al., 2021; Zhang, Mourenas, Artemyev, Angelopoulos, Kurth, Kletzing and Hospodarsky, 2020]. Detailed comparison of wave generation models and spacecraft statistics of whistler-mode wave-packets confirms this scenario of separation of short and long wave-packets [Zhang et al., 2021]. Note that basic empirical models of whistler-mode waves do not include information about wave-packets [Agapitov et al., 2013; Meredith et al., 2012], and this significantly complicates the investigation of nonlinear resonant effects [see discussion in Allanson et al., 2024]. In this review, we utilize wave-packet statistics published in [Zhang, Mourenas, Artemyev, Angelopoulos, Bortnik, Thorne, Kurth, Kletzing and Hospodarsky, 2019; Zhang, Mourenas, Artemyev, Angelopoulos, Kurth, Kletzing and Hospodarsky, 2020].

To model this effect for a monochromatic wave, we may adopt the scheme proposed in [Matsoukis et al., 2000; Wykes et al., 2001] and rewrite the field term in Eq. (2) as

$$B_w(\lambda) \sin \zeta \rightarrow B_w(\lambda) \exp(-A \sin^2(\phi/2\ell)) \sin \zeta$$
where parameter $\mathcal{A}$ determines the depth of wave field modulation, $\ell$ is the wave-packet size in wave length (wave periods) units. Figure 21(c) shows examples of this wave-packet model.

An alternative and more sophisticated wave-packet model for whistler-mode chorus waves (with frequency drift) has been proposed in Furuya et al. [2008] and actively used for investigation of the effect of wave field modulation on efficiency of electron resonant dynamics [Gan et al., 2020; Tao et al., 2013; Tao, Bortnik, Thorne, Albert and Li, 2012]. This model includes the difference between group and phase wave velocities, and thus describes wave-packet modification during the propagation, because different waves contributing to each wave-packet propagate with different velocities. Hiraga and Omura [2020]; Hsieh et al. [2020]; Kubota and Omura [2018] have investigated in detail this chorus wave-packet model and incorporated it into a general scheme of electron flux evaluation in systems with multiple nonlinear resonances based on the Green function approach [Omura et al., 2015].

Resonant interactions with such short wave-packets differ from the simplified picture of interaction with infinitely long packets in four main aspects:

*Nonlinear phase bunching* is affected mostly by the wave-packet propagation in latitude. For an infinite wave-packet, the wave amplitude at the latitude of resonance (for fixed energy and pitch-angle) is determined by a fixed $B_w(\lambda)$ profile that does not depend on time. But for short wave-packets, the wave amplitude depends on the time delay between the particle and the wave-packet arriving at the latitude of resonance with that particle (see scheme in Fig. 22(a)). This significantly randomizes the effect of phase bunching and gives a broader energy distribution in comparison with results for the infinite plane wave. Figure 22(b) shows the probability distribution of energy changes due to a single resonant interaction with a very long wave-packet ($\ell = 300, \ell = 500$) and short wave-packets ($\ell = 10$). In these simulations, the upper limit of the negative energy jump (the phase bunching effect) is the same, because it is determined by the peak wave amplitude, which is the same in all cases. However, for short wave-packets there is a broader distribution of negative energy changes, because particles can encounter various (randomized) instantaneous wave amplitudes at the resonance. These simulations have been set up in such a way that initial electron coordinates are uniformly distributed within the range corresponding to electrons arriving to the resonance when the wave-packet amplitude there exceeds 1/100 of the peak value, i.e., we excluded from the simulations the particles arriving into the resonance between two successive wave-packets and which would not significantly interact with the waves.

*The probability of nonlinear phase trapping* is also mostly affected by the strong wave field gradient at the wave-packet edges. For an infinite wave-packet, the trapping probability is determined by the latitudinal gradients of the wave field and background magnetic field, but for short wave-packets this probability is instead determined by the gradient of the wave-packet envelope [Bortnik et al., 2008]. In contrast to background latitudinal gradients, with typical spatial scales of $\sim LR_E \sim 10^4$km, the typical scale of the wave envelope gradient can shrink to only a few wavelengths, $\sim 10-100$km, which represents a huge increase of
<figure><img src="image_21.png" /><figcaption>**Fig. 21** Panel (a,b) shows four examples of whistler-mode wave-packets propagating with a small wave normal angle: measurements from THEMIS [Le Contel et al., 2008] and Van Allen Probes [Kletzing et al., 2013] search-coil magnetometers (L-shell and MLT are indicated within panels). Panel (c) shows two examples of model wave-packets.</figcaption></figure>

$dS/d\gamma \propto dB_w/d\lambda$ along the resonant trajectory. Figure 22(c) shows how the trapping probability depends on wave-packet size for a given wave-packet model, and increases for shorter packets. In particular, for a short wave-packet, electrons interacting resonantly with the packet leading edge have a larger probability of phase trapping.

*Energy changes due to nonlinear phase trapping* are also strongly affected by the wave-packet size. For such short wave-packets, in principle phase trapping becomes possible at any latitude of resonance, without a strong influence of the background profiles of wave and magnetic fields. However, the resulting strong increase of the probability of trapping is compensated by the brevity of this trapping, because particles rapidly escape from the resonance (and from the trapping) at the edges of wave-packets. For the first cyclotron resonance, the resonant electrons
<figure><img src="image_22.png" /><figcaption>**Fig. 22** Panel (a) shows a schematic view of electron resonant interactions with short wave-packets. The three main parameters are the dependencies of (1) the instantaneous wave amplitude (determining the efficiency of phase bunching) at the resonance and (2) the wave amplitude gradient (determining the phase trapping probability) on the initial electron location relative to the wave-packet generated at the equator, and (3) the wave-packet size (which determines the energy gain due to phase trapping, instead of the background magnetic field gradient in the case of a very long packet). Panel (b) shows the distribution of energy changes due to a single resonant interaction with long ($\ell$ = 300 and $\ell$ = 500) and short ($\ell$ = 10) wave-packets. Panel (c) shows the distribution of the trapping probability, $\Pi$, for a single resonant interaction, as a function of wave-packet length, $\ell$. For numerical simulations shown in panels (b,c) we use field-aligned whistler-mode waves with a peak amplitude of 500 pT and electrons with an initial energy of 100 keV and an equatorial pitch-angle of 30°; other system characteristics are the same as in Fig. 2.</figcaption></figure>
and waves are propagating in opposite directions, corresponding to a maximum
possible trapping time $\Delta t_{trap} \approx 2\pi\ell/k|v_g + v_R|$, where $v_g = \partial\omega/\partial k_\parallel$ denotes the
wave group velocity and $v_R = (\omega - \Omega_0/\gamma)/k_\parallel < 0$ the velocity of the particle in
cyclotron resonance with the wave. Using a simple approximation $v_g \approx \omega/k_\parallel$ [for
parallel propagating waves, the typical ratio of $k_\parallel v_g/\omega$ is 1 - 2, see Stix, 1962],
we obtain $\Delta t_{trap} \approx 2\pi\ell/\Omega_{ce}$. Thus, a finite wave-packet size directly sets an up-
per limit on the electron acceleration. A possible exception has been considered
in [Hiraga and Omura, 2020], where the idea of *successive trapping* into series of
wave-packets (wave-packet train) has been proposed. Such a mechanism allows
to neglect the effect of particle detrapping at the edges of wave-packets, because
such detrapped particles will soon be trapped again by the next packet. However,
*successive trapping* requires that the wave phase remain coherent across the en-
tire wave-packet train, without strong phase stochastization at the wave-packet
edges, which limits its applicability [see examples of observations of both coherent
and non-coherent wave packet trains in Zhang, Agapitov, Artemyev, Mourenas,
Angelopoulos, Kurth, Bonnell and Hospodarsky, 2020].

In Figure 22(b), wave-packets are assumed to be isolated, or separated by ran-
dom jumps of wave frequency and phase as in most observations [Zhang, Agapi-
tov, Artemyev, Mourenas, Angelopoulos, Kurth, Bonnell and Hospodarsky, 2020;
Zhang, Mourenas, Artemyev, Angelopoulos, Kurth, Kletzing and Hospodarsky,
2020]. Compared with long packets, short wave-packets lead to a significantly
broader distribution of small positive energy changes at $\Delta E = 0 - 5$ keV through
phase trapping, due to a randomization of the duration of trapping and of the
related energy change, while the maximum energy changes are strongly reduced
(from 20 keV to 8 keV). Together with the randomization of phase bunching ef-
fects for negative energy changes, short wave-packets therefore lead to a more
symmetrical distribution of energy changes compared with long packets.

The periodical structures present in the distribution of energy changes due to
trapping for long packets ($\ell = 300$ and $\ell = 500$) is due to the discretization of
the trapping time: trapped particles can spend in acceleration regime time inter-
vals divisible to trapping periods $\sim 2\pi/\Omega_{tr}$. As a result, resonant particles may
makes $\sim N_{trap}$ times rotations during their trapping, with $N_{trap} = \Delta t_{trap}\Omega_{tr}/2\pi$.
Changes of this number by $\pm 1$ due to uncertainties of trapping/de-trapping time
result in the formation of several peaks in the $\Delta E$ distribution [see details in
Vainchtein et al., 2017].

Nonlinear anomalous phase trapping effects are also affected, because a very long
packet corresponds to one long phase trapping interval, whereas short wave-packets
result in multiple trapping occurrences, more probable but much more short-lived.
For medium to high pitch-angles, this leads to a modification of the probabil-
ity distribution of energy changes (see an example in Fig. 22(b)), whereas for
small pitch-angles this may prevent an effective particle transport away from the
loss-cone, i.e., it can reduce the efficiency of anomalous trapping [see Appendix
in Mourenas et al., 2022]. Figure 23 illustrates this idea: for long wave-packets a
significant fraction of low pitch-angle particles are transported by anomalous trap-
ping to higher pitch-angles [see also Kitahara and Katoh, 2019], whereas for short
wave-packets the resonant interactions become less regular and anomalous trap-
ping competes with phase bunching, leading to a global transport of initially low
equatorial pitch-angle electrons to both smaller and higher pitch-angles, in roughly
<figure><img src="image_23.png" /><figcaption>Fig. 23 Dynamics of electron equatorial pitch-angles due to a single resonant interaction with a plane wave, corresponding to an infinitely long ($\ell \to \infty$) wave-packet (a) and with a short ($\ell = 10$) wave-packet (b). The system parameters are the same as in Fig. 22.</figcaption></figure>

similar quantities (note that only the first cyclotron resonance is considered here), corresponding to a symmetrization of pitch-angle changes.

## 5.1 Kinetic equation

Let us consider the kinetic equation (19) under the short wave-packet approximation. The effects of phase trapping at arbitrary latitude (arbitrary $\gamma$ or $I$) and the variability of the phase bunching amplitudes cannot be straightforwardly included into this kinetic equation (we shall describe these effects using the mapping technique in Sect. 6.1). Thus, we limit our consideration to the effects in a small energy (or $I$) range with $V(\gamma) \propto S(\gamma) \neq 0$. The maximum energy change due to trapping is equal to $mc^2(\gamma_{max} - \gamma_{min})$, where $\gamma_{min,max}$ are boundary values of the $S(\gamma)$ profile, i.e., $S(\gamma_{min,max}) = 0$. Therefore, for a small range [$\gamma_{min}, \gamma_{max}$] we can expand the nonlocal term of Eq. (19) as:

$$f(\gamma') \approx f(\gamma) - \frac{\partial f}{\partial \gamma} \cdot (\gamma - \gamma') \approx f - 2 \frac{\partial f}{\partial \gamma} (\gamma - \gamma_0)$$

where $\gamma_0$ is the energy where $S(\gamma)$ reaches its peak value. Substituting this expansion into Eq. (19), we obtain

$$\frac{\partial f}{\partial t} = \frac{\partial}{\partial \gamma} \left( D_{\gamma\gamma} \frac{\partial f}{\partial \gamma} \right) - \tilde{V} \frac{\partial f}{\partial \gamma}, \quad \tilde{V} = V_{\gamma} + 2 \frac{dV_{\gamma}}{d\gamma} (\gamma - \gamma_0) \Theta \qquad (29)$$

$$\Theta = 1, \quad \gamma > \gamma_0, \quad \& \quad \Theta = 0, \quad \gamma < \gamma_0$$

Using the variable $x = (\gamma - \gamma_0)/(\gamma_{max} - \gamma_{min})$ with $\gamma_0 = (\gamma_{max} + \gamma_{min})/2$ and $S$ from Eq. (21), we can write

$$V_{\gamma} = -V_{0} (1 - x^{2})^{5/4}, \quad \tilde{V} = -V_{0} (1 - x^{2})^{1/4} (1 - x^{2} (5\Theta + 1))$$
Figure 24(a) shows the profile of $\tilde{V}(x)$: $\tilde{V}$ is positive for $x < 0$ ($\gamma < \gamma_0$, $\Theta = 0$) and negative for $x > 1/\sqrt{6}$ ($\gamma > 0$, $\Theta = 1$). Therefore, electrons with $x > 1/\sqrt{6}$ and $x < 1/\sqrt{6}$ are drifting in opposite directions, corresponding to the opposite effects of trapping-induced acceleration and deceleration by nonlinear scattering.

If we omit the diffusion term, $\sim D_{\gamma\gamma}$, Eq. (29) represents a simple drift equation

$$ \frac{\partial f}{\partial t} = -\tilde{V}(x) \frac{\partial f}{\partial x}, $$

and for an initial condition $f_{t=0}(x)$ this equation has an analytical solution: $f$ is a constant along trajectories $t - \int dx/\tilde{V}(x) = \text{const}$. The function $\int dx/\tilde{V}(x) \to \pm\infty$ at $x = 1/\sqrt{6}$ where $\tilde{V}(x) = 0$, implying that the particles should drift in energy space away from $x = 1/\sqrt{6}$. Indeed, the solution of the drift equation displayed in Figure 23(b) demonstrates the formation of a plateau in energy space around $x = 1/\sqrt{6}$. This effect is essentially the same as what occurs during the well-known quasi-linear evolution of a particle distribution function: resonant wave-particle interactions result in a flattening of the phase space density profile within the range of resonant energies [e.g., Galeev and Sagdeev, 1979; Vedenov and Ryutov, 1975]. However, in contrast with quasi-linear diffusion, this peculiar type of nonlinear interaction is much faster and corresponds to drifts $\sim \tilde{V} \propto (B_w/B_0)^{1/2}$. Figures 23(c,d) confirm this theoretical result of plateau formation within a narrow energy range for particles interacting resonantly with short wave-packets. These figures show the results of a numerical integration of Hamiltonian equations of motion for a large particle ensemble interacting with a monochromatic wave-packet (so that Eq. (4) is well satisfied and we may consider a 1D evolution of the particle distribution function).

Note that the location $x = 1/\sqrt{6}$ of the plateau formation shown in Figure 23(b) is largely determined by the system parameters, while for an ensemble of wave-packets (even if such packets are well separated in space and time) these plateaus will be formed around different $x$ (at different energies). Therefore, the global evolution of the electron distribution function can ultimately lead to a reduction of phase space density gradients over a wide energy range [see Section 6 and Artemyev, Neishtadt, Vasiliev and Mourenas, 2021, for details of such evolution].

## 5.2 Limit of nonlinear diffusion

Shorter wave-packets correspond to smaller energy changes due to trapping and a higher trapping probability (due to the high probability for electrons to resonate with wave-packet edges having a strong $dS/d\gamma \propto dB_w/d\lambda$ gradient). Therefore, as discussed earlier, the decrease of the wave-packet size should lead to a more symmetrical $\Delta E$ distribution (e.g., see Figs. 22(b) and 25, and the works from An, Wu and Tao [2022]; Frantsuzov et al. [2023]; Gan et al. [2022]; Tao et al. [2013]; Zhang, Agapitov, Artemyev, Mourenas, Angelopoulos, Kurth, Bonnell and Hospodarsky [2020]). The natural limit is the fully symmetrical $\Delta E$ distribution with $\langle \Delta E \rangle \sim 0$ and an equal amount of trapped and phase bunched particles. This regime of wave-particle resonant interactions is equivalent to diffusion, but the diffusion rate may not anymore linearly depend on wave intensity, contrary to the case of quasi-linear diffusion, where we have $\langle (\Delta E)^2 \rangle \propto B_w^2$ [Kennel and
<figure><img src="image_24.png" /><figcaption>Fig. 24 Panel (a) shows the profile of $\tilde{V}(x)/V_0$. Panel (b) shows the solution to Eq. (29) at several times $t$, with $t_0 = V_0 \cdot R/c$. Panels (c, d) show two examples of numerical solutions, using test particle tracing. System parameters are the same as in Fig. 2, but we use a wave amplitude $B_w = 100$ pT and a wave packet size $\ell = 15$. Two values of the constant $h$ from Eq. (4) are used: $h = 1.1$ (c) and $h = 1.3$ (d).</figcaption></figure>

Engelmann, 1966]. Let us examine this regime of nonlinear diffusion in more detail. We start with test particle simulations showing the $\Delta E$ distributions for different wave amplitudes $B_w$ and wave-packet sizes. Figure 25 shows $\Delta E$ distributions obtained from a numerical integration of Hamiltonian equations (2) with

$$B_w = B_{w,peak} \frac{1 - \exp(-A \sin^2(\phi/2\ell))}{1 - \exp(-A)}$$

where $B_{w,peak}$ is the peak wave amplitude, $\ell$ defines the number of wave oscillations (periods) within one wave packet, and $A$ controls the depth of amplitude modulations. We also introduce the effective wave amplitude, which corresponds to the wave intensity averaged over the wave-packet size:

$$B_w = \sqrt{\langle B_w^2 \rangle_{\phi \in [0, 2\pi\ell]}} = B_m \frac{\sqrt{1 - 2I_0(A/2) e^{-A/2} + I_0(h) e^{-A}}}{1 - e^{-A}}$$

where $I_n(z)$ is the modified Bessel function of the first kind.
Figure 25(a) shows the detailed distribution of individual (separate) electron energy changes $\Delta E$ for long wave packets. For small wave amplitudes, $B_w/B_0 < 5 \cdot 10^{-4}$, the distribution is symmetric relative to $\Delta E = 0$ with a dispersion scaling as $\langle (\Delta E)^2 \rangle \propto (B_w/B_0)^2$. This is the regime of purely diffusive electron scattering. As the wave amplitude increases, the $\Delta E$ distribution shows the formation of a population of phase trapped particles with $\Delta E \sim 10$ keV (well separated from the main distribution) and of a large population of phase bunched particles (increase of the probability for $\Delta E < 0$ in comparison with the probability of $\Delta E > 0$; the distribution boundary scales with wave amplitude as $\propto \sqrt{B_w/B_0}$). This is the standard picture of nonlinear resonant interactions with plane waves or very long wave-packets, when two distinct populations of phase trapped and phase bunched particles are formed, well separated in the $\Delta E$ space.

Figure 25(b) demonstrates the effect of including a wave-packet modulation, corresponding to electron interaction with a short wave-packet. The main difference with Fig. 25(a) is the absence of a distinct asymmetry in the $\Delta E$ distribution for large wave amplitudes, $B_w/B_0 < 5 \cdot 10^{-4}$. Although some asymmetry of $\Delta E > 0$ vs. $\Delta E < 0$ is present, the $\Delta E$ distribution is much more symmetric when a short wave-packet modulation is included, with a larger population of $\Delta E > 0$. Such symmetrization occurs because the wave amplitude modulation, or limited size, of a short wave-packet, leads to a randomization of the magnitudes and occurrences of nonlinear effects, such that phase trapping provides smaller and more random, but also more probable energy increases, while phase bunching provides more random energy decreases [see also An, Wu and Tao, 2022; Tao et al., 2013; Zhang, Agapitov, Artemyev, Mourenas, Angelopoulos, Kurth, Bonnell and Hospodarsky, 2020]. As a result, the main effect of a wave-packet modulation (i.e., of a limited wave-packet size) is a symmetrization of the $\Delta E$ distribution, such that wave-particle interactions become more diffusive-live, even in the limit of high wave amplitudes. Note that the dispersion of the $\Delta E$ distribution for large amplitudes scales with $B_w$ as $\langle (\Delta E)^2 \rangle^{1/2} \propto \sqrt{B_w/B_0}$, i.e., the diffusion by intense modulated wave-packets does differ from quasi-linear diffusion.

Let us briefly consider the limit of extremely short wave-packets, corresponding to fully incoherent waves. In such a case, the $\Delta E$ distribution can be characterized by its variance alone, and the only system characteristic is the diffusion coefficient. Note that this is only an ideal limiting case: it should not be confused with the case of usual short wave-packets that lead to both energy diffusion and distinct energy drifts due to phase trapping and bunching. This ideal limiting case is interesting, though, because we can derive the corresponding characteristic diffusion coefficient analytically, using Eq. (10), as the resonant scattering is then a purely local process (in contrast with phase trapping that depends on the wave-packet shape determining the $S(\gamma) \neq 0$ energy range). For such local processes, the wave-packet shape can be simply taken into account by introducing an effective wave amplitude parameter $B_w$. Equation $m_e c^2 \Delta \gamma = \omega \Delta I$ and (10) provide the energy change due to a single resonant interaction for an arbitrary wave intensity (exceeding the quasi-linear threshold) and, therefore, can be used to evaluate the variance of energy changes:

$$ (\Delta\gamma)^2 = \left( \frac{\omega\sqrt{2M|A|}}{m_e c^2} \right)^2 F^2(\xi, a) $$
<figure><img src="image_25.png" /><figcaption>**Fig. 25** Distributions of energy changes $\Delta E$ for different $B_{eff}/B_0$ with $A=0$ (a) and $A=1$ and $\ell=20$ (b) waves, with $\Delta E \propto B_w$ fits (white dashed lines). System parameters are $E_0 = 100$ keV, $\alpha_{eq} = 60^{\circ}$, L-shell = 6, whistler-mode waves with a frequency equal to 0.35 of the electron cyclotron frequency at the equator, and a constant plasma frequency equal to 10 of the electron cyclotron frequency at the equator. For each $B_{eff}$, we use $10^4$ trajectories to evaluate the $\Delta E$ distribution, each particle resonating with the wave only once, and $\Delta E$ being the energy change for a single resonance. Initial particle and wave phases are random and thus in the modulated amplitude case the actual wave amplitude is different for different particles having the same energy/pitch-angle. See details in [Frantsuzov et al., 2023].</figcaption></figure>

Figure 9(b) shows $\langle F^2 \rangle_\xi$ and its mean value, $\langle F \rangle_\xi$, evaluated as a function of parameter $a \propto B_w$. The variance has a peak at $a=1$, the boundary value between the weak wave approximation (no phase trapping and bunching for $a < 1$) and the intense wave approximation ($a > 1$). The diffusion coefficient for such intense waves takes the form

$$ \langle (\Delta\gamma)^2 \rangle_{\xi} = \left( \frac{\omega\sqrt{2M|A|}}{m_e c^2} \right)^2 \langle F^2(\xi, a) \rangle_{\xi} \quad (30) $$

and for very large wave amplitudes, $a \gg 1$, this coefficient has an asymptotic form [Frantsuzov et al., 2023]

$$ \langle (\Delta\gamma)^2 \rangle_{\xi} \approx \left( \frac{\omega\sqrt{2M|A|}}{m_e c^2} \right)^2 \frac{32a}{\pi^2} = \frac{64}{\pi^2} \left( \frac{\omega}{m_e c^2} \right)^2 MB \propto B_w $$

Therefore, the diffusion rate corresponding to intense monochromatic waves in the ideal limit of extremely short packets increases with wave intensity weaker than expected based on quasi-linear theory, where $D \propto B_w^2$. Figure 26 illustrates this effect: the diffusion rate scales with $B_w$ as $D \propto B_w^2$ up to the threshold for nonlinear resonant interactions, $B_w^*$, and then significantly decreases relative the $D \propto B_w^2$ expectation. This effect has been noted in numerical verifications of the quasi-linear diffusion equations, which demonstrated that the diffusion rate is lower than the expected quasi-linear level for fully incoherent waves of high intensity [Gan et al., 2022; Tao, Bortnik, Albert and Thorne, 2012]. These results underline
<figure><img src="image_26.png" /><figcaption>**Fig. 26** Schematic view of diffusion rate dependence on wave intensity. Intensity $B_w^*$ corresponds to the threshold for nonlinear resonant interactions, $a = 1$. See details in Frantsuzov et al., 2023.</figcaption></figure>

that even extremely short, strongly modulated wave-packets will not lead to a purely quasi-linear diffusion of electrons for a sufficiently high wave intensity: in such a case, the wave intensity cannot be used anymore as a simple scaling factor [i.e., the procedure of wave intensity averaging over spatial and temporal domains may affect the estimation of the diffusion rate; see also Watt et al., 2021, 2017, for discussion of a similar effect of temporal/spatial variability of diffusion rates].

A natural question arising from the results in Fig. 26 is whether quasi-linear diffusion models systematically overestimate diffusion rates during geomagnetically active periods, characterized by enhanced wave intensity. Previous simulations of radiation belt dynamics during geomagnetic storms rather suggest that quasi-linear diffusion underestimates, or estimates reasonably well, the rates of electron acceleration and losses [Allison and Shprits, 2020; Glauert et al., 2014, 2018; Li et al., 2014; Thorne et al., 2013]. Therefore, there should be some balance between, on one hand, the reduction of the electron pure diffusion rate at high wave amplitude due to the strong modulation of intense wave-packets [Frantsuzov et al., 2023] and, on the other hand, the coexistence of a faster, inherently nonlinear and advective electron transport provided by the main population of short and moderately intense waves-packets [Mourenas et al., 2018] and an even faster nonlinear electron acceleration due to a much smaller but finite population of intense long wave-packets [Hiraga and Omura, 2020]. Over the long run, such electron interactions with various wave-packets should result in a quasi-diffusive transport roughly similar to, or somewhat faster than, diffusion by the global time-averaged wave intensity, whereas over shorter time scales, the finite population of intense and not-too-short wave-packets can occasionally provide a much faster nonlinear resonant electron acceleration [Agapitov et al., 2015b; Foster et al., 2017]. In the next sections, we discuss how the net effect of intense wave-packets can be estimated in terms of diffusion.
# 6 Mapping technique for short wave-packets

The mapping technique is a more flexible approach than the kinetic equation (19), and this allows to make several generalizations for the inclusion of short wave packets. However, one of the main limitations for the description of short wave-packet systems, the absence of an analytical model for $S(\gamma)$, cannot be overcome anyway. Thus, we should either use some synthetically generated $S(\gamma)$ in mapping equations (25) or resort to a numerical integration of electron trajectories to determine the $S(\gamma)$ profiles.

## 6.1 Simplified mapping and typical scalings

We start our consideration with the limit of very short wave-packets, when $S(\gamma)$ can be approximated by Eq. (21), but the dynamics of wave-packets make it possible for electrons to interact resonantly with different wave amplitudes (and different $dB_w/d\lambda$ gradients) at any latitude. Note that the smallest wave-packet size (the shortest timescale of phase trapping) has been considered in Appendix C, while in the present section, we only consider short wave-packets with a size exceeding this minimum threshold size.

First, it is important to recall that nonlinear resonant interactions lead to significant energy drifts or jumps, and that such energy drifts and jumps cannot be accurately (fully) described within a diffusive framework, which intrinsically assumes very small and random individual energy changes, with zero mean value. Therefore, we focus here on a more general characteristic of the dynamic evolution of the electron distribution, common to both the quasi-linear and nonlinear regimes: the characteristic timescale of relaxation of the electron distribution toward its asymptotic steady-state. This asymptotic steady-state was indeed shown in Section 4.2 to be essentially the same state for both regimes. As will be demonstrated below, this common characteristic parameter will provide a way to approximately incorporate nonlinear effects in the diffusive framework.

Let us examine the characteristic timescale of the evolution of the electron distribution in a system with nonlinear wave-particle resonant interactions. Numerical and analytical solutions of Eq. (19) have shown that a plateau will ultimately be formed in the distribution function of electrons within the energy range of nonlinear resonant interactions, $S \neq 0$. This implies that the electron distribution will tend toward a solution identical to the asymptotic solution of the Fokker-Planck diffusion equation. Therefore, we may assume that such evolution can *on average* be approximated by a diffusive transport, albeit with an effective nonlinear diffusion coefficient $D_{NL}$ different from the quasi-linear coefficient $D_{QL} \propto B_w^2$. Such effective diffusion cannot describe all the peculiarities of nonlinear wave-particle interactions (e.g., the formation of transient electron populations at high-energies or the drift due to phase bunching), but it should provide a reasonable description of both the asymptotic (plateau) state of the electron distribution and of the characteristic timescale for reaching this asymptotic state. This effective nonlinear diffusion coefficient $D_{NL}$ should depend on wave intensity $B_w^2$ and wave-packet size $\ell$, and for very small $\ell$ this coefficient should tend toward the scaling $D_{NL} \propto B_w$ derived in Sect. 5.2.
### 6.1.1 Effective nonlinear diffusion

To develop the procedure of $D_{NL}$ derivation, let us start with the simple form of $S(x)$ given by Eq. (21), with a multiplication factor $\delta_0^{5/4}$ (see Appendix C). We introduce the parameter $\kappa$ through the equation $\delta_0 = (B_w/B_0)^{4\kappa/5}$. This parameter $\kappa$ is related to the maximum number of trapped particle periods, $N_{trap}$, during resonant interaction with a wave-packet, and a larger $\kappa$ corresponds to a smaller $N_{trap}$ (see Appendix C, and here below). We rewrite Eq. (21) as

$$S = (B_w/B_0)^{1/2+\kappa} (1-x^2)^{5/4} \qquad (31)$$

Figure 27 shows examples of trajectories given by this mapping and examples of the evolution of the particle distribution function $f(x)$. For $\kappa = 0$, the trajectory is very similar to the trajectory displayed in Fig. 16(a), but for $\kappa = 1/3$ we have smaller $S \propto (B_w/B_0)^{1/2+\kappa} = (B_w/B_0)^{5/6} \ll (B_w/B_0)^{1/2}$, which results in smaller $x$ changes due to bunching (longer intervals of $x$ decrease) and a smaller probability of trapping (rarer jumps with $x$ increase). Despite these differences, the evolution of the particle distribution $f(x)$ is quite similar, but it takes more time for $\kappa = 1/3$. Figure 27(b,c) shows that after 100 resonant interactions (1000 for $\kappa = 1/3$) the initially localized $f(x)$ is transformed into a uniform distribution, in agreement with the asymptotic solution derived in Sect. 4.2 (see also Fig. 24(b)). Therefore, $\kappa$ regulates the time of $f(x)$ evolution, but does not change the main features of this evolution and its asymptotic solution.

To characterize the timescale of the evolution of $f(x)$ as a function of $\kappa$, we numerically integrate a large ensemble of trajectories described by the map (31), for various values of $\kappa$. The initial $x_{0,i}$ values for the $i = 0...N$ trajectories are uniformly distributed within the $[-1, 1]$ range, and we calculate two characteristics:

$$M_1(n) = N^{-1} \sum_{i=0..N} (x_{n,i} - x_{0,i}) \quad M_2(n) = N^{-1} \sum_{i=0..N} (x_{n,i} - x_{0,i})^2 - M_1^2(n)$$

where $x_{0,i}$ are initial $x_i$ values. Figure 28(a) shows $M_1$ and $M_2$ profiles for two $\kappa$ values: $M_1$ is around zero and $M_2$ first increases with the number of resonant interactions, and then saturates. This saturation means a full mixing of particles within the $S(x) \neq 0$ domain (i.e., the formation of a plateau in the corresponding $f(x)$ distribution, without any $df/dx$ gradient). Before the saturation, the growing fragment of $M_2$ can be fitted as

$$M_2(n) = Qn^q$$

and if $q \approx 1$, we may interpret the proportionality coefficient $Q$ as a diffusion rate $D_{NL}$, because this coefficient describes the rate of increase of particle variance $dM_2/dn$ during particle mixing. Note that, instead of being a local diffusion rate (that should depend on $x$), the coefficient $D_{NL}$ determined from the $M_2(n)$ profile is a global system characteristics that is averaged over the entire $x$ range of $S(x) \neq 0$. Because the particle dynamics include very large $x$ jumps (due to phase trapping; see Fig. 27(a)), this dynamics cannot be described locally by a diffusion, but the overall particle mixing within the entire $x$ range stills follows a diffusive-like behavior. Therefore, $D_{NL}$ may be considered as a global characteristic of nonlinear wave-particle interactions, corresponding to the inverse characteristic timescale of
<figure><img src="image_27.png" /><figcaption>Fig. 27 Panel (a) shows two examples of $x_n$ trajectories described by map (26) with $S$ given by Eq. (31). Panels (b,c) show evolution of $f(x)$ distribution for two $\kappa$ values. Results are shown for $B_w/B_0 = 10^{-3}$.</figcaption></figure>

evolution of the electron distribution, and it can be compared with average quasi-linear diffusion rates. Figure 28(b) shows the $D_{NL}$ scaling with $\kappa$ and $(B_w/B_0)$: $D_{NL} \propto (B_w/B_0)^{1/2+\kappa}$, i.e., it is proportional to $S$ magnitude.

To explain the scaling $D_{NL} \propto (B_w/B_0)^{1/2+\kappa}$, let us consider an ensemble of $N$ electrons with the same initial $x$ (i.e., with the same initial energy and $h$). The distribution of $\Delta x$ changes consists of two well separated populations: the most representative population contains almost all particles, $\sim N$, and have a finite $\langle x \rangle \propto S \propto (B_w/B_0)^{1/2+\kappa}$; this is the phase bunched population. A much smaller electron population contains only a tiny fraction $(B_w/B_0)^{1/2+\kappa}$ of all resonant particles, $\sim NII \sim N(B_w/B_0)^{1/2+\kappa}$ where $\Pi = dS/dx \propto S$. This population is characterized by a finite $\langle \Delta x \rangle \sim O(B_w/B_0)$. Such a large $x$ change is due to phase trapping. The dispersion of the entire distribution of $\Delta x$ is about

$$ \langle (\Delta x)^2 \rangle \approx \langle \Delta x \rangle_{\text{bunching}}^2 \frac{N_{\text{bunching}}}{N} + \langle \Delta x \rangle_{\text{trapping}}^2 \frac{N_{\text{trapping}}}{N} \\ \sim \left(\frac{B_w}{B_0}\right)^{1+2\kappa} \cdot O(1) + O(1) \cdot \left(\frac{B_w}{B_0}\right)^{1/2+\kappa} \propto \left(\frac{B_w}{B_0}\right)^{1/2} $$
<figure><img src="image_28.png" /><figcaption>Fig. 28 Panel (a) shows $M_{1,2}$ profiles for two $\kappa$ values and $B_w/B_0 = 10^{-3}$. Panel (b) shows $D_{NL}$ dependence on $\kappa$; fitting $D_{NL} \propto (B_w/B_0)^{1/2+\kappa}$ is shown by black lines.</figcaption></figure>

The parameter $\kappa$ (or equivalently, $\delta_0$) can be expressed through a more physical variable – the number of oscillations of phase trapped particles around the resonance, $N_{trap}$. For a small $S$ in Eq. (31), the period of such oscillations can be written as $N_{trap} \sim (B_w/B_0)^{3\kappa/5-1/2}$ (see Appendix C). We use the normalized number of trapping periods $N_{trap} = N_{trap}/N_{trap}^{\infty}$, where $N_{trap}^{\infty} \sim O(1)/\Omega_{tr} \sim (B_w/B_0)^{-1/2}$ is the number of periods in long wave-packets, where trapping and detrapping are determined by the latitudinal gradients of the background magnetic field and of the wave amplitude. Therefore, $N_{trap}/N_{trap}^{\infty} \sim (B_w/B_0)^{3\kappa/5}$ and $(B_w/B_0)^{\kappa} \sim N_{trap}^{5/3}$. The effective diffusion rate $D_{NL}$ can be then written as

$$D_{NL} \sim (B_w/B_0)^{1/2+\kappa} \sim (B_w/B_0)^{1/2} N_{trap}^{5/3}$$

For $N_{trap} \to 1$, in the limit of long wave-packets, the effective nonlinear diffusion rate is $D_{NL} \sim (B_w/B_0)^{1/2} \propto B_w^{1/2}$. The limit of intense but short wave-packets with $D_{NL} \sim (B_w/B_0) \propto B_w$ can be reached for $N_{trap} \sim (B_w/B_0)^{3/10}$ (or $N_{trap} \sim (B_w/B_0)^{-1/5}$). It is worth emphasizing that the short wave-packet approximation (or, in another words, the approximation of a symmetrical distribution of $\Delta x$ changes; see Fig. 25) provides a full analytical description of the transition from the classical quasi-linear diffusion rate $D_{QL} \propto B_w^2$ to the strongly nonlinear effective diffusion rate $D_{NL} \propto B_w$ (such that $D_{NL} > D_{QL}$), including all the numerical factors quantifying the diffusion rates. Note also that the effective nonlinear diffusion rate $D_{NL} \propto B_w$ produced by intense short wave-packets (which represent the overwhelming majority of the observed chorus wave-packets in the Earth's outer radiation belt) is larger than $D_{QL}$ because it includes a finite, important contribution from nonlinear interactions which hasten the evolution of the electron distribution, contrary to the ideal case of extremely short packets, briefly investigated in Section 5.2, where phase trapping and bunching strictly compensate each other, leaving only a purely diffusive transport at a rate lower than $D_{QL} \propto B_w^2$ even in the high $B_w$ range.
The effective nonlinear diffusion rate for long wave-packets, $D_{NL} \propto B_w^{1/2}$ has a very large magnitude (see numerical test particle simulation results in Section 7), allowing $D_{NL} \propto B_w^{1/2}$ to be much higher than $D_{QL} \propto B_w^2$ over a wide range of $B_w$ values. In a real system containing a large distribution of wave-packets of various sizes and peak amplitudes, the properties of this distribution will determine the actual magnitude of the effective nonlinear diffusion rate $D_{NL}$. However, because the scaling of the larger $D_{NL}$ with $B_w$ in the case of (rare) long packets is slower ($D_{NL} \propto B_w^{1/2}$) than for the weaker $D_{NL}$ due to short packets ($D_{NL} \propto B_w$), which is itself increasing more slowly with $B_w$ than the even weaker quasi-linear diffusion rate $D_{QL} \propto B_w^2$, it is no surprise that the total effective nonlinear diffusion rate may be close to the quasi-linear rate $D_{QL} \propto B_w^2$ extrapolated into the high $B_w$ range. This may explain why quasi-linear simulations often manage to reproduce the observed electron flux dynamics in the radiation belt [Glauert et al., 2018; Li et al., 2014; Ma et al., 2018; Mourenas, Artemyev, Agapitov, Krasnoselskikh and Li, 2014; Thorne et al., 2013] even when the measured wave intensity exceeds the threshold for nonlinear wave-particle interactions [Zhang, Mourenas, Artemyev, Angelopoulos, Bortnik, Thorne, Kurth, Kletzing and Hospodarsky, 2019; Zhang, Thorne, Artemyev, Mourenas, Angelopoulos, Bortnik, Kletzing, Kurth and Hospodarsky, 2018].

### 6.1.2 Diffusion by multiple short wave-packets

The mapping (31) can successfully describe electron resonant interactions with a short wave-packet, but this interaction is assumed to take place within the entire energy range where $S \neq 0$. This is a simplified case. In a more realistic situation, we deal with multiple short wave-packets. At a given time, each of these wave-packets covers only a small range of magnetic latitudes and, thus, will interact resonantly with electrons such that their energy is comprised within only a small portion of all the resonant energies at these latitudes. We consider a bounce period as an elementary time-step between electron resonant interactions that change their energy and pitch-angle. After each time-step, wave-packets will have significantly propagated, and each electron resonant interaction will occur with a new distribution of wave-packets along the magnetic field line. Therefore, we should include into the model (1) the probability for each electron to meet one wave-packet at the resonant latitude, (2) the corresponding limited range of resonant energies for each wave-packet. This limited energy range of actually resonant electrons represents a fraction $\chi < 1$ of the total energy range of particles potentially reaching resonance with a wave packet over the full length of the magnetic field line.

Electron trapping by a short packet leads to a smaller energy change than trapping by an ideal *infinitely long* wave packet, because the electron is released from trapping faster, corresponding to a reduction of its energy change by a factor $\sim \chi$. If short wave packets are rare (not appearing in close succession) and occur only approximately once every bounce period, this corresponds to an additional reduction of the occurrence rate of resonant interaction by a factor $\sim \chi < 1$, as compared with both the case of a close succession of short packets and the ideal case of *infinitely long* packets assumed in the preceding subsection.

Accordingly, let us consider a realistic situation where the magnitude of $S \sim (B_w/B_0)^{1/2}$ remains the same, but the range of nonlinear resonant interactions
shrinks. There are two kinds of systems with an total resonant range $x \in [-1, 1]$, but where $S$ is not equal to zero only for $x \in [-\ell, \ell]$, with $\chi = \sqrt{\delta_0}$. In systems of the first kind, we center $S$ around a $\bar{x}_n$ randomly generated at each map iteration around $x_n$, i.e., for each iteration there is a finite change of $x_n$. This type of mapping mimics electron resonant interactions with a set of short wave-packets that would fill the entire magnetic field lines, but which cannot trap electrons for a long time because of the limited packet duration. Figure 29(a) shows observational examples of such a dense filling of the field-line by short wave packets.

In the system of the second kind, we center $S$ around a randomly generated $\bar{x}_n$, but we do not control the position of $\bar{x}_n$ relative to $x_n$ and some map iterations can occur without any change of $x_n$ because $x_n$ is outside the $S \neq 0$ range. This type of mapping mimics electron resonant interaction with rare short wave-packets propagating with a significant time separation. It corresponds to a situation where, during each bounce period, only one wave packet is present and many electrons reach the latitude of cyclotron resonance without encountering this intense packet there. Figure 29(b) shows observational example of such rarefied filling of a magnetic field-line by short wave packets.

For both systems, the mapping (31) can be rewritten as

$$x_{n+1} = x_n + \begin{cases} \bar{x}_n - 2x_n, & \xi \in [0, \Pi] \\ -S_n(x_n), & \xi \in (\Pi, 1] \end{cases} \quad (32)$$

and

$$S_n = (B_w/B_0)^{1/2} (\chi^2 - (x - \bar{x}_n)^2)^{5/4} \\ \Pi = \frac{dS_n}{dx} = -\frac{5}{2} (x - \bar{x}_n) (B_w/B_0)^{1/2} (\chi^2 - (x - \bar{x}_n)^2)^{1/4}$$

where $\bar{x}_n = x_n + \chi \cdot R$ for the system of the first kind, and $\bar{x}_n = -1 + \chi + 2(1 - \chi) \cdot R$ for the system of the second kind; here $R$ is a random number with a uniform distribution within $[-1, 1]$ (note that for $\bar{x}_n = y_n + \chi \cdot R$ we also control that $S \neq 0$ range does not cross $x \pm 1$ boundaries). An important property of both systems is that for small $\chi$, trapping is possible for any $x_n$ values, whereas for $\chi = 1$ (the initial map given by Eq. (31)) trapping is possible only for $x_n < 0$. Figure 29(c,d) shows a set of sample trajectories for each of the two kinds of systems. The rate of change of $x$ is going down as $\chi$ decreases, and this effect is stronger for the system of the second kind (Fig. 29(d)).

For these two kinds of systems, we set the range of $\chi$ and for each $\chi$ value we calculate $M_2(n)$. Then we fit the growing fragment of $M_2(n)$ by $D_{NL}n$. Fig. 29(e) shows the corresponding $D_{NL}$ dependence on $\chi$. For both systems, $D_{NL}$ scales with $\chi$ as $D_{NL} \sim (B_w/B_0)^{1/2} \chi^n$. But $\eta \approx 7/2$ for the system of the first kind, where resonances occur at each iteration, whereas $\eta \approx 9/2$ for the system of the second kind, where resonance occurrence is decreased by a factor $\chi \ll 1$. Taking into account that $\chi = \sqrt{\delta_0}$, we can write $N_{trap} \sim \chi^{3/2}$. Therefore, for the scaling of $D_{NL}$ as a function of $N_{trap}$ we have $D_{NL} \sim (B_w/B_0)^{1/2} \chi^n \sim (B_w/B_0)^{1/2} N_{trap}^{2\eta/3}$ and $D_{NL} \sim (B_w/B_0)^{1/2} N_{trap}^{2\eta/3}$. For $\eta = 5/2$, this scaling allows to recover the scaling of the ideal case of electron interactions with infinitely long wave packets/waves present all the time along magnetic field lines: $D_{NL} \sim (B_w/B_0)^{1/2} N_{trap}^{5/3}$. In the more realistic situation of short wave packets, two different kinds of systems are
possible, corresponding to $\eta = 7/2$ and $\eta = 9/2$. For the system of the first kind, we get $D_{NL} \sim (B_w/B_0)^{1/2} N_{trap}^{7/3}$. For the system of the second kind, we have $D_{NL} \sim (B_w/B_0)^{1/2} N_{trap}^3$. The parameter $N_{trap}$ depends on wave characteristics (dispersion, wave-packet size, and wave magnitude) and resonant electron characteristics. Although simplified scalings of this parameter with wave and electron characteristics can be derived analytically (see Figs. 28 and 29 and [Artemyev, Neishtadt, Vasiliev and Mourenas, 2021]), a more accurate (likely numerical) approach is needed to complete the quantification of the diffusion regimes for different wave amplitudes and wave coherence. Beside a fully numerical approach through test particle simulations, there are two promising directions for solving the problem of the dynamics of a large ensemble of resonant electrons. The first direction requires a significant modification of the mapping technique, with inclusion of realistic distributions of $B_w^2$ and $\ell$. We discuss this approach in Section 6.2. The second approach requires incorporating a mixture of quasi-linear and nonlinear resonant effects, weighted by the $P(B_w, \ell)$ distribution, into the Fokker-Planck diffusion equation; we discuss this second approach in Section 7.

## 6.2 Synthetic map

The mapping (32) is a simplified version of the realistic map describing electron resonant interaction with multiple different wave-packets. Although this simplified version does reproduce many peculiarities of such interactions, we may still generalize it to incorporate the observed distribution of wave characteristics, $P(B_w, \ell)$, as in examples shown in Figure 30(b) for field-aligned whistler-mode waves. This distribution shows the prevalence of short packets, but also demonstrates the existence of very long packets with $\ell > 100$. Note that the method of wave-packet determination (i.e., the criterion used for identifying the wave-packet edges) can significantly affect the shape of the $P(B_w, \ell)$ distribution [see results in Zhang et al., 2021; Zhang, Mourenas, Artemyev, Angelopoulos, Bortnik, Thorne, Kurth, Kletzing and Hospodarsky, 2019]. Wave-packets for very oblique whistler-mode waves are generally shorter (compare panels (a) and (b) of Fig. 30) because their comparatively lower magnetic fields (corresponding to very high electric field amplitudes) are often closer to the noise level. The dimension of wave characteristics can be further increased by including wave frequency, $P(B_w, \ell) \to P(B_w, \ell, \omega/\Omega_0)$. Each bin of this distribution corresponds to a specific map, but in contrast to the $\ell \to \infty$ limit where such a map can be evaluated analytically (see Section 4.4), for an arbitrary $\ell$ there is no analytical model of wave-particle resonant interactions. Therefore, this approach should rely on test-particle simulations, and the results of such simulations should be either incorporated into the Green function approach [Hsieh et al., 2020, 2022], or somehow used for an improved mapping.

The idea of incorporating such wave-packet distributions into the mapping technique has been proposed in [Shi et al., 2023; Tonoian et al., 2023], and can be named the synthetic map approach. This approach consists in the construction of a combined map

$$ \hat{M} = \sum_{B_w, \ell, \omega} P(B_w, \ell, \omega) \hat{m} \qquad (33) $$
<figure><img src="image_29.png" /><figcaption>**Fig. 29** (a,b) show observational examples of dense and rarefied filling of a magnetic field line by short wave-packets. Data are collected by THEMIS search-coil magnetometer [Le Contel et al., 2008]. Details of these wave events can be found in [Zhang, Agapitov, Artemyev, Mourenas, Angelopoulos, Kurth, Bonnell and Hospodarsky, 2020]. Panels (c,d) show examples of trajectories $x(n)$ for the two types of models of mapping (32). Panel (e) shows the corresponding scaling of $D_{NL} = dM_2/dn$ with $\chi$ for the two models of mapping.</figcaption></figure>
<figure><img src="image_30.png" /><figcaption>Fig. 30 Distributions of wave-packets amplitudes and sizes for very oblique (a) and field-aligned (b) whistler-mode waves (L-shell∈ [6, 10] for very oblique waves, THEMIS statistics; L-shell∈ [4, 6] for parallel waves, Van Allen Probe statistics [see details of dataset in Zhang, Mourenas, Artemyev, Angelopoulos, Bortnik, Thorne, Kurth, Kletzing and Hospodarsky, 2019]). All levels of geomagnetic activity are included.</figcaption></figure>

where $\hat{m}$ is the mapping described by Eq. (32) with the

$$S_{B_w, \ell, \omega} = C_{B_w, \ell, \omega} (B_w/B_0)^{1/2} (\chi^2 - (x_n - x)^2),$$

where $C_{B_w, \ell, \omega}$, $\chi$, and the probability distribution of $P(x_n)$ have to be determined through test particle simulations. Such simulations (for fixed $B_w, \ell, \omega/\Omega_0$) provides the probability distribution function of particle energy changes, and the main elements of $S_{B_w, \ell, \omega}$ model can be derived from this distribution – see schematic view in Fig. 31(a). The maximum energy change due to trapping determines the parameter $\chi$, which should be proportional to the wave-packet size parameter $\ell$ (see Fig. 22(b)). The magnitude of negative energy changes due to phase bunching determines the parameter $C_{B_w, \ell, \omega}$, which denotes the normalized magnitude of the area $S$. The relative number of test particles such that their energy does not change in the first resonant interaction, due to an absence of wave-packet at the latitude of resonance, determines the distribution $P(x_n)$ [see details in Tonoian et al., 2023]. Therefore, even a quite roughly reproduced (due to limited statistics from test particle simulations) probability distribution of electron energy changes may allow to reproduce the main mapping characteristics and to construct a synthetic map. Figure 31(c) shows an example of such a synthetic map, constructed to reproduce the numerically obtained distribution of electron energy changes from Fig. 31(c). Although the synthetic map is simplified and omits many details of this distribution, the application of this map for the simulation of the dynamics of an ensemble of electrons shows very similar results to what can be obtained directly from test particle simulations (e.g., compare panels (d) and (e) in Fig. 31). The general dynamics of the electron energy distribution for this simulation shows an increase of energetic particle population due to trapping, but because we deal with an ensemble of short wave-packets, the formation of the plateau in the distribution function takes a much longer time than for simulations with long wave-packets.

The method of construction of a synthetic map resembles a general decomposition technique where the results of multiple different processes (bunching and
<figure><img src="image_31.png" /><figcaption>**Fig. 31** Panel (a) provides a schematic view of the construction of the $\hat{m}$ mapping using the numerically evaluated distribution of energy changes. Note that we normalize electron energy to the minimum energy, $E_{min}$, of the range where the nonlinear resonant interactions are possible. Panels (b), (c) show the numerically determined 2D probability distribution $P(E/E_{min}, \Delta E/E_{min})$, as fitted by a sum of $\hat{m}$ maps. Panels (d), (e) show the evolution of the electron distribution function obtained directly from test particle simulations and, indirectly, using the mapping technique with the synthetic map shown in panel (c). Details of this demonstration of the usefulness of the synthetic map can be found in [Tonoian et al., 2023].</figcaption></figure>

trapping by wave-packets having different amplitudes and sizes) is expressed as a linear sum of elementary processes, each of them being described by a simple model (e.g., the mapping model (32)). Such a decomposition requires a minimization of the difference between $P(\Delta E)$ derived numerically and obtained from Eq. (33). Numerical schemes for this minimization are not yet well developed ([Tonoian et al., 2023] used a quite simplified form of determination of the parameters of Eq. (32) shown in Fig. 31) and, therefore, the synthetic map approach has not been thoroughly investigated yet. However, this approach looks quite promising for complicated (realistic) systems, and it would be worth generalizing and developing it in the future.
## 6.3 SDE versus mapping

The mapping approach utilizes the probability distribution function of electron resonant energy changes, $P(ΔE)$, and thus, it is quite similar to the more widespread stochastic differential equation (SDE) approach, which is mostly used for simulations of the purely diffusive regime of wave-particle interactions. The basic idea of SDE is that, instead of directly solving the Fokker-Planck equation, one can write and solve the corresponding Ito stochastic differential equations of quasi-particle trajectories [Tao et al., 2008; Zheng et al., 2014]:

$$p(t + Δt) = p(t) + μ(p(t), t)Δt + σ(p(t), t)dW_t \quad (34)$$

where $Δt$ is the time step over which we calculate the change of $p = (p_∥, p_⊥)$, $μ(p(t), t)$ is a 2-dimensional vector of the drift coefficient, $σ(p(t), t)$ is a 2 × 2-dimension matrix related to the diffusion coefficients written in such a way that $\hat{D}_{pp} = \frac{1}{2}σσ^T$, $\hat{D}_{pp}$ is a 2 × 2 matrix of diffusion coefficients[Albert, 2018; Lyons and Williams, 1984; Schulz and Lanzerotti, 1974], $W_t$ is a N-dimension standard Wiener process; $dW_t = \sqrt{Δt}N$, where N is a vector of standard normal random values, $N_i \sim N(0, 1)$.

The term quasi-particles means that we do not directly integrate the equations of motion, but treat the change of $p$ as a stochastic process and approximate it by the equation (34). As a result, two quasi-particles having equal initial conditions $p_0$ may have different trajectories $p(t)$ (in the numerical integration of the equation (34), one can fix the seed of the pseudo random number generator to preserve the sequence of random numbers and make the results repeatable). We can examine electron distributions in the energy space, and thus the Ito equation (34) can be rewritten in the following form

$$E(t + Δt) = E(t) + μ_E(E(t)) Δt + \sqrt{2D_{EE}}(E(t))dW_t \quad (35)$$

This equation describes energy evolution for a fixed $h$ given by Eq. (4), i.e., for a monochromatic wave we may reduce the energy, pitch-angle evolution to the energy evolution only and calculate pitch-angle changes from $h$ conservation. In the more general case of a wave spectrum, 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 energy and pitch-angle evolution [see detail in Tao et al., 2008]. Note also 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 and Horne, 2005], but the Ito equation will still have the same form. Equation (35) is quite similar to the mapping equations (compare with Eqs. (26)), but instead of a more complicated probability distribution function containing phase trapped and bunched populations, this equation utilizes a Gaussian distribution of random energy jumps.

In this Section, we shall compare 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 for small wave amplitudes, when wave-particle resonant interactions are quasi-linear. The test particles approach is expected to be more precise, because it does not rely on a constructed $ΔE$ distribution and it is based on the full set of equations of motion. The main
<figure><img src="image_32.png" /><figcaption>**Fig. 32** Two examples of $P(ΔE)$ distributions for low-amplitude waves with $B_w/B_0 = 10^{-4}$ (a) and intense waves with $B_w/B_0 = 10^{-3}$ (b). In panel (b) the red insert zooms in on $P(ΔE)$ for trapped particles. System parameters are: $ω/Ω_0(0) = 0.3$, $L = 5$, $h = 1.5$, field-aligned waves [see details in Lukin et al., 2024].</figcaption></figure>

advantage of SDE and mapping techniques is their computational efficiency in long-term simulations. Therefore, one can accept a lower accuracy of SDE and mapping, provided that they still describe the main features of the evolution of the electron distribution and that their results statistically reproduce the results from direct test particle simulations. Figure 33(a) shows the evolution of the electron distribution function for a small wave amplitude ($B_w/B_0 = 10^{-4}$): we show the initial distribution and the distributions obtained by three methods after $\sim 100$ bounce periods (see figure caption for details). Without nonlinear resonant interactions, both SDE and mapping technique show results consistent with the test particle simulation. In this case, the evolution is diffusive (as expected from 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 \sim 420$ keV, is due to an overestimation of the diffusion coefficients. We evaluate diffusion coefficients as half of the variance of $ΔE$ distributions, and thus we assume that $ΔE$ distributions are symmetric relative to their mean value. However, even in the case of low-amplitude waves (see, e.g., Figure 32(a)) this assumption may not work, leading to an overestimation of the diffusion rate. The mapping technique does not require any assumptions about $ΔE$ distributions and, therefore, it performs better even in the case of low-amplitude waves.

Figure 33(a) shows that for small wave amplitudes, all three methods, SDE, mapping, and test particle simulations, provide the same evolution of the electron distribution function. This results underlines that the mapping technique is equivalent to SDE approach for the quasi-linear regime. Moreover, for such amplitudes the details of the $P(ΔE)$ distribution are not important: despite the fact that the distributions $P(ΔE)$ from Fig. 32(a) are not Gaussian, their dispersion well describes the electron dynamics due to a very large number of resonant interactions. This is a consequence of the *Central limit theorem*, which states that the mean
<figure><img src="image_33.png" /><figcaption>**Fig. 33** Two examples of evolution of electron distribution function, evaluated using different techniques, for low-amplitude waves with $B_w/B_0 = 10^{-4}$ (a) and for intense waves with $B_w/B_0 = 10^{-3}$ (b). System parameters are given in the caption of Fig. 33, the typical bounce period of electrons is $\tau_b \approx 0.7$s, and results are shown after 1100 s of real time [see details in Lukin et al., 2024].</figcaption></figure>

value of weakly correlated variables is normally distributed as long as the sample size is large enough, whatever the probability distribution of these variables [Dudley, 2014; Elliott, 1980]. Therefore, we should not care too much about describing the fine details of the probability distributions of energy changes, but rather focus on their main characteristics ($\Delta E$ range, mean value, etc.).

Figure 33 (b) shows the evolution of the electron distribution function for a large wave amplitude ($((B_w/B_0) = 10^{-3})$: we show the initial distribution and the distributions obtained by test particle simulations and by the mapping technique after ~ 100 bounce periods (see figure caption for details). For such intense waves, the SDE approach becomes fully inapplicable, but we can still compare the results from test particle simulations and from the mapping technique. The mapping technique accounts for nonlinear resonant interactions (e.g., phase trapping) and describes well the evolution of electron distribution. After several wave-particle resonant interactions, 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 toward smaller energy, the probability of particle trapping increases (see Figure 32(b)) and more particles become trapped and accelerated. Accelerated particles are transported to latitudes of resonance where phase trapping is not possible anymore, and thus, these particles start losing their energy due to phase bunching. Around the time when the main population (at the initial peak of electron phase space density) reaches the left boundary of the allowed energies, the processes of phase bunching and phase trapping statistically compensate each other. This results in the formation of a plateau in the distribution function; this moment is shown in Fig. 33(b). Such an evolution of the electron distribution is consistent with theoretical predictions for the system with multiple nonlinear resonances (see Section 4.2 and [Artemyev, Neishtadt and Vasiliev, 2019]).