ULF Wave Driven Radial Diffusion During Geomagnetic Storms: A Statistical Analysis of Van Allen Probes Observations

The impact of radial diffusion in storm time radiation belt dynamics is well‐debated. In this study we quantify the changes and variability in radial diffusion coefficients during geomagnetic storms. A statistical analysis of Van Allen Probes data (2012–2019) is conducted to obtain measurements of the magnetic and electric power spectral densities for Ultra Low Frequency (ULF) waves, and corresponding radial diffusion coefficients. The results show global wave power enhancements occur during the storm main phase, and continue into the recovery phase. Local time asymmetries show sources of wave power are both external solar wind driving and internal sources from coupling with ring current ions and substorms. Wave power enhancements are also observed at low L values (L < 4). The accessibility of wave power to low L is attributed to a depression of the Alfvén continuum. The increased wave power drives enhancements in both the magnetic and electric field diffusion coefficients by more than an order of magnitude. Significant variability in diffusion coefficients is observed, with values ranging over several orders of magnitude. A comparison to the Kp parameterized empirical model of Ozeke et al. (2014) is conducted and indicates important differences during storm times. Although the electric field diffusion coefficient is relatively well described by the empirical model, the magnetic field diffusion coefficient is approximately ∼10 times larger than predicted. We discuss how differences could be attributed to data set limitations and assumptions. Alternative storm‐time radial diffusion coefficients are provided as a function of L* and storm phase.

anomalies in space due to electrostatic charging and discharging, and the occurrence of these anomalies increases with geomagnetic activity (e.g., Choi et al., 2011;Love et al., 2000). Therefore, it is a key outstanding goal of the Space Weather community to fully understand the dynamics and variability of the radiation belt population during geomagnetic storms (Fry, 2012;Schrijver et al., 2015).
However, the source and loss processes of radiation belt electrons is complex and multifaceted (Reeves, et al., 2003). Radially localized peaks in the distribution of the radiation belts can arise due to local wave-particle interactions Claudepierre et al., 2013;Horne & Thorne, 1998;Horne, 2007;Horne, Thorne, Shprits, et al., 2005;Reeves, et al., 2013;Thorne et al., 2013) and flux dropouts at the outer boundary Turner et al., 2012). In the presence of strong radial phase space density gradients, electrons can then be subjected to significant radial redistribution that acts to smooth the localized peaks. This process is radial diffusion. Radial diffusion of electrons can act to both accelerate electrons as they are transported earthwards and decelerate electrons as they are transported radially outwards, due to the conservation of the first adiabatic invariant. Consequently, radial diffusion has been attributed to playing an important and sometimes key role in radiation belt energization (Mathie & Mann, 2000;Rostoker et al., 1998;Ukhorskiy et al., 2005), rapid electron flux dropouts (e.g., Shprits et al., 2006), and the "slot-penetrating" supply to the inner radiation belt (e.g., Loto'aniu et al., 2006;Zhao & Li, 2013). Understanding the dynamics of radial diffusion during geomagnetic storms is of considerable importance.
In a stationary magnetic field, radiation belt electron motion forms closed paths around the Earth, termed drift shells. However, when these electrons are subjected to random electromagnetic fluctuations occurring on timescales comparable to their drift period, violation of the third adiabatic invariant occurs, and the electron is scattered radially. Radial diffusion describes the average rate at which the trapped population is radially scattered to new drift shells due to a large number of small electromagnetic perturbations (Fälthammar, 1965;Kellogg, 1959;Parker, 1960). Previous work has established that broadband ULF waves are a key driver of radial diffusion in radiation belt electrons, owing to wave frequencies (∼1-10 mHz) that are comparable to the electron drift frequency (e.g., Brito et al., 2020;Elkington et al., 2003).

Diffusion Coefficients
Significant efforts have quantified the magnitude of potential radial diffusion due to ULF wave power into radial diffusion coefficients. In particular, Fei et al. (2006) assumed that the diffusion coefficient for a dipole field, D LL can be separated into magnetic and electric field components,  Lejosne (2019) the assumption made by Fei et al. (2006) that the magnetic and electric fields are independent and uncoupled is incorrect, and using the approach of Fei et al. (2006) can lead to an underestimation of the total radial diffusion coefficients by a factor of 2. We note that this discrepancy is minor compared to the large variability of values observed, which is demonstrated in Section 3.2. Although we opt to use the Fei et al. (2006) formalism in this analysis, the associated assumptions should be kept in mind when interpreting.
Using the Fei et al. (2006) where B E is the equatorial magnetic field strength at the surface of the Earth, R E is the radius of the Earth, L is the L-shell, M is the first adiabatic invariant, q is the electron charge, γ is the relativistic SANDHU ET AL. respectively. The field perturbations have an azimuthal wave number and wave frequency, ω, that satisfies the drift resonance condition (ω − mω d = 0) with an electron with bounce-averaged drift frequency ω d . Note that the magnetic field perturbations considered correspond to the compressional component (parallel to the background magnetic field) and the electric field perturbations to the azimuthal component (perpendicular to the background magnetic field and in the east-west direction).
As demonstrated by Ozeke et al. (2014), these descriptions of the diffusion coefficients can be simplified by assuming that the wave power at a given m value can be represented as a fraction of the total observed wave power, giving where P B (L, f) and P E (L, f) are the power of field perturbations, which is dependent on the frequency, f = 2π/ω d of the perturbations as well as the L value. The angled brackets in Equations 3 and 4 indicate that the average power is taken over all frequencies. The analytical expressions provided in Equations 3 and 4 allow for the estimation of diffusion coefficients directly from measurements of the power spectral density of magnetic and electric field perturbations, at a given time and location. Note that this approach assumes that the power spectral density value is representative of the average power spectral density along an electron drift path. Several studies, such as Ozeke et al. (2014) and Brautigam and Albert (2000), have developed empirical models parameterizing the diffusion coefficients by the Kp index. It should be noted that the parameterization is of the ULF wave power itself, since this is the only observed quantity remaining in Equations 3 and 4. The empirical models provide an easily accessible method to estimate the magnitude of radial diffusion in the inner magnetosphere for a given time, requiring only knowledge of the Kp index and without relying on suitable spacecraft coverage. The results of studies such as Ozeke et al. (2014) and Brautigam and Albert (2000) have informed our understanding of how radial diffusion varies with L value and with geomagnetic activity. Ozeke et al. (2014) show that the electric field component of radial diffusion is consistently greater than the magnetic field component (  E B LL LL D D ) for 1 ≤ L < 7. In contrast, Brautigam and Albert (2000) observe that, although the electric field dominates at low L values, the magnetic field diffusion coefficient is dominant for approximately L > 4. In both the Ozeke et al. (2014) and Brautigam and Albert (2000)  The enhancement in radial diffusion coefficients during geomagnetic storms is further supported by case studies, where event-specific radial diffusion coefficients are estimated from observations of the magnetic and electric field perturbations. Case studies such as Li et al. (2014), Olifer, et al. (2019, Pokhotelov et al. (2016) and Schiller et al. (2017) Ozeke et al. (2014). The magnitude of mis-estimation varied throughout the event, and at times the difference between empirically modeled values and event-specific diffusion coefficients was multiple orders of magnitude. In this study, we aim to identify how and why radial diffusion varies during geomagnetic storms through a statistical analysis of ULF wave fluctuations observed by the Van Allen Probes. The results determine the typical variations during geomagnetic storms and establish whether case studies of single storms are characteristic of typical storm time variations. SANDHU ET AL.

Data and Method
In situ observations of perturbations in the local magnetic and electric field were provided by the Van Allen Probes (Mauk et al., 2013). The Van Allen Probes consisted of two identically instrumented spacecraft (Probe A and Probe B) that commenced observations of the inner magnetosphere in October 2012, with operations ceasing in July 2019 for Probe B and October 2019 for Probe A. The spacecraft were in a 9 h orbit, with an apogee of 5.8 R E geocentric radial distance, and an inclination of 10°. The orbital perigee was ∼600 km during the prime phase, and began a period of perigee lowering maneuvers in January 2019. The orbital apogee precessed in local time, spanning all local times in less than 2 years. Due to the repeated sampling of the inner magnetosphere over a range of L and Magnetic Local Time (MLT) values for almost 7 years the Van Allen Probes present a highly suitable data set for statistical analyses of the outer radiation belt.
This study employs magnetic field observations from the Electric and Magnetic Field Instrument Suite and Integrated Science (EMFISIS) instrument (Kletzing et al., 2013) and electric field observations from the Electric Field and Waves (EFW) instrument (Wygant et al., 2013). The full data set includes observations from both Probe A and Probe B. We use Probe A data from January 01, 2013 to January 01, 2016, and Probe B from January 01, 2013 to June 16, 2017, where these time limits exclude degraded EFW data quality toward the end of the Van Allen Probes lifetimes. In the following subsection we detail how the magnetic and electric field observations are used to estimate radial diffusion coefficients. Figure 1 shows an illustrative example of how the magnetic and electric radial diffusion coefficients are estimated for a typical geomagnetic storm that occurred during March 2013. Figure 1a shows the Sym-H index (Iyemori, 1990) trace during the event that exhibits the characteristic storm time variations. Following the approach of Walach and Grocott (2019), the Sym-H index was used to identify the event, as well as identifying the start and end times of the storm initial, main and recovery phases. The storm phase timings are marked by the vertical dashed lines in Figure 1. Further details about the storm identification are provided in Section 2.2. The L and MLT location of Van Allen Probe A is shown in Figures 1b and 1c.

Estimating Radial Diffusion Coefficients
The magnetic and electric field measurements were then taken for the full event, with a time resolution of 4 s. The spin axis electric field was suitably computed from the EFW observations under the assumption that E ⋅ B = 0. In order to extract perturbations, the background magnetic and electric field was first identified by taking a running average over a 20 min sliding window. The background field was subtracted and the residual field observations were transformed to a magnetic field-aligned coordinate system. The coordinate system was defined as having the parallel component aligned with the background magnetic field unit vector, the unit vector in the azimuthal direction is defined as oriented eastwards and perpendicular to the geocentric position vector, and the unit vector in the poloidal direction completes the coordinate system. Following the coordinate transformation, the parallel component of the magnetic field perturbation, B ∥ , and the azimuthal component of the electric field perturbation, E ϕ , was selected.
The perturbations were then used to estimate the associated power spectral densities using a Morlet wavelet transform. The power spectral density was limited to a frequency range of 1-15 mHz and an L range of SANDHU ET AL.  In chronological order, the dashed lines show the start of the initial phase, the start of the main phase, the start of the recovery phase, and the end of the recovery phase. 3-7.5. Lower L values were excluded as field perturbations cannot be reliably separated from the rapidly varying background field observed by the spacecraft at perigee. Figures 1d and 1e shows the total power spectral density (power summed over frequencies from 1 to 15 mHz) for the magnetic field and the electric field, where values have been binned and averaged for L value and time. The L value bin width is 0.5 and the time bin width is 3 h. For this storm, clear variations in power spectral density during the storm are apparent, with enhancements across all L values during the storm main phase. We note the presence of data gaps in the electric field power spectral densities. During periods of spacecraft charging and Earth eclipse events, electric field data is considered unreliable and omitted. This results in lower sampling of the electric field perturbations throughout the study, but we note that more than sufficient sampling is obtained for reliable statistical analysis.
The power spectral densities of the magnetic field, P B , and electric field, P E , were substituted into Equations 3 and 4, respectively, and provide estimates of the radial diffusion coefficients B LL D and E LL D . The estimation assumes that the observed power spectral densities are constant over all MLT values along the drift path, and the accuracy of this assumption is explored later. Note that throughout the manuscript, we explore the dependences and variations of the total power spectral density in the magnetic and electric field, P B and P E , to understand dynamics across the ULF wave band (e.g., Figures 1d and 1e). However, it is emphasized that the frequency-averaged wave power is used to calculate the corresponding radial diffusion coefficients, as mandated by Equations 3 and 4. Figures 1f and 1g shows the cluded. This database is statistically analyzed to explore how the power spectral density and radial diffusion coefficients vary on average during storm times, the variability across storms, and the physical drivers of the changes.
It is noted that although we have opted to use the method of Ozeke et al. (2014) to convert power spectral densities to diffusion coefficient values, other formalisms exist (see Lejosne and Kollmann (2020) for a thorough review). As mentioned previously, important differences can exist between the different approaches and it is highlighted that the approach adopted by Ozeke et al. (2014) (which follows work by Fei et al. (2006)) can lead to an underestimation of the total radial diffusion coefficients by a factor of 2 (Lejosne, 2019). Although it is shown in the subsequent analysis that this discrepancy is comparatively minor relative to the large variability in the observed values. Furthermore, a key aim of this study is to assess the reliability of the widely used empirical model of Ozeke et al. (2014), particularly focusing on the effect of parameterization by the Kp index. Employing a consistent method of calculation will aid in the comparison, and allow the impact of other factors (e.g., parameterization and datasets) to be isolated.

Identifying Geomagnetic Storms
In order to obtain the storm time database of radial diffusion coefficients, geomagnetic storms were identified for the Van Allen Probes data period. The storms were identified from timeseries of the Sym-H index using an automated algorithm developed and fully detailed by Walach and Grocott (2019). The Sym-H index was employed instead of the traditional Dst index (Sugiura & Poros, 1964) as it can be considered a high time resolution version of the Dst index (Wanliss & Showalter, 2006) and hence will provide higher accuracy of storm phase timings. Figure 1a shows a typical example of a storm identification using the Sym-H index. A geomagnetic storm can generally be separated into three distinct phases: The initial phase, the main phase, and the recovery phase. The start and end of these storm phases are indicated by the dashed vertical lines in Figure 1a. The initial phase is evident from a positive perturbation in the Sym-H index due SANDHU ET AL.

10.1029/2020JA029024
5 of 20 to enhancements in the magnetopause currents and typically lasts ∼1 day (Walach & Grocott, 2019). The inital phase is followed by the main phase, which is characterized by a rapid negative change in the Sym-H index due to significant and dramatic enhancements in the ring current population (e.g., Sandhu et al., 2021). The main phase has a typical duration of a few hours (Walach & Grocott, 2019). The final storm phase is the recovery phase, where the Sym-H index gradually recovers and increases to quiet time levels as the ring current decays (e.g., Sandhu et al., 2021). The recovery phase typically lasts several days (Walach & Grocott, 2019).
The Walach and Grocott (2019) algorithm first identifies a storm from any period where the Sym-H index crosses below a storm time threshold of −80 nT. The start of the recovery phase is marked as the time of the Sym-H index minimum. The start of the main phase and the end of the recovery phase are marked as the times immediately prior to and after the storm peak, where the Sym-H index is less than a quiet time level of −15 nT. The initial phase is defined to include the time of the Sym-H index maximum, and the start of the initial phase is marked as the time immediately prior to the Sym-H index maximum where the Sym-H index is at the quiet time level of −15 nT. Note that the storm time threshold (−80 nT) and the quiet time threshold (−15 nT) follow definitions by Hutchinson et al. (2011). The use of a −80 nT storm time threshold will naturally exclude any weaker storms that have a Sym-H index minima greater than this, and so the storms considered in this analysis will represent relatively large storm events. The −80 nT threshold was chosen to reliably avoid misidentification of other processes (such as large substorms, reconnection events, and magnetopause oscillations) that can manifest as weak depressions in the Sym-H index (Hutchinson et al., 2011). Therefore, the following analysis will reliably represent ULF wave activity and radial diffusion during strong geomagnetic storms. Future work will aim to explore the dependence of radial diffusion on storm intensity and storm drivers.
Overall, the algorithm identifies 45 storms during the Van Allen Probes data set. Over all these storms, the final data set comprises 8 × 10 6 samples of the magnetic field diffusion coefficient, and 5 × 10 6 samples of the electric field diffusion coefficient. This corresponds to total sampling time periods of 350 and 250 days, respectively.

Power Spectral Density
First we explore spatial dependences in the ULF wave power, providing insight into the sources of the wave power and the key physical drivers of the radial diffusion coefficients. We show that important local time dependences exist, implying significant consequences for estimating event-specific radial diffusion coefficients from single spacecraft measurements. sectors. During the main phase, Figure 2b demonstrates that the P B values are enhanced compared to the initial phase, with typical values of ∼10 2 nT 2 Hz −1 . The enhancement is observed across nearly all spatial bins. The enhancement is smallest in the pre-noon MLT sector, imparting a weak local time asymmetry in the spatial distribution of P B . Following the main phase, Figure 2c shows that the median P B values subside to ∼10° nT 2 Hz −1 in the recovery phase. Values peak in the afternoon sector during both the main and recovery phases.
The power spectral density of the electric field also presents noteworthy storm-time variations. Figure 2d shows the median values of P E in the initial phase are ∼10 −1 mV 2 m −2 Hz −1 . No strong local time asymmetries are present in the initial phase. In contrast to the initial phase, Figure 2e shows that P E is significantly enhanced during the storm main phase with values exceeding ∼10 1 mV 2 m −2 Hz −1 . Values are strongly enhanced across all MLT sectors, with an L dependence observed such that P E tends to increase with L. Figure 2f shows that the values subside to P E ∼10° mV 2 m −2 Hz −1 in the recovery phase. A local time asymmetry is apparent with values peaking in the morning sector (06 < MLT < 12; Figure 2e).
Overall, Figure 2 shows that strong local time asymmetries arise during geomagnetic storms. Interestingly, the asymmetries and storm-time variations differ between the magnetic and electric fields. The physical processes that shape and contribute to the observed spatial dependences of P B and P E will be discussed in Section 4.
Although Figure 2 shows that there is important spatial structure in the power spectral density, the effect on a radiation belt electron is dependent on the power spectral density encountered by the electron on its drift path. In order to estimate the P B and P E values along an electron drift path, we bin data for L*, where L* is the third adiabatic invariant (Roederer, 1970;Roederer & Lejosne, 2018). The L* value is calculated using the International Radiation Belt Environment Modeling (IRBEM) code (https://sourceforge.net/projects/ irbem/) with the Tsyganenko and Sitnov (2005) magnetic field model for an equatorially trapped particle. It is noted that multiple methods for calculating L* values exist, and this is a source of additional variability to the statistical analysis (Thompson et al., 2020). It is also stressed that in a distorted and non-dipolar magnetic field configuration, the L value and L* values are not directly comparable. Asymmetries in P B and P E observed along the electron drift path are a combination of local time asymmetries displayed in Figure 2 and the asymmetry of the drift path.  Figure 3a shows that during the main phase the power spectral density is quasi-symmetric in local time across all L* values, with typical values of ∼10 −1 nT 2 Hz −1 . Values are slightly higher for dayside MLT sectors, but there is less than an order of magnitude variation across MLT sectors. During the main phase, the P B values are enhanced across all MLT sectors. A local time asymmetry is prominent, with values peaking above 10 2 nT 2 Hz −1 in the post-noon and nightside MLT sectors (Figure 3b). In the recovery phase, Figure 3c shows that the local time asymmetry is subdued and P B has reduced to ∼10° nT 2 Hz −1 . ). For the outermost L* bin (5 ≥ L* < 6), the values are lowest in the morning sector. In the recovery phase, the observed P E is reduced to approximately initial phase values.
In Figure 3 the outermost L* values are not included as there are limited observations for L* > 6, particularly during the main phase. During storm time conditions, magnetopause shadowing reduces the number of closed drift paths (e.g., Staples et al., 2020) and the last closed drift shell is often within L* = 7 (Olifer, et al., 2018). Consequently, there are often no particles with closed drift paths in 6 ≤ L*<7 during geomagnetic storms. Across the observable L* values, Figure 3 shows that there is limited variation with L* and the profiles are typically within interquartile ranges.
An important consequence of the local time asymmetries is the spatial variations they impart to calculated diffusion coefficients. As radial diffusion is a drift averaged process, the radial diffusion coefficients should SANDHU ET AL.

10.1029/2020JA029024
7 of 20 describe an average over all local times. However, Figures 2 and 3 show that during the main and recovery phases there are strong local time asymmetries observed in power spectral density and hence calculated radial diffusion coefficients. This implies that observations of the ULF waves that control radial diffusion during storm times are highly dependent on the local time position sampled. Figures 2 and 3 show that using spacecraft (single or multiple) or ground based magnetometers that sample a limited MLT sector for a single event can lead to significant under-or over-estimates of radial diffusion, as the approach would neglect spatial asymmetries. For example, event-specific radial diffusion coefficients estimated by Schiller et al. (2017) used Van Allen Probe observations in the pre-dawn sector. If ULF wave activity was significantly enhanced in, for example, the dusk sector, the estimated event-specific radial diffusion coefficients would be unrealistically low.

Radial Diffusion Coefficients
Section 3.1 demonstrates that the power spectral density of both the magnetic field and electric field varies along a drift path. As radial diffusion is a drift-averaged process, the following results bin data by L* to describe the drift path and include contributions from all MLT sectors. Statistical estimates of radial diffusion coefficients are obtained by averaging over all MLT sectors for given L* values. D for each individual sample at a given time. The data have been further binned for storm phase, where the blue profiles correspond to the initial phase, the pink profiles correspond to the main phase and the orange profiles correspond to the recovery phase. Due to the comparatively longer duration of the recovery phase compared to the initial and main phases, the number of samples in the recovery phase is higher. The median value for each distribution is indicated by the position of the open colored circles above each panel, using the same color coding as the distributions. In order to reduce the spatial variability of the values, only samples within an L* range of 4 ≤ L* < 5 are shown here, and it is noted that the features are consistent for other L* ranges. Figure 4 highlights some key features of the diffusion coefficients. First, the electric field diffusion coefficients are typically larger than the magnetic field counterparts (comparing Figures 4a and 4b), with the SANDHU ET AL.   ).
Although Figure 4 shows key features of the storm time diffusion coefficients, information regarding the temporal variability during a storm is lost by binning solely by storm phase. Alternatively, the storm time variations can be further investigated through a normalized superposed epoch analysis, and the results of this analysis are shown in Figure 5. As before, values within an L* range of 4 ≤ L* < 5 are shown to reduce variability due to spatial dependences and focus on storm time variations, although the variations are representative of the other L* values. In Figure 5, a sample is binned according to its sample time relative to the duration of the current phase. For a given storm phase, the duration is scaled to the average duration of that phase across all storm events. These average durations are 37 h for the initial phase, 11 h for the main phase, and 61 h for the recovery phase. This method organizes the temporal variations based on the progression during a storm, and has been demonstrated to be an effective approach in organizing storm time behavior throughout the inner magnetosphere (e.g., Hutchinson et al., 2011;Murphy, et al., 2018;Wharton et al., 2020;Yokoyama & Kamide, 1997). Furthermore, this approach benefits by assigning each storm phase equal importance, as opposed to favoring a particular epoch in the non-normalized superposed epoch approach.  demonstrated to be dominant solar wind parameters in controlling ULF wave power (Bentley et al., 2018), and thus are expected to be important factors in shaping the storm time diffusion coefficients. Figure 5b shows that during the initial phase the solar wind speed is typically 350-400 km s −1 . During the main phase the solar wind speed is progressively enhanced, peaking in the early recovery phase at v SW ∼550 km s −1 . Throughout the remainder of the recovery phase, v SW remains elevated compared to the initial phase with median values of ∼400-500 km s −1 , and large variability in values as indicated by the broad interquartile ranges. Similarly to v SW , Figure 5c shows that the IMF B Z component is at low levels during the initial phase with values close to 0 nT. During the main phase, Figure 5c shows that B Z experiences a dramatic southward excursion and reaches values of ∼−8 nT by the late main phase. In the recovery phase the B Z values return to lower magnitudes with a slight southward bias. The interquartile ranges show significant variability throughout the early recovery phase.
Magnetic field diffusion coefficient,  As Figures 4 and 5 are focused on the L* range of 4 ≤ L* < 5, Figure 6 investigates storm time changes across multiple L* values. Figure 6a shows the median B LL D as a function of L*, where the bars show the interquartile range of the samples within the bin. The blue, pink, and orange profiles correspond to samples in the initial phase, the main phase, and the recovery phase, respectively. Figure 6b shows the E LL D samples in the same format as Figure 6a. Figure 6 shows that both B LL D and E LL D increase with L*, across all storm phases, with an approximate L * 10 dependence. This is a consequence of both the diffusion coefficient formalism, where Equations 3 and 4 are shown to include dependences on L, as well as the slight radial gradients in P B and P E observed in Figures 2 and 3. As highlighted in the presentation of Figures 3 and 6 also shows that the higher L* bins (L* ≥ 5.5) are often unsampled during the main phase of the storm. This is a direct consequence of the Earthwards displacement of the last closed drift shell during geomagnetic storms, as previously discussed.
Finally, Figure 6c shows the L* profile for the sum of the magnetic and electric field diffusion coefficients, in the same format as Figures 6a and 6b. As expected, the profiles are similar to the electric field diffusion coefficients shown in Figure 6b, which dominate over the magnetic field contribution. The values shown in Figure 6c provide an estimate of the average diffusion coefficients specifically for storm time conditions, which may prove useful for the modeling and understanding of radial diffusion of typical geomagnetic storms.

A Comparison to Empirical Models
Widely used empirical models, such as Ozeke et al. (2014), are parameterized by the Kp index and show that during geomagnetically active conditions (high Kp values), both the magnetic and electric field diffusion coefficients are increased compared to quiet time conditions (low Kp values). In order to establish whether the Kp parameterized models are representative of typical storm time conditions, we directly compare the observed values to empirically modeled values using the Ozeke et al. (2014) Figures 7a and 7c shows an L* dependence such that the difference is larger at low L* values and smaller at high L* values. D D line during the initial and recovery phase, although the ratio has considerably smaller values of ∼1.4. It can also be seen that the values are much closer to 1 and lie slightly below the line at low L* bins in the main phase (Figure 7e). Taking into account to variability of the data as shown by the interquartile ranges, the difference between observed and modeled values is not deemed to be significant for the electric field component.
Overall, the results show that although a Kp parameterized model does provide realistic storm time electric field diffusion coefficients, there is a tendency for the magnetic field diffusion coefficients to be underestimated by a factor of ∼10. The implications of this are discussed further in Section 4. It is also highlighted that the discrepancy between the observed and modeled values is comparable to the large degree of variability in values, as shown by Figure 4a and the interquartile ranges in Figure 6a.

Why do Radial Diffusion Coefficients Increase During Storms?
The analysis presented in this manuscript demonstrates key features of the storm time radial diffusion coefficients. We find that both the magnetic and electric field diffusion coefficients exhibit strong and significant enhancements during geomagnetic storms across all L* values observed. The enhancements are observed to SANDHU ET AL.  begin in the late initial phase of the storm, continue through the entirety of the main phase, and only begin to subside in the early recovery phase. The changes are driven by enhancements in the respective power spectral density, which are observed to have strong local time asymmetries (Figures 2 and 3). Here we will discuss potential and likely drivers of the power spectral density enhancements, and thus the enhancements in the radial diffusion coefficients.
It is well established that a key and dominant source of broadband ULF wave power in the inner magnetosphere is directly linked to external solar wind driving. Kelvin-Helmholtz instabilities at the magnetopause flanks generate ULF wave power locally (Chen & Hasegawa, 1974a, 1974bMann et al., 1999;Rae et al., 2005;Southwood, 1974) and buffeting of the magnetopause by solar wind pressure fluctuations enhances ULF wave power across the dayside (Kepko et al., 2002). Previous work supports these sources of ULF wave power by showing strong correlations between storm time wave power and solar wind pressure, IMF B Z , and coupling functions (Koskinen & Tanskanen, 2002;Posch et al., 2003;Takahashi et al., 2012). When only single solar wind parameters are considered, solar wind speed, IMF B Z , and variation in number density are the dominant parameters that account for a large degree of the variability in ULF wave power for all conditions (Bentley et al., 2018). For active conditions, the elevated solar wind speed associated with both CMEs and CIRs is observed to correlate with enhanced magnetospheric ULF wave power Simms et al., 2010;Zong et al., 2009). Due to the characteristics of the external wave power sources, local time asymmetries are imparted to the occurrence and amplitude of ULF wave activity such that they maximize in the dawn-noon MLT sector (Liu et al., 2009;Nosé et al., 1995;Nykyri, 2013;Pahud et al., 2009;Takahashi et al., 2016). The analysis presented in this manuscript shows that during the main phase and into the early recovery phase, where the solar wind driving displays dramatic enhancements (Figures 5b and 5c), corresponding enhancements in wave power are observed in the dawn-noon MLT sector and across the dayside (Figure 2). The wave power is elevated for both the magnetic and electric field components, and contributes to the corresponding enhancements in the diffusion coefficients.
However, Figure 2 also demonstrates that enhancements in ULF wave power occur globally during the main phase and recovery phase compared to the initial phase. Large enhancements are observed in the afternoon sector for the magnetic field fluctuations, and across the nightside sector for the electric field fluctuations. ULF wave power in these sectors is weakly linked to solar wind driving (Takahashi et al., 2012), and the sources of this wave power can instead be attributed to internal sources of ULF wave activity. Coupling of ULF waves and westward drifting ring current ions can drive enhancements in narrowband ULF wave activity predominantly in the afternoon and dusk MLT sectors (Baddeley et al., 2005;Engebretson & Cahill, 1981;Hughes, 1983;James et al., 2013James et al., , 2016Nosé et al., 1998;Woch et al., 1990). Referring to Figure 2, afternoon sector enhancements during the main and recovery phase relative to the initial phase are apparent and are particularly prominent for the magnetic field component. The results suggest that ULF wave power excited by ring current ions contribute to the storm time radial diffusion coefficients, and could be dominant for the magnetic field component.
An additional and noteworthy internal source of ULF wave power is substorm activity. Studies, such as Nosé et al. (1998), show that substorms generate azimuthal ULF fluctuations across the nightside magnetosphere and peaks at 01-02 MLT. During the main phase of geomagnetic storms, substorm size and occurrence is markedly elevated and this enhancement in the source is reflected by main phase enhancements in ULF wave power for both the magnetic and electric field component across the nightside (Figures 2b  and 2e). During the recovery phase the level of substorm activity subsides, and the postmidnight and broader nightside enhancements are reduced in magnitude.
Overall, storm time enhancements in ULF wave power occur across all MLT sectors due to a combination of both external and internal sources. These internal and external sources also exhibit important differences in the resulting azimuthal wave number, m. Studies show that externally driven waves that dominate the dayside MLT sectors typically have low m numbers (Rae et al., 2005;Tu et al., 2012), whereas ULF waves on the nightside have markedly higher m numbers (Barani et al., 2019;Fenrich et al., 1995;Sarris & Li, 2017). The high m ULF waves can originate from substorm-related activity, where they can be driven by drift-bounce resonances with injected ions (James et al., 2013;Walker et al., 1982;Yeoman et al, 2010). The presence of these high m number waves has implications for our estimations of radial diffusion coefficients. underestimates of the diffusion coefficients of more than an order of magnitude (Barani et al., 2019;Tu et al., 2012). This assumption is a focus for future work, where we will more accurately account for variations in the m number and thus refine the estimated storm time diffusion coefficients.
Another key feature of the storm time ULF wave power enhancements is associated with the radial distribution. Figures 2 and 3 show that the enhancements extend across all L and L* values observed, and the wave power penetrates to low radial distances. During quiet and non-storm time conditions, external solar wind driving is a dominant source of ULF wave power and the ULF wave amplitudes are expected to decay exponentially from the source (i.e., the dayside magnetopause) (e.g., Southwood, 1974). Furthermore, wave power propagating from the dayside magnetopause can be effectively reflected at the plasmapause boundary (Lee et al., 2002). In contrast with quiet time conditions and in agreement with the observations presented here, previous studies show that during storm times enhanced wave power is observed at very low radial distances and within L values of ∼4 (Georgiou et al., 2018;Lee et al., 2007;Loto'aniu et al., 2006;Rae et al., 2019). As highlighted by these studies, the accessibility of wave power to low L values can be enabled by changes in the Alfvén continuum. Studies, in particular Wharton et al. (2020), show that Eigen frequencies across a large range of L values are significantly depressed during the main phase of the storm. This is attributed to a weakening of the local magnetic field strength due to the magnetic field perturbations associated with an enhanced ring current population (Jorgensen et al., 2004;Sandhu et al., 2018), as well as an increase in the local plasma mass density due to enhanced concentrations of heavy ions (e.g., Kale et al., 2009;Kistler & Mouikis, 2016;Sandhu et al., 2017;Yue et al., 2019). Furthermore, studies such as Kale et al. (2009), Katus et al. (2015 show that during the main phase of the storm the plasmapause boundary moves to lower L values. This Earthward displacement can exceed several L values for extreme cases. The formation of plasmaspheric plumes in the afternoon sector during the main and recovery phase of the storm can also "trap" fast mode waves within the plumes (Degeling et al., 2018), which may explain the afternoon sector enhancement observed in Figures 2b and 2c. Overall, the suppressed Eigen frequency continuum in conjunction with the Earthwards displacement of the plasmapause allows ULF wave power to access lower L values during the main phase, as evidenced by the enhanced wave power at low L values (Figure 2). Power enhancements at low L values can also be enabled by the compression of the magnetopause, acting to move the external ULF wave power source Earthwards. Murphy et al. (2015) shows that for both CME and CIR driven storms, the compression of the magnetopause during the main phase independently contributes to significant enhancements in power at low L values.
The results demonstrate that both internal and external sources of wave power contribute significantly to enhanced radial diffusion coefficients during geomagnetic storms. It is also suggested that large scale changes in the shape and internal conditions of the magnetosphere result in enhanced diffusion coefficients at low radial distances.

How Reliable is the Kp Parameterized Model During Storms?
A key aim of the study was to identify whether a simple Kp parameterization, often used in empirical models such as Brautigam and Albert (2000) and Ozeke et al. (2014), provide an accurate description of radial diffusion coefficients during geomagnetic storms. As established by Figure 7, the magnetic field diffusion coefficients are systematically higher by ∼10 times compared to the empirically modeled values. The electric field diffusion coefficients are also higher, but by a smaller multiplication factor of ∼2. This offset between observed and empirically modeled values is consistent across all L* values and during all storm phases. It is also interesting to note the large spread of samples in Figure 7, which highlights that the discrepancy between observed and empirically modeled values can be highly variable for single storm samples.
Differences from the empirical model are supported by storm time case studies of radial diffusion. For example, Jaynes et al. (2018) and Olifer et al. (2019) demonstrated that during the main phase of the March 2015 storm the observed diffusion coefficients varied significantly compared to values using the Ozeke et al. (2014) model. The magnetic field diffusion coefficients were larger than the model values, whereas the electric field diffusion coefficients were lower than the model values. This is roughly consistent with the trends shown in Figures 7b and7e. In contrast, a statistical study by Murphy et al. (2015) showed that observed main phase electric field diffusion coefficients were typically larger by more than an order of SANDHU ET AL.

10.1029/2020JA029024
14 of 20 magnitude compared to the values predicted by the Ozeke et al. (2014) model, and the overestimation was most prominent at low L values. However, Murphy et al. (2015) also note the large degree of variability observed in the diffusion coefficients, similarly to that observed in this study. Discrepancies with empirical models are further supported by Murphy et al. (2016), where it was shown that empirical Kp parameterized models can under/over-estimate electric field diffusion coefficients by four orders of magnitude similarly to the results shown in Figures 7d-7f.
The driver of discrepancies between observed and model values may arise through the measurement and analysis techniques. First, we consider differences in the magnetic field component. The Ozeke et al. (2014) empirical model for the magnetic field diffusion coefficient is based on observations from the Geostationary Operational Environmental Satellite (GOES) and the Time History of Events and Macroscale Interactions during Substorms (THEMIS) spacecraft. These observations cover L values from 5-7. Coverage of the non-geosynchronous L values is provided by the THEMIS observations, which are restricted to a time period of 2007-2011. The resulting empirical model will not account for the large enhancements in power at low L values, and the sampled L values are likely to lie outside of the last closed drift shell during storm times. It can be argued that the region assessed in this study lies outside the coverage of the empirical model leading to discrepancies in the magnetic field diffusion coefficient. In addition, the temporal coverage of the THEMIS spacecraft mostly corresponds to the solar minimum phase, whereas the Van Allen Probes measurements cover the solar maximum and declining phases. Studies such as Murphy et al. (2011) show that ULF wave power is enhanced during the solar maximum and declining phases, suggesting that the diffusion coefficients measured in this study will be larger than diffusion coefficients measured in the solar minimum phase. This supports the underestimation of the magnetic field diffusion coefficient by the Ozeke et al. (2014) empirical model (Figures 7a-7c). Finally, a key difference between the Ozeke et al. (2014) model and the results shown here are the choice of drift shell coordinate. The Ozeke et al. (2014) model averages over values binned for L, whereas we choose to average over L* as it is a more realistic description of electron drift trajectories. The difference between L and L* can be significant at large radial distances (approximately L > 5) (Roederer, 1967). However, as the radial dependence of power is small compared to, for example, local time asymmetry and event variability, the impact of using L as opposed to L* is not expected to be the dominant source of differences.
The Ozeke et al. (2014) description of the electric field diffusion coefficients is based on ground magnetometer measurements from stations in the Canadian Array for Realtime Investigations of Magnetic Activity (CARISMA) (Mann et al., 2008) and U.K. Sub-Auroral Magnetometer Network (SAMNET) (Yeoman et al., 1990) arrays. The data is restricted to samples in the dayside sector (06-18 MLT) and covers L values from 2 to 8 (Ozeke et al., 2012). A dipole magnetic field is assumed for the mapping from station magnetic latitude to L value. Although the L coverage is well suited for a comparison with the Van Allen Probes observations, the limitation of the Ozeke et al. (2014) to dayside MLT sectors and the assumption that the observed power is independent of local time may be important. The empirical model would not account for any power enhancements in the electric field in the nightside sector, which can be important in the main and recovery phase (Figures 2e and 2f). This would result in the model underestimating the average power over a drift orbit, and consequently underestimating the electric field diffusion coefficient. The assumption of a dipole magnetic field configuration may also contribute to discrepancies. During the main phase of a geomagnetic storm, the enhanced ring current significantly distorts and inflates the magnetic field configuration in the inner magnetosphere and field lines with a given footprint latitude can map to larger L values than during non-storm periods (Ganushkina et al., 2010;Sandhu et al., 2018). The electric field power on the dayside tends to increase with L so the magnetometer measurements would be overestimating values. This is reflected in the main phase comparison shown in Figure 7e, where observed values tend to be slightly lower than the empirically modeled values for low L* bins. As with the magnetic field component, the L value is also used for drift-averaging instead of the L* values, but as only the dayside sector is considered by Ozeke et al. (2014) the impact of using the L value is not considered to be significant.
A key feature of the Ozeke et al. (2014) model is the choice of parameterization. As shown by Figure 5, the diffusion coefficients exhibit important variations within the storm. In particular, the average values peak in the late main phase. This is a feature that cannot be extracted from Kp parameterizations, which averages together samples from different storm phases, as well as averaging over both storm time and quiet time sam-SANDHU ET AL.

10.1029/2020JA029024
15 of 20 ples. Studies have shown that the storm time ULF wave power is significantly larger than non-storm time values for the same Kp value (Dimitrakoudis et al., 2015), resulting in observed values being larger than the modeled values (Figure 7). These discrepancies support efforts and suggestions to consider parameterization by solar wind driving and storm epoch to more reliably describe radial diffusion coefficients (Bentley et al., 2019;Dimitrakoudis et al., 2015;Lejosne, 2020). Furthermore, probabilistic modeling approaches to space weather effects (e.g., Morley et al., 2018) present a key avenue for more accurate descriptions of storm time radial diffusion coefficients, allowing for the high degree of observed variability in diffusion coefficients to be represented.

Conclusion
This study presents a statistical analysis of Van Allen Probes observations of ULF wave power and estimated radial diffusion coefficients during geomagnetic storms. Dependences on storm phase and L* were explored, and the temporal variation within storms were assessed using a superposed epoch analysis. Crucially, and in contrast to many previous studies, the local time dependence of wave power was accounted for when calculating the diffusion coefficients. We summarize the key findings here: 1. Storm time diffusion coefficients vary significantly and by several orders of magnitude during geomagnetic storms, experiencing enhancements in the late initial phase and maximizing in the late main phase. The results shown here provide storm time specific estimates of radial diffusion coefficients as a function of storm phase and L* 2. Analysis suggests that increases in radial diffusion coefficients is a consequence of both internal and external sources of ULF wave activity, and also that the inner magnetospheric conditions (ring current intensity and mass density distribution) play a key role in the propagation of wave power during geomagnetic storms 3. A comparison to the empirical model of Ozeke et al. (2014) highlights some key differences, most notably that the model tends to underestimate the magnetic field diffusion coefficient during storm times Throughout the analysis we have highlighted the large variability in radial diffusion coefficients during storm times. The range of values for a given epoch time spans several orders of magnitude for both the magnetic and electric components. Future analysis aims to explore the key contributors to variability during storm times. This will include a consideration of storm magnitude on shaping the diffusion coefficient values. The role of storm drivers (e.g., CME and CIR drivers) is also suggested to play an important role (Kalliokoski et al., 2020;Kilpua et al., 2015;Simms et al., 2010), and these effects on the magnitude and temporal variation of diffusion coefficients will also be assessed.

Data Availability Statement
The Van Allen Probes data is publicly available online (https://cdaweb.gsfc.nasa.gov/index.html/). The solar wind data and Sym-H index data are publicly available online (http://wdc.kugi.kyoto-u.ac.jp/index.html). V002554/2. I. J. Rae was supported by NERC Grants NE/P017185/2 and NE/ V002554/2. C. E. J. Watt was supported by STFC grant ST/R000921/1 and NERC grants NE/P017274/1 and NE/ V002759/2. R. B. Horne was supported by Highlight Topic Grant NE/ P01738X/1 (Rad-Sat). M.-T. Walach was supported by Natural Environments Research Council grants NE/P001556/1 and NE/T000937/1. The authors thank the EMFISIS instrument team for data provision. The authors thank the EFW team for data provision and the work by the EFW team was conducted under JHU/APL contract 922613 (RBSP-EFW).