Particle‐in‐Cell Experiments Examine Electron Diffusion by Whistler‐Mode Waves: 2. Quasi‐Linear and Nonlinear Dynamics

Test particle codes indicate that electron dynamics due to interactions with low amplitude incoherent whistler mode‐waves can be adequately described by quasi‐linear theory. However there is significant evidence indicating that higher amplitude waves cause electron dynamics not adequately described using quasi‐linear theory. Using the method that was introduced in Allanson et al. (2019, https://doi.org/10.1029/2019JA027088), we track the dynamical response of electrons due to interactions with incoherent whistler‐mode waves, across all energy and pitch angle space. We conduct five experiments each with different values of the electromagnetic wave amplitude. We find that the electron dynamics agree well with the quasi‐linear theory diffusion coefficients for low amplitude incoherent waves with (Bw,rms/B0)2≈3.7·10−10, over a time scale T of the order of 1,000 gyroperiods. However, the resonant interactions with higher amplitude waves cause significant nondiffusive dynamics as well as diffusive dynamics. When electron dynamics are extracted and analyzed over time scales shorter than T, we are able to isolate both the diffusive and nondiffusive (advective) dynamics. Interestingly, when considered over these appropriately shorter time scales (of the order of hundreds or tens of gyroperiods), the diffusive component of the dynamics agrees well with the predictions of quasi‐linear theory, even for wave amplitudes up to (Bw,rms/B0)2≈5.8·10−6. Quasi‐linear theory is based on fundamentally diffusive dynamics, but the evidence presented herein also indicates the existence of a distinct advective component. Therefore, the proper description of electron dynamics in response to wave‐particle interactions with higher amplitude whistler‐mode waves may require Fokker‐Planck equations that incorporate diffusive and advective terms.


Introduction
Whistler-mode waves play a significant role in the acceleration and loss of electrons within the Earth's radiation belts Thorne, 2010). Physics-based models most typically use the long-established quasi-linear diffusion theory (Kennel & Engelmann, 1966;Lerche, 1968;Lyons, 1974;Summers, 2005) to model particle dynamics due to interactions with whistler-mode waves (e.g., see Albert & Bortnik, 2009;Beutier & Boscher, 1995;Glauert et al., 2014;Su et al., 2010;Subbotin et al., 2010). Observations of whistler-mode waves in the inner magnetosphere (e.g., see Breneman et al., 2011;Cattell et al., 2008;Li et al., 2011;Tyler et al., 2019) have demonstrated not only that very high wave amplitudes (with electric component) ≥100 mV/m exist, but that they are not uncommon (Cully et al., 2008;Kellogg et al., 2011;Watt et al., 2017Watt et al., , 2019Wilson III et al., 2011;Zhang et al., 2019). Theoretical calculations and numerical experiments have demonstrated that such large-amplitude waves should cause electron dynamics to evolve in a manner quite different from that dictated by quasi-linear theory (e.g., see Albert, 2002;Bortnik et al., 2008;Bell, 1986;Karpman et al., 1974;Mourenas et al., 2018;Nunn, 1971;Omura et al., 2007). Indeed, recent observations have also shown that the electron-whistler-mode wave interaction can deviate significantly from that expected by the use of quasi-linear diffusion theory: (i) estimates using data from the THEMIS and the Van Allen Probes satellites suggest that 10-15% of chorus whistler-mode wave packets have wave amplitudes sufficiently high so as to interact nonlinearly with relativistic electrons ; (ii) using data from the Arase satellite, Kurita et al. (2018) showed that deformations in the electron distribution function due to wave-particle interactions occur at a rate that is faster than that expected from quasi-linear theory, such as may be found in .
Much theoretical work on the applicability (or otherwise) of linear and quasi-linear theory has focused on electron interactions with (quasi-)monochromatic whistler-mode wave packets. Inan et al. (1978) studied electron pitch angle diffusion and precipitation due to field-aligned coherent VLF waves at L=4 (the 5 kHz Siple transmitter, Antarctica) using a test particle code and compared the results to a linear analysis of the equations of motion. In their case, it was found that the the linear theory used broke down for wave amplitudes above the relatively low threshold of 3 pT. In particular, they developed a quantitative criterion (based upon an "inhomogeneity parameter") that determined the applicability of linear theory: based upon the conclusion that nonlinear effects should be included if time spent in resonance is greater than the trapping period. This theory was developed by, for example, Bell (1984), to include nonzero wave-normal angles, and further by Solov'ev and Shklyar (1986); by Albert (1993) using a Hamiltonian analysis in a slab magnetic field model; by Omura et al. (2008) to include the effect of time-varying frequency; and by Nunn and Omura (2015) to consider the self-consistent nonlinear interaction in oblique whistlers for any order of resonance. A review of many other relevant works that consider inhomogeneities of the plasma and ambient magnetic field is given in Shklyar and Matsumoto (2009). Test particle numerical experiments conducted by Tao, Bortnik, et al. (2011) have indicated that electron dynamics due to interactions with low amplitude and incoherent field-aligned whistler-mode waves in a uniform background magnetic field can be adequately described by quasi-linear diffusion theory in certain cases. Building on their previous work,  indicated that the quasi-linear diffusion theory may overestimate the value of diffusion coefficients for sufficiently large whistler-mode wave amplitudes (in the case of bounce-averaged diffusion coefficients). We note that there has been much recent work on electron interactions with whistler-mode chorus waves (Omura et al., 2008), in particular, the nonlinear electron interactions that can be difficult or impossible to describe using the standard quasi-linear theory (e.g., see Artemyev et al., 2018;Gan et al., 2020;Mourenas et al., 2018;Vainchtein et al., 2018;Zhang et al., 2018Zhang et al., , 2012. A number of these work discuss dimensionless parameters (e.g., the inhomogeneity parameter or a related parameter) that can predict the dominance of either nonlinear or quasi-linear dynamics, for a given situation (e.g., see a simplified version in Bortnik et al., 2008).
Most of the works that consider numerical experiments of nonlinear wave-particle interactions consider (i) the (in some respects) more physically realistic scenario of background magnetic fields with curvature (or at least spatial gradients) and hence bounce motion and bounce averaging (e.g., see Gan et al., 2020;Silva et al., 2017;; (ii) test particle codes for which the wave spectrum is exactly specified (e.g., see Mourenas et al., 2018;Silva et al., 2018;Tao, Bortnik, et al. 2011);and (iii) interactions with chorus wave packets (e.g., see Artemyev et al., 2018;Vainchtein et al., 2018). In this work we intentionally study a scenario that may be considered more basic in some respects. That is, we study field-aligned broadband and incoherent whistler-mode waves in a uniform field (equatorial) approximation and with a cold plasma. These conditions are in principle very much like the basic assumptions employed in the derivation of quasi-linear theory as applied to the radiation belts. Therefore, we are able to isolate the effect of increasing the wave amplitude on electron dynamics and compare the results to the fundamental quasi-linear theory. Furthermore, we use a particle-in-cell method to accurately capture as much kinetic physics as possible such as is relevant for the fundamentally kinetic wave-particle interaction mechanism. This paper is organized as follows. Section 2 describes the numerical experiments and electromagnetic wave power spectra. Section 3 includes the analysis of electron dynamics, dynamical components, and associated time scales, as well as the extraction of diffusion coefficients and their comparison to the theoretical prediction. Sections 4 contains discussion and summary.

Outline of the Numerical Experiments
In a recent work, Allanson et al. (2019) benchmarked a boundary value problem particle-in-cell method to analyze electron dynamics due to interactions with whistler-mode wave spectra. In particular, Allanson et al. (2019) (hereafter referred to as "Paper 1") presented novel experimental and analysis techniques (using the EPOCH particle-in-cell code) to (a) excite electromagnetic whistler-mode waves within the interior of the domain, by perturbing a boundary; (b) extract diffusive characteristics of electrons across all energies and pitch angles, through the use of a distribution of noninteracting test particles embedded within the experiment. In this paper we use the method established in Paper 1 to examine the dependence of the electron response to incoherent whistler-mode waves on the electromagnetic wave amplitude. In particular, we intend to compare the dynamics and directly extracted diffusion coefficients with those obtained from quasi-linear theory. In contrast to some other works (e.g., see Silva et al., 2018; that invoke the spatial gradient of the background magnetic field, we consider a uniform background field approximation appropriate for the magnetic equator, in a similar way to Tao, Bortnik, et al. (2011). In particular, the conditions used in this paper are chosen to isolate the effect of increasing wave amplitude on electron dynamics, with all other conditions chosen to be appropriate for testing the basic quasi-linear theory (which is formulated for a cold plasma with uniform background magnetic field).
In this paper we use version 4.17 of EPOCH to conduct our numerical experiments. In particular we use the one-dimensional version (EPOCH1D), which only has spatial gradients (and hence electromagnetic wave propagation) in the x-direction. All particle and field quantities do have y and z components however.
The five runs to be discussed are listed in Table 1. All runs use n x = 3,587 grid points with spacing Δx ≈ 235 m, giving a total domain length of L x ≈ 843 km. These values of (Δx, n x , L x ) are the same as that used and justified in Paper 1. As described fully therein, we utilize an electromagnetic source at the left-hand boundary (with root-mean-squared value of the magnetic wave component as listed in column 2 of Table 1) to excite a propagating and incoherent spectrum of right-hand polarized electromagnetic whistler-mode waves within the computational domain. These electromagnetic waves have the overwhelming majority of their power situated in the frequency range 0.2f ce < f < 0.4f ce , for f ce the ordinary electron gyrofrequency. These wave-amplitude spectral profiles are similar to those used in Tao, Bortnik, et al. (2011) and Paper 1, for example. Whistler-mode waves are natural eigenmodes of the cold plasma bulk (Stix, 1992) and propagate through the domain with negligible damping.
The Fourier transforms of the E y and B y components for all five runs are plotted in Figure 1. These amplitude spectra are averaged over all run time and space, and we see that electromagnetic wave power remains well-localized to the 0.2f ce < f < 0.4f ce range (more than 96% of the total magnetic wave power in all cases). Each pair of red and blue curves marks the B y and E y power spectrum for each run, in ascending order (Run 1 is the "lowest" pair and Run 5 is the "highest" pair). Note that we get near-identical results for the spectra of the B z and E z components, as should be the case with circularly polarized waves that propagate in the x-direction only. The vertical black lines mark 0.2 and 0.4 of the electron gyrofrequency (which is the vertical green line). Each pair of horizontal red and blue lines marks the average of the B y and E y power within the 0.2f ce < f < 0.4f ce range. For all runs, dominant power for E y , E z , B y , B z is localized to the whistler-mode dispersion curve (Stix, 1992), as in Figure 4 of Paper 1. Further note that the electromagnetic wave amplitude measured within the domain is not equal to that of the electromagnetic source at the left-hand boundary. The reason for this effect is discussed in detail in Section 2.2 of Paper 1 and is due to the coupling mechanism of the electromagnetic source to the bulk plasma. We mainly discuss the magnetic components in the following text for the sake of brevity (and as is customary when discussing quasi-linear theory). However, note that the electromagnetic waves discussed here have self-consistent electric wave field components following |E w | ≈ (c/η)|B w | (for η = η(ω) the frequency dependent refractive index), and the phase relationship as necessary for right-hand polarized whistler-mode waves (Stix, 1992). The squared ratio of magnetic wave power to background magnetic field amplitude, (B w,rms /B 0 ) 2 , is listed in column 4 of Table 1. These values span the wave amplitude threshold discussed by  that marked a transition from validity to invalidity of the quasi-linear diffusion coefficients (for bounce-averaged diffusion coefficients). For example,  found this wave amplitude threshold to be (B w,rms /B 0 ) 2 ≥ 2 × 10 −7 for 10 keV electrons, and approximately (B w,rms /B 0 ) 2 ≥ 5 × 10 −7 for 215.4 keV electrons (see Section 3.2 and Figure 3 of that paper, respectively).
The domain length L x corresponds to approximately 40-65 wavelengths, depending on the frequency of interest. All experiments allow for wave-particle interactions to occur for T ≈ 1,008t ce ≈ 0.26 s. This corresponds to roughly 200-400 wave periods, once again depending on the frequency of interest. All runs are performed at 500 particles per cell per species. We use electron and ion species each with a number density of n e ¼ n i ¼ 10 7 m −3 ; the correct mass ratio of m i /m e = 1,836.2; and isotropic Maxwellian distributions at 0.1 eV. There is a uniform background magnetic field of B 0 = (B 0x ,0,0) = (140 nT, 0, 0). These parameters are chosen to represent conditions in the magnetosphere close to L ≈ 6, except for the absence of a minority warm electron component. The evolution of wave-particle interactions in the presence of a warm anisotropic component that prevents wave damping will be studied in future works. These parameters give the ratio of electron plasma frequency to gyrofrequency as f pe /f ce ≈ 7.2. Each run required approximately 25,000-35,000 cpu hours (depending on the number of nodes used), and these were performed using the Science and Technology Facilities Council "DiRAC" HPC Facility (www.dirac.ac.uk). We also ran smaller test runs on the Reading Academic Computing Cluster at the University of Reading and the Natural Environment Research Council "ARCHER" HPC facility (www.archer.ac.uk).

Electron Dynamics Extracted from Particle-in-Cell Numerical Experiments
For all runs, the electron dynamics are extracted from the experiments using the specific techniques introduced in Paper 1, and so, we only give a limited description here. At the beginning of the wave-particle interaction, all N total ≈ 10 8 electrons to be considered are binned in (relativistic kinetic) energy and pitch angle space. We use 125 bins in energy space spanning the resonant energies from approximately 3 keV < E < 300 keV and 45 bins in pitch angle space spanning 0°≤ α ≤ 90°. Each of these individual 5,625 bins contains N bin = 17,777 electrons. Electron data are dumped roughly every 11 or 22 gyroperiods. Each electron remains identified with its "initial" bin for the duration of the experiment; that is, electrons are not re-binned at each data dump.
In Paper 1 it was shown that in some regions of (E,α) space, the diffusion can proceed at a rate that is not constant in time (sometimes known as anomalous diffusion Bouchaud & Georges, 1990; del Figure 1. Fourier amplitude of the B y (red) and E y components (blue) of the waves within the PiC domain for each of the five runs. Each pair of red/blue curves marks the B y /E y power spectrum for each run, in ascending order (Run 1 is the "lowest" pair and Run 5 is the "highest" pair). Each pair of horizontal red and blue lines marks the average of the B y and E y power within the 0.2f ce < f < 0.4f ce range. Vertical black lines mark the lower and upper bounds of the driven wave spectrum. The vertical green line marks f ce . Red and blue curves have scales/units defined by the left-hand and right-hand axes, respectively. Castillo-Negrete et al., 2004;Metzler & Klafter, 2000;Perrone et al., 2013;Zaslavsky, 2002;Zimbardo et al., 2015), according to (1) in the case of pitch angle diffusion (for example), for a≠1 and {α l } the set of pitch angle values of the l = 1,2,…,N bin electrons at time t, identified with a given bin α 0 . To be clear, these are the electrons with α ≈ α 0 at t = 0 and ⟨…⟩ the ensemble average. In order not to overburden the notation, we shall omit the {…} set brackets for the rest of the paper and take the subscript l on α l or E l to imply the full set of values of the l = 1,2,…,N bin electrons in a given bin.
However, it was also shown that this anomalous diffusion usually appears as an initial transient feature, such that a ≈ 1 after some short time has passed, t anom. . Figure 2 shows the evolution of ⟨(α l −α) 2 ⟩ for all five runs (asterisks) for electrons initially in a given bin (E 0 ,α 0 ) ≈ (7 keV,30°). As the run number (and hence the wave power) is increased, visual inspection of the plots suggests that the time spent in the anomalous diffusion phase decreases, and so t anom. decreases. This feature is common across (E, α) space. As such, we consider linear fits (solid black line) to the data for each run once t anom. has passed, as follows: (a) 250t ce < t < T; and (e) 50t ce < t < T (T ≈ 1,008t ce for all plots). The gradient of this line can be used to determine a value for the "PiC diffusion coefficient" in a given bin, (E,α) ≈ (7 keV,30°), using Equation A3. To be exact, we extract the gradients of straight line fits to the curves (after t anom. has passed) defined by the right hand sides of the following equations: for (E l ,α l ) the electrons in bin (E,α). The gradient of the dashed black line corresponds to the evolution of ⟨(α l −α) 2 ⟩ as is predicted using the PADIE code  based upon quasi-linear theory and obtained using the values of (B w,rms /B 0 ) 2 as listed in Table 1 and the wave-power spectrum as shown in Figure 1. The midpoint of the dashed black line is shifted so that it coincides with the solid black line, in order to allow a comparison. Figures 2a-2d show very good agreement between the quasi-linear diffusion theory and the diffusion coefficient extracted from the PiC experiments, for this bin. However, Figure 2e shows a less good agreement between the theory and the PiC data. This feature will be further discussed in the following subsection.

Effective Extraction of Particle Dynamics-The Localization Problem
Paper 1 stressed the point that one can only extract a meaningful diffusion coefficient (from particle data) for particular values E 0 = ⟨E l (t = 0)⟩, α 0 = ⟨α l (t = 0)⟩, and then describe it as a function of (E 0 ,α 0 ), while the given subset of N bin particles remain well "localized" to (E 0 ,α 0 ). One clear test of "localization" is to count the number of particles that leave their given initial bin. We can simply consider the number of particles, # out,α or # out,E , that have left the pitch angle or energy bin that they were initially identified with. Figure 3 presents for all five runs, calculated over the time 0 < t < T ≈ 1,008t ce : (a) and (b) describe Run 1; (c) and (d) describe Run 2; (e) and (f) describe Run 3; (g) and (h) describe Run 4; and (i) and (j) describe Run 5. A value of P α > 0.5 or P E > 0.5 indicates that more than half of the particles initially identified with a given bin (E 0 ,α 0 ) at t = 0 are now located "outside" of the bin at t = T. In such cases, it would be inappropriate to describe the dynamics of those particles as a function of (E 0 ,α 0 ). The overplotted black curves mark the values of energy and pitch angle that are in "n = −1 resonance" (see Allanson et al., 2019) with waves of frequency 0.2f ce ("dash"), 0.3f ce ("solid"), and 0.4f ce ("dash-dot"). We can observe the following: (i) runs with higher amplitude waves give higher values of P α and P E ; (ii) P α > 0.5 is observed for Runs 3-5 and P E > 0.5 is observed for Runs 4-5; (iii) significant values of P α and P E are largely confined to the resonant regions, except for Run 5, which also interestingly includes drift in nonresonant regions; and

Journal of Geophysical Research: Space Physics
(iv) higher values of P α are more common than P E , indicating that pitch angle drift is more significant than energy drift in this case.

Journal of Geophysical Research: Space Physics
where α bin,max and α bin,min represent the upper and lower bounds of the given pitch angle bin. Equation 7 states that ⟨Δα⟩=⟨Δα⟩(E,α) is equal to the change over time 0 < t < T ≈ 1,008t ce in the mean value of pitch angle of the l=1,2,…,N bin electrons initially identified with bin (E 0 ,α 0 ). Equation 8 then normalizes this quantity by the size of that pitch angle bin, to give the dimensionless quantity ⟨Δα⟩. If j⟨Δα⟩ðE; αÞj exceeds 0.5, then the mean value of the electron pitch angles has changed by more than half of the initial bin width, and so, ⟨α l ⟩ no longer remains within the initial bin (with pitch angle value α 0 ). One can construct similar quantities ⟨ΔE⟩ and ⟨ΔE⟩ in order to track the dynamics in energy space.
The dashed red line in Figures 2a-2e is the evolution of ⟨ΔE⟩, and the dotted red line is the evolution of ⟨Δα⟩. We see that these values remain very small in Figures 2a-2d. However, in Figure 2e, ⟨α l ⟩ is seen to rapidly grow beyond 0.5. Therefore, the behavior of those particular electrons cannot fairly be described as a function of the bin (E 0 ,α 0 ) ≈ (7 keV,30°), over the entire time frame 0 < t < T ≈ 1,008t ce . We argue that this is a significant cause of the poor agreement shown between the theoretical and obtained diffusion coefficients in Figure 2e. and (f) describe Run 3; (g) and (h) describe Run 4; and (i) and (j) describe Run 5. The magnitude of ⟨Δα⟩ is typically seen to be larger than ⟨ΔE⟩, with the magnitude of both increasing with run number, that is, wave amplitude. Furthermore, it is clear that the dominant changes are observed within the region outlined by the resonance curves. Therefore, we can deduce that these ⟨Δα⟩, ⟨ΔE⟩ dynamics are due to resonant wave-particle interactions. j⟨Δα⟩j can exceed 0.5 for Run 3, while both j⟨Δα⟩j and j⟨ΔE⟩j can exceed 0.5 for both Runs 4 and 5.

Interpreting the Electron Dynamics: Advection and Diffusion
At first glance, the nonzero values of ⟨Δα⟩ and ⟨ΔE⟩ would seem to indicate the existence of an "advective" component to the electron dynamics, in energy and pitch angle space (see the appendix for more discussion of the statistical treatment of both diffusive and advective processes). However, as we shortly discuss, it is not straightforward to disentangle the advective and diffusive components of the dynamics, when analyzing experimental data. Therefore, we term these nonzero values of ⟨Δα⟩ and ⟨ΔE⟩ as "apparently advective" at present, pending further investigation. Furthermore, the common features observed between Figures 3 and 4 suggest that these "apparently advective" dynamics play a significant role in the de-localization of electrons from their initial bins.
At this point it is worth discussing the mathematical tools most commonly used in this field. A statistical description of radiation belt electron dynamics is most usually considered by using a Fokker-Planck equation of the form in Equation A7 (e.g., see , Schulz & Lanzerotti, 1974. Such equations fundamentally describe diffusive dynamics, characterized by the diffusion matrix D ij . However, while the underlying dynamics that are consistent with Equation A7 are diffusive in nature, it does not necessarily mean that one might not observe "apparently advective" dynamics, such as in our experiments. One can have, for example, net acceleration/energization of a given particle subpopulation. In Equation A7, such "apparently advective" dynamics originate from nonzero gradients in the diffusion matrix, that is, terms of the form for F the electron distribution function and J i the action integrals (Roederer & Zhang, 2013) (see the appendix for more detail). Therefore, it is important for us to try and better understand the nature of any dynamics that appear to be nondiffusive, that is, nonzero values of ⟨Δα⟩ and ⟨ΔE⟩. Do these Figure 5. The pitch angle diffusion coefficient D αα calculated from numerical experiment Run 1 (R1) over the time duration 250t ce < t < T ≈ 1,008t ce (as also considered for a particular bin in Figure 2a). The overplotted black curves mark wave-particle resonance curves as in Figure 4.
nonzero values have their origin in fundamentally diffusive dynamics, in a manner as per Equation 9, and/or are they due to fundamentally advective dynamics (e.g., see Mourenas et al., 2018;Gan et al., 2020;Vainchtein et al., 2018;Zhang et al., 2018), that are likely nonlinear in origin? We will now analyze the results from our numerical experiments in order to (i) try to answer these questions and (ii) extract meaningful diffusion coefficients and then compare with those predicted using the standard quasi-linear formalism. < t < T ≈ 1,008t ce ; (c) shows for Run 2 over the time duration t anom. < t < 664t ce ; (e) shows for Run 3 over the time duration t anom. < t < 321t ce ; (g) shows for Run 4 over the time duration 0 < t < 46t ce ; (i) shows for Run 5 over the time duration 0 < t < 46t ce . The overplotted black curves mark wave-particle resonance curves as in Figure 4.

Diffusion Coefficients
First, we use the method outlined above and in Paper 1 to construct pitch angle, energy, and mixed diffusion coefficients from all five runs over the full-time t anom. < t < T ≈ 1,008t ce , in each of the 125 × 45 = 5,625 (E,α) bins. In Figure 5 we plot D αα,PiC obtained for Run 1 as an example. The overplotted resonance curves indeed outline the regions of (E,α) space in which we see the dominant effects, as should be expected for resonant interactions. Empty (or white) regions of the plot indicate regions for which there was negligible diffusion, as defined in Paper 1. In Figure 6 we present comparisons of D αα,PiC with the diffusion coefficient obtained using the PADIE code by plotting For reasons that will be justified shortly, we present the calculation of D αα,PiC over a variety of time scales. Figures 6a-6i present D αα,ratio as follows: (a), (b), (d), (f), and (h) show D αα,ratio for Runs 1-5 over the time duration t anom. < t < T ≈ 1,008t ce ; (c) shows D αα,ratio for Run 2 over the time duration t anom. < t < 664t ce ; (e) shows D αα,ratio for Run 3 over the time duration t anom. < t < 321t ce ; (g) shows D αα,ratio for Run 4 over the time duration 0 < t < 46t ce ; (i) shows D αα,ratio for Run 5 over the time duration 0 < t < 46t ce . When possible, the anomalous diffusion time scales that are used are the same as for Figure 2, that is, t anom. ∈{250tce ,200t ce ,150t ce ,100t ce ,50t ce } for Runs 1-5, respectively. However, the diffusion coefficients that are plotted in Figures 6g and 6i do not allow for the treatment of these anomalous diffusion time scales, since they are calculated over a total time 46t ce , which is less than the smallest anomalous time scale. All plots are saturated to highlight the level 0:5 ≤ D αα;ratio <2: We consider that there is good agreement between the PiC and PADIE diffusion coefficients only for regions in (E,α) space for which this (fairly standard) "factor-of-two" condition is satisfied. We remind that only Runs 1 and 2 allow for both ⟨E⟩; ⟨α⟩<0:5 and P α ,P E < 0.5 across all considered (E,P) space over all 0 < t < T ≈ 1,008t ce . The different timeframes considered in Figures 6c, 6e, 6g, and 6i present alternate calculations of D αα,PiC , as performed over 0 < t < T end , for T end < T, and T end a different value in each case. In particular, the timeframes considered in Figures 6e and 6g (for Runs 3 and 4, respectively) are chosen since they are the maximum possible timeframes that allow for both ⟨E⟩; ⟨α⟩<0:5, and P α ,P E < 0.5 across all considered (E,α) space. Therefore, we can confidently describe electron dynamics as a function of the "initial bin," (E 0 ,α 0 ), over those timeframes. For Run 5 it is not possible to both satisfy these localization conditions and have enough data points in time to extract a diffusion coefficient. However, we present Figure 6i as a best possible estimate.
When considered over the "entire" time t anom. < t < T ≈ 1,008t ce (i.e., Figures 6a, 6b, 6d, 6f, and 6h), we can make the following observations: 1. We see excellent agreement for the lowest amplitude case (Run 1), with almost the entire region satisfying Equation 11. 2. Runs 2,3, and 4 show a mixture of very good and less good agreement, with different regions of energy and pitch angle space not satisfying Equation 11. 3. There is a mixture of very good, less good, and poor agreement for Run 5, with significant regions of (E,α) space such that D αα,ratio < 0.5 and even < 0.25 and < 0.125. 4. If there is not good agreement between the PiC and PADIE diffusion coefficients, then this is almost always due to D αα,PiC < 0.5D αα,PADIE .  also observed the feature that diffusion coefficients extracted from numerical experiments tended to be smaller than those obtained from quasi-linear theory (see Figure 2in that paper), although they studied bounce averaged diffusion coefficients with a test particle code The important and interesting point is that when considered over these appropriately shorter time scales (i.e., Figures 6c, 6e, 6g, and 6i), the agreement with the quasi-linear theory is markedly improved. These alternate calculations of the diffusion coefficient show very good agreement with PADIE across the vast 10.1029/2020JA027949

Journal of Geophysical Research: Space Physics
majority of resonant (E,α) space. This result might be considered as unexpected, since we see good agreement in particular for values of (B w,rms /B 0 ) 2 from 3.73·10 −10 to 1.77·10 −6 . This range spans values for which one might typically expect the inclusion of nonlinear dynamics to prevent diffusive dynamics from being accurately described by quasi-linear theory. Note that while the time scale considered for Figure 6b does not break the de-localization constraints (described in Section 3.1), we present an alternate calculation in Figure 6c over a shorter timeframe, and we see a better agreement between the < t < T ≈ 1,008t ce ; (c) shows for Run 2 over the time duration t anom. < t < 664t ce ; (e) shows for Run 3 over the time duration t anom. < t < 321t ce ; (g) shows for Run 4 over the time duration 0 < t < 46t ce ; (i) shows for Run 5 over the time duration 0 < t < 46t ce . The overplotted black curves mark wave-particle resonance curves as in Figure 4.

Journal of Geophysical Research: Space Physics
PiC and PADIE diffusion coefficients. This further highlights that improvements can be obtained when considering dynamics over shorter time scales. Figures 7 and 8 present D Eα,ratio and D EE,ratio , respectively (defined analogously to Equation 10, for completeness). These plots are saturated as in Figure 6, and the comparisons yield very similar conclusions to those made for the pitch angle diffusion coefficients. < t < T ≈ 1,008t ce ; (c) shows for Run 2 over the time duration t anom. < t < 664t ce ; (e) shows for Run 3 over the time duration t anom. < t < 321t ce ; (g) shows for Run 4 over the time duration 0 < t < 46t ce ; (i) shows for Run 5 over the time duration 0 < t < 46t ce . The overplotted black curves mark wave-particle resonance curves as in Figure 4.

10.1029/2020JA027949
Journal of Geophysical Research: Space Physics Table 2 presents some important quantities that we have extracted from our numerical experiments, over the reduced time scales 0 < t < T end , with T end ≤ T and T end different in each case. For each run, we list the maximum obtained value of D αα,PiC , D EE,PiC , as well as the maximum value of the following quantities for T end the end time of the analysis in each case. μ α,max and μ E,max characterize the maximum values of nondiffusive dynamics across the (E,α) domain. The final column in Table 2 lists the shortest dynamical time scale present in the five runs, defined by In all cases the shortest dynamical time scale is defined by 1/D ααPiC . However, it is clear from Table 2 that 1/μ α,max also defines a short time scale and therefore is an important component of electron dynamics. Figure 9presents the limiting time scales defined by each of the quantities considered above, as a function of the magnetic wave power (listed for each run in Table 2). We see that there is a clear separation between the time scales defined by dynamics in pitch angle (blue) and energy space (red), with pitch angle dynamics clearly the dominant feature, and as is expected for the plasma conditions that we consider, f pe ≫f ce (Summers, 2005). Pitch angle diffusion (blue asterisks) seems to dominate over pitch angle advection (blue triangles). However, it is interesting to see that energy advection (red triangles) seems to dominate over energy diffusion (red asterisks). Furthermore, it is evident that all time scales appear to scale as ∝ 1= B 2 w;rms . This is an expected feature for the diffusive component of the dynamics, at least for the standard quasi-linear theory (Kennel & Engelmann, 1966).
It is known that in general the time scales associated with electron dynamics are also significantly influenced by the background magnetic field inhomogeneity. Coherent large amplitude (chorus-type) waves in the presence of a nonuniform background magnetic field have been shown to cause advective energy dynamics, ⟨ΔE⟩, to be proportional to ðB 2 w;rms Þ 1=4 (e.g., see Albert et al., 2013;Mourenas et al., 2018;Zhang et al., 2018), that is, a time scale that is proportional to 1=B 1=2 w;rms . Furthermore, the diffusive energy dynamics, ⟨(ΔE) 2 ⟩, in such intense wave cases have been shown to alter such that they have associated time scales proportional to 1/B w,rms . Figure 9 does not display these trends, rather it shows the nonlinear behavior and time scales associated purely with the presence of large amplitude waves, that is, the relaxation of the quasi-linear assumption. Future work will determine the additional changes to the wave-particle interaction in the presence of magnetic field inhomogeneities.

Summary and Discussion
In this paper we have used the novel experimental and analysis methods established in Allanson et al. (2019) to analyze electron dynamics due to wave-particle interactions with broadband and incoherent whistler-mode waves, as a function of root-mean-square wave amplitude. Our findings are intended to provide a better understanding of the applicability of quasi-linear diffusion under the usual hypotheses of sufficiently low wave amplitude and broad wave spectrum, and in the simple case of a constant background magnetic field, for which the theory is actually derived (Kennel & Engelmann, 1966). It is the focus of this work to isolate these circumstances and explore the applicability of quasi-linear theory specifically for the case of increasing wave amplitude. In terms of direct physical applications, the results may be considered relevant to energy and pitch angle diffusion over short time scales at (or close to) the magnetic equator (i.e., not considering bounce motion or bounce averaging), and the results are likely to change when more complicated configurations are included.
Before proceeding, we note that in an inhomogeneous ambient magnetic field (and with plasma conditions consistent with the outer radiation belt), the main nonlinear effects appear to be phase trapping and phase bunching (e.g., see Albert, 2002;Artemyev et al., 2018;Bortnik et al., 2008;Inan et al., 1978;Omura et al., 2008) that can strongly modify particle energy within a very short time scale. We do not consider these effects and instead focus on the advective processes that are solely due to the large amplitude waves. We also note that in general the value of any nonlinear wave amplitude threshold is determined by a combination of the properties of the wave spectrum (we consider incoherent waves) and by the gradients of the background magnetic field (e.g., see Albert, 1993;Karpman, 1974). We have described why we ignored the magnetic field inhomogeneity, but it should be expected that any wave amplitude thresholds obtained here would be different for systems with background magnetic field gradients (e.g., see Solov'ev & Shklyar, 1986).
The main conclusions of this paper are as follows: 1. For the lowest amplitude waves, we see excellent agreement between the electron dynamics extracted from the numerical experiment (over the full experimental time 0 < t < T ≈ 1,008t ce ≈ 0.26 s) and the theoretical prediction of quasi-linear theory. 2. In particular, we observe minimal nondiffusive dynamics for the lowest amplitude case and therefore would consider a Fokker-Planck equation containing only diffusive terms (Schulz & Lanzerotti, 1974) as an appropriate model to describe the observed dynamics. 3. As the wave amplitude is increased, we observe that the diffusive component of the electron dynamics agrees less well with the quasi-linear theory. This less good agreement occurs when we naively extract dynamics over 0 < t < T ≈ 1,008t ce . 4. We argue that the apparent less good agreement is in fact due to a naive treatment of the dynamics, in which electron dynamics are not able to be meaningfully ascribed to the intended region of energy and pitch angle space. Specifically, as the wave power is increased, a variety of dynamics can result in electrons no longer being "localized" to the region of energy and pitch angle space that they were originally identified with. 5. We observe that this de-localization correlates well with increasing wave amplitude. As the wave amplitude is increased, we observe that there is a significant component to the dynamics that appears to be advective in nature, as well as enhanced diffusion. 6. The separation of "fundamentally diffusive" and "fundamentally advective" processes is very difficult, since gradients in the (theoretical) diffusion coefficient can lead to "apparently advective" dynamics, that is, nonzero movement of the mean pitch angle or energy of a given electron subpopulation, in energy and/or pitch angle space. 7. When we constrain the time scale over which we extract electron dynamics to 0 < t < T end , with T end ≤ T and different in each case (such that electrons remain "localized"), we see much better agreement between the diffusive component of the electron dynamics and the theoretical prediction of quasi-linear theory. These reduced time scales are of the order of tens or hundreds of gyroperiods. 8. However, we still observe nonzero and sometimes rather significant contributions to the electron dynamics that are of a nondiffusive nature, even over these shorter time scales. 9. These results suggest that nondiffusive (i.e., advective) processes may be important to include when modeling electron dynamics due to wave-particle interactions with higher amplitude whistler-mode waves. The inclusion of truly advective dynamics could require a more complete version of the Fokker-Planck equation, such as that discussed in the Appendix A (following the treatment discussed in Zheng et al., 2019, as applied to EMIC waves). 10. Furthermore, even if "only" considering diffusive dynamics, we suggest that it may be necessary to further investigate the (formal) limiting time scales imposed by large values of the diffusion coefficients, when using numerical Fokker-Planck models of electron dynamics in the radiation belts.
One of the most interesting (and perhaps unexpected) results is that the quasi-linear theory can describe the diffusive component of the dynamics surprisingly well even when in the regime for which one might expect nonlinear effects, given the results found for nonuniform background fields (Tao, Bortnik, Albert, et al., 2012a), but crucially only when one carefully considers important time scales. Of course, a complete description of the dynamics would also require any nondiffusive (i.e., advective) components of the electron motion included. Without including both the diffusive and advective components, one would have only a partial description.
In order to make conclusions about the relative importance of the nonlinear wave-particle interaction for incoherent and broadband waves in radiation belt dynamics over longer time scales than considered herein, the fundamental effect of the background magnetic field inhomogeneity should be properly considered. The work presented here provides very helpful insight into the fundamental nature of electron dynamics in energy and pitch angle space due to the wave-particle interaction under the uniform background-field assumption, used in the basic quasi-linear theory. In particular, the work presented here has allowed us to study the changing importance of advection (which may be described using a Fokker-Planck equation) as a fundamental process when increasing the wave amplitude to nonlinear levels, in the absence of other types of nonlinear dynamics (phase trapping and phase bunching), which may not be described using a standard Fokker-Planck equation (e.g., see Artemyev et al., 2018;Mourenas et al., 2018;Omura et al., 2015;Vainchtein et al., 2018;and references therein).
This study motivates a variety of future numerical experiments and theoretical investigations on the applicability (or otherwise) of quasi-linear diffusion theory to electron dynamics in the Earth's radiation belts, using the particle-in-cell method established in Allanson et al. (2019).

Appendix A: Markov Processes and the Fokker-Planck Equation
Consider a particle variable (e.g., an action integral), X, with a "dynamical history" H n−1 of values X 1 , …, X n −1 at times t 1 , …, t n−1 : H n−1 =(X 1 , …, X n−1 ;t 1 , …, t n−1 ). For a given history, H n−1 , a Markov process is defined by the following: the conditional probability, P n (H n−1 |X n ,t n ), that X lies within (X n ,X n +dX n ) at time t n only depends on the "current" and "previous" states, that is, (X n ,t n ) and (X n−1 ,t n−1 ) (Wang & Uhlenbeck, 1945). This statement can be summarized as P n ðH n−1 jX n ; t n Þ ¼ P 2 ðX n−1 ; t n−1 jX n ; t n Þ: If these stochastic changes in X for a Markov process are "sufficiently small" (e.g., see Reif, 2009;Wang & Uhlenbeck, 1945;Zheng et al., 2019), then the time evolution of the distribution function that describes the ensemble of particles is found to satisfy the Fokker-Planck equation, Particle dynamics under the influence of resonant interactions with waves are commonly treated as Markov processes. If we consider that a particle may be described by a collection of action integrals J i (i=1,2,3), then Equation A1 can be generalized to have the form with the underlying advective processes in "direction" J i governed by components of the vector h i and diffusion described by the symmetric matrix D ij . A particularly clear and thorough presentation of the physical constraints and mathematical background is given in the recent work by Zheng et al. (2019) based upon the derivations in (Reif, 2009;Wang & Uhlenbeck, 1945). The diffusion coefficients and "advection coefficients" in Equation A2 are formally defined according to where ⟨…⟩ stands for the ensemble average over the primed coordinates. Note that Zheng et al. (2019) use the letter B instead of μ. For small stochastic changes, the expression in Equation A3 reduces to the slightly different expression that is perhaps more well known in the literature, such as was used in Tao, Bortnik, et al. (2011), for example.
Due to the manner in which it is constructed, quasi-linear diffusion theory as applied to the Earth's radiation belts is based upon a specific variant of the Fokker-Planck equation (Equation A2), but with h i identically set to zero,

Data Availability Statement
The supporting information provides: (S1) basic instructions on how to run the same experiments that are presented in the main article. Furthermore, data is provided at https://doi.org/10.6084/m9. figshare.12355964. This data file provides: (D1) the contents of the input text files, used for the numerical experiments that are presented in the main article, as further discussed in S1. Therefore, 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.