Compacting the description of a time-dependent multivariable system and its multivariable driver by reducing the state vectors to aggregate scalars: the Earth's solar-wind-driven magnetosphere

Using the solar-wind-driven magnetosphere–ionosphere–thermosphere system, a methodology is developed to reduce a state-vector description of a time-dependent driven system to a composite scalar picture of the activity in the system. The technique uses canonical correlation analysis to reduce the time-dependent system and driver state vectors to time-dependent system and driver scalars, with the scalars describing the response in the system that is most-closely related to the driver. This reduced description has advantages: low noise, high prediction efficiency, linearity in the described system response to the driver, and compactness. The methodology identifies independent modes of reaction of a system to its driver. The analysis of the magnetospheric system is demonstrated. Using autocorrelation analysis, Jensen–Shannon complexity analysis, and permutation-entropy analysis the properties of the derived aggregate scalars are assessed and a new mode of reaction of the magnetosphere to the solar wind is found. This state-vector-reduction technique may be useful for other multivariable systems driven by multiple inputs.

1 Introduction

In this report a methodology is described that can produce a compact description of the behavior of a time-dependent, multivariable system driven by a time-dependent, multivariable driver or by multiple drivers. The system used to develop this methodology is the Earth's magnetosphere–ionosphere–thermosphere system driven by the time-dependent solar wind. The spatial domain wherein the Earth's magnetic field dominates over the solar wind is known as the magnetosphere. The interaction between the solar wind, and the magnetosphere is surprisingly complex and the magnetosphere's evolution in response to the time-varying solar wind is rich and diverse. The magnetospheric system is characterized by multiple subsystems that interact with each other (cf. Lyon, 2000; Otto, 2005; Siscoe, 2011; Eastwood et al., 2015; Borovsky and Valdivia, 2018): almost 6 orders of magnitude of spatial scales are involved in the global behavior of the magnetosphere, from ∼1 to $\sim \mathrm{6}\times {\mathrm{10}}^{\mathrm{5}}$ km. This system is highly coupled, dynamic, with memory and with feedback loops. Multiple physical processes act to couple the various subsystems, with the strength of the couplings evolving with time as the subsystems evolve owing to the couplings. Even after a half of a century of measurements and analysis, its subsystems and the couplings between its subsystems are not fully understood (Stern, 1989, 1996; Denton et al., 2016). It has been argued that the system adjectives “adaptive”, “nonlinear”, “dissipative”, and “complex” apply to the magnetospheric system (Borovsky and Valdivia, 2018). (See also the earlier systems analyses by Horton et al., 1999, Chapman et al., 2004, Valdivia et al., 2005, 2013, and Sharma, 2010.) The magnetospheric system is well measured: there are hundreds of thousands of hours of simultaneous measurements of various aspects of the magnetospheric system and its solar-wind driver over the five decades of the “space age” (cf. Stern, 1989, 1996; King and Papitashvili, 2005).

The methodology developed here creates a compact description of a time-dependent, multivariable system driven by a time-dependent, multivariable driver. The diverse variables describing the system may be intercorrelated, and the variables describing the driving may be intercorrelated. The methodology was developed to gain an understanding of the Earth's magnetosphere–ionosphere–thermosphere system as driven by the solar wind. To utilize the methodology the system and its driver are conceptualized by a time-dependent, multidimensional system state vector S(t) and a time-dependent, multidimensional driver state vector D(t), with the assumption that the driver vector D affects the system vector S, but not vice versa, written D→S. The individual time-dependent scalar variables making up the state vector S(t) are time-dependent measures of various forms of activity in the system and various properties of the system, and the individual time-dependent scalar variables making up the driver state vector D(t) are various time-dependent measures of the properties of the drivers of the system. We will utilize the correlation properties between the components (individual time-dependent variables) of S and the components of D. Canonical correlation analysis (CCA) will be used to derive scalar projections (dot products) of S(t) and scalar projections of D(t) that have the highest Pearson linear correlation coefficient between them. The derived scalar projections S₍₁₎(t), S₍₂₎(t), S₍₃₎(t), … of the vector S(t) will be composite (aggregate) measures of activity in the system, and the derived scalar projections D₍₁₎(t), D₍₂₎(t), D₍₃₎(t), … of the vector D(t) will be the composite drivers of S₍₁₎(t), S₍₂₎(t), S₍₃₎(t), …, respectively. In essence, the aggregate variables S₍₁₎(t), S₍₂₎(t), S₍₃₎(t), … are “latent variables” of the system constructed from the “manifest variables” in the system state vector. This reduced scalar picture D_(i)→S_(i) of the system driven by the driver focuses on the time-dependent properties of the system that react to the driver. By maximizing the correlations, the predictability of the system from knowledge of the state of the driver is also maximized.

The solar-wind-driven magnetospheric system very cleanly follows the D→S picture where the driver affects the system, but the system does not affect the driver. The Earth's magnetosphere has no influence whatsoever on the properties of the solar wind that passes the Earth. Measurements of this magnetospheric system will be used in Sects. 2 and 3 to explore the mathematical reduction of the state-vector D(t)→S(t) picture to the composite-scalar D_(i)(t)→S_(i)(t) picture. Table 1 lists the nine time-dependent measurements of the magnetosphere in the system state vector S and the eight time-dependent measurements of the solar wind in the driver state vector D. The individual variables in the system state vector and in the driver state vector are described in the Appendix.

Table 1 The nine time-dependent variables going into the system state vector S(t) of the magnetosphere and the eight time-dependent variables going into the driver state vector D(t) of the solar wind.

This report is organized as follows. In Sect. 2 the CCA approach is applied to the magnetospheric system driven by the solar wind to derive the first three time-dependent sets of composite variables S₍₁₎(t) and D₍₁₎(t), S₍₂₎(t) and D₍₂₎(t), and S₍₃₎(t) and D₍₃₎(t) from the state vectors S(t) and D(t). In Sect. 3 the three sets of composite variables S_(i) and D_(i) for the magnetospheric system are explored, and the complexity–entropy properties of the aggregate variable S₍₁₎(t) are analyzed. In Sect. 4 the advantages of the reduced D_(i)→S_(i) scalar description are examined: these advantages include (a) a compact description of global system-wide reactions to variations in the driver, (b) increased predictability of the system from knowledge of the driver, (c) linearity in the description of the system's response to the driver, and (d) lower noise in correlations between the system variables and the driver variables. The reduced scalar picture can also reveal independent modes of reaction of the system to the driver, providing insight into the behavior of the system in reaction to complexities in the driver. The variables of the magnetospheric and solar-wind state vectors are described in the Appendix.

2 Creation of composite (aggregate) variables from the state vectors

Using the Earth's magnetosphere–ionosphere–thermosphere system as driven by the solar wind, the reduction of a time-dependent state-vector picture D(t)→S(t) to the time-dependent composite-variable-pair picture D_(i)(t)→S_(i)(t) will be performed. The nine measured variables chosen for the nine-dimensional magnetospheric system state vector S appear in the first column of Table 1, and the eight measured variables chosen for the eight-dimensional solar-wind driver state vector D appear in the second column of Table 1, with explanations of those measures deferred to the Appendix.

One-hour averages of all magnetospheric and solar-wind variables are used in the years 1991–2007. No time lags are used between the solar-wind measurements and the magnetospheric measurements: most expected time lags will be about 1 h (e.g., Clauer et al., 1981; Smith et al., 1999), which is the time resolution of the data set.

Canonical correlation analysis (CCA) is applied to the time-dependent state vectors S(t) and D(t). CCA finds correlation patterns between two multivariable data sets (Nimon et al., 2010; Hair et al., 2010). It yields pairs of composite (aggregate) variables (a) that are linear combinations of the variables of the two data sets and (b) that have maximal correlations with each other. Each pair of composite variables is called the “Nth canonical correlation”. From the data sets of S(t) and D(t) the first pair of composite variables yielded (the first canonical variates) is S₍₁₎(t) and D₍₁₎(t): these two variables are projections of S and D given by $S_{\left(\mathrm{1}\right)}\left(t\right)={\mathit{C}}_{S\mathrm{1}}\cdot \mathit{S}\left(t\right)$ and $D_{\left(\mathrm{1}\right)}\left(t\right)={\mathit{C}}_{D\mathrm{1}}\cdot \mathit{D}\left(t\right)$, where C_S1 and C_D1 are time-independent coefficient (weight) vectors. S₍₁₎ and D₍₁₎ are the composite variables from S and D that have the highest Pearson linear correlation coefficient with each other. Here, CCA is in a sense creating the system function S₍₁₎(t) that is most reactive to the driver vector D(t) and creating the driver scalar function D₍₁₎(t) that describes that driving. CCA then yields other pairs of composite variables S₍₂₎ and D₍₂₎ (the second canonical correlation), S₍₃₎ and D₍₃₎ (the third canonical correlation), etc. S₍₂₎ and D₍₂₎ are the projections of S and D that have the highest correlation with each other, provided that S₍₂₎ and D₍₂₎ are uncorrelated with S₍₁₎ and D₍₁₎. S₍₃₎ and D₍₃₎ are the projections of S and D that have the highest correlations with each other, provided they are uncorrelated with S₍₁₎, S₍₂₎, D₍₁₎, and D₍₂₎. S₍₁₎(t), S₍₂₎(t), and S₍₃₎(t) represent three independent modes of reaction of the global system to the driver D(t). The CCA process will identify these modes (and their respective drivers).

CCA is a matrix equation solution, non-iterative, that yields a single unique solution (Johnson and Wichern, 2007). CCA operates on standardized variables (with the mean value subtracted and the values then divided by the standard deviation), denoted with an asterisk. (For each variable the mean value and standard deviation are calculated for the entire data set.) CCA operates most efficiently on variables that are Gaussian distributed: hence the logarithms of some variables are used to yield more-Gaussian-like distributions. All standardized variables v^* have a mean value of zero, a standard deviation unity, and no units.

When CCA is applied to the 1991–2007 S(t) and D(t) data sets (see Table 1), the first canonical pair of time-dependent variables is

$$\begin{array}{}\text{(1a)}& {\displaystyle }\begin{array}{rl}{S}_{\left(\mathrm{1}\right)}& =\mathrm{0.0260}{\log }_{\mathrm{10}}{\left(\mathrm{1}+\left|\mathrm{AL}\right|\right)}^{*}\\ & +\mathrm{0.1151}{\log }_{\mathrm{10}}{\left(\mathrm{1}+\left|\mathrm{AU}\right|\right)}^{*}+\mathrm{0.2160}|\mathrm{PCI}{|}^{*}\\ & +\mathrm{0.1451}{K_{\mathrm{p}}}^{*}+\mathrm{0.2881}{\log }_{\mathrm{10}}{\left(\mathrm{1}+\mathrm{am}\right)}^{*}\\ & +\mathrm{0.0201}d\left|D_{\mathrm{st}}\right|/\mathrm{d}{t}^{*}+\mathrm{0.0492}{\log }_{\mathrm{10}}{\left(\mathrm{0.01}+{\mathrm{mP}}_{\mathrm{e}}\right)}^{*}\\ & +\mathrm{0.2531}{\log }_{\mathrm{10}}{\left(\mathrm{0.01}+{\mathrm{mP}}_{\mathrm{i}}\right)}^{*}\\ & +\mathrm{0.0854}{\log }_{\mathrm{10}}{\left(\mathrm{0.01}+P_{\mathrm{ips}}\right)}^{*}\end{array}\text{(1b)}& {\displaystyle }\begin{array}{rl}{D}_{\left(\mathrm{1}\right)}& =\mathrm{0.8378}{\log }_{\mathrm{10}}{\left(v_{\mathrm{sw}}\right)}^{*}+\mathrm{0.6876}{\log }_{\mathrm{10}}{\left(n_{\mathrm{sw}}\right)}^{*}\\ & +\mathrm{0.1018}{\log }_{\mathrm{10}}{\left(F_{\mathrm{10.7}}\right)}^{*}-\mathrm{0.1676}{\left(-B_z\right)}^{*}\\ & +\mathrm{0.3547}f{\left(M\right)}^{*}+\mathrm{0.3844}{〈{\sin }^{\mathrm{2}}\left({\mathit{\theta }}_{\mathrm{clock}}/\mathrm{2}\right)〉}_{\mathrm{3}}^{*}\\ & +\mathrm{0.0960}{〈{\mathit{\theta }}_{\mathrm{Bn}}〉}_{\mathrm{3}}^{*}+\mathrm{0.0943}{\log }_{\mathrm{10}}{\left(\mathrm{0.1}+\left|\mathrm{\Delta }B\right|\right)}^{*}.\end{array}\end{array}$$

S₍₁₎ and D₍₁₎ have mean values of zero and standard deviations of unity. The derived composite variables given by expressions (1) are robust and reproducible: applying the CCA process to various subsets of the full 1991–2007 data set, the CCA process repeatedly yields essentially the same coefficients that are in Eqs. (1a) and (1b) (cf. Borovsky and Denton, 2018).

CCA applied to the time-dependent state vectors S(t) and D(t) for the 1991–2007 data set yields the second canonical pair of time-dependent scalar variables as

$$\begin{array}{}\text{(2a)}& {\displaystyle }\begin{array}{rl}{S}_{\left(\mathrm{2}\right)}& =-\mathrm{0.2628}{\log }_{\mathrm{10}}{\left(\mathrm{1}+\left|\mathrm{AL}\right|\right)}^{*}\\ & -\mathrm{0.0874}{\log }_{\mathrm{10}}{\left(\mathrm{1}+\left|\mathrm{AU}\right|\right)}^{*}-\mathrm{0.1302}|\mathrm{PCI}{|}^{*}\\ & -\mathrm{0.0556}{K_{\mathrm{p}}}^{*}+\mathrm{0.1928}{\log }_{\mathrm{10}}{\left(\mathrm{1}+\mathrm{am}\right)}^{*}\\ & +\mathrm{0.0028}d\left|D_{\mathrm{st}}\right|/\mathrm{d}{t}^{*}-\mathrm{0.8506}{\log }_{\mathrm{10}}{\left(\mathrm{0.01}+{\mathrm{mP}}_{\mathrm{e}}\right)}^{*}\\ & +\mathrm{0.9218}{\log }_{\mathrm{10}}{\left(\mathrm{0.01}+{\mathrm{mP}}_{\mathrm{i}}\right)}^{*}\\ & +\mathrm{0.3493}{\log }_{\mathrm{10}}{\left(\mathrm{0.01}+P_{\mathrm{ips}}\right)}^{*}\end{array}\text{(2b)}& {\displaystyle }\begin{array}{rl}{D}_{\left(\mathrm{2}\right)}& =\mathrm{0.1195}{\log }_{\mathrm{10}}{\left(v_{\mathrm{sw}}\right)}^{*}+\mathrm{0.8874}{\log }_{\mathrm{10}}{\left(n_{\mathrm{sw}}\right)}^{*}\\ & +\mathrm{0.1202}{\log }_{\mathrm{10}}{\left(F_{\mathrm{10.7}}\right)}^{*}-\mathrm{0.1138}{\left(-B_z\right)}^{*}\\ & +\mathrm{0.2669}f{\left(M\right)}^{*}-\mathrm{0.5079}{{〈{\sin }^{\mathrm{2}}\left({\mathit{\theta }}_{\mathrm{clock}}/\mathrm{2}\right)〉}_{\mathrm{3}}}^{*}\\ & -\mathrm{0.0186}{〈{\mathit{\theta }}_{\mathrm{Bn}}〉}_{\mathrm{3}}^{*}+\mathrm{0.0260}{\log }_{\mathrm{10}}{\left(\mathrm{0.1}+\left|\mathrm{\Delta }B\right|\right)}^{*}.\end{array}\end{array}$$

For the 1991–2007 data set CCA yields the third canonical pair of time-dependent scalar variables as

$$\begin{array}{}\text{(3a)}& {\displaystyle }\begin{array}{rl}{S}_{\left(\mathrm{3}\right)}& =-\mathrm{0.1796}{\log }_{\mathrm{10}}{\left(\mathrm{1}+\left|\mathrm{AL}\right|\right)}^{*}\\ & -\mathrm{0.2220}{\log }_{\mathrm{10}}{\left(\mathrm{1}+\left|\mathrm{AU}\right|\right)}^{*}-\mathrm{1.0351}|\mathrm{PCI}{|}^{*}\\ & +\mathrm{0.8265}{K_{\mathrm{p}}}^{*}+\mathrm{0.5809}{\log }_{\mathrm{10}}{\left(\mathrm{1}+\mathrm{am}\right)}^{*}\\ & -\mathrm{0.2169}d\left|D_{\mathrm{st}}\right|/\mathrm{d}{t}^{*}+\mathrm{0.3856}{\log }_{\mathrm{10}}{\left(\mathrm{0.01}+{\mathrm{mP}}_{\mathrm{e}}\right)}^{*}\\ & -\mathrm{0.6100}{\log }_{\mathrm{10}}{\left(\mathrm{0.01}+{\mathrm{mP}}_{\mathrm{i}}\right)}^{*}\\ & +\mathrm{0.1064}{\log }_{\mathrm{10}}{\left(\mathrm{0.01}+P_{\mathrm{ips}}\right)}^{*}\end{array}\text{(3b)}& {\displaystyle }\begin{array}{rl}{D}_{\left(\mathrm{3}\right)}& =\mathrm{0.4241}{\log }_{\mathrm{10}}{\left(v_{\mathrm{sw}}\right)}^{*}-\mathrm{0.1985}{\log }_{\mathrm{10}}{\left(n_{\mathrm{sw}}\right)}^{*}\\ & -\mathrm{0.1404}{\log }_{\mathrm{10}}{\left(F_{\mathrm{10.7}}\right)}^{*}-\mathrm{0.6704}{\left(-B_z\right)}^{*}\\ & -\mathrm{0.1008}f{\left(M\right)}^{*}+\mathrm{0.0572}{〈{\sin }^{\mathrm{2}}\left({\mathit{\theta }}_{\mathrm{clock}}/\mathrm{2}\right)〉}_{\mathrm{3}}^{*}\\ & -\mathrm{0.3134}{〈{\mathit{\theta }}_{\mathrm{Bn}}〉}_{\mathrm{3}}^{*}+\mathrm{0.3055}{\log }_{\mathrm{10}}{\left(\mathrm{0.1}+\left|\mathrm{\Delta }B\right|\right)}^{*}.\end{array}\end{array}$$

The properties of S₍₁₎(t) and D₍₁₎(t), S₍₂₎(t) and D₍₂₎(t), and S₍₃₎(t) and D₍₃₎(t) as given by Eqs. (1)–(3) are explored in Sect. 3.

3 Properties of the scalar reduced picture for the magnetospheric system

The three sets of composite variables S₍₁₎ and D₍₁₎, S₍₂₎ and D₍₂₎, and S₍₃₎ and D₍₃₎ for the magnetospheric system are explored and the advantages of the reduced D_(i)→S_(i) scalar description are investigated.

3.1 The primary mode of system response as represented by D₍₁₎→S₍₁₎

In Fig. 1 the composite system variable S₍₁₎ (as given by Eq. 1a) is plotted for the years 1991–2007 as a function of the composite driver variable D₍₁₎ (as given by Eq. 1b). Each black point in Fig. 1 represents 1 h of data. The Pearson linear correlation coefficient between S₍₁₎ and D₍₁₎ for the 1991–2007 data set is r_corr=0.921. Accordingly, $r_{\mathrm{corr}}^{\mathrm{2}}=\mathrm{84.8}$ % of the variance of the system function S₍₁₎(t) is described by the driver function D₍₁₎(t), and so 15.2 % of the variance of S₍₁₎ is unaccounted for by D₍₁₎. The blue line in Fig. 1 is a linear-regression fit to S₍₁₎, and the red curve is a 50-point vertical running average of the black points. Note the approximate linearity of system variable S₍₁₎ as a function of driver variable D₍₁₎, indicated by the manner in which the running average tracks the linear-regression line.

Figure 1 The aggregate system scalar S₍₁₎ is plotted as a function of the driver scalar D₍₁₎ for the 1 h resolution 1991–2007 data set. Each black point is 1 h of data.

Note that whereas the correlation coefficient between S₍₁₎(t) and D₍₁₎(t) is r_corr=0.921, the maximum correlation coefficient between any single variable in the system state vector S(t) and any single variable in the driver state vector D(t) is only r_corr=0.586 (between ${〈{\sin }^{\mathrm{2}}({\mathit{\theta }}_{\mathrm{clock}}/\mathrm{2})〉}_{\mathrm{3}}$ and ${\log }_{\mathrm{10}}(\mathrm{1}+|\mathrm{AL}\left|\right))$.

As a further note, the Pearson linear correlation coefficients between S₍₁₎ and various “physics-based” solar-wind driver functions from the literature are the following: +0.378 for −v_swB_z, +0.557 for v_swB_south (Eq. 2 of Holzer and Slavin, 1979), +0.679 for the Newell function dΦ∕dt (Eq. 1 of Newell et al., 2007), +0.723 for the quick reconnection function R_q (Eq. 8 of Borovsky and Birn, 2014), and +0.761 for the nonlinear reconnection-coupled MHD generator with Bohm viscosity (Eq. 65 of Borovsky, 2013). All of these driver functions have poor correlations with S₍₁₎ in comparison with the +0.921 correlation of D₍₁₎ with S₍₁₎.

Figure 2 Plots of the nine components of the coefficient vectors used to project the system state vector S into the aggregate variables S₍₁₎ (a), S₍₂₎ (c), and S₍₃₎ (e) and plots of the eight components of the coefficient vectors used to project the driver state vector D into the driver scalar variables D₍₁₎ (b), D₍₂₎ (d), and D₍₃₎ (f).

In the six panels of Fig. 2 the coefficients of the six vectors C_S1, C_D1, C_S2, C_D2, C_S3, and C_D3 are plotted. (These are the coefficients in Eqs. 1–3.) Examining these six panels enables the reaction modes represented by S₍₁₎, S₍₂₎, and S₍₃₎ to be interpreted as well as their drivers D₍₁₎, D₍₂₎, and D₍₃₎. Figure 2a indicates that all coefficients of S₍₁₎ are positive: this indicates a mode of the magnetospheric system in which all measures of activity in the system vector S increase or decrease in unison, with S₍₁₎ representing a “global activity index”. Figure 2b indicates that all of the coefficients of D₍₁₎ are positive. The variables in the driver state vector D (Table 1) and their signs were all chosen so that a positive increase in each variable would result in a generally accepted increase in magnetospheric activity. The individual variables on the right-hand side of Eq. (1b) have all been correlatively associated with the driving of magnetospheric activity (Berthelier, 1976; Borovsky and Funsten, 2003; Newell et al., 2007; Borovsky and Denton, 2014; Borovsky and Birn, 2014; Osmane et al., 2015). S₍₁₎ is selected by the CCA process to have highest correlation with solar-wind variability: S₍₁₎ is focused on activity that reacts to the solar-wind driver.

Using the linear-regression curve in Fig. 1 as a “prediction” of the value of S₍₁₎ from knowledge of the value of D₍₁₎ yields

$$\begin{array}{}\text{(4)}& S_{\left(\mathrm{1}\right)\mathrm{pred}}=\mathrm{0.9209}{D}_{\left(\mathrm{1}\right)}-\mathrm{4.4}\times {\mathrm{10}}^{-\mathrm{5}}.\end{array}$$

In Fig. 3 the autocorrelation functions of S₍₁₎(t) (red curve), D₍₁₎(t) (blue curve), and S₍₁₎(t)−S_(1)pred(t) (green curve) are plotted. In Fig. 3a it is seen that the autocorrelation functions of S₍₁₎ and D₍₁₎ are very similar, with 1∕e autocorrelation times of 23.3 h for S₍₁₎ and 22.7 h for D₍₁₎. In Fig. 3b the three autocorrelation functions are plotted for time shifts up to 40 d. Note the 27 d peak in the autocorrelation functions of D₍₁₎(t) and S₍₁₎(t): this is associated with the 27 d rotation period of the Sun as viewed from the Earth and the persistence of features on the solar surface that give rise to solar wind with characteristic properties. This causes the driver D(t) properties to have a 27 d periodicity, which drives the system S(t) with a 27 d periodicity.

Figure 3 The autocorrelation functions for the system scalar S₍₁₎ (red), the driver scalar D₍₁₎ (blue), and the unaccounted-for variance E₍₁₎E_(1)pred (green) are plotted. In panel (a) the plot extends to 50 h and in panel (b) the plot extends to 40 d.

The quantity S₍₁₎−S_(1)pred is the portion of S₍₁₎(t) that is not accounted for by D₍₁₎(t), i.e., the unaccounted-for variance of S₍₁₎(t). S₍₁₎(t)−S_(1)pred(t) is completely uncorrelated with D₍₁₎(t). Further, S₍₁₎(t)−S_(1)pred(t) is completely uncorrelated with each of the eight individual solar-wind variables on the right-hand side of Eq. (1b). Since S₍₁₎ is so similar to D₍₁₎, the standard analyses of the S₍₁₎(t) time series (e.g., determining the correlation dimension, examining the state space, or Fourier analyzing; Sharma et al., 2005a; Vassiliadis, 2006) would largely be an analysis of the properties of the solar-wind time series D₍₁₎(t) – not so for S₍₁₎(t)−S_(1)pred(t), which is uncorrelated with D₍₁₎. The autocorrelation function of S₍₁₎(t)−S_(1)pred(t) in Fig. 3a is very different from the autocorrelation function of D₍₁₎: the 1∕e autocorrelation time of S₍₁₎(t)−S_(1)pred(t) is 2.4 h. Determining what the unaccounted-for variance S₍₁₎−S_(1)pred originates from is of great interest. Four suggestions of what contributes to S₍₁₎−S_(1)pred are made here. First, some fraction of S₍₁₎−S_(1)pred may be associated with noise in the various measurements of the magnetospheric system and of the solar wind. Shot noise (random noise in the values of the variables) would have an autocorrelation time of less than 1 h, the autocorrelation function of the shot-noise going from 1 to 0 in one data-resolution time shift (cf. Sect. 2.4 of Borovsky et al., 1997). Second, some fraction of S₍₁₎(t)−S_(1)pred(t) may be owed to errors in the measurement values in the state vectors S(t) and D(t). Errors in the values of the variables of D could be caused by the spatial structure of the solar wind and the measuring spacecraft upstream of the Earth not intercepting the exact solar-wind structures that hit and drive the Earth (cf. Weimer et al., 2003; Borovsky, 2018a): this could affect all of the variables of D. Extrapolating local measures to estimate global properties can also lead to errors: this might affect the hemispheric particle-precipitation variables mP_e and mP_i (Emery et al., 2008) in S and also the magnetospheric pressure values P_ips (Borovsky, 2017) in S. Variables reacting to more than one physical process (such as $\mathrm{d}\left|D_{\mathrm{st}}\right|/\mathrm{d}t$ and P_ips) could also appear to have error in the values when relating the values to D₍₁₎. Third, unaccounted-for time lags between solar-wind variables and magnetospheric variables may be resulting in weakened correlations: most time lags are 1 h or less, but measurements of magnetospheric particle populations can have lags of several hours (e.g., Denton and Borovsky, 2009; Borovsky, 2017). The fourth suggestion is that some fraction of S₍₁₎(t)−S_(1)pred(t) might be associated with system variations that are not directly associated with the solar-wind driver as measured by D. The autocorrelation time of S₍₁₎(t)−S_(1)pred(t) is approximately the 2–3 h time duration of a magnetospheric substorm (Borovsky et al., 1993; Weimer, 1994; Chu et al., 2015). Substorms are large transients in the reaction of the magnetospheric system to solar-wind driving. (Substorms have been described as self-organized criticality events in the driven magnetospheric system; Klimas et al., 2000.) The occurrence of a substorm is notoriously difficult to predict from solar-wind data (Freeman and Morley, 2004; Hsu and McPherron, 2009; Newell and Liou, 2011). The timing of substorm occurrence would be particularly difficult to infer from the 1 h resolution variables going into D because of the 3 h smoothing used on the clock-angle term ${〈{\sin }^{\mathrm{2}}({\mathit{\theta }}_{\mathrm{clock}}/\mathrm{2})〉}_{\mathrm{3}}$ in Eq. (1b) for D₍₁₎, with the clock angle being critical for substorm occurrence (Newell and Liou, 2011). The occurrence of a substorm would produce signatures in many of the variables used in S₍₁₎, typically an enhancement in the variable's amplitude lasting 2–3 h (Weimer, 1994).

To investigate this substorm hypothesis for S₍₁₎(t)−S_(1)pred(t), the variables D₍₁₎(t), S₍₁₎(t), and S₍₁₎(t)−S_(1)pred(t) are superposed-epoch averaged in Fig. 4 for a collection of 2155 substorm events; the collection is from Borovsky and Yakymenko (2017). The zero epoch in Fig. 4 is the onset time of each of the 2155 substorms. Substorms are associated with intervals of driving of the magnetosphere (e.g., Caan et al., 1977; Morley and Freeman, 2007); this is indicated by the increase in the superposed average of D₍₁₎ beginning prior to the onset time in Fig. 4. However, substorms also represent a transient release of stored energy in the magnetosphere (Birn et al., 2006); this is indicated in Fig. 4 by the superposed average of S₍₁₎ exceeding the superposed average of D₍₁₎ after the substorm onset and by the positive perturbation of the S₍₁₎−S_(1)pred curve after onset. The S₍₁₎−S_(1)pred curve indicates a transient in S₍₁₎ that is unaccounted for by D₍₁₎ associated with the occurrence of substorms. The autocorrelation time of the green superposed-average S₍₁₎−S_(1)pred time series plotted in Fig. 4 is 2.6 h, similar to the Fig. 1 autocorrelation time of the full 1991–2007 time series of S₍₁₎−S_(1)pred.

Figure 4 Superposed epoch averages of S₍₁₎ (red), D₍₁₎ (blue), and S₍₁₎−S_(1)pred (green) for 2155 substorms. The epoch time (t=0) is the time of onset of each substorm.

Compacting the description of the system from a multidimensional state vector to a few variables S₍₁₎, S₍₂₎, S₍₃₎, … is a form of dimensional reduction to a small set of fundamental latent variables: in that dimensional reduction a potential question is whether the reduced (more-fundamental) variables themselves exhibit a reduction of their embedding dimension from the embedding dimensions of the manifest variables in the state vector. Additionally, it would be valuable to differentiate S₍₁₎ from other indices commonly used to characterize magnetospheric activity. In order to achieve this task we use the methodology of Rosso et al. (2007) based on the combined use of permutation entropy (Bandt and Pompe, 2002) and Jensen–Shannon complexity mapping. This mapping developed by Rosso et al. (2007) is particularly useful to disentangle deterministic and stochastic time series. Additionally, one can use the Jensen–Shannon complexity alone to extract correlational structures in time series. In other words, one can seek timescales upon which coherent structures/fluctuations/modes arise. In some instances, coherent modes are signatures of deterministic chaos (Maggs and Morales, 2013) or a reduction to the number of degrees of freedom (Osmane et al., 2019). Readers with little familiarity with these two information theoretic measures can consult the reviews of Riedl et al. (2013) and Zanin et al. (2012) or the pedestrian methodology section found in Osmane et al. (2019). In Figs. 5 and 6, we map the values of AL (red), am (blue), S₍₁₎ (black) and D₍₁₎ (pink) on the complexity–entropy plane for an interval of 500 h duration with no data gaps, embedding dimensions of d=4, and embedding delay T ranging between 2 and 40 h. Because there are a total of d!=24 ordinal patterns we are limited to embedding delays of the order of 2 d. For embedding delays greater than T=48 h the number of segments $N-(d-\mathrm{1}{)}^{*}T$ becomes too small to ascertain the likelihood of forbidden ordinal patterns. Error bars for the Jensen–Shannon complexity, shown for the zoomed panels of the complexity–entropy planes, are computed as the square root of the number of ordinal patterns divided by the number of segments available. Hence, larger embedded dimensions d, require a larger number of segments $N-(d-\mathrm{1}{)}^{*}T$ to determine whether the Jensen–Shannon complexity lies significantly above the stochastic boundary (see below).

Figure 5 Jensen–Shannon complexity and permutation entropy, AL (red) and am (blue) for 500 points with no gaps sampled at 1 h interval, embedding dimension d=4, and embedding delays ranging between 2 and 40 h.

Figure 6 Jensen–Shannon complexity and permutation entropy for S₍₁₎ (black), D₍₁₎ (pink) for 500 points with no gaps sampled at 1 h interval, embedding dimension d=4, and embedding delays ranging between 2 and 40 h.

The bottom left panels in Figs. 5 and 6 show the value of the permutation entropy for AL, am, D₍₁₎ and S₍₁₎ as a function of embedding delay. Similarly, the bottom right panels show the value of the Jensen–Shannon complexity for AL, am, D₍₁₎ and S₍₁₎ as a function of embedding delay. What we notice is that all four signals are highly stochastic since the normalized permutation entropy is very close to 1. However, we see that the Jensen–Shannon complexity for S₍₁₎ is of comparable magnitude as for am, and it is significantly larger than for AL. This is not a surprise because the construction of S₍₁₎ was based on am, and the Jensen–Shannon complexity indicates that the former preserved the correlated structures of the latter on timescales ranging between a few hours to a few days. The top left panel of Figs. 5 and 6 shows the complexity–entropy plane, and the top right panel is a zoom of the right corner where most of the data for AL, am, and S₍₁₎ lies. In both figures the blue line curves represent the maximum and minimum value of complexity for a fixed entropy value, and the dashed curve represents the complexity–entropy mapping of fractional Brownian motion (fBm) with Hurst exponent ranging between 0 and 1, which is a stochastic process that also contains correlated structures. The fBm curve is a boundary between deterministic (above) and stochastic (below) fluctuations. We note that AL is effectively stochastic, whereas am and S₍₁₎ lie above the fBm boundary for a few tens of hours. The explanation for this behavior from am lies in its construction: it is repeated for 3 h at a time. Hence, ordinal patterns of size d=4 and embedding delays of a few hours will register the repetition as correlated structures. Since S₍₁₎ is constructed in part with am, it also contains part of its correlated structure.

For longer embedding delays of the order of seasonal variations ranging from 27 to 45 d (not shown), all four time series overlap and are indistinguishable from stochastic fluctuations. In terms of complexity–entropy plane it translates into a permutation entropy of approximately 1 and a Jensen–Shannon complexity of almost 0. Hence, the system variable S₍₁₎, based on various magnetospheric indices, preserves the stochastic and correlational structures of its individual components. The comparable values of the permutation entropy (and therefore Jensen–Shannon complexity) with the system variable with the indices for long times are not fortuitous. The permutation entropy is invariant under any monotonic transformations (for instance, if one scales the time series by a positive real number, or if one takes the logarithm). However, if one used a linear combination of non-monotonic functions, for instance trigonometric functions, then the permutation entropy would not be invariant. Since the Jensen–Shannon complexity is a function of the permutation entropy, it is also invariant under monotonic transformations. Additionally, if one takes an average around the mean of some time series over a time TAU, one will reduce the noise level for fluctuations with timescales less than TAU. Thus, the stochastic nature of the signal will be reduced, and the permutation entropy and Jensen–Shannon complexity would move up in the plane towards the chaotic and/or periodic regions. The equivalence mapping of the information theoretic measures for the system variables and geomagnetic indices is a consequence of the monotonic transformation linking the former to the latter and the absence of coarse-graining of the indices.

3.2 The secondary modes of reaction represented by S₍₂₎ and S₍₃₎

In Fig. 7a and b the second and third scalar pairs are plotted, S₍₂₎ as a function of D₍₂₎ and S₍₃₎ as a function of D₍₃₎, respectively (Eqs. 2 and 3). The correlation coefficient for the second pair is still quite high (r_corr=0.775) but lower than that of the first pair (Fig. 1). This correlation coefficient r_corr=0.775 for the secondary mode is better than correlations obtained in most studies of solar-wind–magnetosphere coupling using single measures of the magnetospheric system (e.g., Table 3 of Newell et al., 2007; Table 1 of Borovsky, 2013). D₍₂₎ describes $r_{\mathrm{corr}}^{\mathrm{2}}=\mathrm{60.0}$ % of the variance of S_(2). In Fig. 7b the correlation coefficient for the third pair is low (r_corr=0.456); D₍₃₎ only describes $r_{\mathrm{corr}}^{\mathrm{2}}=\mathrm{20.8}$ % of the variance of S₍₃₎. Canonical pairs beyond the third pair have even weaker correlations.

Figure 7 For the 1 h resolution 1991–2007 data set, the aggregate system scalar S₍₂₎ is plotted as a function of the driver scalar D₍₂₎ in panel (a) and the aggregate system scalar S₍₃₎ is plotted as a function of the driver scalar D₍₃₎ in panel (b). Each black point is 1 h of data.

Figure 2c shows that mode S₍₂₎ (Fig. 7a) is dominated by opposite-signed coefficients for mP_i and mP_e, which respectively are measures of the global ion precipitation into the atmosphere versus the global electron precipitation into the atmosphere. In this S₍₂₎ mode the intensity of ion and electron precipitation reacts oppositely. Figure 2d shows that D₍₂₎ (the driver of S₍₂₎) is dominated by the solar wind number density n_sw opposite to the clock angle sin ²(θ_clock∕2) of the solar-wind magnetic field, with the solar wind density increasing while the clock angle decreases resulting in more ion precipitation and less electron precipitation. This ion-versus-electron precipitation mode is a newly uncovered mode of reaction of the magnetosphere–ionosphere–thermosphere system to the solar wind.

Figure 2e shows that S₍₃₎ (Fig. 7b) is characterized by the polar cap index (PCI) and mP_i acting oppositely to K_p and am. PCI is a measure of high-latitude electrical currents in the magnetosphere, and mP_i is a measure of high-latitude ion precipitation; K_p and am are measures of global magnetospheric convection. This S₍₃₎ mode is very similar to a high-latitude versus convection mode uncovered by Borovsky (2014) and by Holappa et al. (2014). Figure 2f indicates that the driver D₍₃₎ for this mode is the solar wind velocity acting oppositely to the magnetic field clock angle: the wind velocity increasing while the clock angle is reduced producing more convection and less high-latitude activity, or the wind slowing down while the clock angle increases producing less convection and more high-latitude activity.

4 Advantages of the reduced (aggregate-variable) representation of the system

The aggregate variable S₍₁₎ acts as a global activity index for the magnetospheric system: S₍₁₎ is new and unfamiliar, and experience using S₍₁₎ is needed to gain an understanding of the full usefulness of this measure. S₍₁₎ could be thought of as a next-generation magnetospheric index. In Earth systems science global aggregate variables are familiar: the global warming index (Hasselman, 1997; Haustein et al., 2016), the global mean sea level (Vermeer and Rahmstorf, 2009), the mean global temperature (Hansen et al., 2006), the Palmer Drought Severity Index (Wells et al., 2004), and Sea Surface Temperature indices (Kaplan et al., 1998). And stock-market aggregate indices are well known and are more-meaningful gauges of an economy than the price of a single stock (e.g., Pan and Mishra, 2018). Here the aggregate variable S₍₁₎ is mathematically derived. The individual variables of S that go into the definition of S₍₁₎ represent familiar and identifiable aspects of activity in the magnetospheric system. The composite variable S₍₁₎ is a mix of these understood measurements, the mix reflecting some global properties of the system's reaction to the solar wind. Unfamiliar as it is, the composite-scalar D₍₁₎→S₍₁₎ reduction of the state-vector D→S picture exhibits some outright advantages for the magnetospheric system. This is particularly true in comparison with the standard method of analysis of magnetospheric driving by the solar wind that uses only a single measurement of magnetospheric activity and a single function of solar-wind variables. Four advantages are discussed in the following four paragraphs.

• Linearity. The plotted points in Fig. 1 indicate that there is a linear response of the composite system variable S₍₁₎(t) to the composite driver D₍₁₎(t). Usually, single measures of the magnetosphere tend to have a nonlinear response to the solar wind (e.g., Voros, 1994; Valdivia et al., 1996; Sharma et al., 2005b; Borovsky, 2013; Stepanova and Valdivia, 2016), with the individual activity variables saturating (becoming anomalously weak) when solar-wind driving becomes strong (e.g., Fig. 3 of Reiff and Luhmann, 1986; Fig. 17 of Lavraud and Borovsky, 2008; Fig. 6 of Borovsky, 2013). Such a saturation is not seen for S₍₁₎ driven by D₍₁₎. Undoubtedly, the linearity of the result is in part owed to the maximizing of the “linear” correlation coefficient in the CCA process. The linearity of the S₍₁₎-vs.-D₍₁₎ relation has a great advantage: the same mathematical relationship between S₍₁₎ and D₍₁₎ (i.e., Eq. 4) holds for weak driving of the system (small D₍₁₎) (e.g., Kerns and Gussenhoven, 1990) and for strong driving of the system (large D₍₁₎) (e.g., Sharma and Veeramani, 2011).

• Low noise. The high correlation between S₍₁₎ and D₍₁₎ (cf. Fig. 1) indicates that there is a relatively low level of noise in the linear-regression fit to S₍₁₎: the activity in the system as described by S₍₁₎ responds directly to the solar-wind driving as described by D₍₁₎. For example, the unaccounted-for variance of S₍₁₎ is only 15.2 %. Single measures of the magnetospheric system have much weaker Pearson linear correlation coefficients with solar-wind variables than S₍₁₎ and D₍₁₎ do. Examples can be found in Table 3 of Newell et al. (2007) and Table 1 of Borovsky (2013): the maximum correlation coefficient in those tables is 0.860 (for the D_st index), but usually it is much lower. The lower noise is also confirmed by the Jensen–Shannon complexity analysis of S₍₁₎: the points for S₍₁₎ and D₍₁₎ sit closer to the maximum complexity curve than AL and other indices. The lower noise (and higher r_corr) reduces “regression dilution bias” (Bock and Petersen, 1975; Hutcheon et al., 2010) when the system activity is fit by the driver strength. Regression dilution bias can lead to spurious interpretation of trends in the data when subsets of the data are compared, particularly when a subset with systematically weaker driving is compared with a subset with systematically stronger driving.

• High prediction efficiency. In magnetospheric physics, predicting what the reaction of the magnetospheric system will be to measured upstream solar-wind conditions is very important: i.e., the prediction of “space weather” (Singer et al., 2001). The high correlation between S₍₁₎ and D₍₁₎ means that there will be a high prediction efficiency when the value of S₍₁₎ is predicted from knowledge of the value of D₍₁₎. Note that this is a high prediction of S₍₁₎(t) without using past values of S₍₁₎, just using the present value of D₍₁₎(t). By optimizing the Pearson linear correlation coefficient between S and D, S₍₁₎ was designed to focus on aspects of the magnetospheric system that are responsive to the conditions of the solar wind. Internal dynamics of the system that are not dependent on the time-varying state of the driver are de-emphasized in S₍₁₎.

• Compactness of the description. Reductionist analysis has concluded that the magnetosphere–ionosphere–thermosphere system is extremely complex (e.g., Siscoe, 2011; Eastwood et al., 2015; Borovsky and Valdivia, 2018). As it is driven by the solar wind, there are major outstanding issues as to how the system functions (e.g., Denton et al., 2016). Having a single scalar variable S₍₁₎(t) that describes a universal global reaction of the system to its driver promises to yield insight as to how the combined system operates.

• Uncovering new modes of reaction. In the CCA analysis of the system and driver state vectors, two additional aggregate variables S₍₂₎(t) and S₍₃₎(t) were generated (Eqs. 2a and 3a). Analysis in Sect. 3.2 showed these two variables to represent two modes of reaction of the system to the driver that are independent of (uncorrelated with) the global-activity mode represented by S₍₁₎(t). The mode represented by S₍₃₎ is known (having been independently discovered by this CCA methodology in Borovsky, 2014, and by a principle-components methodology in Holappa et al., 2014), but the mode represented by S₍₂₎ has until now been unknown. The CCA methodology used here also identifies the aggregate driver variable that drives each of the independent modes. In the future, expanding the system state vector to include a larger number of measurements in the diverse magnetospheric system should enable this state-vector-reduction methodology to uncover more unknown modes of reaction of the system to the driver. Once a reaction and its driver are uncovered, reductionist analysis can be applied to determine the physical reasons why the mode arises.

For a system measured by multiple time-dependent variables (which are collected into a time-dependent system state vector S(t)), with that system driven by multiple time-dependent factors (inputs) (which are collected into a time-dependent driver state vector D(t)), canonical correlation analysis (CCA) can be used to reduce the D(t)→S(t) state-vector picture to a D_(i)(t)→S_(i)(t) composite-scalar picture. The reduction will work, even if there is influence on the driver by the system (i.e., D(t)↔S(t)). The advantageous properties of this reduction that were examined for the magnetospheric system should apply to systems in general.

Future developments of this methodology will focus on the introduction of time lags between the driver and the system, on the introduction of integro-differential correlations rather than algebraic correlations (e.g., Borovsky, 2017), and on the use of dynamic canonical correlation analysis (e.g., Dong and Qin, 2018a, b).

Appendix A: The variables comprising the magnetospheric system state vector and the solar-wind driver state vector

The time-dependent variables of the magnetospheric system state vector and the solar-wind driver state vector are listed in Table 1.

The magnetospheric variables measure various aspects of activity in the magnetosphere. The auroral upper index AU (Davis and Sugiura, 1966) measures the electrical current in the high-latitude ionosphere: this variable is taken to be a measure of electrical currents in the dayside magnetosphere (Goertz et al., 1993). The auroral lower index AL (Davis and Sugiura, 1966) measures the electrical current in the high-latitude nightside ionosphere: this variable is taken to be a measure of auroral activity in the nightside magnetosphere (Goertz et al., 1993). The polar cap index PCI is a measure of the strength of cross-polar-cap electrical current in the ionosphere (Troshichev et al., 1988). The planetary K index K_p is a measure of the strength of global convection in the magnetosphere (Thomsen, 2004). The range index am (Mayaud, 1980) is another measure of the strength of global magnetospheric convection. The disturbance storm-time index D_st measures plasma pressure in the inner magnetosphere (Dessler and Parker, 1959); D_st also reacts to the currents on the dayside boundary of the magnetosphere and to the cross-magnetotail currents in the nightside magnetosphere. The time derivative of the magnitude of the D_st index $d\left|D_{\mathrm{st}}\right|/\mathrm{d}t$ is a compound measure of magnetospheric activity: when $d\left|D_{\mathrm{st}}\right|/\mathrm{d}t$ is positive, hot plasma is being convected from the magnetotail into the dipolar portion of the magnetosphere, and when $d\left|D_{\mathrm{st}}\right|/\mathrm{d}t$ is negative, convection has recently subsided. The variables mP_e and mP_i are estimates of the full-Earth power in magnetospheric electron precipitation into the atmosphere and magnetospheric ion precipitation into the atmosphere (Emery et al., 2008, 2009), with these estimates coming from observations on only a few spacecraft in orbit around the Earth. The average of the ion-plasma-sheet particle pressure P_ips around the Earth (Borovsky, 2017) is obtained from three to five spacecraft.

The variables going into the solar-wind driver state vector are various measures of the time-dependent solar wind at Earth. The solar wind speed v_sw ranges from 244 to 1045 km s⁻¹ in the 1991–2007 data set. The solar wind number density ranges from 0.3 to 98.2 particles cm⁻³ in the data set. B_z is the magnetic-field component in the solar-wind plasma that is approximately aligned with the Earth's magnetic-dipole orientation. The function f(M) (Borovsky and Birn, 2014) is a function of the solar-wind Mach number M that accounts for the properties of the bow shock that forms upstream of the Earth in the supersonic solar-wind flow. The clock angle θ_clock measures the angular alignment of the solar-wind magnetic-field vector with the Earth's magnetic-dipole orientation. The angle θ_Bn measures the orientation of the solar-wind magnetic-field vector with respect to the Sun–Earth line. F_10.7 is the 10.7 cm radio flux from the Sun, a proxy for the ionization of the upper atmosphere of the Earth by solar photons.