The “SafeSpace” database of ULF power spectral density and radial diffusion coefficients: dependencies and application to simulations

Radial diffusion has been established as one of the most important mechanisms contributing to both the acceleration and loss of relativistic electrons in the outer radiation belt, as well as to the supply of particles to the inner radiation belt. In the framework of the “SafeSpace” project, we have used 9 years (2011–2019) of multi-point magnetic and electric field measurements from THEMIS A, D and E satellites to create a database of radial diffusion coefficients (D_LL) and ultra-low-frequency (ULF) wave power spectral densities (PSDs) spanning an L^∗ range from 3 to 8. In this work we investigate the dependence of the D_LL on the various solar wind parameters, geomagnetic indices and coupling functions, as well as the L-shell, during the solar cycle 24. Moreover, we discuss the uncertainties introduced on the estimation of D_LL time series by the partial azimuthal coverage provided by in situ measurements. Furthermore, we investigate, via a superposed analysis, the dependence of the D_LL on solar wind drivers. We show, for the first time to the best of our knowledge, that the interplanetary coronal mass ejection (ICME)-driven disturbances accompanied by high solar wind pressure values combined with intense magnetospheric compression can produce $D_{\mathrm{LL}}^{\mathrm{B}}$ values comparable to or even greater than the ones of $D_{\mathrm{LL}}^{\mathrm{E}}$. This feature cannot be captured by semi-empirical models and introduces a significant energy dependence on the D_LL. Finally, we show the advantages of using D_LL time series by means of numerical simulations of relativistic electron fluxes performed with the Salammbô code and significant deviations in the predictions of several semi-empirical models depending on the level of geomagnetic activity and L-shell.

1 Introduction

The outer radiation belt consists of electrons at energies from a few hundred kiloelectronvolts (keV) to several megaelectronvolts (MeV) (Daglis et al., 2019). Radial diffusion has been established as one of the most important mechanisms that contributes to this broad energy range of electrons since it can lead to both energization (Jaynes et al., 2015; Li et al., 2016; Katsavrias et al., 2019b; Nasi et al., 2020) and loss of relativistic electrons (Morley et al., 2010; Turner et al., 2012b; Katsavrias et al., 2015, 2019a).

Ultra-low-frequency (ULF) waves in the Pc4–5 band (1–25 mHz) can violate the third adiabatic invariant L^∗ of the energetic electrons. This drives radial diffusion by conserving the first two adiabatic invariants under the drift resonance condition ω=mω_d, where ω is the wave frequency, m is the azimuthal wave mode number and ω_d is the electron drift frequency (Elkington et al., 2003). Most often, radial transport is described as a stochastic process, the result of incoherent transport of particles by electromagnetic fields that vary irregularly on timescales of the drift period of radiation belt electrons (of the order of minutes). The radial diffusion coefficient, D_LL, has been defined to represent the mean square change of L^∗ for a large number of particles over time.

Currently there are two widely used formalisms in order to derive radial diffusion coefficients. Fälthammar (1965) distinguished the contribution of single-mode fluctuations in Earth's magnetic field and induced electric fields ($D_{\mathrm{LL}}^{\mathrm{M}}$) and perturbations in convection electric fields ($D_{\mathrm{LL}}^{\mathrm{E}}$) to derive a mathematical formulation for D_LL. On the other hand, Fei et al. (2006) included the contributions from all azimuthal wave modes. Nevertheless, the latter authors made the additional assumption that the magnetic field perturbations and the inductive electric field perturbations are independent, something that assumes that the two perturbations are uncorrelated. As discussed by Lejosne (2019), such an assumption is inconsistent with Faraday's law ($\mathrm{\nabla }x\mathit{E}=-\partial \mathit{B}/\partial t$).

Specifically, Fei et al. (2006) assumed radial diffusion coefficients to be the sum of the effects of perturbations in the azimuthal electric field and the parallel magnetic field:

$$\begin{array}{}\text{(1)}& D_{\mathrm{LL}}=D_{\mathrm{LL}}^{\mathrm{B}}+D_{\mathrm{LL}}^{\mathrm{E}}.\end{array}$$

These two components of the radial diffusion coefficients are given by

$$\begin{array}{}\text{(2)}& {\displaystyle }{D}_{\mathrm{LL}}^{\mathrm{B}}={\displaystyle \frac{{\mathit{\mu }}^{\mathrm{2}}{L}^{\mathrm{4}}}{\mathrm{8}{q}^{\mathrm{2}}{\mathit{\gamma }}^{\mathrm{2}}{B}_{\mathrm{E}}^{\mathrm{2}}{R}_{\mathrm{E}}^{\mathrm{4}}}}\cdot \sum _m{m}^{\mathrm{2}}{P}_m^{\mathrm{B}}\left(m{\mathit{\omega }}_{\mathrm{d}}\right)\text{(3)}& {\displaystyle }{D}_{\mathrm{LL}}^{\mathrm{E}}={\displaystyle \frac{L^{\mathrm{6}}}{\mathrm{8}{B}_{\mathrm{E}}^{\mathrm{2}}{R}_{\mathrm{E}}^{\mathrm{2}}}}\cdot \sum _m{P}_m^{\mathrm{E}}\left(m{\mathit{\omega }}_{\mathrm{d}}\right),\end{array}$$

where μ is the first adiabatic invariant, L is Roederer's L^∗, q is the charge of the diffused electrons, γ is the Lorentz factor, R_E is Earth's radius and B_E is the strength of the equatorial geomagnetic field on the Earth's surface. Moreover, P corresponds to the wave power at a specific drift frequency (ω_d) for all the azimuthal mode numbers (m). Note that $D_{\mathrm{LL}}^{\mathrm{B}}$ includes contributions only from the magnetic field oscillations, while $D_{\mathrm{LL}}^{\mathrm{E}}$ contains contributions from the total (inductive and convective) electric field.

It is clear, from the aforementioned formulation, that in order to have accurate calculations of the radial diffusion coefficients, we need accurate magnetic and electric field measurements, which, of course, are not always available. To that end, efforts have been devoted to provide empirical relationships of D_LL for radiation belt simulations, parameterizing the diffusion coefficients by the Kp index and L^∗ parameter. These empirical models have the advantage of providing estimations and predictions of the D_LL without the dependence on the in situ measurements. Nevertheless, it is also obvious (see also Table 1) that the use of a single input parameter (Kp index) is an oversimplification for a complex process such as the radial diffusion of electrons. Moreover, Kp is a global geomagnetic index, which is a proxy for the global changes in the geomagnetic field (Mayaud, 1980). On the other hand, two of the most important (external) sources for ULF waves are (a) solar wind pressure pulses and (b) Kelvin–Helmholtz instabilities powered by the increased solar wind speed (Claudepierre et al., 2008). Since the Kp index does not present significant correlation with either of these two solar wind parameters, it cannot account well for the mechanism of radial diffusion that enhances or depletes the electron population in the outer radiation belt.

Brautigam and Albert (2000)Boscher et al. (2018)Liu et al. (2016)Ozeke et al. (2014)Ali et al. (2016)

Table 1 Widely used semi-empirical models for the estimation and prediction of the radial diffusion coefficients, their mathematical formulation, trained datasets and limitations.

In addition, the observed D_LL has been shown to be highly event-specific (Jaynes et al., 2018), and physics-based models, such as the Versatile Electron Radiation Belt, cannot simulate the dynamics of the outer radiation belt observed during every storm using these empirically estimated coefficients (Drozdov et al., 2021). Several case studies have demonstrated deviations of the event-specific diffusion coefficients from the Kp-parameterized models. The recent study of Liu et al. (2018) suggests that the difference between the various models is negligible for low levels of geomagnetic activity at an equatorial distance of L-shell = 7.5 R_E but can be different by orders of magnitude at high levels of geomagnetic activity. At the same extent, Olifer et al. (2019) observed that, during the March 2015 geomagnetic storm, $D_{\mathrm{LL}}^{\mathrm{B}}$ was consistently underestimated and $D_{\mathrm{LL}}^{\mathrm{E}}$ was consistently overestimated by the empirical model of Ozeke et al. (2014). Furthermore, the magnitude of misestimation varied throughout the event, and, at times, the difference between empirically modelled values and time series of diffusion coefficients was multiple orders of magnitude.

In this work we present a new database of ULF power spectral density (PSD) and the derived radial diffusion coefficients, which has been developed in the framework of the “SafeSpace” project funded by Horizon 2020. The SafeSpace project aims at advancing space weather nowcasting and forecasting capabilities and, consequently, at contributing to the safety of space assets through the transition of powerful tools from research to operations. To that end, a database of radial diffusion coefficients derived from in situ magnetic and electric field measurements, coupled with solar wind and geomagnetic parameters, as well as the accompanied analysis, is of outmost importance, not only for statistical purposes, but also for any future efforts to develop accurate models for nowcasting and forecasting the D_LL. The rest of this paper is organized as follows: Sect. 2 describes the datasets used as input in the D_LL database as well as the roadmap towards its creation, Sect. 3 reports statistics which are important for future modelling efforts and Sect. 4 presents examples of the importance of the use of D_LL time series in radiation belt simulations.

2 ULF PSD and D_LL database

2.1 Data and methods

The radial diffusion coefficients were calculated directly from in situ measurements using the approach based on the Fei et al. (2006) formulation. As mentioned before, this approach considers the compressional component of the magnetic field and the toroidal component of the electric field. To that end, we used 4 s resolution measurements of the magnetic field vector from the THEMIS A, D and E fluxgate magnetometers (Auster et al., 2008) as well as electric field measurements from the EFI instrument (Bonnell et al., 2008) covering solar cycle 24 (2011–2019). Complementary measurements of solar wind and geomagnetic parameters were obtained from the NASA OMNIWeb database populated by NASA's Space Physics Data Facility with propagated values at the bow shock nose (https://omniweb.gsfc.nasa.gov/, last access: 8 June 2022). These complementary data were used not only for the parameterization of the database, but also for the statistical analysis performed in this study.

Figure 1 shows the steps followed in order to create the ULF power spectral density (PSD) and D_LL database, from the collection of the input data to the final data products. In detail, THEMIS magnetic and electric field vector data were preprocessed by transforming them into a mean field-aligned (MFA) coordinate system, similar to Balasis et al. (2013). The MFA is a local coordinate system defined by the ambient magnetic (or electric) field. The z axis is aligned with the unperturbed field (compressional or parallel component). The unperturbed field is obtained by a 30 min running average on the fields magnitude, and, then, the compressional component is calculated as follows:

$$\begin{array}{}\text{(4)}& B_{\text{com}}=\mathrm{\Delta }B\cdot {\displaystyle \frac{B}{\left|B\right|}}=(B-\stackrel{\mathrm{‾}}{B})\cdot {\displaystyle \frac{B}{\left|B\right|}},\end{array}$$

where $\stackrel{\mathrm{‾}}{B}$ is the unperturbed field. The y axis is perpendicular to the field's meridian pointing predominantly eastward (toroidal or azimuthal component), while the x axis completes the triad having an outward component (poloidal or radial component). After the transformation, the toroidal and compressional components of the electric and magnetic field, respectively, were de-trended using a 20 min moving average in order to eliminate the slow field variations. This is quite similar with a high-pass filtering with cut-off frequency at ≈0.83 mHz in order to focus on the Pc4–5 ULF frequencies.

Figure 1 Work logic towards the creation of the SafeSpace radial diffusion coefficient database.

The next step was to estimate the power spectral density (PSD) of the waves in the 2–25 mHz frequency range for the two time series, which corresponds to the drift periods of near-equatorial mirroring electrons roughly in the 0.4–13 MeV. For the spectral analysis of the electric and magnetic field measurements we made use of the continuous wavelet transform (CWT – see also Torrence and Compo, 1998) with the Morlet wavelet as the wavelet basis function (Morlet et al., 1983). Subsequently, the 4 s PSDs obtained were averaged in 1 min windows, which were considered one sample. For each time stamp, we also estimated the geomagnetic coordinates L^∗ and magnetic local time (MLT) using the International Radiation Belt Environment Modelling (IRBEM) library (Bourdarie and O'Brien, 2009) and the Olson–Pfitzer 1977 (Olson and Pfitzer, 1977) external magnetic field model. Note that, as the satellites move inbound and outbound with high velocities at low L-shells, the magnetic field measurements exhibit, not only orders of magnitude increase, but very large gradients as well. Therefore, after the application of the filtering, we perform a visual inspection of the time series, and the steep gradients that have not been removed by 20 min averaging are removed by hand.

Finally, using the estimated PSDs, the $D_{\mathrm{LL}}^{\mathrm{B}}$ and $D_{\mathrm{LL}}^{\mathrm{E}}$ were calculated following Eqs. (2) and (3), respectively. For the magnetic component, we have calculated $D_{\mathrm{LL}}^{\mathrm{B}}$ for 39 values of the first adiabatic invariant (μ) in the 50–20 000 MeV G⁻¹ range.

The PSDs of both the toroidal electric and the compressional magnetic field as a function of time, L^∗ and magnetic local time (MLT), were stored in daily CDF files for each THEMIS probe, separately. Then, $D_{\mathrm{LL}}^{\mathrm{B}}$ and $D_{\mathrm{LL}}^{\mathrm{E}}$ were calculated as hourly values from the averaged PSD of all three THEMIS probes and grouped in bins $\mathrm{d}{L}^{\ast }=\mathrm{0.1}$ (the process and assumptions are discussed in detail in the following sections). These grouped radial diffusion coefficients were also stored in daily CDF files. The database of both the ULF PSD and the D_LL can be found at https://synergasia.uoa.gr/modules/document/?course=PHYS120 (last access: 8 June 2022).

2.2 Weighted average power

The wave power included in Eqs. (2) and (3) corresponds to the power at a specific drift frequency for all m values, which essentially means that particles are radially transported via stochastic acceleration with various frequency waves (main frequency and harmonics). Nevertheless, to calculate the power at various m values, one would need at least 2m observations simultaneous in time, which is not trivial. To address this issue, it is often assumed that power at high m values is consistently lower than power at m=1 and, subsequently, that all power is contained in the lowest m=1 wave mode of ULF waves driving diffusion (Ozeke et al., 2014). This assumption can lead to underestimation of the radial diffusion coefficient, since higher m values are shown to be often significant (e.g. m=2 up to m=5 at recovery phase of storms; see also Sarris et al., 2013). To address this issue, we have opted to use, in the place of power at a specific frequency, the weighted average power over the whole frequency range under study (in our case Pc4 and Pc5 frequency range). This weighted average power is given by Torrence and Compo (1998) as follows:

$$\begin{array}{}\text{(5)}& P_{\text{total}}=\left({\displaystyle \frac{\mathrm{d}j\cdot \mathrm{d}t}{\text{Cdelta}}}\right)\cdot \sum _f\text{PSD}\left(f\right),\end{array}$$

where Cdelta is a smoothing factor (which for the Morlet wavelet is empirically derived as 0.776) and

$$\mathrm{d}j=-{\displaystyle \frac{{\log }_{\mathrm{2}}\left(\frac{f_{\text{min}}}{f_{\text{max}}}\right)}{\mathrm{1}/f_{\text{min}}}}$$

is the sampling scale (for more details, see also Sect. 3 of Torrence and Compo, 1998, and Sect. 4 of Katsavrias et al., 2022a).

2.3 Assumptions

Even though we have followed a well-established methodology in order to calculate – as accurately as possible – the ULF PSD and the corresponding D_LL, there are still assumptions worth mentioning, which are based on the theoretical approach we have used as well as on the inherent limitations of the in situ data.

As already discussed, important differences can exist between the two approaches by Fei et al. (2006) and by Fälthammar (1965). It is estimated that the former can underestimate (compared to the latter) the total D_LL by a factor of 2 (Lejosne, 2019), given that the magnetic field disturbances are described by the simple model introduced by Fälthammar (1965) and that there is no electric potential disturbance. Nevertheless, Fei's approach is the more widely used due to the fact that it is very difficult to separate the total measured electric field from single point measurements in space (Brautigam et al., 2005) into its convective and inductive components (Lejosne and Kollmann, 2020). Furthermore, it has been shown that this discrepancy by a factor of 2 is comparatively minor relative to the large variability in the observed values (Sandhu et al., 2021).

In addition, the theoretical approach of Fei et al. (2006) formulas applies for equatorially mirroring particles only, while THEMIS satellites do not necessarily sample the magnetic equator. Nevertheless, they remain very close to the magnetic equator throughout their trajectories in the heart of the outer belt (Angelopoulos, 2008; Turner et al., 2012a), something that allows us to assume that the uncertainty in the D_LL calculation will be rather small. Moreover, the distribution of both the electric and magnetic field power with respect to the B ratio ($B_{\text{eq}}/B_{\text{local}}$, not shown here) shows that the uncertainty in ULF power is up to a factor of 2, at least for B ratio values larger than 0.8. We emphasize that this 0.8 threshold in the B ratio is applied in our database in order to minimize the corresponding uncertainty; thus all the following results correspond to D_LL values at points with B ratio>0.8. We also note that it is difficult to differentiate between spatial and temporal field variations. In our calculations we have assumed that since the spatial variations, especially at the radial distances corresponding to the outer belt, are usually slow variations, they are filtered out by the filtering process. This is of course usually not possible at very low L-shells; thus the corresponding calculations are removed from the dataset.

Finally, we have to mention the uncertainties introduced by the use of the Olson–Pfitzer quiet model for the estimation of the magnetic ephemeris data, as well as the uncertainties of the instruments used (e.g. Califf et al., 2015).

3 Uncertainties generated by the azimuthal dependence of ULF power

Equations (2) and (3) implicitly assume a uniform distribution of wave power in azimuth. In reality, the azimuthal distribution of the wave power in the Pc4–5 range depends on their generation mechanism; e.g. the wave power due to the Kelvin–Helmholtz instability is expected to be greater near dawn and dusk sectors, while due to the pressure, pulses from the solar wind are expected to be greater near noon. Furthermore, radial diffusion coefficients are by definition drift-averaged, which is opposite in paradigm to these existing physical anisotropies. However, the in situ calculation of radial diffusion coefficients with full azimuthal coverage, and therefore averaging, would require a large spacecraft constellation with appropriate positioning providing concurrent and intensity coherent measurements. Obviously this is not currently possible. Our efforts have been focused here on quantifying the magnitude of radial diffusion due to ULF waves observed by the THEMIS spacecraft, which provide a maximum MLT coverage that rarely exceeds 6 h per hour and per L^∗.

In order to discuss the possible uncertainties generated by the limited azimuthal coverage, we calculate a 1 min resolution proxy of the D_LL at each point of the spacecraft orbit, for each spacecraft separately. Since the D_LL proxy at each L^∗ value has been calculated as the product of the weighted averaged power with a simple multiplication factor, it is expected to reflect directly the azimuthal distribution of wave power for both the magnetic and the electric component. Figure 2 shows the logarithms of the mean D_LL proxy as a function of MLT and L^∗ for three levels of geomagnetic activity: Kp<3 (left column panels), $\mathrm{3}<\text{Kp}<\mathrm{5}$ (middle column panels) and Kp>5 (right column panels).

Figure 2 Logarithms of the mean D_LL proxy as a function of MLT (dMLT=1 h) and L^∗ ($\mathrm{d}{L}^{\ast }=\mathrm{0.1}$) for three levels of geomagnetic activity: (a, d) Kp<3, (b, e) $\mathrm{3}<\text{Kp}<\mathrm{5}$ and (c, f) Kp>5. Panels (a)–(c) and (d)–(f) correspond to the electric and the magnetic component (for μ=1000 MeV G⁻¹) of the D_LL proxy, respectively.

As shown, there are obvious differences between the two components. During quiet times, the $D_{\mathrm{LL}}^{\mathrm{E}}$ proxy (top left panel) exceeds the value of 10 d⁻¹ outside the geosynchronous orbit and is approximately equal to 1 d⁻¹ at the 4.5–6 L^∗ range, while there is a significant MLT asymmetry at the 00:00–03:00 MLT range. As we move to higher geomagnetic activity levels ($\mathrm{3}<\text{Kp}<\mathrm{5}$ – top middle panel), the $D_{\mathrm{LL}}^{\mathrm{E}}$ proxy intensifies, and, in addition, this asymmetry becomes stronger and reaches even lower L^∗ values. During intense geomagnetic activity levels (top right panel), $D_{\mathrm{LL}}^{\mathrm{E}}$ proxy values range between 10 and 100 d⁻¹ at $L^{\ast }>\mathrm{5}$, and they reach approximately the value of 1 d⁻¹ even down to at $L^{\ast }=\mathrm{3.5}$. On the other hand, the $D_{\mathrm{LL}}^{\mathrm{B}}$ proxy exhibits a different behaviour. During quiet times (bottom left panel), the $D_{\mathrm{LL}}^{\mathrm{B}}$ proxy values reach 1 d⁻¹ at $L^{\ast }>\mathrm{7}$ and only at the dayside sector (approximately in the 09:00–15:00 MLT range). As we move to higher geomagnetic activity levels, the $D_{\mathrm{LL}}^{\mathrm{B}}$ proxy exceeds the value of 10 d⁻¹, even inside the geosynchronous orbit $L^{\ast }<\mathrm{6}$. Furthermore, the MLT asymmetry becomes more intense and wide (approximately in the 05:00–18:00 MLT range during Kp>5 periods).

The aforementioned feature of the $D_{\mathrm{LL}}^{\mathrm{B}}$ proxy is in agreement with the correlation results shown in Fig. 3 and indicates that the magnetic field, and consequently the magnetic component of the diffusion coefficient, is linked with ULF waves generated through solar wind pressure pulses (Kepko et al., 2002). On the other hand, the observed asymmetry in the $D_{\mathrm{LL}}^{\mathrm{E}}$ proxy is not only linked with solar wind speed, but also with internal mechanisms such as substorm activity, especially during quiet or moderate magnetospheric activity. This is supported by the remarkable agreement of the $D_{\mathrm{LL}}^{\mathrm{E}}$ proxy with the findings of Nosé et al. (1998), who showed that substorms generate azimuthal ULF fluctuations at the nightside which peak at 01:00–02:00 MLT. Furthermore, this is also in agreement with the results of Fig. 3 and the significant correlation of $D_{\mathrm{LL}}^{\mathrm{E}}$ with the Auroral Electrojet (AE) index. These results are also in good agreement with Sandhu et al. (2021), who used Van Allen Probe data to infer the radial diffusion coefficients. This agreement indicates that the uncertainty introduced by the magnetic latitude (and already discussed in Sect. 2) is insignificant, even though there is no straightforward comparison with the dataset used by the latter authors.

Figure 3 Pearson correlation coefficients between the natural logarithms of the hourly mean values of $D_{\mathrm{LL}}^{\mathrm{E}}$ (a) and $D_{\mathrm{LL}}^{\mathrm{B}}$ (for μ=1000 MeV G⁻¹ – b) with various solar wind parameters, geomagnetic indices and coupling functions as a function of L^∗ (with $\mathrm{d}{L}^{\ast }=\mathrm{0.1}$).

Our D_LL calculations employ a fraction of the full azimuthal coverage and, thus, are expected to have some uncertainty. As shown in Fig. 2, the asymmetry in the electric component of the D_LL proxy is limited in the 00:00–03:00 MLT range for low to moderate geomagnetic conditions, which means that the uncertainty in the $D_{\mathrm{LL}}^{\mathrm{E}}$ will not be significant given the up to 6 h coverage in MLT from the three THEMIS spacecraft. In addition, the magnetic component of the D_LL proxy during low geomagnetic activity does not exhibit significant variation with the MLT, which also corresponds to insignificant uncertainties in the $D_{\mathrm{LL}}^{\mathrm{B}}$. On the other hand, the $D_{\mathrm{LL}}^{\mathrm{B}}$ proxy during moderate and intense geomagnetic activity exhibits a difference of approximately 1 order of magnitude between dayside and nightside. This means that the partial azimuthal coverage provided by the three THEMIS spacecraft could lead to an uncertainty of up to 1 order of magnitude in the estimation of the $D_{\mathrm{LL}}^{\mathrm{B}}$ for particular spatial configurations of the three THEMIS spacecrafts, e.g. when all three are located in the nightside or all three are located in the dayside. We note that such uncertainties are present in all D_LL time series estimated by in situ measurements, e.g. Jaynes et al. (2018), Olifer et al. (2019) and Sandhu et al. (2021) as well as in semi-empirical models; e.g. Ozeke et al. (2014) have used only dayside measurements from ground-based magnetometers to infer the electric component of the D_LL.

4 Dependence on solar wind and geomagnetic parameters

Figure 3 shows the Pearson correlation coefficients (henceforward CCs) between the logarithm of the hourly mean values of D_LL with various solar wind parameters, geomagnetic indices and coupling functions in the 3–8 L^∗ range. Shown, from left to right, are the interplanetary magnetic field (IMF), solar wind velocity, pressure and number density, plasma β parameter, the geomagnetic indices SYM-H, AE and Kp, the ϵ parameter (Akasofu, 1981), the southward solar wind field (here we show the exponential of Bs), the Half-Wave Rectifier (Burton et al., 1975) and Newell's function (Newell et al., 2007).

Generally, the CCs of the magnetic component exhibit greater values than the ones of the electric component with maxima at ≈0.7 and ≈0.4, respectively. Note that we only show the $D_{\mathrm{LL}}^{\mathrm{B}}$ at 1000 MeV G⁻¹, but the CCs do not change at all if we account for the μ value. In detail, both D_LL components exhibit their best correlation with the geomagnetic indices AE and Kp. Nevertheless there is a pronounced difference concerning the L^∗ location of the maximum CC. For the electric component, the maximum CC (≈0.4 for both AE and Kp) is located roughly in the 4.5–6.5 L^∗ range. For the magnetic component, the maximum CC with AE (≈0.65) is located roughly at the 4.5–8 L^∗ range and the maximum CC with Kp (≈0.7) covers approximately the whole L^∗ range. The latter is in agreement with Dimitrakoudis et al. (2015), who found that the Kp index provides the best parameterization of the $D_{\mathrm{LL}}^{\mathrm{B}}$. Our results indicate that this parameterization may not work equally for the electric component, especially for L^∗ values higher than 6.5 and lower than 4.5.

Furthermore, the CC between solar wind speed and D_LL is at ≈0.4 and ≈0.5 for the electric and magnetic component, respectively, but both in the 4.5–6.5 L^∗ range. The importance of magnetopause instabilities – induced by the increased solar wind velocity – has been well established before (Bentley et al., 2018), but here we show that it can similarly affect both D_LL components. Another interesting feature is exhibited by the correlation between the D_LL and solar wind dynamic pressure, even though there is no significant correlation with number density. For the electric component, the CC does not exceed the 0.2 value, but for the magnetic component, it is larger than 0.5 at $L^{\ast }>\mathrm{4.5}$. A possible explanation of this feature could be that, since solar wind pressure pulses produce mainly global magnetospheric oscillations (Kepko et al., 2002; Takahashi et al., 2012), they do not affect the azimuthal electric field variations and thus the electric D_LL component.

It is worth mentioning that the only parameter which exhibits an anti-correlation with the D_LL is the plasma β parameter at all L^∗ values. Nevertheless, the maximum CC at both components does not exceed −0.2. Finally, the CCs between the $D_{\mathrm{LL}}^{\mathrm{B}}$ component with Newell's function and Akasofu's ϵ parameter exhibit a similar trend to the AE index but with lower CC maxima (≈0.4). This is expected since these parameters are known to be well correlated with substorm activity (Katsavrias et al., 2021).

Figure 4 shows the cross-correlation between $D_{\mathrm{LL}}^{\mathrm{E}}$ (left panels) and $D_{\mathrm{LL}}^{\mathrm{B}}$ (for μ=1000 MeV G⁻¹ – right panels) with (top to bottom) solar wind speed, dynamic pressure, plasma β, AE and Kp index. Note that in this figure we are showing only the parameters which, according to Fig. 3, exhibited noteworthy correlations. Similar to Fig. 3, the CCs of the magnetic component are systematically higher than the ones of the electric component, at least for time lags of up to 12 h, with the exception of plasma β. As shown, the maximum CCs for the magnetic component (right panels) are exhibited at zero time lag, while they become negligible for time lags greater than 9 h. A similar trend is exhibited for the CCs of the electric component with solar wind speed and AE index. On the contrary, the CC of the electric component with Kp index exhibits a maximum at the 0–3 h time lag.

Figure 4 Pearson correlation coefficients between the natural logarithms of the hourly mean values of $D_{\mathrm{LL}}^{\mathrm{E}}$ (left panels) and $D_{\mathrm{LL}}^{\mathrm{B}}$ (for μ=1000 MeV G⁻¹ – right panels) with (top to bottom) solar wind speed, dynamic pressure, plasma β, AE and Kp index. CCs are presented – both colour-coded and with black contours – as a function of L^∗ (with $\mathrm{d}{L}^{\ast }=\mathrm{0.1}$) and the time lag from 0 to 36 h.

5 ICME- vs. SIR-driven geospace disturbances

The role of solar wind drivers (e.g. interplanetary coronal mass ejections, ICMEs, and stream interaction regions, SIRs) has been suggested to play an important role in the generation of ULF waves and, consequently, in the evolution of radial diffusion coefficients (Simms et al., 2010; Kilpua et al., 2015). In order to investigate the dependence of the $D_{\mathrm{LL}}^{\mathrm{E}}$ and $D_{\mathrm{LL}}^{\mathrm{B}}$ on the solar wind driver we have selected 25 ICME- and 46 SIR-driven geospace disturbances (71 events in total) in the 2011–2019 time period, following the criteria of Katsavrias et al. (2019a). More specifically, we have chosen events that include a single driver and have no pre-conditioning in solar wind parameters for at least 12 h before the arrival of the ICME or SIR. Since we have applied no criteria depending on the Dst index (non-storm events are also included), we have used the time of the maximum compression of the magnetopause (Lmp_min)as zero-epoch time (t₀) as it is given by the empirical model of Shue et al. (1998).

Figure 5 shows the results of the superposed epoch analysis. As shown, both groups exhibit several differences. During ICME-driven disturbances, the maximum increase in $D_{\mathrm{LL}}^{\mathrm{E}}$ takes place on t₀ at all $L^{\ast }>\mathrm{4}$ and reaches a median value of 1000 d⁻¹ at $L^{\ast }>\mathrm{5}$, while significant activity reaches down to L≈3.5 up to 12 h. After these 12 h, the activity is still significant at $L^{\ast }>\mathrm{5}$ and lasts up to 96 h (4 d). During SIR-driven disturbances, the $D_{\mathrm{LL}}^{\mathrm{E}}$ exhibits quite a similar trend (it lasts up to 4 d after t₀), but both its magnitude and the penetration to inner L^∗ are lower compared to the ICME-driven disturbances. On the other hand, the $D_{\mathrm{LL}}^{\mathrm{B}}$ exhibits much more pronounced differences. During ICME-driven disturbances, the maximum increase in $D_{\mathrm{LL}}^{\mathrm{B}}$ takes place on t₀, and the penetration of the activity reaches down to $L^{\ast }\approx \mathrm{4}$. The overall enhancement occurs on $-\mathrm{8}<t_{\mathrm{0}}<\mathrm{12}$ h. During SIR-driven disturbances, the $D_{\mathrm{LL}}^{\mathrm{B}}$ hardly reaches $L^{\ast }\approx \mathrm{4}$, and the maximum increase reaches a value of 10. Nevertheless, the overall activity lasts up to approximately 30 h after t₀. Furthermore, the enhancement, as well as the penetration of $D_{\mathrm{LL}}^{\mathrm{B}}$ to low L^∗, is very well correlated with the enhancement in both solar wind dynamic pressure and Kp index and, consequently, is in agreement with the findings of Fig. 3. This result is also in agreement with Simms et al. (2010), who indicated that ground Pc5 power was greater during coronal mass ejection (CME) storms, especially during the main and recovery phase. One step further, Kalliokoski et al. (2020) studied 37 ICME-driven sheath regions in the Van Allen Probe era and linked the increased Pc5-ULF activity at GEO with the increased pressure during the sheath.

Figure 5 Superposed epoch analysis of the 25 ICME-driven (left column panels) and 46 SIR-driven (right column panels) geospace disturbances. Top to bottom: median (black line), 25th and 75th quantiles (red lines) of the magnetopause location predicted by the Shue et al. (1998) model, solar wind speed, Kp index, solar wind dynamic pressure, the logarithm of the median values of $D_{\mathrm{LL}}^{\mathrm{E}}$, $D_{\mathrm{LL}}^{\mathrm{B}}$ (for μ=1000 MeV G⁻¹) and their ratio. The binning is performed with dt=1 h and $\mathrm{d}{L}^{\ast }=\mathrm{0.1}$.

Finally, a very important feature is exhibited by the ratio of the electric over the magnetic component, which generally spans the 0.1–100 range. As shown in the bottom panels of Fig. 5, the electric component is mostly dominant – up to 2 orders of magnitude compared with the magnetic component. This feature changes dramatically during ICME-driven disturbances and around ±6 h from the maximum compression of the magnetopause, where the D_LL ratio decreases below 1 at all L^∗ values. Furthermore, at $L^{\ast }>\mathrm{6}$, the D_LL ratio is approximately 1 up to 12 h after t₀. The relative strength of the two D_LL has been discussed before by Olifer et al. (2019), who studied the component ratio during the St. Patrick's event of 2015. These authors indicated that during the main phase of this ICME-driven storm, the magnetic component exceeded the electric by approximately 1 order of magnitude, something that semi-empirical models cannot reproduce. Here we replicate this result using a statistical sample of 25 ICME-driven disturbances independent of the magnitude of Dst index. Note that this feature, even though it is not that obvious, may be important during SIR disturbances as well. As shown in the bottom right panel of Fig. 5, the D_LL ratio at $L^{\ast }>\mathrm{5.5}$ is decreased from approximately 100 to approximately 1 at ±3 h from t₀. We suggest that this difference in the D_LL ratio between ICME- and SIR-driven disturbances is probably attributed to the existence (or not) of shocks, which produce a significant increase of the dynamic pressure and accompany, more often, the ICME-driven events.

We must emphasize the fact that this feature introduces a significant energy dependence on the D_LL, since the magnetic component is energy-dependent, that may be of great importance to radiation belt simulations. Furthermore, this feature is expected to be dependent on the first adiabatic invariant as well, since greater values of μ produce greater values of $D_{\mathrm{LL}}^{\mathrm{B}}$, which will consequently lead to changes in the D_LL ratio. It is also expected that, except the magnitude, the change in μ will affect both the duration and the L^∗ coverage of this feature as well. In a future study, we intend to investigate these changes in greater detail.

6 The use of D_LL time series in physics-based models

6.1 Comparison with semi-empirical models

As already discussed in the Introduction, even though the semi-empirical Kp-parameterized models have the advantage of providing estimations and predictions of the D_LL without the dependence on the in situ measurements, they can significantly deviate from the calculated D_LL time series. In order to statistically establish these deviations we directly compare the calculated D_LL values ($D_{\mathrm{LL}}^{\mathrm{B}}$ is always at μ=1000 MeV G⁻¹) from the SafeSpace database to the empirically modelled values of Table 1 for the whole 2011–2019 time period.

We note that the comparison of radial diffusion coefficients among multiple methods is anything but a straightforward process since the details of each method are different (e.g. different datasets, different time periods in a solar cycle or different theoretical approaches). In addition, here we attempt a comparison of a dataset inferred from in situ measurements with the estimations of semi-empirical models, which are by definition two different things. Nevertheless, we attempt a comparison in order to show a clearer picture of the statistical behaviour of each model compared to our D_LL inferred from in situ data.

Figure 6 shows this comparison parameterized by the Kp index. As shown, there is a general trend with all empirical models (and their components), where the D_LL has lower values at low levels of geomagnetic activity at all L^∗ and higher values at high levels of geomagnetic activity at high L^∗ values, compared to the SafeSpace D_LL. In detail, concerning the Ozeke et al. (2014) model, the $D_{\mathrm{LL}}^{\mathrm{E}}$ exhibits lower values at all L^∗ for Kp<4 (panel a), while there is a relatively good agreement with SafeSpace D_LL at $\mathrm{4}<\text{Kp}<\mathrm{6}$. The observed higher values at $L^{\ast }>\mathrm{7}$ and Kp>6 cannot be discussed since they exceed the limits of the Ozeke model (see also Table 1). These features are in good agreement with the results of Sandhu et al. (2021) and Murphy et al. (2015), even though the latter authors performed a statistical comparison during storm time only and with D_LL values calculated using RBSPs (Radiation Belt Storm Probes) and ground-based data, respectively. On the other hand, the $D_{\mathrm{LL}}^{\mathrm{B}}$[OZ] (panel b) exhibits persistently lower values of at least a factor of 10 at all L^∗ and Kp values, which is in agreement with the results of Olifer et al. (2019). Also note that this feature is expected to be even more pronounced with increasing μ values, since the $D_{\mathrm{LL}}^{\mathrm{B}}$[OZ] is not energy-dependent. A possible source for the aforementioned disagreement between the Ozeke model and our calculated D_LL may come from the assumption that ULF wave power is concentrated in the lowest mode (m=1), thus underestimating the true ULF wave power, while in this study, a weighted averaging over all frequencies has been used. Furthermore, the $D_{\mathrm{LL}}^{\mathrm{E}}$[OZ] is based on ground magnetometer measurements restricted in the dayside sector with the assumption that the observed power is independent of the MLT. This is in contrast with our results in Fig. 2, where the contribution from internal mechanisms (i.e substorm activity) are shown to be important. The most striking difference is the persistent disagreement of the $D_{\mathrm{LL}}^{\mathrm{B}}$[OZ], even though it is also derived by THEMIS magnetic field measurements. We suggest that the primary sources of this feature are the m=1 assumption and the fact that its training dataset is restricted to the 2007–2011 time period, thus including mostly quiet magnetospheric conditions.

Figure 6 Comparison of the SafeSpace D_LL values with the five semi-empirical models listed in Table 1 binned in Kp index values (dKp=0.5) and L^∗ ($\mathrm{d}{L}^{\ast }=\mathrm{0.1}$). $D_{\mathrm{LL}}^{\mathrm{B}}$ corresponds to μ=1000 MeV G⁻¹.

It is worth mentioning that the Brautigam and Albert (2000) model (panel e) exhibits a very similar trend with the electric component of the Ozeke et al. (2014) model, even though the two models have been developed using completely different theoretical approaches and datasets.

Concerning the Ali et al. (2016) model, both components exhibit significantly lower values compared to the SafeSpace D_LL that reaches approximately 2 orders of magnitude (panels c and d), with the exception of $D_{\mathrm{LL}}^{\mathrm{E}}$[ALI], which appears higher at high L^∗ and Kp values (top right corner of panel c). Nevertheless, this area is outside the limits of the model as described in Table 1. The overall behaviour of the Ali et al. (2016) model presented in this figure is in agreement with the results of Drozdov et al. (2021), who showed that simulations performed with the Versatile Electron Radiation Belt (VERB) code using this D_LL model exhibited significantly lower flux levels. Similarly with Ozeke et al. (2014), the D_LL[ALI] values are calculated with the assumption that the observed power is independent of the MLT and that all power falls into m=1. Moreover, the Ali et al. (2016) model considers ULF wave power in a narrow frequency range (1.67–6.67 mHz), while in this study we have considered the full Pc4–5 range up to 25 mHz. We suggest that the aforementioned differences possibly account for the lower values of the D_LL[ALI] in comparison with our calculated ones.

To the same extent, the Liu et al. (2016) model for the $D_{\mathrm{LL}}^{\mathrm{E}}$ (panel g) exhibits mostly lower values, compared to our calculated ones, up to a factor of 10, even though it also uses THEMIS electric field measurements to derive the electric component of the D_LL. In addition, the Liu model exploits the entire Pc3–5 frequency range to derive the $D_{\mathrm{LL}}^{\mathrm{E}}$, which is quite similar to the range used in this study, but the same assumption that all power falls into m=1 is used here as well. Nevertheless, a significant difference is that the dataset used to derive the model spans the 2008–2014 time period. This time period includes the extended minimum of 2008–2009 (which is not included in our dataset) and misses the recovery phase of solar cycle 24, where several intense events occurred.

The Boscher et al. (2018) model (panel f) exhibits very similar results compared to the SafeSpace D_LL. In detail, the modelled D_LL is in good agreement for $\mathrm{4}<\text{Kp}<\mathrm{7}$ at $L^{\ast }<\mathrm{6}$ and for $\mathrm{3}<\text{Kp}<\mathrm{6}$ at $L^{\ast }>\mathrm{6}$. Nevertheless, it exhibits significantly lower values of the D_LL up to a factor of 10 during quiet times at all L^∗ and significantly higher values (at least a factor of 10) for Kp>6 approximately outside the geosynchronous orbit, which nevertheless is outside the limits of the model.

Finally, we have to consider the uncertainties in the SafeSpace-calculated D_LL that may also be a cause of disagreement versus the aforementioned semi-empirical models. As discussed in Sect. 2, it has been shown that the Fei et al. (2006) approach can underestimate the radial diffusion coefficient by a factor of 2 compared with the Fälthammar (1965) approach. This is sufficient to explain the difference exhibited by SafeSpace D_LL and the Brautigam and Albert (2000) and Boscher et al. (2018) models at $L^{\ast }>\mathrm{4}$, but it cannot explain the difference of up to a factor of 10 at lower L-shells. Another uncertainty, also discussed in Sect. 2, comes from the limited azimuthal coverage of THEMIS satellites used in this study. Nevertheless, the results of Fig. 6 are averaged values of the SafeSpace D_LL for specific values of L^∗ and Kp over a 9 years time period, thus including all MLT values. Last but not least, we have to consider the uncertainties introduced by the use of the OP77 model, especially at high L-shells (outside GEO) and high Kp values.

6.2 St. Patrick's 2015 event

In the previous section, we presented an extended comparison of the various semi-empirical models with the calculated D_LL from the SafeSpace database, showing that all of them (more or less) exhibit significant deviations at different L^∗ and Kp values range. These deviations correspond to the cause (D_LL) and not the effect (electron radial diffusion). In order to evaluate the actual effect of these calculated radial diffusion coefficients on the outer belt dynamics, we have performed simulations without the energy diffusion term using the Salammbô model (Beutier and Boscher, 1995). Figure 7 shows the results of this simulation for two electron energies at 500 (left column panels) and 1500 keV (right column panels) during the March 2015 time period, which includes the St. Patrick's event of 17 March. Note that the magnetospheric model used in the simulations is the Olson–Pfitzer quiet model. In addition, the dataset gaps due to the partial azimuthal coverage of THEMIS constellation were filled using power law interpolation and extrapolation schemes in L^∗, and then a time interpolation at each L^∗ value was performed.

Figure 7 Simulation of the outer radiation belt dynamics during the St. Patrick's event of 2015 using the Salammbô model for two electron energies: (left column panels) 500 keV and (right column panels) 1.5 MeV. From bottom to top the Kp index, the AL index, the electron flux measured by the MagEIS instrument on board the Van Allen Probes, the simulation results using the Boscher et al. (2018) semi-empirical model and the simulation results using the SafeSpace D_LL values are shown.

As shown in the 500 keV electron energy, simulation results exhibit more intense radial transport at the outer edge of the outer radiation belt ($\mathrm{4}<L^{\ast }<\mathrm{5.5}$) both during the relatively quiet period on early March and during the intense St. Patrick's storm when using the calculated D_LL compared to the Boscher et al. (2018) model. This is in agreement with the results shown in Fig. 6, where the semi-empirical models underestimate the D_LL at high L^∗ values during active geomagnetic conditions. Moreover, as shown in the 1500 keV electron energy, the simulation captures not only the re-distribution of the relativistic electron population, but also generally the dynamics and the magnitude of the 1500 keV electron fluxes more realistically. The latter is particularly important since it has been reported that during the St. Patrick's event of 2015, radial diffusion contributed not only to the enhancement of 1–2 MeV electrons but also to further acceleration to ultra-relativistic energies (Jaynes et al., 2018). We note that the magnitude of the flux in the Salammbô simulation is not expected to agree with the MagEIS data due to the lack of the energy diffusion term (in situ acceleration by VLF chorus waves), which for the St. Patrick's event of 2015 has been shown to be crucial, especially for 1–2 MeV electrons (Li et al., 2016).

We must emphasize the fact that the aforementioned comparison is performed between the calculated μ-dependent D_LL from the SafeSpace database and the Boscher et al. (2018) model, only. This is done in accordance to the results discussed in the previous section (see also Fig. 6), where we showed that the Boscher model exhibited the best comparison with the SafeSpace calculated diffusion coefficients. It is also notable that the simulation results indicate that the uncertainties introduced by the limited azimuthal coverage provided by the THEMIS spacecraft are rather small compared to the uncertainties of using medians or other statistical quantities over long-term datasets.

7 Conclusions

In the framework of the SafeSpace project, we have used 9 years (2011–2019) of multi-point magnetic and electric field measurements from THEMIS A, D and E satellites to create a database of ULF power spectral density and the estimated radial diffusion coefficients. We have further exploited this database in order to investigate the dependence of these calculated D_LL values on several solar wind and geomagnetic parameters, and on solar wind drivers (ICMEs and SIRs), with respect to L^∗.

The results of this analysis can be summarized as follows:

1. $D_{\mathrm{LL}}^{\mathrm{B}}$ exhibits higher correlation coefficients, with both solar wind and geomagnetic parameters, compared to $D_{\mathrm{LL}}^{\mathrm{E}}$. In addition, $D_{\mathrm{LL}}^{\mathrm{E}}$ exhibits its highest coefficients in the 4.5–6.5 L^∗ range and $D_{\mathrm{LL}}^{\mathrm{B}}$ at $L^{\ast }>\mathrm{4.5}$.

2. Both D_LL components (magnetic and electric) exhibit good correlation with the Kp and AE index. Furthermore, $D_{\mathrm{LL}}^{\mathrm{E}}$ exhibits good correlation with solar wind speed, while $D_{\mathrm{LL}}^{\mathrm{B}}$ exhibits good correlation with both solar wind speed and pressure with zero time lag.

3. The superposed epoch analysis reveals important differences between the evolution of D_LL during ICME- and SIR-driven disturbances. During the former, the high solar wind pressure values combined with the intense magnetospheric compression produce $D_{\mathrm{LL}}^{\mathrm{B}}$ values comparable to or even greater than the ones of $D_{\mathrm{LL}}^{\mathrm{E}}$. This feature cannot be captured by semi-empirical models and introduces a significant energy dependence on the D_LL.

Furthermore, the comparison of the semi-empirical models with the D_LL from the SafeSpace database reveals significant deviations depending on the level of geomagnetic activity and the drift shell. Generally, all models exhibit lower values of the D_LL during quiet times at low L^∗ values and higher values during high levels of geomagnetic activity at high L^∗ values, compared to the SafeSpace D_LL. Finally, we have evaluated these calculated D_LL values through simulations of relativistic electrons using the Salammbô code. The simulations show that the model captures not only the re-distribution of the relativistic electron population, but also generally the dynamics and the magnitude of the electron fluxes more realistically when using the SafeSpace D_LL compared with the well-established Boscher semi-empirical model.

We believe that these results may offer significant insight for future modelling efforts in order to develop an accurate nowcasting and forecasting model for radial diffusion coefficients.