Using ARMAX Models to Determine the Drivers of 40-150 keV GOES Electron Fluxes.

We investigate the drivers of 40-150 keV hourly electron flux at geostationary orbit (GOES 13) using autoregressive moving average transfer functions (ARMAX) multiple regression models which remove the confounding effect of diurnal cyclicity and allow assessment of each parameter independently. By taking logs of the variables, we create nonlinear models. While many factors show high correlation with flux in single variable analysis (substorms, ULF waves, solar wind velocity (V), pressure (P), number density (N) and electric field (E y ), IMF Bz, Kp, and SymH), ARMAX models show substorms are the dominant influence at 40-75 keV and over 20-12 MLT, with little difference seen between disturbed and quiet periods. The Ey influence is positive post-midnight, negative post-noon. Pressure shows a negative influence, strongest at 150 keV. ULF waves are a more modest influence than suggested by single variable correlation. Kp and SymH show little effect when other variables are included. Using path analysis, we calculate the summed direct and indirect influences through the driving of intermediate parameters. Pressure shows a summed direct and indirect influence nearly half that of the direct substorm effect. N, V, and B z , as indirect drivers, are equally influential. While simple correlation or neural networks can be used for flux prediction, neither can effectively identify drivers. Instead, consideration of physical influences, removing cycles that artificially inflate correlations, and controlling the effects of other parameters gives a clearer picture of which are most influential in this system.

Introduction

Geostationary/geosynchronous orbit (GEO) is highly populated with active satellites (http://www.unoosa.org/oosa/osoindex/) that can experience damaging surface charging due to high energy electrons present in the radiation belts (e.g., Choi et al., 2011; Koons et al., 2000; Lam et al., 2012; Loto'aniu et al., 2015; Matéo‐Vélez et al., 2018). These and other studies suggest that surface charging is a function of factors in the space environment, including solar and geomagnetic activity, electron and ion flux magnitudes, and particle energy spectrum hardness. Charging events may also be more likely when the satellite is in the Earth's shadow (eclipse). Surface charging events often occur when there are increased electron fluxes at 10–50 keV (kilo electronvolt), and <100 keV electrons may be more responsible for the most rapid surface charging events than electrons at higher energies (Matéo‐Vélez et al., 2018; Thomsen et al., 2013). The abundance of these electrons fluctuates on time scales of minutes and also shows high spatial variability over the magnetosphere. For this reason, daily/orbit averaging misses much of the behavior of these electrons. Moreover, even moderate storms are not necessary for electron enhancements in this energy range, with many surface charging events detected during low to moderate substorm activity and no direct dependence on substorm strength (Ganushkina et al., 2021; Matéo‐Vélez et al., 2018).

A better understanding of keV electron flux behavior is needed, including details of how fluxes are driven and by what parameters. While a prediction model may hint at the drivers and mechanisms, no matter how well it may forecast, it is not a valid tool for effectively testing hypotheses about physical drivers. Hypothesis testing is best done with statistical tools developed specifically for this. Regression is one such tool, with multiple regression being the more appropriate test when multiple drivers are being considered. (Confusingly, the least squares method used in regression can also be used to create prediction models, but this should not be confused with hypothesis testing). The autoregressive moving average transfer functions (ARMAX) method, which we discuss below, is a refinement of regression that allows the modeling of time series behavior before the testing of input parameters. This will reduce possible spurious correlations that can occur if both dependent and independent variable time series cycle or trend simultaneously.

MeV (mega electronvolt) electron fluxes at GEO have been more extensively studied and may show high overall correlations with solar wind parameters when daily averaged (e.g., Blake et al., 1997), although the hourly response may be much lower (Simms et al., 2022). Solar wind speed is often cited as the most important driver (Li et al., 2001; Paulikas & Blake, 1979), although the relationship is complex (Reeves et al., 2011) and, for example, Lyatsky and Khazanov (2008) and Balikhin et al. (2011) have shown that the solar wind density is more associated with MeV electron variations. The solar wind electric field (E ) also shows an association with MeV flux (Baker et al., 2019) (but see Pulkkinen et al. (2016) for a discussion of how well the solar E maps to the magnetospheric convective currents). However, the direct influence of many solar wind drivers on even MeV electron flux is still unclear, both because much of the solar wind influence may not be direct but instead mediated by waves and electron injections following substorms (e.g., Simms, Engebretson, Clilverd, Rodger, Lessard, et al., 2018), and because simple correlations of solar wind parameters with electrons may be inflated by common cycles and trends if these commonalities are not removed via such methods as differencing transformation or ARMAX modeling (Simms et al., 2022). For keV electrons, there are even fewer simple answers as to which of the solar wind parameters drive their variations.

Fluxes of low energy electrons have been modeled with a first principle kinetic approach in several ring current simulations (e.g., Chen et al., 2015; Fok et al., 2014; Ganushkina et al., 2014; Harel et al., 1981; Jordanova et al., 2016). These models are driven by different sets of solar wind, Interplanetary Magnetic Field (IMF) parameters and geomagnetic indices but the drivers are predetermined. The first principle models cannot define the driving parameters themselves.

Empirical models can determine correlates of electron flux energies from eVs to several MeVs using a variety of fitting techniques. Among them, (a) one of the earliest models, the NASA (National Aeronautics and Space Administration) radiation belt models for electrons such as AE8 traditionally used to specify the average charged particle flux for space missions (Vette, 1991), (b) the improved AE9/SPM models (Ginet et al., 2013) derived from measurements made over an extended period of time by particle detectors and dosimeters on board many satellites in a variety of orbits (see table 3 in Ginet et al. (2013)), (c) a Particle ONERA (Office National d’Etudes et de Recherches Aérospatiales/French Aeronautics and Space Research Center)‐LANL Electron (POLE) model (Boscher et al., 2003) of energetic electron flux developed using 25 years of LANL data with input as the year in the solar cycle, (d) the extended POLE model known as the new International Geostationary Electron model (IGE‐2006) (Sicard‐Piet et al., 2008) created by adding the data from the Japanese spacecraft Data Relay Test Satellite (DRTS), and (e) the electrons model (Roeder et al., 2005) based on Polar HYDRA (Hot Plasma Analyzer) data proving the average flux as a function of the position in the Earth's magnetosphere. The models above were not parameterized on geomagnetic conditions and did not capture the magnetic local time (MLT) dependence and variations on time scales of less than a day.

The Kp (Planetarische Kennziffer) index, a simple 0–9 index as compared to the more complex variations of solar wind and IMF parameters, has been used to organize keV electron fluxes (e.g., Korth et al., 1999). Using LANL satellites data in the range from 1 eV to 40 keV at GEO, Denton et al. (2015); Denton et al. (2016) developed a model which predicts electron flux values based on energy and local time for given values of the 3‐hr Kp‐index and ‐V B (the electric field of the solar wind, where V is the solar wind speed, B is the z‐component of the IMF), under the assumption that both Kp and the solar wind electric field are correlated with magnetosphere activity, for example, for Kp (Freeman, 1974; Thomsen, 2004); for ‐V B (Akasofu, 1964; Burton et al., 1975). The Kp version of the model also provides flux values for given values of the daily F10.7 index (solar radio flux at 10.7 cm). However, while the Kp index may correlate well with flux (at least in daily averaged data), it is neither the best predictive parameter, nor what we would consider to be a physical driver of electron flux variations. Kp, as it Earth‐based (measured at ground magnetometers), may not represent conditions in the magnetosphere well. It is most likely a proxy measure, representing a combination of both relevant and non‐relevant correlated factors, which tells us little about which specific processes drive flux. While the ease of obtaining it might offset this drawback in prediction models, it may be nearly useless in models seeking instead to explain what drives electrons. Its 3 hr time cadence may also make it unsuitable even for prediction models, given that electron fluxes fluctuate much more rapidly. The ‐V B measure could be an improvement over Kp as it can be obtained at hourly or faster cadence and each is a specific physical parameter rather than a possible conglomeration of generalized response (as the Kp is). However, this measure, alone, only accounts for two possible driving parameters.

Several studies have examined the response of geosynchronous keV electron flux measured at LANL satellites to solar wind parameters. For example, Shi et al. (2009) found electron flux increases due to solar wind dynamic pressure enhancements and Li et al. (2005) and Kellerman and Shprits (2012) concluded that higher solar wind speed results in higher electron fluxes. Hartley et al. (2014) have found an effect of solar wind speed on the 30–600 keV electron density, temperature and energy density from the MAGED (MAGnetospheric Electron Detector) instrument onboard Geostationary Operational Environmental Satellites (GOES) 13–15.

Sillanpää et al. (2017), using 5 years of GOES 13 MAGED electron flux data, fit an empirical model using both solar wind and IMF B to predict electron fluxes at 40, 75 and 150 keV energies, after concluding that the other two IMF components and solar wind density, temperature, and pressure were of less importance. This is in line with earlier studies (e.g., Ganushkina et al., 2019; Kellerman & Shprits, 2012; Li et al., 2005). The effects of multiplicative combinations of parameters such as ‐V B (Denton et al., 2016) were not studied and it is possible that not a single parameter but the combined effect of multiple driving parameters that result in the observed fast variations of the keV electrons.

Ganushkina et al. (2021) discovered that the AE/AL (Auroral Electrojet/Auroral Lower) indices, together with solar wind speed, provide a better model of the severe environments related to surface charging of satellites by keV electrons measured by LANL (1990–2005) than do IMFB , Kp, and solar wind number density. Based on integral electron fluxes, among 400 events of worst‐case severe environments (categorized based on four criteria (Matéo‐Vélez et al., 2018) of the solar wind and IMF parameters and geomagnetic indices), 100 were in one criterion based on the measured spacecraft potential and 300 in the other 3 criteria based on these electron flux measurements.

In recent years, multivariate approaches have been explored to refine and complement physical and single variable empirical models, and to determine the main driving parameters of keV electrons. Some techniques used for predictions of mainly MeV radiation belt electrons include linear prediction filters (e.g., Baker et al., 1990; Castillo Tibocha et al., 2021; Rigler et al., 2004), dynamic linear models (e.g., Osthus et al., 2014), conditional mutual information (Wing et al., 2022), multiple regression (e.g., Sakaguchi et al., 2013; Simms et al., 2014, 2016, 2018a, 2018b; Stepanov et al., 2021), neural networks (e.g., Freeman et al., 1998; Koons & Gorney, 1991; Ling et al., 2010; Simms & Engebretson, 2020), and Nonlinear AutoRegressive Moving Average with eXogenous (NARMAX) inputs (e.g., Balikhin et al., 2011, 2016; Boynton et al., 2015; Boynton et al., 2016).

GOES 13–15 40 keV electron flux data were used by Boynton et al. (2019) to develop a model of time series of the electron flux for each of 24 MLTs employing NARMAX methodology. They found that the IMF factor, a combination of IMF B and B component (Balikhin et al., 2010; Boynton et al., 2011), B (t) = B (t) sin6(θ(t)/2), where and θ = tan−1(B (t)/B (t)), controls most of the output variance. Another important variable was determined to be the solar wind velocity. The square root of the solar wind pressure and solar wind density were also chosen by the algorithm but their contributions are small. Boynton et al. (2019) stressed that the time resolution of the parameters used in the model development influences the importance of these parameters. For comparison, the earlier study by Boynton et al. (2013), in which daily averaged 10–100 keV electron fluxes measured at LANL satellites were used, the role of southward IMF was found to be insignificant.

In the present study, we test the influence of several possible drivers of low energy electron flux (40–150 keV) observed by GOES 13 and GOES 16 satellites: solar wind velocity (V), number density (N), pressure (P), and the electric field (E ), IMF B , and substorms (as measured by the AE index). We use ARMAX (autoregressive moving average transfer function) models both to measure the cumulative effects and to remove common cycles and trends that may inflate correlations between variables (Simms et al., 2022). These parameters may act in combination, with influence accumulating over time. It is also possible that some variables may not influence electron flux directly but indirectly via other parameters. For the latter case, we develop subset models showing postulated direct and indirect effects.

Regression can be a powerful tool for testing which drivers could have a possible controlling influence on electron flux levels. However, regression on time series data, because it often violates the assumption of uncorrelated errors, can result in highly inflated hypothesis test statistics, giving the impression that certain factors may be strong drivers of flux when they are only cycling or trending in common (Simms et al., 2022). While this may not be a problem if we are using a regression model to forecast flux, it will invalidate the hypothesis tests that allow us to determine if solar wind, IMF, and substorm factors are meaningfully correlated with flux. We may also find that using more of the information present in the data (i.e., the time behavior) results in more accurate predictions.

There are several approaches to modeling the periodic behavior of a time series. We will do so with autoregressive (AR) and moving average (MA) terms (Hyndman & Athanasopoulos, 2018; Pankratz, 1991). When chosen well, these reduce the autocorrelation in the errors of the model and fully describe the cycling behavior of the series. With this behavior effectively removed (by the introduction of these terms) the remaining variability in the data can be tested for its response to external factors (the independent variables). This last step results in a transfer function model (X), giving the acronym ARMAX. The ARMA terms therefore partition out the part of the variance due to common cycling. What we will be left with is the actual relationship between the predictor and flux.

Further assumptions of regression models, in general, are a linear relationship between response and predictor variables, homoscedasticity (equality of variances in the residual errors), and normally distributed errors. To achieve linearity, we take the logs of the variables (excepting those with both positive and negative values). This not only allows the use of the linear model technique (regression) to describe what may be a nonlinear process, it can also normalize the residuals and reduce the heteroscedasticity (non‐equality of variances) of the residual error (Neter et al., 1990).

Other studies have used ARMAX modeling to predict higher energy electron fluxes in geostationary orbit, and these provide further information on this approach of describing the underlying cyclical behavior of flux with AR and MA terms (Balikhin et al., 2011, 2016; Simms & Engebretson, 2020; Simms et al., 2022). We note that nonlinearity was introduced into the models of Balikhin et al. (2011) with polynomial terms (square and cubic terms) instead of the logs we use here. Polynomial terms often result in a similar description of nonlinearity as merely taking logs of the variables. Additive polynomial terms may be a better choice if the slope of the relationship cannot be made constant with the log transformation, but there are potential disadvantages. The exponent of a polynomial term is chosen by the researcher, while linearizing the relationship using a log transformation allows the model itself to choose the exponent of the power function. The polynomial approach may also not have the same ability to fix problems of heteroscedasticity and non‐normality of the errors. However, although these appear to be significantly different approaches, most relationships in space weather data are such that either approach will give a similar description of the data.

We also note that ARMAX models may sometimes be called ARIMAX models, with the additional I conveying that the data is differenced at some time step n with a y  − y transformation. However, as we did not find it necessary to difference the present data set for the full models, ARMAX without the additional I is the more descriptive acronym.

In this study, we extend this previous work by using the ARMAX technique to determine the most influential drivers of lower energy electron flux behavior. While previous studies (e.g., Balikhin et al., 2011) may choose an optimal, parsimonious set of predictors that describe the variance in the data set (e.g., through the Error Reduction Ratio technique), using least squares regression (as applied to an ARMA model) we are able to show the statistically significant, relative contributions of each parameter rather than reducing the model to only highlight the most essential variables. In other words, we are able to test for the inutility of certain parameters in describing flux, rather than just choosing those parameters that have the strongest correlation. This provides more information on the additive influence of parameters, even if the influence of some is not as strong as others. This results in a deeper understanding of the ensemble effects. We also explore a reduced model consisting of just those parameters we hypothesize are the direct physical drivers of flux: AE (as a measure of electron injections from substorms), pressure, and the solar wind electric field (E , or − VB ).

The description of the data is given in Section 2. Section 3 presents the results for drivers for 40–150 keV. The findings are discussed and the conclusions are drawn in Section 4.

Data for Defining the keV Electron Drivers

For electron fluxes, we use hourly averaged data from the geostationary GOES‐13 satellite. We analyze the measurements from the MAGED instrument consisting of the nine collimated solid state telescopes (e.g., Rowland & Weigel, 2012), each with a 30° full‐angle conical field of view. All nine telescopes measured the directional differential electron fluxes in units of cm−2 s−1 sr −1 keV −1. We use the fluxes in the first three energy channels where the fluxes are defined at the midpoints of the energy ranges, that is, at 40, 75, and 150 keV. We compute one omnidirectionally averaged flux (flight direction‐integrated differential electron fluxes) for each of the energies using pitch angles calculated from the onboard magnetometer data following the method presented in Sillanpää et al. (2017) and Ganushkina et al. (2019). The GOES‐13 MAGED data of electron fluxes and the data for the pitch angles of each telescope with 5‐min averaging are available at https://www.ngdc.noaa.gov/stp/satellite/goes/dataaccess.html.

The time interval of this study is 10 June 2013–6 August 2016. There were minimal data gaps of only several hours during these time periods. For the time‐dependent analyses (ARMAX models) these gaps were filled with linear interpolation between the existing observations.

Solar wind parameters (solar wind velocity V, number density N, pressure P, IMF B and B (including only the southward component of B ), and the electric field E ) and magnetic indices (Kp, AE and SymH) were obtained from OMNIWeb (https://omniweb.gsfc.nasa.gov/ow.html) with 1 hr resolution with data time‐shifted to the bow shock nose. We use an hourly ground ULF wave index (ULF) as a global ULF activity proxy reconstructed from 1‐min data from the world‐wide array of magnetic stations in the Northern hemisphere (data available at: https://doi.org/10.2205/ULF-index) (Kozyreva et al., 2007; Pilipenko et al., 2017).

Analyses based on the least squares regression methodology assume that the relationship between predictor and response variables be linear, with the residual errors (that variance i.e., unexplained by the model) being random, normally distributed, and with equal variance over the range of predicted values. This requirement applies even to such analyses as simple correlation. However, the relationship between flux and predictor parameters is often nonlinear and inspection of the residual errors of these analyses performed on non‐transformed data shows this nonlinearity, as well as non‐normality and an inequality of variances at different levels of the predictors. Fortunately, these problems can usually be fixed by taking the log of at least electron flux, with further improvements obtained by taking the log of transformable predictor variables as well. We therefore take log 10 of all variables ≥0. Variables containing zero values which cannot be logged without creating missing values (i.e., Kp) were transformed by adding 1 to all values before the log transformation. B and E , as they have both positive and negative values, were not logged. Examination of residual plots of the ARMAX models (not shown) showed that this transformation fixed all three problems.

Because the dependent variable (electron flux) is log‐transformed, this results in nonlinear models between flux and all the variables, a power function relationship for those predictor variables that are also log‐transformed, and an exponential function relationship for those predictor variables that are not logged.

Subsequent to the log transformation, all variables were standardized by subtracting that series mean and dividing by its standard deviation. This creates unitless variables (Z‐scores) for which regression coefficients (slopes) can be directly compared. Although it makes no difference to the outcome of the correlations, we also used the Z‐scores for the correlation analysis for consistency. We note, however, that neither the log nor the Z‐score transformation reduces either the serial autocorrelation or common cycles seen in these time series datasets. This autocorrelation inflates the simple correlations and must be further dealt with by describing/removing the autocorrelation and common trends and cycles via the introduction of AR and MA terms and/or differencing, as described below in Section 3.2 (Granger & Newbold, 1974; Simms et al., 2022).

ARMAX models were developed in IBM SPSS Statistics (formerly known as the Statistical Package for the Social Sciences), with additional statistical analysis in MATLAB.

Drivers of 40–150 keV Electrons at Geostationary Orbit

Cross Correlations of Electron Fluxes With Solar Wind and IMF Parameters and Geomagnetic Indices

Simple cross correlations of hourly measured parameters (Figure 1) show values near 0.5 for some parameters, most notably and in keeping with previous studies, V, ULF, and AE (e.g. (Hartley et al., 2014; Kellerman & Shprits, 2012; Li et al., 2005; Simms et al., 2014),). Positive correlations are shown with solid lines, negative with dashed red lines. Correlations are performed between electron flux and individual parameters from each hour (0–48 hr) before the flux measurement. At higher electron energies, the AE and ULF correlations are lower, with peak correlations at earlier times. The correlation with V may be somewhat higher, but there is also a tendency for its peak correlation with electron flux to occur earlier at higher energies. The correlation of flux with N is less than that with V, although it does become more prominent at 150 keV, if negative.

Figure 1

Cross correlations between Geostationary Operational Environmental Satellites electron flux and possible drivers (hourly averages) for (a) 40 keV, (b) 75 keV, (c) 150 keV. Solid lines are positive correlations; dashed lines are negative correlations. Note that most correlations are <0.5 in magnitude.

B and B correlations with flux are similar to each other. There appears to be no particular advantage to using the B parameter over B . (Correlations of By with flux at all hours and energies were <0.05. This parameter was therefore not studied further). The negative correlations of B and SymH with flux are as expected, as each of these parameters are measured on a negative scale indicating increasing strength at more negative values. While the B strength shows less association with flux, SymH and Kp show similar patterns of correlation to each other, likely because both are generalized measures of disturbance at ground magnetometers. These parameters also show an increased correlation at earlier time steps at higher flux energy.

P and E are somewhat different from the other variables in that they are mathematical combinations of other measured parameters (V 2 and N, and V and B in the cases of P and E , respectively), but, at the same time, they may have more physical interpretability. That the P‐flux correlation is similar to that of the flux correlation with V or N can be seen where the P correlation drops off in a manner similar to the N correlation, albeit, with some tempering of this decrease as the V correlation rises at the same point in time. The E ‐flux correlation follows the pattern of the B ‐flux correlation nearly exactly.

Interpretation Problems With Simple Correlations: Poorly Defined Variables, Autocorrelation, and Spurious Correlations

Most of the parameters of Figure 1 show more association with flux in the few hours just prior to a flux measurement at the lower energies, but with maximum correlations at the higher energies occurring further back in time. However, it is difficult to interpret a single peak or even a rise in correlation at a given hour as a physical process that happens at that particular time, given that all these parameters are strongly auto correlated in time. A variable strongly correlated with itself in previous time steps will show a similar correlation with flux at every one of those time steps, making it impossible to determine the exact time of physical action from simple correlation analysis.

Another difficulty with simple correlation analysis is that correlations between predictor variables may distort the apparent association between a predictor and flux by confounding the true relationship. The well known correlation between V and N, for example, even if it is negative, will result in both predictors showing a correlation with flux, even if only one of them has an actual association. Besides this, any co‐cycling variables will show a strong correlation even when there is no association other than a similar response to time. This is a particular difficulty in space weather data where both diurnal cycles and longer cycles are common.

Although we find reasonable correlations of SymH and Kp with flux, to justify including these in a model attempting to find the physical drivers of flux, there must be some basis for thinking there is a physical connection between these particular indices and electrons. While Kp, derived from midlatitude stations, may be sensitive to variations at the inner edge of the electron plasma sheet (Freeman, 1974; Thomsen, 2004), there is no guarantee that this is all or even most of what Kp measures. As the measure itself is merely the maximum geomagnetic disturbance recorded in a 3 hr period, it may not be specific to that particular area of the magnetosphere, nor temporally fine tuned enough to be of much use. The discrete nature of the index values would also work to obscure much of the information it could carry. That there are high correlations between electron flux and Kp (see Figure 1) is not an argument in favor of its necessary inclusion in a meaningful physical model, but may more likely only indicate that Kp is a proxy that represents a large number of processes that we would, instead, prefer to know the effects of individually. In addition, as parameters that are averaged over longer periods of time tend to show higher statistical correlations without any meaningful increase in association (Simms et al., 2022), this alone could explain the Kp, at a 3 hr cadence, having a higher correlation with flux than that of many other parameters. SymH may be an indirect measure of the free energy available for local wave acceleration of keV electrons up to MeV energies, but is perhaps more representative of inner magnetospheric plasma pressure, about 12% of which is keV electron pressure (Kumar et al., 2021). SymH may be worth testing as a representation of these processes, but the applicability to electron flux in the outer radiation belts appears weak. While the AE index can be interpreted as a measure of the substorm activity that may result in electron injections, we do not have the a similarly meaningful physical interpretation of Kp and SymH other than that they measure the overall level of disturbance in the magnetosphere. But if “disturbance” is a meaningful concept, it is more accurately measured by such parameters as V, B , etc., which also have a physical meaning in the system. In previous work it has also been found that indices from magnetometers tend to correlate highly with each other, meaning that it may only be useful, or possible, to include one index in a multivariate analysis without reaching problematic levels of multicollinearity that make it impossible to determine which variables are most associated with flux (Simms et al., 2016). Therefore, we need to use care in deciding which index to use and not include every one possible. Instead, we should settle on the one that best describes the physical processes we suspect are occurring.

However, these arguments are somewhat moot. If we do include all three indices (Kp, SymH, and AE in a full regression ARMAX model (see below; Table 1), Kp and SymH are not strong candidates, as their influence can be up to an order of magnitude below that of AE. Although Kp and SymH have high simple correlations with flux, and even if we were to believe they represented physical drivers, when variables are tested simultaneously, these two indices do not perform well. In the subset models, we therefore use the AE index both because it is representative of substorm activity and because it is a stronger correlate, at least at 40 keV. In future work, if we planned to create prediction models only, but not to identify physical drivers, this restriction would not apply and all three indices could be included (with the caveat that this did not result in overfitting and, therefore, poor predictive ability).

Table 1

Autoregressive Moving Average Transfer Function Standardized Regression Coefficients of the Full Models (One for Each Electron Energy) Including All Variables Except B

	40 keV	75 keV	150 keV
Intercept	−0.057 a	−0.054 a	−0.03 a
AR1	0.836*	0.845*	0.855*
MA1	0.204*	0.207*	0.069*
MA2	0.302*	0.217*	0.202*
DailyAR1	0.999*	1.000*	1.000*
DailyMA1	0.986*	0.993*	0.994*
V lag 1 hr	0.632n.s.	0.888*	−0.196*
Decay	0.270	0.822	−0.147
N lag 1 hr	1.087*	1.358*	−0.087*
Decay	0.126	0.811	0.854
Bz lag 1 hr	0.265*	0.386*	0.306*
Decay	0.283	0.429	0.673
Kp lag 1 hr	0.017 a	0.041*	0.023*
Decay	−0.563	0.937	0.967
SymH lag 1 hr	−0.028*	−0.004*	0.056*
Decay	0.726	0.975	−0.447
AE lag 1 hr	0.170*	0.131*	0.019*
lag 2 hr	–	0.050*	0.062*
Decay	0.379*	−0.055	0.551
ULF lag 1 hr	0.021*	0.001 a	0.003 a
Decay	0.959	−0.988	0.984
P lag 1 hr	−0.992*	−1.274*	0.035 a
Decay	0.177	0.813	0.849
Ey lag 1 hr	0.257*	0.352*	0.263*
lag 2 hr	−0.131	−0.040*	–
Decay	−0.046	0.414	0.731
R 2	67.4%	69.2%	78.1%

Not Statistically Significant p‐value >0.10.

Note. *: Statistically Significant p‐value <0.05. n.s.: not statistically significant, p‐value > 0.10).

Although simple correlations can suggest possible drivers, further work must be done to elucidate these relationships. Below, in our ARMAX models, we address these issues, performing multivariate analyses to account for spurious simple correlations due to the confounding of variables, adding autoregressive (AR) and moving average error (MA) terms to account for serial autocorrelation and co‐cycling of variables, and choosing predictors that have a reasonable basis for some physical relationship with flux. In regards to this latter issue, we also choose 4 variables (AE, ULF, P, and E ) as possible direct physical drivers of flux (direct effects) and explore their relationship with the other solar wind and IMF parameters (indirect effects).

Additionally, below, we explore whether certain parameters are more correlated during geomagnetically disturbed periods and at different times of the day. For the former, we must use a differencing transformation (y  − y ) to reduce serial autocorrelation as, without a complete time series, we are unable to remove this with ARMA terms. To study varying influences by time of day, we add indicator variables to the ARMAX model to identify each hour (MLT: magnetic local time).

ARMAX Models

As noted in the previous section, simple cross correlations of time series variables may be highly inflated by common cycles and trends often seen in time series data (Granger & Newbold, 1974). These correlations may, therefore, not say anything useful about the relationship between variables. In addition, analyzing the effect of each predictor individually gives us no information about the relative importance of each, or the effect of each when the others are held constant. Multiple regression analysis would assess the strength of the relationship between each predictor with the effects of the other predictors eliminated. Additionally, as regression gives us the slope of the relationship between predictor and flux (the coefficients of the regression equation), there will be more information about the form of the relationship. We can further improve on a multiple regression model by introducing terms to specifically describe the cycling, trends, and autocorrelation that may be present in time series data. These terms may take the form of an autoregressive component (regressing on previous values of the dependent variable: an AR term), or a moving average component (regressing on the errors of the model at preceding time steps: an MA term). A difference term, which subtracts a previous value from each observation, may also be used to fit an overall trend, but we found this was not needed for this full set of hourly averaged flux data. For data that cycles “seasonally” (at a set time period) it may be helpful to also fit seasonal AR and MA terms (Hyndman & Athanasopoulos, 2018).

We fit ARMAX models, using AR and MA terms, along with “seasonal” (daily) AR and MA terms, to describe the cycling behavior of the dependent variable. We are then able to test input variables for their possible correlation separate from these common cycles. The “seasonality” we incorporate is the daily variation in flux seen as the observing satellite passes between drift shells due to the asymmetric dipole of the Earth's magnetic field. Typically, higher energy (MeV) electron flux data collected at geosynchronous orbit shows higher levels on the dayside where the field is compressed and lower flux levels on the night side where the fields are stretched (e.g., Boynton et al., 2019; O’Brien & McPherron, 2003). For keV electrons, fluxes are highest in the morning hours and lowest in the evening hours due to their trajectories and losses (e.g., Korth et al., 1999; Sillanpää et al., 2017).

As all variables were standardized by subtracting that series mean and dividing by its standard deviation, we are able to compare these unitless regression coefficients between variables. Note that these are not correlation coefficients, but slopes. A 1 unit increase in a predictor variable is thus associated with a certain increase in the dependent variable. Taking log 10 of those variables for which it made sense (i.e., not B , e.g., which has both positive and negative values) effectively creates a non‐linear model, despite how we are using the linear model technique of ARMAX regression.

For each electron flux energy (40, 75, and 150 keV), we fit an AR1, MA1,2, seasonal AR1, seasonal MA1 model. More specifically, each regression contained two flux autoregressive terms (from 1 hr previous and 24 hr previous) and the moving average of the errors of the model from 1,2, to 24 hr previous as predictors, in addition to the exogenous AE, Kp, SymH, ULF, and solar wind and IMF variables. The 24 hr AR and MA terms represent the “seasonality” terms that model the diurnal fluctuations in flux due to the movement of the satellite through field lines (in other words, the “seasons” are days (Table 1). This reduced all terms of the partial autocorrelation function (PACF) to non‐statistically significant levels.

Full ARMAX Model Including All Variables

V, N, IMF B , AE, ULF, P, E , Kp, and SymH were first entered as numerator (influence) terms at 1 and 2 hr delays, with a denominator (decay) term at 1 hr (Table 1). Influence terms with p‐value >0.10 were dropped from the model. The p‐value is the probability that the null hypothesis of no association is true. A p‐value <0.05 is generally considered to be statistically significant, or, put another way, that the null hypothesis of no association has been rejected. (The calculation of the p‐value is dependent on the presumed sampling distribution of the test statistic under the null hypothesis, which is itself dependent on the standard error of the estimate. Simply put, dividing the parameter (coefficient) by its standard error gives a t‐statistic which can be compared to standard tables giving the probability of that value given the relevant sample sizes. While it can be calculated by hand, it is generally left to the statistical package to do the calculation and table look up (Neter et al., 1990). Due to this p‐value restriction, not all influence and decay terms are retained, however, at least one influence and the decay term are retained for each predictor, even if statistical significance fell above a p‐value >0.10, in order to describe the relative influence of each term. (The constant term is not significantly different from zero because all variables were standardized and therefore centered around zero. However, we retain it for the small amount of explanatory value it adds to the model). We report standardized regression coefficients which describe the slope of the relationship between predictor and response variables on a standard (unitless) scale. Due to this standardization we are able to directly compare the influences of each predictor with all the others. (These are slopes, not correlation coefficients, so are not constrained to lie between −1 and 1).

The R 2, or coefficient of determination, measures the percent of variation in the data that is explained by the model. (Note that the R 2 is mathematically equivalent to the prediction efficiency used by some other authors when applied to a training data set). In simple correlation, the R 2 is equivalent to the square of the correlation coefficient (r 2). The highest simple correlations (e.g., AE and V of Figure 1) around r = 0.5 would therefore have an R 2 of 25%, explaining 25% of the variation in the data. Thus, the multiple regression ARMAX models which use both ARMA terms and more than one predictor variable, explain more of the variation than any of the simple correlations. Much of the increase in R 2 is due to the introduction of the ARMA terms, but the ARMAX models do also tell us which independent variables are most important and how they compare in influence with each other. This addition of predictor variables would also allow the ARMAX model to be used for prediction. If there are no exogenous (independent) variables in the model, predictions would quickly revert to the mean value of zero, the constant of the ARMAX equation.

The predictor coefficients can be represented with an empirical prediction equation (Equation 1). For the 40 keV electrons:

40 keV electron flux at time t is predicted by the other variables at previous times steps (t − 1, etc), the model predicted value of flux at t − 1 and t − 24 (“daily”), and the error between model and observation (ɛ) at t − 1, t − 2, and t − 24. Each influence term is represented in a numerator, with decay terms in the denominator. For clarity, we do not label the variables that have been logged (flux, V, N, Kp, AE, ULF, and P) in the empirical prediction equation, however, due to this transformation, this is effectively a non linear model in the terms for which we have taken logs. For example, the partial influence coefficient for V of 0.632 would be interpreted as a 1% change in the Z score of V resulting in a 0.632% change in the Z score of electron flux. However, given the number of differences between our model and that of previous models (including more variables, using Z scores, including decay terms) a direct comparison of these coefficients to other studies using somewhat similar techniques is not particularly meaningful. Instead, it makes more sense to compare which predictor variables are most influential and not the details of exactly how much flux changes in response to a unit change in a predictor.

The influence (numerator) and decay (denominator) terms of Equation 1 give us the tools to calculate the cumulative effects of each input variable. An influence that appears at t‐1 dissipates at a rate given by the decay term. Thus, although there may only be 1 hr in which a variable input appears, the exponential decay over time means influence may spread from previous time periods. The influence at a given forward time step from some time step (t) in the past will be that influence × (1 − decayfactor) . Graphically, this results in a time delay of influence that appears similar to a cross correlation, however, the transfer function gives regression coefficients (i.e., slopes), not correlations. While a correlation can be interpreted as the strength of a relationship between two variables, a regression coefficient can be interpreted as the magnitude of the impact of one variable on another. We use the predictor coefficients of Table 1 to create the cumulative influence bar charts of Figure 2. It should be remembered that these regression coefficients represent the influence of each variable with the others held constant, unlike the simple correlations of Figure 1. Each panel of this figure shows the response of an electron energy (40, 75, and 150 keV) to the influence of each of the 9 exogenous variables when the other 8 predictor variables are held constant. The influence of each begins from the hour previous to the flux measurement. The decay term describes the fall off in influence over time.

Figure 2

Cumulative effects of all possible drivers of electron flux. For each flux energy, variables are entered simultaneously into an autoregressive moving average transfer function regression model as a predictor at a delay of 1 and 2 hr. Only statistically significant time steps are retained, along with a decay factor. Standardized regression coefficients may be compared within each model (a). 40 keV, (b). 75 keV, (c). 150 keV to determine the relative influence of each variable on flux. Note that each row has the same scale, but scales vary between rows, in order to compare more effectively between the strongest associations (V, N, and P) and between the indices (AE, Kp, and SymH) and other variables with lower influence (ULF, B , and E ).

These ARMAX models incorporating all nine possible predictors show little influence of Kp and SymH. AE has the highest influence of the geomagnetic indices, but it is weaker than the strong and lasting effects of V, N, and P, particularly at 75 keV. The V, N, and P influences are superficially similar to those seen in the simple cross correlations (Figure 1) but the sign of influence of N and P have switched. B and E also superficially show the same influence as in the cross correlations, but, again, the sign of influence of B is switched.

What are we to make of these losses of influence (particularly Kp and SymH) and the changes in sign? It is obvious that simple correlations are highly unreliable as each parameter is highly correlated with all the other parameters of interest, and because any one of them may show a spurious correlation with electron flux due to common cycling behavior.

Second, geomagnetic indices (particularly Kp and SymH) do not even appear to influence electron flux when other variables are present. In this full model, Kp and SymH have little influence. However, even if they were the most “influential” parameters in these models, for the reasons mentioned above would we be justified in calling them drivers? Or are they merely correlated proxies? Is SymH a predictor variable at all? Or just another measure of our response variable, the electron flux? These questions can only be answered from a consideration of what information these indices actually contain. As we have discussed above, while Kp and SymH may roughly represent disturbance in the magnetosphere, we don't know exactly which processes and how much of each process they might represent. AE is a different case. First, AE does show more influence than the other two indices, and second, we know that this index measures substorm activity which can lead to electron injection. For this latter reason, we will retain AE in further models.

Both P and N act more as a pointed shock to the system with less long term influence, however, the opposite sign of these two predictors, at similar magnitudes, suggests that there is some degree of multicollinearity occurring between these two. This is not surprising, as P, partially calculated from N, is highly correlated with N and the amount of information about the influence of each on flux is almost wholly contained in the other. Unfortunately, this can result in a pattern of presumed “influence” (as seen here) that reflects a competition for explanatory power rather than actual opposing effects on flux, and the inclusion of both in the model is misleading. It is unclear why their influences appear larger at 75 keV than at 40 or 150 keV, but this could simply be the result of only a small difference in flux related to the N and P variables being amplified by the competition between the two. It is unlikely to represent anything tangible and only demonstrates that coefficients of highly correlated variables in the same model are not trustworthy.

B and E have more modest influences on flux. Despite the high ULF‐flux correlation seen in the simple correlations, the ULF influence on flux is very low. This is likely due to two factors. First, when other variables are included in the model any proxy correlation ULF may have represented is removed from the ULF influence. Second, the high simple correlation may be simply due to this ULF index and satellite‐measured flux both showing a diurnal cycle. When this cycling is removed (via the AR and MA terms) the correlation between these variables disappears (Simms et al., 2022). (The occasional oscillating pattern of influence in several of the variables is the result of a negative decay term found by the regression. It is often unclear whether this has any real physical meaning).

As these are standardized regression coefficients, we can calculate the impact of a predictor on flux. For example, as we are using standardized coefficients, a 1 standard deviation increase in log 10(AE) 1 hr previous would result in 0.17 standard deviation increase in log 10(40keV flux), holding all the other predictors constant.

Choice of Variables

Pressure (P) and number density (N) are difficult to incorporate into a model simultaneously. As pressure is the product of the V 2 and N, the strong correlation between pressure and N can lead to unexpected and puzzling behavior. In the models of Figure 2 and Table 1, there is a strong initial influence of P, and an opposing strong influence of N in the same time period. As we know that P and N are highly correlated with each other, it is difficult to interpret this as each having a strong, opposing, and, most importantly, independent influence. It is more likely that these opposing effects are merely the result of the two terms acting counter to each other in an effort to explain the same small bit of variation. The same is true of E with V and IMF B , as E is the product of V and B . In fact, it likely makes little sense to draw firm conclusions about the physical drivers based on this full model. A more plausible model could be achieved by dropping one of either P and N, and one of E and IMF B . For example, dropping the two derived parameters (E and P) would allow us to more accurately see the effects of V, N, and B .

However, we may be able to do better by separating out just those parameters we believe could be influencing flux directly. These direct parameters would be AE (as a measure of substorms which inject electrons), ULF (waves in this frequency are thought to drive electrons to higher energies), E (with the solar wind electric field plausibly having some influence on electron behavior), and pressure (which could influence flux levels through acceleration, through magnetopause shadowing, and by compression of the magnetosphere at the altitude of the satellite, bringing the satellite into higher drift shells with lower electron density). We note that there is little theoretical work to suggest either N or V alone drive flux. Their physical action, instead,is thought to derive through pressure. We similarly assume that the physical action of B is likely through the electric field rather than simply the B itself. We fit a reduced ARMAX model using only these 4 parameters. The coefficients of this reduced model are presented in Table 2.

Table 2

Autoregressive Moving Average Transfer Function Standardized Regression Coefficients of the Three Reduced Models Using AE, ULF, P, and Ey as Predictors

	Log 40 keV flux	Log 75 keV flux	Log 150 keV flux
Constant	−0.090n.s.	−0.093n.s.	−0.056n.s.
AR1	0.825*	0.843*	0.86*
MA1	0.197*	0.195*	0.055*
MA2	0.293*	0.212*	0.201*
Daily AR1	0.998*	0.998*	0.999*
Daily MA1	0.981*	0.987*	0.993*
Log(AE) 1h lag	0.216*	0.130*	0.004n.s.
2h lag	0.154*	0.091*
Decay 1h	0.882	0.542	0.053
Decay 2 hr		0.349
Log(ULF) 1h lag	0.017*	0.021*	0.03*
Decay 1h	0.965	0.969	0.97
Log(P) 1h lag	−0.025*	−0.039*	−0.055*
Decay 1h	0.717	0.728	0.801
Ey 1h lag	−0.018*	−0.014*	−0.03*
2h lag	—	−0.022	—
Decay 1h	−0.763	−0.381	0.412
R 2	67.10%	68.50%	76.90%

Note. *; statistically significant, p‐value <0.05, n.s.; not statistically significant

From the coefficients of this table, we once again calculate the cumulative effects of each variable on flux (Figure 3). At 40 keV (Figure 3a), this simpler model of the presumed direct effects alone shows a strong effect of AE, peaking at 2 hr before the flux and with influence over many hours. Pressure, E , and ULF, while still statistically significant effects, are much lower in magnitude. The effect of pressure is negative, presumably as most of its effect is due to the compression of the magnetosphere which positions the satellite into a less populated drift shell and to magnetopause shadowing. The small E association cycles between positive and negative. A similar pattern is seen for the 75 keV electrons (Figure 3b), although the AE influence is slightly lower and the P and ULF effects somewhat stronger. The 150 keV electrons (Figure 3c) show a much lower response to AE, and, again, a somewhat stronger response to P and ULF.

Figure 3

Cumulative effects of the possible direct drivers of electron flux. For each flux energy, AE, P, E , and ULF are simultaneously entered into an autoregressive moving average transfer function regression model as predictors at 1 and 2 hr, but only significant time steps are retained, along with a decay factor. Standardized partial regression coefficients may be compared within each model to determine the relative influence of each variable on flux: (a). 40 keV, (b). 75 keV, (c). 150 keV.

But what of the strong influence of V we saw in the full model of Figure 2? Although our direct effects model (of Figure 3) may make more physical sense, we still would like to understand the correlation of V with flux. We can do this by using the other, indirect parameters to predict our set of more physically interpretable variables, decomposing each correlation into components. In other words, we can use N, V, and IMF B to predict AE, ULF, P, and E , which we subsequently use to predict flux.

To accomplish this, we presume a causal model (Figure 4) and run a series of regressions to determine the coefficients of the paths. In this figure, we present the standardized regression coefficients obtained by predicting 40 keV flux from AE, ULF, P, and E . We then predict both AE and ULF using P, E , N, V, and IMF B from 1 hr previous. (These models are not shown explicitly as the input parameter coefficients are all that we need here, but these are simply the exogenous coefficients from an ARMAX model also incorporating AR and MA terms. For this particular model, we use only a lag 1 hr influence term and no decay term to simplify the effects of each input variable). Similarly, we show the exogenous variable coefficients for predicting P from N and V, and E from V and B , using N, V, and B , but from the same hour as P and E . (There are no paths from V to either ULF or AE because it was not a statistically significant direct influence on either). In this figure, green arrows run to and from AE, gold arrows to and from ULF, and blue arrows to and from P and E .

Figure 4

Postulated direct drivers of 40 keV Geostationary Operational Environmental Satellites electron flux (green arrows to and from AE, gold arrows to and from ULF, blue arrows to and from P and E ) may be influenced by solar wind and Interplanetary Magnetic Field parameters (V, N, and B ). Standardized coefficients of the influence of AE, ULF, P, and E on flux (from an the autoregressive moving average transfer function [ARMAX] model with predictors measured 1 hr before flux) are given. ULF and AE are postulated to be driven by P, E , V, N, and B (coefficients from ARMAX models with predictors measured 1 hr before). P and E , being mathematically dependent on N, V, and B , are predicted from ARMAX models with all variables measured at the same hour. Influences of V, N, and B on P and E are from the same hour. These paths break down the overall correlations into components, attributable to the various associations between variables. Only statistically significant links between variables are retained. As a consequence, there is no direct link from V to either ULF or AE.

These standardized regression coefficients from this series of regression models are known as path coefficients (Wright, 1934). The path coefficients can be multiplied (through connecting arrows, or paths), then summed to show the full cumulative effect of each of the indirect drivers (V, N, and B ) on the direct drivers (AE, ULF, P, and E ) and, subsequently, on flux.

The maximum direct effect of each variable is shown by arrows leading directly to flux. Simple correlations between the exogenous, or indirect, variables (N, V, and B ) are shown (in black curved arrows). This decomposition allows the correlation between a pair of variables to be broken down into direct effects, indirect effects, and spurious correlation due to associations between the exogenous variables. We are interested in the direct and indirect effects and will ignore spurious correlations due to the associations between N, V and B . For example, the direct effect of pressure on flux is represented by the arrow from pressure to flux (−0.04 coefficient). This is rather low, but to this we can add the indirect effect of pressure: the path from P through AE to flux (coefficients 0.52 and 0.25). This indirect effect of P via its influence on AE (which subsequently influences flux) is the product of the steps in the path: 0.52 × 0.25 = 0.13. The contribution of several indirect paths can be calculated by summing these products (Table 3). In the first column of this table, we show the direct effect of AE, ULF, P, and E on flux (coefficients on the arrows leading directly to flux of Figure 4). In the second column we show the results of the calculations for the indirect effects of each variable through AE, in the third column, these indirect effects through ULF, in the fourth, indirect effects through P, and in the fifth column, these indirect effects through E . (Details of example calculations are given in the footnote). The last column is the sum of the first 5 columns, showing the total influence of each variable, both through its direct influence (if any) and its indirect influence via other parameters.

Table 3

Calculating the Sum of Direct and Indirect Influences on 40 keV Flux

	Direct	Via AE	Via ULF	Via P	Via Ey	Sum
						Direct + Indirect Influence
AE	0.25	–	–	–	–	0.25
ULF	0.02	–	–	–	–	0.02
P	−0.04	0.13 a	0.018	–	–	0.11
Ey	−0.01	−0.055	−0.005	–	–	−0.070
N	–	−0.12	−0.015	0.12 b	–	−0.014
V	–	0	0	0.137	−0.0007	−0.024
Bz	–	−0.13	−0.0088	–	0.070	−0.071

As an example, the indirect path of P influence through AE = (effect of P on AE) × (effect of AE on flux) = 0.52 × 0.25 = 0.13, using coefficients from the paths in Figure 4.

The more complicated paths of N through P are summed: (N on P) × (P on flux) + (N on P) × (P on AE) × (AE on flux) + (N on P) × (P on ULF) × (ULF on flux) = 1.1 × (−0.04) + 1.1 × 0.52 × 0.25 + 1.1 × 0.88 × 0.02 = 0.12.

The result of these calculations are that we can now see a clearer picture of which variables are most influential on flux and through which processes that influence is mediated (given this particular, hypothesized, causal structure). Predictors not postulated to directly influence flux, such as V, still show an overall moderate degree of influence when paths connecting it indirectly to flux are considered (mainly, in this case, via P). However, N, which has a moderate (if negative) simple correlation with 40 keV flux, has less influence than V when all influences are added. N appears to drive several competing processes: reducing AE and ULF while simultaneously (through P) increasing flux. Thus, the lower correlation of N with flux is not an indication that it does not influence flux, but that it does so through several opposing processes that cancel out each other's effects in an overall correlation.

Certain parameters, such as ULF, which show a strong simple correlation with flux (Figure 1), are not influential. So why does the simple correlation appear so high in comparison? This is due to several factors which we have now accounted for: inflated correlations due to common cycles and trends (accounted for by the AR and MA terms of the ARMAX regression), correlations with confounding variables (now accounted for by the use of multivariate regression instead of single correlations), and the possibility that ULF over the short term (hourly, in this case) has little influence.

For parameters such as V and N, influence has been diminished by their relegation to indirect driver status in the path analysis. This is a choice made based on the hypothesis that neither is postulated to have the physical ability to directly drive electron flux. If there were reason to believe they did, these could be moved up the hierarchy in the path analysis, allowing them to have more influence in that correlational structure.

We can do these calculations for each of the electron energies, giving the summed influence of each parameter on flux (Table 4). AE appears only as a direct effect, and is thus comparable directly between electron energies, with the strongest effect at 40 keV (0.25) but a lower effect above this range (−0.001–0.15). The summed influence of P is generally larger and positive compared to its weak negative direct effect, particularly at 40 keV. The summed E effect is similar in magnitude to P. The summed effects of V, N, and B are all somewhat equal to each other, with somewhat more effect of V at 40 keV and a higher influence of N at 150 keV. For the most part, these three indirect drivers are negative in influence overall.

Table 4

Summed Direct and Indirect Influences on 40, 75, and 150 keV Flux

	a. AE	b. ULF	c. P	d. Ey	e. N	f. V	g. Bz
40 keV	0.25	0.020	0.11	−0.070	−0.014	−0.024	−0.071
75 keV	0.15	0.005	0.021	−0.052	−0.051	−0.049	−0.029
150 keV	−0.001	−0.008	−0.090	−0.023	−0.092	−0.065	0.027

MLT Dependence of 40–150 keV Electron Flux Response to , , , and

Electrons at geostationary orbit show different flux levels at different magnetic local times (MLT) (Boynton et al., 2019). With geostationary satellites, which orbit synchronously with MLT, it is unclear whether these are spatial or temporal variations, however, electron injection has been observed in the hours around local midnight (Birn et al., 1997; Thomsen et al., 2001). Using ARMAX models, we investigate not only whether flux differs at varying MLT, but also whether the identified drivers show different influences (i.e., a different coefficient slope) at each MLT. We do not subset the data into MLT bins and analyze them separately, but identify each MLT in the data set and calculate a different slope coefficient for each. This is done by creating a set of 23 indicator variables spanning the MLT hours: each is set to 1 for a different, particular MLT and 0 at all other times. The interaction term between each of these indicator variables and each predictor variable (obtained by multiplying each indicator variable by each predictor) gives the slope of the relationship between flux and predictor at each MLT (Neter et al., 1990). By not splitting the data set by MLT (i.e., by identifying MLT by indicator variables instead), we are able to analyze the data set as a continual ARMA process. We report these slopes (standardized regression coefficients) for each MLT (Figure 5).

Figure 5

Varying effects of AE, ULF, P, and E over magnetic local time. Each variable is entered into an autoregressive moving average transfer function regression model as a predictor at 1 hr (a). 40 keV, (b). 75 keV, (c). 150 keV.

At 40 and 75 keV, AE is the most influential (positive) parameter, but it is most effective over 3–11 MLT (40 keV) and 6–17 MLT (75 keV). Not only is the flux higher at these times (Boynton et al., 2019), but the effect of the strongest driver (AE) is also at its highest level.

The other direct drivers (ULF, P, and E ) are, as demonstrated above, less influential, but there are MLT differences in their effects. ULF has somewhat more positive effect at 19‐0 MLT on the 40 keV electrons. P shows a stronger negative effect over 16‐4 MLT, with the most effect being seen at 150 keV. E , at 40 and 75 keV, shows a positive effect over 23–8 MLT, with a negative effect over 9–22 MLT. The E switch in influence from positive to negative likely accounts for its overall lack of effect in the analyses above that are not broken down by MLT. Although less dramatic, the switch in ULF from positive to slightly negative or near zero also results in an overall lack of influence when MLT is not considered, even though ULF does show a modest positive influence at some times.

Disturbed Versus Quiet Response

To produce an ARMAX model, a continuous time series is needed. This means that disturbed and quiet periods must be combined in the same analysis. However, it may be that the flux response to each predictor varies depending on conditions. A simpler multiple regression model could be used to explore the response between quiet and disturbed periods, however, this can result in spurious correlations if variables are cycling together (e.g., a diurnal cycle) or show a common trend (Simms et al., 2022). A regression model that accounts for these co‐occurring cycles and trends can be produced by differencing the data: subtracting the previous value from each observation (y  − y ). This results in regression coefficients that describe the change in flux as predicted by the change in the independent variables, rather than in the original units, but tests of significant influence and comparisons of relative influence can still be made.

To pinpoint those periods when predictors may have differing influence on electron flux due to a change in geomagnetic conditions, we compare disturbed versus quiet conditions by assembling a subset of “disturbed” hours (a day before and a week following each Dst dip to −100 nT) and a “quiet” set (>2 weeks after a Dst dip below −30 nT). Within the generalized disturbance periods, we also create a subset (“recovery”) with the usually short period of the Dst drop, or storm main phase, removed. The period of time when the Dst is dropping is both short and of less interest to the question of what drives electron flux changes as most of the observable change in this period is due to the compression of the radiation belts below the altitude of the satellite. Following the storm main phase, Dst rises slowly to −30 nT and higher, the recovery and after recovery periods when electron fluxes increase. The rise in electron flux can take up to a week to occur following the main phase, particularly after the strong storms we use in this data subset (Simms et al., 2014). We first perform a multiple regression on the differenced data with AE, ULF, P, and Ey as predictors in order to compare their relative effects via the standardized regression coefficients (Figure 6). We then compare this to the same analyses performed on undifferenced data to show the effect of removing spurious correlations that are the result of common cycles and trends.

Figure 6

Standardized regression coefficients (AE, ULF, P, and E ) from multiple regression (not autoregressive moving average transfer function) models. 1. All data differenced by subtracting the previous hour's observation: during disturbed periods (a, d, and g), quiet periods (b, e, and h), and storm recovery periods (c, f, and i). 2. The same for undifferenced data. Note the difference in scale between 1. and 2. Significant effects (p‐value < 0.05) are shown in blue.

With differenced data (Figure 6.1), the AE effect is consistent over these three periods (strongest effect on the 40 keV flux, least effect on 150 keV flux). No matter the geomagnetic conditions, substorms (as measured by AE) show a statistically significant positive influence on flux, with the most effect at the lower electron energies. P does not contribute significantly at most periods or energy levels (the exception being at 150 keV during disturbed periods). E shows a negative effect in the quiet periods but a positive effect in recovery. ULF has little or a negative influence, even when periods are selected that would be expected to show a strong effect such as recovery following storms.

We present an analysis of undifferenced data in Figure 6.2 to show the danger of correlating variables with common diurnal cycles. In the undifferenced data, we do find the “expected” strong ULF effect (Figure 6.2; note the larger scale compared to the differenced data), but this is only a demonstration of the spurious nature of this high correlation. High correlations between ULF wave activity and electron flux in hourly data are likely only describing a common diurnal cycle and say little about physical driving mechanisms (Simms et al., 2022). Note that it is not so much that the correlation is “wrong” but that the differencing or ARMA modeling removes the portion of the correlation that is irrelevant to the questions we are interested in. ULF waves may be a more long term driver of flux, with positive influences only appearing after 24 hr (Simms et al., 2021). The other predictors also show stronger effects when not differenced (Figure 6.2), likely also due to common diurnal cycles in the data.

Discussion and Conclusions

A number of variables show high simple (single variable) correlations with keV electron flux, but by using an ARMAX analysis, which removes the confounding effect of diurnal cyclicity and allows assessment of each parameter independently, we show more definitively that substorms (measured by AE) are the most influential process at 40 and 75 keV. This accords with previous work that found substorms to be an important correlate with both keV (Ganushkina et al., 2021) and MeV electrons (Simms, Engebretson, Clilverd, Rodger, Lessard, et al., 2018).

There are major differences between the single variable correlations and the full (all variable) ARMAX results. Certain variables lose most of their apparent influence when all possible drivers are considered at once. Some parameters, for example, AE, Kp and SymH, may show a decrease in correlation because they all essentially describe the same geomagnetic perturbations. Although all three geomagnetic indices (Kp, SymH, and AE) show high simple (single variable) correlations with electron flux, the influences of Kp and SymH disappear in a full regression model where other variables are included. It is likely that these two indices mostly measure generalized disturbance in the magnetosphere which is better described using solar wind and IMF variables. However, the AE index, as it is better positioned to measure substorms and subsequent electron injections, is more representative of the physical processes that drive flux.

Other variables may appear to lose influence because most of their correlation was due to sharing a common diurnal cycle with the satellite‐observed electron flux data. If certain of these variables were measured by the same satellite, this would be an obvious source of error. However, even the ULF ground index that we use shows its own diurnal cycle and this will inflate the correlation with flux (Simms et al., 2022). Adding ARMA terms to partition out the variance due to common cycling leaves us with the portion of the variance that actually describes the relationship between predictor and flux.

Other variables are simply highly correlated with each other (e.g., N and V), while others are derived from these same variables (P and E ). This means that single variable analyses, with no correction for cycles and trends, are highly unreliable for pinpointing the actual correlates of electron flux, and also that variable sets should be chosen to minimize these intercorrelations and to maximize the variables of possible direct influence.

In the overall ARMAX models, there is a somewhat lesser effect of Ey (calculated as ‐V B ) in contrast to previous single‐variable studies (Denton et al., 2016) (but see below, as this may be due to varying influence of this factor over the 24 hr period: positive over 23‐8 MLT, negative effect over 9–22 MLT). P is more influential at 150 keV, acting to decrease electron flux. The contrast to previous findings, where pressure showed a positive association with flux (Shi et al., 2009), is due to our present study incorporating more predictors at one time. ULF shows little influence on keV electrons in these hourly, full variable models, despite its influence on MeV electrons (Simms, Engebretson, Clilverd, Rodger & Reeves et al., 2018; Simms, Engebretson, Clilverd, Rodger, Lessard, et al., 2018; Simms et al., 2021) and its strong positive correlation when it is the only predictor and the confounding effects of the diurnal cycle are not removed (Figure 1). This is not to say that ULF waves have no influence, but rather that the single variable correlation mischaracterizes these wave effects as much more dominant than they really are.

However, it should be noted that all correlational analysis of this system, whether single variable, multivariable, or correcting for cycles, is observational in nature, not experimental. As treatments cannot be randomly assigned (e.g., ULF waves cannot be increased or decreased to explore their effect), the correlations found can only be evidence of an association, not definitive proof of driving. However, despite this, we continue to do the best we can by correcting the issues that we are able to address.

As electron flux is log‐transformed in our analyses, all the relationships we find here are nonlinear even though they are tested with the linear model method of ARMAX regression. As B and E are not log‐transformed, they show an exponential relationship with electron flux. All other predictors, which are log‐transformed, are described by a power function relationship. The use of transformations such as the log‐transform or polynomial terms (Balikhin et al., 2011) are both able to create intrinsically linear regression models from many, but not all, data distributions. However, while both serve to linearize the relationship, the log‐transform is thought to also help the data to meet other conditions of the linear model framework such as equality and normality of the residual error (Neter et al., 1990).

While there are sizable simple correlations of some parameters with electron flux, single variable correlations can misrepresent the actual relationships. If neither common cycles and trends, nor confounding variables are accounted for, simple correlational analysis may show large associations between variables that have no physical relationship. This has been demonstrated before, where removal of common cycles results in either a complete elimination of a correlation between some space weather parameters (e.g., the commonly observed ULF wave correlation with solar wind velocity or with electron flux (Simms et al., 2022)) or a reduction in correlation (Simms et al., 2021)). An ARMAX model, used in this study, can account for common cycles in time series data (and trends, if necessary) by the use of AR and MA terms (and differencing, if needed). Entering several predictor variables into the same analysis then allows each variable's influence to be calculated while the others are held constant.

The Reduced Model: , , , and as Predictors

However, adding all possible explanatory variables to a model may not correctly identify the most important physical parameters but only those that correlate best, for whatever reason. While a reasonable predictive model may be achieved by using all available variables in a regression or neural network, leaving an algorithm to choose the model with the highest validation correlation, this is unlikely to identify actual drivers in the system. This approach, instead, can lead to several problems: (a) “opposing” variables may appear extremely influential as they compete to explain the same small bit of variation, (b) theoretical considerations of physical influence tend to be ignored in favor of factors that happen to correlate well, (c) coefficient estimates may be biased if extraneous variables are included or if important variables are excluded (Smith, 2018; Whittingham et al., 2006). In the worst case, a model may report that factors that cannot physically influence the dependent variable are the only factors that have any effect at all. For this reason, to determine whether a factor has an actual driving influence, care must be taken to choose only those for which a likely physical effect can be postulated and not just all that are available. This is why we have chosen to do further analyses on a set of presumed direct drivers (substorms, ULF waves, pressure, and electric field), as well as analyses that show the relative correlations of all possible variables.

Using the ARMAX method on such a reduced model, we find that the influence of substorms (AE) on hourly electron flux remains substantial over the 40–75 keV range at geostationary orbit (approximately L6) although of less importance at 150 keV. This influence is strongest after midnight into the mid‐morning hours MLT, coinciding with the post midnight to dawn injection of electrons from the magnetotail (Birn et al., 1997; Thomsen et al., 2001). The AE influence is slightly higher during storm recovery periods than during either disturbed or quiet periods. Substorms, therefore, are the dominant driver within our postulated “direct driver” set (substorms, ULF waves, solar wind pressure, and electric field) and presumably show the influx of electrons injected from the magnetotail.

The hourly E parameter (electric field of the solar wind) shows little influence when MLT is ignored. However, introducing MLT into the model results in a positive effect of E over 20‐8 MLT, corresponding with the observation that IMFB influence is strongest during this period as well (Dubyagin et al., 2016). There is, however, a mostly negative effect of E at other times of day. These opposing influences cancel each other out in a model that does not account for variations over MLT. The E influence also varies by geomagnetic conditions, with no influence during disturbed periods, a negative influence during quiet periods, and a positive influence during recovery after storms. This is likely a reflection of the magnetospheric electric convective field (which reflects the E ) breaking up during disturbed periods and becoming enhanced during storm recovery.

Overall, P shows a moderately negative direct effect on flux. When the analysis accounts for MLT, this negative influence is strongest over 20‐12 MLT.

ULF waves, thought to accelerate electrons to higher energies, show lower immediate (hourly) influence than the AE. A strong correlation of ULF waves with high energy electron flux (>1.5 MeV) found in previous studies may be a consequence of correlating two variables with a common diurnal cycle, or a reflection of only long term (at least day long) physical driving (with no short term influence), or both. We find here that the most significant short term driving of 40–75 keV electrons by ULF appears to be negative and only during quiet or recovery periods, while there is little short term effect at 150 keV. The negative short term effect in quiet or recovery periods accords with previous work on higher energy electrons that shows a negative influence of the ULF index on the first day, with more positive influences on the next day (Simms et al., 2021), but is not characteristic of the analyses using all hours. While it is possible that the ULF index does not effectively capture the immediate, localized, and frequency specific ULF wave power that would be most responsible for the resonant interaction driving electrons to these energies (e.g., Claudepierre et al. (2013)), leading to a lower apparent effect, we note that the ULF index influence is still statistically significant at all three energies in analyses using all hours (see Table 2). We have thus captured the importance of ULF waves as a driver, even if they appear to have less influence than AE once the diurnal cycle is removed.

While all four of these factors in the reduced model (AE, P, E and ULF) are statistically significant, the AE has up to 10 times greater influence on the 40–75 keV flux. This would seem to be the result of electron injection in the lower energy ranges being the major driving factor of increased electron flux. At 150 keV, P and ULF influences are somewhat stronger than at the lower energies, with their effects being similar to that of AE. Therefore, electron injection is no longer the most important factor at 150 keV. Flux, at the higher energy, is more dependent on acceleration due to ULF waves and losses due to magnetopause shadowing induced by pressure. As noted above, the E influence appears low because it switches in sign from a positive influence post‐midnight to a negative influence post‐noon.

The amount of variation explained by these reduced (4 factor) models (R 2 = 67%–76%) is high, but much of this is due to the AR and MA terms that describe the cycling behavior of flux. In other words, most of what is being explained is the background rise and fall of flux at the satellite as it orbits the Earth. The explanatory value due to AE and the other possible drivers pales in comparison, but this is not as discouraging as it may at first seem. Most of the hourly fluctuation in flux is the rather uninteresting observation that flux varies widely over MLT as the satellite passes through various field lines. We are more interested in the response of electron flux when this behavior is eliminated. Once this is removed, we have a clearer picture of which variables are most influential. The relevant comparison between correlations, then, is not to the single variable, uncorrected correlations of Figure 1 but between the coefficients of the tested possible drivers of Table 2, and as previously stated, all of these are statistically significant, although the AE is the strongest influence at 40–75 keV. (Note that the coefficients of the ARMAX models are not correlations. They cannot be compared directly to the single variable correlations of Figure 1).

The response of electron flux to our identified possible direct drivers (AE, ULF, P, and E ) varies only somewhat between disturbed, quiet, and storm recovery periods. AE is a stronger influence during recovery, for example, than during quiet or disturbed periods.

At 150 keV, there is the least response of hourly averaged flux to the presumed physical drivers. This may represent the longer time frame of action required from these processes to bring electrons to higher energies. Even the cross correlations (Figure 1) show higher effects from 24 to 48 hr previous, with ULF and AE showing their least influence in the 12 hr preceding a flux measurement and the E influence peaking at 12 hr.

Indirect Drivers ( , , and )

In addition to these variables that we label direct, physical drivers of flux, we consider several other parameters as possible indirect drivers (solar wind N and V and IMF Bz) which show fairly equivalent influences on flux via their effects on the direct drivers. This supports previous findings concerning these three solar wind and IMF influences (Ganushkina et al., 2019; Hartley et al., 2014; Kellerman & Shprits, 2012; Li et al., 2005; Sillanpää et al., 2017). Stepanov et al. (2021) when controlling for other variables, also found solar wind velocity and a magnetospheric convection variable (the dayside merging electric field, similar to the E we use) to be the strongest correlates of keV flux near the plasmasheet midplane (A similar multiplicative variable, the IMF factor (Balikhin et al., 2010; Boynton et al., 2011), and solar wind velocity appear to control hourly averaged 40 keV electrons. However, these last studies did not include a test of AE influence for comparison).

We note again that we are not trying to predict flux. We are trying to understand what drives it. We do not test N, V, and B for their direct effect on flux in the reduced model for the simple reason that we have not hypothesized a direct physical connection between these variables and flux. There is little theory that would suggest this. There is theoretical work suggesting direct influences of pressure (the combination of N and V 2), AE (through substorm injection of electrons), and ULF waves (either accelerating electrons to these energies or causing electron loss). Therefore, we test those particular hypotheses with an ARMAX regression model (the upper part of Figure 4 and the description in Section 4.1). Although V or N may show a strong single variable correlation with electron flux, this is either via their role in pressure, or via their influence on ULF waves. These second hypotheses about the indirect V or N influence on flux are tested in the lower part of Figure 4. By only including hypothesized direct drivers in the top part of Figure 4, we are better able to test these particular hypotheses. If we include N and V, as in the full regression model, then we are no longer testing hypotheses about direct drivers and, in fact, are making those particular hypothesis tests impossible because of the interference between these variables.

We are able to compare effects of the other correlates by summing their indirect influence through the presumed physical drivers. We are able to calculate that at 40 keV, P shows a summed influence (both direct and indirect) nearly half that of the most influential parameter, AE, with E having about a fourth the influence of AE. Of the postulated indirect drivers, N, B , and V show nearly equal effects. The N and V influences are negative, while the B influence switches sign above 75 keV.

The Problem of Prediction Versus Driver Identification

If the purpose of a model is accurate prediction, then a simple validation correlation of observation with prediction on a withheld test set is the statistic of interest. In that case, predictor variables can be chosen simply on the basis of availability and ability to correlate well with the response. Alternatively, the ARMAX‐regression models we present here address the question of what parameters drive flux changes. We use hypothesis testing within the ARMAX‐regression framework to determine whether certain parameters show an association with electron flux. As our questions concern the science of the system (i.e., which variables are drivers), we consider, first, which variables most justifiably have a physical association with flux and which are only highly correlated because they are proxies. Of particular concern is the removal of the diurnal cycle due to satellite orbit. This factor alone may increase apparent correlations tremendously without having any bearing on the physical driving relationships.

A model such as this, developed for determining the actual relationships, statistically tests the slope of association with flux for each identified variable, with less concern attached to either a validation correlation or the overall R 2 (proportion of variation explained). That these more focused hypothesis tests appear to explain less than the unspecific correlations only reflects the removal of spurious correlations that lead to an incorrect understanding of the system.