Magnetospheric chaos and dynamical complexity response during storm time disturbance

In this study, we examine the magnetospheric chaos and dynamical complexity response to the disturbance storm time (D_st) and solar wind electric field (VB_s) during different categories of geomagnetic storm (minor, moderate and major geomagnetic storm). The time series data of the D_st and VB_s are analysed for a period of 9 years using non-linear dynamics tools (maximal Lyapunov exponent, MLE; approximate entropy, ApEn; and delay vector variance, DVV). We found a significant trend between each non-linear parameter and the categories of geomagnetic storm. The MLE and ApEn values of the D_st indicate that chaotic and dynamical complexity responses are high during minor geomagnetic storms, reduce at moderate geomagnetic storms and decline further during major geomagnetic storms. However, the MLE and ApEn values obtained from VB_s indicate that chaotic and dynamical complexity responses are high with no significant difference between the periods that are associated with minor, moderate and major geomagnetic storms. The test for non-linearity in the D_st time series during major geomagnetic storm reveals the strongest non-linearity features. Based on these findings, the dynamical features obtained in the VB_s as input and D_st as output of the magnetospheric system suggest that the magnetospheric dynamics are non-linear, and the solar wind dynamics are consistently stochastic in nature.

1 Introduction

The response of chaos and dynamical complexity behaviour with respect to magnetospheric dynamics varies (Tsurutani et al., 1990). This is due to changes in the interplanetary electric fields imposed on the magnetopause and those penetrating the inner magnetosphere and sustaining convection, thereby initiating a geomagnetic storm (Dungey, 1961; Pavlos et al., 1992). A prolonged southward turning of interplanetary magnetic field (IMF, B_z), which indicates that solar-wind–magnetosphere coupling is in progress, was confirmed on many occasions for which such a geomagnetic storm was driven by co-rotating interaction regions (CIRs), by the sheath preceding an interplanetary coronal mass ejection (ICME), or by a combination of the sheath and an ICME magnetic cloud (Gonzalez and Tsurutani, 1987; Tsurutani and Gonzalez, 1987; Tsurutani et al., 1988; Cowley, 1995; Tsutomu, 2002; Yurchyshyn et al., 2004; Kozyra et al., 2006; Echer et al., 2008; Meng et al., 2019; Tsurutani et al., 2020; Echer et al., 2006). The sporadic magnetic reconnection between the southward component of the Alfvén waves and the earth's magnetopause leads to isolated substorms or convection events such as the high-intensity long-duration continuous AE activity (HILDCAA, where AE represents auroral electrojet) which are shown to last from days to weeks (Akasofu, 1964; Tsurutani and Meng, 1972; Meng et al., 1973; Tsurutani and Gonzalez, 1987; Hajra et al., 2013; Liou et al., 2013; Mendes et al., 2017; Hajra and Tsurutani, 2018; Tsurutani and Hajra, 2021; Russell, 2001). Notably, the introduction of the disturbance storm time (D_st) index (Sugiura, 1964; Sugiura and Kamei, 1991) unveiled a quantitative measure of the total energy of the ring current particles. Therefore, the D_st index remains one of the most popular global indicators that can precisely reveal the severity of a geomagnetic storm (Dessler and Parker, 1959).

The D_st fluctuations exhibit different signatures for different categories of geomagnetic storm. Ordinarily, one can easily anticipate that fluctuations in a D_st signal appear chaotic and complex. These may arise from the changes in the interplanetary electric fields driven by the solar-wind–magnetospheric coupling processes. At different categories of geomagnetic storm, fluctuations in the D_st signals differ (Oludehinwa et al., 2018). One obvious reason is that as the intensity of the geomagnetic storm increases, the fluctuation behaviour in the D_st signal becomes more complex and non-linear in nature. It has been established that the electrodynamic response of the magnetosphere to solar wind drivers are non-autonomous in nature (Price and Prichard, 1993; Price et al., 1994; Johnson and Wings, 2005). Therefore, the chaotic analysis of the magnetospheric time series must be related to the concept of input–output dynamical process (Russell et al., 1974; Burton et al., 1975; Gonzalez et al., 1989, 1994). Consequently, it is necessary to examine the chaotic behaviour of the solar wind electric field (VB_s) as input signals and the magnetospheric activity index (D_st) as output during different categories of geomagnetic storms.

Several works have been presented on the chaotic and dynamical complexity behaviour of the magnetospheric dynamics based on an autonomous concept, i.e. using the time series data of magnetospheric activity alone such as auroral electrojet (AE), amplitude lower (AL) and D_st index (Vassiliadis et al., 1990; Baker and Klimas, 1990; Vassiliadis et al., 1991; Shan et al., 1991; Pavlos, 1994; Klimas et al., 1996; Valdivia et al., 2005; Mendes et al., 2017; Consolini, 2018). Authors found evidence of low-dimensional chaos in the magnetospheric dynamics. For instance, the report by Vassiliadis et al. (1991) shows that the computation of Lyapunov exponent for AL index time series gives a positive value of Lyapunov exponent, indicating the presence of chaos in the magnetospheric dynamics. Unnikrishnan (2008) studied the deterministic chaotic behaviour in the magnetospheric dynamics under various physical conditions using AE index time series and found that the seasonal mean value of Lyapunov exponent in winter season during quiet periods (0.7±0.11 min⁻¹) is higher than that of the stormy periods (0.36±0.09 min⁻¹). Balasis et al. (2006) examined the magnetospheric dynamics in the D_st index time series from pre-magnetic storm to magnetic storm period using fractal dynamics. They found that the transition from anti-persistent to persistent behaviour indicates that the occurrence of an intense geomagnetic storm is imminent. Balasis et al. (2009) further reveal the dynamical complexity behaviour in the magnetospheric dynamics using various entropy measures. They reported a significant decrease in dynamical complexity and an accession of persistency in the D_st time series as the magnetic storm approaches. Recently, Oludehinwa et al. (2018) examined the non-linearity effects in D_st signals during minor, moderate and major geomagnetic storms using recurrence plots and recurrence quantification analyses. They found that the dynamics of the D_st signal is stochastic during minor geomagnetic storm periods and deterministic as the geomagnetic storm increases.

Also, studies describing the solar wind and magnetosphere as a non-autonomous system have been extensively investigated. Price et al. (1994) examine the non-linear input–output analysis of AL index and different combinations of interplanetary magnetic field (IMF) with solar wind parameters as input functions. They found that only a few of the input combinations show any evidence whatsoever for non-linear coupling between the input and output for the interval investigated. Pavlos et al. (1999) presented further evidence of magnetospheric chaos. They compared the observational behaviour of the magnetospheric system with the results obtained by analysing different types of stochastic and deterministic input–output systems and asserted that a low-dimensional chaos is evident in magnetospheric dynamics. Prabin Devi et al. (2013) studied the magnetospheric dynamics using AL index and the southward component of IMF (B_z). They observed that the magnetosphere and turbulent solar wind have values corresponding to non-linear dynamical system with chaotic behaviour. The modelling and forecasting approach have been applied to magnetospheric time series using non-linear models (Valdivia et al., 1996; Vassiliadis et al., 1999; Vassiliadis, 2006; Balikhin et al., 2010). These efforts have improved our understanding that the concept of non-linear dynamics can reveal some hidden dynamical information in the observational time series. In addition to these non-linear effects in D_st signals, a measure of the exponential divergence and convergence within the trajectories of a phase space known as maximal Lyapunov exponent (MLE), which has the potential to depict the chaotic behaviour in the D_st and VB_s time series during a minor, moderate or major geomagnetic storm, has not been investigated. In addition, to the best of our knowledge, computation of approximate entropy (ApEn) that depicts the dynamical complexity behaviour during different categories of geomagnetic storm has not been reported in the literature. The test for non-linearity through delay vector variance (DVV) analysis that reveals the non-linearity features in D_st and VB_s time series during minor, moderate and major geomagnetic storms is not well known. It is worthy to note that understanding the dynamical characteristics in the D_st and VB_s signals at different categories of geomagnetic storms will provide useful diagnostic information to different conditions of space weather phenomenon. Consequently, this study attempts to carry out comprehensive numerical analyses to unfold the chaotic and dynamical complexity behaviour in the D_st and VB_s signals during minor, moderate and major geomagnetic storms. In Sect. 2, our methods of data acquisition are described. Also, the non-linear analysis that we employed in this investigation is detailed. In Sect. 3, we unveil our results and engage in the discussion of results in Sect. 4.

2 Description of the data and non-linear dynamics

The D_st index is derived by measurements from ground-based magnetic stations at low-latitude observatories around the world and depicts mainly the variation of the ring current, as well as the Chapman–Ferraro magnetopause currents, and tail currents to a lesser extent (Sugiura, 1964; Feldstein et al., 2005, 2006; Love and Gannon, 2009). Due to its global nature, the D_st time series provides a measure of how intense a geomagnetic storm was (Dessler and Parker, 1959). In this study, we considered D_st data for the period of 9 years from January to December between 2008 and 2016 which were downloaded from the World Data Centre for Geomagnetism, Kyoto, Japan (http://wdc.kugi.kyoto-u.ac.jp/dst_final/index.html, last access: 15 May 2021). We use the classification of geomagnetic storms as proposed by Gonzalez et al. (1994) such that D_st index value in the ranges $\mathrm{0}\le D_{\mathrm{st}}\le -\mathrm{50}$ nT, $-\mathrm{50}\le D_{\mathrm{st}}\le -\mathrm{100}$ nT and $-\mathrm{100}\le D_{\mathrm{st}}\le -\mathrm{250}$ nT are classified as minor, moderate and major geomagnetic storms, respectively, and each time series is being classified based on its minimum D_st value. The solar wind electric field (VB_s) data are archived from the National Aeronautics and Space Administration, Space Physics Facility (https://omniweb.gsfc.nasa.gov/form/dx1.html, last access: 15 May 2021). The sampling time of D_st and VB_s time series data was 1 h. It is well known that the dynamics of the solar wind contribute to the driving of the magnetosphere (Burton et al., 1975). Furthermore, we took the solar wind electric field (VB_s) as the input signal (Price and Prichard, 1993; Price et al., 1994). The VB_s was categorized according to the periods of minor, moderate and major geomagnetic storms. Then, the D_st and VB_s time series were subjected to a variety of non-linear analytical tools explained as follows.

2.1 Phase space reconstruction and observational time series

An observational time series can be defined as a sequence of scalar measurements of some quantity, which is a function of the current state of the system taken at multiples of a fixed sampling time. In non-linear dynamics, the first step in analysing an observational time series is to reconstruct an appropriate state space of the system. Takens (1981) and Mañé (1981) stated that one time series or a few simultaneous time series are converted to a sequence of vectors. This reconstructed phase space has all the dynamical characteristic of the real phase space provided the time delay and embedding dimension are properly specified.

$$\begin{array}{}\text{(1)}& X\left(t\right)={\left[x\left(t\right),x\left(t+\mathit{\tau }\right),x\left(t+\mathrm{2}\mathit{\tau }\right),\mathrm{\dots },x\left(t+\left(m-\mathrm{1}\right)\mathit{\tau }\right)\right]}^T,\end{array}$$

where X(t) is the reconstructed phase space, x(t) is the original time series data, τ is the time delay and m is the embedding dimension. An appropriate choice of τ and m is needed for the reconstruction of phase space, which is determined by average mutual information and false nearest neighbour, respectively.

2.2 Average mutual information (AMI)

The method of average mutual information (AMI) is one of the non-linear techniques used to determine the optimal time delay (τ) required for phase space reconstruction in observational time series. The time delay mutual information was proposed by Fraser and Swinney (1986) instead of an autocorrelation function. This method takes into account non-linear correlations within the time series data. It measures how much information can be predicted about one time series point, given full information about the other. For instance, the mutual information between x_i and x_(i+τ) quantifies the information in state x_(i+τ) under the assumption that information at the state x_i is known. The AMI for a time series x(t_i), $i=\mathrm{1},\mathrm{2},\mathrm{\dots },N$ is calculated as

$$\begin{array}{}\text{(2)}& \begin{array}{rl}I\left(T\right)& =\sum _{x\left(t_i\right),x(t_i+T)}^{.}P\left(x\left(t_i\right),x\left(t_i+T\right)\right)\\ & \times {\log }_{\mathrm{2}}\left[{\displaystyle \frac{P\left(x\left(t_i\right),x\left(t_i+T\right)\right)}{P\left(x\left(t_i\right)\right)P\left(x\left(t_i+T\right)\right)}}\right],\end{array}\end{array}$$

where x(t_i) is the ith element of the time series, T=kΔt ($k=\mathrm{1},\mathrm{2},\mathrm{\dots },k_{max}$), P(x(t_i)) is the probability density at x(t_i), and $P\left(x\right(t_i),x(t_i+T\left)\right)$ is the joint probability density at the pair xt(t_i) and x(t_i+T). The time delay (τ) of the first minimum of AMI is chosen as optimal time delay (Fraser and Swinney, 1986; Fraser, 1986). Therefore, the AMI was applied to the D_st and VB_s time series, and the plot of AMI versus time delay is shown in Fig. 3. We notice that the AMI showed the first local minimum at roughly τ=15 h. Furthermore, the values of τ near this value of ∼ 15 h maintain constancy for both VB_s and D_st. In the analysis, τ=15 h was used as the optimal time delay for the computation of maximal Lyapunov exponent.

2.3 False nearest neighbour (FNN)

In determining the optimal choice of embedding dimension (m), the false nearest neighbour method was used in the study. The method was suggested by Kennel et al. (1992). The concept is based on how the number of neighbours of a point along a signal trajectory changes with increasing embedding dimension. With increasing embedding dimension, the false neighbour will no longer be neighbours; therefore, by examining how the number of neighbours changes as a function of dimension, an appropriate embedding dimension can be determined. For instance, suppose we have a one-dimensional time series. We can construct a time series y(t) of D-dimensional points from the original one-dimensional time series x(t) as follows:

$$\begin{array}{}\text{(3)}& y\left(t\right)=\left(x\left(t\right),x(t+\mathit{\tau }\right),\mathrm{\dots },x\left(t+\left(D-\mathrm{1}\right)\mathit{\tau }\right),\end{array}$$

where τ and D are time delay and embedding dimension. Using the formula from Kennel et al. (1992) and Wallot and Monster (2018), if we have a D-dimensional phase space and denote the rth nearest neighbour of a coordinate vector y(t) by y^(r)(t), then the square of the Euclidean distance between y(t) and the rth nearest neighbour is

$$\begin{array}{}\text{(4)}& R_D^{\mathrm{2}}\left(t,r\right)=\sum _{k=\mathrm{0}}^{D-\mathrm{1}}{\left[x\left(t+k\mathit{\tau }\right)-x^{\left(r\right)}(t+k\mathit{\tau })\right]}^{\mathrm{2}}.\end{array}$$

Now applying the logic outlined above, we can go from a D-dimensional phase space to D+1-dimensional phase space by time-delay embedding, adding a new coordinate to y(t), and ask what is the squared distance between y(t) and the same rth nearest neighbour:

$$\begin{array}{}\text{(5)}& R_{D+\mathrm{1}}^{\mathrm{2}}\left(t,r\right)=R_D^{\mathrm{2}}\left(t,r\right)+{\left[x\left(t+D\mathit{\tau }\right)-x^{\left(r\right)}(t+D\mathit{\tau })\right]}^{\mathrm{2}}.\end{array}$$

As explained above, if the one-dimensional time series is already properly embedded in D dimensions, then the distance R between y(t) and the rth nearest neighbour should not change appreciably by some distance criterion R_tol (i.e. R<R_tol). Moreover, the distance of the nearest neighbour when embedded into the next higher dimension relative to the size of the attractor should be less than some criterion A_tol (i.e. $R_{D+\mathrm{1}}<A_{\mathrm{tol}}$). Doing this for the nearest neighbour of each coordinate will result in many false nearest neighbours when embedding is insufficient or in few (or no) false neighbours when embedding is sufficient. In the analysis, the FNN was applied to the D_st and VB_s time series to detect the optimal value of embedding dimension (m). Figure 4 shows a sample plot of the percentage of false nearest neighbour against embedding dimension in one of the months under investigation (other months show similar results; thus, for brevity we depict only one of the results). We notice that the false nearest neighbour attains its minimum value at m≥5, indicating that embedding dimension (m) values from m≥5 are optimal. Therefore, m=5 was used for the computation of maximal Lyapunov exponent.

2.4 Maximal Lyapunov exponent (MLE)

The maximal Lyapunov exponent (MLE) is one of the most popular non-linear dynamics tools used for detecting chaotic behaviour in a time series. It describes how small changes in the state of a system grow at an exponential rate and eventually dominate the behaviour. An important indication of chaotic neighbour of a dissipative deterministic system is the existence of a positive Lyapunov exponent. A positive MLE signifies divergence of trajectories in one direction or expansion of an initial volume in this direction. On the other hand, a negative MLE exponent implies convergence of trajectories or contraction of volume along another direction. The algorithm proposed by Wolf et al. (1985) for estimating MLE is employed to compute the chaotic neighbour of the D_st and VB_s time series at minor, moderate and major geomagnetic storms. Other methods of determining MLE include Rosenstein's method, Kantz's method and so on. In this study, the MLE for minor, moderate and major geomagnetic storm periods was computed with m=5 and τ=15 h as shown in Figs. 5 and 6 (bar plots) for D_st and VB_s. The calculation of MLE is explained as follows: given a sequence of vector x(t), an m-dimensional phase space is formed from the observational time series through an embedding theorem as

$$\begin{array}{}\text{(6)}& x\left(t\right),x\left(t+\mathit{\tau }\right),\mathrm{\dots },x\left(t+\left(m-\mathrm{1}\right)\mathit{\tau }\right),\end{array}$$

where m and τ are as defined earlier; after reconstructing the observational time series, the algorithm locates the nearest neighbour (in Euclidean sense) to the initial point $x\left(t_{\mathrm{0}}\right),\mathrm{\dots },x(t_{\mathrm{0}}+(m-\mathrm{1}\left)\mathit{\tau }\right)$ and denotes the distance between these two points L(t₀). At a later point t₁, the initial length will have evolved to length L^′(t₁). Then the MLE is calculated as

$$\begin{array}{}\text{(7)}& \mathit{\lambda }={\displaystyle \frac{\mathrm{1}}{t_M-t_{\mathrm{0}}}}\sum _{k=\mathrm{1}}^M{\log }_{\mathrm{2}}{\displaystyle \frac{L^{\prime }\left(t_k\right)}{L\left(t_{k-\mathrm{1}}\right)}},\end{array}$$

where M is the total number of replacement steps. We look for a new data point that satisfies two criteria reasonably well: its separation, L(t₁), from the evolved fiducial point is small. If an adequate replacement point cannot be found, we retain the points that were being used. This procedure is repeated until the fiducial trajectory has traversed the entire data.

2.5 Approximate entropy (ApEn)

Approximate entropy (ApEn) is one of the non-linear dynamics tools that measures the dynamical complexity in observational time series. The concept was proposed by Pincus (1991), who provide a generalized measure of regularity, such that it accounts for the logarithm likelihood in the observational time series. For instance, a dataset of length N that repeat itself for m points within a boundary will again repeat itself for m+1 points. Because of its computational advantage, ApEn has been widely used in many disciplines to study dynamical complexity (Pincus and Kalman, 2004; Pincus and Goldberger, 1994; McKinley et al., 2011; Kannathan et al., 2005; Balasis et al., 2009; Shujuan and Weidong, 2010; Moore and Marchant, 2017). The ApEn is computed using the formula below:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}\mathrm{ApEn}\left(m,r,N\right)& ={\displaystyle \frac{\mathrm{1}}{N-m+\mathrm{1}}}\sum _{i=\mathrm{1}}^{N-m+\mathrm{1}}\log C_i^m\left(r\right)\\ & -{\displaystyle \frac{\mathrm{1}}{N-m}}\sum _{i=\mathrm{1}}^{N-m}\log C_i^m\left(r\right),\end{array}\end{array}$$

where $C_i^m\left(r\right)=\frac{\mathrm{1}}{N-m+\mathrm{1}}\sum _{j=\mathrm{1}}^{N-m+\mathrm{1}}\mathrm{\Theta }\left(r-∥x_i-x_j∥\right)$ is the correlation integral, m is the embedding dimension and r is the tolerance. To compute the ApEn for the D_st and VB_s time series classified as minor, moderate and major geomagnetic storms from 2008 to 2016, we choose m=3 and τ=1 h. We refer interested readers to the works of Pincus (1991), Kannathal et al. (2005) and Balasis et al. (2009) where all the computational steps regarding ApEn were explained in detail. Figures. 5 and 6 depict the stem plot of ApEn for D_st and VB_s from 2008 to 2016.

2.6 Delay vector variance (DVV) analysis

The delay vector variance (DVV) is a unified approach in analysing and testing for non-linearity in a time series (Gautama et al., 2004; Mandic et al., 2007). The basic idea of the DVV is that if two delay vectors of a predictable signal are close to each other in terms of the Euclidean distance, they should have a similar target. For instance, when a time delay (τ) is embedded into a time series x(k), $k=\mathrm{1},\mathrm{2},\mathrm{\dots },N$, then a reconstructed phase space vector is formed which represents a set of delay vectors (DVs) of a given dimension.

$$\begin{array}{}\text{(9)}& X\left(k\right)={\left[X_{k-m\mathit{\tau }},\mathrm{\dots },X_{k-\mathit{\tau }}\right]}^T\end{array}$$

Reconstructing the phase space, a set (λ_k) is generated by grouping those DVs that are within a certain Euclidean distance to DVs (X(k)). For a given embedding dimension (m), a measure of unpredictability σ^∗2 is computed over all pairwise Euclidean distances between delay vector as

$$\begin{array}{}\text{(10)}& d\left(i,j\right)=∥x\left(i\right)-x\left(j\right)∥(i\ne j).\end{array}$$

Then, sets λ_k(r_d) are generated as the sets which consist of all delay vectors that lie closer to x(k) than a certain distance r_d.

$$\begin{array}{}\text{(11)}& {\mathit{\lambda }}_k\left(r_{\mathrm{d}}\right)=\left\{x\left(i\right)∥x\left(k\right)-x\left(i\right)∥\le r_{\mathrm{d}}\right\}\end{array}$$

For every set λ_k(r_d), the variance of the corresponding target σ^∗2(r_d) is

$$\begin{array}{}\text{(12)}& {\mathit{\sigma }}^{\ast \mathrm{2}}\left(r_{\mathrm{d}}\right)={\displaystyle \frac{\frac{\mathrm{1}}{N}{\sum }_{k=\mathrm{1}}^N{\mathit{\sigma }}_k^{\mathrm{2}}\left(r_{\mathrm{d}}\right)}{{\mathit{\sigma }}_k}},\end{array}$$

where σ^∗2(r_d) is target variance against the standardized distance, indicating that Euclidean distance will be varied in a manner standardized with respect to the distribution of pairwise distance between DVs. The iterative amplitude-adjusted Fourier transform (IAAFT) method is used to generate the surrogate time series (Kugiumtzis, 1999). If the surrogate time series yields a DV plots similar to the original time series and the scattered plot coincides with the bisector line, then the original time series can be regarded as linear (Theiler et al., 1992; Gautama et al., 2004; Imitaz, 2010; Jaksic et al., 2016). On the other hand, if the surrogate time series yields DV plot that is not similar to that of the original time series, then the deviation from the bisector lines indicates non-linearity. The deviation from the bisector lines grows as a result of the degree of non-linearity in the observational time series.

$$\begin{array}{}\text{(13)}& t^{\mathrm{DVV}}=\sqrt{〈\left({\mathit{\sigma }}^{\ast \mathrm{2}}\left(r_{\mathrm{d}}\right)-{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N{\mathit{\sigma }}_{\mathrm{s},\mathrm{i}}^{\ast \mathrm{2}}}{N_{\mathrm{s}}}}\right)〉},\end{array}$$

where ${\mathit{\sigma }}_{\mathrm{s},\mathrm{i}}^{\ast \mathrm{2}}(r_{\mathrm{d}}$) is the target variance at the span r_d for the ith surrogate. To carry out the test for non-linearity in the D_st signals (m=3, n_d=3), the number of reference DVs = 200 and the number of surrogates N_s=25 were used in all of the analysis. Then we examined the non-linearity response at minor, moderate and major geomagnetic storms.

3 Results

In this study, D_st and VB_s time series from January to December were analysed for the period of 9 years (2008 to 2016) to examine the chaotic and dynamical complexity response in the magnetospheric dynamics during minor, moderate and major geomagnetic storms. Figures 1 and 2 display the samples of fluctuation signatures of D_st and VB_s signals classified as (a) minor, (b) moderate and (c) major geomagnetic storms. The plot of average mutual information against time delay (τ) shown in Fig. 3 depicts that the first local minimum of the AMI function was found to be roughly at τ=15 h. Furthermore, we notice that the values of τ near this value of ∼ 15 h maintain constancy for both VB_s and D_st. Also, in Fig. 4, we display the plot of the percentage of false nearest neighbour against embedding dimension (m). It is obvious that a decrease in false nearest neighbour when increasing the embedding dimension drops steeply to zero at the optimal dimension (m=5); thereafter, the false neighbours stabilize at that m=5 for VB_s and D_st. Therefore, m=5 and τ=15 h were used for the computation of MLE at different categories of geomagnetic storm, while m=3 and τ=1 h are applied for the computation of ApEn values.

Figure 1 Samples of D_st signals classified as (a) minor, (b) moderate and (c) major geomagnetic storms.

Figure 2 Samples of solar wind electric fields (VB_s) during (a) minor, (b) moderate and (c) major geomagnetic storms.

Figure 3 The plot average mutual information against embedding time delay (τ).

Figure 4 The plot of percentage of false nearest neighbours against embedding dimension (m).

The results of MLE (bar plot) and ApEn (stem plot) for D_st at minor, moderate and major geomagnetic storms are shown in Fig. 5. During minor geomagnetic storms, we notice that the value of MLE ranges between 0.07 and 0.14 for most of the months classified as minor geomagnetic storms. Similarly, the ApEn (stem plot) ranges between 0.59 and 0.83. It is obvious that strong chaotic behaviour with high dynamical complexity are associated with minor geomagnetic storms. During moderate geomagnetic storms, (see Fig. 5b), we observe a reduction in MLE values (0.04–0.07) compared to minor geomagnetic storm periods. Within the observed values of MLE during moderate geomagnetic storms, we found a slight rise of MLE in the following months: March in 2008; April in 2011; January, February, and April in 2012; July, August, September, October, and November in 2015; and November in 2016. Also, the ApEn revealed a reduction in values between 0.44 and 0.57 during moderate geomagnetic storms. The lowest values of ApEn were noticed in the following months: May 2010, March 2011 and January 2016. During major geomagnetic storms as shown in Fig. 5, the minimum and maximum values of MLE are respectively 0.03 and 0.04, implying a very strong reduction of chaotic behaviour compared with minor and moderate geomagnetic storms. The lowest values of MLE were found in the months of July 2012, June 2013 and March 2015. Interestingly, further reduction in ApEn value (0.29–0.40) was as well noticed during this period. Thus, during major geomagnetic storms, chaotic behaviour and dynamical complexity subside significantly.

Figure 5 The MLE (bar plot) and ApEn (stem plot) of D_st at (a) minor, (b) moderate and (c) major geomagnetic storms.

Figure 6 The MLE (bar plot) and ApEn (stem plot) of solar wind electric field (VB_s) during: (a) minor, (b) moderate and (c) major geomagnetic storms.

We display in Fig. 6 the results of MLE and ApEn computation for the VB_s which has been categorized according to the periods of minor, moderate and major geomagnetic storms. The values of MLE (bar plot) were between 0.06 and 0.20 for VB_s. The result obtained indicates strong chaotic behaviour with no significant difference in chaoticity during minor, moderate and major geomagnetic storms. Similarly, the results obtained from computation of ApEn (stem plot) for VB_s depict a minimum value of 0.60 and peak value of 0.87 as shown in Fig. 6. The ApEn values of VB_s indicates high dynamical complexity response with no significant difference during the periods of the three categories of geomagnetic storm investigated.

The test for non-linearity in the D_st signals during minor, moderate and major geomagnetic storms was analysed through the DVV analysis. Shown in Fig. 7 is the DVV plot and DVV scatter plot during minor geomagnetic storms for January 2009 and January 2014. We found that the DVV plots during minor geomagnetic storms reveals a slight separation between the original and surrogate data. Also, the DVV scatter plots show a slight deviation from the bisector line between the original and surrogate data which implies non-linearity. Also, during moderate geomagnetic storms, we notice that the DVV plot depicts a wide separation between the original and the surrogate data. Also, a large deviation from the bisector line between the original and the surrogate data was also noticed in the DVV scatter plot as shown in Fig. 8 thus indicating non-linearity. In Fig. 9, we display samples of DVV plot and DVV scatter plot during major geomagnetic storm for October 2011 and December 2015. The original and the surrogate data showed a very large separation in the DVV plot during major geomagnetic storm. The DVV scatter plot depicts the greatest deviation from the bisector line between the original and the surrogate data which is also an indication of non-linearity. The DVV analysis of the VB_s time series during minor, moderate and major geomagnetic storms shown in Figs. 10–12 revealed a separation between the original and surrogate data with no significant difference between the periods of minor, moderate and major geomagnetic storms.

Figure 7 The DVV plot and scatter plot for D_st during minor geomagnetic storm for January 2009 and January 2014.

Figure 8 The DVV plot and scatter plot for D_st during moderate geomagnetic storm for March 2011 and January 2015.

Figure 9 The DVV plot and scatter plot for D_st during major geomagnetic storm for October 2011 and December 2015.

Figure 10 The DVV plot and scatter plot for VB_s during minor geomagnetic storm for January 2009 and January 2014.

Figure 11 The DVV plot and scatter plot for VB_s during moderate geomagnetic storm for March 2011 and January 2015.

Figure 12 The DVV plot and scatter plot for VB_s during major geomagnetic storm for October 2011 and December 2015.

4 Discussion of results

4.1 The chaotic and dynamical complexity response in D_st at minor, moderate and major geomagnetic storms

Our result shows that the values of MLE for D_st during minor geomagnetic storm are higher, indicating significant chaotic response during minor geomagnetic stormy periods (bar plot, Fig. 5). This increase in chaotic behaviour for D_st signals during minor geomagnetic storms may be as a result of asymmetry features in the longitudinal distribution of solar source region for the corotating interaction regions (CIRs) signatures responsible for the development of geomagnetic storms (Turner et al., 2006; Kozyra et al., 2006). CIR-generated magnetic storms are generally weaker than ICME- or MC-generated storms (Gonzalez et al., 1994; Tsurutani et al., 1995; Feldstein et al., 2006; Richardson and Cane, 2011). Therefore, we suspect that the increase in chaotic behaviour during minor geomagnetic storms is strongly associated with the asymmetry features in the longitudinal distribution of solar source region for the corotating interaction region (CIR) signatures. For most of these periods of moderate geomagnetic storms, the values of MLE decreases compared to minor geomagnetic storms. This revealed that as geomagnetic stormy events build up, the level of unpredictability and sensitive dependence on initial condition (chaos) begin to decrease (Lorentz, 1963; Stogaz, 1994). The chaotic behaviour during major geomagnetic storms decreases significantly compared with moderate geomagnetic storms. The reduction in chaotic response during moderate and its further declines at major geomagnetic storms may be attributed to the disturbance in the interplanetary medium driven by sheath preceding an interplanetary coronal mass ejection (ICME) or combination of the sheath and an ICME magnetic cloud (Echer et al., 2008; Tsurutani et al., 2003; Meng et al., 2019). Notably, the dynamics of the solar-wind–magnetospheric interaction are dissipative and chaotic in nature (Pavlos, 2012), and the electrodynamics of the magnetosphere due to the flux of interplanetary electric fields had a significant impact on the state of the chaotic signatures. For instance, the observation of strong chaotic behaviour during minor geomagnetic storms suggests that the dynamics was characterized by a weak magnetospheric disturbance. The reduction in chaotic behaviour at moderate and major geomagnetic storm period reveals the dynamical features with regards to when a strong magnetospheric disturbance begins to emerge. Therefore, our observation of chaotic signatures at different categories of geomagnetic storm has a potential capacity to give useful diagnostic information about monitoring space weather events. It is important to note that the features of D_st chaotic behaviour at different categories of geomagnetic storm has not been reported in the literature. For example, a previous study of Balasis et al. (2009, 2011) investigated dynamical complexity behaviour using different entropy measures and revealed the existence of low dynamical complexity in the magnetospheric dynamics and attributed it to ongoing large magnetospheric disturbance (major geomagnetic storm). The work of Balasis et al. (2009, 2011), where a certain dynamical characteristic evolved in the D_st signal was revealed, was limited to 1 year of data (2001). It is worthy to note that the year 2001, according to sunspot variations, is a period of high solar activity during solar cycle 23. It is characterized by numerous and strong solar eruptions that were followed by significant magnetic storm activities. This confirms that on most of the days in year 2001, the geomagnetic activity is strongly associated with major geomagnetic storms. The confirmation of low dynamical complexity response in the D_st signal during major geomagnetic storms agree with our current study. However, the idea of comparing the dynamical complexity behaviour at different categories of geomagnetic storms and revealing its chaotic features was not reported. This is the major reason why our present investigation is crucial to the understanding of the level of chaos and dynamical complexity involved during different categories of geomagnetic storms. As an extension to the single-year investigation done by Balasis et al. (2009, 2011) during a major geomagnetic storm, we further investigated 9 years of data of D_st that covered minor, moderate and major geomagnetic storms (see Fig. 5, stem plots) and unveiled their dynamical complexity behaviour. During major geomagnetic stormy periods, we found that the ApEn values decrease significantly, indicating reduction in the dynamical complexity behaviour. This is in agreement with the low dynamical complexity reported by Balasis et al. (2009, 2011) during a major geomagnetic period. Finally, based on the method of DVV analysis, we found that a test of non-linearity in the D_st time series during major geomagnetic storms reveals the strongest non-linearity features.

4.2 The chaotic and dynamical complexity behaviour in the VB_s as input signals

The results of the MLE values for VB_s revealed a strong chaotic behaviour during the three categories of geomagnetic storms. Comparing these MLE values during minor to those observed during moderate and major geomagnetic storms, the result obtained did not indicate any significant difference in chaoticity (bar plots, Fig. 6). Also, the ApEn values of VB_s during the periods associated with minor, moderate and major geomagnetic storms revealed high dynamical complexity behaviour with no significant difference between the three categories of geomagnetic storms investigated. These observation of high chaotic and dynamical complexity behaviour in the dynamics of VB_s may be due to interplanetary discontinuities caused by the abrupt changes in the interplanetary magnetic field direction and plasma parameters (Tsurutani et al., 2010). Also, the indication of high chaotic and dynamical complexity behaviour in VB_s signifies that the solar wind electric field is stochastic in nature. The DVV analysis for VB_s revealed non-linearity features with no significant difference between the minor, moderate and major geomagnetic storms. It is worth mentioning that the dynamical complexity behaviour for VB_s is different from what was observed for D_st time series data. For instance, our results for D_st time series revealed that the chaotic and dynamical complexity behaviour of the magnetospheric dynamics are high during minor geomagnetic storms, reduce at moderate geomagnetic storms and further decline during major geomagnetic storms. The VB_s signal revealed a high chaotic and dynamical complexity behaviour at all the categories of geomagnetic storm period. Therefore, these dynamical features obtained in the VB_s as input signal and the D_st as the output in describing the magnetosphere as a non-autonomous system further support the finding of Donner et al. (2019) that found increased or unchanged behaviour in dynamical complexity for VB_s and low dynamical complexity behaviour during storms using the recurrence method. This suggests that the magnetospheric dynamics are non-linear, and the solar wind dynamics are consistently stochastic in nature.

5 Conclusions

This work has examined the magnetospheric chaos and dynamical complexity behaviour in the disturbance storm time (D_st) and solar wind electric field (VB_s) as input during different categories of geomagnetic storms. The chaotic and dynamical complexity behaviour at minor, moderate and major geomagnetic storms for solar wind electric field (VB_s) as input and D_st as output of the magnetospheric system were analysed for the period of 9 years using non-linear dynamics tools. Our analysis has shown a noticeable trend of these non-linear parameters (MLE and ApEn) and the categories of geomagnetic storm (minor, moderate and major). The MLE and ApEn values of the D_st have indicated that the chaotic and dynamical complexity behaviour are high during minor geomagnetic storms, low during moderate geomagnetic storms and further reduced during major geomagnetic storms. The values of MLE and ApEn obtained from VB_s indicate that chaotic and dynamical complexity are high with no significant difference during the periods of minor, moderate and major geomagnetic storms. Finally, the test for non-linearity in the D_st time series during major geomagnetic storms reveals the strongest non-linearity features. Based on these findings, the dynamical features obtained in the VB_s as input and D_st as output of the magnetospheric system suggest that the magnetospheric dynamics are non-linear, and the solar wind dynamics are consistently stochastic in nature.