Disentangling nonlinear geomagnetic variability during magnetic storms and quiescence by timescale dependent recurrence properties

Understanding the complex behavior of the near-Earth electromagnetic environment is one of the main challenges of Space Weather studies. This includes both the correct characterization of the different physical mechanisms responsible for its configuration and dynamics as well as the efforts which are needed for a correct forecasting of several phenomena. By using a nonlinear multi-scale dynamical systems approach, we provide here new insights into the scale-to-scale dynamical behavior of both quiet and disturbed periods of geomagnetic activity. The results show that a scale-dependent dynamical transition occurs when moving from short to long timescales, i.e., from fast to slow dynamical processes, the latter being characterized by a more regular behavior, while more dynamical anomalies are found in the behavior of the fast component. This suggests that different physical processes are typical for both dynamical regimes: the fast component, being characterized by a more chaotic and less predictable behavior, can be related to the internal dynamical state of the near-Earth electromagnetic environment, while the slow component seems to be less chaotic and associated with the directly driven processes related to the interplanetary medium variability. Moreover, a clear difference has been found between quiet and disturbed periods, the former being more complex than the latter. These findings support the view that, for a correct forecasting in the framework of Space Weather studies, more attention needs to be devoted to the identification of proxies describing the internal dynamical state of the near-Earth electromagnetic environment.


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 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 Topical Issue -Space Weather research in the Digital Age and across the full data lifecycle 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, lowlatitude 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 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 spatiotemporal 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 methodsempirical 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., 2017Alberti et al., , 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 purposerecurrence 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 Mendes et al., 2017;Donner et al., 2018Donner et al., , 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.

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:// omniwebgsfc.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-geomagneticstorms/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 groundbased geomagnetic observatories (see http://www.stce.be/ newsletter/pdf/2018/STCEnews20180831.pdf).
T. Alberti et al.: J. Space Weather Space Clim. 2020, 10, 25 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 component (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).

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Þ . 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: 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Þ (Huang et al., 1998). Then, the steps from 1 to 5 must be repeated on the signal 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., with r c being typically between 0.2 and 0.3 (Huang et al., 1998). For our analysis a threshold value r c ¼ 0:3 is selected.
At the end, we can write where 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Þ and phase u 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Þ and u k ðtÞ by introducing a complex signal z k ðtÞ ¼ a k ðtÞe iu k ðtÞ ¼ c k ðtÞ þ i Hfc k ðtÞg where is the Hilbert Transform of the k-th IMF and P is the Cauchy principal value (e.g., Huang et al., 1998;Carbone et al., 2018). Thus, we derive such that c k ðtÞ ¼ a k t ð Þ cos u k ðtÞ ½ (Huang et al., 1998;Huang & Wu, 2008). Having derived the instantaneous phase u k ðtÞ, its timederivative allows us to introduce the concept of instantaneous frequency x k ðtÞ ¼ du k ðtÞ dt (e.g., Huang et al., 1998), from which a typical timescale of oscillation s k ¼ 2phx k ðtÞi À1 , where h. . .i 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).

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Þ 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 À 2dÞ; . . .Þ. 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.
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 e ðkÞ is applied to all pairwise distances between state vectors obtained by embedding the 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 where Hð. . .Þ denotes the Heaviside function and jj. . .jj 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 ) (which is directly related to a generalized version of a fractal dimension; Donner et al., 2011) and the average path length (L). For details on both measures, we refer to Donner et al. (2015) and Zou et al. (2019).

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 s k and associated standard deviations r s k , obtained from instantaneous frequencies (e.g., r s k = std f2 p x À1 k ðtÞg), are reported in Table 2. 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., s 3~1 2 min, s 7~1 10 min, s 11~1 .5 days, and s 15~1 3 days) and the residue 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Þ in comparison with that of 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 .
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 x 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 .

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Þ and 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Þ) exhibits a much richer dynamics than 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 e, 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).
Although visually appearing less densely populated than for c 7 ðtÞ, the recurrence plots for 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Þ 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Þ, 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Þ basically reflects the considerably higher degree of smoothness and associated larger decorrelation time of c 11 ðtÞ.

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Þ, 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Þ 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.

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., 2017Alberti et al., , 2018Consolini 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.
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 s 7~1 10 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 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.

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 . 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., 2017Alberti et al., , 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 .
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., 2018Lekscha & 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 magnetosphereionosphere-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.
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.

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 techniquethe EMDand two flavors of RA, we have been able to disentangle the different features of both quiet and disturbed periods. Specifically, we have found that, (i) 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 , are confirmed by both recurrence quantification analysis and recurrence network analysis. (ii) 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 . 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 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 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(Tozzi et al., , 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.