Diffusive scattering of energetic electrons by intense whistler-mode waves in an inhomogeneous plasma

Abstract

Electron resonant interactions with electromagnetic whistler-mode waves play an important role in electron flux dynamics in various space plasma systems: planetary radiation belts, bow shocks, solar wind and magnetic reconnection regions. Two key wave characteristics determining the regime of wave–particle interactions are the wave intensity and the wave coherency. The classical quasi-linear diffusion approach describes well electron diffusion by incoherent and low-amplitude waves, whereas the nonlinear resonant models describe electron phase bunching and trapping by highly coherent intense waves. This study is devoted to the investigation of the regime of electron resonant interactions with incoherent but intense waves. Although this regime is characterized by electron diffusion, we show that diffusion rates scale linearly with the wave amplitude, $D\propto B_w$, in contrast to the quasi-linear diffusion scaling $D_{QL}\propto B_w^2$. Using observed wave amplitude distributions, we demonstrate that the quasi-linear diffusion model significantly overestimates electron scattering by incoherent, but intense whistler-mode waves. We discuss the results obtained in the context of simulations of long-term electron flux dynamics in space plasma systems.

1. Introduction

Resonant electron interactions with whistler-mode waves are one of the main drivers of electron pitch-angle scattering and acceleration in various space plasma systems, e.g. solar flares (Bespalov, Zaitsev & Stepanov 1991; Filatov & Melnikov 2017; Melnikov & Filatov 2020), solar wind (Tong et al. 2019; Cattell et al. 2020, 2021; Mozer et al. 2021), shock waves (Hull et al. 2012; Wilson et al. 2013; Oka et al. 2017; Page et al. 2021), planetary radiation belts (Li et al. 2021; Menietti et al. 2021; Thorne et al. 2021) and magnetic reconnection regions (Le Contel et al. 2009; Breuillard et al. 2016; Zhang et al. 2018a). The basic theoretical framework for description of such interactions is the quasi-linear model (Drummond & Pines 1962; Vedenov, Velikhov & Sagdeev 1962; Andronov & Trakhtengerts 1964; Kennel & Engelmann 1966) that is based on the assumption of weak perturbation of particle dynamics by each single resonance. This assumption reduces the Vlasov equation to the Fokker–Planck diffusion equation (Drummond & Pines 1962; Vedenov et al. 1962) where the main characteristics of wave–particle resonant interactions are diffusion rates. The requirement of a weak perturbation of particle trajectories for a single resonance is equivalent to the requirement that each individual wave–particle resonant interaction should not last for a long time (so particle energy/pitch-angle change for a single resonance is sufficiently small), and there are several mechanisms responsible for particle escape from the resonance.

The original quasi-linear diffusion model assumes the broad spectrum of waves resonating with charged particles (Drummond & Pines 1962; Vedenov et al. 1962), when the resonance width in velocity space $\Delta v_R$ is equal to the difference of resonance $v_R$ velocity (determined for the cyclotron resonant conditions) and wave group velocity $v_g=\partial \omega /\partial k$ (where $\omega$ and $k$ are wave frequency and wavenumber). The estimate for the resonance width can be derived from the condition that a change of the resonant particle velocity, $\Delta v_R \sim |v_R-v_g|\Delta k/k$, will remove the particle from the cyclotron resonance (Karpman 1974). The small factor $\Delta k/k$ is determined by the wave spectrum width in wavenumber space, $\Delta k$. This mechanism determines the shortness of an individual resonance and justifies the applicability of the diffusion approximation for modelling the dynamics of the charged particle ensemble (Karpman 1974; Le Queau & Roux 1987; Shapiro & Sagdeev 1997). This description works well for low-amplitude whistler-mode waves resonating with electrons in homogeneous systems (without spatial gradients of the background plasma and magnetic field), e.g. in the solar wind (see review by Verscharen et al. (2022) and references therein).

The assumption of background magnetic field homogeneity, however, does not work for many space plasma systems. Resonant electron scattering by whistler-mode waves is often observed in magnetic field traps, regions with a spatially localized minimum of the magnetic field magnitude, where charged particles can be trapped and bouncing. Important examples of such traps are the radiation belt dipole field (Lyons & Williams 1984; Schulz & Lanzerotti 1974) and magnetic holes generated by compressional perturbations on a bow shock (Oka et al. 2019; Hull et al. 2020; Yao et al. 2021). Bouncing within magnetic traps, electrons periodically resonate with whistler-mode waves, and the resonance width for cyclotron resonance in an inhomogeneous field is determined from the condition that a change of the resonant particle velocity (due to the field spatial gradient) $\Delta v_R \sim |\partial v_R/\partial s|/k$ will remove the particle from the cyclotron resonance (Trakhtengerts & Rycroft 2008). If $\Delta v_R$ is finite, the quasi-linear diffusion model works even for monochromatic waves ($\Delta k \to 0$) resonating with electrons in magnetic traps (Karpman & Shklyar 1977; Albert 2001, 2010; Shklyar 2021). Thus, the only condition required for the application of the quasi-linear diffusion model is that the mirror force due to the background magnetic field gradient should be stronger than the Lorentz force of the wave field (Karpman 1974).

The small wave intensity approximation, however, is often violated for whistler-mode waves observed in the highly unstable plasma of shock waves (Zhang et al. 1999; Artemyev et al. 2022) and plasma injections (Zhang et al. 2018b, 2019). Such intense waves may resonate with electrons in a nonlinear regime, including effects of phase trapping and phase bunching (e.g. Nunn 1971, 1974; Karpman, Istomin & Shklyar 1974; Inan, Bell & Helliwell 1978). Although phase bunching is a strongly nonlinear effect (Albert 1993; Bortnik, Thorne & Inan 2008), due to the smallness of the electron energy and pitch-angle changes in a single resonant phase bunching, it can be incorporated as a drift term into the Fokker–Planck equation (see discussion in Artemyev et al. (2014), Gan et al. (2020) and Allanson et al. (2021)). Changes of electron energy and pitch angle due to phase trapping are comparable with the initial energies/pitch angles (Omura, Furuya & Summers 2007; Summers & Omura 2007), and thus it is not clear how to include this effect into the Fokker–Planck equation. Several approaches with different integral operators describing the phase trapping contribution to the electron flux dynamics have been proposed (e.g. Omura et al. 2015; Artemyev et al. 2018b; Vainchtein et al. 2018; Hsieh, Kubota & Omura 2020), but the evaluation of such operators is computationally expensive and significantly changes the Fokker–Planck equation. Thus, it is important and practically useful to propose an approach for incorporation of nonlinear effects without significantly altering models based on the Fokker–Planck equation.

The principal possibility for such an approach has been proposed in Solovev & Shkliar (1986): namely, the total contribution of trapping and bunching may compensate each other. This idea has been reinvestigated in Mourenas et al. (2018), where effects of wave modulations were taken into account. Spacecraft observations (e.g. Oka et al. 2019; Zhang et al. 2019, 2020b; Foster, Erickson & Omura 2021; Artemyev et al. 2022) and numerical simulations (e.g. Nunn & Omura 2012; Katoh & Omura 2016; Demekhov, Taubenschuss & Santolık 2017; Tao et al. 2020; Zhang et al. 2021) show that intense whistler-mode waves mostly propagate in a form of short, modulated wave packets. Typical wave packets include only a few wave periods (see figure 1), which can be an effect of sideband instability of wave generation (Nunn 1986) or overlapping of several waves with close wave frequencies (Zhang et al. 2020b; Nunn et al. 2021). Such modulation reduces the efficiency of phase trapping (Tao et al. 2012b, 2013), and can make the net effect of electron resonant interactions with waves more diffusive (Zhang et al. 2020a; Allanson et al. 2020, 2021; Gan et al. 2020; An, Wu & Tao 2022; Mourenas et al. 2022). Thus, the derivation of diffusion rates is the main question for theoretical description of such regime of wave–particle interaction.

Figure 1. Examples of typical wave packets of whistler-mode waves captured by THEMIS spacecraft (Angelopoulos 2008) in the Earth bow shock (a), foreshock transient (b), outer radiation belt (c) and plasma injection region (d). These events are picked up from statistics published in Artemyev et al. (2022), Shi et al. (2020), Zhang et al. (2018b) and Zhang et al. (2018a).

In this paper, we propose an approach for the evaluation of diffusion rates including nonlinear effects for intense, but strongly modulated waves. First, in § 2 we describe the concept of the diffusion coefficient model. Then, in § 3 we provide the basic model equations for the diffusion rate and evaluate diffusion rates for arbitrary wave intensity. Finally, in § 4 we show the contribution of nonlinear effects to diffusion rates averaged over wave intensity distributions and discuss the obtained results.

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 1984) 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. (2022). 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 1984). 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 1962). 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 1988):

(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 1988):

(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. 2019, 2021). 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 2017a,b; Artemyev et al. 2021b; Zhang et al. 2022). 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. 2013; Zhang et al. 2020a; An et al. 2022; Mourenas et al. 2022).

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. 2011, 2012a; Allanson et al. 2020). After the wave amplitude reaches a certain threshold (depending on the electron energy, pitch angle and system characteristics; see Omura et al. 1991; Shklyar & Matsumoto 2009), 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 2010), 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 (2006) and Artemyev et al. (2018a) 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 (2000) and Albert et al. (2022)):

(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 2010). 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$.

4. Discussion and conclusions

In this study, we derive the diffusion rate for electrons resonantly interacting with intense whistler-mode waves. Although such intense waves may resonate with electrons nonlinearly, the efficiency of this interaction would be significantly reduced by wave modulation (Zhang et al. 2020a; Tao et al. 2012b; Allanson et al. 2021; An et al. 2022; Gan et al. 2020, 2022). Indeed, most of the intense whistler-mode waves that are observed in space plasma systems (like Earth's radiation belts, bow shock, foreshock transients and plasma injections) are presented in the form of short well-modulated wave packets (see examples in figure 1 and references in the figure caption).

Using test particle simulations, we show that such modulation will reduce the difference between energy changes of phase trapping and phase bunching electrons, and make the energy change distribution more symmetric (see figure 3). In the limit of total symmetrization of energy change distribution (extremely modulated wave packets), the main (and the only one) characteristic of electron resonant scattering will be the diffusion rate describing the dispersion of this distribution. Such a diffusion rate can be evaluated analytically: the dispersion of the energy changes $\langle (\Delta \gamma )^2\rangle _\xi$ tends to $\langle \Delta \gamma \rangle _\xi ^2$ for large wave amplitudes. That is to say, $\langle (\Delta \gamma )^2\rangle _\xi \to \langle \Delta \gamma \rangle _\xi ^2 \propto B_w$ for $a\gg 1$ (i.e. when wave amplitude is sufficiently large).

Using this theoretical result, we can calculate the ratio of such a nonlinear diffusion rate and the quasi-linear rate, $\langle (\Delta \gamma )^2\rangle _\xi ^{QL} \sim B_w^2$. We can extrapolate to the large-wave-amplitude limit using a renormalization:

(4.1)\begin{equation} \langle(\Delta\gamma)^2\rangle_\xi^{QL}(B_w) =\left(B_w/B_{w,\min}\right)^2\langle(\Delta\gamma)^2\rangle_\xi^{QL}(B_{w,\min}), \end{equation}

where $B_{w,\min }$ corresponds to the $a\ll 1$ limit. The ratio $\langle (\Delta \gamma )^2\rangle _\xi /\langle (\Delta \gamma )^2\rangle _\xi ^{QL}$ should show how quasi-linear diffusion models overestimate the diffusion rates for intense waves, because such models scale $\langle (\Delta \gamma )^2\rangle _\xi ^{QL}$ with wave intensity $B_w^2$. Both the numerator and the denominator of this ratio should be weighted with the actual distribution of observed wave intensities, ${\mathcal {F}}(B_w/B_0)$, and we use two distributions of whistler-mode wave packets collected in the inner magnetosphere.

Figure 7(a) shows two examples of ${\mathcal {F}}(B_w/B_0)$: the main difference between these distributions is in the definition of wave packets used in the two statistics (see details in Zhang et al. (2019, 2020b)). Using these distributions, we plot the ratio $\langle \langle (\Delta \gamma )^2\rangle _{\xi }/\langle (\Delta \gamma )^2\rangle _{\xi }^{QL} \rangle _{\varepsilon }$ as a function of energy and pitch angle in figure 7(b,c). The region with $\langle \langle (\Delta \gamma )^2\rangle _{\xi }/\langle (\Delta \gamma )^2\rangle _{\xi }^{QL} \rangle _{\varepsilon } \approx 1$ corresponds to the dominant contribution of waves with insufficiently large wave amplitude, where $a<1$ for most part of $\varepsilon$, and the diffusion rate is $\langle (\Delta \gamma )^2\rangle _{\xi } \propto B_w^2$. The region with $\langle \langle (\Delta \gamma )^2\rangle _{\xi }/\langle (\Delta \gamma )^2\rangle _{\xi }^{QL} \rangle _{\varepsilon } < 1$ corresponds to the dominant contribution of high-intensity waves, where $a>1$ for a significant fraction of the ${\mathcal {F}}(B_w/B_0)$ distribution, and the diffusion rate $\langle (\Delta \gamma )^2\rangle _{\xi } \propto B_w$. Note that electrons of smaller energy/larger pitch angle resonate with waves closer to the equator, where $a<1$ for the larger part of the ${\mathcal {F}}(B_w/B_0)$ distribution.

Figure 7. (a) Probability distribution functions ${\mathcal {F}}(B_w/B_0)$ of whistler-mode wave intensities from Zhang et al. (2019) (blue) and Zhang et al. (2020b) (red). (b,c) Two-dimensional distributions of $\langle \langle (\Delta \gamma )^2\rangle _{\xi }/\langle (\Delta \gamma )^2\rangle _{\xi }^{QL} \rangle _{\varepsilon }$ in $(E_0,\alpha _{0})$ space for the two ${\mathcal {F}}(B_w/B_0)$ distributions from (a).

Figure 7(b,c) demonstrates clearly that the quasi-linear diffusion model significantly overestimates the real diffusion of electrons of smaller energy/larger pitch angle. Note that a similar effect of diffusion rate reduction relative to the quasi-linear theory predictions has been obtained for broadband waves (see Tao et al. 2011, 2012a). This overestimation will be stronger for active geomagnetic conditions with higher wave intensity (Meredith et al. 2003, 2012; Agapitov et al. 2013, 2018), because $\langle (\Delta \gamma )^2\rangle _{\xi }/\langle (\Delta \gamma )^2\rangle _{\xi }^{QL} \propto 1/B_{w}$. Figure 7(b,c) is plotted under the assumption that there is no net contribution of nonlinear resonant phase bunching and phase trapping effects, and that all wave–particle interactions can be solely described by diffusion because the wave field is dominated by well-modulated short wave packets. The opposite assumption consists of a dominant role of phase trapping and phase bunching effects of electron resonant interactions with coherent long wave packets. For long-term electron dynamics, such effects also can be fitted by diffusion, but the diffusion rate will be much larger than the quasi-linear one (Artemyev et al. 2021a).

The schematic of figure 8 generalizes both regimes of wave–particle resonant interactions for large-amplitude waves: $\langle (\Delta \gamma )^2\rangle _{\xi } \propto B_{w}$ for well-modulated short wave packets and $\langle (\Delta \gamma )^2\rangle _{\xi } \propto B_{w}^{1/2}$ for highly coherent long wave packets. For Earth's radiation belts, the quasi-linear simulations of long-term electron flux dynamics generally describe observed electron fluxes with a reasonable tuning of the averaged wave intensity (e.g. Thorne et al. 2013; Li et al. 2014; Drozdov et al. 2015; Ma et al. 2018; Allison et al. 2021). This will not be possible if either of the two limiting cases shown in figure 8 would work. Therefore, we conclude that wave–particle resonant interactions include both diffusion by short wave packets and nonlinear phase trapping/bunching by rare long wave packets, and a fine balance of these two regimes results in electron diffusion that can be mimicked by the quasi-linear diffusion models. However, for each specific event, such mimicking would require a proper tuning of wave intensity. This underlines the importance of investigations of nonlinear resonant interactions for accurate inclusion of the net effects of phase trapping and phase bunching into wave–particle interaction models.

Figure 8. A schematic view of diffusion rate scaling with wave amplitude normalized to the amplitude threshold for nonlinear resonant interactions, $B_w^{*}$. Quasi-linear diffusion rate $D_{QL}\propto B_w^2$ works only for $B_w/B_w^{*}<1$. For $B_w/B_w^{*}>1$ and incoherent waves should work $D\propto B_w$, whereas for coherent waves should work $D_{NL}\propto \sqrt {B_w}$ (see Artemyev et al. 2021a).

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.