Electron Diffusion and Advection During Nonlinear Interactions With Whistler‐Mode Waves

Radiation belt codes evolve electron dynamics due to resonant wave‐particle interactions. It is not known how to best incorporate electron dynamics in the case of a wave power spectrum that varies considerably on a “sub‐grid” timescale shorter than the computational time‐step of the radiation belt model ΔtRBM, particularly if the wave amplitude reaches high values. Timescales associated with the growth rate of thermal instabilities are very short, and are typically much shorter than ΔtRBM. We use a kinetic code to study electron interactions with whistler‐mode waves in the presence of a thermally anisotropic background. For “low” values of anisotropy, instabilities are not triggered and we observe similar results to those obtained in Allanson et al. (2020, https://doi.org/10.1029/2020JA027949), for which the diffusion roughly matched the quasilinear theory over short timescales. For “high” levels of anisotropy, wave growth via instability is triggered. Dynamics are not well described by the quasilinear theory when calculated using the average wave power. Strong electron diffusion and advection occur during the growth phase (≈100 ms). These dynamics “saturate” as the wave power saturates at ≈ 1 nT, and the advective motions dominate over the diffusive processes. The growth phase facilitates significant advection in pitch angle space via successive resonant interactions with waves of different frequencies. We suggest that this rapid advective transport during the wave growth phase may have a role to play in the electron microburst mechanism. This motivates future work on macroscopic effects of short‐timescale nonlinear processes in radiation belt modeling.

theory (Drummond & Pines, 1962;Kennel & Engelmann, 1966;Lerche, 1968;Lyons, 1974;Summers, 2005;Vedenov et al., 1962). Two of the core assumptions of this theory are that: (i) the wave amplitude is a small perturbation to the background field; (ii) and that the wave spectrum is constant in time over the duration for which the diffusion coefficient is calculated -the so-called "limit of resonant diffusion" (Kennel & Engelmann, 1966) (i.e., electromagnetic waves have zero growth rate).
It is now well known that large amplitude whistler-mode waves are not rare in the outer radiation belt Cattell et al., 2008;Cully et al., 2008;Kellogg et al., 2011;W. Li et al., 2011;Tyler et al., 2019;Watt et al., 2017Watt et al., , 2019Wilson III et al., 2011). For example, up to 30% of long-duration chorus waves possess properties that indicate that the quasilinear treatment should fail (Zhang et al., , 2019. These are termed "nonlinear wave-particle interactions." Despite the fact that the waves evident from these observations appear to violate the above assumption, diffusion codes based on the quasilinear theory can yield very good results (e.g., see Glauert et al., 2018;Thorne et al., 2013).
Electron dynamics due to nonlinear interactions with large-amplitude chorus waves (and other large amplitude electromagnetic waves) are in principle markedly different to those that are consistent with the quasilinear theory (Albert, 2002;Albert & Bortnik, 2009;Bell, 1986;Bortnik et al., 2008;Karpman et al., 1974;Mourenas et al., 2018;Nunn, 1971;Omura et al., 2007). One of the major unresolved problems in radiation belt physics is therefore to understand the relative importance of nonlinear wave-particle dynamics in system-scale modeling of the radiation belts, and how to implement these effects. Theoretical and modeling studies indicate that a full treatment of the evolution of electron distributions during nonlinear interactions may require the inclusion of extra "advective" terms in the Fokker-Planck equation in energy and pitch-angle space (Albert & Bortnik, 2009;Allanson et al., 2020;Artemyev et al., 2018;Gan et al., 2020;Lee et al., 2018;Liu et al., 2012;Vainchtein et al., 2018;Zheng et al., 2019), as well as the "diffusive" terms that are the mainstay of numerical models of the radiation belts (Glauert et al., 2014;Schulz & Lanzerotti, 1974).
Test particle and particle-in-cell simulations have indicated that electron dynamics closely follow the quasilinear result for low amplitude whistler-mode and electromagnetic ion cyclotron (EMIC) waves, but that the diffusion coefficient itself "saturates"/"falls-off" as the wave amplitude increases above a certain threshold (Allanson et al., 2020;Liu et al., 2010;Tao et al., 2011Tao et al., , 2012, i.e., above this wave amplitude threshold, the extracted diffusion coefficient is smaller than the quasilinear result. These numerical experiments were performed in the context of a background cold plasma (either directly in the particle-in-cell case, or indirectly via the cold plasma dispersion relation in the test-particle case), such that there is no wave generation via thermal anisotropy instabilities. While this is often a good approximation since the inner magnetosphere is overwhelmingly composed of cold and thermally isotropic plasma, it is well known that "hot" and/or thermally anisotropic electron populations are prevalent (L. Chen et al., 2012;Gao et al., 2014;W. Li et al., 2010) and play a key role in wave and particle dynamics within the inner magnetosphere.
Recent theoretical work has in fact demonstrated the significant effect of hot and anisotropic plasma populations (via modifications to the dispersion relation) on the calculation of diffusion coefficients due to plasmaspheric hiss waves (Cao et al., 2020). Particle-in-cell experiments investigating electron dynamics in energy and pitch-angle space have also demonstrated the strong time-dependence of electron diffusion during the growth and saturation phase of self-consistently generated whistler-mode waves (Camporeale & Zimbardo, 2015). Intriguingly, it was shown that electron diffusion was strongest during the growth phase of the waves, proceeding at a rate that was more rapid than one might expect by naively applying the quasilinear results. Strictly speaking one should not apply the "resonant diffusion" version of the quasilinear theory result during a wave growth phase, but for reasons to be discussed in Section 4, it can be a useful comparison. The authors compared their particle-in-cell results to test-particle codes, since test-particle codes are expected to reproduce quasilinear theory in the low-amplitude case (Tao et al., 2011). Furthermore, the diffusion was weaker during saturation (when wave amplitude is nearly constant in time), and proceeding at a rate slower than that the comparison extracted from the test-particle codes. These results are interesting for two reasons: (i) diffusion was weakest during the wave saturation period, that is, the period of strongest wave amplitude (quasilinear theory predicts diffusion that scales with the square of wave amplitude); (ii) the observed 'saturation' of the diffusion during wave saturation somewhat resembles the aforementioned results obtained by Liu et al. (2010) and Tao et al. (2012). To be specific, whilst the experiments considered in Camporeale and Zimbardo (2015) are quite different to those in Liu et al. (2010) and Tao et al. (2012), the interesting resemblance regards the saturation of the diffusion as the wave amplitude increases/is increased.
Kinetic physics numerical codes are the only means available to study self-consistent particle dynamics during nonlinear wave growth and saturation. Here, we use a particle-in-cell code to study the effects of: (i) increasing background thermal anisotropy; (ii) and increasing wave amplitude on electron dynamics due to interactions with whistler-mode waves. In particular, we wish to discover: (i) whether or not the quasilinear theory comprehensively describes the electron dynamics; (ii) under what circumstances this breaks down; (iii) the nature of the dynamics if not quasilinear diffusion. The experiments here are designed to complement and build upon other recent studies, in particular, the test-particle experiments by Tao et al. (2011), and the particle-in-cell experiments by Camporeale and Zimbardo (2015) and Allanson et al. (2020), that were briefly discussed above. We combine some features from each of these works and analyze the electron response to a combination of "driven" and self-consistently generated whistler-mode waves with a particle-in-cell code. In particular, we are interested to discover in which cases the wave-driving or the self-consistent internally generated waves dominate.
This study is organized as follows. In Section 2, we outline the setup of the 10 numerical experiments and their wave properties. In Sections 3 and 4, we analyze the results from experiments with "lower" and "higher" thermal anisotropy respectively. Sections 5 and 6 conclude with a discussion, summary, and an indication of possible future work.

Outline of the Numerical Experiments
We analyze electron dynamics in response to interactions with electromagnetic waves using the EPOCH particle-in-cell (PiC) code (Arber et al., 2015). We use the boundary-value-problem experimental method as introduced and benchmarked in Allanson et al. (2019). This method was developed to analyze the electron response to whistler-mode waves in the context of a kinetic plasma physics code, via the extraction of electron "tracer" particle dynamics in energy and pitch-angle space. In Allanson et al. (2020), this method was used to analyze the particular effect of increasing whistler-mode wave amplitude on electron dynamics, in the presence of a thermally isotropic background cold plasma and uniform background magnetic field. The nature (and justification) of the set-up of the numerical experiments to be used here almost exactly replicates those used in Allanson et al. (2020) (e.g., boundary conditions, number of particles per cell, total domain length, run-time, spatial and temporal resolution, wave-driving mechanism), and so we shall omit most of this discussion for the sake of brevity. The most important physical parameters for the wave-particle interaction are chosen as follows: electron and proton number density (n e = n p = 10 7 m −3 ), appropriate for L ≈ 6; equatorial approximation of the uniform background magnetic field   ( , , ) 140 x y z B B B nT  ,0,0 ; ordinary electron gyrofrequency f ce ≈ 3919 Hz; electron plasma to gyrofrequency ratio f pe /f ce ≈ 7.2 (note that these are the same as the parameters used in Tao et al. [2011]). We use a realistic proton-electron mass ratio m p /m e = 1, 836.2. Each run required approximately 20,000-30,000 cpu hours (depending on the number of nodes used), and these were performed either using the UK Research and Innovation (UKRI) Science and Technology Facilities Council "DiRAC" facility (www.dirac.ac.uk), or the UKRI Natural Environment Research Council "ARCHER" facility (www.archer.ac.uk).
As is a very common approach in such numerical experiments (e.g., see Allanson et al., 2019Allanson et al., , 2020Camporeale, 2015;Camporeale & Zimbardo, 2015;H. Chen et al., 2017H. Chen et al., , 2018Gao et al., 2017;Kuzichev et al., 2019;Liu et al., 2010Liu et al., , 2012Tao et al., 2011Tao et al., , 2017, we consider one field-aligned dimension (x), which is a reasonable approximation near the magnetic equator (i.e., uniform fields and field-aligned wave propagation) (Tsurutani & Smith, 1977), i.e., the region that is most important for whistler-mode wave generation (Gao et al., 2014;W. Li et al., 2010Teng et al., 2018). We therefore use the one-dimensional version of the EPOCH code, such that all quantities have spatial gradients in the x-direction only. This is often referred to as a "3V1D" configuration (i.e., three dimensions in velocity space and one dimension in real space). We note that on the one hand, the limitations imposed by our uniform magnetic field idealization mean that we can't realistically simulate the generation and subsequent development of nonlinear rising-/falling-tone chorus waves in an inhomogeneous field appropriately. Nonlinear chorus generation is usually (Hikishima et al., 2020;Omura et al., 2008), but not always (Wu et al., 2020), considered to require background field inhomogeneities. However, the idealized nature of our numerical experiments does allow us to more accurately compare the experiment results to the original fundamental quasilinear theory (e.g., Kennel & Engelmann, 1966), developed under the uniform magnetic field assumption. We hope that any lessons learned from these experiments will guide us in future works that may include background gradients, and therefore electron bounce motion.
Whilst the cold plasma approximation is often a good first-order treatment in the radiation belt context, plasma within the Earth's inner magnetosphere is known not to be exclusively composed of cold and thermally isotropic components. A better approximation to reality would be to consider the addition of some higher temperature and thermally anisotropic electron components (e.g., see L. Chen et al., 2012;Gao et al., 2014;W. Li et al., 2010), which can be thermally unstable and susceptible to self-consistent generation of whistler-mode waves (Gary, 1993). The main difference between this study and that presented in Allanson et al. (2020) is that we consider the effect of allowing the electron component of the background plasma to have minority populations characterized by higher temperatures and thermal anisotropy. This approach allows us to analyze the changing nature of the electron response in different experiments for different values of both: (a) the background anisotropy; (b) the driving wave-amplitude. This scheme builds upon and blends the approaches used in Camporeale and Zimbardo (2015) and Allanson et al. (2020), to give an even more detailed understanding of nonlinear electron dynamics.
The exact initial conditions for the electron species used in the experiments are fully detailed in Table 1. All distribution functions are initialized in our experiments as relativistic bi-Maxwellians if anisotropic, and as relativistic Maxwellians if isotropic. Parallel and perpendicular temperatures are defined using the second-order moments of the distribution (e.g., see Xiao et al., 1998 for a discussion). The parameter that describes the temperature anisotropy is defined by A T = T ⊥ /T ∥ − 1 such that a thermally isotropic species has A T = 0. The electron cold, "cool," "warm," and "hot" species have densities according to 0.9778n p , 0.02n p , 0.002n p , and 0.0002n p respectively (these sum up to n p to ensure quasineutrality). These density parameters and the temperatures shown in Table 1 are chosen to resemble "typical" (Runs I-V) and "higher" (Runs VI-X) levels of thermal anisotropy, respectively. The parameters in Runs I-V were chosen to resemble typical values that are consistent with the findings of W. Li et al. (2010), L. Chen et al. (2012), and the parameters for Runs VI-X were chosen to represent correspondingly "high" values (i.e., two times the value of T ⊥ /T ∥ of Runs I-V). The ion component is initialized as a cold and isotropic Maxwellian species in every run with A T = 0 and T ⊥ = T ∥ = 1 eV.
In the same manner as described in Allanson et al. (2019) and utilized in Allanson et al. (2020), we excite incoherent right-hand polarized and field-aligned whistler-mode waves within the computational domain by perturbing the left-hand boundary with electromagnetic field oscillations with frequencies uniformly distributed in the range 0.2f ce < f < 0.4f ce . Table 2 presents the rms amplitude of the magnetic component of these waves: (a) at the left-hand boundary source; (b) spatially and temporally averaged over the entire spatial domain and run time. This choice of frequency spectrum is intended to be a generic example of "broadband and ALLANSON ET AL.   incoherent waves," and is the same as chosen in the test-particle study by Tao et al. (2011) and Allanson et al. (2019Allanson et al. ( , 2020. Figure 1 presents the time-dependent wave power normalised to the background magnetic field (with   2 2 2 w y z B B B ) for each run. The blue lines depict Runs I-V, and the pink-red lines show the power for Runs VI-X. The vertical black line at t ≈ 0.03s represents the time at which the traveling waves excited within the domain first reach the right-hand boundary. A clear separation in resultant wave power is evident in Runs I-V (which have the lower anisotropy). These runs are therefore dominated by the value of the incident wave amplitude, and no significant thermal instability is triggered. The situation is quite different for Runs VI-X. Despite the fact that these runs are driven with the same incident wave amplitudes as Runs I-V, respectively, they exhibit much higher wave amplitudes within the domain. These runs exhibit clear linear and nonlinear wave growth, and subsequent saturation (reached by t ≈ 0.11s in the "latest" case). Runs VI-X are therefore dominated by the thermal instability. Interestingly, wave growth and saturation occur most quickly for Run X, which is the case with the highest driving wave amplitude. These results demonstrate the well-known fact that for an instability to promote the growth of waves quickly to high amplitudes and nonlinear saturation, one requires a "sufficient" supply of free energy, as determined by the number density and value of A T for the minority anisotropic components (e.g., see Tang & Summers, 2019).
Figure 2 presents various spectral information regarding the B y component of the electromagnetic waves within the domain, for different runs. Since we have a one-dimensional domain, we only have field-aligned wave propagation. Therefore the whistler-mode waves are right-hand polarised, and the spectral information regarding the B y component is qualitatively the same as for the B z component. We use B y (or we could equivalently use B z ) for the Fourier analysis, because this quantity is readily available from the experiment output and contains all information required for field-aligned waves, that is, a field-aligned whistler-mode wave with frequency f that propagates in the x-direction is a sum of two linearly (y and z) polarised sinusoidal waves (that themselves propagate in the x-direction) with the appropriate phase difference. Note that the "Fourier amplitude" that we shall refer to for B y at a given frequency in Figures 2a and 2d is the amplitude A of one linearly polarized wave mode ∼ A sin(kx − 2πft). The Fourier amplitude values at each individual frequency f are (as should be expected) significantly lower than the value of B w , since B w is formally calculated by integrating wave power (in B y and B z ) over all relevant frequencies. The Fourier amplitude at a given frequency in our experiments is observed to be about one order of magnitude lower than the value of B w , corroborating with previous experiments conducted with similar broadband spectra (e.g., see discussions in Tao et al. [2011Tao et al. [ , 2012 Figure 2a show the average power in the 0.2f ce < f < 0.4f ce range. Runs I-V once again demonstrate clear separation in wave power, as determined by the amplitude of the wave-driving mechanism. Furthermore, we see that wave power within the domain is well-localized to the 0.2f ce < f < 0.4f ce region, matching the frequency range of the driver. However, we see from Figure 2d that the wave power for Runs VI-X: (i) do not exhibit significant separation; (ii) do display significant wave power over a much wider frequency band than for Runs I-V. The thermal instability in Runs VI-X has clearly triggered wave power over the majority of the whistler frequency domain. Figures 2b and 2e show time-frequency contour plots of the power in the B y component for Runs III and VIII, respectively, as examples. Figure 2b further confirms the localization of power to the 0.2f ce < f < 0.4f ce band, whereas Figure 2e shows that the wave power spectra rapidly "spreads" over wider regions of frequency space over a time period that roughly matches that of the instability.
Whistler-mode waves in the magnetosphere frequently (but not always, e.g., see Teng et al., 2019) have a "power gap" at 0.5f ce , separating the "lower-band" (f lh < f < 0.5f ce , for f lh the lower-hybrid frequency) and "upper-band" (0.5f ce < f < f ce ) (Burtis & Helliwell, 1969;Tsurutani & Smith, 1974). However, we do not see that phenomenon in our experiment. This is not unexpected for the following reasons. It is now accepted ALLANSON ET AL.  that Landau damping by whistler-mode waves plays an important role in establishing the power gap (e.g., see H. Chen et al., 2021;J. Li et al., 2019;Ratcliffe & Watt, 2017). Landau damping by whistler-mode waves in this context requires a parallel component of the whistler-mode wave field perturbation vector, and this can only exist if the wave normal angle is non-zero (i.e., obliquely propagating modes with k × B 0 ≠ 0). We note that, whilst this Landau damping mechanism necessarily requires the existence of obliquely propagating whistler-mode waves to occur, the power-gap effect can be observed in both the field-aligned and oblique modes (e.g., see Taubenschuss et al., 2014). However, we only have electromagnetic wave modes that propagate (anti-)parallel to the field-aligned direction in our experiments (k × B 0 = 0). Therefore there is no Landau damping by whistler-mode waves in our experiments.
Once the thermal instability saturates (at roughly t ≈ 0.08s for Run VIII), we see that the wave spectrum stops widening in frequency space, and slowly begins to become more narrow (in a manner somewhat similar to that in Tao et al. [2017], H. Chen et al. [2017Chen et al. [ , 2018, and Kuzichev et al. [2019]). This narrowing occurs as the instability saturates and the anisotropy decreases, reducing the maximum frequency with positive growth rate (i.e., the condition on marginal instability is satisfied in a smaller ranger of frequencies) (Kennel & Petschek, 1966;Ossakow et al., 1972;Tao et al., 2017).  dispersion relation (Stix, 1992) overplotted in cyan. We see that wave power is localized to the whistler band in both cases. These are presented as examples, and this result is common for the other runs. As is now to be expected, wave power is once again localized to the 0.2f ce < f < 0.4f ce band of the whistler-curve for Run III, but exhibits a wider spread in frequency space for Run VIII.
Coherent chirping (rising-tone or falling-tone) chorus waves  are the whistler-mode wave modes that are most usually considered (and modeled) to reach high amplitudes and be responsible for nonlinear wave-particle interactions (see discussion in Section 1). A background magnetic field inhomogeneity is considered to be fundamentally important to promote rising-tone and falling-tone chorus waves to high amplitudes (Omura et al., 2008) (however note recent experiments that demonstrate rising-tone and falling-tone chorus in a uniform magnetic field; Wu et al., 2020). However, less structured whistler-mode waves have been observed by a number of studies to be prevalent in the equatorial/near-equatorial regions outside of the plasmasphere (Gao et al., 2014;Shumko et al., 2018;Tsurutani & Smith, 1974;Tsurutani et al., 2013). These wave modes have been termed "hiss-like" or even "hiss-like chorus" (hereafter we use the term "hiss-like"). They can reach high amplitudes, and they have frequencies in the lower-band chorus range, that is, above the range of plasmaspheric hiss.
For example, W. Li et al. (2012) conducted a survey of whistler-mode waves using near-equatorial Time History of Events and Macroscale Interactions during Substorms (THEMIS, Angelopoulos, 2008) data, and found many instances of quasi-field aligned broadband hiss-like whistler-mode waves with 0.1-0.5f ce lower-band frequencies outside the plasmasphere. They found that the hiss-like bands predominantly in the regions of higher f pe /f ce , including high-amplitude waves (magnetic component > 300pT). Furthermore, Gao et al. (2014) similarly considered THEMIS whistler-mode wave data outside the plasmapause, and correlated occurrence statistics of hiss-like waves with the relative proportion of "hot" plasma ("N h ") to total plasma ("N t "). They found that hiss-like waves were more common than rising tones or falling tones. They also found that ≈ 0.5−1 nT hiss-like waves could be observed, and quoting: "With the increase of N h /N t , the occurrence rate of hiss-like emissions generally increases, and the corresponding wave amplitudes become larger". We also highlight the work of Katoh and Omura (2013), using hybrid numerical experiments of whistler-mode wave generation (both chorus and "broadband hiss-like" emissions). They found that the excitation of broadband hiss-like emission with amplitudes comparable to discrete chorus elements was found in the experiment with lowest field-line inhomogeneity. Similarly, unstructured lower-band whistler-mode spectra have been observed in kinetic physics (spatially one-dimensional) numerical experiments that studied whistler-mode wave generation via anisotropic instability (H. Chen et al., 2017Chen et al., , 2018Kuzichev et al., 2019;Tao et al., 2017).

Lower Anisotropy: Dynamics Dominated by Incident Waves
Theoretical radiation belt quasilinear diffusion coefficients are derived using complicated formulae that depend on the background plasma, number density, and electromagnetic wave spectral properties Kennel & Engelmann, 1966;Lerche, 1968;Lyons, 1974;Summers, 2005). Radiation belt numerical models then solve diffusion equations that utilize these diffusion coefficients. These Fokker-Planck equations describe the statistical evolution of the particle distribution functions, in a manner that is consistent with Markovian stochastic dynamics (Wang & Uhlenbeck, 1945). For a given ensemble of particles that are initially localized at t = t 0 to some initial value of energy and pitch angle (E 0 , α 0 ), the formal (Markovian) definition of the pitch-angle diffusion coefficient (for example) at the point (E 0 , α 0 ) in energy pitch-angle space is given by where 〈…〉 is the ensemble average (i.e., the mean) (Allanson et al., 2020;Liu et al., 2010;Zheng et al., 2019). Note that this definition does not diverge as Δt → 0, since the numerator is itself of order (Δt) 1 (Reif, 2009;Wang & Uhlenbeck, 1945;Zheng et al., 2019). This equation defines the diffusion coefficient as being proportional to the rate of increase of the ensemble mean of the "squared pitch-angle change" (the change is with respect to the initial value, that is, α 0 = 〈α(t = t 0 )〉. This formula is slightly different to the one that is actually more commonly used and discussed, which is given by where 〈α〉 is the mean of the pitch angle values at t = t 0 + Δt. This equation defines the diffusion coefficient as being proportional to the rate of increase of the variance of the ensemble of values. We can observe that the expression in Equation 2 is different to that in Equation 1, by virtue of the fact that and therefore D αα,formal ≥ D αα,approx , since the final term on the RHS of Equation 3 is non-negative. As explained in Zheng et al. (2019) and references therein, for a system with well-defined small stochastic changes, D αα,formal → D αα,approx as Δt → 0, since   as Δt → 0.
In Figure 3, we present the evolution of the quantities      2 0 ( ) ( Figure 3a) and 〈(α −〈α〉) 2 〉 (Figure 3b), for ensembles of particles that are initially localized to bins in energy and pitch-angle space that are resonant with whistler-mode waves at f = 0.3f ce . Each ensemble of electrons has 17,777 members. We start tracking the electron data from the time (t ≈ 0.03s) at which the whistler-mode waves first reach the right-hand side of the physical domain (as discussed in Section 2). The plots are labeled and color-coded for electrons with initial values of pitch angle α 0 = 5°, 15°, 25°, 35°, 45°, 55°, 65°, 75°, 85°, and 89°. The asterisks show data from the numerical experiment, and the straight lines are those that are consistent with the D αα as predicted using the PADIE software  for a Gaussian fit to the average of the wave spectrum. We see that the two figures show near-identical results to each other (despite the different measures used), and very good agreement with the PADIE results. This indicates that treatment with a "standard" quasilinear Fokker-Planck diffusion code (i.e., one with only diffusive terms, e.g., see Schulz & Lanzerotti, 1974) would be expected to yield good results (at least for values of Δt roughly of the order of that used in this experiment, but possibly longer). We have performed similar comparisons with the results from PADIE for all Runs I-V, and find qualitatively similar results to those obtained for the 5 runs conducted with a cold and isotropic plasma background in Allanson et al. (2020) (e.g., see Figure 6 in that paper). Therefore we do not repeat that presentation for the sake of brevity.
To gain further insight into the nature of the particle dynamics, we present scatter plots that track the movement of ensembles of particles in (E, α) space for Runs I and V. Figures 4a-4d each present the trajectories of 17,777 electrons scattering from a given initial (E, α) bin at four snapshots in time from Run I. In each case, the blue rectangle outlines the domain (i.e., the bin) in (E, α) space of the electrons at beginning of the tracking period t ≈ 0.03s (this is labeled as "t = 0" in this Figure  ) ; (b) 〈(α −〈α〉) 2 〉 for electrons with initial values of pitch angle α 0 = 5°, 15°, 25°, 35°, 45°, 55°, 65°, 75°, 85°, and 89°, and energies that correspond to particles resonant with f = 0.3f ce . The asterisks show data from the numerical experiment, and the straight lines are those that are consistent with the D αα as predicted using the PADIE software for a Gaussian fit to the average of the wave spectrum. Each ensemble of electrons has 17,777 members. We start tracking the electron data from the time (t ≈ 0.03s) at which the whistler-mode waves first reach the right-hand side of the physical domain (as discussed in Section 2). defines resonant energies and pitch angles with f = 0.3f ce waves. Underneath each scattering plot, we present a histogram of those particles in pitch-angle space. Each bin represents electrons with energies initially resonant with f = 0.3f ce waves, for pitch angles according to (a) 15°; (b) 35°; (c) 55°; (d) 75°. The vertical black/blue line marks the mean value of the distribution at the initial/current time. The scatter plots show that diffusion is very weak in this case, since particles barely move outside of their initial bin. Likewise, the histograms do not demonstrate significant changes in the individual pitch-angle distributions.
Figures 5a-5d present the same information but for Run V. The overplotted curves are the lines that define resonant energies and pitch angles with waves at frequency f = 0.2f ce , 0.3f ce , 0.4f ce (such that lower energy particles, at a fixed pitch angle, resonate with higher frequency waves). The vertical red line in Figures 5a and 5b are drawn at α = 3°, which is the approximate value of the equatorial loss cone at L = 6. Of course, our experiments are performed in the context of a uniform background magnetic field and we have no particle loss, and so this line is drawn purely to provide an indication of possible dynamics into/out of the equatorial loss cone. In contrast to the results from Run I, we see much more significant diffusion. Within the total run-time, we see that significant proportions of electrons are able to move in (E, α) space so that they are in resonance with waves of a frequency ± 0.1f ce different to that with which they were initially resonant. The histograms show the clear development of initially narrow distributions into broader distributions, as is generally expected in the diffusive paradigm. Furthermore, Figures 5a and 5b (the ensembles initially located at 15°, 35°) do show some transport to pitch angle values within the loss cone.
In Figure 6, we plot the trajectories of 30 randomly chosen electron trajectories for particles with energies and pitch angles initially resonant with waves of frequency f = 0.3f ce for Runs I ((a) and (b), 35° and 75°, respectively) and V((c) and (d), 35° and 75°, respectively). The thin "central" overplotted black line represents the mean value of the 30 electrons, and the thicker black lines represent the mean ±1 standard deviation, which is a proxy measure for diffusion according to Equation 2. The results from Run I further demonstrate the minimal diffusion present over the timescale considered. Furthermore, the mean value of the electrons effectively stays constant. For Run V, we see that an efficient diffusive process is taking place. The electron ensemble initially localized to 75° demonstrates some small level of advection in the mean value, but the diffusive process still appears to dominate.
It is helpful to construct a proxy quantity that compares the "rate of diffusion" with the "rate of advection." We use dimensionless quantities that were also used by Gan et al. (2020) when analyzing nonlinear electron-whistler interactions, namely in energy space and in pitch-angle space, for (E 0 , α 0 ) the initial value that defines a given ensemble. Both the absolute value and the time-derivative of these quantities can yield useful information. Figure 7 plots these quantities to indicate the ratio of advection to diffusion for Run V in energy space (a) and pitch angle space (b). We do not show the results for Run I since the ratios are essentially nil in that case, corroborating with the previous results to show that advective processes are negligible for that case, on this timescale. The color-coded plots show data for electrons with initial values of pitch angle α 0 = 5°, 15°, 25°, 35°, 45°, 55°, 65°, 75°, 85°, and 89°, and energies that correspond to particles resonant with f = 0.3f ce . The results from Run V show that: ALLANSON ET AL.

10.1029/2020JA028793
Figure 4. Tracking the particle scattering from a given initial (E, α) bin at four snapshots in time from Run I In each case the blue rectangle outlines the domain (i.e., the bin) in (E, α) space of the particles at t = 0. The overplotted curve is the line that defines resonant energies and pitch angles with f = 0.3f ce waves. Underneath each scattering plot, we present a histogram of those particles in α space. The vertical black/blue line represents the mean value of the distribution at the initial/current time. Each bin represents electrons with energies initially resonant with f = 0.3f ce waves, for pitch angles according to (a) 15°; (b) 35°; (c) 55°; (d) 75°. Each ensemble of electrons has 17,777 members. "t = 0" in these plots refers to the time (t ≈ 0.03s) at which the whistler-mode waves first reach the right-hand side of the physical domain (as discussed in Section 2), and "t = T" refers to the end of the experiment.
(i) advective processes are present; (ii) pitch-angle advection appears to be stronger than energy advection (as might be expected for a plasma with f pe ≫ f ce in which pitch-angle dynamics are dominant, e.g., see Summers, 2005); (iii) in all cases the results are bounded by ± 1, thus indicating that diffusive dynamics are the dominant process. Interestingly, it seems that electrons with pitch-angle values < 45° experience net positive advection, and that electrons with pitch-angles > 45° experience net negative advection. As a phenomenological interpretation, these results make sense when we observe the evolution of the mean of the histograms in Figure 5, as the ensemble pitch-angle distributions rapidly spread across pitch-angle space. We note that similar results were found in the test-particle studies by Liu et al. (2012) and Lee et al. (2018) for the case of electron interactions with EMIC waves in a uniform background field.

Higher Anisotropy: Dynamics Dominated by Instability
In this section, we consider the results from Runs VI-X which have the higher levels of background anisotropy (see Table 1 for all details). The wave power and spectral profiles are all qualitatively similar for each Runs VI-X since they are dominated by the instability and saturate at roughly the same amplitude. Therefore we only analyze Run VIII in detail as a representative example.   (d)). The thin "central" overplotted black line represents the mean value of the 30 electrons, and the thicker black lines represent the mean ± 1 standard deviation. We start tracking the electron data from the time (t ≈ 0.03s) at which the whistler-mode waves first reach the right-hand side of the physical domain (as discussed in Section 2). electrons with the same color coding as for Figure 3, and also resonant with waves of frequency f = 0.3f ce . The asterisks show data from the numerical experiment, and the straight lines are once again those that are consistent with the D αα as predicted using the PADIE software for a Gaussian fit to the average of the wave spectrum. In contrast to the results from Runs I-V, we see that the quasilinear diffusion prediction does not effectively describe the dynamics over the timescale considered, since they overestimate the amount of diffusion in both cases. We note here that, formally speaking, the variety of quasilinear theory that we compare to (the limit of "resonant diffusion;" Kennel & Engelmann, 1966) should only be applied to time-invariant wave spectra. This is clearly not the case for Runs VI-X, as demonstrated in Section 2, and so it would seem that we should not compare our experimental data in this manner. However, the total run-time of the dynamics considered is ≈ 0.26s, and this is much shorter than the typical timesteps used in radiation belt diffusion codes, which are usually of the order of tens of seconds, minutes, or even longer. Numerical diffusion codes therefore do not model rapid variations in the wave power spectral densities such as those presented here. Instead, they consider the wave power averaged over longer timescales. We are therefore comparing the results from our numerical experiments with the predictions made using the averaged wave power in the quasilinear theory, to more effectively compare with the philosophy and approach that is used in current diffusion codes.
We can immediately make three main observations from considering Figures 8a and 8b. First, the two different measures of diffusive dynamics (Equations 1 and 2) give very different results. This is unlike the case presented in Figure 3 for which they gave nearly identical results. As discussed in Section 3, this difference must be determined by the following quantity ALLANSON ET AL.

10.1029/2020JA028793
13 of 22 Figure 7. Ratio of advection to diffusion for Run I in energy space (a) and pitch angle space (b); and also for Run V in energy space (c) and pitch angle space (d). The plots show data for electrons with initial values of pitch angle α 0 = 5°, 15°, 25°, 35°, 45°, 55°, 65°, 75°, 85°, and 89°, and energies that correspond to particles resonant with f = 0.3f ce . We start tracking the electron data from the time (t ≈ 0.03s) at which the whistler-mode waves first reach the right-hand side of the physical domain (as discussed in Section 2). ) ; (b) 〈(α −〈α〉) 2 〉 for electrons with initial values of pitch angle α 0 = 5°, 15°, 25°, 35°, 45°, 55°, 65°, 75°, 85°, and 89°, and energies that correspond to particles resonant with f = 0.3f ce . The asterisks show data from the numerical experiment, and the straight lines are those that are consistent with the D αα as predicted using the PADIE software for a Gaussian fit to the average of the wave spectrum. We start tracking the electron data from the time (t ≈ 0.03s) at which the whistlermode waves first reach the right-hand side of the physical domain (as discussed in Section 2). Therefore we can deduce that the mean value of electron ensembles change significantly over the timescale considered in pitch-angle space. The difference between these two measures strongly suggests that one needs to be careful when extracting diffusion data from numerical experiments, interpreting diffusion coefficients, and when considering the results in the context of the standard diffusion equation used in radiation belt modeling (Schulz & Lanzerotti, 1974). In particular, with the experiment timescale that we consider (≈0.26s), and also therefore the typical time-step values (Δt RBM ≫ 0.26s) of the magnitude commonly used in radiation belt diffusion models, we cannot appropriately construct a diffusion coefficient. This timescale/ time-step size is too large (with respect to the rapid diffusion and advection) to be used in a numerical approximation of the mean square limit in the definition of the diffusion coefficient in Equation 1 or 2.
Second, and despite the fact that the two plots present very different results, we observe that the nature of the diffusion follows a similar pattern in both cases (in all but the two highest presented pitch angle samples). Namely, a growth period and subsequent saturation period that roughly correlates with the observed wave growth and saturation timescales. As described in Section 1, this feature bears an interesting resemblance to results obtained in Liu et al. (2010), Tao et al. (2012), and Camporeale and Zimbardo (2015). Note that each of those three papers treated different problems in a different context (see Section 1 for a full description), but they all observed a "saturation" in the diffusion due to wave-particle interactions, either for: (i) the case of individual experiments with a constant, but high wave amplitude (Liu et al. [2010] in the case of EMIC waves and Tao et al. [2012] in the case of whistler-mode waves); (ii) or for an experiment in which a whistler wave instability driven by thermal anisotropy underwent growth and saturation phases, in a manner somewhat similar to the experiments that we present (Camporeale & Zimbardo, 2015).
Third, most of the particle ensembles considered in Figure 8b seem to saturate/asymptote toward similar values within the range ≈ 400−460deg 2 . This feature was also noted in the case of the experiment run by Camporeale and Zimbardo (2015), and can correspond to the pitch-angle variance of a 90°-peaked pitch-angle distribution f(α) ∼ sin α.
In Figure 9, we track particle scattering for Run VIII, in the same manner as presented in Figures 4 and 5. Each bin once again represents electrons with energies initially resonant with f = 0.3f ce waves, for pitch angles according to (a) 15°; (b) 35°; (c) 55°; (d) 75°. The overplotted curves are the lines that define resonant energies and pitch angles with waves at frequency f = 0.1f ce , 0.2f ce , 0.3f ce , 0.4f ce , 0.5f ce , 0.6f ce , 0.7f ce , 0.8f ce , 0.9f ce (such that lower energy particles, at a fixed pitch angle, resonate with higher frequency waves). As before, the vertical black, blue and red lines in the histograms depict the mean value at the initial and current times, and the loss cone respectively. In all cases, we see that electrons travel a much greater distance in pitch-angle space, as compared to Runs I and V. This rapid transport (diffusion and advection) allows electrons to cover regions of energy and pitch-angle space that are resonant with waves throughout the majority of the whistler-mode wave frequency range, as demonstrated by the overplotted resonance curves for different frequencies. These dynamics are facilitated by the wideband frequency range observed in the numerical experiment (see Figure 2). Furthermore, all cases demonstrate electron transport within the loss cone.
The pitch-angle distributions presented in Figure 9 of the individual ensembles demonstrate not only that each initial localized ensemble rapidly spreads to cover all pitch-angle space, but that the distribution is 90°-peaked. This is in stark contrast to Runs I-V, in which electron ensembles only spread locally and with a broad distribution. This 90°-peaked feature is common to all four presented cases in Figures 9a-9d. This corroborates with the results presented in Figure 8b (i.e., the asymptotic of the variance toward ≈ 400−460deg 2 ), and those results in Camporeale and Zimbardo (2015).
In Figure 10, we plot the trajectories of 30 randomly chosen electron trajectories from Run VIII for particles with energies initially resonant with waves of frequency f = 0.3f ce with pitch angle 15° (a); 35° (b); 55° (c); 75° (d). Once again the thin "central" black line represents the mean value, and the thicker black lines represent the mean ± 1 standard deviation. For the 15°, 35°, 75° cases we see a significant advective motion of the mean value. We also observe the expected diffusive spreading (i.e., the expansion of the standard deviation "width") in all cases. These diffusive and advective dynamics are seen to be strongest during the wave growth-phase, and slow down during the saturation phase.
Finally, in Figure 11, we present the ratios of energy ( Figure 11a) and pitch-angle (Figure 11b) advective and diffusive dynamics for Run VIII, as was done for Runs I and V in Figure 7, and using the formulas given by Equations 4 and 5. This is done for electrons with initial values of pitch angle α 0 = 5°, 15°, 25°, 35°, 45°, 55°, 65°, 75°, 85°, and 89°, and energies that correspond to particles resonant with f = 0.3f ce . The plots show similar results as discussed for Run V in Figures 7c and 7d), namely that: (i) advective processes are present; (ii) pitch-angle advection is larger than energy advection; (iii) and electrons with smaller pitch-angle values experience net positive advection, while electrons with larger pitch-angles experience net negative advection. However, in contrast to Run V, we see that: (i) the significant dynamics occur during the wave-growth phase and saturate thereafter; (ii) the values reached are significantly larger, with             2 0 | ( ) / ( ) | bounded by 3 for Run VIII, as opposed to 1 for Run V. This indicates that the wave-growth phase has facilitated significant advective dynamics to occur and to dominate on short timescales.

Discussion
Analysis of the electromagnetic wave properties has demonstrated the well-known and crucial role that the thermal composition of a background plasma can have on the resultant wave power and its evolution. In particular, our experiments demonstrate/confirm that for "high" levels of anisotropy: (i) the wave power can reach "high" amplitudes ≈ 1 nT such that B w ∼ O(B 0 /100); (ii) the high amplitude of the final saturation state is essentially unaffected by the lower amplitude incident wave power; (iii) but that saturation is reached more quickly for a higher amplitude of incident wave; (iv) wave power can rapidly spread over wider regions of frequency space over the wave-growth period, but then stop spreading as the wave power saturates (in a manner that is similar to other recent 1D PiC experiments that study whistler-wave growth via anisotropic instabilities (H. Chen et al., 2017Chen et al., , 2018Kuzichev et al., 2019;Tao et al., 2017).
However, for the lower/"typical" levels of anisotropy, we see that wave propagation and spectral properties are relatively unaffected by the background thermal conditions, and thermal instabilities need not occur. When considering the electron dynamics in these low anisotropy cases (Runs I-V), we observe that the quasilinear approach seems appropriate over the short timescales considered. To be precise, the diffusive dynamics are observed to dominate over advective dynamics, and the extracted diffusion coefficients agree well with the predictions of the quasilinear theory obtained using PADIE software, when considered over "appropriately short timescales" (as discussed in much more detail in Allanson et al. [2020]). Similar results are found in other numerical (test-particle) experiments that are consistent with a cold and isotropic plasma background (e.g., Tao et al. [2011] for whistler-mode waves and Liu et al. [2010Liu et al. [ , 2012 for EMIC waves).
In the cases with higher anisotropy, we observe markedly different dynamics. Significant diffusive and advective responses are observed, and these correlate with the wave-growth phase. As the rms wave-amplitude rapidly grows in time, we observe that the spectral width (bandwidth) also increases. By tracking the individual electron dynamics during this time we can see that this allows rapid transport to regions of energy and pitch-angle space that are resonant with very different frequencies to those with which they were initially resonant.
Interestingly, these rapid dynamics contribute to the emergence of 90°-peaked pitch-angle distributions for individual ensembles of electrons. Ensembles that are initially localized to a given pitch-angle, scatter across all 0° < α < 90° pitch-angle space, to give "wedge-shaped" 90°-peaked distributions: (i) that have a functional profile and variance roughly equivalent to that of a sin n α distribution (Camporeale & Zimbardo, 2015); (ii) and such as are known occur in the radiation belt plasma environment (Gannon et al., 2007;Wrenn et al., 1979;Zhao et al., 2020). Note that the distributions that we present describe a given ensemble, and do not necessarily resemble the pitch angle distribution across all energy space.
ALLANSON ET AL.
10.1029/2020JA028793 16 of 22 Figure 9. Tracking the particle scattering from a given initial (E, α) bin at four snapshots in time from Run VIII. In each case, the blue rectangle outlines the domain (i.e., the bin) in (E, α) space of the particles at t = 0. The overplotted curves are the lines that define resonant energies and pitch angles with f = 0.1f ce , 0.2f ce , 0.3f ce , 0.4f ce , 0.5f ce , 0.6f ce , 0.7f ce , 0.8f ce , 0.9f ce waves (such that higher frequency waves resonate with lower energies, for a given pitch angle). Underneath each scattering plot, we present a histogram of those particles in α space. The vertical black/blue line represents the mean value of the distribution at the initial/ current time, and the vertical red line indicates the equatorial loss cone value α ≈ 3°. Each bin represents electrons with energies initially resonant with f = 0.3f ce waves, for pitch angles according to (a) 15°; (b) 35°; (c) 55°; (d) 75°. Each ensemble of electrons has 17,777 members. "t = 0" in these plots refers to the time (t ≈ 0.03s) at which the whistler-mode waves first reach the right-hand side of the physical domain (as discussed in Section 2), and "t = T" refers to the end of the experiment.
Electron microbursts (Anderson & Milton, 1964) are localized and short timescale (of the order of 100 ms) electron precipitation events, with energy spectra from a few keV and up to the MeV range. Microburst events have been shown to correlate with equatorial and near-equatorial whistler-mode waves activity in a general sense, that is, both: (i) coherent rising-tone chorus elements (Breneman et al., 2017;Lorentzen et al., 2001;Mozer et al., 2018); (ii) and also less structured lower-band whistler-mode spectra, including hiss-like spectra (Shumko et al., 2018;Tsurutani & Smith, 1974;Tsurutani et al., 2013). It has also been demonstrated that whilst quasilinear theory results can sometimes happen to characterize the precipitation rate, it is likely that nonlinear wave-particle processes play a fundamental role, especially when wave amplitudes reach ≈ 1 nT values (Breneman et al., 2017;Mozer et al., 2018;Shumko et al., 2018). Strong whistler-mode wave activity (including the less structured hisslike variety) has been observed to correlate with raised levels of energetic electron density and anisotropy (Gao et al., 2014;W. Li et al., 2010), as according to the standard theory (Davidson et al., 1972;Kennel & Petschek, 1966;Ossakow et al., 1972;Tang & Summers, 2019). Here, we show that significant diffusion and advection occur during short-timescale (≈100 ms) dynamics due to interactions with high-amplitude hisslike whistler-mode waves (≈100 nT at saturation), which are generated by strong (thermal-anisotropy driven) instabilities. Of particular note for comparison is the work by Shumko et al. (2018), in which Van Allen Probe data of near-parallel-propagating and near-equatorial hiss-like whistler-mode waves at L ⋆ ≈ 6 (with background electron number density ≈ 10 7 m −3 ) indicated electron microbursts with a duration 150−500 ms. These parameters and characteristics are very similar to those presented in this study. Collecting together all of these discussed results leads us to hypothesize that the wave-growth phase itself may play an important role in the electron microburst mechanism, for the few tens of keV energy range. Since our experiments are conducted in the equatorial/uniform B approximation, this hypothesis will need to be confirmed via future studies that include electron bounce motion due to background magnetic field line inhomogeneities, and therefore the existence of a loss cone.
The results here demonstrate the pivotally important role that anisotropic components can play in electron dynamics, despite the fact that such components may make up a very small portion of the total plasma density. This presents the further motivation for more detailed statistical models of the anisotropy and number density of the fractional higher temperature anisotropic species in the radiation belts. In particular, it would be interesting to know, on a statistical basis: (i) how extreme the levels of anisotropy can be; (ii) how often such "large" values occur; (iii) and what proportion of the total electron density they constitute.

Summary and Future Work
The current limits on computational power dictate that any operational (physics-based) forecasting model of the Earth's radiation belts must be a reduced-physics model, such as Fokker-Planck diffusion codes. However, it is clear that nonlinear wave-particle interaction effects, and more general nonlinear kinetic effects (such as wave growth and saturation) are not only interesting to study, but important to understand for modeling of high energy electron dynamics. The important and underlying ALLANSON ET AL.
10.1029/2020JA028793 17 of 22 Figure 10. 30 randomly chosen electron trajectories from Run VIII for particles with energies initially resonant with waves of frequency f = 0.3f ce with pitch angle 15° (a); 35° (b); 55° (c); 75° (d). The thin "central" overplotted black line represents the mean value of the 30 electrons, and the thicker black lines represent the mean ± 1 standard deviation. We start tracking the electron data from the time (t ≈ 0.03s) at which the whistlermode waves first reach the right-hand side of the physical domain (as discussed in Section 2). question could be summarized: "how can we best incorporate very short-timescale effects (that may include the influence of very high-amplitude waves) into longer timescale global modeling of the radiation belts?" The study here tries to stimulate further discussion on this important topic, and highlight some particularly interesting features of this large and challenging problem.
Using the method established in Allanson et al. (2019), we analyze (in the radiation belt context) the electron response to a combination of externally driven and internally generated whistler-mode waves with realistic background thermal anisotropy, and with a kinetic plasma physics code. We observe that: 1. For our experiments run with "lower/typical" levels of anisotropy, wave properties are relatively unaffected by the background thermal conditions. Furthermore, the quasilinear approach seems appropriate over the short timescales considered (with some caveats in a similar manner to those that are discussed in Allanson et al. (2020), and similarly by Liu et al. (2010) in the case of EMIC waves). 2. For experiments run with the "higher" levels of anisotropy, an instability is triggered with growth and saturation phases. Electron dynamics during this instability are radically different to that predicted using the typical average wave power method in the quasilinear theory. 3. During the wave growth phase the dynamics include both strong diffusive and advective components which both saturate as wave power saturates, with the advective component often dominant. 4. The growth phase facilitates rapid advective changes in electron pitch angle, via successive electron interactions with the (widening) resonance regions in energy and pitch-angle space. 5. We hypothesize that the wave-growth phase for whistler-mode waves that saturate at high amplitudes may itself play an important role in the electron microburst mechanism for the few tens of keV energy range.
The study here is intended to complement other recent studies that explicitly compare the efficacy of the quasilinear diffusion theory with physics-based numerical experiments to describe electron dynamics when it may not be expected to be formally appropriate (e.g., see Allanson et al., 2020;Camporeale & Zimbardo, 2015;Gan et al., 2020;Liu et al., 2010Liu et al., , 2012Tao et al., 2011Tao et al., , 2012Tao et al., , 2013. In particular, this study motivates future work on the macroscopic impact of very short timescale and nonlinear processes in radiation belt modeling, and the effects of anisotropic background plasma components on electron scattering. Future work could include rigorous and detailed comparisons between results from kinetic experiments and numerical diffusion models in energy and pitch-angle space. Furthermore, a more realistic treatment of radiation belt particle and wave dynamics should include particle mirror/bounce motion due to background magnetic field-line inhomogeneities (e.g., see Tao et al. [2012] for example experiments using a test-particle method). We aim to consider this problem using the particle-in-cell method in future works. In particular, it will be interesting to consider the suggestion raised here that the nonlinear electron dynamics (diffusion and advection) during the rapid (≈100 ms) wave-growth phase for whistler-mode waves that saturate at high amplitude (≈1 nT) may contribute to the electron microburst mechanism.
ALLANSON ET AL.
10.1029/2020JA028793 18 of 22 Figure 11. Ratio of advection to diffusion for Run VIII in energy space (a) and pitch angle space (b). The plots show data for electrons with initial values of pitch angle α 0 = 5°, 15°, 25°, 35°, 45°, 55°, 65°, 75°, 85°, and 89°, and energies that correspond to particles resonant with f = 0.3f ce . We start tracking the electron data from the time (t ≈ 0.03s) at which the whistler-mode waves first reach the right-hand side of the physical domain (as discussed in Section 2).

Data Availability Statement
The Supporting Information provides (S1) basic instructions on how to run the same experiments that are presented in the main article. Data file D1 is available at https://doi.org/10.6084/m9.figshare.13653065.v1. This data file provides the contents of the input text files, used for the numerical experiments that are presented in the main article. A combination of the information provided in S1 and D1 will enable readers to locally generate the same experimental data as was considered in the main article.