Issue
J. Space Weather Space Clim.
Volume 10, 2020
Topical Issue - Space Weather research in the Digital Age and across the full data lifecycle
Article Number 25
Number of page(s) 12
DOI https://doi.org/10.1051/swsc/2020026
Published online 29 June 2020

© T. Alberti et al., Published by EDP Sciences 2020

Licence Creative CommonsThis is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1 Introduction

The dynamics of the near-Earth electromagnetic environment exhibits a nonlinear character which manifests in the emergence of complex variability at a wide range of spatial and temporal scales (Tsurutani et al., 1990; Consolini, 2018). The origin of this multiscale character of the near-Earth space (the magnetosphere–ionosphere system) is related to both, the interplanetary medium variations (e.g., solar wind variability), which is a turbulent plasma medium, and the internal dynamics associated with the variability of a set of different electric currents flowing in the magnetosphere and the ionosphere (e.g., Alberti et al., 2017). This environment is a complex (with different components that may interact with each other), highly dynamic (given its multi-scale nature) and chaotic (reflected by its unpredictable behavior) nonlinear system which exhibits scale-invariant features (Vassiliadis et al., 1990; Consolini et al., 1996; Wanliss, 2005).

The dynamical configuration and variability of this complex environment can be monitored and studied by means of time series of several geomagnetic indices, which are capable of monitoring the different current systems. For instance, the Auroral Electrojet (AE) indices (Davis & Sugiura, 1966) monitor the evolution of the current systems in the high-latitude ionosphere, whose dynamical changes can be related to different source mechanisms like the dynamics of the tail current in the magnetosphere, the electric convection field driving electric currents within the polar dynamo regions, or the inherent enhancements of solar wind particles entering the magnetosphere in the high-latitude regions (e.g., Ahn et al., 1983; Daglis et al., 1994; Baker et al., 1995; Consolini et al., 1996). These electric currents are the main drivers responsible for the occurrence of high-latitude enhancements in geomagnetic field variations, such as the development of magnetospheric substorms (e.g., Kamide & Kokubun, 1996). On the other hand, low-latitude geomagnetic field variations are related to the occurrence of storms which are mainly driven by the enhancements of the magnetospheric ring current activity (e.g., Weygand & McPherron, 2006), whose variability is monitored, among others, by the SYM-H index (Iyemori, 1990; Joselyn & Tsurutani, 1990; Wanliss & Showalter, 2006; Consolini, 2018).

Nowadays, the ongoing development of our technological society results in an increased variety of both data and phenomena which need to be properly analyzed and investigated in the framework of Space Weather research (Camporeale et al., 2018). This obviously requires, on the one hand, an increase in the computational and data science resources for managing what we know as big data, and, on the other hand, the development and/or combination of both established and novel data analysis and modeling tools for advancing our understanding of the near-Earth electromagnetic environment and processes (Alberti et al., 2017; Camporeale, 2019).

The understanding of the multi-scale and complex dynamics of the magnetosphere–ionosphere (MI) system in response to the changes of the physical conditions of the solar wind and interplanetary medium is a crucial point to develop a reliable modeling of the MI dynamics, which is at the basis of every attempt for correctly forecasting the effects of a solar disturbance (e.g., Alberti et al., 2018; Camporeale, 2019). In other words, the investigation of the nature of the MI dynamics is a central issue in Space Weather research. For instance, questions like what are the features of the MI system at different spatio-temporal scales, or to what extent the dynamics of the solar wind driving linearly affects the MI response and/or triggers internal processes, have to be considered central for a correct modeling of several Space Weather processes and phenomena, such as geomagnetic storms and substorms (e.g., Camporeale et al., 2018; Consolini et al., 2018).

In this paper, we demonstrate that a combination of two state-of-the-art nonlinear data analysis methods – empirical mode decomposition (EMD) and recurrence analysis (RA) – applied to the SYM-H index (Iyemori, 1990; Alberti et al., 2018) can provide a better understanding of the multi-scale dynamical properties of a geomagnetic storm (as opposed to periods of magnetospheric quiescence). First, the EMD is used to decompose the recorded SYM-H time series into a finite number of intrinsic components at different timescales without any specific a priori mathematical requirements (e.g., Huang et al., 1998; Alberti et al., 2017, 2018). Then, RA is applied to each derived intrinsic component to investigate different aspects of the dynamical complexity of temporal fluctuations captured in those different variability modes. Specifically, representative measures originating from two complementary quantitative concepts are used for this purpose – recurrence quantification analysis (RQA) (Marwan et al., 2007) and recurrence network analysis (RNA) (Donner et al., 2010; Zou et al., 2019). Although both, EMD and RQA/RNA have been recently employed in studies on geomagnetic variability (Alberti et al., 2017, 2018; Mendes et al., 2017; Donner et al., 2018, 2019), the combination of these concepts is novel and has not been used before for time series analysis in the context of space physics. Thereby, insights are to be expected that go beyond those provided by previous studies which have utilized more traditional analysis concepts.

2 Data description

We focus our attention on the low-latitude geomagnetic disturbances by analysing the SYM-H geomagnetic index (Iyemori, 1990), which is a proxy for the geomagnetic activity at mid-to-low latitudes as it monitors the symmetric (SYM) variations of the horizontal component (H) of the magnetic field recorded at mid-to-low latitude geomagnetic observatories. This index is provided at a 1 min temporal resolution and allows to investigate the dynamical changes of the magnetospheric equatorial ring current due to the solar wind parameters’ changes, as for example, the effects related to variations of the solar wind dynamic pressure and the interplanetary magnetic field orientation, which enhance plasma convection toward the inner magnetospheric regions (Wanliss & Showalter, 2006). Data can be freely retrieved from the NASA OMNI website (http://omniweb.gsfc.nasa.gov).

We have selected two different time intervals during the months of August and September 2018 according to the different mean values of the SYM-H index: a quiet period corresponding to the time interval between 1 and 10 August, and a disturbed storm period between 24 August and 3 September, as reported in Figure 1. During the quiet period the average value of the SYM-H index was ~3 nT reaching a minimum value of −16 nT, while during the disturbed period the average and minimum values were −52 nT and −206 nT, respectively. The latter contains the geomagnetic storm that occurred on 26 August, which was characterized by a value of Kp = 7 (https://www.spaceweatherlive.com/en/auroral-activity/top-50-geomagnetic-storms/year/2018), ranking it as the third strongest geomagnetic storm of the current solar cycle and the strongest storm of the year 2018. It was related to a large-scale quiet-region solar filament located around the central meridian, which gradually erupted around 19:00 UT and produced a partial halo CME detected by SoHO/LASCO around 21:12 UT (see http://www.stce.be/newsletter/pdf/2018/STCEnews20180831.pdf). This structure was originally expected to have only minor consequences on the near-Earth electromagnetic environment; however, its effects were not minor. Indeed, the associated interplanetary coronal mass ejection (ICME), which impacted the Earth’s magnetosphere on 26 August, was characterized by a strong southward component of the interplanetary magnetic field which persisted throughout most of the day, and by an increase of the plasma density, thus enabling reconnection and causing the geomagnetic storm observed by ground-based geomagnetic observatories (see http://www.stce.be/newsletter/pdf/2018/STCEnews20180831.pdf).

thumbnail Fig. 1

The SYM-H time series of the period considered in this work. The green and red lines refer to the quiet and the disturbed period, respectively. The dashed gray line is used as a reference to SYM-H = 0 nT.

During the disturbed period the Earth’s magnetosphere was initially strongly compressed by the impact of the interplanetary (IP) shock associated with the ICME, which caused the sudden impulse (SI) before the storm main phase (Joselyn & Tsurutani, 1990), with the SYM-H index reaching the positive value of 27 nT (see Fig. 1). Then, a sudden decrease of the SYM-H values was observed, reaching the minimum value of −206 nT on 26 August at 07:11 UT during the main phase of the geomagnetic storm associated with the southward turning of the interplanetary magnetic field B z $ {B}_z$ component ( B z $ {B}_z$ ~ −15 nT). This main phase lasted for many hours during which the equatorial ring current became intensified, and was followed by a very long recovery phase of several days during which the SYM-H index moved towards normal quiet-time values (see Fig. 1).

3 Methods

3.1 Empirical mode decomposition

In Huang et al. (1998) developed a new tool for decomposing time series with a fully data-adaptive approach, allowing to reduce mathematical assumptions and artifacts. The new tool, known as empirical mode decomposition (EMD), is based on an algorithmic procedure through which the decomposition basis can be derived, while assuming neither stationarity nor linearity, directly from the local properties of an observed signal s ( t ) $ s(t)$ (Alberti et al., 2018). Thus, the basis is not fixed a priori and formed by a finite set of empirical modes, called intrinsic mode functions (IMFs), satisfying two properties: (i) the number of zero crossings and local extrema must be equal (or differ at most by one) and (ii) at any point, the sum of the upper and lower envelopes should be zero. The IMFs are derived by the so-called sifting process, which consists of the following steps:

  1. Derive a zero-mean signal γ ( t ) = s ( t ) - s ( t ) $ \gamma (t)=s(t)-\left\langle s(t)\right\rangle$, with $ \left\langle \dots \right\rangle$ being the time average.

  2. Identify the local extrema of γ ( t ) $ \gamma (t)$ (i.e., { t [ 1 , T ] : d γ ( t ) d t | t = t = 0 } ) $ \{{t}^{\star }\in [1,T]:\frac{\mathrm{d}\gamma (t)}{\mathrm{d}t}{|}_{t={t}^{\star }}=0\})$.

  3. Interpolate local maxima (minima) by using a cubic spline interpolation to derive the upper (lower) envelope u ( t ) $ u(t)$ ( l ( t ) $ l(t)$).

  4. Evaluate the mean envelope m ( t ) = u ( t ) + l ( t ) 2 $ m(t)=\frac{u(t)+l(t)}{2}$.

  5. Evaluate the detail h ( t ) = γ ( t ) - m ( t ) $ h(t)=\gamma (t)-m(t)$.

The main purpose of the sifting process is to eliminate riding waves from time series. It requires that the sifting process must be repeated several times until an IMF is obtained and stored as c k ( t ) = h ( t ) $ {c}_k(t)=h(t)$ (Huang et al., 1998). Then, the steps from 1 to 5 must be repeated on the signal s ( t ) - c k ( t ) $ s(t)-{c}_k(t)$ until a new IMF is found, iterating this process until no more IMFs can be obtained. As suggested by Huang et al. (1998), in order to guarantee that the IMFs retain enough physical meaning of both amplitude and frequency modulations, a stopping criterion needs to be defined, taking into account to limit the size of the standard deviation between two consecutive sifting results, e.g., σ = t = 1 T | h j - 1 ( t ) - h j ( t ) | 2 h j - 1 2 ( t ) < σ c , $$ \sigma =\sum_{t=1}^T \frac{|{h}_{j-1}(t)-{h}_j(t){|}^2}{{h}_{j-1}^2(t)} < {\sigma }_c, $$(1)with σ c $ {\sigma }_c$ being typically between 0.2 and 0.3 (Huang et al., 1998). For our analysis a threshold value σ c = 0.3 $ {\sigma }_c=0.3$ is selected. At the end, we can write s ( t ) = k = 1 N c k ( t ) + r ( t ) , $$ s(t)=\sum_{k=1}^N {c}_k(t)+r(t), $$(2)where r ( t ) $ r(t)$ is the final residue of the decomposition, a monotonic function from which no more empirical modes (i.e., IMFs) can be extracted (Huang et al., 1998). Thus, a completely adaptive and a posteriori decomposition method is built, particularly suitable to investigate various aspects of both known and unknown physical processes in the framework of nonlinear geosciences (e.g., Huang et al., 1999; Huang & Wu, 2008; De Michelis et al., 2012; Consolini et al., 2017).

Another interesting property of the EMD is that each IMF can be written as an oscillating function with time-dependent amplitude a k ( t ) $ {a}_k(t)$ and phase φ k ( t ) $ {\phi }_k(t)$, respectively, by using the so-called Hilbert–Huang Transform (HHT) (Huang et al., 1998). Indeed, the HHT allows us to derive both a k ( t ) $ {a}_k(t)$ and φ k ( t ) $ {\phi }_k(t)$ by introducing a complex signal z k ( t ) = a k ( t ) e i φ k ( t ) = c k ( t ) + i   H { c k ( t ) } $ {z}_k(t)={a}_k(t){e}^{i{\phi }_k(t)}={c}_k(t)+i\enspace \mathcal{H}\{{c}_k(t)\}$ where H { c k ( t ) } = π - 1   P 0 c k ( t ) t - t d t $$ \mathcal{H}\{{c}_k(t)\}={\pi }^{-1}\enspace \mathcal{P}{\int }_0^{\mathrm{\infty }} \frac{{c}_k({t}^\mathrm{\prime})}{t-{t}^\mathrm{\prime}}\mathrm{d}{t}^\mathrm{\prime} $$(3)is the Hilbert Transform of the k-th IMF and P $ \mathcal{P}$ is the Cauchy principal value (e.g., Huang et al., 1998; Carbone et al., 2018). Thus, we derive a k ( t ) = c k ( t ) 2 + H { c k ( t ) } 2 $$ {a}_k(t)=\sqrt{{c}_k(t{)}^2+\mathcal{H}\{{c}_k(t){\}}^2} $$(4) φ k ( t ) = tan - 1 [ H { c k ( t ) } c k ( t ) ] $$ {\phi }_k(t)={\mathrm{tan}}^{-1}\left[\frac{\mathcal{H}\{{c}_k(t)\}}{{c}_k(t)}\right] $$(5)such that c k ( t ) = a k ( t ) cos [ φ k ( t ) ] $ {c}_k(t)={a}_k(t)\mathrm{cos}\left[{\phi }_k(t)\right]$ (Huang et al., 1998; Huang & Wu, 2008).

Having derived the instantaneous phase φ k ( t ) $ {\phi }_k(t)$, its time-derivative allows us to introduce the concept of instantaneous frequency ω k ( t ) = d φ k ( t ) d t $ {\omega }_k(t)=\frac{\mathrm{d}{\phi }_k(t)}{\mathrm{d}t}$ (e.g., Huang et al., 1998), from which a typical timescale of oscillation τ k = 2 π ω k ( t ) - 1 $ {\tau }_k=2\pi \left\langle {\omega }_k(t){\right\rangle^{-1}$, where $ \left\langle \dots \right\rangle$ again refers to the time average, can be derived for each IMF (e.g., Huang et al., 1998; Alberti et al., 2019).

Thus, some limitations of classical decomposition methods are overcome by the EMD, such as the requirement of a priori assumptions on the decomposition basis or the use of fixed scale eigenfunctions not allowing for a suitable description of nonstationarities. Obviously, some outstanding features need to be properly addressed in order to deal with specific problems like end or boundary effects and to avoid misleading results. For more details, we refer the reader to previous works (e.g., Huang et al., 1999; Huang & Wu, 2008; Wu & Huang, 2009; Alberti et al., 2018).

3.2 Recurrence analysis

Our nonlinear characterization of the time series provided by the different IMFs of the EMD is based on the concept of recurrences in the state space of a dynamical system. Since univariate time series commonly cannot fully capture the dynamical complexity of a higher-dimensional complex system (in our case, the externally driven Earth’s magnetosphere), additional unobserved components are first qualitatively reconstructed by the concept of time delay embedding (Takens, 1981). Here, we replace each value c k ( t ) $ {c}_k(t)$ of a given IMF by a pattern of m values associated with distinct points in time, which are mutually separated by some time delay d, i.e. c k ( t ) = ( c k ( t ) , c k ( t - d ) , c k ( t - 2 d ) , ) $ {{c}}_k(t)=({c}_k(t),{c}_k(t-d),{c}_k(t-2d),\dots )$. This so-called reconstructed state vector or embedding vector does not only characterize the present value of the time series, but rather the dynamical evolution of the system culminating in the associated observation.

The two parameters of the described embedding process, embedding dimension m and delay d, can be selected by some standard criteria like the false nearest neighbor approach and average mutual information method (Donner et al., 2018). In the present analysis, a trade-off is taken into consideration for both parameters.

On the one hand, a choice of m = 3 balances between the possibly high-dimensional nature of the observed dynamics especially in the faster EMD modes (which would call for some large m) and the fact that suitable statistics can only be expected up to a given maximum embedding dimension depending on the number of available state vectors (requiring a low m). Specifically, since we will close our present study with some sliding window analysis with a window size of only 200 min (see Sect. 4.5), using a larger embedding dimension could lead to some parts of the reconstructed state space not being sufficiently populated by state vectors to allow for a proper statistical characterization. Although we may expect different intrinsic dimensionality of different IMFs, we will therefore keep this value fixed throughout all analyses.

On the other hand, the intrinsic timescales of variability differ strongly among the different IMFs, so that we select the embedding delays individually for each IMF according to the first zero crossing of the associated autocorrelation function as a generally accepted criterion in nonlinear time series analysis (see Table 1). As we will show in Section 4.1, each IMF is associated with a relatively narrow range of instantaneous timescales, so that this first root is always of the order one fourth of the mean instantaneous timescale as expected for oscillatory signals. Notably, although the choice of d may well affect the recurrence characteristics of a time series quantitatively, the general results typically do not change qualitatively under minor variations of this embedding parameter.

Table 1

Embedding delays d $ d$ (in min) for each IMF during the quiet and storm period, respectively. Only the first 13 IMFs are studied in terms of their recurrence properties, since the remaining four do not exhibit sufficient dynamics over the respective study periods to allow for a meaningful statistical analysis.

A more detailed discussion on how to choose embedding parameters for applications to geomagnetic activity time series can be found in the Supplementary Material of Donner et al. (2018).

In the abstract metric space spanned by the coordinates of the aforementioned embedding vectors, the recurrence of a state vector is defined according to the spatial closeness of two vectors corresponding to different points in time (Eckmann et al., 1987; Marwan et al., 2007). As a common criterion for proximity, we select the 5% mutually closest pairs of state vectors among all pairs of such vectors. That is, a proximity threshold ε ( k ) $ {\epsilon }^{(k)}$ is applied to all pairwise distances between state vectors obtained by embedding the k $ k$-th IMF that corresponds to the empirical 5th percentile of all those distances. The location of recurrent state vectors in terms of their associated observation times is summarized in the binary recurrence matrix, accordingly defined as R ij ( k ) ( ε ( k ) ) = Θ ( ε ( k ) - c k ( t i ) - c k ( t j ) ) , $$ {R}_{{ij}}^{(k)}({\epsilon }^{(k)})=\mathrm{\Theta }({\epsilon }^{(k)}-\Vert {{c}}_k({t}_i)-{{c}}_k({t}_j)\Vert ), $$(6)where Θ ( ) $ \mathrm{\Theta }(\dots )$ denotes the Heaviside function and $ \Vert \dots \Vert $ is the maximum norm.

The dynamical patterns encoded in this recurrence matrix can be characterized in several complementary ways. On the one hand, recurrence quantification analysis (RQA) makes use of the length distributions of diagonal (and vertical) line structures arising in the graphical representation of the recurrence matrix (the so-called recurrence plot; Eckmann et al., 1987; Marwan et al., 2007), emphasizing that diagonal lines indicate the parallel evolution of two trajectory segments over a considerable amount of time, which reflects a certain degree of predictability of the dynamics. Here, we focus on the so-called degree of determinism (DET), defined as the fraction of recurrent state vectors contained in diagonal line structures consisting of at least two points, and the diagonal line length entropy (ENT), which is simply the Shannon entropy of the frequency distribution of observed diagonal line lengths (Marwan et al., 2007).

In addition to the line-based RQA measures that capture temporal dependency structures, we also employ two characteristics from recurrence network analysis (RNA), which reinterprets the recurrence matrix equation (6) as the adjacency matrix of a random geometric graph in the metric space of the embedding vectors (Donner et al., 2010; Zou et al., 2019). Since network characteristics are invariant under time-reordering of the state vectors, the recurrence networks solely capture geometric properties of the observed dynamical system’s trajectory in the reconstructed state space. Nonetheless, RNA measures have been found recently to sensitively distinguish between geomagnetic storm and quiescence phases (as reflected in the disturbance storm-time index (Dst) (Donner et al., 2018), which is closely related to the SYM-H index studied in the present work). Here, we restrict ourselves to two RNA measures, the network transitivity ( T $ \mathcal{T}$) (which is directly related to a generalized version of a fractal dimension; Donner et al., 2011) and the average path length ( L $ \mathcal{L}$). For details on both measures, we refer to Donner et al. (2015) and Zou et al. (2019).

4 Results and discussion

4.1 Timescale variability of the SYM-H index

As a first step, the timescale variability of the SYM-H index during the whole study period (e.g., between 1 August and 1 October, 2018) is investigated by applying the EMD procedure. A set of N = 17 empirical modes (IMFs) is derived whose characteristic mean timescales τ k $ {\tau }_k$ and associated standard deviations σ τ k $ {\sigma }_{{\tau }_k}$, obtained from instantaneous frequencies (e.g., σ τ k $ {\sigma }_{{\tau }_k}$ = std { 2   π   ω k - 1 ( t ) } $ \{2\enspace \pi \enspace {\omega }_k^{-1}(t)\}$), are reported in Table 2.

Table 2

Characteristic timescales τ k $ {\tau }_k$ and standard deviations σ τ k $ {\sigma }_{{\tau }_k}$ for each IMF estimated for the whole time period between 1 August and 1 October, 2018.

The SYM-H index variability clearly shows a well-defined multiscale variability over a wide range of timescales, e.g. from few minutes up to several days. The short-timescale modes are related to the fast dynamics of the SYM-H index, which is mainly due to the internal dynamical state of the near-Earth electromagnetic environment (e.g., Kamide & Kokubun, 1996; Consolini & DeMichelis, 2005), while the long-timescale variability is related to the occurrence of geomagnetic storms and to the effects of the interplanetary medium variability on the magnetosphere-ionosphere system (e.g., Alberti et al., 2017). These features can be clearly evidenced by looking at the behavior of some selected IMFs at different timescales (i.e., τ 3 ~ 12 min, τ 7 ~ 110 min, τ 11 ~ 1.5 days, and τ 15 ~ 13 days) and the residue r ( t ) $ r(t)$ of the decomposition reported in Figure 2. As expected, empirical modes typically show amplitude enhancements during the geomagnetic storm (see the red line in Fig. 2), although irregular amplitude variations are also found during the quiet period (see the green line in Fig. 2), especially in the short timescale IMFs (see, for example, the behavior of c 7 ( t ) $ {c}_7(t)$ in comparison with that of c 11 ( t ) $ {c}_{11}(t)$). These features completely support the existence of different dynamical components characterizing the SYM-H index variability: a fast component, which is less affected by the solar wind variability, and a slow component, which is strongly coupled with solar wind variations (Alberti et al., 2017, 2018).

thumbnail Fig. 2

The SYM-H time series (upper panel), the empirical modes c 3 ( t ) $ {c}_3(t)$, c 7 ( t ) $ {c}_7(t)$, c 11 ( t ) $ {\mathrm{c}}_{11}(t)$, and c 15 ( t ) $ {c}_{15}(t)$ (middle panels), and the residue r ( t ) $ r(t)$ (lower panel). The green and red lines refer to the quiet and the disturbed period, respectively. The dashed gray line is used as a reference to SYM-H = 0 nT.

The different multiscale characteristics of SYM-H are also highlighted by looking at the behavior of the energy content during both quiet and disturbed periods, the latter being mainly controlled by the solar wind variability as opposed to the former (e.g., Baker et al., 1995; Wanliss, 2005; Consolini, 2018). This can be simply done by evaluating the mean variances of all empirical modes during the two different time intervals separately, e.g. during the quiet and the disturbed period. For this purpose, we also evaluated the mean timescales during both periods separately by averaging the inverse of the instantaneous frequencies ω k ( t ) $ {\omega }_k(t)$ over the two different time intervals. The obtained mean timescales do not significantly differ from those evaluated by considering the whole period of two months, especially for the first IMFs associated with the short-term variability of the SYM-H index, thus suggesting that the observed multiscale variability is a common feature of both quiet and disturbed periods. Moreover, by looking at the variance-timescale distribution of empirical modes (see Fig. 3), we clearly note that a larger energy content (e.g., larger variance) is stored in the long-timescale IMFs during the disturbed period as compared with the quiet one (up to one order of magnitude), while a similar energy content is found at timescales below 200 min during both the quiet and the disturbed periods. This implies that the role of the long-term geomagnetic variability significantly increases during disturbed periods suggesting a close connection with the solar wind variability associated with the effects of a solar disturbance on the near-Earth electromagnetic environment (Alberti et al., 2017).

thumbnail Fig. 3

The variance-timescale distribution of the different empirical modes ( σ k 2 $ {\sigma }_k^2$ vs. τ k $ {\tau }_k$). The green and red symbols refer to the quiet and the disturbed period, respectively. The mean timescales are evaluated by considering only the quiet or the disturbed time interval, respectively.

4.2 Dynamical complexity of the overall SYM-H variability

In order to further characterize the differences between the short- and long-term variability of the SYM-H index, we investigated their complexity by firstly looking at the behavior of the resulting recurrence plots for two different empirical modes, i.e., c 7 ( t ) $ {c}_7(t)$ and c 11 ( t ) $ {c}_{11}(t)$, as reported in Figure 4, during both the quiet and the storm periods, respectively. It is evident that the IMF capturing faster timescales (e.g., c 7 ( t ) $ {c}_7(t)$) exhibits a much richer dynamics than c 11 ( t ) $ {c}_{11}(t)$, resulting in the presence of diagonal and vertical line structures arranged in larger blocks. While diagonal lines indicate a certain degree of predictability, with the typical line length corresponding to the prediction horizon with respect to a maximum tolerable prediction error determined by the recurrence threshold ε $ \epsilon $, vertical lines indicate slowly varying states and, hence, temporary laminar dynamics. Taking both together, the emerging large-scale block structures are commonly indicative of multistable and/or intermittent dynamics (Marwan et al., 2007). The latter interpretation matches very well with the known presence of intermittency not only in the solar wind as the main external driver of magnetospheric dynamics, but also geomagnetic activity itself (e.g., Bruno et al., 2003; Wanliss et al., 2004; Vörös et al., 2005), which is expected to manifest particularly in the fast component of geomagnetic variability (e.g., Alberti et al., 2018; Consolini et al., 2018).

thumbnail Fig. 4

An example of recurrence plots for two empirical modes c 7 ( t ) $ {c}_7(t)$ and c 11 ( t ) $ {c}_{11}(t)$. The green and red colors refer to the quiet and the disturbed period, respectively.

Although visually appearing less densely populated than for c 7 ( t ) $ {c}_7(t)$, the recurrence plots for c 11 ( t ) $ {c}_{11}(t)$ exhibit the same total number of recurring pairs of state vectors. This can be understood when zooming into the respective plots (not shown): while the recurrence patterns for c 7 ( t ) $ {c}_7(t)$ have complex fine-scale structures, the slower ones form rather thick diagonal lines especially during the quiet period. This behavior originated from the slow variability encoded in c 11 ( t ) $ {c}_{11}(t)$, necessarily appearing more regular and predictable due to its clearer oscillatory character (cf. Fig. 2), also calling for the emergence of long diagonal line structures in the recurrence plot. At the same time, the overall amplitude of variations is rather large as compared to the amplitude differences of subsequent points, which implies that subsequent state vectors (appearing next to each other in the recurrence plot) are very likely similar and, hence, most often identified as recurrent pairs. Hence, the more regular appearance of the obtained recurrence plot as compared to c 7 ( t ) $ {c}_7(t)$ basically reflects the considerably higher degree of smoothness and associated larger decorrelation time of c 11 ( t ) $ {c}_{11}(\mathrm{t})$.

4.3 Dynamical complexity during quiet and disturbed periods

While the overall differences between the recurrence plots for different IMFs can be explained in a straightforward manner as discussed above, comparing the recurrence structures of the same mode during quiet and disturbed periods is more informative. Despite the fact that IMFs are data-adaptive features and as such can obey marked changes in the mean instantaneous frequency as time proceeds, the qualitative differences of the shown plots between both considered time intervals can not only be explained by a change in the intrinsic frequency, but also reflect marked amplitude changes. This is particularly evident for c 11 ( t ) $ {c}_{11}(t)$, for which the secondary diagonal lines (beside the main diagonal) visible during the quiet period almost completely vanish along with the storm interval. Instead, we find an entirely different morphology of the recurrence plot structures involving bended recurrence structures especially along with the recovery phase after the storm. While the absence of recurrences over a large portion of the recurrence plot simply reflects the uniquely large amplitude of the IMF during parts of this time interval, bended lines are commonly indicative of changing frequencies of oscillations (Marwan et al., 2007), which are also clearly visible in Figure 2.

Conversely, the differences between the quiet and the disturbed period for c 7 ( t ) $ {c}_7(t)$ are more subtle, with the quiet period being characterized by block structures in the recurrence plot that are more pronounced than the storm one. This observation suggests that the effect of intermittency is dominant in the quiet-time magnetosphere, which could appear contradictory. However, due to the complex nature of the magnetosphere (Vassiliadis et al., 1990; Consolini, 2018) this result needs to be interpreted in terms of the concurrent effects of fast and slow dynamical processes. Indeed, a geomagnetic storm can be understood as a coherent large-amplitude pattern embedded in a multiscale chaotic signal which reflects in the slow dynamical component, with the faster ones carrying less information on relevant dynamical processes. Thus, since fast dynamical processes are dominated by the slow ones during the main and the recovery phases of a geomagnetic storm, this potentially reflects in a higher degree of observational stochasticity (e.g., due to observing a spatially extended system only at some selected sites) relative to what is contained in the same IMF during magnetospheric quiescence.

4.4 Recurrence measures at different timescales: quiet vs. storm periods

The dynamical patterns of the SYM-H index variability at different timescales can be characterized by evaluating some recurrence measures for each empirical mode, thus providing an independent way to investigate the existence of a timescale separation with respect to previous works (e.g., Alberti et al., 2017, 2018; Consolini et al., 2018). For this purpose, we estimated the degree of determinism, the network transitivity, the diagonal line length entropy, and the average shortest path length for each empirical mode as shown in Figure 5 during both the quiet and the storm period. As highlighted in Figure 5, all measures manifest a clear scale-dependent behavior.

thumbnail Fig. 5

Recurrence characteristics of the different IMFs of the SYM-H index in dependence on the associated average instantaneous periodicity of each mode: degree of determinism (upper left), network transitivity (upper right), diagonal line length entropy (lower left), and average shortest path length (lower right). The green and red symbols refer to the quiet and the disturbed period, respectively.

First of all, the degree of determinism exhibits a steep rise to values close to 1 (which are expected for a regular, completely deterministic dynamics). Notably, those values are already reached for timescales of the order of 100 min. This behavior can be considered typical for signals with pronounced oscillations that are densely sampled (i.e., in our case for a sampling interval of 1 min for oscillations with an instantaneous periodicity of 100 min or even far longer). In such situations, the degree of determinism loses its ability to sensitively distinguish between more or less predictable dynamics, as it is based on line structures in the recurrence plot, which naturally appear in case of slow oscillations even if those were completely unpredictable. Note that predictability needs to be understood as relative to the typical timescales of variability embedded in a given signal. In this regard, RQA provides a methodology that is intrinsically tailored to studying short-term fluctuations. As a consequence, the timescale of 100 min at which DET saturates should not be confused with the scale separation between fast and slow modes, since it most likely depends on the time resolution of the underlying data.

Turning to the diagonal line length entropy as a second RQA measure, we note that this property exhibits a steady increase as the timescales increase. In order to understand this, it is important to note that this entropy is conceptually different from other common entropic characteristics used in nonlinear time series analysis, exhibiting increases if the observed dynamics becomes more regular (Letellier, 2006). Other variants of recurrence plot based entropy measures have been developed recently to account for this conceptual ambiguity (Letellier, 2006; Eroglu et al., 2014; Corso et al., 2018), but have not yet become standard tools in applications of RQA and therefore not been used in the present work.

As an interesting feature, we find that the recurrence network transitivity generally increases with increasing IMF timescale, but only very slowly. Recalling the interpretation of this measure as a proxy for the fractal dimension of the system, this implies that the corresponding dimensionality of the different IMF components does not change markedly over a range of time scales up to the order of, say, one day. Only for the slowest modes, transitivity increases more steeply, indicating a successive reduction in dimensionality.

Concerning the aforementioned three recurrence characteristics, it is remarkable that the overall behavior with increasing timescale does not exhibit marked differences between quiet and disturbed periods. This is, however, distinctively different for the average path length of the obtained recurrence networks. Similar to network transitivity and diagonal line length entropy, this measure also first increases with increasing timescale, before suddenly changing its behavior at time scales of the order of 100 min (for the quiet period) and 1000 min (for the storm period), respectively (see Fig. 5, lower right panel).

According to the former observation, the average path length is the only studied recurrence measure that reveals a clear difference between quiet and storm periods when considering timescales beyond 100 min. Specifically, while for the disturbed period the average shortest path length continues to increase with rising timescale up to a scale comparable with the typical duration of a magnetic storm, the quiet period exhibits a transition near τ 7 ~ 110 min, beyond which the dynamics is characterized by an almost constant value of the average path length of about 5. A likely explanation for this qualitative change would be the decomposition of the recurrence network into two or more mutually disconnected parts (which could be avoided by considering a higher number of recurrent pairs of state vectors than the 5% used here), which prohibits the further rise of long graph distances among the nodes. In general, unless such a splitting of the network takes place, a successive increase in average path length with rising smoothness of the signal would be the expected behavior. This is because for an oscillatory dynamics, the topology of the network would increasingly resemble that of a ring structure. Notably, for a given number of nodes and links in a network as in the present case, such an almost regular ring structure would maximize the average path length (except for structures resembling a linear lattice with marked end nodes).

We emphasize that especially the network transitivity does not seem to show strong differences between quiet and storm periods when considering its variation with the different IMFs. This does not necessarily contradict the results of previous studies for the Dst index (Donner et al., 2018, 2019) which have demonstrated that this measure can sensitively distinguish between magnetospheric variability in both types of conditions. The former studies have notably focused on the integral variability aggregating over all time scales, while the present analysis explicitly accounts for the different scales separately by utilizing the corresponding IMFs. Besides having used different yet closely related indices, an absence of differences in network transitivity for individual IMFs but presence in the aggregated variability could point to changes in the interdependence structure between different variability modes, a hypothesis that should be further studied in future work.

The presented analysis has necessarily been limited by the choice of only four recurrence characteristics. While the selection made here has been mainly guided by previous successful applications and the complementary nature of the four measures, we do not rule out that other statistical properties of the IMFs’ recurrence structures could more sensitively trace differences between quiet and storm periods. Future in-depth studies should further examine this idea, along with the consideration of different quiet and storm periods.

4.5 Network transitivity for fast and slow dynamical components

The above results support the previous findings of a timescale separation between the fast and the slow dynamical components of the SYM-H variability (Alberti et al., 2017). In order to examine the corresponding dynamics in more detail, we finally reconstructed both the fast and the slow component by considering partial sums of equation (2) involving empirical modes with timescales below (fast component) and above (slow component) 200 min (e.g., Alberti et al., 2017, 2018), respectively, as reported in Figure 6. We note that the slow component clearly reproduces the large-scale features of a geomagnetic storm: it fluctuates around zero when a quiet period is considered, and is characterized by a positive value (the so-called sudden impulse; Tsurutani et al., 1990) before exhibiting negative values corresponding to the main phase of the storm. Conversely, the fast dynamics is overall characterized by fluctuations around the common values during a quiet period, with larger amplitude enhancements during the main phase of the geomagnetic storm, although significant deviations from the background level can also be observed during the quiet period (see Fig. 6, green line in the upper panel). This is clearly the reflection of the occurrence of coherent intermittent activity bursts due to loading-unloading dynamics of the Earth’s magnetotail (Kamide & Kokubun, 1996; Uritsky & Pudovkin, 1998; Consolini, 2002; Consolini & De Michelis, 2005), usually occurring on timescales of 100–200 min, characterizing the internal state of the Earth’s magnetosphere (Alberti et al., 2017, 2018).

thumbnail Fig. 6

Reconstruction of the fast (upper panel) and slow (lower panel) components of SYM-H based on the obtained empirical modes. The green and red lines refer to the quiet and the disturbed period, respectively. The dashed gray line is used as a reference to SYM-H = 0 nT.

To investigate the dynamical changes of recurrence properties in both fast and slow SYM-H components, during the quiet and the disturbed periods, we finally evaluated the resulting network transitivity following a sliding window approach (Donner et al., 2018, 2019; Lekscha & Donner, 2019), using embedding delays of 100 and 200 min for the fast and slow component, respectively. The results thus obtained are reported in Figure 7. What is interesting here is that we observe more dynamical anomalies (i.e., periods with transitivity values that do not match the expected range of values when simply taking random state vectors from the time series for constructing a recurrence network) in the quiet than in the storm period (without attempting here any physical interpretation of the individual highlighted episodes, which appears challenging without considering additional characteristics of the coupled magnetosphere–ionosphere–solar wind system). This at first apparently surprising result can be explained by the origin of processes involved in the fast and the slow dynamics. Indeed, the magnetosphere is continuously interacting with the solar wind, which determines its shape but not exactly its internal configuration. Thus, all the currents inside the magnetosphere are stable also if the solar wind would not exist. During a storm, a large number of solar wind (and solar, e.g., suprathermal) particles enter the magnetosphere, thereby changing the intensity of the current systems, transferring more information to these currents by directly driving their changes and configurations. During a quiet period, the configuration and dynamics is less driven and, hence, less regular, since it is mostly the result of the information flow between the different current systems.

thumbnail Fig. 7

Recurrence network transitivity T $ \mathcal{T}$ estimated for sliding windows in time with a width of 200 min for the fast (top) and slow (bottom) SYM-H components as described in the text, using data from the quiet (left, green) and storm period (right, red). The grey horizontal shades indicate the respective ranges of expected values of T $ \mathcal{T}$ if 200 state vectors had been randomly drawn from the whole series for generating a recurrence network.

Again, we note that an additional consideration of further complementary recurrence characteristics could potentially help better understanding the observed differences in the frequency of dynamical anomalies between quiet and storm period. Although it might at first glance appear discouraging that individual anomalies in the network transitivity cannot be directly related with marked disruptive events affecting the Earth’s magnetosphere, our findings rather motivate further in-depth studies to uncover the dynamical differences between magnetic field variations associated with “normal” and “extraordinary” transitivity values along with their possible (internal or external) origins.

5 Summary and perspectives

In this paper, we have dealt with the characterization of the nonlinear recurrence properties of geomagnetic variability during periods characterized by a different geomagnetic disturbance level by using a timescale dependent approach. By a combined use of a nonlinear decomposition technique – the EMD – and two flavors of RA, we have been able to disentangle the different features of both quiet and disturbed periods. Specifically, we have found that,

  1. Previous findings by Alberti et al. (2017) on the existence of two characteristic ranges of timescales, which have been identified as fast and slow dynamics (Alberti et al., 2018), are confirmed by both recurrence quantification analysis and recurrence network analysis.

  2. Both fast and slow components are characterized by different recurrence properties, with more temporary dynamical anomalies during the quiet than the storm period.

The above results can be interpreted in the framework of the complex and nonlinear interaction between the solar wind and the Earth’s magnetosphere as a signature of the different origins of processes occurring at different timescales, due to different sources located in the interplanetary and near-Earth environment (Kamide & Kokubun, 1996; Consolini & De Michelis, 2005; Alberti et al., 2017). Indeed, the more regular behavior, which can be observed in terms of its recurrence properties as well as other nonlinear characteristics (e.g., Donner et al., 2018, and references therein), during the disturbed period is a direct reflection of the strong degree of dependency of geomagnetic activity on solar wind variability (Alberti et al., 2017). This implies that the magnetospheric response to interplanetary changes can be quite reasonably forecasted during a disturbed period in terms of an enhancement of convection processes, as well as, that the slow component, i.e. the dynamics at timescales longer than 200 min, is characterized by a lower degree of complexity, which is an indication of the external origin of processes at these timescales (Alberti et al., 2017; Consolini, 2018). Conversely, the more complex and irregular behavior at short timescales, i.e., below 200 min, can be related to impulsive energy releases occurring in the tail central plasma sheet as well as in the turbulent state of the tail neutral sheet, being the fingerprint of turbulent energy releases associated with the unloading mechanisms and less affected by the external changes of interplanetary conditions (Alberti et al., 2017; Consolini, 2018).

The results of this investigation show that more efforts are needed to correctly deal with the complex nature of the near-Earth electromagnetic environment. As also shown by Consolini et al. (2018), the short-term dynamics is characterized by a high-dimensional, low-predictability behavior, while the long-term one shows a large forecast horizon (~50 min). This can be interpreted as a sort of topological phase transition (Chang et al., 2003) confirmed by the behavior of the recurrence network transitivity (see Fig. 5) revealing a dimensionality reduction when approaching large timescales.

Moreover, this manuscript offers new insights for Space Weather purposes of forecasting the behavior of geomagnetic indices, not only in terms of their response to the interplanetary medium changes. Indeed, our analysis, in a complementary way with that based on traditional fractal measures (e.g., Consolini et al., 2018), allows to strengthen the need of providing high-resolution proxies of the internal dynamical state of the near-Earth electromagnetic environment for correctly dealing with the characterization of fast dynamical processes, generally related to the explosive dynamics of the tail current (Consolini, 2002; Vörös et al., 2005). The typical features of the recurrence plots for the short-term component during both the quiet and the disturbed periods point towards the primary role of temporal laminar dynamical features which are associated with multistable and/or intermittent dynamics of the magnetotail plasma (Wanliss et al., 2004; Vörös et al., 2005). This type of dynamics has profound consequences on reliable estimates of several Space Weather phenomena, such as, for example, the generation of ground-induced currents (Tozzi et al., 2019a, 2019b).

Future studies will be devoted to a statistical characterization of the multiscale nature of geomagnetic indices, the dynamical interplay between variations at different timescales, the comparison of a larger variety of recurrence-based measures applied to different geomagnetic indices monitoring both magnetospheric and ionospheric current systems (e.g., ASY-H, AE, AU, AL), as well as the investigation of higher-resolution time series which could help in a deeper characterization of the short-term fluctuations of geomagnetic variability.

Acknowledgments

The OMNI data were obtained from the GSFC/SPDF OMNIWeb interface at http://omniweb.gsfc.nasa.gov. The authors also acknowledge World Data Center for Geomagnetism (Kyoto) for making available the geomagnetic index data. All the derived data products presented in this paper can be provided by the authors upon request. TA acknowledges the Institute for Space Astrophysics and Planetology (IAPS) for its support within the framework “Bando Nuove Idee IAPS 2019”. GC and PDM thanks the financial support by Italian MIUR-PRIN grant 2017APKP7T on “Circumterrestrial Environment: Impact of Sun-Earth Interaction”. RVD and JL have been financially supported by the German Federal Ministry for Education and Research (BMBF) via the Young Investigators Group CoSy-CC2 (grant no. 01LN1306A) and the German Academic Scholarship Foundation. All recurrence analyses reported in this work have been obtained using the Python package pyunicorn (Donges et al., 2015), which is freely available at GitHub. TA, GC and RVD acknowledge fruitful discussions within the scope of the International Team “Complex Systems Perspectives Pertaining to the Research of the Near-Earth Electromagnetic Environment” at the International Space Science Institute in Bern, Switzerland. The editor thanks two anonymous reviewers for their assistance in evaluating this paper.

References

Cite this article as: Alberti T, Lekscha J, Consolini G, De Michelis P & Donner RV. 2020. Disentangling nonlinear geomagnetic variability during magnetic storms and quiescence by timescale dependent recurrence properties. J. Space Weather Space Clim. 10, 25.

All Tables

Table 1

Embedding delays d $ d$ (in min) for each IMF during the quiet and storm period, respectively. Only the first 13 IMFs are studied in terms of their recurrence properties, since the remaining four do not exhibit sufficient dynamics over the respective study periods to allow for a meaningful statistical analysis.

Table 2

Characteristic timescales τ k $ {\tau }_k$ and standard deviations σ τ k $ {\sigma }_{{\tau }_k}$ for each IMF estimated for the whole time period between 1 August and 1 October, 2018.

All Figures

thumbnail Fig. 1

The SYM-H time series of the period considered in this work. The green and red lines refer to the quiet and the disturbed period, respectively. The dashed gray line is used as a reference to SYM-H = 0 nT.

In the text
thumbnail Fig. 2

The SYM-H time series (upper panel), the empirical modes c 3 ( t ) $ {c}_3(t)$, c 7 ( t ) $ {c}_7(t)$, c 11 ( t ) $ {\mathrm{c}}_{11}(t)$, and c 15 ( t ) $ {c}_{15}(t)$ (middle panels), and the residue r ( t ) $ r(t)$ (lower panel). The green and red lines refer to the quiet and the disturbed period, respectively. The dashed gray line is used as a reference to SYM-H = 0 nT.

In the text
thumbnail Fig. 3

The variance-timescale distribution of the different empirical modes ( σ k 2 $ {\sigma }_k^2$ vs. τ k $ {\tau }_k$). The green and red symbols refer to the quiet and the disturbed period, respectively. The mean timescales are evaluated by considering only the quiet or the disturbed time interval, respectively.

In the text
thumbnail Fig. 4

An example of recurrence plots for two empirical modes c 7 ( t ) $ {c}_7(t)$ and c 11 ( t ) $ {c}_{11}(t)$. The green and red colors refer to the quiet and the disturbed period, respectively.

In the text
thumbnail Fig. 5

Recurrence characteristics of the different IMFs of the SYM-H index in dependence on the associated average instantaneous periodicity of each mode: degree of determinism (upper left), network transitivity (upper right), diagonal line length entropy (lower left), and average shortest path length (lower right). The green and red symbols refer to the quiet and the disturbed period, respectively.

In the text
thumbnail Fig. 6

Reconstruction of the fast (upper panel) and slow (lower panel) components of SYM-H based on the obtained empirical modes. The green and red lines refer to the quiet and the disturbed period, respectively. The dashed gray line is used as a reference to SYM-H = 0 nT.

In the text
thumbnail Fig. 7

Recurrence network transitivity T $ \mathcal{T}$ estimated for sliding windows in time with a width of 200 min for the fast (top) and slow (bottom) SYM-H components as described in the text, using data from the quiet (left, green) and storm period (right, red). The grey horizontal shades indicate the respective ranges of expected values of T $ \mathcal{T}$ if 200 state vectors had been randomly drawn from the whole series for generating a recurrence network.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.