The Implications of Temporal Variability in Wave‐Particle Interactions in Earth's Radiation Belts

Changes in electron flux in Earth's outer radiation belt can be modeled using a diffusion‐based framework. Diffusion coefficients D for such models are often constructed from statistical averages of observed inputs. Here, we use stochastic parameterization to investigate the consequences of temporal variability in D. Variability time scales are constrained using Van Allen Probe observations. Results from stochastic parameterization experiments are compared with experiments using D constructed from averaged inputs and an average of observation‐specific D. We find that the evolution and final state of the numerical experiment depends upon the variability time scale of D; experiments with longer variability time scales differ from those with shorter time scales, even when the time‐integrated diffusion is the same. Short variability time scale experiments converge with solutions obtained using an averaged observation‐specific D, and both exhibit greater diffusion than experiments using the averaged‐input D. These experiments reveal the importance of temporal variability in radiation belt diffusion.

Fokker-Planck model, the diffusion coefficients contain all the important subgrid physics; hence, our choice of methods to model the diffusion coefficients becomes a key part of the model. Models of diffusion coefficients largely employ two different strategies: parameterized "average" models and event-specific types. Parameterized models offer comprehensive spatial coverage, and are often parameterized by a geomagnetic activity index. For example, the energy and pitch-angle diffusion due to whistler-mode hiss or chorus has been parameterized by Kp (e.g., Spasojevic et al., 2015;Wang et al., 2019) or AE (e.g., Horne et al., 2013;Meredith et al., 2018Meredith et al., , 2020.
The availability of new data sets (e.g., from the NASA Van Allen Probes and the JAXA Arase mission) have encouraged the construction of "event-specific" models of diffusion coefficients. For example, event-specific models of whistler-mode waves can include not only the variation of wave characteristics, but the variation of the magnetic field and number density that controls the resonant condition between waves and electrons (e.g., Ripoll et al., 2017Ripoll et al., , 2019Ripoll et al., , 2016Tu et al., 2014;Zhao et al., 2018) or indeed the wave amplitudes themselves (e.g., Malaspina et al., 2016Malaspina et al., , 2018. Event-specific models capture more of the variability in wave-particle interactions (see e.g., Ripoll et al., 2017), and perhaps importantly, more of the extremes in this process than averaged parameterized models. In what follows, we investigate the effects of including temporal variability in a more general diffusion model, where we build a stochastic parameterization (Berner et al., 2017) from a statistical database of observations. We focus on pitch-angle diffusion due to plasmaspheric hiss (e.g., Hartley et al., 2018;Li et al., 2015;Malaspina et al., 2017Malaspina et al., , 2018.

Constructing D αα (t)
In Figure 1, we show the probability distribution function of observation-specific bounce-averaged pitch-angle diffusion coefficients D αα (X i , Y i ) constructed using 4 years of observations from Van Allen Probe A. Individual pairs of observations X i , Y i of ratio of plasma frequency to electron gyrofrequency (X i = ω pe /Ω e ) and wave intensity at f = 252 Hz (Y i = δB 2 ) are used as input to the calculation of each D αα (α). Diffusion coefficients are shown for E = 0.5 MeV. Full details of other inputs and method of calculation are described in Watt et al. (2019). Figure 1a demonstrates a large amount of variability in D αα (α). Where the diffusion coefficient is weakest (60° ≲ α ≲ 75°), the variability is highest. There is a region at α ∼ 75° where the resonant condition is not satisfied, and so there is no diffusion across this pitch-angle range. The diffusion coefficient is rarely greater than 10 −5 s −1 (similar to a diffusion time scale of 1 day) across all pitch angles.
It is possible to construct an estimate of the average diffusion in two ways (Horne et al., 2013;Shprits et al., 2009;Wang et al., 2019;Zhu et al., 2019): (i) using the average values of input conditions D αα (〈X i 〉, 〈Y i 〉) where 〈…〉 indicates an arithmetic mean (cf. methods of (Horne et al., 2013;Shprits et al., 2009;Wang et al., 2019;Zhu et al., 2019) or (ii) constructing an average of each observation-specific diffusion coefficient 〈D αα (X i , Y i )〉. In Figure 1a, the white solid line indicates the value of 〈D αα (X i , Y i )〉, and the dashed line Two example probability density functions for α = 30° and α = 60° are shown in Figure 1b on a logarithmic scale, indicating the non-Gaussian nature of the distributions. Figure 1c shows the quantile-quantile plot for the probability distribution function of D αα (α = 45°), representative of all the probability distribution functions shown in Figure 1a. The deviations from standard normal quantiles (red line) toward the righthand side of the plot indicate the presence of a statistically heavy tail.
We construct a variable D αα (t) that preserves the underlying distribution and range of D αα as displayed in Figure 1. We choose a simple method where the function D αα (α) is randomly selected from the 2,377 calculated values and kept constant for a period t = Δt (note the importance of retaining the functional dependence of D αα on the pitch angle). At the end of this period, another D αα (α) is chosen and kept constant for the same time period. The modeled time series is denoted by D rand (Δt), and our choices of Δt are motivated in the following section.

Time Scales
To construct D rand (Δt), an estimate of the temporal scales of variability of D αα is required. Since orbiting spacecraft that sample the outer radiation belt (e.g., NASA Van Allen Probes, JAXA Arase, CRRES) traverse this region relatively quickly, distinguishing spatial and temporal variations is difficult. We obtain estimates for the temporal scale of variability by studying the temporal evolution of inputs to the diffusion coefficients (magnetic field strength B 0 , number density n e , and wave intensity δB 2 at f = 252 Hz), and the diffusion coefficients themselves, calculated in our previous work (Watt et al., 2019). It takes around 2 min for the spacecraft to traverse the bin, and around 9 h to return to a similar spot. Note that since location bins are very small in L*, MLT, and magnetic latitude λ, the spacecraft is not guaranteed to traverse the same bin each orbit; it can be >9 h from one data point shown in Figure 2 to the next.
Panel (a) shows the variability of the ambient magnetic field B 0 . The standard deviation of the magnetic field is due to the spatial variation of B 0 as the spacecraft traverses the radial distance 2.95 < L* < 3.05. The variation between passes is very small compared to the variation during a pass, with the exception of the pass on November 18, 2012. For n e (panel (b)), the variation from one pass to the next is much greater than the variation seen within each pass. There is a moderate geomagnetic storm with minimum Dst reaching ∼110 nT that commences near midnight on November 13, 2012 (see gray trace in Figure 2b). Prior to this time, n e mainly lies between 1.5 and 2 × 10 9 m −3 , but after the storm starts, n e is depressed to between 0.5 and 1.5 × 10 9 m −3 , and remains so at least until the end of the month. The variations in wave intensity δB 2 are demonstrated in panel 2(c). The wave intensity is significantly more variable than B 0 or n e , varying over at least 3 orders of magnitude over the interval presented (cf. Figure 5 in Ripoll et al. (2017)). Interestingly, WATT ET AL. the variation of hiss power during the storm appears opposite to that reported by Malaspina et al. (2018) (i.e., that plasmaspheric hiss power decreases with decreasing density) although storms may prove to display different wave power dependencies than in quiet times. The variation in δB 2 between spacecraft traverses through the bin is usually much bigger than the variation seen during each pass. Figure 2d shows the variability in the calculated D αα (E = 0.5 MeV, α = 30°) during November 2012 (calculations of mean and standard deviation are for log 10 (D αα )). The variability of D αα during each pass should be interpreted as an estimate of the uncertainty in its calculation. The variation in diffusion coefficient largely tracks the variation in the wave intensity, although the increase in diffusion is much more pronounced in the middle of the month than the increase in wave intensity. This is due to the systematic decrease in electron number density in this bin as a result of the storm as pitch-angle diffusion increases when the density decreases when all other parameters are kept constant (e.g., Glauert & Horne, 2005).
We conclude from Figure 2d that temporal scales of variability of D αα in this observation bin lie between 2 min and 9 h. We therefore choose two temporal scales Δt = 2 min and Δt = 6 h for our study to span this range. We stress that these choices serve an illustrative purpose in the following numerical experiments and reflect the constraints provided by spacecraft coverage (cf. Ripoll et al., 2017). We wish to investigate whether pitch-angle diffusion depends on Δt; a full characterization of the dependence of diffusive processes on the full range of possible Δt is a larger task left for future work.
WATT ET AL.

Numerical Diffusion Experiments
For the illustration presented here, we assume that pitch-angle diffusion dominates, and ignore energy diffusion due to plasmaspheric hiss. One-dimensional diffusion experiments at a single energy are performed with different choices of bounce-averaged diffusion coefficients D αα (t). The energy E = 0.5 MeV is chosen because the original calculations of the distribution of D αα were performed at this energy. Results may also vary with energy, which will be a focus of future work.
The time evolution of the phase space density at each pitch angle, f, can be found by solving the 1-D Fokker-Planck equation for pitch-angle diffusion, given by where T(α) is given by The second term on the right-hand side of Equation 1 accounts for losses due to atmospheric collisions and the loss time scale τ L is taken to be a quarter of the bounce period inside the loss cone and infinite outside (Shprits et al., 2008).
In the following analysis, Equation 1 is solved using an explicit time stepping scheme in steps of 0.1 s. The pitch-angle grid has a resolution of 1° and boundary conditions at α = 0° and α = 90° are required to define the calculation domain. We assume that far into the loss cone, collisions are sufficient for the phase space density distribution to be isotropic, and use f    = 0 at 0° and 90° (see e.g., Glauert et al. 2014). All of the following experiments initialize the simulation with an isotropic pitch-angle distribution, assuming an electron flux of 5 × 10 3 cm −2 s −1 sr −1 keV −1 for all pitch angles; this distribution is then allowed to evolve over a 30-day period.
We run a series of numerical experiments using the diffusion coefficients displayed in Figure 1. Specifically, each numerical experiment employs the same initial and boundary conditions, and they differ through the choice of D αα (t): 1. Ensemble experiments with D αα (t) = D rand (Δt) with Δt = 2 min and 6 h. The ensembles each contain 60 individual scenarios where a different random selection of D αα (α) is drawn from the observed distribution (see Figure 1). (Ensemble convergence is demonstrated in supporting information). 2. Constant diffusion experiment where D αα (t) = D αα (〈X i 〉, 〈Y i 〉) = constant (represented by dashed line in Figure 1a). 3. Constant diffusion experiment where D αα (t) = 〈D αα (X i , Y i )〉 = constant (represented by solid line in Figure 1a).
The number of unique diffusion coefficients in the distribution is 2,377, defined in each case for 0 < α < 90°, and so given uniform sampling of the distribution, the Δt = 2 min ensemble is likely to have sampled all coefficients after ∼3.3 days and the Δt = 6 h ensemble after ∼400 days. Experiments are run to 30 days, with the caveat that the Δt = 2 min ensemble may show signs of oversampling of the discrete distribution. The experiments presented here do not distinguish between the amount or nature of temporal variability on different time scales, even though in reality variations on ∼2 min time scales are likely to have different causes than those on ∼6 h time scales. Our aim is to investigate the consequences of variability, and how the solutions to the diffusion equation depend on time scale with all other factors treated equally.

Results
In general, each member of the ensemble experiments evolves less smoothly compared to the case where a constant diffusion coefficient is used. An example is shown in Figure 3, where panel (a) displays the evolution of f during a single ensemble member from the numerical experiment where the diffusion coefficient is WATT ET AL.

10.1029/2020GL089962
randomly varied D αα (t) = D rand (Δt = 6 h), and panels (b) and (c) show the evolution of f during the averaged experiments with D αα (t) = 〈D αα (X i , Y i )〉 = constant and D αα (t) = D αα (〈X i 〉, 〈Y i 〉) = constant, respectively. Five days of evolution are shown. Given the initial and boundary conditions of the numerical experiments, and important features of the diffusion coefficients, it is expected that f(α) will rapidly approach 0 as α tends to 0, and the values of f(α ≲ 75°) will decrease as the diffusion progresses.
In the constant diffusion coefficient case, the variations in f(α) as a function of time are very smooth. There is more diffusion in panel (b) than in panel (c), reflecting that 〈D αα (X i , Y i )〉 > D αα (〈X i 〉, 〈Y i 〉). The scenario where the diffusion coefficient is randomly varied often seems quite flat, with sporadic sudden changes in the evolution of f across all α, notably at t ∼ 1.7, 2.5, 3.5, and 4.8 days. The sudden changes indicate times when a large D αα has been chosen from the distribution, and the solution after 30 days in each member of the ensemble is dependent upon the number of large jumps experienced. The times where the solution is practically constant (e.g., from t ∼ 3.6 days to t ∼ 4.8 days) indicate extended periods where mainly very small values of D αα (α) have been chosen. There are also times where the value of f at a constant pitch angle appears to experience brief increases with time (see e.g., α = 20° for 2.6 < t < 3.6 days and for 20 < α < 75 at t = 4.8 days). These occur during times when D αα (α) significantly changes shape from one Δt period to the next, and are most obvious in numerical experiments where Δt = 6 h (i.e., the variation time is slow). During experiments with fast variations of D αα (α), small increases in f are seen only for a very short time.
All solutions tend toward a picture similar to the top of Figure 3a where f(α > 75°) remains at the initial condition, and is elsewhere reduced to small values due to pitch-angle scattering and removal of phase space density in the loss cone. Snapshots of the process after 2, 10, and 30 days can be seen in Figure S1.  Figure 1, we can see that 〈D αα (X i , Y i )〉 > D αα (〈X i 〉, 〈Y i 〉), so it is always expected that the black solid line will lie below the black dashed line in Figure 4a. That is, the diffusion due to the average of all the diffusion coefficients is more than the diffusion due to the diffusion coefficient constructed from WATT ET AL.
10.1029/2020GL089962 6 of 10 The Δt = 2 min ensemble exhibits a time history very similar to the constant D αα (t) = 〈D αα (X i , Y i )〉 experiment, with very little spread in solutions. The median of the Δt = 2 min experiment tracks the 〈D αα (X i , Y i )〉 experiment very well, until t ∼ 20 days. At this point, the median of the ensemble indicates additional diffusion from the constant diffusion experiment and we can see from Figure 4c that the variability of the solutions increases. For a constant diffusion coefficient, the solution is likely to asymptotically approach a limit that takes into account the boundary conditions and is defined by the diffusion lifetime. However, when the shape and strength of the diffusion coefficient varies in time, there is no such limit and no associated lifetime, and diffusion is enhanced.
The time history of the Δt = 6 h ensemble shows slightly slower diffusion initially in the median of the ensemble in Figure 4a, and there is a large spread in the solutions (see Figure 4d). For Δt = 6 h, the median f(α = 30°) decreases at a roughly constant rate until around 25 days and the final values of f(α = 30°) are on average much lower than in the other experiments.

Discussion
We have demonstrated that the evolution of the phase space density in an idealized diffusion experiment depends not only on whether the diffusion coefficient varies with time, but also on the time scale of that variation. There are three notable differences in the ensemble results for Δt = 2 min and Δt = 6 h: differences in the evolution of median f, i.e., the trend in the behavior of the ensemble, differences in the variance of the WATT ET AL.
10.1029/2020GL089962 7 of 10 ensembles, and differences in the final values of f between the two ensembles. It is important to note that both the ensemble experiments, with Δt = 2 min and Δt = 6 h, experience almost exactly the same time-integrated diffusion (see Figure 4b). Yet both the evolution of each ensemble, and the final state after 30 days, is markedly different and depends upon the time scale Δt. Most importantly, for Δt = 6 h, most members of the ensemble experiment experience much more diffusion and reach lower values of f(α = 30°) than the two similar solutions (i.e., the ensemble result with Δt = 2 min, and the averaged result with D αα (t) = 〈D αα (X i , Y i )〉). We have begun investigations into why there are such differences in the two ensembles. We note that a key feature of our ensemble experiments is that the shape of the diffusion coefficient D αα (α) changes with time, in addition to the strength of the diffusion coefficient and remind the reader that with a temporally varying D αα (α), there is no longer a well-defined lifetime.
It is important to note that the averaged 〈D αα (X i , Y i )〉 result is a good approximation of the ensemble median when time scales of variation are short. The diffusion coefficient constructed from averaged inputs, which we have denoted D αα (〈X i 〉, 〈Y i 〉), is much slower than the average diffusion rate estimated by averaging all of the individual diffusion coefficients together, which we denoted 〈D αα (X i , Y i )〉. In our experiments, we noted that the difference in f(α = 30°) after 30 days using these two different "averages" could be more than a factor of 10 at large pitch angles (see Figure S1), and around a factor of 3 at lower pitch angles. These differences might depend sensitively on the location of diffusion "gaps" in pitch angle, or on the size of the variance in the D αα (α) distribution, or a combination of both. Our results strongly suggest that databases of many individual diffusion coefficients should be constructed from colocated and simultaneous measurements. Diffusion coefficient models could then be constructed using appropriate averages of these individual diffusion coefficients, or stochastic methods such as suggested here. Preliminary studies of these new methods of averaging suggest that they are more effective than previous methods (Ross et al., 2020).
Other evidence suggests it is important to understand the underlying distribution of diffusion coefficients.
Idealized numerical experiments using a radial diffusion equation (Thompson et al., 2020) noted that the amount of diffusion depends upon the nature of the underlying distribution of diffusion coefficients. Those distributions with statistically heavier tails experienced greater diffusion, even when the distributions of diffusion coefficients had the same statistical average value.
The ensemble experiments are examples of probabilistic models that can yield rich information about the potential behavior and uncertainty in the physical system. This uncertainty may exist due to a lack of knowledge, indicating that further parameterization of the diffusion coefficients is merited, perhaps using geomagnetic indices or other magnetospheric conditions. Or it may be true that there is inherent natural variability in the system that cannot be parameterized away.

Conclusions
We have presented the results from a series of idealized numerical experiments that highlight the response of the pitch-angle diffusion equation to temporally varying diffusion coefficients that reproduce the full range of observation-specific wave-particle interactions observed over a 4-year period by NASA Van Allen Probes. We present evidence to show that both the wave intensity, and number density observed in the same region of L*, MLT, and magnetic latitude over a period of around 30 days varies significantly, causing changes in D αα of orders of magnitude on time scales of <9 h. We perform idealized numerical experiments of the resulting pitch-angle diffusion that could result from different methods of averaging the diffusion coefficient, as well as two ensemble experiments that deploy stochastic parameterization techniques. If the time scale of variability is very short (Δt = 2 min), then the ensemble result is very similar to the result using the average of many observation-specific coefficients. Where inputs to the diffusion coefficient are averaged prior to its calculation, then the amount of diffusion experienced is much less than in any other numerical experiment. Most interestingly, in the ensemble experiments where diffusion is varied on different time scales, the phase space density solution of the experiment with longer variability time scales (Δt = 6 h) reaches lower values than in the experiment with faster time scales, even though the total time-integrated diffusion in each experiment is the same.
Both this paper and Thompson et al. (2020) highlight that the distribution and variability time scales of diffusion coefficients are important for the evolution of electron phase space density due to diffusion. In other words, key details of the microphysical wave-particle interaction are important for accurate modeling of the macroscale radiation belt system, and the evolution of phase space density is not solely reliant on the average properties of the diffusion coefficients. Our preliminary results isolate pitch-angle scattering due to plasmaspheric hiss, but the concepts illustrated in this paper are likely to be important for all wave-particle interactions in the inner magnetosphere.

Data Availability Statement
Diffusion coefficient data that are shown in Figure 1 can be found at http://dx.doi.org/10.17864/1947.212. Ensemble and averaged numerical experiment results can be found at http://dx.doi.org/10.5281/zenodo.4290006.