A New Approach to Constructing Models of Electron Diffusion by EMIC Waves in the Radiation Belts

Electromagnetic ion cyclotron (EMIC) waves play an important role in relativistic electron losses in the radiation belts through diffusion via resonant wave‐particle interactions. We present a new approach for calculating bounce and drift‐averaged EMIC electron diffusion coefficients. We calculate bounce‐averaged diffusion coefficients, using quasi‐linear theory, for each individual Combined Release and Radiation Effects Satellite (CRRES) EMIC wave observation using fitted wave properties, the plasma density and the background magnetic field. These calculations are then combined into bounce‐averaged diffusion coefficients. The resulting coefficients therefore capture the combined effects of individual spectra and plasma properties as opposed to previous approaches that use average spectral and plasma properties, resulting in diffusion over a wider range of energies and pitch angles. These calculations, and their role in radiation belt simulations, are then compared against existing diffusion models. The new diffusion coefficients are found to significantly improve the agreement between the calculated decay of relativistic electrons and Van Allen Probes data.


Introduction
The electrons that make up the outer electron radiation belt exist over a wide range of energies, extending up to well over 10 MeV, where curvature and gradient drift dominates and they complete closed paths around the Earth. Electrons can be rapidly accelerated up to relativistic energies during geomagnetic storms through a combination of wave-particle interactions with whistler mode chorus waves Horne, Thorne, Shprits, et al., 2005), shocks (Blake et al., 1992), and radial transport (Schulz & Lanzerotti, 1974). Relativistic electron fluxes can be increased by several orders of magnitude (Baker et al., 1994;Reeves et al., 2003) on the timescale of hours to days and can pose a risk to satellite instrumentation (Wrenn et al., 2002). These fluxes then decay on timescales of 5-10 days (Claudepierre et al., 2020;Meredith et al., 2006), and it is believed that their removal is strongly influenced by pitch angle diffusion by electromagnetic ion cyclotron (EMIC) waves (Albert, 2003;Horne & Thorne, 1998). Coincident observations of EMIC waves and relativistic electron precipitation support the idea that EMIC waves are influential at scattering MeV electrons (Miyoshi et al., 2008;Rodger et al., 2008;Usanova et al., 2014;Yuan et al., 2018). The effects of EMIC waves can also be seen in the narrowing of ultrarelativistic electron pitch angle distributions Shprits et al., 2016;Usanova et al., 2014) and deepening of local minima in phase space density .
Electron resonance with EMIC waves, and hence diffusion, is sensitive to both the wave spectrum and the plasma properties (Meredith et al., 2003;Summers et al., 2007). For instance, increasing the electron density decreases the minimum resonant energy and allows diffusion at larger pitch angles, providing a mechanism by which near-equatorial mirroring electrons can be removed. Similarly, for wave frequencies that approach the respective ion cyclotron frequency from below, the minimum resonant energy can be reduced down to ∼400 keV (Ukhorskiy et al., 2010). In fact, studies using helium band EMIC waves with Gaussian spectra with central frequencies of 0.9f cHe , where f cHe is the helium ion gyrofrequency, have been found to improve the agreement between radiation belt models of the ≳1 MeV electron flux and observations Ma et al., 2015). It is unclear, however, whether such a spectrum is representative of the EMIC waves in the radiation belts given that mean wave frequency to gyrofrequency ratios observations are on average substantially below these adopted values. More specifically, are EMIC waves with characteristics in the tails of the wave frequency and plasma density distributions sufficiently efficient at diffusing near-equatorially mirroring MeV energy electrons to provide observed levels of electron losses?
In this paper we present a new approach for calculating the average quasi-linear diffusion coefficients for the diffusion of relativistic electrons by averaging over observation-specific bounce-averaged diffusion coefficients. We begin by illustrating the variability of EMIC observations in section 2. In section 3 we then present an alternative method of calculating statistical bounce-averaged EMIC diffusion coefficients by performing bounce-averaged diffusion calculations for each EMIC wave observation before averaging over the observations. These new diffusion coefficients are then compared to existing models using average wave and plasma properties. In section 4, the new coefficients are included into 3-D Fokker-Plank based radiation belt simulations to assess the effect of the new diffusion coefficients on the relativistic electron population.

EMIC Wave Spectra
For this study, we use data from the Combined Release and Radiation Effects Satellite (CRRES) (Johnson & Kierein, 1992). The magnetic field was measured by the 16 Hz fluxgate magnetometer (Singer et al., 1992) and provides the ambient magnetic field from which we determine the gyrofrequencies. Additionally, the field measurements are fast Fourier transformed using 1,600 points with a 400 point step to obtain wave spectra extending from 0 to 8 Hz in steps of 10 mHz at a 25.6 s resolution. Each spectral profile is analyzed, and two Gaussian functions are fitted, one for the hydrogen band and the other for the helium band . Noise spikes and bands that do not have Gaussian profiles are omitted. For each Gaussian profile, we have the peak power spectral density, PSD m , the corresponding peak frequency f m , and width, f w , defined by (1) The plasma densities are derived from the plasma wave spectra measured by the Plasma Wave Experiment  as described in Meredith et al. (2002). We store the central frequency of the fitted Gaussian f m /f cHe , the Gaussian width f w /f cHe and the wave intensity B 2 w , and the corresponding electron plasma frequency to electron gyrofrequency ratio, f pe /f ce , together with L * , the magnetic latitude m , and the associated activity indices (AE, Kp, Dst, and solar wind pressure, P dyn ). Here f cHe is the helium gyrofrequency. The geomagnetic coordinates are calculated using the ONERA-DESP library and the International Geomagnetic Reference Field at the middle of the corresponding year combined with the Olson-Pfitzer quiet time model (Olson & Pfitzer, 1977).
Key parameters in determining wave-particle interactions include the wave intensity, B 2 w (nT 2 ), the proton gyrofrequency, f cp , and the plasma density. In Figures 1a and 1b we show scatter plots of B 2 w (nT 2 ) against the determined f pe /f ce for H + and He + EMIC bands, respectively. For comparison, the constant value of pe ∕ ce = 10 that was used in Kersten et al. (2014) (henceforth Kersten2014) is overlaid. To aid with interpretation we over plot the survival function of f pe /f ce , where the survival function S of parameter x with probability density function, f (x), is defined as (2) Intuitively, for a given X, S(X) can be interpreted as the probability that x exceeds X. The survival function is one minus the cumulative distribution function. Strong wave power is observed over a wide range of f pe /f ce , with pe ∕ ce > 20 ≳10% of the time and reaches ∼30 for both EMIC bands on occasion. Considerable wave power is present over a range f m /f cHe (panels c and d), with the distribution particularly heavy tailed for helium band EMIC waves (panel d). Again for comparison, f m /f cHe and the frequency cutoffs (f m ± 2f w )/f cHe of the Gaussian spectra that are assumed by Kersten2014 are shown in blue. It is clear that the spectrum used by Kersten2014 does not overlap with the f m /f cHe of a substantial fraction of observations. Additionally, each of the EMIC observations has a finite spectral width, further increasing the wave power coverage to higher and lower ratios. Significant wave power is therefore present close to the bounding upper gyrofrequency, which may substantially affect electron diffusion. Finally, the central frequency and frequency range assumed by Ma et al. (2015) (henceforth Ma2015) and Drozdov et al. (2017) for He + EMIC waves is shown in Figure 1d. Their power spectrum lies within the tail of the distribution of fitted f m /f cHe with only ∼1% exceeding their adopted value, suggesting their choice is unrepresentatively close to the helium ion gyrofrequency.

Diffusion Coefficients
In order to determine the effects of EMIC waves on radiation belt electrons, we calculate bounce-averaged diffusion coefficients for each observation using the concurrent spectral and plasma measurements (Watt et al., 2019). These results are then combined to create statistical bounce and drift-averaged diffusion coefficients. Although simulations suggest that EMIC waves can become highly nonlinear (Omura & Zhao, 2012), here we assume that the wave amplitude is sufficiently low and that the waves are sufficiently broad band that quasi-linear theory can be applied. The concurrent measurement bounce-averaged diffusion coefficients for each CRRES observation are calculated using the code PADIE . PADIE solves the cold plasma dispersion relation and so does not include warm plasma effects. We assume there is no EMIC wave power at absolute magnetic latitudes greater than 20 • within L * = 5.5 (Meredith et al., 2003). However, Allen et al. (2015) found EMIC wave at larger absolute latitudes. Those at larger absolute latitudes are mainly linearly polarized and therefore more likely to be oblique. We are assuming the main contribution to diffusion comes from lower latitudes. Hence, our diffusion coefficients are a conservative estimate rather than an over estimate.
For each observation, we perform latitudinally restricted bounce-averaged diffusion coefficient calculations using 10 latitude bins evenly spaced in absolute latitude between the equator and 20 • . The wave power is assumed to be 0 outside the bin while the wave spectrum and plasma properties inside are determined from the local observation mapped to that bin latitude. These latitudinally restricted bounce-averaged diffusion coefficients are then summed to calculate the bounce-averaged diffusion coefficient for that observation. The procedure is explained in more detail below.
For a given observation, the local observed ratio f pe /f ce , and the fitted ratios f m /f cHe and f w /f cHe are mapped to other latitudes assuming a dipole magnetic field and constant cold plasma density along a magnetic field line. We also assume that f m and f w are constant as the waves propagate along the magnetic field lines from their source regions close to the magnetic equator. EMIC waves are thought to become more oblique as they propagate to higher latitudes from the equator ; however, Ni et al. (2015) investigated the sensitivity of electron diffusion to EMIC wave normal angle and found that the main differences where at eq < 40 • and E > 5 MeV. Therefore, following Kersten et al. (2014) to allow a direct comparison, we assume that the waves have a Gaussian distribution in X, X = tan , where is the wave normal angle, centered on = 0 with a width X w = tan 15 • . The cutoff in X is taken to be 2X w . For the ion composition we assume abundances of (H + , He + , O + ) = (0.94, 0.05, 0.01) .
In a cold plasma, EMIC waves of band i will only exist for f ci + 1 < f < f ci . When such a wave propagates to higher latitudes, f /f ci decreases and f may approach f ci + 1 from above, where i + 1 is the next ion species in the series {H + , He + , O + }. At which point the waves will become unguided and may reflect, refract or damp (Horne & Thorne, 1993 and so we assume they no longer take part in wave-particle interactions. Since the ion-ion hybrid frequency usually lies very close to the heavier ion cyclotron frequency we introduce the condition that f > f ci + 1 into the calculation. Similarly, at low latitudes, the high-frequency tail of the power spectrum may approach f ci from below. These waves will be damped by warm plasma effects due to wave-particle interactions with substantial populations of warm particles. As we are working in the cold plasma framework, we must not include wave power at frequencies where this theory breaks down, that is, when the refractive index tends to infinity at the waves as the wave frequency approaches the gyrofrequency. We therefore add the constraint that f /f ci < 0.97.
After summing over magnetic latitude to calculate a bounce-averaged diffusion coefficient for each observation, the results are then binned by activity and L * to construct average values. In this calculation, the individual bounce-averaged diffusion coefficients are summed within each bin and divided by the total number of observations, including those with no EMIC activity, within that bin, in order to obtain the average diffusion experienced at a given time within this bin. This results in bounce-and drift-averaged EMIC diffusion coefficients for the given activity and L * binning. For the L * bins we use a step size of 0.5L * .
In Figures 2d-2f we show the resulting pitch angle diffusion coefficients for our Kp parametrization at L * of 4, 4.5, and 5 during moderate activity, 2 ≤ Kp< 4. As expected, the diffusion is strongest at lower pitch angle and at higher energies but extends up to nearly 90 • equatorial pitch angle, eq . Stronger diffusion is found at L * = 5 as a result of more frequent EMIC wave activity. Electron diffusion by EMIC waves increases with Kp (not shown) but significant diffusion occurs even at relatively low Kp, Figures 2d-2f. This is in contrast to diffusion coefficients parametrized by P dyn , Figures 2a-2c, where a much stronger activity dependence is found and there is significantly lower diffusion at P dyn < 5 nPa, suggesting that electron scattering by EMIC waves is strongly influenced by P dyn .
In the lowest two panels of Figure 2, we compare the new Kp parametrized coefficients against previous results obtained by Kersten2014 and Ma2015. Note that we have reproduced Ma2015's calculations with

Effects of EMIC Waves on Electron Fluxes
In the previous section we calculated diffusion coefficients for EMIC wave-particle interactions using concurrent measurements. We will now assess their implications on radiation belt electron modeling. For this study, we perform 3-D simulations of the Earth's radiation belts using the British Antarctic Survey Radiation Belt Model (BAS-RBM) (Glauert et al., 2014a(Glauert et al., , 2014b, which is based on the phase-averaged Fokker-Plank equation that calculates the evolution of the electron phase space density. Primarily, BAS-RBM simulations are used for reproducing historical events and short-term forecasting. Currently, Kp is forecasted and used as a wave parametrization index when observational data is unavailable. In addition to Kp with bin infima of {0, 2, 4, 6} (henceforth Kp Ross), we present results using solar wind pressure, P dyn , as the EMIC driving index (henceforth P dyn Ross) with bin infima of {0, 2, 5} nPa. This choice is motivated by the observed link between solar wind pressure and EMIC activity (Anderson & Hamilton, 1993;Lessard et al., 2019;Saikin et al., 2016;Usanova et al., 2012). Additionally, Drozdov et al. (2017) found that a P dyn EMIC wave parametrization minimized the mean absolute error of 4.2 MeV electrons in radiation belt simulations. For comparison, we also present two otherwise identical simulations where the EMIC waves are from Kersten et al. (2014) and Ma et al. (2015), parametrized by Kp (henceforth Kp Kersten and Kp Ma).  electron fluxes are used for both the initial condition and the boundary conditions at L min , L max , and E min . The inner L * boundary is set to 2.5 and the outer L * boundary to 5.5. The minimum energy boundary is determined by assuming constant first adiabatic invariant with E = 150 keV at the outer L * boundary. Note that BAS-RBM assumes a dipole magnetic field configuration and so L * = L, although the wave measurements are converted to L * (section 2). The chorus diffusion coefficients are parametrized by AE (Horne et al., 2013), and the hiss diffusion coefficients are also parametrized by AE (Meredith et al., 2018). The chorus and hiss waves are spatially separated by an activity dependent mask as given in Meredith et al. (2018). For radial diffusion we use the Kp-driven magnetic component of the (Brautigam & Albert, 2000) formalism.
The location of the last closed drift shell is accounted for following Glauert et al. (2014a).
In Figure 3 we show the results of our simulations at 2.6 MeV. In panel (e) we compare the eq = 85 0 simulated flux at L * = 5 to Van Allen Probe data. P dyn Ross (purple), Kp Ross (orange), and Kp Ma (green) significantly improve the agreement between the modeled and observed (blue) 2.6 MeV equatorially mirroring electron fluxes compared to Kp Kersten (red). This is most apparent during quasi-steadily decaying times and during rapid loss. For instance, in mid-October 2015, a rapid order of magnitude decease occurs followed by ∼10 days of steady decay. P dyn Ross, Kp Ross, and Kp Ma reproduce the rapid drop and slow decay; however, the losses during both of these phases is insufficient in Kp Kersten, leading to approximately an order of magnitude discrepancy between modeled and observed fluxes by the 1 November. At lower L * , 10.1029/2020GL088976 panels (c) and (d), P dyn Ross and Kp Ross again perform better than Kp Kersten, with decay consistent with observations. Contrastingly, Kp Ma underestimates the electron fluxes at L * = 4. During electron acceleration phases, such as early October and early November, the flux in Kp Ross and Kp Ma are a factor of ∼2 less than is obtained in the P dyn Ross and Kp Kersten, leading to a greater discrepancy with observations. This is due to additional losses in this run coincident with the acceleration as a result of large Kp, as seen in panel (f). Note that Kp Ma has constant EMIC activity above Kp = 2. From panels (a) and (b), it is clear that the agreement at L * < 4 between the models and the observations is significantly worse that at L * ≳ 4.0.
In particular, relativistic electrons are transported into the slot region in the simulation, while this is not found in the observational data. This is likely a result of a combination of factors. First, our EMIC diffusion coefficients do not extend below L * = 3.5 due to lack of EMIC observations in the CRRES data set as a result of the magnetometer switched into low sensitivity mode in this region and EMIC waves were unable to be identified at these times. At L * < 4.0 the CRRES observational coverage is unlikely to be sufficient to calculate reliable diffusion coefficients in this way. Contrastingly, the statistical coverage of EMIC waves peaks at L * ∼ 5 where the discrepancy between modeled fluxes and observations is at its minimum . Second, rapid injection of electrons via radial transport is not always well captured by a simple radial diffusion parametrization. The most notable difference being a significantly larger discrepancy between observed fluxes and those in P dyn Ross. This again suggests that the observational coverage EMIC coverage of CRRES is insufficient, at L * ≲ 4, to fully reproduce relativistic electron decay.

Discussion
In this paper we have presented a new approach for calculating, using quasi-linear theory, average diffusion coefficients representing the diffusion of electrons by EMIC waves by averaging over observation-specific diffusion coefficients. In our EMIC diffusion model, we are assuming that electrons interact with EMIC waves sufficiently frequently on the diffusion timescale for an average representation to be appropriate. We have shown that 3-D Fokker-Planck radiation belt models using these new concurrent measurement EMIC diffusion coefficients considerably improve agreement between modeled and observed decay of relativistic electrons in the outer radiation belt at L * ≳ 4 when compared to models using EMIC diffusion coefficients calculated from average spectral and plasma properties, with the strongest agreement when the EMIC diffusion coefficients are parametrized by P dyn . In particular, relativistic near-equatorially mirroring electrons are scattered in pitch angle by the EMIC waves and are eventually precipitated into the atmosphere. This is in contrast to equivalent simulations using EMIC diffusion coefficients constructed using averaged spectral and plasma properties. It is necessary to capture the effects of the full range of wave-particle interactions that occur, not only those of an interaction under average conditions and with average spectral properties. Diffusion is possible over a wider domain in ( eq , E, L * ) space than is otherwise found using the simpler prescription. More precisely, diffusion extends nearly to 90 • and so electron losses into the atmosphere are possible over a much wider range of pitch angles than before. The differing behavior is a result of the wave-particle interactions in high-density regions and when the wave frequency is close to the EMIC wave bands upper gyrofrequency. At E = 4.2 MeV, the results of our Kp driven model are comparable to a Kp-driven model using an EMIC spectral profile where the central wave frequency is chosen to be close to the upper-bounding gyrofrequency following Ma et al. (2015). Contrastingly, at E = 2.6 MeV, our Kp-driven EMIC model gives closer agreement to observations at L * = 4.0. Our analysis of CRRES EMIC observations suggest that the central spectral frequency to helium gyrofrequency adopted in Ma et al. (2015) is unrepresentatively high and is only reached or exceeded in ∼1% of EMIC observations.
The L * range of the EMIC diffusion coefficients calculated here is limited due to the CRRES 16 Hz magnetometer data set having no identifiable EMIC observations at L * ≲ 3.5. Contrastingly, Qin et al. (2019) found EMIC waves at L * as low as ∼2.5 using the 64 Hz fluxgate magnetometer, EMFISIS, on board the Van Allen Probes. Additionally, owing to the brevity of the CRRES mission and the infrequency of EMIC waves, the statistical coverage of L * ≲ 4 is insufficient to create representative diffusion coefficients from the limited number of measurements. New Van Allen Probe observations of EMIC waves down to L * ∼ 2.5, and thus, the method shown here may help reduce the uncertainty in modeling the decay of the ultrarelativistic storage ring (Baker et al., 2013;Shprits et al., 2016;Usanova et al., 2014).
Given the improvement in the global simulations found when using our new diffusion coefficients, we propose that a similar approach should be explored for hiss and chorus waves as suggested in Watt et al. (2019). In this way plasma density and wave spectral variation will be included implicitly and thus enable a much better calculation of electron acceleration to relativistic energies which becomes very effective in low-density regions.