On the Variability of EMIC Waves and the Consequences for the Relativistic Electron Radiation Belt Population

The interactions between electromagnetic ion cyclotron (EMIC) waves and relativistic electrons are influential in diffusing radiation belt electrons into the loss code from which the electrons are lost into the atmosphere. These wave‐particle interactions between EMIC waves and electrons with energies of a few MeV or more depend strongly on wave spectra and plasma properties. Here we study the variability in wave spectra and plasma properties as a function of L* found during Van Allen Probe EMIC observations. These results are used to calculate statistical bounce and drift average diffusion coefficients that include the variation in wave spectra and plasma density as a function of L* and activity by averaging observation‐specific diffusion coefficients. The diffusion coefficients are included in global radiation belt simulations and the effect of the EMIC waves is explored. The distribution in the plasma frequency to electron gyrofrequency ratio decreases to lower values as L* decreases. As a result, few EMIC waves are able to resonate with 2–3 MeV electrons at L* ≤ 3.75 while electrons of the same energy at larger L* are diffused by EMIC waves in high density regions. In comparison, a sufficient number of EMIC waves are able to resonate with higher energy electrons, ≥4.2 $\ge 4.2$ MeV, at L* ≥ 3.25 to significantly affect the decay in electron flux. EMIC wave parametrisations of electron diffusion by EMIC waves are compared and solar wind dynamic pressure is found to give the best agreement with Van Allen Probe observations.

The effect of magnetospheric waves on radiation belt electrons is included in global radiation belt models via bounce and drift averaged diffusion coefficients. Recent studies have highlighted the importance of including the variability in plasma properties and wave spectra in radiation belt modeling, for EMIC waves (Ross et al., 2020) as well as whistler-mode waves (Watt et al., 2019(Watt et al., , 2021 and Ultra Low-Frequency wave (Thompson et al., 2020). Ross et al. (2020) developed a new approach for calculating electron diffusion coefficients by averaging observation specific diffusion coefficients (each calculated using the wave spectra and plasma density of each individual EMIC wave observation), rather than averaging the wave spectra and density measurements before calculating diffusion coefficients (e.g., Kersten et al. [2014]). This approach, applied to CRRES EMIC observations, significantly improved the agreement between global radiation belt models and Van Allen Probe observations of relativistic and ultra-relativistic electron flux compared to EMIC diffusion models using average wave spectra and plasma density. In particular, higher pitch-angled particles are diffused by the waves, removing a greater proportion of the electron population.
The CRRES EMIC observations, used in Ross et al. (2020), are limited to L* ≥ 4.0 due to the magnetometer switching to a lower sensitivity mode at lower L* (Meredith et al., 2003). However, EMIC waves have been observed by the Van Allen probes at lower L* (Gamayunov et al., 2020;Qin et al., 2019). Furthermore, the majority of EMIC observations during storms by the CRRES mission were in the storm main phase (Halford et al., 2010), while Saikin et al. (2016) reported more storm-time EMIC observations during the recovery phase during the early years of the Van Allen Probes mission. Additionally, Saikin et al. (2016) noted that EMIC observations observed by the Van Allen Probes have an activity distribution significantly different from those observed by CRRES. The CRRES mission coincided with a considerably more active solar cycle than the Van Allen Probes mission and therefore observed EMIC activity under different solar activity (Wang et al., 2016).
In this paper, we build on the work from Ross et al. (2020) to calculate statistical EMIC diffusion coefficients, by applying a similar approach which includes the variation in wave spectra and plasma properties, by using the more extensive Van Allen Probes data set allowing us to improve the statistics and extend our modeling to lower L*. The method presented in Ross et al. (2020) for calculating the diffusion coefficients is computationally expensive for extensive wave observation datasets due to the calculation of diffusion coefficients for each observation individually. In this paper, we present an approximate method that reduces the number of calculations necessary while agreeing well with the full calculation. A similar method may be used for modeling electron 3 of 25 diffusion by other wave modes in the radiation belts. These EMIC diffusion coefficients are then included in global radiation belt simulations to explore the effects of EMIC waves on ultrarelativistic electrons.
Our identification and analysis of the Van Allen Probes EMIC wave observations is presented in Section 2, with a focus toward the variation in wave spectra and plasma density and the implication for resonant energies. In Section 3 we calculate statistical bounce and drift averaged EMIC diffusion coefficients by averaging observation-specific diffusion coefficients parameterized by P dyn , Dst and Kp. In Section 4 we include these matrices into a global radiation belt model to asses the importance of EMIC waves and to compare the EMIC parameterizations. The results are discussed in Section 5 and our conclusions are presented in Section 6.

EMIC Wave Observations
In this study, we use data from the Van Allen Probe satellites that were launched on the 30 August 2012 into a near-equatorial (10° inclination) elliptical orbit with a perigee of 1.1R E and apogee of 5.8R E (Mauk et al., 2013). For the EMIC wave observations, we use the 64 Hz magnetometer data from the Electric and Magnetic Field Instrument Suite and Integrated Science (EMFISIS) between 11 September 2012 and 13 October 2019 (Kletzing et al., 2013). We do not use data from the Van Allen Probe B as we cannot reliably determine if the two sets of measurements are independent; the same EMIC waves observed by one instrument may also be observed by the other due to the temporal and spatial coherence of the waves (Blum et al., 2017). Additionally, the magnetometer data from Van Allen Probe A is generally less noisy than that from B, and therefore EMIC wave identification is more reliable (Wang et al., 2015).

Data Analysis and Extraction
The wave spectra are calculated using the magnetic field measurements rotated into a mean-field-aligned coordinate system where the coordinate transformation is performed based on the method from Jun et al. (2019), using a 32 s moving average for the background field. The mean-field-aligned orthonormal basis is defined by =̄∕|̄| , the y-component is in the direction ̄× and completes the right hand system via = × . We apply for a 3rd order Butterworth bandpass filter with a lower cut-off of 0.1 Hz and an upper cut-off of 10 Hz along each component (Wang et al., 2015). The power spectral density (PSD) of the waves is calculated using a fast Fourier transform applied to each magnetic field component with a Hanning time window of 64s with a stepsize of 16s. The wave properties are then calculated from the co-variance matrix in the frequency domain ( ) = ( ) * ( ) , where X i (f) are the Fourier transforms of the magnetic field components, following Means (1972). The wave properties are then averaged over 4 windows, giving a 64s measurement cadence.
At L* < 3, the proton gyrofrequency increases to ≥10 Hz which becomes comparable to the Nyquist frequency of the 64 Hz magnetometer on EMFISIS. We, therefore, limit our analysis to L* > 3, otherwise, we would be artificially restricting EMIC wave detections to frequencies significantly lower than the proton gyrofrequency at lower L* (Wang et al., 2015).
Wave features are extracted by the following algorithm based on Bortnik et al. (2007). In order to remove ambient noise, the rolling 3 min average PSD must exceed 10 times the daily average (Usanova et al., 2012), and therefore long lasting background features, likely resulting from the satellite itself, are excluded. For example, the narrow band wave power at ∼1.5 Hz (horizontal line in Figure 1a) does not satisfy this criterion and is therefore excluded. Occasionally, short-duration bursts of wave power occur which do not correspond to EMIC activity (Usanova et al., 2012), we exclude these by imposing the constraint that wave features must persist for a minimum of 5 min. . Panel (a) shows the power spectrum with the hydrogen, helium, and oxygen gyrofrequencies plotted in white. The power over the frequency-dependent threshold is shown in panel (b), where 2 events are identified and illustrated by the red and green shading. The average power spectra over the duration of the events, used in the events validation, are shown in the inset to panel (b). Panel (c) shows the result of the identification when the wave power in the helium band has been removed as it has not been classified as EMIC.
10.1029/2021JA029754 4 of 25 In order to capture most of the wave power of bursty events, shorter duration bursts in wave power that exceeds the threshold that is separated by less than 5 min from an existing feature are identified as the same feature (Wang et al., 2015), as illustrated in Figure 1b. As we require simultaneous density and spectral measurements for our diffusion coefficient calculations, we excluded observations when there is no f pe /f ce provided in the Level 4 EM-FISIS data  leaving ∼60% of the Van Allen Probe A data which we analyze.
After performing the feature extraction, we verify that the features are EMIC waves and remove broad band signals and magnetosonic waves (Anderson & Hamilton, 1993). For this, we use a combination of visual identification and supervised machine learning. First, we visually validate the features during 2013 and 2014 and label the features appropriately as EMIC or non-EMIC, for each band separately. For each feature, we then average the power spectral density over time on a frequency grid normalized by the proton gyrofrequency (see the subplot in Figure 1b). A Random Forest is then trained and tested on the average power spectral density profiles and validation labels from the visually inspected 2014 and 2013 observations, respectively. The prediction accuracy is ∼95%, where the differences in classification are low amplitude EMIC waves that are excluded by the Random Forest. The trained Random Forest is then applied to the remaining dataset of candidate features to identify the EMIC waves. An example of EMIC wave activity identified by the Random Forest is shown Figure 1c. The broadband wave power below the helium gyrofrequency that is picked out by our extraction algorithm is classified as non-EMIC by the Random Forest method and therefore rejected from the EMIC dataset.
To aid our analysis of the EMIC wave spectra and to only select the most intense harmonic within the band, we fit a Gaussian distribution to the power spectral density at each instance in time to each EMIC band. The corresponding upper and low cut-off's are determined by the highest and lowest frequency within the band that was identified as being above the threshold. Note that each EMIC wave feature will be made up of multiple line spectra at different instances in time, resulting in multiple fitted distributions for a given event, each representing 64 s of the EMIC wave, rather than averaging over the duration of the duration as is often done. For each of the spectral profiles, we retain the central Gaussian frequency, f m , the Gaussian width, and the upper and lower cutoffs. The wave intensity, 2 ( 2 ) , is then calculated by integrating the Gaussian distribution between the the cut-offs. Additionally, we retain f pe /f ce , MLT, the magnetic latitude and L* along with solar wind parameters and geomagnetic indices from the OMNI web. For consistency with our previous diffusion coefficient calculations, we calculate L* using OP77Q (W. P. Olson & Pfitzer, 1977).

Activity and Spatial Dependence
The distribution of EMIC waves depends on the solar and geomagnetic activity as well as the spatial location in the radiation belts. We, therefore, perform a statistical analysis of the EMIC wave activity which we can use to interpret later results. For brevity, we focus on the distributions in MLT and L* parameterized by solar wind pressure, P dyn , as similar statistical studies have been performed before, for example, see Chen et al. (2019), Jun et al. (2019), Saikin et al. (2016).
During the weak solar activity, P dyn ≤ 2 nPa shown in Figure 2 a and b, the occurrence rates of hydrogen and helium band EMIC waves peak at ≲1% just prior to noon. In comparison, the helium band waves have a broader distribution in MLT with slightly enhanced rates in the afternoon sector and predawn. The occurrence rate increases with increasing L* for both bands. The MLT and L* dependence of the average wave intensity, Figures 3a and 3b are well correlated with the occurrence EMIC rate.
At high solar wind activity, P dyn > 5 nPa shown in Figures 2e and 2f, the EMIC occurrence rates are significantly increased from the levels found during weak solar activity. Again the rates increase with L* and by L* ∼ 5, the occurrence rates reach 5%-10%, a factor of ∼10 increase compared to during the weak solar activity. There are enhanced hydrogen band EMIC occurrence rates at L* ≥ 5 over a broad range of MLT, from 2-19. Contrastingly, the wave intensity is peaked in the afternoon sector, Figure 3e. The hydrogen band waves outside of the afternoon sector are typically low amplitude, resulting in a lower average power spectral intensity. On the other hand, the helium band waves are largely found between 10 and 20 MLT, with few observations in the predawn sector. Enhanced EMIC wave activity on the dusk side is consistent with previous statistical studies Saikin et al., 2015) and self-consistent global models of EMIC wave excitation (Jordanova et al., 2008).
Before we calculate drift and bounce averaged diffusion coefficients, we need to understand the L* and activity dependence. We therefore average observations from each MLT and bin by L*, and P dyn , Kp, and Dst in turn.

Figures 4a and 4b
show that the MLT averaged occurrence rates increase with P dyn for both hydrogen and helium band EMIC waves. The strongest P dyn dependence is found at larger L*, with the occurrence rate reaching over 4% for hydrogen and over 3% for helium at L* = 5.25. Contrastingly, while P dyn < 2 nPa, the occurrence rates at all L* are less than 0.5% for each band. However, there are more EMIC wave observations during these quiet times than during high solar wind pressure, as shown in Figures 4c and 4d, despite the low occurrence rate as a result of longer satellite dwell times. The total duration of the EMIC wave observations increases with L* until L* = 5.25. The decrease at L* = 5.75 is due to lower dwell times and that the low-frequency limit of the Butterworth band pass filter approaches f cHe in the weaker magnetic field making identification of helium band EMIC waves more difficult. The average EMIC wave intensity for both bands, including times with EMIC activity and without EMIC waves identified that is, the average wave intensity over all time in the bin, increases significantly with P dyn even at low L* with 2 ( 2 ) > for helium bands waves reaching 10 −2 2 at L* = 3.75 and P dyn ≥ 5 nPa, compared to ∼10 −4 nT 2 when P dyn < 1 nPa. For our activity binning for our diffusion coefficient calculations, we want the binning to reflect the changes in the EMIC wave occurrence and wave intensity while retaining good statistical coverage. At L* = 4.75 and 5.25, the occurrence rate increases substantially between the 4 nPa ≤ P dyn < 5 nPa and P dyn ≥ 5 nPa bins, and therefore, for our diffusion coefficient calculations, we choose 5 nPa as the lower threshold for as our largest P dyn bin. The minimum total observation duration, with or without EMIC activity, when P dyn > 5 nPa within each L* and P dyn bin is 72 hr, although the duration of the observed EMIC activity at L* = 3.25 within this bin is low. For the lower bins, we adopt P dyn < 1 nPa, 1 nPa ≤ P dyn < 2 , and 2 nPa ≤ P dyn < 5 nPa, each of which has greater than 650 hr of observations. Occurrence rates for hydrogen and helium band electromagnetic ion cyclotron waves as a function of magnetic local time and L* for P dyn < 2 nPa, 2 nPa ≤ P dyn < 5 nPa and P dyn ≥ 5 nPa. The gray shading indicates that there are at least 2 hr of observations in a bin.
Figures 5a and 5b show that the occurrence rates of both hydrogen and helium band EMIC waves increase with Kp at all L* between Kp < 1 and 4 ≤ Kp < 5. The dependence is again strongest at larger L*, particularly for helium band EMIC, with occurrence rates reaching ∼4% at L* ≥ 4.75 during active conditions. The average wave power of both bands shows a strong dependence on Kp at all L* and therefore we take Kp < 2, 2 ≤ Kp < 4 and Kp > 4 as our bins for diffusion coefficient modeling. At Kp > 4, each L* bin has greater than 63 hr of observations while the lower Kp bins have at least ∼630 hr.
Finally, Figures 6a and 6b show that the highest occurrence rates are found when Dst decreases at all L*. At L* = 3.25, the majority of the EMIC wave observations are at Dst < −20 nT with few at higher values of Dst, Figures 6c and 6d. However, at L* ≥ 4.25, the dependence of EMIC occurrence rates of Dst is significantly weaker, particularly at larger L*. The average wave intensity also increases with decreasing Dst at all L*'s although the dependence is again stronger at lower L*. We note that there is a substantial difference between the strength of the geomagnetic storms with in the −200 nT ≤ Dst < −50 nT and −50 nT ≤ Dst < −20 nT bins, however the average wave intensity within these bins are very similar, and therefore we combine them. We discuss this further in the discussion. Our activity binning is then Dst < − 20 nT, −20 nT ≤ Dst < 5 nT, and Dst > 5 nT, with each bin having more than 350 hr of observations. . Average wave intensity for hydrogen and helium band electromagnetic ion cyclotron waves as a function of magnetic local time and L* for P dyn < 2 nPa, 2 nPa ≤ P dyn < 5 nPa and P dyn ≥ 5 nPa. The gray shading indicates that there are at least 2 hr of observations in a bin. 7 of 25

Variability in Spectral and Plasma Properties
In this section, we discuss the variability of the EMIC wave spectra and plasma density as these two properties strongly influence the wave-particle interactions. For this, we map the observations to the equator, as equatorial values are necessary for the diffusion coefficient calculations, and diffusion by EMIC waves is dominated by low latitude wave-particle interactions (see Section 3.2). For the mapping we assume a dipole magnetic field and that the wave power and plasma density remain constant along the field. In an average sense, close to the magnetic equator (|MLAT| < 20°) the plasma density deviates little from the equatorial value, for example, see Ozhogin et al. (2012) Figure 9.  (d) show the corresponding EMIC duration within the P dyn and L* bins. The average wave intensities are given in (e) and (f). The round brackets indicate strict inequality while the square brackets allow for equality for example, (0, 1) corresponds to 0 nPa ≤ P dyn < 1 nPa.
As mentioned above, the resonance between EMIC waves and electrons is sensitive to the wave spectrum and f pe / f ce . This can be seen in Figure 7, when the wave frequency approaches the upper bounding gyrofrequency of the band, the minimum resonant energy, E min , significantly decreases. Similarly, in high-density regions where f pe /f ce is large, E min is also decreased. Therefore, knowledge of how the wave spectra and f pe /f ce relate to each other and vary with L* is necessary for accurately modeling wave-particle interactions.
In order to represent the wave spectra, we can consider the Gaussian central frequency, f m , as this gives an indication of the frequency at which the EMIC wave power is centered. The distribution of ∕ remained fairly constant between L* = 3.2 and 5.7 for both hydrogen and helium band waves, Figures 8a and 8b, with little variation in average (crosses) and the 20% and 80% percentiles (black lines). There is a second small population of hydrogen band EMIC waves close to the proton gyrofrequency with ∕ > 0.9 found at L* ≥ 4.0 (Teng et al., 2019). At L* ≤ 4, the Butterworth upper cut-off frequency may be limiting the hydrogen distribution at high frequencies (red line), however, Teng et al. (2019) found that the number of high frequency EMIC waves decrease with deceasing L* and few were found at L* ≤ 3.5. The distribution of ∕ for hydrogen band EMIC waves is peaked at ∼0.42, consistent with the CRRES EMIC observations Meredith et al., 2014), with the 80% of observations with ∕ < 0.45 . The helium band waves have a broader distribution in f m relative to it's gyrofrequency with a significant number of observations close to the helium frequency. For instance, 20% of helium EMIC waves have f m above ∼0.75f cHe and ∼20% have f m less than ∼0.5f cHe , while the average value is ∼0.65f cHe , allowing helium band waves to resonate with a broad range of electron energies and pitch-angles. For both bands, the average ∕ , shown by the crosses in Figures 9a and 9b, increases with L* although the increase is much shallower for hydrogen band waves. Between the two wave bands, there is a clear distinction in distribution in ∕ . Helium band EMIC waves are found in predominantly higher density regions, such as a plasmasphere and plasmaspheric plumes, while hydrogen band EMIC waves are largely associated with lower density regions, such as the plasma trough, where the 80th percentile of hydrogen band densities is close to the 20th percentile of helium band densities at L* ≥ 4.0 (black lines) Horne & Thorne, 1994;Zhang et al., 2016). However, there is significant variation in ∕ at fixed L* with the highest obtained ∕ increasing with L*. As a result of the lower ∕ values found at low L*, these EMIC waves will typically resonate with higher energies than those at high L*, unless there is wave power very close to the upper bounding gyrofrequency, Figure 7. Note that while an electron drifts around the Earth, at fixed L*, it may interact with EMIC waves in a wide range in plasma densities as it changes MLT, this variation must be captured when calculating drift averaged diffusion coefficients.
We have over-plotted the L dependent density models from Sheeley et al. (2001) (red dotted line) and Ozhogin et al. (2012) (blue dotted line) in Figure 9, calculated assuming a dipole magnetic field and so L = L*. At L* ≥ 4.0, the L* dependent average ∕ found in the helium band EMIC observations exceeds both plasmaspheric density models and, by L* ∼ 5.0, 80% of the observations have higher densities than these models. At L* ≤ 3.6, the average ∕ for a helium band observation drops below the plasmaspheric density models and the majority of the observations are at lower densities. For hydrogen band EMIC waves, the L* dependent average ∕ lies between the plasmaspheric models and the Sheeley et al. (2001)  Finally, we show the distribution in ∕ as a function of ∕ in Figure 10. From this combination of parameters, assuming an ion composition of 94% hydrogen, 5% helium, and 1% oxygen  and that the EMIC waves are field-aligned, we can calculate the minimum resonant energies which are shown by the black contours over-plotted in Figure 10. Note that the wave power of the waves extends to higher and lower frequencies as the wave spectra have a nonzero width. There are a considerable number of both hydrogen and helium band waves that can resonate with ≲3 MeV electrons but significantly fewer that can resonate with ≲1 MeV electrons, and these are largely hydrogen band waves close the proton gyrofrequency, for instance, the high frequency Figure 7. Profiles of the minimum resonant energy assuming field-aligned waves for various f pe /f ce and ion composition of 94% hydrogen, 5% helium, and 1% oxygen. The red shaded regions give the prohibited frequency ranges used in the diffusion coefficient calculation, specified by f up . hydrogen population (Teng et al., 2019). Note that the cold plasma dispersion relation, used here, breaks down close to the gyrofrequency when the refractive index tends to infinity. As a result, the lowest resonant energies (<1 MeV) are likely to be underestimated by this approximation.

EMIC Diffusion Coefficients
In order to include electron diffusion by EMIC waves in global radiation belt models, we calculate statistical drift and bounce averaged diffusion coefficients, > , using the PADIE code . PADIE solves the cold plasma dispersion relation and uses quasi-linear theory to calculate the diffusion coefficients. In this work, we use a modified version of the PADIE code that takes an arbitrary wave power spectral density as an input rather than Gaussian inputs. We calculate statistical bounce averaged diffusion coefficients for both hydrogen and helium band EMIC waves by averaging observation-specific diffusion coefficients (Ross et al., 2020).

Calculation
The orbit of the Van Allen Probe satellites is inclined with respect to the geomagnetic equator, we, therefore, map each of the EMIC observations to the magnetic equator assuming a dipole, and that the number density and magnetic power spectra do not change with latitude along the magnetic field line. The average EMIC wave power observed by the Van Allen probes is approximately constant with latitude for absolute magnetic latitudes of ≤17°. We use the superscript eq to identify equatorial values and prevent ambiguity. We restrict our diffusion coefficient calculations to magnetic latitudes less than 20° assuming that the EMIC wave spectral profile does not vary with geomagnetic latitude. Summers et al. (2007) found that electron diffusion, by EMIC waves, at the loss cone is dominated by waves at |λ mag | < 20°: diffusion at larger equatorial pitch angles is restricted close to the magnetic equator. By neglecting wave power at higher magnetic latitudes we underestimate diffusion at high energies typically above 10 MeV. We discuss this further below in Section 3.2. For our nominal ion abundances, we use 94% hydrogen, 5% helium, and 1% oxygen following Kersten et al. (2014). For the local upper cut-off, f up , we take a value of 0.97 following Ross et al. (2020) to avoid warm plasma effects close to the upper bounding gyrofrequency  where the refractive index tends to infinity in the cold plasma approximation for oblique propagation. We discuss the sensitivity to this choice in the discussion section. Jun et al. (2021) performed a comprehensive analysis of EMIC wave normal angles using both Van Allen Probe and Arase data. They found that the majority of hydrogen and helium band EMIC wave modes have wave normal angles of <30 • . More oblique waves are detected near dawn at L > 8 Jun et al., 2021;Min et al., 2012). EMIC waves are thought to become more oblique as they propagate to higher latitudes from the equator Horne & Thorne, 1994). Gamayunov et al. (2018) noted a discrepancy between observations at low L and theory, they suggested that the superposition of multiple waves can lead to underestimation in the wavenormal angle when calculated using Fast Fourier Transforms. Ni et al. (2015) investigated the sensitivity of electron diffusion to EMIC wave normal angle and found that at L ≤ 5.5 the decay timescales were nearly independent of the wave normal angle model at E ≥ 2 MeV unless a very oblique wave is adopted over all latitudes. We performed calculations using the latitudinal varying wave normal angle model from Ni et al. (2015) and found minimal difference in the diffusion coefficients in the domain of interest to the fixed wave normal angle model used by Kersten et al. (2014) and Ross et al. (2020). Therefore 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 cut-off in X is taken to be 2X w . With this wave normal angle prescription, we are omitting the contribution from oblique waves near dawn in the outer magnetosphere, however, the average EMIC wave intensity is generally much lower in this region compared to the noon and the afternoon sectors (Figures 2-5). Ross et al. (2020) demonstrated a new method of calculating diffusion coefficients for electron diffusion by EMIC waves. In these calculations, observation-specific diffusion coefficients are calculated for each observation using the observed wave spectra and plasma frequency to electron gyrofrequency ratio. The drift and bounce averaged diffusion coefficients were then calculated by averaging the bounce averaged diffusion coefficients within an L* and activity bin. In this way, the full range of wave particle interactions resulting from the variability of the plasma environment was captured in the diffusion coefficients. However, their method is numerically expensive on large datasets as it requires the calculation of bounce averaged diffusion coefficients for each observation. Here we approximate the method from Ross et al. (2020) by binning the observations in ∕ and L* allowing multiple observations to be combined into one diffusion coefficient calculation while retaining the spectral information from each observation. The EMIC observations presented in the previous section are first binned by activity (which in this study is either P dyn , Kp, or Dst) and L* (adopting the bin centered L* value for all observations in the bin). Within each bin, we further bin the results by ∕ , again assigning the central bin value to ∕ . Given that we have point measurements, we make the assumption that the latitudinal variation of the wave is the same for each observation. As a result, the observations within a ∕ , activity and L* bin only differ through their equatorial spectra. For each observation, we reconstruct the EMIC power spectral density using the Gaussian parameters and cut-offs determined from the observations (see Section 2.1). We can add the power spectral densities of each measurement in a bin and calculate the bounce averaged diffusion coefficients using this summed profile. Note that this summed profile is rarely a Gaussian but has multiple peaks, coming from individual observations. The resulting diffusion coefficients are then a close approximation of the sum of the bounce averaged diffusion coefficients calculated from the observation specific diffusion coefficients within the ∕ , activity and L* bin. Finally, after summing the diffusion coefficients over the ∕ bins, while still retaining the L* and activity binning, and dividing by the number of observations within the L* and activity bin, including those with no EMIC observed activity, we have an approximation of the average observation specific diffusion coefficient. A more formal description is given Appendix A.
The error in the approximation is dependent on the granularity of the L* and x binning, with smaller bins giving better agreement and in the limit of infinitely small bin sizes is exact. After exploring different bin sizes, we adopt L* bins of width 0.5 L* and ∕ bin widths of 2.5, which give very good agreement with the full calculation ( Figure 11). Figure 11. Comparison of full method (average of observation specific diffusion coefficients) and reduced method (wave spectra of observations summed within ∕ bins before calculating the diffusion coefficients in each bin and summing over bins, normalized by the number of observations) for helium band electromagnetic ion cyclotron waves at |λ m | < 5° for P dyn > 5 nPa and L* = 4.75.

Results
Here we present our diffusion coefficients parameterized by P dyn , our results for the other parameterizations are somewhat similar and so are not shown here.
The bounce and drift averaged pitch-angle diffusion coefficients parameterized by P dyn are shown in Figure 12. Electron pitch-angle diffusion by EMIC waves dominates over the energy and cross term diffusion ) and therefore we focus on pitch-angle diffusion here. Diffusion increases with increasing solar wind pressure at all energies. When P dyn < 1 nPa at L* = 3.75, there is rapid drop off in pitch-angle diffusion below 5.6 MeV, Figure 12a, even at 10 MeV diffusion is limited to α eq ≤ 60° and therefore can not remove higher pitch-angled electrons. At high solar wind pressures, electrons with lower energies and higher α eq are diffused. The diffusion coefficients increase with increasing L* and extend to larger α eq . For instance, at L* = 3.75 when 2 nPa ≤ P dyn < 5 nPa (Figure 12b), diffusion of 1.78 MeV electrons (orange line) is limited to α eq ≤ 40° but at L* = 4.75 (Figure 12e) diffusion occurs up to near 90°, allowing for near equatorially mirroring particles to be removed from the belt.
The pitch-angle profiles of the diffusion coefficients look more irregular in profile than those calculated using average spectra which assume a single-peaked, Gaussian power distribution. This is because, in our model, there are contributions from many EMIC wave spectra, giving rise to more complex diffusion coefficient profiles.
In Section 3.1, we discussed the restricted latitude range used in this study and the rapid decrease in contribution with increasing latitude. Figure 13 shows that the diffusion coefficients decrease with latitude and the minimum resonant energy increases. The electron energies that resonate with the EMIC waves increase with increasing latitude as a result of the decreasing local wave frequency to proton gyrofrequency and plasma to electron gyrofrequency ratios, Figure 7. Additionally, the increased local pitch-angle of the particle at higher latitudes further increases the resonant energies.

Long Term Global Radiation Belt Simulations
In this section, we explore the effect of these new diffusion coefficients on the relativistic and ultra-relativistic populations using year long global radiation belt simulations. The interval between 1 March 2015 and 1 March 2016 is chosen as it exhibits a range of dynamics including multiple active periods where there is significant acceleration up to ultra-relativistic energies followed by quiescent periods with gradual decay, allowing us to study the effects of EMIC activity during a range of magnetospheric conditions.

Method
We perform 3D 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 solves the phase-averaged Fokker-Plank equation that calculates the evolution of the electron phase space density. The evolution of the phase space density, f, is given in terms of equatorial pitch-angle, α eq , Roederer L* and energy, E. (3) ( ) = (1.3802 − 0.3198(sin + sin 1∕2 )).
Where E 0 is the electron rest mass, τ c and τ M are the timescales of electron loss into the atmosphere and to the magnetopause, respectively. Here, μ = p 2 sin 2 α/(2m e B 0 ) and J = ∮ p ‖ ds are the first and second adiabatic invariants (Schulz & Lanzerotti, 1974), where m e is the electron rest mass, B 0 is the mean value of the geomagnetic field at the Earth's surface and p's the electron's momentum.
The fundamental equations solved by BAS-RBM are similar to other global radiation belt models such as VERB (Shprits et al., 2009), DREAM (Reeves et al., 2012), and SALAMBO (Beutier & Boscher, 1995). The models differ through their modeling of the wave-particle interactions by diffusion coefficients. Each of the global radiation belts models uses statistical models of the average wave and plasma properties to calculate the diffusion coefficients. The EMIC wave diffusion coefficients presented here are calculated in a significantly different manner where the variation in observed wave-spectra and plasma properties are included in the diffusion coefficient calculations. This is a necessary extension as the operation of calculating the diffusion coefficients and the process of averaging do not commute that is, the electron diffusion by an average wave is not equal to the average diffusion by waves (Watt et al., 2019(Watt et al., , 2021. By taking our approach a more representative effect of electrons on the radiation belts is included. In all of the runs below, we include hiss diffusion coefficients from Glauert et al. (2014a) which are parameterized AE*, defined as the maximum AE over the previous 3 hr, and restricted to inside the plasmapause. Chorus wave diffusion coefficients are parameterized by AE (Horne et al., 2013) and are restricted to outside the plasmapause. The plasmapause location is given by Carpenter and Anderson (1992). For radial diffusion, we use the Kp-driven magnetic component of the Brautigam and Albert (2000) formalism. The location of the last closed drift shell is calculated following the method of Glauert et al. (2014b) and effects included in the simulations.
Van Allen Probe A and B background corrected MagEIS electron flux (Claudepierre et al., 2015) and REPT electron flux (Baker et al., 2012) 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 be inside the slot region at 2.5, and the outer L* boundary to 5.3 which is the largest L* consistently sampled throughout the simulation period. The minimum energy boundary is determined by assuming constant first adiabatic invariant with E = 150 keV at the outer L* boundary in order to exclude lower electron energies that are influenced by physics not accounted for in the model. Note that BAS-RBM assumes a dipole magnetic field configuration internally as is usual for all global radiation belt models.

Contributions to Relativistic Electron Flux Decay
Both plasmaspheric hiss waves Pinto et al., 2019;Thorne et al., 2013) and EMIC waves (Drozdov et al., 2017;Kersten et al., 2014;Ma et al., 2015;Shprits et al., 2016;Usanova et al., 2014) are believed to be important for the decay of relativistic and ultrarelativistic electrons. Additionally, Wang and Shprits (2019) showed that chorus waves at high latitude can also contribute to the loss of 0.9 MeV MeV electrons. In order to assess the relative importance of these waves as a function of L* and electron energy, we compare three simulations with the following combinations of EMIC and plasmaspheric hiss diffusion coefficients: EMIC and plasmaspheric hiss; EMIC and no plasmaspheric hiss; plasmaspheric hiss and no EMIC. Note that radial diffusion and chorus waves are included in each of these simulations. For these runs, we use our P dyn parameterized EMIC wave diffusion coefficients. We focus on the near equatorially mirroring electron flux for the comparison as the decay of this flux is dependent on electron diffusion at lower pitch-angles and therefore requires an accurate model of diffusion at all pitch-angles. For BAS-RBM we take the α eq = 85° and for the Van Allen Probe observations the central pitch-angle bin. Figure 14 shows the results for 4.2 MeV electron flux. For L* ≥ 4.25, losses due to EMIC and hiss waves (red curve) and EMIC waves alone (green) are much higher than losses due to hiss alone (blue). This indicates that hiss is not very effective in removing electrons at these energies. At L* = 3.75 losses due to hiss and EMIC (red) remain very effective whereas losses due to hiss alone (blue) or EMIC waves alone (green) are far less effective. This indicates that the combination of the two waves modes is very important for removing electrons at large pitch angles in this region, as suggested earlier by the width of the pitch angle distribution in Figure 12. Hiss provides diffusion at pitch-angle close to 90° while EMIC waves diffuse the particles at lower pitch-angles. At L* = 3.25 initially, the best agreement between the data and the model is for hiss and EMIC (red) which give the largest losses for the whole period compared to those wave modes on their own. Between 1 March and the middle of August the combined hiss and EMIC model does better than the other two models, but after that, all three models show several rapid increases in flux whereas there is a large reduction and drop out in the data in early September followed by a more gradual increase that is not reproduced. It is interesting to note that the data show a significant increase in flux at L* = 3.75 compared to a reduction at 3.25. Since EMIC and hiss waves contribute mainly to electron loss and not acceleration, the increase in flux in the model is due to inward radial diffusion as shown in the discussion. This suggests that the radial diffusion coefficients used in this model, and in many other global models, needs to be substantially improved for low L*. Note that the divergence of the model from the observation does not correspond to the largest Kp when D LL is at its largest. We return to this in the discussion. Figure 15 shows the results for 2.6 MeV. At L* = 4.75 the model results for EMIC and hiss (red) and EMIC alone (green) are very similar and reproduce the flux variations quite well indicating that EMIC waves are the most important in this region. In contrast, hiss alone (blue) is able to reproduce some but not all of the decay. A similar result is obtained at L* = 4.25. However at lower L* ≤ 3.75 the EMIC and hiss model performs the best and the model for EMIC waves alone (green) does not capture the decay. The results illustrate that EMIC waves acting alone are not very effective at causing loss at energies of a few MeV, rather the combined action of hiss waves and EMIC waves is necessary to explain the observed decay.

Activity Parameterization
As described in Section 3, we have calculated 3 separate EMIC diffusion coefficients, parameterized by P dyn (as used in the previous section), Kp and Dst. Here we compare these parameterizations and their agreement with Van Allen Probe observations to determine which parameterizations work the best. In the previous section, we found that the impact of EMIC waves is easily identifiable at 4.2 MeV and L* ≥ 3.25, in which case we focus our analysis of the EMIC driving parameter on this energy channel. Figure 16 shows the 4.2 MeV electron flux with EMIC waves parameterized by P dyn (red), Kp (blue), and Dst (green). At L* ≥ 4.25, all three simulations perform well and there is little difference between the parameterizations, although the simulation Kp driven EMIC waves typically under estimates the peak electron flux by a greater amount than the other simulations. At lower L* ≤ 3.75 (Figures 16a and 16b) the Kp and P dyn parameterizations perform significantly better than the Dst parameterization, with the Dst parameterization systematically overestimating the flux. The P dyn and Kp driven models at L* = 3.75 agree fairly well with the observations. Again, at L* = 3.25, there is good agreement between the P dyn and Kp model runs compared to the observations until August 2015, with decay rates consistent with those observed.

Metrics
In order to asses the EMIC models, we make use of the metrics from Morley et al. (2018) and Glauert et al. (2018), namely, the median symmetric accuracy (MSA) and the symmetric signed bias (SSPB). Both of these metrics are straight-forward to interpret, symmetric when penalizing under and over prediction, and robust to outliers. The median symmetric accuracy is given by with Q i = X i /Y i with Y i denoting the model value and X i the observation. It is informative to note that 50% of results are with a factor exp(Median {|ln |}) = 1 + MSA∕100 of the data, referred to as the median error factor. The SSPB gives the percentage overestimate or underestimate by the given percentage of the median error and is defined as The metrics are applied to near equatorially mirroring electron flux, namely, the differential α eq = 85° electron flux modeled by BAS-RBM is compared to the 90° flux measurement observed by the Van Allen Probes.
Figures 17a and 17b show that the median error factors and SSPBs for 2.6 MeV electrons at L* ≤ 3.75 are smallest when hiss waves are included (blue crosses) compared to when they are omitted (orange dots). In contrast, including EMIC waves has little effect on the metrics (blue crosses compared to green dots). At larger L*, both EMIC waves and hiss waves are important with the smallest errors when both waves are included (blue crosses compared to orange and green dots). For 4.2 MeV electrons, Figures 17c and 17d, including EMIC waves significantly improve the metrics at L* ≤ 3.75 (blue crosses compared to green dots), although hiss waves are also very important for reducing the errors and biases (blue crosses compared to orange dots). At L* ≥ 4.25, omitting hiss waves (orange dots) has minimal effect on the metrics and achieves similar scores to an identical simulation where they are included (blue crosses), while if EMIC waves are left out (green dots) the errors and bias are much greater. The importance of both EMIC waves and hiss waves for the decay of multi-MeV electron fluxes is consistent with the conclusions of Drozdov et al. (2020) based on long-term simulations, with the EMIC waves providing rapid losses close to the loss cone and hiss waves contributing to electron diffusion at larger pitch-angles.
At L* ≥ 4.25, the simulations with EMIC waves parameterized with P dyn (blue crosses), Kp (red crosses), and Dst (purple crosses) perform comparably well, although Kp has a significant negative SSPB for 4.2 MeV electrons at L* = 4.25 which is related to the lower flux shown by the blue line in Figure 16d. At L* ≤ 3.75, the P dyn and Kp simulations perform better than the Dst simulation with substantially lower median error factors and SSPBs.  Figure 9 clearly shows EMIC waves are organized by plasma density, with hydrogen band waves in low density regions such as the plasma trough while helium band waves are in high-density regions such as the plasmasphere and plasma plumes Zhang et al., 2016). However, we do not parameterize the waves by density, instead, we parameterize by solar or geomagnetic activity, and the spread in plasma environment is captured in the model, due to the way we calculate our drift-averaged DC from observations specific diffusion coefficients which sample from all MLTs, which is not captured in previous average statistical models.

Discussion
Results from numerical diffusion experiments (Watt et al., 2021) indicate that averaging observation-specific diffusion coefficients may be appropriate if the wave activity varies sufficiently rapidly, even given a large amount of variability in observations of wave intensity and plasma conditions and no further parameterization. The parameterizations investigated here remove some sources of variability and can be used to construct more accurate models. Further investigations into the temporal and spatial variability of diffusion coefficients for EMIC waves using the observations shown here and results from for example, Blum et al. (2017) would also confirm whether the averaging strategy used here is most accurate. There is no denying that the creation of a deterministic model of averaged diffusion coefficients is much more tractable and appealing than using multiple stochastic parameterizations (Watt et al., 2021) to capture all variability in the wave-particle interactions important to Earth's radiation belts.
EMIC waves that approach their upper bounding gyrofrequency are affected by warm plasma. Cold plasma theory, as used in this study, breaks down in this limit as the refractive index tends to infinity when the wave frequency approaches the gyrofrequency. To avoid invalidating the approximation, we have set the upper cut-off to be f up = 0.97f ci . However, the exact threshold at which cold plasma theory breaks down depends on the plasma properties. We, therefore, test the sensitivity of our results to the imposed threshold and perform an identical set of diffusion coefficient calculations parameterized by P dyn with f up = 0.95f ci . There is a negligible difference between the two model results (not shown) giving us confidence that results robust to this assumption. Shprits et al. (2013) considered the energy range of electrons that interact with whistler mode waves compared to EMIC waves during a long lasting third radiation belt in September 2012. They concluded that at L ∼ 3.3 whistler mode waves were able to resonate with electrons at relativistic energies (<2 MeV) and but were unable to do so with ultrarelativistic electrons (>2 MeV) near the equator. Instead, they suggested that EMIC waves are important for the decay of the third belt at ultrarelativistic energies (Shprits et al., 2016(Shprits et al., , 2018. Our results here are consistent with Shprits et al. (2013) at L ∼ 3.3, but also show that the electron energy range that are significantly affected by EMIC waves depends on L* as a result of the increasing f pe /f ce with L*, allowing for substantial electron diffusion at ≲2 MeV and L* ≥ 4.0.
The global radiation belt simulations with EMIC waves parameterized by Kp underestimate peak electron flux more than when Dst and P dyn are used. This may be due to overestimating EMIC losses outside of the plasmasphere, removing the electrons that are being accelerated by chorus waves and radial diffusion. Kp is typically enhanced for longer than P dyn while Dst recovers more quickly. Therefore EMIC waves in the Kp driven simulation are able to counteract electron acceleration by chorus waves and radial diffusion for longer, suppressing the build-up of relativistic electrons. Alternatively, as the chorus diffusion model does not include variability in plasma density, other than through variation in the MLT, acceleration up to MeV energies may be under represented .
At L* ≤ 3.75, the Dst simulation performs significantly worse than the other 2 simulations, while the Kp simulation performs the best (Figure 17). The greater overestimation of the flux in the Dst simulation is likely partly due to the larger reservoir of flux at higher L*, compared to the Kp and P dyn simulations, that is diffused inward. In Figure 6 we showed that Dst is a weaker indicator of EMIC wave activity, and hence electron diffusion, than both Kp and P dyn , particularly at L* ≥ 4.25. The weaker dependence of EMIC wave activity on Dst may be related to the inability of the Dst index to capture high-speed-stream-driven geomagnetic storms (Borovsky & Denton, 2010;Borovsky & Shprits, 2017) which can drive enhanced EMIC activity (Gamayunov et al., 2020) and are associated with a daylong increase in solar wind pressure at storm onset. Furthermore, Borovsky (2017) found that the correlation between radiation belt flux and Kp is much stronger than the correlation with Dst suggesting that the Dst index is insufficient to characterize the radiation belts.
We obtain good agreement with observations in our long-term simulations with EMIC waves parameterized by P dyn . Drozdov et al. (2017) performed long-term Fokker-Planck simulations of the radiation belts using the VERB code comparing EMIC waves parameterized by Kp, Dst, and AE indices, solar wind velocity, and solar wind pressure. In their study, they used a representative EMIC spectral profile and an average plasma density for their EMIC diffusion coefficient calculations. They minimized the absolute mean error by optimizing the wave amplitude and activity threshold for the presence of EMIC waves. They found that solar wind dynamic pressure provides the best parameterization of EMIC waves. Our Van Allen Probe EMIC data set shows significantly increased EMIC activity during magnetospheric compression (Figure 4), consistent with previous statistical results (Anderson & Hamilton, 1993;J. V. Olson & Lee, 1983;Saikin et al., 2016;Usanova et al., 2012). Furthermore, the average wave intensity increases with solar wind pressure, even at L* < 4, which in turn leads to faster rates of electron diffusion. EMIC waves at low L*s have previously been linked with magnetospheric compressions; Qin et al. (2019) found EMIC waves at 1.9 < L* < 3.2 during extreme magnetospheric compression on 22 June 2015. However, it is unlikely that solar wind pressure can generally directly lead to EMIC waves at L* < 4.0 as moderate magnetospheric compression will only affect the outer magnetosphere. Rather, it is likely that solar wind pressure is correlated with geomagnetic conditions which are favorable for EMIC wave activity at lower L* In Figures 14-16 we showed that the modeled >2.6 MeV electron flux at L* ≥ 3.75 agree well with the Van Allen Probe observations. However, at L* = 3.25 (e.g., Figure 16a) we obtain sudden increases in >2.6 MeV electron flux during August 2015 that are not observed in the Van Allen Probe data. Instead, the observations show a gradual increase in flux on a day to week timescale. Contrastingly, observations of lower energy flux (<1 MeV) do show a sudden increase in flux (not shown) which agrees well with the modeled flux. This suggests that we are missing or poorly capturing the physical processes governing the MeV population during this period. In Figure 18a we show results (blue line) from a simulation with radial diffusion reduced by a factor of 100 from the 1 August 2015. The flux no longer shows the sudden increases found previously (red line) implying these increases are a result of transport by radial diffusion. Similarly, later in the simulation with radial diffusion, there are additional increases in the modeled flux L* = 3.25 which are also not observed. Two possibilities explain these differences. The radial diffusion model may break down at several MeV energies in this region (Lejosne et al., 2013) or losses by local wave-particle interactions during this period are sufficiently strong and localized to L* ∼ 3.25 to counter act the flux increase by radial diffusion. The ratio of f pe /f ce is generally low at L* ∼ 3.25 and exceeds 12 only 20% of the time at (Figure 9c) and, therefore, unless there is substantial wave power very close to the gyrofrequency, the 2.6 MeV electrons will not be diffused by EMIC waves at these low densities ( Figure 10). Furthermore, there is no evidence for strong EMIC waves close to an upper bounding gyrofrequency during these times. If the differences were due to hiss waves then the increase in flux at lower energies would also be reduced. It is therefore more likely the discrepancy is a result of radial diffusion rather than diffusion by local wave-particle interactions.
The pitch-angle and energy of the electrons that can resonate with an EMIC wave depend on the cold ion composition. Measurements of the cold ion composition are particularly difficult due to satellite charging (Olsen et al., 1985) preventing reliable direct measurements. Warm ion measurements can be made but they may differ significantly from the cold ion population that dominates the mass. Therefore, we are unable to use the ion composition observed during each observation to complement the observation-specific wave spectra and plasma density. In this study we have adopted the ion composition from Kersten et al. (2014) however a range of fixed values have previously been adopted (Ma et al., 2015;Meredith et al., 2003;. A higher He + percentage leads to electron diffusion at lower energies for helium EMIC waves but the converse for hydrogen EMIC waves. By the very nature of a statistical model, we are averaging over many events within an activity and L* bin. In a statistical approach parameterized by a geomagnetic activity index or solar wind measurement, we are assuming that the index captures the variation and captures the average response of the electron population under those conditions. The true state of the system however may not always be well represented by a statistical model. For the case of EMIC waves, the growth of the EMIC waves depends on the ion temperature anisotropy and electron plasma density. Similarly, the wave-particle interactions depend on electron plasma density and ion composition. In turn, the ion population and electron plasma density are affected by the time history of the magnetosphere through electric fields and disturbances in the magnetic field. Under different magnetospheric driving, these populations will evolve differently leading to different EMIC wave behavior (Gamayunov et al., 2020) which will not necessarily be well captured in a statistical model of this form. Coupled ring current-plasmaspheric models (Jordanova et al., 2008) provide a means to include these effects more consistently however it is challenging to constrain such models due to the number of degrees of freedom.

Conclusions
We have analyzed the L* dependence of the wave spectra and the corresponding plasma frequency to electron gyrofrequency ratio of EMIC waves observed by the Van Allen Probe A satellite. These observations have then been used to calculate bounce and drift averaged EMIC diffusion coefficients from averaging observation specific diffusion coefficients. When calculated in this way, the variability in wave spectra and plasma density is included in the diffusion coefficients. These have been included in a global radiation belt model, BAS-RBM, and year long simulations compared against Van Allen Probe observations to determine at what energies and L*s EMIC waves are important and which activity parametrization gives the best agreement. Our principle results are as follows: 1. For L* > 4 the average plasma frequency to electron gyrofrequency ratio for helium band EMIC waves is much higher than statistical plasmaspheric models, while hydrogen band waves are typically in low density regions. 2. EMIC waves do not contribute significantly to the decay of the ultrarelativistic storage ring at 2.6 MeV due to the low values of f pe /f ce in the region, instead diffusion by hiss waves is the dominant loss process. In comparison, EMIC wave and hiss waves are both necessary for the decay of the 4.2 MeV component of the storage ring. 3. At L* ≥ 4.25, electron flux decay is largely controlled by EMIC waves at ≥4.2 MeV, but a combination of hiss waves and EMIC waves determine the electron losses at 2.6 MeV.

Appendix A: Reduced Diffusion Coefficient Method
For observation i, let the bounce averaged diffusion coefficients be denoted by where we have dropped the α, E assumed dependence and we have defined = ∕ for brevity. D is the bounce average diffusion operator assuming that ion composition is constant between observations and that the latitudinal variation of the wave normal angle is the same between observations. Now consider a specific activity bin and L* bin, and label the set of point in this bin as  and  , respectively. The average of the observation specific diffusion coefficients within that bin is then given by (A2) where the summation includes all observations, even those without observed EMIC waves. Note that calculating the numerator involves calculating bounce averaged diffusion coefficients for each observation. For large datasets, this becomes computationally heavy and therefore we make a series of approximations to reduce the number of calculations that are necessary.
First, we adopt the central L* bin value, ̄ * , for observations in a L* bin. The bounce average diffusion coefficient then only depends on x i and ( ) 2 . Furthermore, if the observations are now binned by x, into bins labeled  , and the bin center values, , are adopted for all x i in each  , then the bounce averaged diffusion coefficients within that density bin only depend on ( ) 2 . We can now use the fact that operator D is directly proportional to ( ) 2 to rewrite the numerator of Equation A2 as (A3)