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

1 Introduction

The near-Earth space environment, which comprises the Earth’s ionosphere and magnetosphere, is a highly dynamical system which responds to the changes of the physical conditions of the interplanetary medium and to the solar magnetic activity displaying a very complex dynamics. Since the early 90s, it was realised that in the course of magnetic storms and substorms the Earth’s magnetosphere does not passively respond to the changes of the solar wind conditions but conversely shows a very rich dynamics characterised by nonlinearity (Ahn et al., 1983; Tsurutani et al., 1990; Vassiliadis et al., 1990, and references therein). This nonlinear dynamics appears in a nearly chaotic and critical behaviour, manifesting in scale-invariant features of geomagnetic indices time series (Consolini et al., 1996; Klimas et al., 1996; Consolini, 1997, 2002; Uritsky & Pudovkin, 1998; Uritsky et al., 2002). Currently, the best scenario for the Earth’s magnetospheric dynamics is that of complex, dissipative dynamical system (Consolini, 2018), which is the result of multi-scale nonlinear processes both of purely-internal and externally-driven origin in an out-of-equilibrium configuration.

The overall magnetospheric dynamics, as monitored by geomagnetic indices during magnetic storms and substorms (Daglis et al., 1994; Davies & Sugiura, 1966), is due to both purely-internal and externally-driven processes, which are generally characterised by different spatial and temporal scales. For instance, the Earth’s magnetotail dynamics displays fast intermittent coherent relaxation processes (such as local plasma energization and acceleration occurring in the near- and mid-tail) which manifest in localized auroral emissions in the polar ionosphere during geomagnetic substorms. These coherent intermittent processes (activity bursts) occur on short timescales (fast processes), typically smaller than 100 min, and are due to the loading-unloading dynamics of the Earth’s magnetotail (Kamide & Kokubun, 1996; Consolini & De Michelis, 2005). Conversely, the directly-driven processes show dynamical features occurring on longer timescales (τ > 100 − 200 min, slow processes). In order to avoid any confusion we remark that both the fast and slow processes, taking place during quiet and disturbed geomagnetic conditions and responsible for the observed changes in the geomagnetic indices, are internal to the magnetosphere, although they can be directly-driven and/or triggered by external processes (i.e., changes of the interplanetary conditions).

In a recent paper, Alberti et al. (2017a), using non-traditional methods of analysis based on the empirical mode decomposition (EMD) and delayed mutual information (DMI), have shown that the magnetospheric dynamics is inherently multi-scale during geomagnetic storms, and that fluctuations on short and long timescales are characterized by a different degree of correlation (coupling) according to interplanetary and solar wind parameters changes. In particular, it has been shown that geomagnetic indices fluctuations at scales longer than 200 min share a meaningful information content with a set of external driving functions (e.g., Perrault-Akasofu ε coupling function, vB _s, etc.), while the fluctuations at shorter timescales (below 200 min) do not seem to share a meaningful information with the same interplanetary medium quantities. This result suggests that there is a different degree of correlation between fluctuations of the external drivers and those of the magnetospheric dynamics which depends on the considered timescales. The different degree of correlation between short and long timescale fluctuations and the solar wind changes has been interpreted in terms of a different degree of driving of fast and slow magnetospheric dynamics; the fast magnetospheric dynamics (i.e., that related to short timescale fluctuations) is mainly influenced by the internal magnetospheric state, conversely the slow dynamics, although related to the internal magnetospheric state and processes, is more directly driven by external changes. Thus, the fast and slow dynamics refer to the dynamics of the fluctuations at the different timescales and do not refer to the typical timescales of the magnetospheric response, which can be of the order of some tens of minutes has shown in several previous works (see, e.g., Bargatze et al., 1985; Baker et al., 1995). Furthermore, we remark that the typical timescale of about 200 min is very well in agreement with the timescale found by Tsurutani et al. (1990) where a clear break between the linear and the nonlinear response of the magnetosphere to external solar wind drivers has been found in correspondence of this typical timescale. Here, the terminology fast and slow dynamics refers to the observed timescale separation of the fluctuations which is more or less (linearly or nonlinearly) directly related to the solar wind changes, i.e., more or less affected by externally-driven or purely-internal dynamics.

A correct characterization of the fast and slow dynamics and of the processes at their origin represents a critical issue in the framework of Space Weather studies. Indeed, the understanding of the features of the fluctuations occurring at different timescales is of fundamental importance to correctly forecast the magnetospheric dynamics starting from the measurement of the solar wind conditions. Several attempts to forecast geomagnetic activity in terms of geomagnetic indices have clearly shown the difficulty of reproducing the short timescale dynamics (typically below 1 h) both in terms of fluctuation amplitude and timescales only using the interplanetary magnetic field and solar wind plasma parameters (see, e.g., Pallocchia et al., 2007). The correct modelling of such short timescale certainly represents a central issue to correctly forecast a large series of Space Weather disturbances, such as for instance the geomagnetic induced currents (GICs), which may cause great damages to many anthropic systems.

In this work we present a further step in advancing the understanding and characterization of the fast and slow dynamics involved in the Earth’s magnetospheric response to interplanetary and solar wind changes. We propose a dynamical system approach based on a stochastic Langevin equation to describe both short and long timescale fluctuations observed in the course of geomagnetic storms as monitored by SYM-H index. We attempt an investigation of the state functions governing the dynamics of the fluctuations above and below 200 min and of their evolution/changes between quiet and disturbed periods, using a Fokker-Planck approach to the stochastic Langevin equation. The results clearly demonstrate the different nature of the processes at the basis of the slow and fast dynamics, the former governed by the formation of metastable out-of-equilibrium states, the latter related to a quasi-invariant single-well state function.

The paper is organised as follows: in Section 2 we discuss the data and the methods used in our analysis, providing a detailed introduction to the empirical mode decomposition (EMD) method and to the stochastic Langevin approach along with some introductory examples; Section 3 presents the analysis of SYM-H in the selected time periods along with a detailed discussion of the obtained state functions and their evolution; in Section 4 we traces the conclusions and the relevance of our results in the framework of Space Weather studies.

2 Data and methods

2.1 Data

We focus our attention on three intense geomagnetic storms, characterized by a single minimum in the Disturbance Storm Time (Dst) index, with values lower than −200 nT. Table 1 reports the selected time periods, each of which has a length of 1 month or more centered around the principal storm. The three selected time intervals refer to the famous Bastille’s day storm on 14 July 2000, to the 24 August 2008 and to the St. Patrick’s day on 15 March 2015 storms, the latter being the most prominent storm of the last solar cycle. The selected storms are all associated with fast moving Halo CMEs (on 14 July 2000 at 10:54 UT, 22 August 2005 at 01:31 UT, and 15 March 2015 at 01:48 UT, respectively¹), having a high transit time velocity (>750 km/s), and are related to three forward interplanetary shocks² (observed at the ACE location on 15 July 2000 at 14:15 UT, 24 August 2005 at 05:45 UT, and 17 March 2015 at 04:05 UT, respectively). Moreover, for the Bastille day and the 24 August 2005 storms, the associated solar eruptions produced the emission of solar energetic particles, observed at the Earth as intense events (Laurenza et al., 2009; Alberti et al., 2017b) both in terms of flux intensity (being of class >S2, according to the NOAA scale; Laurenza et al., 2018) and energy content (as ground level enhancements-GLEs were observed by neutron monitors). Finally, these geomagnetic storms were accompanied by Forbush decreases in the cosmic ray intensity, greater than 4% as observed by the SVIRCO Observatory in Rome (http://webusers.fis.uniroma3.it/svirco/pag_1.html). A fourth storm was detected on 2 September 2005, associated with a non halo CME (on 31 August 2005 at 11:30 UT) and possibly a forward interplanetary shock (on 2 September 2005 at 14:50 UT at ACE location).

To correctly monitor the evolution of the three geomagnetic storms and to ensure a good statistics for our stochastic analyses, instead of using the Dst index, which has a time resolution of 1 h, we consider its high resolution version, the low-latitude SYM-H index (Iyemori, 1990). This index, which describes the geomagnetic disturbances at mid-latitudes in terms of symmetric (SYM) disturbances for the horizontal component (H) of the magnetic field near the equator, has a time resolution of 1 min and is considered as a version of the Dst index at higher resolution. Compared to Dst, SYM-H gives the opportunity to better visualize the effects on geomagnetic field due to the variations of the solar wind dynamic pressure and is capable of resolving higher-frequency variations, which are averaged out in the calculation of Dst and that instead can be useful in the studies related to space storm triggers. However, SYM-H differs from Dst not only on the distinct time resolution but also on the number and location of the geomagnetic observatories used and on a different method of convolution applied to produce it (Iyemori, 1990). Nevertheless, it has been found that during the quiet periods and the occurrence of small magnetic storms the values of the two indices differ at most 10 nT and during intense storms the difference is less than 20 nT (Wanliss & Showalter, 2006). Using this index, whose values can be freely retrieved at OMNI website (https://omniweb.gsfc.nasa.gov), each dataset describing the three selected storm periods consists of more than 43 k-points.

Figure 1 shows the SYM-H time series relative to the selected periods. The storm periods can be identified, having SYM-H < −200 nT, along with a series of less relevant disturbances.

Fig. 1. The SYM - H time series of the three periods considered in this work and reported in Table 1.

2.2 The empirical mode decomposition (EMD): A brief introduction

Theoretically, data analysis methods should require minimal assumptions and mathematical artefacts since physical processes are usually not known beforehand. This point is particularly crucial also to avoid misleading results and extract local features from time series. In this framework, a useful tool for analysing both stationary and nonstationary time series data with the above prescriptions is the empirical mode decomposition (EMD), an adaptive data analysis technique developed by Huang et al. (1998).

This algorithmic procedure decomposes a set of observed data, in our case the SYM-H time series Φ_SYM-H(t), into a finite number N of embedded structures having common local features C _k (t), named intrinsic mode functions (or empirical modes), so that$${\mathrm{\Phi }}_{S\mathrm{YM}-\mathrm{H}}\left(t\right)=\sum _{k=1}^N\mathrm{}{C}_k\left(t\right)+\mathrm{res}\left(t\right),$$(1)where res(t) is the final residue of the decomposition, a monotonic function from which no more empirical modes can be extracted (Huang et al., 1998). Each empirical mode has the same number of local extrema (i.e., local maxima and minima) and zero-crossings (or they differ at most by one) and the average of the mean envelope obtained from local maxima and minima envelopes is zero.

The EMD is a completely adaptive and a posteriori decomposition method where the basis functions (i.e., the empirical modes) are derived from the data by an iterative process, namely the sifting process, representing the core of the decomposition procedure. Thus, the number of empirical modes is not “a priori” chosen but it is related to the local features of the signal and its complexity (Huang et al., 1998).

Table 2 summarizes the main steps of the sifting process (see, e.g., Huang et al., 1998; Flandrin et al., 2004; Huang & Wu, 2008; Wu & Huang, 2004; Wu et al., 2008, for more details). By using the so-called Hilbert-Huang Transform, each intrinsic mode function can be seen as an oscillating function with both time-dependent amplitude a _k (t) and frequency ω _k (t), respectively, which are one of the most important novelties introduced by this technique. Indeed, while the sifting process is useful to have no a priori assumptions on the decomposition basis (i.e., carrying out nonlinearities embedded into the time series), the concepts of instantaneous amplitude and frequency are suitable for investigating nonstationary features (Huang et al., 1998). In this way, each mode is characterized by a typical timescale τ _k which can be derived by using different methods based on the auto-correlation analysis, on the local extrema distance and on spectral methods (see, e.g., De Michelis et al., 2012; Alberti et al., 2017a, and references therein). The latter is used in the present work such that the characteristic mean frequency ${\overline{\omega }}_k$ of each empirical mode is estimated by means of the associated Fourier power spectral density S _k (ω) as$${\overline{\omega }}_k=\frac{{\int }_0^{\mathrm{\infty }}\mathrm{}\omega S_k\left(\omega \right)\mathrm{d}\omega }{{\int }_0^{\mathrm{\infty }}\mathrm{}{S}_k\left(\omega \right)\mathrm{d}\omega },$$(2)and, thus, the characteristic timescale is ${\tau }_k=2\pi {\overline{\omega }}_k^{-1}$.

Moreover, since amplitude and frequency of each empirical mode are time-dependent, both local and global time series energy contents can be estimated by the so-called Hilbert-Huang spectrum, a frequency-time distribution of signal amplitude (or energy) permitting the identification of localized features, and by its time integration, the marginal Hilbert spectrum, showing the total energy that each frequency value contribute with (Huang et al., 1998). Finally, since the set of empirical modes, forming the decomposition basis, empirically satisfies mathematical requirements of both completeness and local orthogonality (Huang et al., 1998; Huang & Wu, 2008), the EMD can be used as a filter by reconstructing partial sums of Equation (1) in a chosen frequency range (see, e.g., Alberti et al., 2014, 2017a; Laurenza et al., 2014; De Michelis et al., 2015).

Summarizing, the EMD is particularly suitable to overcome several limitations of other decomposition analysis techniques since it does not require a priori assumptions on the functional form of the basis of the decomposition (as for Fourier or Wavelet analysis), avoiding misleading results and carrying out local (nonstationary and nonlinear) features from time series that cannot usually be obtained by using fixed eigenfunctions. Nevertheless, as for each analysis technique there are some outstanding open problems which need to be outlined out as end effects and/or stopping criteria for the sifting process. Particularly, end effects can produce misleading empirical modes, propagating into the decomposition process through the sifting steps, since they are not the extreme values of the time series, causing divergence problems of the local extrema envelopes. Several methods have been proposed to avoid this problem, mostly based on mirror and/or data extending methods (Huang & Wu, 2008), particularly suitable to get a better spline fit at the ends.

Another frequent drawback of the EMD is the appearance of mode mixing, which is defined as a single empirical mode either consisting of signals of widely and clearly separate scales, or a signal of a similar scale residing in different empirical mode components. This is a direct consequence of signal intermittency, aliasing the time-frequency distribution and devoid empirical modes of physical meaning (Huang et al., 1998, 1999). The first attempt to reduce this problem was the intermittency test proposed by Huang et al. (1999) introducing a sifting process that sets an upper limit of the distance between extrema before each sifting. Although this reduces the mixing of different scale extrema, the EMD ceases to be totally adaptive. A new approach, consisting of a noise-assisted sifting process and known as Ensemble EMD, was proposed by Wu & Huang (2009). A white noise series is added to the original data and the EMD is used to decompose the new data. Then, these two steps are repeated different times by using several white noise series added at each time. Finally, the true empirical mode is obtained as the ensemble means of the corresponding intrinsic mode functions of the decomposition. In this way, a uniform reference frame in the time-frequency space is provided with the added noise collating comparable scales into one empirical mode and separating scales without assumptions, as for the intermittence test. However, this algorithm assumes that the noise content of a time series is reproducible by white noise processes which cannot be really true when natural signals are analyzed (as, e.g., in climate research, Ditlevsen, 1999, but also in magnetospheric one, Silverman & Shapiro, 1983).

The EMD has been recently applied also to study Sun-Earth relationships in the framework of solar and space physics (Stangalini et al., 2014; Vecchio et al., 2012, 2017; Consolini et al., 2017; Piersanti et al., 2017), as well as in the characterization of the magnetospheric and ionospheric dynamics and configuration (Balasis & Egbert, 2006; De Michelis et al., 2012, 2015; De Michelis & Consolini, 2015; Alberti et al., 2016, 2017a; Alberti, 2018) providing evidence for new insights.

2.3 Stochastic Langevin equation and the Fokker-Planck approach

The study of a dynamical system can be approached in terms of a non-linear stochastic model based on a generalized Langevin equation. This approach, that is, for instance, extensively used in the climate framework (see, e.g., Ditlevsen, 1999; Kwasniok & Lohmann, 2009; Livina et al., 2010; Alberti et al., 2014), can be summarised as follows.

Let us consider an N-dimensional dynamical system described by the system variables {X _t } = {X ₁(t), X ₂(t),…, X _N (t)}, and be the time evolution of the ith system variable described by a discrete Langevin dynamics (Langevin, 1908; Lemons & Gythiel, 1997),$$X_i\left(t+1\right)=\mu \left(X_i\left(t\right),t\right)+\sqrt{2D\left(X_i\left(t\right),t\right)}{W}_i\left(t\right)$$(3)where μ(X _i (t), t) and D(X _i (t), t) are the drift and diffusion coefficients, respectively, and W _i (t) represents a Wiener process,$$W_i\left(t\right)=\sqrt{2}\sum _{k=1}^{\mathrm{\infty }}\mathrm{}{\xi }_k\frac{\sin \left[\left(n-\frac{1}{2}\right)\pi t\right]}{\left(n-\frac{1}{2}\right)\pi }$$(4)where ξ _k are independent Gaussian variables having zero mean and unit variance, so that the expectation of W _i (t) is zero (E[W _i (t)] = 0) and its variance is t (Var[W _i (t)] = t) (Wiener, 1923). Then, the statistical properties of the ith system variable can be described by its probability distribution function (PDF) and, in particular, the probability of finding X _i (t) in a state x at a time t, namely p(x, t), is given by the Fokker-Planck equation,$$\frac{\mathrm{\partial }p\left(x,t\right)}{\mathrm{\partial }t}=\frac{\mathrm{\partial }}{\mathrm{\partial }x}\left[\mu \left(x,t\right)p\left(x,t\right)\right]+\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{x}^2}\left[D\left(x,t\right)p\left(x,t\right)\right].$$(5)

This equation provides the time evolution of the probability distribution function p(x, t). From the above equation it is very simple to recover the stationary solution by simple posing ∂ _t p(x, t) = 0.

As an example, we discuss the simple case of a 1-dimensional system, described by a continuous system variable x(t), whose slow dynamics is driven by a forcing term F(x) (i.e., μ(X _i (t), t) = F(x)), while the fast dynamics is described in terms of a noise process η(t) of amplitude σ (i.e., D(X _i (t), t) = σ ²/2, W _i (t) = η(t)), where σ will be the standard deviation of x. In this case the Langevin equation becomes,$$\mathrm{d}x=-U^{\mathrm{\prime }}\left(x\right)\mathrm{d}t+\sigma \mathrm{d}\eta \left(t\right),$$(6)where we assume the existence of the state function U(x) such that $F\left(x\right)=-\frac{\mathrm{\partial }U\left(x\right)}{\mathrm{\partial }x}=-{U\mathrm{\prime }}^{}\left(x\right)$.

From a dynamical system point of view, U(x) captures all the statistical properties of the system, since it can be linked to the probability distribution function of the system variable x. Indeed, in this case the Fokker-Planck equation now reads$$\frac{\mathrm{\partial }p\left(x,t\right)}{\mathrm{\partial }t}=-\frac{\mathrm{\partial }}{\mathrm{\partial }x}\left[U^{\mathrm{\prime }}\left(x\right)p\left(x,t\right)\right]+\frac{{\sigma }^2}{2}\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{x}^2}\left[p\left(x,t\right)\right].$$(7)

In this case, the stationary solution for the probability distribution, p _s(x), is a direct function of the state function U(x), being$$p_{\mathrm{s}}\left(x\right)\sim \exp \left[-\frac{2U\left(x\right)}{{\sigma }^2}\right].$$(8)

Then, once the stationary distribution function p _s(x) is known, by inverting the relation (8) it is possible to get the state function U(x) responsible for the slow dynamics, being U(x) given by the follow expression$$U\left(x\right)=-\frac{{\sigma }^2}{2}\log p_s+c,$$(9)where c is an arbitrary constant, which does not affect the dynamics of x(t).

The stationary solution obtained from Equation (8) is also an equilibrium solution (Ichimaru, 1973), so that the state function U(x) can be used to investigate the number and the nature of the available system states. Indeed, having even-order state functions (i.e., with positive curvature at both minus and plus infinity), the number of zero crossings of the first derivative U′(x) corresponds to the number of system states, while the sign of the second derivative evaluated at the zeros of U′(x) characterises their nature (positive: stable, negative: unstable) (Ditlevsen, 1999; Kwasniok & Lohmann, 2009; Livina et al., 2010; Alberti et al., 2014). When discrete time series are considered, stationary solution refers to ∂p/∂t ~ 0 and it can be seen as a mean of the instantaneous (i.e., at a fixed time t = t _i ) solution over a discrete time interval Δt.

To illustrate the potentiality of such a method we have done a very simple simulation of the stochastic processes described in Equation (6), using a fourth-order double well state function,$$U\left(x\right)=x^4+\alpha x^2+\beta$$(10)and for the stochastic noise η(t) a delta-correlated white noise (〈η(t′)η(t)〉 = δ(t′ − t)). The noise amplitude σ is chosen to ensure transitions between the two wells. Figure 2 shows an example of the simulation results for the dynamics of the stochastic process, x(t), in the case of the fourth-order double well state function. x(t) is characterised by stochastic jumps between the two state function wells (states) plus noisy fluctuations.

Fig. 2. The simulated behaviour of the stochastic process x(t) relative to a motion in a fourth-order double well state function.

Using the Fokker-Planck approach described above, we have computed the stationary probability distribution p _s(t) and the corresponding state function U(x). The results are reported in Figure 3. The obtained state function U(x) is very well in agreement with the characteristic shape of a fourth-order double well.

Fig. 3. The stationary probability distribution p s(t) (left panel) and the state function U(x) (right panel) obtained by applying Equation (9) for the stochstic process reported in Figure 2.

3 Analysis and results

We start our analysis of the three selected periods by decomposing the corresponding time series into empirical modes using the EMD. For each period we find a number of 17 or 18 empirical modes plus a residue. Figure 4 shows the results of EMD method applied to the SYM-H time series of the period B (10 August 2005–10 September 2005), which is characterised by the occurrence of a large geomagnetic storm on 24 August, and two other minor storms on 31 August, and 2 September, respectively. In this specific case we find 17 empirical modes, each of which characterised by a variability on a different timescale, and a residue. As expected, the amplitudes of empirical modes show an increase in correspondence of the three geomagnetic storms. However, the empirical modes associated with the shortest timescales, i.e., C ₁(t) – C ₈(t), show irregular amplitude enhancements also before/after the storms (see, e.g., the time period from 10 to 15 August 2005), suggesting that at these timescales the dynamics of SYM-H may be not directly controlled by the external driving, i.e., not linearly correlated to interplanetary changes. Indeed, according to Alberti et al. (2017a) the fluctuations at these timescales do not show a one-to-one coupling with solar wind and interplanetary parameters, indicating the occurrence of a nonlinear response at these timescales (see, e.g., Tsurutani et al., 1990). Using the Hilbert transform (see, e.g., Huang et al., 1998), we compute the corresponding Hilbert marginal power spectral density (PSD) for each of the three selected periods, which are reported in Figure 5. All the Hilbert marginal PSDs show a quasi-power law behaviour characterised by a spectral slope of approx. −2 over a wide range of frequencies, [2 × 10⁻⁴, 1 × 10⁻¹] min⁻¹ (i.e. for timescales in the range from 10 min to ~1.5 days). Furthermore, there is no evidence for any spectral break at timescales similar to those for which a break is found in AE-indices [3.3 × 10⁻³ min⁻¹;Tsurutani et al., 1990). Our Hilbert marginal PSDs are very well in agreement with previous findings by Wanliss & Showalter (2006), that observed a quasi-Lorentzian PSD with a typical spectral break at about 2 × 10⁻⁶ Hz (1.2 × 10⁻⁴ min⁻¹, ~6 days). As a first step, we apply the stochastic analysis to the raw time series of SYM-H index for the three selected periods in order to compute the average structure of the state function. Figure 6 reports the obtained average state functions U for the three selected periods. The state functions U show very similar shapes characterised by a very deep well, which is associated with the most probable state, and secondary (one or more, less probable) wells.

Fig. 4. 17 empirical modes (C k (t)) and the residue (res(t)) extracted from the SYM-H time series for the time period B, from 10 August 2005 to 10 September 2005. The empirical modes are ordered by increasing characteristic timescale τ k (from top to bottom, see text for the definition).

Fig. 5. Marginal Hilbert power spectral density for the three different selected periods: period A (black dots), period B (red dots), and period C (green dots), respectively. The black dashed line is used as a reference to the spectral slopes.

Fig. 6. The average state function U (in log-scale) as dependent by SYM-H during the three time intervals: A (left panel), B (middle panel), and C (right panel).

From a dynamical system point of view, the most probable state is associated with the minimum value of the state function U. Looking at our results we notice that deep minima in U are all located near values of SYM-H close to 0 nT, suggesting that these minima identify the quiet-time nearly-stationary state, which corresponds to the stable state of the ring current dynamics during quiet-time periods. In the following, we will refer to this ground state in terms of Non-Equilibrium Stationary State (NESS).

Apart from the ground state, the state functions U show a series of local minima for values of SYM-H < 0, which correspond to the emergence of temporary dynamical states during the development of geomagnetic storms. For instance, in the case of the Bastille’s day storm (case study A) a very pronounced minima is observed near SYM-H ~ −300 nT. These secondary minima, which correspond to less probable states, can be identified as the formation of metastable states during the development of a geomagnetic storm. These metastable states have to be associated with the multi-stage processes occurring during the recovery phase of a geomagnetic storm, when in the decay of the ring current several loss mechanisms take place, or with a long-standing solar wind perturbation transferring energy, mass and momentum to the magnetosphere. These metastable states are characterised by a finite lifetime during which the state-describing parameter (i.e., the SYM-H index) reaches and holds for a certain time interval near a fixed value.

The emerging scenario for the geomagnetic storm is, thus, that of an out-of-equilibrium phase transition among ground NESS and time-local metastable states, which happens as a consequence of the changes of the system boundary conditions. In other words, both nodes and saddles occur during a geomagnetic storm with the system moving between the minima of the state functions, from the stable point, classified as a node, to one or more metastable node points, passing through one or more unstable saddle points (local maxima of the state functions). The number and the depth of these metastable states clearly depend on the characteristics of the geomagnetic storm, that means they depend on the behavior of the main and the recovery phases and on the value of the minimum of the SYM-H index. In particular, apart from the ground NESS, associated with the quiet condition, the emergence of secondary minima, i.e., metastable states, depends on both the external drivers (solar wind and interplanetary magnetic field conditions) and the specific internal instantaneous magnetospheric configuration. Indeed, looking at the average state functions reported in Figure 6 we can notice that they present different structures. All are characterised by a deeper minimum corresponding to SYM-H ~ 0 nT but while the state function of period B is characterised by a multi-well structure due to the occurrence of two similar and consecutive geomagnetic storms, those associated with periods A and C are clearly characterized by a single minimum due to the occurrence of a single storm. This clearly suggests that different state functions, i.e., characterized by different multi-well shapes with different depth, can be found for different storm intensities and behaviors.

3.1 Scale-to-scale stochastic analysis

Considering that the Earth’s magnetosphere dynamics during geomagnetic storms and substorms involves both fast and slow dynamical processes on different timescales (Kamide & Kokubun, 1996; Consolini & De Michelis, 2005; Alberti et al., 2017a), we investigate the behavior of the state function U in dependence on the timescale. In order to do this, we reconstruct the state function U for all the timescales identified by the EMD method, i.e. we compute the state function U _k associated with each empirical mode C _k .

Figure 7 reports the obtained empirical mode state functions, U _k , for different characteristic timescales, τ _k . The dynamical configuration changes (or a bifurcation occurs) on a timescale of ~200 min, in agreement with previous findings (Alberti et al., 2017a), at which fast dynamical components can be distinguished from the slow ones. Interestingly, the former are characterised by single-well state functions, suggesting the existence of a single stable state, while the latter seems to present different dynamical states. This feature is typical of dynamical systems where fast and slow dynamics can be ascribed to different sources, as for example the stochastic resonance (Benzi et al., 1981, 1982) or the paleoclimate changes (Ditlevsen, 1999; Alberti et al., 2014). Furthermore, Figure 7 suggests that the slow dynamics of the ring current, which can be directly related to the external driver (i.e., changes in the interplanetary conditions) (Alberti et al., 2017a), is characterised by a multi-state configuration. Indeed, although the most probable state corresponds to SYM-H ~ 0 nT, the enhancement of the ring current density during geomagnetic storms manifests in the appearance of metastable states associated with both positive (during sudden impulse, SI+; Joselyn & Tsurutani, 1990) and negative (during the main phase of the storm) values of the SYM-H index. Conversely, the shape of the state functions characterising the fast dynamics persists in a single-well form which means that processes operating on these timescales are not significantly affected by external sources (Kamide & Kokubun, 1996; Consolini & De Michelis, 2005; Alberti et al., 2017a).

Fig. 7. The state functions U k for different empirical modes with timescales τ k moving from 8 min up to 4 days for the three different periods: A (left panel), B (middle panel), and C (right panel). The emergence of metastable states occurs at τ 200 min.

This scenario is confirmed by the state functions extracted during the time intervals under investigation (see Fig. 8), where fast/slow dynamics is captured by summing empirical modes with characteristic timescales shorter/longer than 200 min, as described in Alberti et al. (2017a). In particular, we can see that the single-well state function associated with the fast dynamics is quasi-invariant. This evidence is substantiated by rescaling the state function U(x) according to the following transformation$$\{\begin{array}{l}x \to \frac{x}{{\sigma }_x}\\ U\left(x\right) \to U\left(x\right)+c\end{array},$$(11)and investigating the collapsing of the state functions. Figure 9 shows the state functions collapsing, which confirms our statement on the fast dynamics.

Fig. 8. State functions of both fast (upper panels) and slow (lower panels) dynamical components during the time periods A (left panels), B (middle panels), and C (right panels).

Fig. 9. Collapsing of the state functions associated with the fast dynamics for the three time periods A (black line), B (red line), and C (green line).

3.2 Time-windowed analysis

As a last step of our analysis, we investigate the evolution of the state functions associated with fast and slow dynamics during the different phases of the geomagnetic storms under investigation. This is done using a moving time-window technique. The size of the time window is 5 days (7200 data points), and state functions are evaluated by using a constant binsize of 10 nT to correctly calculate the statistical quantities and the shift of moving time window is equal to half day (720 data points).

We start our analysis by investigating the fast dynamics features. Figure 10 shows the state functions related to the fast component (characteristic timescales τ < 200 min) of the SYM-H index evaluated over the different sliding windows.

Fig. 10. Time-windowed contour of the state functions of the fast dynamics of SYM - H during the different storm phases for the selected time periods: A (left panel), B (middle panel), and C (right panel). Here, QT = quiet time, SSC = storm sudden commencement, ST = storm time, RP = recovery phase.

The SYM-H fast dynamics shows a single-state configuration (single-well state function) during the selected periods, even during storm-time intervals, suggesting that this component, and consequently the associated physical processes, is not substantially affected by external-driving sources (Kamide & Kokubun, 1996; Consolini & De Michelis, 2005; Alberti et al., 2017a). It is evident in the case of the Bastille’s and August 2005 storms. The only observable changes are related with the additive constant that is a measure of the relative contribution in probability of the fluctuations associated with these short timescales. In other words, during the storm time the state function values tend to increase suggesting that, although fast dynamics is not directly driven by the external drivers, it triggers external variations (Alberti et al., 2017a) capable of moving the state function towards a less probable (metastable) state. To avoid any confusion, we again remark that with the term triggering we mean the action of a perturbation which does not imply a linear response, conversely the term directly driven refers to a (quasi-)linear response to a perturbation.

The situation is completely different in the case of the slow dynamics state functions which show a clear evolution in the course of the geomagnetic storm phases. In Figure 11, we report the evolution of the state functions related to the slow dynamics component of SYM-H index. The state function shows a significant evolution with the different storm phases, displaying the formation of short-living metastable states which form during the main and the recovery phase, respectively. This is the effect of the competition between the external driving and the internal ring current dynamics. To stress this point we report in Figure 12 the evolution of the state function during the different storm phases for the Bastille event. This plot clearly shows how during the storm time the state function is characterized by many local minima (multi-well state function), which can be interpreted as short-living metastable states. The multi-well nature of the storm time state function also suggests that different internal processes, characterized by different decay times, can be involved in the evolution of the ring current during the main and recovery phases of the storm. This has been widely documented in literature, where different mechanisms have been claimed to work for ring current decay such as wave-particle scattering, O⁺/H⁺ charge-exchange, Coulomb scattering and ion-convection out-flows at the magnetopause frontside. All these mechanisms act during the recovery phase of a storm and are characterized by different timescales (see, e.g., Hamilton et al., 1988; Kozyra et al., 2002;Weygand & McPherron, 2006, and references therein). This supports the common view according to which the magnetospheric dynamics is governed by internal processes, although the slow dynamics is directly driven by (i.e., correlated to) the external driver (Alberti et al., 2017a).

Fig. 11. Time-windowed contour of the state functions of the slow dynamics of SYM-H during the different storm phases for the selected time periods: A (left panel), B (middle panel), and C (right panel). As for Figure 10, QT = quiet time, SSC = storm sudden commencement, ST = storm time, RP = recovery phase.

Fig. 12. The evolution of the state function (bottom panels) associated with the SYM-H slow dynamics (upper panels) for the different phases of the Bastille’s day storm during the period from 01 July to 25 July 2000.

Another interesting feature is that the quiet time dynamics of both fast and slow dynamics is characterized by a single well state function U(x) (see Fig. 13), suggesting that during quiet conditions the ring current, as monitored by SYM-H index, is in a quasi-stable state of minimum energy configuration. Consequently, the timescale separation emerges only during disturbed periods, suggesting that the external disturbances directly drive only a certain range of timescales.

Fig. 13. Comparison of scaled state functions for the two different ranges of timescales during the interval of quiet time conditions (01–10 July 2000).

4 Discussion and conclusions

Before discussing the possible implications of our results, we start summarizing the main features that emerge from our analysis:

1. the analysis, resulting from a combination of the EMD method and the Langevin/Fokker-Planck approach, confirmed the previous findings by Alberti et al. (2017a) on the existence of two characteristic ranges of timescales, which have been identified as fast and slow dynamics;

2. these two ranges of timescales emerge during disturbed periods, i.e., during geomagnetic storms, when the state function associated with slow dynamics undergoes to significant changes.

This suggests that the two ranges of timescales are characterized by a different dynamical behavior during geomagnetic storms: the fast dynamics is practically unperturbed by the changes in the solar wind conditions (the state function remains the same without showing the emergence of secondary or metastable state), the slow dynamics is conversely strongly affected by external conditions showing dynamical changes in the state function with the formation of short-living metastable states.

The above points move toward the existence of a very well defined separation of timescales in the response of the Earth’s magnetosphere to the solar wind changes, which clearly emerges during perturbed periods. This separation can be due to different physical processes which are/are not affected by external dynamics during geomagnetic storms.

Indeed, during quiet periods, i.e., when there is no plasma entry from the interplanetary medium due to long-standing periods of northward interplanetary magnetic field, SYM-H is characterized by a single quasi-stable state of minimum energy configuration. Conversely, during geomagnetic storms a clear separation of timescales emerges in SYM-H.

In terms of a stochastic Langevin equation the observed behavior of the fast and slow dynamics can be described by two different equations. The fast dynamics is essentially characterized by a quasi-stable single-well state function which is nearly quadratic, so that the associated Langevin equation takes the form,$$\mathrm{d}x\left(t\right)=-\gamma x\left(t\right)\mathrm{d}t+\sigma \mathrm{d}\eta \left(t\right),$$(12)where γx(t) = U′(x), σ is the noise term amplitude and η(t) is the stochastic noise process. The slow dynamics shows a state function which undergoes to dynamical changes during disturbed periods and under the effect of external driving, so that a good modeling of slow dynamics in terms of Langevin equation is the following,$$\mathrm{d}x\left(t\right)=U^{\mathrm{\prime }}\left(x\left\{{\alpha }_i\left(t\right)\right\}, t\right)\mathrm{d}t+\sigma \mathrm{d}\eta \left(t\right),$$(13)where U(x|{α _i (t)}, t) is the state function which depends on both time $t$ and a set of parameters {α _i (t)}, reasonably representing the changes of the interplanetary magnetic field and solar wind conditions (Alberti et al., 2017a). In other words, the topology of the state function U(x|{α _i (t)}, t) contains valleys and hills depending on the set of parameters {α _i (t)}.

The fact that a correct description of the slow dynamics is represented by a state function U(x|{α _i (t)}, t), which undergoes to dynamical changes as a function of the external driving, is the evidence for the occurrence of a first-order-like topological phase transition (Chang et al., 1992, 2003). In this scenario, the set of parameters {α _i (t)} can be associated with a sort of critical parameter vector $\vec{c}\left(t\right)$, which governs the occurrence of dynamical transitions. Here, this dynamical transition takes the form of a order-disorder transition (Chang et al., 2003).

The above scenario for the slow dynamics of SYM-H confirms the previous results by Balasis et al. (2006) and Wanliss & Dobias (2007) on the occurrence of a dynamical phase transition in proximity of and during a geomagnetic storm (see also Wanliss, 2005). Similar concepts have also proposed to explain the dynamics of geomagnetic substorms (see, e.g., Sitnov et al., 2000, 2001).

Regarding the fast dynamics the quasi-independent shape of the state function, which acquires a near quadratic shape, suggests that this range of fluctuation timescales (τ < 200 min) may be directly connected to processes of internal origin, which are less affected by the external changes of interplanetary magnetic field and solar wind conditions (Alberti et al., 2017a). A possible explanation of this short-timescale fluctuations could be related to spatio-temporal turbulence occurring in the equatorial plasma sheet regions. A more detailed investigation of the features of these short timescale fluctuations is required in order to better clarify the meaning of the two terms of the corresponding stochastic Langevin equation. This point is demanded to a further work, strictly devoted to the features of the fast dynamics.

In conclusions, we have presented an extended analysis of the fluctuations occurring during geomagnetic storms as monitored by SYM-H index in terms of a Langevin/Fokker-Planck approach, showing how the features at timescales longer/shorter than 200 min are completely different. These results are important and crucial to a correct forecasting of geomagnetic indices, such as SYM-H, during magnetic storms, especially in relationship with the forecasting of the internal fast dynamics which is fundamental in the framework of Space Weather.

Acknowledgments

This research has been supported by the Italian MIUR-PRIN grant 2012P2HRCR on “The Active Sun and Its effects on Space and Earth Climate”. The OMNI data were obtained from the GSFC/SPDF OMNIWeb interface at https://omniweb.gsfc.nasa.gov/. The authors also acknowledge World Data Center for Geomagnetism (Kyoto) for the use of the geomagnetic indices data. All the derived data products presented in this paper are available upon request by email to the Author (tommaso.alberti@inaf.it). The editor thanks two anonymous referees for their assistance in evaluating this paper.

References

All Tables

All Figures

	Fig. 1. The SYM - H time series of the three periods considered in this work and reported in Table 1.
In the text

	Fig. 2. The simulated behaviour of the stochastic process x(t) relative to a motion in a fourth-order double well state function.
In the text

	Fig. 3. The stationary probability distribution p s(t) (left panel) and the state function U(x) (right panel) obtained by applying Equation (9) for the stochstic process reported in Figure 2.
In the text

	Fig. 4. 17 empirical modes (C k (t)) and the residue (res(t)) extracted from the SYM-H time series for the time period B, from 10 August 2005 to 10 September 2005. The empirical modes are ordered by increasing characteristic timescale τ k (from top to bottom, see text for the definition).
In the text

	Fig. 5. Marginal Hilbert power spectral density for the three different selected periods: period A (black dots), period B (red dots), and period C (green dots), respectively. The black dashed line is used as a reference to the spectral slopes.
In the text

	Fig. 6. The average state function U (in log-scale) as dependent by SYM-H during the three time intervals: A (left panel), B (middle panel), and C (right panel).
In the text

	Fig. 7. The state functions U k for different empirical modes with timescales τ k moving from 8 min up to 4 days for the three different periods: A (left panel), B (middle panel), and C (right panel). The emergence of metastable states occurs at τ 200 min.
In the text

	Fig. 8. State functions of both fast (upper panels) and slow (lower panels) dynamical components during the time periods A (left panels), B (middle panels), and C (right panels).
In the text

	Fig. 9. Collapsing of the state functions associated with the fast dynamics for the three time periods A (black line), B (red line), and C (green line).
In the text

	Fig. 10. Time-windowed contour of the state functions of the fast dynamics of SYM - H during the different storm phases for the selected time periods: A (left panel), B (middle panel), and C (right panel). Here, QT = quiet time, SSC = storm sudden commencement, ST = storm time, RP = recovery phase.
In the text

	Fig. 11. Time-windowed contour of the state functions of the slow dynamics of SYM-H during the different storm phases for the selected time periods: A (left panel), B (middle panel), and C (right panel). As for Figure 10, QT = quiet time, SSC = storm sudden commencement, ST = storm time, RP = recovery phase.
In the text

	Fig. 12. The evolution of the state function (bottom panels) associated with the SYM-H slow dynamics (upper panels) for the different phases of the Bastille’s day storm during the period from 01 July to 25 July 2000.
In the text

	Fig. 13. Comparison of scaled state functions for the two different ranges of timescales during the interval of quiet time conditions (01–10 July 2000).
In the text