Multifractal and Chaotic Properties of Solar Wind at MHD and Kinetic Domains: An Empirical Mode Decomposition Approach

Abstract

Turbulence, intermittency, and self-organized structures in space plasmas can be investigated by using a multifractal formalism mostly based on the canonical structure function analysis with fixed constraints about stationarity, linearity, and scales. Here, the Empirical Mode Decomposition (EMD) method is firstly used to investigate timescale fluctuations of the solar wind magnetic field components; then, by exploiting the local properties of fluctuations, the structure function analysis is used to gain insights into the scaling properties of both inertial and kinetic/dissipative ranges. Results show that while the inertial range dynamics can be described in a multifractal framework, characterizing an unstable fixed point of the system, the kinetic/dissipative range dynamics is well described by using a monofractal approach, because it is a stable fixed point of the system, unless it has a higher degree of complexity and chaos.

1. Introduction

The interplanetary space is permeated by a supersonic and super-Alfvénic plasma known as solar wind which develops a strong turbulent character during its expansion phase [1]. Due to the presence of a mean magnetic field, solar wind low-frequency fluctuations are usually described within the magnetohydrodynamic (MHD) framework [2,3,4]. These fluctuations show turbulent properties that are characterized by a quasi-Kolmogorov energy scaling [5,6,7,8]. Indeed, magnetic energy density seems to follow a spectral decay as $E\left(k\right)\sim k^{-5/3}$, although the theoretical scaling derived from MHD equations suggests a slightly different spectral exponent, e.g., $E\left(k\right)\sim k^{-3/2}$ [9,10] for Alfvènic turbulence as it should be in the case of solar wind. Turbulence is a phenomenon showing the presence of small scale fluctuations in the velocity and pressure fields (for fluids), as well as in the magnetic field (for plasmas), and an increased rate of mixing of mass and momentum [1,11]. Turbulent flows exhibit characteristic phenomena like coherent structures in the flow and intermittency. Coherent structures are usually defined as regions of concentrated vorticity where phase correlation exists with a typical lifetime larger than that of the stochastic fluctuations surrounding them, while intermittency is the manifestation of sudden field changes, modifying the shape of the probability distribution functions of field gradients (e.g., velocity and temperature in fluids, magnetic in plasmas) [12,13], and resulting in an anomalous scaling of the field increments [14]. As the analytic and numerical solution of such flows is expensive, investigators rely on models to simulate and simplify their dynamics. Such turbulence models include two-equation models (like the k-$ϵ$ model and the k-$\omega$ one [15]), Reynolds stress models (like the Speziale-Sarkar-Gatski model [16] and the Mishra-Girimaji model [17]), along with models in Large Eddy Simulations [18]. Similar models have been also developed for describing turbulent features in plasmas like two-dimensional hybrid-Vlasov simulations [19], compressible Hall MHD direct numerical simulations [20], and shell models [21].

In the framework of turbulence, several phenomena inside the MHD/inertial domain are described, by using the nonlinear energy cascade via the Yaglom law, an exact relation for the scaling of the third-order moment of fluctuations [22,23], or by analyzing the role of the intermittency in changing scaling properties of magnetic fluctuations within a multifractal approach [14,24,25,26]. Both previous findings are derived from the structure function analysis, through which scaling laws can be investigated (e.g., [27]), exploiting Kolmogorov’s universality assumptions (e.g., [28]). More specifically, a turbulent flow is sustained by a persistent source of energy which is rapidly dissipated via the so-called nonlinear energy cascade [29], converting the kinetic energy into internal energy through viscous processes. Indeed, turbulence causes the formation of eddies at different scales and energy is transferred from large- to small-scale structures through an inertial and inviscid mechanism [5,28], i.e., the nonlinear energy cascade. According to Kolmogorov’s theory, if all possible symmetries of the Navier-Stokes equation are restored in a statistical sense, then the turbulent flow is self-similar at small scales and has a finite mean energy dissipation rate $ϵ$ such that, at very high (but not infinite) Reynolds numbers, all small-scale statistical properties are uniquely and universally determined by the length scale ℓ, the mean energy dissipation rate $ϵ$, and the kinematic viscosity $\nu$ (e.g., [28]). By simply exploiting dimensional arguments, these assumptions imply that the energy spectrum at large wavenumbers assumes a universal form as

$$E\left(k\right)=F\left(\nu \right){ϵ}^{2/3}{k}^{-5/3}.$$

(1)

In the limit of infinite Reynolds numbers, Equation (1) becomes independent by the viscosity $\nu$ such that

$$E\left(k\right)=C{ϵ}^{2/3}{k}^{-5/3}.$$

(2)

being C a universal dimensionless constant [28]. The above assumptions are only valid for all those scales which are smaller than the integral scale L, where long-range correlations between particles are found, and are greater than the dissipative scale ${\ell }_D$, where viscosity dominates and the turbulent kinetic energy is dissipated, i.e., ${\ell }_D\ll \ell \ll L$ (e.g., [5,28]). Particularly, the dissipation of kinetic energy mostly takes place at the so-called Kolmogorov microscale defined as

$$\eta ={\left(\frac{{\nu }^3}{ϵ}\right)}^{1/4}$$

(3)

which is well separated from the integral scale L, corresponding to the size of the eddies when they are formed. These two scales mark the extrema of the energy cascade: since eddies with size L are much larger than the dissipative eddies with size $\eta$, kinetic energy is not dissipated at large scales but it is essentially transferred to smaller scales where viscous effects become dominant. Within this range, where nonlinear interactions between eddies take place, inertial effects are larger than viscous ones such that it is usually named “inertial range”. Due to the large separation between L and $\eta$, the dissipation rate is primarily determined by the large scales since viscous effects at Kolmogorov scales rapidly dissipate energy. Thus, the overall rate of dissipation is only controlled by the nonlinear scale-to-scale transfer such that the dissipation rate is approximately given as

$$ϵ=\frac{u^3}{L}$$

(4)

being u the bulk velocity of the flow, and consequently

$$\frac{\eta }{L}=R{e}^{-3/4}$$

(5)

where $Re$ is the Reynolds number. Equation (5) can be used as a measure of the number of scales within the inertial range (i.e., the extension of this range of scales) which only depends by the Reynolds number of the flow.

Similar assumptions and scalings can be also derived by using MHD equations for describing plasma dynamics and, particularly, within the same Kolmogorov’s assumptions of isotropy, homogeneity, and stationarity, and exploiting dimensional arguments, the energy spectrum at large wavenumbers behaves as

$$E\left(k\right)=C{\left(ϵv_0\right)}^{1/2}{k}^{-3/2}.$$

(6)

where C a universal dimensionless constant and $v_0$ is the rms component of the total turbulent velocity [9,10].

Kolmogorov assumptions break down beyond a scale ${\ell }_b\sim {\ell }_i$, being ${\ell }_i$ the ion inertial length, where the MHD/inertial description cannot be used [30,31], and for which a steeper slope of the energy density spectrum is found [32]. The second power-law domain at small scales can be explained by invoking several possible dispersive/kinetic phenomena as wave-wave coupling, Landau damping, Kinetic Alfvén Waves (KAWs), and so on [33,34,35].

Solar wind magnetic field fluctuations, including spectral features like the Kolmogorov energy scaling, are usually investigated by following a statistical approach of magnetic increments

$$\delta \mathbf{b}(\mathbf{x},\ell )\doteq \mathbf{b}(\mathbf{x},\ell )-\mathbf{b}\left(\mathbf{x}\right)$$

(7)

requiring that the statistics of these increments is invariant under arbitrary translations

$$\langle \delta \mathbf{b}(\mathbf{x}+\mathbf{r},\ell )\rangle =\langle \delta \mathbf{b}(\mathbf{x},\ell )\rangle .$$

(8)

In this way, it is possible to derive an exact law for the scaling of the mean-square magnetic field increments between two points separated by a distance ℓ according to which $\langle {\left(\delta \mathbf{b}(\ell )\right)}^2\rangle \sim {\ell }^{1/2}$, from which the Iroshnikov-Kraichnan $k^{-3/2}$ scaling law can be simply derived by noting that $\langle {\left(\delta \mathbf{b}(\ell )\right)}^2\rangle$ is related to the magnetic energy [9,10]. Moreover, from dimensional analysis it is easy to recover the second-order structure function which is equal to $S_2\left(l\right)=\langle {\left(\delta \mathbf{b}(\ell )\right)}^2\rangle =C_{IK}{\left(ϵc_A\right)}^{1/2}{\ell }^{1/2}$, where $c_A$ is the Alfvén speed and $C_{IK}$ is a universal constant. More generally, generalized structure functions can be calculated from any arbitrary and finite order q from increments such that

$$S_q(\ell )\doteq \langle \mid \delta \mathbf{b}(\mathbf{x},\ell ){\mid }^q\rangle .$$

(9)

They are extensively used to investigate scaling properties of both velocity and magnetic field fluctuations (e.g., [27]), mostly devoted to the characterization of intermittency and self-similarity properties of solar wind turbulence [14,36,37,38]. In this context, different multifractal models (e.g., $\beta$-model, p-model, and their variations (e.g., [24,26])) have been proposed to explain the evolution of intermittency across the heliosphere, with the solar wind becoming more multifractal in nature when leaving the Sun [1]. In particular, the analysis of scaling exponents $\zeta \left(q\right)$ of the $q$-th order structure function shows that different turbulent scenarios can develop, being mainly characterized by a nonlinear dependence of $\zeta \left(q\right)$ with q (e.g., [39,40,41,42]).

This paper approaches the study of the self-similarity properties of solar wind magnetic field fluctuations at different timescales by using a novel method to evaluate structure functions at different orders q. It is based on the Empirical Mode Decomposition (EMD) which allows to correctly derive the timescales embedded into the analyzed time series, as well as to better evaluate structure functions by using local (non-constant) timescales [43,44]. The results evidence a different behavior of magnetic field fluctuations at MHD/inertial and kinetic/dissipative scales. While the former are characterized by a multifractal character, the latter show a monofractal scaling. In a dynamical system framework, both behaviors can be seen as corresponding to two different fixed points: an unstable saddle for the MHD/inertial domain, and a stable node for the dissipative one. The solar wind magnetic field fluctuations undergo a saddle-node bifurcation when moving from the MHD/inertial down to kinetic/dissipative scales.

2. Data

We consider solar wind magnetic field measurements from ESA-Cluster mission on 10 January 2004 from 05:30 UT to 06:30 UT. This period is characterized by a high speed (i.e., it is a fast stream, $v\sim$ 540 km/s), a mean magnetic field intensity ($B\sim$ 11 nT), and a high density ($n\sim$ 14 cm${}^{-3}$). We used combined magnetic field data from the fluxgate magnetometer (FGM) and the experiment called “spatiotemporal analysis of field fluctuations” (STAFF) onboard Cluster 3 spacecraft to obtain a resolution of data equals to 450 Hz. For computational purposes, the time resolution has been reduced of a factor 4, moving it to $\Delta t$ = 8.9 ms.

Figure 1 shows the magnetic field intensity (B) and the three magnetic field components ($B_x,B_y,B_z$) for the selected time interval in the GSE reference system.

3. Methods

3.1. The Empirical Mode Decomposition (EMD): A Brief History

During the past decades, several decomposition procedures have been suggested to investigate scale variability of time series. Most common methods rely on Fourier-based techniques, Wavelet transforms and/or eigenfunction analysis (e.g., [45]). These methods, by choosing a decomposition basis in a mathematical space with requirements of completeness and orthogonality, allow to obtain oscillating components, with fixed scales and amplitudes, embedded inside time series (e.g., [45]). However, neither stationarity nor linearity is assured when natural phenomena are investigated, unless scaling law theory is mostly derived by exploiting these two requirements (e.g., [27,28]). Recently, the Empirical Mode Decomposition (EMD) has been developed to provide a suitable decomposition method for time series by exploiting their local properties, allowing us to reduce mathematical assumptions by using a completely adaptive and a posteriori decomposition procedure where the number of the extracted empirical modes depends on the signal complexity [43]. The EMD carries out a finite set of embedded modes, usually named Intrinsic Mode Functions (IMFs), from a given time series $x\left(t\right)$ by using an iterative process known as sifting process. The main steps of this process can be summarized as follows:

• evaluate the mean of a signal $x\left(t\right)$ and subtract from it to produce a zero-mean signal $x_m\left(t\right)=x\left(t\right)-\langle x\left(t\right)\rangle$;
• find local maxima and minima of $x_m\left(t\right)$;
• use a cubic spline to evaluate the upper ($e_{max}\left(t\right)$) and lower ($e_{min}\left(t\right)$) envelopes from local maxima and minima, respectively;
• evaluate the mean envelope $e_m\left(t\right)$ and subtract from $x_m\left(t\right)$ to have $h\left(t\right)=x_m\left(t\right)-e_m\left(t\right)$;
• check if $h\left(t\right)$, often called detail or “candidate” IMF, is an IMF that is, check if the number of zero crossings and local extrema differs at most by one and if the local mean is zero;
• if $h\left(t\right)$ is an IMF, then store it ($c_k\left(t\right)=h\left(t\right)$), else repeat steps from 1 to 5 on the signal $x_h\left(t\right)=x_m\left(t\right)-h\left(t\right)$ until an IMF is obtained.

Once the decomposition is complete, i.e., when no more IMFs can be extracted from $x\left(t\right)$, the time series $x\left(t\right)$ can be written as

$$x\left(t\right)=\sum _{k=1}^N{c}_k\left(t\right)+r\left(t\right)$$

(10)

where $r\left(t\right)$ is the residue of the decomposition, a non-oscillating function of time [43]. Mathematically, the sifting process stops only when the number of iterations $n\to \infty$; numerically, it can be stopped after $n^{\prime }$ iterations by defining a stopping criterion [46] like the Cauchy convergence criterion [43], according to which the sifting algorithm stops when ${\sigma }_{n^{\prime }}<{\sigma }_0$, being ${\sigma }_{n^{\prime }}={\sum }_{j=1}^T\frac{|h_{n^{\prime }}\left(t_j\right)-h_{n^{\prime }-1}\left(t_j\right){|}^2}{h_{n^{\prime }-1}^2\left(t_j\right)}$, where $h_{n^{\prime }}$ is the $n^{\prime }$ detail and T the length of the time series $x\left(t\right)$, and ${\sigma }_0$ is a threshold value which usually varies between 0.2 and 0.3 [43], or by the threshold method proposed in Reference [47] in which two thresholds, ${\theta }_1$ and ${\theta }_2$, are chosen to guarantee globally small fluctuations and, in the meanwhile, to take into account locally large excursions. In this way, by defining $\sigma \left(t\right)=\left|\frac{2h_{n^{\prime }}\left(t\right)}{e_{max}\left(t\right)-e_{min}\left(t\right)}\right|$, the sifting process is iterated until $\sigma \left(t\right)<{\theta }_1$ for a prescribed fraction $1-\alpha$ of the total duration, and $\sigma \left(t\right)<{\theta }_2$ for the remaining fraction, being typically ${\theta }_1=0.05$ and ${\theta }_2=10{\theta }_1$ [47,48]. More details about the sifting process and its features can be found in several previous works (e.g., [43,47,48,49]).

The EMD is a fundamental step for providing non-stationary oscillating components which can be used as inputs for the Hilbert Spectral Analysis (HSA), which permits us to investigate amplitude-frequency modulation embedded in time series (e.g., [43,50]). Through the Hilbert Transform (HT), which is a linear mathematical operator that takes each IMF $c_k\left(t\right)$ and produces a function $H\left[c_k\right]\left(t\right)$ by convolution with the function $\frac{1}{\pi t}$, each empirical mode can be written as modulated both in amplitude and in frequency

$$c_k\left(t\right)=a_k\left(t\right)\Re \left\{exp\left[i2\pi {\int }_0^t{f}_k\left(t^{\prime }\right)d{t}^{\prime }\right]\right\}$$

(11)

where $a_k\left(t\right)$ and $f_k\left(t\right)$ are the instantaneous amplitude and frequency of the k-th empirical mode, respectively, and ℜ is the real part of the exponential. The HT allows to investigate non-stationary features of the time series, being $f_k\left(t\right)$ a function of time, and also its nonlinear behavior, due to the time-dependence of $a_k\left(t\right)$ (e.g., [43,51]). Derived from both $a_k\left(t\right)$ and $f_k\left(t\right)$, the instantaneous local energy content $E(t,f)$ is studied by contouring the squared-amplitude in a time-frequency plane, i.e., by defining the so-called Hilbert-Huang spectrum [43]. Then, an intermittency measure, similar to that defined by using wavelet analysis, can be introduced as

$$DS\left(f\right)=\frac{1}{n\Delta t}{\int }_t{\left[1-\frac{H(t^{\prime },f)}{h\left(f\right)}\right]}^{2}d{t}^{\prime }$$

(12)

where $h\left(f\right)={\langle H(t^{\prime },f)\rangle }_t$ and $n\Delta t$ is the length of the time series (e.g., [43]). It is often called Degree of Stationarity (DS) (e.g., [43]) and a time series is statistically stationary if DS = 1.

Figure 2 reports the degree of stationarity for the three magnetic field components. A clear increase in the stationary character of time series is found when approaching the frequency $f_b=$ 0.4 Hz $\sim f_i$, being $f_i$ the Doppler-shifted ion cyclotron frequency. This suggests that a high non-stationary behavior characterizes the inertial regime, where MHD processes govern the dynamics of the system, while dissipative processes are characterized by a nearly-stationary dynamics as also previously observed (e.g., [52]). The non-stationary character observed in the MHD/inertial domain could be a counterpart of the intermittent nature of fluctuations in the inertial range.

3.2. The EMD-Based Multifractal Analysis

Recently, a method capable of detecting the fractal dimension of a time series by partitioning the time and scale domain of a signal into fractal dimension regions has been proposed. This method, which is similar to the Wavelet Transform Modulus Maxima (WTMM), is an EMD-based multifractal analysis. It is named EMD-based dominant amplitude multifractal formalism (DAMF) [44] and it has been proposed to investigate singularities and (multi)fractal behavior of time series. The EMD-DAMF method can be summarized in the following steps:

• derive instantaneous amplitude $a_k\left(t\right)$ and mean timescale ${\tau }_k={\langle f_k\left(t\right)\rangle }_t^{-1}$ of each empirical mode;
• determine the dominant amplitude coefficients $u_{j,k}$ over a time support $I_{j,k}$ around the j-th local maximum

  $$u_{j,k}\doteq \underset{k^{\prime }\le k}{sup}\left\{max\left\{|a_{k^{\prime }}(t\in I_{j,k})|\right\}\right\}$$

  (13)

  with $j=1,\cdots ,N_k$, being $N_k$ the number of local maxima of $a_k\left(t\right)$, and $k=1,\cdots ,N$;
• evaluate the q-th-order structure function $S_q\left({\tau }_k\right)$

  $$S_q\left({\tau }_k\right)=\frac{1}{N_k}\sum _{j=1}^{N_k}{\left\{u_{j,k}\right\}}^q;$$

  (14)
• estimate the scaling exponent $\zeta \left(q\right)$ as the linear slope, in a log-log space, of $S_q\left({\tau }_k\right)$ vs. ${\tau }_k$, such that

  $$S_q\left({\tau }_k\right)\sim {\tau }_k^{\zeta \left(q\right)};$$

  (15)
• derive the singularity strengths $\alpha$ and spectrum $f\left(\alpha \right)$ by using the Legendre transform of the scaling exponents $\zeta \left(q\right)$ as usual

  $$\alpha =\frac{d\zeta \left(q\right)}{dq}\&f\left(\alpha \right)=\alpha q-\zeta \left(q\right).$$

  (16)

The main novelty introduced by this method is that structure functions $S_q\left({\tau }_k\right)$ are derived by exploiting the local features of empirical modes such that local extrema can be used to correctly calculate differences/increments between two points, instead of considering a fixed timescale as for canonical structure function analysis. Moreover, timescales are not fixed a priori but they are derived from the EMD analysis of time series considering a finite set of multiresolution coefficients $u_{j,k}$. This allows us to have a limited (and small) number of points in the scaling range such that the scaling exponents can be better evaluated and visually inspected.

4. Results from the EMD-Based Multifractal Analysis

Figure 3 and Figure 4 report the EMD-DAMF results at MHD/inertial and kinetic/dissipative scales for each magnetic field component, respectively. In each figure, the second-order structure function $S_2\left(\tau \right)$ is shown in the upper panel (multiplied by ${\tau }^{1/2}$ and ${\tau }^{3/2}$ to have a compensated structure function), the scaling exponents $\zeta \left(q\right)$ are reported in the middle panel, and the singularity spectrum $f\left(\alpha \right)$ is displayed in the lower panel.

The EMD-DAMF analysis at MHD/inertial scales (Figure 3), i.e., corresponding to the inertial range, can be carried out by considering empirical modes with mean timescales in the range 2–500 s (or $f\in ({10}^{-3},0.4)$ Hz). Indeed, as shown by the second-order structure function $S_2\left(\tau \right)$ a scale-break is found when $f=f_b=$ 0.4 Hz $\sim f_i$. As expected from structure function theory of the MHD/inertial domain (e.g., [9,10]), the second-order structure function behaves as ${\tau }^{1/2}$, suggesting that the Fourier energy spectral density decays as $f^{-3/2}$ (or $k^{-3/2}$ assuming Taylor’s hypothesis) (e.g., [8,9,10]).

This result supports the common view according to which energy is injected at large scales (i.e., larger than a typical injection scale L) and transferred to small scales (i.e., smaller than a dissipative scale ${\ell }_D$) through nonlinear interactions and phenomena taking place at scales ℓ, being ${\ell }_D\ll \ell \ll L$ (e.g., [3,4,24,26]). Here, ${\ell }_D$ stands for the dissipation scale (equivalent to Kolmogorov’s scale in fluid turbulence). This result has been obtained by considering the “true” timescales which are embedded in the raw time series and extracted via an adaptive procedure, with no assumptions on the stationarity of oscillating components.

By considering structure functions $S_q\left(\tau \right)$ with 2 s $<\tau <$ 500 s, the scaling exponents are derived and shown in the middle panel of Figure 3 for the three magnetic field components. From a theoretical point of view, assuming homogeneity, isotropy and scale-invariance of the time series we should obtain $\zeta \left(q\right)=q/4$ in the case of Alfvènic MHD turbulence [10]. Our results show a different behavior with scaling exponents characterized by a nonlinear convex trend with the moment order q like $\zeta \left(q\right)\sim q/4+\phi (q/4)$ [26]. This deviation is the fingerprint of the occurrence of anomalous scaling features, i.e., of an intermittency phenomenon, in the nonlinear energy cascade of the magnetic field, suggesting nonlinear two-point correlations in the real space [39,40,41]. Interestingly, different scaling exponents are obtained for the different magnetic field components indicating the existence of an anisotropy of the scaling features in the different directions which may reflect the anisotropic nature of the fluctuation field [7]. Moreover, the observed nonlinear scaling of $\zeta \left(q\right)$ suggests that the probability distribution functions (PDFs) of increments at MHD/inertial scales are characterized by multifractal scaling features. This aspect can be clearly seen by looking at the singularity spectrum $f\left(\alpha \right)$ (see Figure 3, lower panel) which shows a wide range of singularities (0.1 $<\alpha <$ 0.8) for all components. Wider singularities are found for the $B_y$ component, while a narrower spectrum is found for $B_z$.

A clear different behavior is found when approaching the dissipative range (see Figure 4), i.e., moving towards higher frequencies ($f>f_b=0.4$ Hz). A different scaling law is recovered, moving to a greater scaling exponent (${\tau }^{3/2})$ and, consequently, a steeper slope for the energy spectral density, decaying as $\sim f^{-5/2}$, is recovered. This suggests that different physical processes operate inside this dynamical regime occurring on small scales. From a fractal point of view, magnetic field fluctuations at kinetic/dissipative scales seem to behave as a monofractal system with a Hurst exponent (i.e., $\zeta \left(1\right)$) ∼ 0.8 [52,53], quite similar for all magnetic field components (Figure 4, middle panel). This behavior is confirmed by the singularity spectrum $f\left(\alpha \right)$ (see Figure 4, lower panel) which collapses near the point $(\alpha ,f(\alpha \left)\right)=(1,1)$. Our findings suggest the absence of intermittency at dissipative scales, which is well in agreement with previous works where generalized Hilbert spectra were used [52,53,54].

The difference in the scaling properties between the two ranges of scales can be linked to the different physical processes operating in both the inertial and dissipative domains. On one hand, the inertial range is characterized by the nonlinear interactions between eddies of different size, causing their fragmentation to smaller and smaller ones until viscous effects become dominant (it is worthwhile to remark that eddies must not be thought of as real vortices, but as a description of the triadic interaction between modes). Conversely, when approaching the Kolmogorov scale $\eta$ wave-particle mechanisms and small-scale structures (like current sheets) become the most prominent features which characterize the dissipative processes [31]. Indeed, while the inertial range physics is mostly dominated by large-scale phenomena like plasma instabilities and it is characterized by an inhomogeneous nonlinear transfer of energy, resulting in the generation of localized small-scale structures with scale-dependent features [8,14,36], the dissipative range physics is mainly characterized by several dispersive phenomena generated by velocity-space effects and electron dynamics, driven by wave-wave coupling, scattering processes, and damping mechanisms [1,31,34,35].

Moreover, our results seem to confirm the robustness of the EMD-based method in investigating scaling features of solar wind fluctuations. In addition, by using the EMD-DAMF approach we are able to carry out structure function analysis on both positive and negative q, allowing us to derive the whole singularity spectrum $f\left(\alpha \right)$ such that accurate intermittency measures can be found. Conversely, generalized Hilbert spectra, unless based on the EMD and HSA procedures, cannot be evaluated for $q<0$, thus permitting only a partial detection of singularities (only the increasing branch of $f\left(\alpha \right)$ can be obtained) [55]. Although the difference in the intermittent properties between MHD/inertial and kinetic/dissipative domains remains an open question (e.g., [35]), our results can help to accurately measure scaling exponents and singularities with fewer a priori mathematical assumptions with respect to previous analysis, thus providing useful constraints for modeling purposes.

5. Chaotic Measures and Phase-Space Analysis

A dynamical system, like the solar wind, can be also investigated following a chaotic approach, mostly based on looking at the dimensionality of its phase-space. A system is defined to be chaotic if its dimension is a non-integer value [56]. Different measures have been introduced to quantify the presence and degree of chaos [57]. Particularly, the correlation dimension $D_2$, useful for determining the fractional dimensions of fractal objects, is estimated by embedding a time series $x\left(t\right)$ in a time-delayed m-component state vector as

$${\mathbf{X}}_k=\{x_1\left(t_k\right),x_2\left(t_k\right),\cdots ,x_m\left(t_k\right)\}$$

(17)

where $x_l\left(t_k\right)=x(t_k+(l+1)\Delta )$, m is usually named embedding dimension, and $\Delta$ is a time delay. Then, the correlation integral can defined as

$$C(\rho ,m)=\underset{N_s\to \infty }{lim}\frac{1}{N_s^2}\sum _{i=1}^{N_s}\sum _{j=1}^{N_s}\mathsf{\Theta }(\rho -|{\mathbf{X}}_i-{\mathbf{X}}_j\left|\right)$$

(18)

where $N_s$ is the number of considered phase-space states, $\mathsf{\Theta }$ is the Heavyside step function, and $\rho$ is the phase-space threshold distance between two points. If $\rho \to 0$, a power-law behavior is found for the correlation integral as $C(\rho ,m)\sim {\rho }^{D_2}$, where $D_2$ is defined as

$$D_2=\underset{\rho \to 0}{lim}\frac{logC(\rho ,m)}{log\rho }.$$

(19)

As the embedding dimension m increases, the correlation dimension will converge to its true value. Specifically, if $D_2=m$ then the system will explore the whole phase-space; conversely, if $D_2<m$ a strange attractor will characterize the phase-space dynamics. Of course, both m and $\Delta$ need to be properly chosen. Their choice is crucial for a correct estimation of the correlation dimension in the case of chaotic systems [56,57,58]. Generally, the choice of the time delay $\Delta$ corresponds to the first minimum of the autocorrelation function of the time series, while the choice of the embedding dimension m falls on the lowest value at which $D_2$ approaches from a constant value [57].

Figure 5 shows the behavior of the correlation dimension $D_2\left(\tau \right)$ as a function of the mean frequency of each empirical mode, derived as the inverse of the mean timescale $\tau$. This allows us to investigate how the dynamical behavior changes when moving from MHD/inertial to kinetic/dissipative scales.

The dimensionality of the system clearly exhibits a scale-dependent behavior characterized by an increase in the values of $D_2\left(\tau \right)$ with the mean frequency, approaching from a constant value $D_2\left(\tau \right)\sim$ 2.7 for $f>f_b=0.4$ Hz. This suggests that magnetic field fluctuations are characterized by a superposition of processes working on different timescales and with different dimensionality. While MHD processes can be described by using a low-dimensional dynamical system, since $D_2<$ 2, the kinetic/dissipative domain dynamics cannot be represented as a linear system since at least three system variables ($D_2>$ 2) are needed to describe processes (perhaps dissipation) occurring at these scales.

An interesting result is the continuous change of the correlation dimension moving from MHD/inertial to kinetic/dissipative domains, which suggests that a single correlation dimension is not capable of describing the complexity features of solar wind magnetic field fluctuations at the MHD/inertial scales, while a single correlation dimension seems to describe kinetic features. This can be interpreted as the signature of the intermittent nature of fluctuations in the MHD/inertial domain, where a hierarchy of dimensions is necessary to describe the complex nature of the nonlinear energy cascade. Conversely, at the kinetic/dissipative scales where dissipation may occur, the correlation dimension seems to converge to a single value of $D_2\sim$ 2.7. This is the temporal counterpart of the multifractal nature of turbulence in the MHD/inertial domain and of the monofractal nature of the dissipative regime, as also shown in Section 4.

We can characterize a dynamical system by looking at its phase-space dynamics in order to recover the existence of fixed points and their nature, as well as to investigate the presence of (strange) attractors [57]. Since by using the EMD we are able to decompose our time series into oscillating functions [48,59,60,61], we choose to reconstruct empirical modes according to the different dynamical regimes investigated. We can investigate the dynamics at the MHD/inertial and kinetic/dissipative scales in a separate way by defining

$$\begin{array}{ccc}\hfill R_I\left(t\right)& =& \sum _{f_k\in f_I}{c}_k\left(t\right)\hfill \end{array}$$

(20)

$$\begin{array}{ccc}\hfill R_D\left(t\right)& =& \sum _{f_k\in f_D}{c}_k\left(t\right)\hfill \end{array}$$

(21)

as the reconstructions of empirical modes with characteristic mean frequencies inside the intertial ($f_I$) and kinetic/dissipative ($f_D$) domains. In detail,

$$\begin{array}{ccc}\hfill f_I& \doteq & \left\{f_k|{10}^{-3}\mathrm{Hz}<f_k<f_b\right\}\hfill \end{array}$$

(22)

$$\begin{array}{ccc}\hfill f_D& \doteq & \left\{f_k|f_k>f_b\right\}\hfill \end{array}$$

(23)

being $f_b=0.4$ Hz.

Figure 6 reports the phase-space portraits for the two different dynamical regimes, i.e., inertial scales (left panels) and kinetic/dissipative ones (right panels). The different symbols identify different phase-space trajectories starting at different phase-space positions, identified by a black symbol, and ending with a magenta one. The results look quite interesting and can be interpreted in dynamical system framework.

The dynamics seems to be characterized by an unstable orbit at inertial scales, so that the associated fixed point can be classified as a saddle. Indeed, starting from different phase-space positions each trajectory moves along an unstable manifold such that the system will approach the (unstable) fixed point being repelled on different (and opposite) phase-space points. Thus, the set is a repeller. This hyperbolic equilibrium point does not have any center manifolds, and, near it, the orbits of the system resemble hyperbolas. Conversely, the dynamics at kinetic/dissipative scales is characterized by a set which is an attractor since all phase-space trajectories tend to move towards the stable fixed point, which can be identified as a node. This fixed point, due to its fractal dimension and structure (see Section 4), is a chaotic strange attractor, extremely sensitive to initial conditions. By considering two arbitrarily close initial phase-space positions near the attractor, after several time steps they will move on phase-space positions far apart, and after other several time iterations will lead to phase-space positions which are arbitrarily close together. Thus, the dynamics never depart from the attractor [56].

The obtained results seems to suggest a new view of the dynamics of the solar wind at different scales from the MHD/inertial domain down to the kinetic/dissipative one. The system undergoes a saddle-node bifurcation, a local bifurcation in which two fixed points collide and annihilate each other, with an unstable fixed point (saddle) and a stable one (node). This means that both the inertial and kinetic/dissipative ranges can be seen as fixed points of the governing system equations, one unstable and the other stable. In this way, the phenomenological model of the Richardson cascade [24,26,28,29] can be interpreted in the different context of the dynamical system theory. Energy is injected at a scale L, which represents a stable fixed point of the system; then, when nonlinear interactions develop, corresponding to changes in one or more dynamical bifurcation parameters, the dynamics of the system changes, moving towards an unstable fixed point (i.e., the inertial regime) which, due to its repeller nature, forces the system to explore the available phase-space until a stable fixed point (i.e., the kinetic/dissipative domain) is reached (Figure 7). In a simple conceptual model, a bifurcation parameter could be the timescale of the different processes operating inside the MHD/inertial and kinetic/dissipative domains such that the dynamics of the system, represented in our case by the magnetic field components $B_i$, can be seen as solely dependent on $\tau$

$$\begin{array}{c}\hfill {\dot{B}}_i=g(B_i,B_j,\tau ).\end{array}$$

(24)

The system is surely characterized by a chaotic dynamics due its dimension (our system is described by three variables) with nonlinear interactions between the different variables (i.e., $g(B_i,B_j,\tau )$ is a nonlinear function of $B_i,B_j$) as required when describing turbulent features (e.g., [1,24,26]).

6. Conclusions

Solar wind magnetic field fluctuations at different scales have been investigated by employing both a multifractal and a chaotic approach. The multifractal analysis has been performed by using a novel formalism, the EMD-based dominant amplitude multifractal formalism, through which increments are derived by using local properties of fluctuations at different scales obtained by using the Empirical Mode Decomposition method. The results suggest that MHD fluctuations show an intermittent character, well described in the framework of classical multifractal models (like the p-model (e.g., [26])); conversely, magnetic field fluctuations at kinetic scales (i.e., beyond the ion inertial length) show a monofractal behavior, in agreement with previous findings (e.g., [52,53,54]).

The phase-space dynamics of the two ranges of scales, i.e., inside the MHD/inertial and kinetic/dissipative domains, is characterized by a different degree of chaos, because the system is more chaotic when moving from the MHD down to the kinetic scales. An unstable manifold is recovered for the MHD scales, characterizing an unstable saddle for the magnetic field dynamics. Conversely, a stable manifold, corresponding to a stable node, is found at kinetic scales, suggesting the occurrence of a saddle-node bifurcation passing from MHD down to kinetic scales. These results can open the way to new perspectives in approaching scale-to-scale dynamics of solar wind magnetic field fluctuations as well as in deriving conceptual models to explain the observed dynamical regimes.