2. Basic concept

To propose the approach for the evaluation of such diffusion rates, let us illustrate the wave modulation effect on nonlinear wave–particle interactions. We consider electrons bouncing in a magnetic trap modelled by a curvature-free dipole field (Bell Reference Bell1984) and their resonant interaction with a monochromatic and intense whistler-mode wave. To evaluate a set of test particle trajectories resonating once with whistler-mode waves, we use the approximation of a monochromatic field-aligned wave. The wave field distribution along the magnetic field lines and the concept of description of wave packets are taken from Mourenas et al. (Reference Mourenas, Zhang, Nunn, Artemyev, Angelopoulos, Tsai and Wilkins2022). We start with the Hamiltonian of a relativistic electron (rest mass is $m_e$ and charge is $-e$) bouncing in a magnetic trap and interacting with a field-aligned whistler-mode wave:

(2.1)\begin{equation} H = \sqrt{\left({\boldsymbol p} + \frac{e}{c} {\boldsymbol A}\right)^2 c^2 + m^2 c^4}, \end{equation}

where ${\boldsymbol p}$ is a canonical momentum and ${\boldsymbol A}$ is a vector potential that can be derived from ${\boldsymbol B} = \boldsymbol {\nabla }\times {\boldsymbol A}$ with ${\boldsymbol B} = {\boldsymbol B}_0 + {\boldsymbol B}_w$. Here ${\boldsymbol B}_0$ is the background magnetic field of Earth's dipole and ${\boldsymbol B}_w$ describes the wave field. As the electron gyroradius is significantly smaller than the dipole magnetic field line curvature, $\sim LR_E$ where $R_E$ is Earth's radius and $L$-shell, we consider a curvature-free magnetic field with a pair $(z, p_z)$ of Cartesian coordinates and momentum playing a role of field-aligned coordinate and momentum $(s,p_{\parallel })$ (see Bell Reference Bell1984). The equatorial magnetic field is determined by a dipole model, $B_0(0)\propto L^{-3}$. To define a $B_0(z)$ function, we introduce a geomagnetic latitude $\lambda$:

(2.2a,b)\begin{equation} \frac{{\rm d}z}{{\rm d}\lambda} = \sqrt{1+3 \sin^2\lambda} \cos\lambda,\quad \frac{B_0(\lambda)}{B_0(0)} = \frac{\sqrt{1+3 \sin^2\lambda}}{\cos^6\lambda}. \end{equation}

The smallness of electron gyroradius allows us to write ${\boldsymbol A}_0=x\tilde {B}_0(z) {\boldsymbol e}_y$ (${\boldsymbol B}_0 = \boldsymbol {\nabla }\times {\boldsymbol A}_0$), where ${\boldsymbol e}_{r_i}$ is the unit vector along the $r_i$ axis, $r_i = (x, y, z)$, and $x$ is the cross-field coordinate. As the magnetic field ${\boldsymbol B}_0$ is mainly oriented along the $z$ axis, we use the approximation $\tilde {B}_0(z) \approx B_0(z)$. To derive the equation for the wave vector potential ${\boldsymbol A}_w$, we introduce the wave phase $\phi$ and define the magnetic vector ${\boldsymbol B}_w$ as $B_w {\boldsymbol e}_x \cos \phi - B_w {\boldsymbol e}_y \sin \phi$. Note that $B_w$ may be a function of $z$ (or/and $\phi$), but for the following derivations we assume that $|{\rm d}B_w/{\rm d}z| \ll |(\partial \phi /\partial z)B_w|$. The wave frequency $\omega$ and the wave vector ${\boldsymbol k}$ are defined by equations $\omega = -\partial \phi /\partial t$ and ${\boldsymbol k} = \boldsymbol {\nabla } \phi \approx (\partial \phi /\partial z) {\boldsymbol e}_z$ for field-aligned waves. We consider $\omega ={\rm const.}$ and determine ${\boldsymbol k}$ from a cold plasma dispersion: $k c / \varOmega _{pe} = ( \varOmega _{ce}/\omega - 1)^{-{1}/{2}}$, where $\varOmega _{pe} = \sqrt {4 {\rm \pi}n_e e^2 / m_e}$ is the plasma frequency and $\varOmega _{ce} = e B_0 / m_e c$ is the cyclotron frequency (Stix Reference Stix1962). Therefore, ${\boldsymbol A}_w\approx B_w {\boldsymbol e}_x \cos \phi / k - B_w {\boldsymbol e}_y \sin \phi / k$ and Hamiltonian (2.1) can be presented as

(2.3)\begin{align} H = m_e c^2 \sqrt{1 + \left( \frac{p_z}{m_e c} \right)^2 + \left( \frac{p_x}{m_e c} + \frac{\varOmega_{ce}}{kc} \frac{B_w}{B_0} \cos\phi \right)^2 + \left( \frac{p_y}{m_e c} + \frac{x\varOmega_{ce}}{c} - \frac{\varOmega_{ce}}{kc} \frac{B_w}{B_0} \sin\phi \right)^2}. \end{align}

This Hamiltonian does not depend on $y$, and thus canonical momentum $p_y$ is a constant: $\dot {p}_y = -\partial H/\partial y = 0$. Without loss of generality, we can set $p_y = 0$.

In the absence of a wave, this Hamiltonian describes fast $(x,p_x)$ oscillations and slow $(z,p_z)$ oscillations. Thus, we can introduce an adiabatic invariant $I_x=(2{\rm \pi} )^{-1}\oint \,{{{\rm d} x} p_x}$ as an area surrounded by a closed trajectory in the $(x,p_x)$ plane (Landau & Lifshitz Reference Landau and Lifshitz1988):

(2.4)\begin{equation} I_x = \frac{1}{2 {\rm \pi}} \oint\,{{{\rm d} x} p_x} = \frac{m_e c^2}{2 \varOmega_{ce}} \left( \frac{H^2}{m_e^2 c^4} - 1 - \left( \frac{p_z}{m_e c} \right)^2 \right) \end{equation}

with

(2.5a,b)\begin{equation} H = m_e c^2 \sqrt{1 + \left( \frac{p_z}{m_e c} \right)^2 + \frac{2 I_x \varOmega_{ce}}{m_e c^2}},\quad \frac{2 I_x \varOmega_{ce}}{m_e c^2} = \left(\frac{p_x}{m_e c}\right)^2 + \left(\frac{x\varOmega_{ce}}{c}\right)^2. \end{equation}

We consider the canonical transformation $(x, p_x) \to (\psi, I_x)$ given by a generating function $F_2(x, I_x)=(2{\rm \pi} )^{-1}\int \,{{{\rm d} x} p_x}$ (Landau & Lifshitz Reference Landau and Lifshitz1988):

(2.6)\begin{gather} F_2(x, I_x) ={\pm} \frac{ m_e }{2 {\rm \pi}}\int \,{{\rm d} x} \sqrt{\frac{2 I_x \varOmega_{ce}}{m_e} - x^2 \varOmega_{ce}^2} ={\pm} I_x \left[ \sqrt{\frac{m_e \varOmega_{ce}}{2 I_x}} x \sqrt{1 - \frac{m_e \varOmega_{ce}}{2 I_x} x^2} \right.\nonumber\\ + \left. {\rm arcsin}\left( \sqrt{\frac{m_e \varOmega_{ce}}{2 I_x}} x \right) \right]. \end{gather}

The corresponding variable transformations are

(2.7a,b)\begin{equation} p_x ={\pm} m_e \sqrt{\frac{2 I_x \varOmega_{ce}}{m_e} - x^2 \varOmega_{ce}^2},\quad \psi ={\pm} {\rm arcsin}\left( \sqrt{\frac{m_e \varOmega_{ce}}{2 I_x}} x \right). \end{equation}

Thus, equations for $x$ and $p_x$ are

(2.8a,b)\begin{equation} x = \frac{c}{\varOmega_{ce}} \sqrt{\frac{2 I_x \varOmega_{ce}}{m_e c^2}} \sin\psi,\quad p_x = m_e c \sqrt{\frac{2 I_x \varOmega_{ce}}{m_e c^2}} \cos\psi. \end{equation}

A new equation for $\varGamma = H / m_e c^2$ in terms of variables $(z, p_z)$ and $(\psi, I_x)$ can be written as

(2.9)\begin{equation} \varGamma = \sqrt{1 + \left( \frac{p_z}{m_e c} \right)^2 + \frac{2 I_x \varOmega_{ce}}{m_e c^2} + 2 \sqrt{\frac{2 I_x \varOmega_{ce}}{m_e c^2}} \frac{\varOmega_{ce}}{k c} \frac{B_w}{B_0} \cos(\phi + \psi) + \left( \frac{\varOmega_{ce}}{k c} \frac{B_w}{B_0} \right)^2}. \end{equation}

The Hamiltonian $H=m_ec^2\varGamma$ describes the dynamics of two pairs of conjugated variables, $(z, p_z)$ and $(\psi, I_x)$. The system of Hamiltonian equations can be solved with respect to these variables and time $t$. The general approach consists of expanding $H$ over $B_w/B_0$, keeping only the linear $\sim B_w/B_0$ term:

(2.10a,b)\begin{align} H = m_ec^2\gamma + m_e c^2\sqrt{\frac{2 I_x \varOmega_{ce}}{m_e c^2}} \frac{\varOmega_{ce}}{\gamma k c} \frac{B_w}{B_0} \cos(\phi + \psi),\quad \gamma=\sqrt{1 + \left( \frac{p_z}{m_e c} \right)^2 + \frac{2 I_x \varOmega_{ce}}{m_e c^2}}. \end{align}

We introduce the wave modulation through the $B_w$ dependence on the wave phase $\phi$ ($B_w(\phi )$ periodicity mimics effect if multiple wave packets):

(2.11)\begin{equation} B_w = B_{m} \frac{1 - e^{{-}h \sin^2(\phi / 2l)}}{1 - e^{{-}h}}, \end{equation}

where $B_m$ is the peak wave amplitude, $l$ defines the number of wave oscillations (periods) within one wave packet and $h$ controls the intensity of modulations. The effective wave amplitude $B_{{\rm eff}}$ can be determined by averaging $B_w$ over the period of modulations:

(2.12)\begin{equation} B_{{\rm eff}} = \sqrt{\left\langle B_w^2 \right\rangle}_{\phi\in[0,2{\rm \pi} l]}= B_{m} \frac{\sqrt{1 - 2 {\rm I}_0\left( h/2 \right) e^{ {-}h/2 } + {\rm I}_0\left( h \right) e^{ {-}h }}}{1 - e^{{-}h}}, \end{equation}

where ${\rm I}_n(z)$ is the modified Bessel function of the first kind. The average intensity of the plane wave with $B_w = B_{{\rm eff}} = {\rm const.}$ equals the intensity of the modulated wave with $B_w$ given by (2.11). Thus, to describe the electron diffusion theoretically, the approximation $B_w \approx B_{{\rm eff}} = {\rm const.}$ can be used. However, this assumption neglects nonlinear effects and, therefore, to determine whether the theory is applicable, we perform numerical simulations of modulated waves. Parameter $h$ in (2.11) determines the depth of the modulation: at $h \to \infty$ and $l\to \infty$ the wave packet reduces to a plane wave. Parameter $l$ determines the wave packet size with typical values $l\in [10,30]$ (see Zhang et al. Reference Zhang, Mourenas, Artemyev, Angelopoulos, Bortnik, Thorne, Kurth, Kletzing and Hospodarsky2019, Reference Zhang, Demekhov, Katoh, Nunn, Tao, Mourenas, Omura, Artemyev and Angelopoulos2021). For numerical simulations we use $h=1$ and $l=20$. Note that these parameters well satisfy the condition $|{\rm d}B_w/{\rm d}z| \ll |(\partial \phi /\partial z)B_w|$, because ${\rm d}B_w/{\rm d}z=(\partial \phi /\partial z) ({\rm d}B_w/{\rm d}\phi )$ and $({\rm d}B_w/{\rm d}\phi )B_w^{-1}\sim h/l \ll 1$.

We consider particles having the same initial energy $E_0$ and pitch angle $\alpha _0$, but random wave phase $\phi$ and gyrophase $\psi$. Figure 2 shows typical examples of electron resonance interactions with whistler-mode waves obtained by a numerical integration of Hamiltonian equations: diffusive electron scattering by low-amplitude wave (figure 2a), nonlinear resonant interactions with intense coherent wave (figure 2b) and nonlinear resonant interactions with well-modulated intense wave (figure 2c). The resonant interaction occurs once per simulation interval (half of the bounce period), i.e. there is only one point along the unperturbed particle trajectory where the resonant condition $\dot {\phi } + \dot {\psi } = 0$ is satisfied. The diffusive scattering is characterized by a symmetric (relative to zero) distribution of energy changes, and thus this process should be described by a diffusion rate $\sim \langle (\Delta E)^2\rangle$. The nonlinear resonances with an intense coherent wave are characterized by a small population changing the energy significantly ($\Delta E>0$, the phase trapping effect) and a large population with a small energy change, but almost identical for all particles ($\Delta E<0$, the phase bunching effect). Thus, nonlinear resonances with a coherent wave should be described separately for trapped and bunched particle populations, and due to the large energy change of the trapped population it is not thought to be possible to include this into the Fokker–Planck equation (see Hsieh & Omura Reference Hsieh and Omura2017a,Reference Hsieh and Omurab; Artemyev et al. Reference Artemyev, Neishtadt, Vasiliev, Zhang, Mourenas and Vainchtein2021b; Zhang et al. Reference Zhang, Artemyev, Angelopoulos, Tsai, Wilkins, Kasahara, Mourenas, Yokota, Keika and Hori2022). The nonlinear resonances with modulated waves are characterized by: (i) an increase of probability for the $\Delta E>0$ changes; (ii) but also with a decrease in size of $|\Delta E|$ itself; and (iii) further a randomization of energy change for this population. Thus, for modulated waves, the resonant wave–particle interaction is closer to diffusive scattering (Tao et al. Reference Tao, Bortnik, Albert, Thorne and Li2013; Zhang et al. Reference Zhang, Agapitov, Artemyev, Mourenas, Angelopoulos, Kurth, Bonnell and Hospodarsky2020a; An et al. Reference An, Wu and Tao2022; Mourenas et al. Reference Mourenas, Zhang, Nunn, Artemyev, Angelopoulos, Tsai and Wilkins2022).

Figure 2. A set of trajectories obtained by the numerical integration of Hamiltonian equations for a system (2.10a,b). Results are shown for (a) low-amplitude wave with diffusive scattering, (b) coherent high-amplitude wave with trapping and bunching and (c) modulated high-amplitude wave with diffusive-like scattering. System parameters are: electron energy $E_0=100$ keV, equatorial pitch angle $\alpha _0 = 40^\circ$, number of particles $N = 100$, wave amplitude $B_w = 5$ pT (a), $B_w = 500$ pT (b,c). System parameters correspond to $L\text {-shell} = 6$, whistler-mode waves with frequency equal to $0.35$ of the electron cyclotron frequency at the equator, and the constant plasma frequency equal to $10$ times the electron cyclotron frequency at the equator.

In Figure 3 we show how wave intensity and wave modulation control the efficiency of the nonlinear interactions. Figure 3(a–c) shows distributions of energy changes $\Delta E$ depending on the normalized wave amplitude $\varepsilon =B_w/B_{0}(0)$. Below $\varepsilon \sim 10^{-4}$–$10^{-3}$, the resonant interaction is diffusive with a symmetric $\Delta E$ distribution and $\langle {(\Delta E)^2}\rangle ^{1/2}$ linearly growing with $\varepsilon$. This dependence $\langle {(\Delta E)^2}\rangle ^{1/2}\propto \varepsilon$ demonstrates the applicability of the unperturbed trajectory approximation, a core assumption of the quasi-linear diffusion theory (Tao et al. Reference Tao, Bortnik, Albert, Liu and Thorne2011, Reference Tao, Bortnik, Albert and Thorne2012a; Allanson et al. Reference Allanson, Watt, Ratcliffe, Allison, Meredith, Bentley, Ross and Glauert2020). After the wave amplitude reaches a certain threshold (depending on the electron energy, pitch angle and system characteristics; see Omura et al. Reference Omura, Matsumoto, Nunn and Rycroft1991; Shklyar & Matsumoto Reference Shklyar and Matsumoto2009), the resonant interaction becomes nonlinear with a clear formation of a population of trapped electrons (a separate group of large positive $\Delta E$ in figure 3a–c) and a highly asymmetric $\Delta E$ distribution (most of the electrons experience phase bunching and form a large maximum at $\Delta E<0$). For such nonlinear resonant interactions the $\Delta E$ distribution with phase-trapped and phase-bunched populations cannot be characterized by a diffusion $\langle {(\Delta E)^2}\rangle$ only, and thus it is not known how to include this regime of resonant interactions into the Fokker–Planck equation.

Figure 3. Distributions of the energy change $\Delta E$ for different $B_{{\rm eff}}$ with $h=1$ and $l=20$ for coherent (a–c) and modulated (d–f) waves with $\Delta E \propto B_w$ fitting (black dashed lines) and with $\langle {(\Delta E)^2}\rangle ^{1/2}$ profiles (black dots): (a,d) $E_0 = 100$ keV, $\alpha _0 = 40^{\circ }$; (b,e) $E_0 = 100$ keV, $\alpha _0 = 60^{\circ }$; (c, f) $E_0 = 300$ keV, $\alpha _0 = 60^{\circ }$. For each $B_w$ we use $10^4$ trajectories to evaluate the $\Delta E$ distribution, each particle resonates with the wave only once and $\Delta E$ is the energy change for a single resonance, initial particle phases and wave phases are random and thus for a modulated case the actual wave amplitude is different for different particles having the same energy/pitch angle. For each $B_{{\rm eff}}$ the $\Delta E$ distribution is normalized to one. System parameters correspond to $L\text {-shell} = 6$, whistler-mode waves with frequency equal to $0.35$ of the electron cyclotron frequency at the equator, and the constant plasma frequency equal to $10$ times the electron cyclotron frequency at the equator.

Figure 3(d–f) shows $\Delta E$ distributions as a function of wave amplitude for strongly modulated waves. For small wave intensity, $\varepsilon <10^{-4}$, there is the same diffusive regime of wave–particle resonant interactions as for non-modulated waves: a symmetric $\Delta E$ distribution with $\langle {(\Delta E)^2}\rangle ^{1/2}\propto \varepsilon$. With amplitude increasing, the regime of wave–particle resonant interaction changes. However, the wave modulation does not allow strong trapping acceleration (there is no population with large positive $\Delta E$), but increases the number of trapped particles (the populations of particles with $\Delta E>0$ and with $\Delta E<0$ are comparable even for $\varepsilon >10^{-3}$). Thus, for intense modulated waves, the $\Delta E$ distribution remains almost symmetric and can be characterized by $\langle {(\Delta E)^2}\rangle$. Despite such a symmetric $\Delta E$ distribution, the resonant interaction for large $\varepsilon$ is nonlinear, and the wave field alters electron dynamics in the resonance. This results in inapplicability of the unperturbed trajectory approximation, and thus $\langle {(\Delta E)^2}\rangle ^{1/2}\propto \varepsilon ^A$ with $A<1$. Therefore, for intense modulated waves we deal with diffusion of resonant electrons, but it is not a quasi-linear diffusion. In this study, we aim to derive the diffusion coefficient $\sim \langle {(\Delta E)^2}\rangle$ as a function of $\varepsilon$ for a wide $\varepsilon$ range. Figure 3(d–f) shows that $\langle {(\Delta E)^2}\rangle$ for small $\varepsilon$ should be similar to that for the quasi-linear diffusion (see also Albert Reference Albert2010), whereas for large $\varepsilon$ the dispersion $\langle {(\Delta E)^2}\rangle$ should be about $\langle {\Delta E}\rangle ^2$ of the phase-bunched particle population. We have checked this assumption with an analytical approach, which we further introduce.

3. The diffusion coefficient model for an arbitrary wave intensity

The Hamiltonian $H=m_ec^2\gamma$ with $\gamma$ given by (2.9) determines the electron dynamics. For this Hamiltonian system, we derive an analytical approximation for the energy change $\Delta E = m_e c^2 \Delta \gamma$ due to a single resonant interaction. Such a change depends on the initial electron energy and pitch angle, and on the initial phase $\phi +\psi$. Thus, we finally aim to find a mean value $V_\gamma =\langle \Delta \gamma \rangle$ and variance $D_{\gamma \gamma }=\langle (\Delta \gamma )^2 \rangle$ with the averaging over initial phases.

To derive the $\Delta \gamma$ equation for $H=m_ec^2\gamma$, we follow the approach from Neishtadt & Vasiliev (Reference Neishtadt and Vasiliev2006) and Artemyev et al. (Reference Artemyev, Neishtadt, Vainchtein, Vasiliev, Vasko and Zelenyi2018a) for the perturbation theory application to a resonant system containing a phase and a small wave amplitude, $B_w / B_0 \ll 1$ (more precisely, $\varepsilon = B_w / B_0(0) \ll 1$). For analytical consideration we assume that $B_w$ is a constant or a function of the magnetic latitude with spatial gradient much weaker than wave phase gradients $\partial \phi /\partial s=k$.

The Hamiltonian from (2.10a,b) is time-dependent as $\phi = \phi (t)$. To find the invariant equation, we define a generation function of the second kind $F_2(s, \psi, P_z, I)$ for $(s, p_z, \psi, I_x) \to (\tilde {z}, P_z, \zeta, I)$:

(3.1)\begin{equation} F_2(s, \psi, P_z, I) = (\phi + \psi)I + P_z z \quad \begin{cases} p_z = \dfrac{\partial F_2}{\partial z} = P_z + k I ,\\ I_x = \dfrac{\partial F_2}{\partial \psi} = I ,\\ \tilde{z} = \dfrac{\partial F_2}{\partial P_z} = z ,\\ \zeta = \dfrac{\partial F_2}{\partial I} = \phi + \psi ,\\ {\mathcal{H}} - H = \dfrac{\partial F_2}{\partial t} ={-}\omega I. \end{cases} \end{equation}

As a result, a modified Hamiltonian ${\mathcal {H}}$ is

(3.2)\begin{equation} \left. \begin{aligned} {\mathcal{H}} & = m_e c^2 \gamma - \omega I + m_e c^2 \sqrt{\frac{2 I \varOmega_{ce}}{m_e c^2}} \frac{\varOmega_{ce}}{\gamma k c} \frac{B_w}{B_0} \cos\zeta, \\ \gamma & = \sqrt{1 + \left( \frac{P_z + k I}{m_e c} \right)^2 + \frac{2 I \varOmega_{ce}}{m_e c^2}}. \end{aligned} \right\} \end{equation}

Note that $\partial _t {\mathcal {H}} = 0$. Thus, in the zeroth order of $\varepsilon$ the first invariant of this system can be written as

(3.3)\begin{equation} m_e c^2 \gamma - \omega I = {\rm const.} \end{equation}

Without a wave perturbation, $I$ is invariant ($I = \textrm {const.}$) by its definition, and $\gamma$ has a constant value, considering (3.3). Thus, all derivatives of those variables are first-order terms of $\varepsilon$ or higher. Resonant wave–particle interactions change the particle energy $E = m_e c^2 \gamma$: $\gamma _0 \to \gamma _0 + \Delta \gamma$, where $\gamma _0$ is the initial Lorentz factor (at $t = 0$). The resonance equation can be written as $\dot {\zeta } = \partial {\mathcal {H}}/\partial I = 0$. This equation characterizes the dynamics of the system near the resonance point and, considering (3.2), can be written as (subscript $R$ means that the function is evaluated at resonance)

(3.4a,b)\begin{equation} m_e c^2 \left. \frac{\partial \gamma}{\partial I} \right\vert_R - \omega = 0, \quad \gamma_R= \frac{k_R c}{\omega} \frac{P_z + k_R I_R}{m_e c} + \frac{\varOmega_{ce,R}}{\omega}. \end{equation}

The Hamiltonian ${\mathcal {H}}$ has to be expanded near the resonance point $\dot {\zeta } = 0$ with values $\gamma _R = \gamma _0 + O(\varepsilon )$, $I_R = I_0 + O(\varepsilon )$:

(3.5)\begin{equation} \left. \begin{aligned} {\mathcal{H}} & \approx {\mathcal{H}}_R + \frac{m_e c^2}{2} \left.\frac{\partial^2 \gamma}{\partial I^2}\right\vert_{R} (I - I_R)^2 + m_e c^2 \sqrt{\frac{2 I_R \varOmega_{ce,R}}{m_e c^2}} \frac{\varOmega_{ce,R}}{\gamma_R k_R c} \frac{B_w}{B_{0,R}} \cos\zeta, \\ {\mathcal{H}}_R & = m_e c^2 \gamma_R - \omega I_R, \end{aligned} \right\} \end{equation}

where

(3.6)\begin{equation} \left.\frac{\partial^2 \gamma}{\partial I^2}\right\vert_{R} = \frac{k_R^2 c^2 - \omega^2}{m_e^2 c^4 \gamma_R}. \end{equation}

Thus, the final form of ${\mathcal {H}}$ is

(3.7)\begin{equation} \left. \begin{aligned} {\mathcal{H}} & = {\mathcal{H}}_R + \frac{(I - I_R)^2}{2 g} + m_e c^2\sqrt{\frac{2 I_R \varOmega_{ce,R}}{m_e c^2}} \frac{\varOmega_{ce,R}}{\gamma_R k_R} \frac{B_w}{B_{0,R}} \cos\zeta, \\ g^{{-}1} & \equiv m_e c^2 \left.\frac{\partial^2 \gamma}{\partial I^2}\right\vert_{R} =\frac{k_R^2 c^2 - \omega^2}{m_e c^2 \gamma_R}. \end{aligned} \right\} \end{equation}

The second term of ${\mathcal {H}}$ from (3.7) is an analogue of the kinetic energy of the particle motion near the resonance, with $I - I_R$ being the canonical momentum and $g$ playing the role of mass. Thus, we introduce a generating function of the third kind $F_3(\tilde {z}, \tilde {\zeta }, P_z, I)$ for $(z, P_z, \zeta, I) \to (\tilde {z}, \tilde {P}_z, \tilde {\zeta }, P_{\zeta })$:

(3.8)\begin{equation} F_3(\tilde{z}, \tilde{\zeta}, P_z, I) ={-}(I - I_R) \tilde{\zeta} - P_z \tilde{z} \quad \begin{cases} z ={-}\dfrac{\partial F_3}{\partial P_z} = \tilde{z} - \dfrac{\partial I_R}{\partial P_z} \tilde{\zeta},\\ \zeta ={-}\dfrac{\partial F_3}{\partial I} = \tilde{\zeta} ,\\ \tilde{P}_z ={-}\dfrac{\partial F_3}{\partial \tilde{z}} = P_z - \dfrac{\partial I_R}{\partial \tilde{z}} \tilde{\zeta} ,\\ P_{\zeta} ={-}\dfrac{\partial F_3}{\partial \tilde{\zeta}} = I - I_R . \end{cases} \end{equation}

As $I$ is an invariant in the unperturbed system, $\partial _{P_z} I_R \sim \varepsilon$ and $\partial _{\tilde {z}} I_R \sim \varepsilon$. This means that the first term of the Hamiltonian ${\mathcal {H}}$ can be expanded with respect to $\varepsilon$:

(3.9)\begin{align} {\mathcal{H}}_R & = {\mathcal{H}}_R\left(z, P_z\right) = {\mathcal{H}}_R\left(\tilde{z} - \frac{\partial I_R}{\partial P_z} \zeta, \tilde{P}_z + \frac{\partial I_R}{\partial \tilde{z}} \zeta\right) \nonumber\\ & \approx {\mathcal{H}}_R\left(\tilde{z}, \tilde{P}_z\right) - \frac{\partial {\mathcal{H}}_R}{\partial \tilde{z}} \frac{\partial I_R}{\partial \tilde{P}_z} \zeta + \frac{\partial {\mathcal{H}}_R}{\partial \tilde{P}_z} \frac{\partial I_R}{\partial \tilde{z}} \zeta \nonumber\\ & = {\mathcal{H}}_R\left(\tilde{z}, \tilde{P}_z\right) + \left\{{\mathcal{H}}_R, I_R \right\}_{\tilde{z}, \tilde{P}_z} \zeta . \end{align}

Substituting the expanded form of ${\mathcal {H}}_R$ into (3.7), we obtain two separate parts of the Hamiltonian ${\mathcal {H}}={\mathcal {H}}_R+{\mathcal {H}}_\zeta$, ${\mathcal {H}}_R = m_e c^2 \gamma _R - \omega I_R$ with canonical variables $(\tilde {z}, \tilde {P}_z, \zeta, P_{\zeta })$:

(3.10)\begin{equation} {\mathcal{H}}_{\zeta} = \frac{P_{\zeta}^2}{2g} + \left\{{\mathcal{H}}_R, I_R \right\}_{\tilde{z}, \tilde{P}_z} \zeta + m_e c^2 \sqrt{\frac{2 I_R \varOmega_{ce,R}}{m_e c^2}} \frac{\varOmega_{ce,R}}{\gamma_R k_R c} \frac{B_w}{B_{0,R}} \cos\zeta. \end{equation}

Hamiltonian ${\mathcal {H}}_R$ describes $(\tilde {z}, \tilde {P}_z)\approx (z, P_z)$ dynamics in the resonance, and this Hamiltonian does not depend on the fast phase $\zeta$. Hamiltonian ${\mathcal {H}}_{\zeta }$ is a $\zeta$-dependent pendulum Hamiltonian that describes fast phase and conjugated momentum dynamics around the resonance. Coefficients of Hamiltonian ${\mathcal {H}}_\zeta$ depend on $(\tilde {z}, \tilde {P}_z)$ and slowly change along the resonant trajectory.

In the Hamiltonian ${\mathcal {H}}_\zeta$, an effective potential energy, ${\mathcal {H}}_\zeta -P_\zeta ^2/2g$, contains two terms: the first term $\sim \{{\mathcal {H}}_R,I_R\}_{\tilde {z}, \tilde {P}_z}$ describes the impact of the background magnetic field gradient and the second term $\sim B_w$ describes the effect of the wave's field. The important system parameter is the ratio of magnitudes of these two terms:

(3.11)\begin{equation} a = m_e c^2 \sqrt{\frac{2 I_R \varOmega_{ce,R}}{m_e c^2}} \frac{\varOmega_{ce,R}}{\gamma_R k_R c} \frac{B_w}{B_{0,R}} \left\{{\mathcal{H}}_R, I_R \right\}_{\tilde{z}, \tilde{P}_z}^{{-}1}. \end{equation}

Taking into account ${\mathcal {H}}_R=m_ec^2\gamma _R-\omega I_R$, we can rewrite the Poisson bracket:

(3.12)\begin{equation} \left\{{\mathcal{H}}_R, I_R \right\}_{\tilde{z}, \tilde{P}_z} = m_e c^2 \left\{ \gamma_R, I_R \right\}_{\tilde{z},\tilde{P}_z} \approx m_e c^2 \left\{ \gamma_R, I_R \right\}_{z, P_z}, \end{equation}

where $\gamma _R$ and $I_R$ are given by (3.2) and (3.4a,b). Combining these equations, we obtain (subscript $R$ is omitted in equations below)

(3.13a,b)\begin{equation} \gamma = \frac{k c}{\omega} \frac{P_z + k I}{m_e c} + \frac{\varOmega_{ce}}{\omega},\quad I = \frac{m_e c^2}{2 \varOmega_{ce}} \left[ \gamma^2 - \left(\frac{\omega}{k c} \right)^2 \left( \gamma - \frac{\varOmega_{ce}}{\omega} \right)^2 - 1 \right]. \end{equation}

Therefore, the Poisson bracket $\{\cdots \} = \{\cdots \}_{s,P_{\parallel }}$ can be rewritten as

(3.14)\begin{align} \left\{ \gamma, I \right\} & = I \varOmega_{ce} \left\{ \gamma, \frac{1}{\varOmega_{ce}} \right\} - \frac{m_e c^2}{2 \varOmega_{ce}} \left\{ \gamma, \left(\frac{\omega}{k c} \right)^2 \left( \gamma - \frac{\varOmega_{ce}}{\omega} \right)^2 \right\} ={-}I \frac{\partial \gamma}{\partial P_z} \partial_z \left[ \ln \varOmega_{ce} \right] \nonumber\\ & \quad - \frac{m_e c^2}{\varOmega_{ce}} \frac{\omega}{k c} \left( \gamma - \frac{\varOmega_{ce}}{\omega} \right) \left\{ \gamma, \frac{\omega}{k c} \left( \gamma - \frac{\varOmega_{ce}}{\omega} \right) \right\} ={-}I \frac{\partial \gamma}{\partial P_z} \partial_z \left[ \ln \varOmega_{ce} \right] \nonumber\\ & \quad - \frac{m_e c^2}{\varOmega_{ce}} \frac{\omega}{k c} \left( \gamma - \frac{\varOmega_{ce}}{\omega} \right) \left[ \left( \gamma - \frac{\varOmega_{ce}}{\omega} \right) \frac{\partial \gamma}{\partial P_z} \partial_z \left[ \frac{\omega}{k c} \right] - \frac{\omega}{k c} \frac{\varOmega_{ce}}{\omega} \frac{\partial \gamma}{\partial P_z} \partial_z \left[ \ln \varOmega_{ce} \right] \right] \nonumber\\ & = \frac{m_e c^2}{\omega} \frac{\partial \gamma}{\partial P_z} \left[ \left(\left(\frac{\omega}{kc}\right)^2 \left( \gamma - \frac{\varOmega_{ce}}{\omega} \right) - \frac{\omega I}{m_e c^2} \right) \partial_z \left[ \ln \varOmega_{ce} \right] - \frac{\omega}{\varOmega_{ce}} \left(\frac{\omega}{kc}\right)^2 \right. \nonumber\\ & \quad \times \left. \left( \gamma - \frac{\varOmega_{ce}}{\omega} \right)^2 \partial_z \left[ \ln \frac{\omega}{k c} \right] \right]. \end{align}

To determine $\partial _z [ \ln ({\omega }/{k c}) ]$ we use the cold plasma approximation, $k c / \omega = (\varOmega _{pe} / \omega ) ( \varOmega _{ce}/\omega - 1)^{-{1}/{2}}$, with constant plasma frequency $\varOmega _{pe}$. Therefore, for partial derivatives in $\{ \gamma, I \}$ we can write

(3.15)\begin{equation} \left. \begin{aligned} \frac{\partial \gamma}{\partial P_z} & = \frac{1}{m_e c} \frac{k c}{\omega} \left( 1 + k \frac{\partial I}{\partial P_z} \right),\quad \frac{\partial I}{\partial P_z} = \frac{m_e c^2}{\varOmega_{ce}} \frac{\partial \gamma}{\partial P_z} \left[ \gamma - \left(\frac{\omega}{k c} \right)^2 \left( \gamma - \frac{\varOmega_{ce}}{\omega} \right) \right],\\ \frac{\partial \gamma}{\partial P_z} & ={-}\frac{1}{\gamma m_e c} \frac{\varOmega_{ce}}{\omega} \frac{k c}{\omega} \left[ \left(\frac{k c}{\omega}\right)^2 - 1 \right]^{{-}1},\quad \frac{\partial }{\partial s} \left[ \ln \frac{\omega}{k c} \right] = \frac{1}{2} \frac{\partial_z \left[ \ln \varOmega_{ce} \right]}{1 - \omega/\varOmega_{ce}} \end{aligned} \right\} \end{equation}

and

(3.16)\begin{gather} \left\{ \gamma, I \right\} = \frac{\partial_z \left[ \ln \varOmega_{ce} \right]}{\gamma k} \left[\frac{\left( \gamma - \varOmega_{ce}/\omega \right)^2}{2 \left( 1 - \omega/\varOmega_{ce} \right)} - \varOmega_{ce}/\omega \left(\gamma - \frac{\varOmega_{ce}}{\omega} - \frac{\omega I}{m_e c^2} \left(\frac{k c}{\omega}\right)^2 \right) \right] \nonumber\\ \times \left[ \left(\frac{k c}{\omega}\right)^2 - 1 \right]^{{-}1}. \end{gather}

Having $\{ \gamma, I \}$, we can determine $a$ from (3.11) at the resonance. As we show, this is the main parameter controlling the energy change due to the resonant wave–particle interaction.

To write an equation for the resonant energy change, $\Delta \gamma$, we use the invariant from (3.3):

(3.17)\begin{align} \Delta \gamma & = \frac{\omega}{m_e c^2} \Delta I = \frac{2 \omega}{m_e c^2} \int_{-\infty}^{t_R} \,{\rm d}t \ \dot{I} ={-}\frac{2 \omega}{m_e c^2} \int_{-\infty}^{t_R} \,{\rm d}t \frac{\partial {\mathcal{H}}}{\partial \zeta} \nonumber\\ & \approx 2 \omega \sqrt{\frac{2 I_R \varOmega_{ce,R}}{m_e c^2}} \frac{\varOmega_{ce,R}}{\gamma_R k_R c} \frac{B_w}{B_{0,R}} \int_{-\infty}^{\zeta_R} \,{\rm d}\zeta \frac{\sin\zeta}{\dot{\zeta}}, \end{align}

where $\varOmega _{ce,R}=eB_{0,R}/m_ec$ and $B_{0,R}$ is defined at the resonant $z=z(t_R)$ for given initial energy and pitch angle. For $\dot \zeta$ we use the Hamiltonian equation $\dot \zeta =\partial {\mathcal {H}}_\zeta /\partial P_\zeta$. The resonance is defined by $\dot {\zeta } = P_{\zeta } / g = 0$ with the solution $\zeta _R$ according to $\zeta _R + a\cos \zeta _R=\xi$ with the resonant energy $\xi = {\mathcal {H}}_\zeta / \{ {\mathcal {H}}_R, I_R \}$. Thus, (3.17) can be written as

(3.18)\begin{equation} \left. \begin{aligned} \Delta \gamma & = 2 \omega \sqrt{\frac{g I_R \varOmega_{ce,R}}{m_e c^2}} \frac{\varOmega_{ce,R}}{\gamma_R k_R c} \frac{B_w}{B_{0,R}} \left\{ {\mathcal{H}}_R, I_R \right\}_{\tilde{z}, \tilde{P}_z}^{{-}1/2} \int_{-\infty}^{\zeta_R} \,{\rm d}\zeta \frac{\sin\zeta}{\sqrt{\xi - \zeta - a \cos\zeta}} \\ & = \left( \frac{2 I_R \varOmega_{ce,R}}{m_e c^2} \right)^{\frac{1}{4}} \sqrt{\frac{2 \varOmega_{ce,R} / k_R c}{k_R^2 c^2 / \omega^2 - 1} \frac{B_w}{B_{0,R}}} f(a,\xi), \\ f(a,\xi) & = \int_{-\infty}^{\zeta_R} \,{\rm d}\zeta \frac{\sqrt{a} \sin\zeta}{\sqrt{\xi - \zeta - a \cos\zeta}} . \end{aligned} \right\} \end{equation}

Figure 4(a) shows the $f(a,\xi )$ function. This function is periodic with period $2{\rm \pi}$ for $\xi$ (this can be shown analytically; see Appendix A).

Figure 4. (a) Profiles of $f(a,\xi )$ for three $a$ values. (b) Profiles of $\langle f(a,\xi )\rangle _{\xi }$ and $\langle (f(a,\xi ))^2\rangle _{\xi }$.

The value of $\xi$ varies with the initial conditions, i.e. with wave phase $\phi$, gyrophase $\psi$ and electron position on the trajectory in $(s,P_\parallel )$ plane far from the resonance. Assuming a uniform distribution of these parameters, we numerically integrate an ensemble of trajectories as determined by (2.10a,b), and so determine the probability function for $\xi$. Figure 5 shows such distributions of $\xi$ for three sets of the system parameters. These distributions are very close to uniform distributions with $\xi \in [0,2{\rm \pi} ]$, which suggests that averaging over initial parameters $\phi$, $\psi$ and $s$ can be substituted with averaging over $\xi$ with constant weights (see also discussion in Itin, Neishtadt & Vasiliev (Reference Itin, Neishtadt and Vasiliev2000) and Albert et al. (Reference Albert, Artemyev, Li, Gan and Ma2022)):

(3.19)\begin{equation} \iiint_\varPi \,{\rm d}\phi \, {\rm d}\psi \, {\rm d}\lambda \to \int_{\xi_0}^{\xi_0 + 2 {\rm \pi}} \,{\rm d}\xi, \quad \begin{cases} \xi_0 - \zeta_0 - a \cos\zeta_0 = 0 ,\\ 1 = a \sin \zeta_0, \end{cases} \end{equation}

where $\varPi$ is the parametric range ($3D$ uniform distribution) and $\xi _0 = \xi _0(a)$ determines the case when the integral diverges near the resonance point (see Appendix A). Thus, the equation for the variance $D_{\gamma \gamma }$ can be written as

(3.20)\begin{equation} D_{\gamma \gamma}= \langle (\Delta \gamma)^2 \rangle_{\xi} = \frac{2 \varOmega_{ce,R} / k_R c}{k_R^2 c^2 / \omega^2 - 1} \frac{B_w}{B_{0,R}} \sqrt{ \frac{2 I_R \varOmega_{ce,R}}{m_e c^2} } \langle f^2(a,\xi) \rangle_{\xi}, \end{equation}

where $D_{\gamma \gamma }$ can be considered as a diffusion rate for a unit time interval between two resonant interactions (the actual diffusion rate is the ratio of $\langle D_{\gamma \gamma } \rangle$ and a fraction of electron bounce period).

Figure 5. Probability distributions of $\xi$ for $10^5$ trajectories and (a) $E_0 = 100$ keV, $\alpha _0 = 40^{\circ }$, $B_w = 50$ pT, (b) $E_0 = 100$ keV, $\alpha _0 = 40^{\circ }$, $B_w = 500$ pT, (c) $E_0 = 300$ keV, $\alpha _0 = 60^{\circ }$, $B_w = 500$ pT.

The function $f(a,\xi )$ determines the difference of the diffusion coefficient $\sim D_{\gamma \gamma }$ and the quasi-linear model. For the case of $a \ll 1$, $D_{\gamma \gamma }$ asymptotically tends to the quasi-linear equation $D_{\gamma \gamma } = \langle (\Delta \gamma )^2 \rangle \sim (B_w / B_0)^2$ (Albert Reference Albert2010). To verify that, we have to expand $f(a,\xi )$ in a Taylor series:

(3.21)\begin{equation} \left. \begin{aligned} \lim_{a\to0} \frac{f(a,\xi)}{\sqrt{a}} & = \int_{-\infty}^{\xi} \,{\rm d}\zeta \frac{\sin\zeta}{\sqrt{\xi - \zeta}} = \sqrt{\rm \pi} \sin(\xi - {\rm \pi}/ 4),\\ \lim_{a\to0} \langle f(a,\xi) \rangle_{\xi} & = 0, \quad \lim_{a\to0} \frac{\langle f^2(a,\xi) \rangle_{\xi}}{a} = \frac{\rm \pi}{2}, \end{aligned} \right\} \end{equation}

and thus $D_{\gamma \gamma } \propto (B_w/B_0)\langle f^2(a,\xi ) \rangle _{\xi } \propto (B_w/B_0)^2$, because $a\propto B_w/B_0$ (see (3.11)).

Figure 4(b) shows that there are two $a$ ranges with different wave–particle resonant effects: for $a<1$ the mean value of $\langle f(a,\xi ) \rangle _{\xi }$ equals zero and there is only a particle diffusion $\propto a \langle f^2(a,\xi ) \rangle _{\xi }$, whereas for $a>1$ there is a non-zero negative $\langle f(a,\xi ) \rangle _{\xi }$. The diffusion $\propto a \langle f^2(a,\xi ) \rangle _{\xi }$ has a local maximum at $a\approx 1.0392$, a local minimum at $a\approx 1.5923$ and then increases with $a$ as $\propto a$. The mean value $\langle f(a,\xi ) \rangle _{\xi }$ has an asymptote $4 \sqrt {2} / {\rm \pi}$ for $a\gg 1$. There is an important property of $f(a,\xi )$: $\langle f^2(a,\xi ) \rangle _{\xi }-\langle f(a,\xi ) \rangle _{\xi }^2 \to 0$ for $a\to \infty$ (see Appendix A for details). This property defines the behaviour of the diffusion rate $\propto a\langle f^2(a,\xi ) \rangle _{\xi }$ for $a\gg 1$.

As shown in figure 3, the strong wave modulation should result in a symmetric distribution of $\Delta \gamma$ with a zero mean value and with the dispersion $\langle (\Delta \gamma )^2 \rangle$ about $\langle (\Delta \gamma )^2 \rangle _{\xi }$ where $\xi$ averaging is performed for population of phase-bunched particles (i.e. particles with $\Delta \gamma <0$). Therefore, for such strongly modulated waves, we consider $\langle (\Delta \gamma )^2 \rangle _{\xi }$ as a diffusion rate for both $a<1$ and $a>1$ parametric ranges. Figure 6 shows $\langle (\Delta \gamma )^2 \rangle _{\xi }$ distributions in energy and pitch-angle space for three typical wave intensities. Electrons of lower energy/higher pitch angle resonate with waves closer to the equatorial plane, where $a$ is large because $\{ \gamma _R, I_R \}_{s, P_{\parallel }} \propto \partial \varOmega _{ce}/\partial s$ and tends to zero around the equator.

Figure 6. Two-dimensional energy/pitch-angle maps of $(E_0,\alpha _{0})$ for (a) $B_w = 100$ pT, (b) $B_w = 500$ pT and (c) $B_w = 1000$ pT. Black curve shows $a=1$.

5. Summary

We have derived the diffusion rate for relativistic electron scattering by intense whistler-mode waves of arbitrary amplitude. This diffusion rate repeats the $D\propto B_w^2$ scaling of the quasi-linear diffusion rate $D_{QL}$ for small amplitudes, and tends to $D\propto B_w$ scaling for amplitudes exceeding the threshold of nonlinear wave–particle interactions. Therefore, under the assumption of the absence of main nonlinear resonant effects (phase trapping and phase bunching) due to low wave coherence, the quasi-linear diffusion model will overestimate the diffusion rate for large amplitudes, because $D/D_{QL}\propto B_w^{-1}$. This result demonstrates that the extrapolation of quasi-linear diffusion models should not work for high wave amplitudes, whereas the approximation of the total destruction of nonlinear resonant effects (phase trapping and phase bunching) due to low wave coherence/strong wave modulation will underestimate the rates of electron flux dynamics. Thus, for accurate inclusion of statistics of high wave amplitudes into electron flux models, we must account for the contribution of nonlinear resonant effects.