© ESO 2021

1. Introduction

Investigations of the turbulent nature of solar wind fluctuations have been ongoing for more than half a century (see, e.g., Bruno & Carbone 2016). Advances have been made consistently thanks to the increasingly accurate measurements of several dedicated space mission as well as to the enormous improvement of numerical calculation, new detailed models and theoretical frameworks, and the development of specific data analysis techniques. Nevertheless, the extremely complex nature of the system and the coexistence of multiple actors, scales, and dynamical regimes have led to a number of questions that remain open (Viall & Borovsky 2020). Among these, the very nature of the turbulent cascade of the solar wind flow and its relationship with the small-scale processes still need to be described in full (Tu & Marsch 1995; Bruno & Carbone 2013; Matthaeus & Velli 2011; Chen 2016). Magnetic field fluctuations have been characterized with great detail at magnetohydrodynamic and kinetic scales, for example, through spectral and high-order moments analysis (Tu & Marsch 1995; Bruno & Carbone 2013). The anisotropic nature of magnetic turbulence has also been addressed, and is still being debated, due to the limited access to three-dimensional measurements in space (see, e.g., Horbury et al. 2008, 2012; Sorriso-Valvo et al. 2010; Yordanova et al. 2015; Verdini et al. 2018; Telloni et al. 2019a; Oughton & Matthaeus 2020). Velocity fluctuations have been studied thoroughly (see, e.g., Sorriso-Valvo et al. 1999; Bruno & Carbone 2013), although the kinetic scales still remain quite unexplored for instrumental limitations, most notably in the sampling time resolution. Both the velocity and magnetic field show highly variable turbulence properties, with well developed spectra, strong intermittency (Sorriso-Valvo et al. 1999), anisotropy, and linear third-order moments scaling (Sorriso-Valvo et al. 2007; Carbone et al. 2011). The level of Alfvénic fluctuations (mostly but not exclusively found in fast streams, see e.g., D’Amicis et al. 2011; Bruno et al. 2019) are believed to be associated with the state of the turbulence. In particular, solar wind samples containing more Alfvénic fluctuations are typically associated with less developed turbulence, as inferred from both shallower spectra and reduced intermittency (see Bruno & Carbone 2013, and references therein). This is consistent with the expectation that uncorrelated Alfvénic fluctuations contribute to reduce the nonlinear cascade by sweeping away the interacting structures (Dobrowolny et al. 1980), as also confirmed by the observed anticorrelation between the turbulent energy cascade rate and the cross-helicity (Smith et al. 2009; Marino et al. 2011a,b).

Conversely, density fluctuations have been only partially explored, due in part to the unavailability of high-frequency time series and in part to their supposedly secondary relevance in the nearly incompressible solar-wind dynamics. In recent years, studies have shown that proton density is also turbulent and intermittent (Hnat et al. 2003, 2005), with the characteristics of the turbulence often depending of the Alfvénic nature of each specific solar wind interval. In particular, it is understood that in the nearly incompressible Alfvénic solar wind, the turbulence of density fluctuations is strongly similar to that of magnetic field magnitude, as both behave as scalar quantities passively advected by the turbulent fluctuations of velocity and magnetic field components (Goldreich & Sridhar 1995; Chen et al. 2012). Conversely, in the more compressive non-Alfvénic solar wind (more typically associated with slow streams), they are actively contributing to the nonlinear transfer of energy (Schekochihin et al. 2009; Boldyrev et al. 2013). The radial evolution of proton density turbulence has been examined, providing evidence of complex, unexpected behaviours (Bruno et al. 2014). A modeling of the radial evolution was attempted based on the parametric instability, expected to generate increasingly compressive fluctuations as the solar wind expands away from the sun. However, despite the analysis provided some insight, the lack of statistical description still prevents the validation of models that could provide a prediction for the radial evolution of the density fluctuations. Finally, sub-ion scale density turbulence has been investigated using high-resolution proton density measurements from Spektr-R (Chen et al. 2014; Riazantseva et al. 2019) and electron density from MMS (Roberts et al. 2020), but the results were not be fully conclusive regarding the nature of the multifractal and intermittency properties of the fluctuations (Sorriso-Valvo et al. 2017; Carbone et al. 2018). Therefore, a deeper analysis of solar wind density fluctuations is necessary in order to constrain modelings and to understand the role of density fluctuations in the magnetohydrodynamic turbulent cascade. We note that although the above studies may refer either to proton or electron density, depending on the instrument used for the measurements, at the scales of interest for turbulence, the two can be safely considered equal for the quasi-neutrality condition of solar wind plasma.

Recent developments in studying the electromagnetic wave payload have allowed for the density fluctuations to be obtained from the measurement of the spacecraft potential (Pedersen 1995; Roberts et al. 2017). This has provided access to higher-frequency and more accurate density measurements. In particular, the Solar Orbiter spacecraft (Müller et al. 2020) was launched in February 2020, equipped with both remote and in situ instrumentation, aimed at investigating the Sun and the solar wind from and within the inner heliospehre (Zouganelis et al. 2020). The Radio and Plasma Waves (RPW) instrument (Maksimovic et al. 2020) measures the spacecraft potential with unprecedented accuracy, allowing the estimation of the high-cadence solar wind electron density (Khotyaintsev et al. 2021). Together with the enhanced data quality provided by Solar Orbiter instruments, novel data analysis techniques are emerging to improve the delicate measurement of turbulence parameters in solar wind data. Among these, Empirical Mode Decomposition (EMD, Huang et al. 2008) has been successfully applied to mitigate short-sample and large-scale structures effects (Carbone et al. 2018; Alberti et al. 2019). This paper is aimed at providing the first description of the properties of turbulence as obtained applying EMD-based analysis techniques to the Solar Orbiter RPW density measurements taken over the month of June 2020 (Khotyaintsev et al. 2021).

The results are discussed in relation to solar wind parameters and magnetic field turbulence, in the context of a study applying the same EMD-based techniques to the MAG instrument (Horbury et al. 2020). Moreover, after introducing a novel parameter to describe the presence of ion-scale waves, we analyzed the relationship between the properties of turbulence and such waves.

Section 2 provides a description of the data used for the analysis. In Sect. 3, the techniques used for the analysis of turbulent fluctuations and the corresponding results are described. In Sect. 4, the wave parameter is introduced and the presence of ion-scale waves is discussed. Section 5 gives a detailed discussion on the existing correlations between solar wind, turbulence, and wave parameters along with their implications. Finally, our conclusions are summarized in Sect. 6.

2. Description of data

In order to study the properties of turbulence of solar wind density fluctuations, we made use of Solar Orbiter measurements taken from 7 to 29 June 2020, when the spacecraft was orbiting the Sun between 0.52 AU and 0.55 AU. During that time interval, the Solar Wind Analyser (SWA) plasma instruments (Owen et al. 2020) were not operational, so that direct measurements of proton and electron moments are not available. Here, we use the 16 Hz electron density, n_e, accurately estimated from the RPW probe-to-spacecraft potential measurements (Maksimovic et al. 2020). The equilibrium floating potential of the spacecraft is reached when the current due to photo-electrons emitted from the spacecraft is balanced by the plasma current collected by the spacecraft. The equilibrium is reached instantaneously (0.1–1 ms) on the time scales of interest for the turbulence studies. From the current balance, we can find that the electron density has an exponential relation to the spacecraft potential. Then, by making an exponential fit of the spacecraft potential to the electron density obtained from the high-frequency measurement of plasma quasi-thermal noise, we can find the relation for approximate conversion of the probe-to-spacecraft potential to electron density. Details on the density estimation technique are given in Khotyaintsev et al. (2021).

Localized estimates of the solar wind speed from the deHoffmann-Teller (HT) analysis of electromagnetic fluctuations are used to provide general context, most notably discriminating between fast and slow solar wind streams. The HT analysis is used to find the velocity of the frame in which the electric field is zero. In the solar wind, where current sheets and magnetohydrodynamic (MHD) turbulence are ubiquitous, the HT velocity is in general a good estimate of the solar wind velocity. In order to get reasonable coverage, estimates of the solar wind speed V_sw were obtained every ten minutes by applying the HT analysis on one-hour running intervals of E and B data. Details about the HT analysis can be found in Steinvall et al. (2021). Finally, magnetic field vector, B, measurements from the magnetometer (MAG; Horbury et al. 2020) were studied for the purposes of comparison.

Due to the presence of large gaps in the electron density data and to avoid any strong violation of homogeneity and stationarity, a number of intervals of variable length were extracted from the electron density time series n_e(t). All intervals were chosen as relatively stationary, covering at least 1 h of data and were limited to include only a few data gaps that are shorter than 5 s. The remaining missing points have been interpolated linearly. Such criteria resulted in the selection of 36 sub-intervals, whose list and macroscopic details (starting and ending time, interval duration, mean distance from the Sun, estimated mean solar wind bulk speed, and angle between the ambient magnetic field and the radial direction averaged over each interval) are given in Table 1. Given the estimated solar wind speed, the Taylor hypothesis (Taylor 1938) is considered valid throughout the whole data set, allowing us to interpret the time series as equivalent to a space ensemble, therefore enabling the standard tools for turbulence analysis.

Figure 1 shows an overview of the above parameters for the whole month of June 2020. The overall solar wind conditions are variable, with an alternation of fast and slow streams, as well as complex coronal mass ejection structures (as studied in detail by Telloni et al. 2021). However, the selected sub-intervals, whose durations vary between one and six hours, are typically embedded in homogeneous solar wind conditions.

Fig. 1. Overview of the data. From top to bottom: solar wind electron density ne (RPW), deHoffmann-Teller solar wind velocity estimate Vsw (RPW) and interplanetary magnetic field components, Bi, and magnitude, |B| (MAG), measured by Solar Orbiter during the whole month of June 2020.

3. Analysis of solar wind electron density turbulence

The properties of turbulent fields are usually studied through standard statistical analysis techniques. Among others, these may include (Dudok de Wit et al. 2013): the autocorrelation function, used to determine specific scales of the data; the power spectral density E(f), providing information on the self-similar energy redistribution among scales; and the structure functions (SF) $S_q({\mathcal{l}}_t)\equiv \langle {[n_{\text{e}}(t+{\mathcal{l}}_t)-n_{\text{e}}(t)]}^q\rangle \sim {\mathcal{l}}_t^{\zeta (q)}$, with ℓ_t representing the time scale of the field increments, whose anomalous scaling exponents ζ(q) are able to quantitatively describe the effects of intermittency; namely, the inhomogeneous nature of the turbulent cascade (Frisch 1995). Related to the structure functions, the kurtosis of the distribution function of the fluctuations $K({\mathcal{l}}_t)=S_4({\mathcal{l}}_t)/S_2^2({\mathcal{l}}_t)$ is often used to quantify the intermittency effects. However, such techniques may be sensitive to the data sample characteristics, resulting in undesired effects not attributed to the turbulent energy cascade, but rather due, for example, to finite sample size, limited stationarity, or presence of superposed structures larger than the typical turbulence scales (for example non-turbulence related current sheets or velocity shears). Such issues often occur in ecliptic solar wind intervals, where instrumental performance and the intrinsic wind variability may prevent ideal experimental conditions for the study of turbulence (Carbone et al. 2018). The magnification shown in the top two panels Fig. 2 shows an example of selected intervals (specifically, sub-intervals 26 and 27, both on June 20). The presence of large-scale modulation of the density profile, sporadic sharp gradients similar to ramp-cliff structures, and the general lack of strict stationarity appear evident (Matthaeus & Goldstein 1982; Perri & Balogh 2010). In order to mitigate the effects of such large-scale features, we made use of the Hilbert-Huang transform to obtain more precise estimators of the generalized high-order spectra.

Fig. 2. Selected intervals. Two examples of adjacent intervals during day 20 June 2020. Left: interval 26; right panel: interval 27. From top to bottom: interplanetary magnetic field components, Bi, and magnitude, |B| (MAG); electron density ne (RPW); magnetic filed spectrogram; perpendicular magnetic field components coherence; perpendicular magnetic field components phase angle; wave parameter Qw(t) (see Sect. 4).

3.1. Empirical mode decomposition

The technique used here is based on EMD (Huang et al. 1998; Jánosi & Müller 2005; Carbone et al. 2016b). This is a self-consistent, data-driven projection of a time series (as in this case the solar wind electron density n_e(t)) on a finite number n of empirical basis functions ϕ_j(t), called intrinsic mode functions (IMFs), so that

$$\begin{array}{c}\hfill n_{\mathrm{e}}(t)={\displaystyle \sum _{j=1}^n{\varphi }_j(t)+r_n(t).}\end{array}$$(1)

The additive residual function r_n(t) describes the mean trend. Each IMF can be characterized by the instantaneous timescale τ_j(t). The decomposition is based on a recursive procedure, consisting of two main stages (Huang et al. 1998): (i) the local extrema of n_e(t) are interpolated through cubic spline to provide superior and inferior envelops of the time series; and (ii) the average between the two envelops is subtracted from the original data. The resulting field is accepted as an IMF if it satisfies the following specific criteria: the number of local extrema and zero crossings does not differ by more than one, and the average between the IMF superior and inferior envelopes is zero at all times. Otherwise, the procedure is repeated on the remaining IMF until the criteria are met according to the so called 3-thresholds stoppage criterion (with the following standard choice of parameters: δ = 0.05, ξ₁ = 0.05, and ξ₂ = 10ξ₁, see Rilling et al. 2003).

Some examples of IMFs and the associated residual, as obtained through the above EMD decomposition of the solar wind electron density, n_e, are shown in the left panels of Fig. 3 for sample 1. For clarity, only odd IMFs have been plotted. Although the instantaneous frequency of the modes is variable, the figure highlights that each mode, j, has a characteristic narrow frequency band, so that a mean period, ⟨τ_j⟩, and an associated variance, Var(ϕ_j), can be properly defined. Indeed, EMD acts intrinsically as a dyadic filter bank (Wu & Huang 2004; Flandrin et al. 2004; Flandrin & Goncalves 2004; Huang & Shen 2005), each IMF effectively capturing a narrow frequency band. However, the general features of the various IMFs depend on the specific process under analysis. For a turbulent field (Huang et al. 2008; Carbone et al. 2016a, 2018, 2020a) or for a multifractal process (Carbone et al. 2010; Sorriso-Valvo et al. 2017), the characteristic mean period grows exponentially as ⟨τ_j⟩=αγ^j, where ⟨ ⋅ ⟩ represents an ensemble average (in this case, average over time). The basis, γ, can be evaluated empirically from the IMFs; for an exact dyadic decomposition, γ = 2. Additionally, the variance of the IMF scales as a power of the mean timescale, Var(ϕ_j)∼⟨τ_j⟩^2ℋ. The scaling exponent, ℋ, is the Hurst number, a parameter describing the persistence or anti-persistence of the fluctuations of the process under analysis (Nava et al. 2016; Carbone et al. 2019).

Fig. 3. Results of the EMD analysis. Left panel: example of IMFs ϕj(t) for sample 1 (black lines). The bottom plot (red line) represents the residual rn(t). For better readability, only odd IMFs have been plotted. Central panel: average timescale ⟨τj⟩ of each IMF of sample 1 as a function of the mode j. Error bars represent the 95% confidence bounds. The dashed line is a least square fit obtained from the relation ⟨τj⟩=αγj. Right panel: IMF variance Var(ϕj) as a function of the average period ⟨τj⟩ for sample 1, with the associated 95% confidence bounds. The dashed line represents the relation Var(ϕj)∼⟨τj⟩2ℋ, being ℋ the Hurst number (see the list for all intervals in Table A.1).

The scaling of the mean period for sample 1 is shown in the central panel of Fig. 3. The value γ = 1.88 ± 0.11, obtained through a least-square exponential fit, is compatible with the expected value for a dyadic decomposition, showing that EMD is correctly decomposing the data. Moreover, the Hurst number, ℋ = 0.35 ± 0.04, is obtained from the regression of the IMFs variance versus the average period, as shown in the right panel of Fig. 3. In this example, ℋ is compatible with the standard value for classical fluid or magnetic turbulence, ℋ = 0.37, (Benzi et al. 1993; Arneodo et al. 1996; Bruno & Carbone 2016) and indicates persistence of the fluctuations typical of intermittency.

Using the above procedure, the Hurst number has been evaluated for all 36 intervals. The results, listed in Table A.3, show some variability that is discussed and compared to other parameters in Sect. 5.

3.2. Arbitrary-order Hilbert Spectral analysis

The Hilbert Spectral Analysis (HSA) is an extension of the basic EMD designed to characterize scale-invariant properties directly in the amplitude-frequency space (Huang et al. 2008). It provides the equivalent of the power spectral density and high order moments of a field fluctuations (the structure functions), therefore representing a viable alternative to those standard tools. After decomposing the field under study in its IMFs, the Hilbert transform of each mode is computed as:

$$\begin{array}{c}\hfill {\varphi }_j^{\star }(t)=\frac{1}{\pi }P{\displaystyle {\int }_{-\infty }^{+\infty }\frac{{\varphi }_j(t^{\prime })}{t-t^{\prime }},}\end{array}$$(2)

where P is the Cauchy principle value. The Hilbert representation allows to extract a time-dependent instantaneous frequency $f_j(t)\equiv {\tau }_j^{-1}(t)$ and a time-dependent amplitude modulation 𝒜_j(t), by constructing the so called analytical signal ${\mathrm{\Phi }}_j(t)={\varphi }_j(t)+i{\varphi }_j^{\star }(t)\equiv {\mathcal{A}}_j(t)e^{i{\theta }_j(t)}$ (Cohen 1995). Here, ${\mathcal{A}}_j(t)=|{\mathrm{\Phi }}_j(t)|\equiv \sqrt{{\varphi }_j^2(t)+{{\varphi }_j^{\star }}^2(t)}$ and ${\theta }_j(t)=\arctan [{\varphi }_j^{\star }(t)/{\varphi }_j(t)]$ are the instantaneous amplitude modulation and the instantaneous phase oscillation, respectively (the instantaneous frequency being defined as f_j(t) = 2π⁻¹dθ_j(t)/dt) (Long et al. 1995; Cohen 1995; Flandrin 1999).

After rewriting the original signal in terms of 𝒜_j and θ_j, n_e = Real[∑_j𝒜_j(t)exp(i∫f_i(t)dt)], the energy as a function of the instantaneous frequency f and time, can be defined as $h(f)={\int }_0^{\infty }H(f,t)\text{d}t$, or the marginal integration of the Hilbert spectrum H(f, t)≡𝒜²(f, t) (being H(f, t) a representation of the original signal at the local level) (Huang et al. 1998, 2009). An equivalent definition of H(f, t) can be obtained from the joint probability density function of the instantaneous frequency and amplitude P(f_j, 𝒜_j) (Long et al. 1995), extracted from the IMFs. In this case, the Hilbert marginal spectrum is the second statistical moment of such distribution, analog to the Fourier spectral energy density:

$$\begin{array}{c}\hfill h(f)\equiv {\mathcal{L}}_2(f)={\displaystyle {\int }_0^{\infty }P(f,\mathcal{A}){\mathcal{A}}^2\mathrm{d}\mathcal{A}.}\end{array}$$(3)

The above definition can be then generalized to any arbitrary moment q ≥ 0, representing the analogous of the standard structure functions of the fluctuations:

$$\begin{array}{c}\hfill {\mathcal{L}}_q(f)={\displaystyle {\int }_0^{\infty }P(f,\mathcal{A}){\mathcal{A}}^q\mathrm{d}\mathcal{A}.}\end{array}$$(4)

Equations (3) and (4) are used here to estimate the spectral and intermittency properties of the electron density turbulence of the 36 sub-intervals of this study.

The left panel of Fig. 4 shows one example of the equivalent spectrum ℒ₂(f) (blue circles), obtained through the HSA described above, for sample 1. For comparison, the classical power spectral density E(f) (red line), evaluated through the Welch’s method (Welch 1967), and the second-order structure function $S_2({\mathcal{l}}_t^{-1})$ are also shown. Both spectra clearly display power-law scaling E(f),ℒ₂(f)∼f^−β₂, in the frequency range of f ∈ [3 × 10⁻³, 10⁻¹], which is compatible with the typical inertial range of timescales. Similarly, the second-order structure function also show a very clear power law scaling. In all three cases, the power-law scaling of ℒ₂(f) is compatible with the standard Kolmogorov predictions ($S_2({\mathcal{l}}_t)\sim {\mathcal{l}}_t^{2/3}\to E(f)\sim f^{-5/3}$), also shown in the left panel of the figure.

Fig. 4. Results of the Hilbert-Huang analysis. Left panel: second-order Hilbert spectrum ℒ2(f) (blue circles), the classical Fourier PSD (Welch 1967) E(f) (red line), and the second-order structure function S2(ℓt) (green squares, plotted as a function of the inverse timescale ℓt), for sample 1. Power scaling is present in the same frequency range for all methods. The dashed line represents the classical Kolmogorov scaling ℒ2(f)∼f−5/3, while the dotted line represent the expected scaling for the second-order structure function $S_2({\mathcal{l}}_t)\sim {\mathcal{l}}_t^{\zeta (2)}$ (ζ(2) = β2 − 1). Central panel: generalized Hilbert spectra ℒq(f) for q ∈ [1, 5], obtained for sample 1. The curves have been vertically shifted for clarity. The shaded area represents the frequency range of the bootstrapping least-square fit. Right panel: same as in the central panel but for sample 9. The power behavior is still present, but the power law exponents are considerably different.

The HSA equivalent spectrum ℒ₂(f) displays a slightly better power-law scaling than the traditional Fourier spectral density, which has some weak amplitude modulation. Thanks to the local nature of the EMD analysis, the sources of such modulation can be removed, allowing to obtain a more precise determination of the spectral scaling exponents. This corresponds to isolating the properties of the turbulent cascade from the possible effects of the instrumental noise, and of the larger-scale energy inhomogeneity (Huang et al. 2010; Carbone et al. 2018; Telloni et al. 2019a). We note that the power-law scaling range can vary for the various samples, but always includes at least one decade of scales. The two spectral estimators are not always nicely superposing, with the Fourier spectrum occasionally presenting stronger modulations (not shown). This suggests the possible presence of large-scale modulations, which may affect the Fourier spectra but are controlled well by the HSA. In order to achieve robust estimate, the scaling exponent of the qth Hilbert spectrum ℒ_q(f) was evaluated via the residual resampling (bootstrapping) procedure (Bradley & Robert 1994; Carbone et al. 2020b). First, a least square fit is performed on each ℒ_q(f), then the residuals are randomly resampled and added to the fit, generating a new data-set (replica). The replica is then fitted again and the procedure is repeated a number of times; in this case, N_boot = 10⁴ (Boos & Stefanski 2010; Wilcox 2010). Such a large number of replications is necessary for correctly evaluating confidence or prediction interval, whereas, in general, simple statistical tests require a smaller number, namely, N_boot ∈ [50, 100] (Dale 2012). The probability distribution of the exponents β_q (or the scaling exponents β_q − 1) obtained from the N_boot least-square fits is then used to estimate the 50th percentile (median of the distribution), applied as the best estimate of the exponents, along with the statistical error (95% confidence interval) (Efron et al. 2015; Wilcox 2010). Figure 5 shows two examples of the distribution P(β₂) estimated for samples 1 and 9 using the bootstrapping technique. The red dotted vertical line represents the median, and the black dashed lines represents the statistical error. In the examples shown in the figure, the median of the distribution provides the scaling exponent β₂ = 1.72 ± 0.14 for sample 1 (error indicating the 95% confidence interval), in excellent agreement with the classical scaling for the fully developed hydrodynamic turbulence β₂ ≈ 1.7 (Benzi et al. 1993), and a shallower β₂ = 1.51 ± 0.10 for sample 9, consistent with the Iroshnikov-Kraichnan spectrum for Alfvénic turbulence (Iroshnikov 1964; Kraichnan 1965). Values of β₂ were obtained for all 36 samples. A discussion of these values, presented in Table A.1, is provided in next section and in Sect. 5.

Fig. 5. Probability distribution function P(β2) of the scaling exponents β2, constructed via bootstrap resampling, for sample 1 (upper panel) and sample 9 (lower panel). In both panels, the vertical dotted line represents the median of the distribution (50th percentile), while the vertical dashed bars indicate the 95% confidence interval (enclosed between 2.5th and 97.5th percentile).

3.3. Scaling exponents and intermittency analysis

In the framework of the standard Kolmogorov turbulence, a direct link exists between the Fourier spectral exponent E ∼ f^−β and the scaling exponent of the second-order structure function, $S_2\sim {\mathcal{l}}_t^{\zeta (2)}$, so that β − 1 = ζ(2). This relationship can be extended to any moment order q of the generalized Hilbert spectra ℒ_q ∼ f^−β_q, yielding the generalized scaling exponents ξ(q)≡β_q − 1 (Huang et al. 2010; Carbone et al. 2016a). These are analogous to the scaling exponents, ζ(q), obtained using the standard structure functions (Frisch 1995; Benzi et al. 1993; Arneodo et al. 1996), and can be used to retrieve quantitative information on the properties of turbulence. Additionally, the structure function scaling exponents are linked to the Hurst number via the relation ζ(q) = qℋ (in absence of intermittency corrections). This allows an alternative estimate of the Hurst number using, for example, the first-order exponent ℋ = β₁ − 1 ≡ ξ(1).

The central and right panels of Fig. 4 show two examples of ℒ_q(f) (for intervals 1 and 9), obtained from Eq. (4) up to the fifth order, for the electron density, n_e, in samples 1 and 9. All curves presents good power-law scaling for all orders, approximately in the inertial range of frequencies (shaded areas). The associated generalized Hilbert spectra scaling exponents, β_q, and, hence, the equivalent structure-function scaling exponents, ξ(q), were obtained through the bootstrapping procedure described above and were used to determine the intermittency properties of the electron density. In order to check the quality of the procedure, the first-order exponents were initially used to obtain the alternative estimate of the Hurst number. For the example of sample 1, the value ℋ = 0.32 ± 0.06 was obtained, which is in good agreement with the value obtained through the regression of the IMF variance versus the average timescale, illustrated in Fig. 3. This was consistently observed for all the 36 samples. The power-law exponents of ℒ_q(f) are visibly different in the two intervals shown in the central and right panels of Fig. 4. For q = 2, this can be also more quantitatively noticed by comparing the distributions and median values shown in Fig. 5 for the same intervals. The power-law exponents of ℒ_q(f) are visibly different in the two intervals shown in the central and right panels of Fig. 4.

The scaling exponents, ξ(q), shown in the top-left and bottom-left panels of Fig. 6 for most of the samples, are separated in two groups, as described in the following. The scaling exponents for the magnetic field magnitude, |B|, (central panels) and radial component, B_r, (right panels) are shown for comparison. For all cases included in the figure, the curvature of the exponents with respect to the linear prediction of the non-intermittent Kolmogorov phenomenology is evident. This universal behaviour of turbulence is due to the effects of the intermittency or any anomalous dissipation (Kolmogorov 1962; Schmitt et al. 1994; Schmitt 2003; Bruno & Carbone 2016; Carbone et al. 2019). These effects are related to the multifractal nature of the turbulent cascade (Meneveau & Sreenivasan 1991; Davis et al. 1994; Sorriso-Valvo et al. 2017). In all panels, the scaling exponents from a classical measure of fluid intermittent turbulence are also shown for a comparison with the results of Benzi et al. (1993). In order to describe their intermittency behaviour, the 36 samples were then separated in three groups, according to the behaviour of the generalized scaling exponents and of the Hurst number, with respect to the standard fluid turbulence reference. We note that for each sample, the same behaviour is consistently observed for density and magnetic field.

Fig. 6. Results of the equivalent structure function analysis. Top-left panel: HSA equivalent scaling exponents ξ(q) for group 1 (blue circles) and their average (black crosses). The exponents ζ(q) for velocity fluctuations in the inertial range of hydrodynamic turbulence (red squares) are shown for reference (Benzi et al. 1993). The dashed line represents the theoretical expectation ζ(q) = q/3 (Kolmogorov 1941). Top-center and top-right panels: equivalent scaling exponents for the magnetic field radial component Br and magnitude |B|, for group 1. Bottom panels: same as the top panels, but for the intervals of group 2.

The first group (group 1, including 21 intervals) displays the typical statistical features of fully developed turbulence (top panels of Fig. 6). For these samples, the scaling exponents are consistent with the reference values from fluid turbulence and are well described (not shown) by models of intermittent turbulence (e.g., the p-model by Meneveau & Sreenivasan 1991, not shown). In particular, the equivalent spectral exponent β₂ > 1.55 is always compatible with the Kolmogorov scaling (β₂ = 5/3). Furthermore, for samples in this group, ℋ ∈ [0.30, 0.39], which is compatible with the classical value obtained for ordinary fluid turbulence ℋ = 0.37. From these observations, we conclude that the samples in group 1 are characterized by a standard turbulence, with the expected power-law spectra, presence of small-scale intermittent structures and anti-persistent fluctuations.

In the second group (group 2, featuring 11 samples), the scaling exponents are characterized by more extreme deviation from the linear prediction and by much smaller values (bottom panels of Fig. 6). These exponents also deviate considerably from the fluid reference and their order dependence cannot be described by standard models of turbulence (not shown). The equivalent spectral exponents β₂ < 1.55 are shallower than the Kolmogorov scaling, and in some cases they are compatible with the Iroshnikov-Kraichnan scaling. The Hurst number is also consistently smaller than for standard turbulence, ℋ ∈ [0.16, 0.26]. These observations suggest that, unlike in group 1, the density and magnetic fluctuations in these samples may not be generated by a standard turbulent cascade. Some other processes might coexist, modifying the statistics. We note that the EMD-based analysis constrains the effects of finite-size sample, poor stationarity, and large-scale structures effects. Therefore, it can be claimed that the observed features might be related to the presence of small- or inertial-scale fluctuations that are not generated only by a turbulent cascade.

The third group (group 3, not shown) includes four samples that do not show a clear power-law scaling of the generalized Hilbert spectra for all orders, so that not all the scaling exponents are available. These intervals, associated with small Hurst number and spectral exponent, are therefore not representative of any turbulence.

The multifractal nature of the fluctuations can be quantitatively described fitting the scaling exponents ξ(q) to a log-normal model (Schmitt 2003; Medina et al. 2015):

$$\begin{array}{c}\hfill \xi (q)=q\mathcal{H}-\frac{\mu }{2}(q^2-q).\end{array}$$(5)

The model is able to describe standard intermittent turbulence when the curvature parameter μ ≈ 0.02. For other multifractal processes, not generated by a nonlinear turbulent cascade, different values can be obtained. One example of log-normal model fit of the equivalent scaling exponents ξ(q) is shown in Fig. 7 for group 1 (sample 10, red circles), giving μ = 0.019 ± 0.004. The model was also fitted to the exponents obtained ensemble-averaging all samples of group 1 (stars), providing μ^⋆ = 0.023 ± 0.002. Alternatively, the average parameter computed using the results of the fit of all samples of group 1 is ⟨μ⟩₁ = 0.028 ± 0.010 (the error representing the standard deviation). All the above values are in good agreement with those observed for standard turbulence. Similar values were obtained for the magnetic field magnitude |B| (e.g., ⟨μ⟩₁ = 0.027 ± 0.010). The exponents for group 2 were also fitted to relation 5. In that case, the resulting parameters were generally smaller, with average ⟨μ⟩₂ = 0.013 ± 0.01. The parameters from all 36 intervals are plotted in one of the panels of Fig. 9 and will be discussed in Sect. 5. While there is a considerable spread in both groups, the parameters for group 2 appear generally smaller, confirming that the fluctuations have peculiar, strongly multifractal structures that do not simply originate from a nonlinear turbulent cascade.

Fig. 7. Least-squares fit of the scaling exponents ξ(q) via the log-normal cascade model (5) (lines), for sample 10 (red circles) and for the average of all exponents (black stars). The fitting parameter μ is in agreement with that of classical hydrodynamic turbulence μ ≈ 0.02.

Finally, in order to include one example of standard data analysis technique for intermittent turbulence, we have estimated the kurtosis of the fluctuation distribution, namely the ratio between the fourth-order and the squared second-order structure functions $K({\mathcal{l}}_t)=S_4/S_2^2\sim {\mathcal{l}}_t^{-\kappa }$ (Frisch 1995; Dudok de Wit et al. 2013). The kurtosis provides information on the shape of the distribution of the scale-dependent fluctuations. At large scales (comparable with the system correlation scale) the Gaussian value K = 3 is typically observed. As the scale decreases, the inhomogeneous turbulent cascade generates intermittent structures, so that the distribution deviates from Gaussian, corresponding to increasing K. The scaling properties of turbulence result in the power-law scaling of K in the inertial range. The scaling exponent, κ, is a good measure of the efficiency of the cascade, namely how rapidly the small-scale structures are generated. This depends on the nature of the nonlinear interactions and can be used as a quantitative measure of intermittency (Castaing et al. 1990; Carbone & Sorriso-Valvo 2014). The kurtosis was estimated for all intervals, and a power-law fit was performed whenever a long enough scaling range was observed. An example of K with power-law fit is shown in Fig. 8. The resulting fitting parameters of κ are collected in Table A.3.

Fig. 8. Scale-dependent kurtosis K(ℓt) for samples 26 (group 1) and 27 (group 2). Power-law fits $K({\mathcal{l}}_t)\sim {\mathcal{l}}_t^{-\kappa }$ give the indicated exponents κ. The Gaussian value K = 3 is also shown (horizontal dotted line).

4. Observation of ion cyclotron waves

Observations of solar wind data often reveals the presence of wave activity near the end of the MHD inertial cascade range, and close to the kinetic plasma range. These are typically identified as kinetic Alfvén waves (KAW), or ion-cyclotron waves (ICW), among other modes (see e.g., Bale et al. 2005; Kiyani et al. 2012).

In order to explore the relationship between the observation of ion-scale waves and the properties of the inertial range turbulence, we introduced a quantity that enables the identification of waves and quantitatively assesses their presence in the time series. The technique for the identification of the waves is described in detail in Khotyaintsev et al. (2021), where it is used to determine that the observed fluctuations are most likely ion cyclotron waves. A brief description of the technique is given in the following. The first step is to rotate the magnetic field into the field-aligned coordinates using B, which is then low-pass filtered at 0.01 Hz as the background magnetic field. The power-spectrum of the resulting transverse component, δB_⊥, is shown in the third panels from the top of Fig. 2. Subsequently, the coherence between the two perpendicular magnetic field components is computed (fourth panel). This is expected to result in large values if circularly-polarized ion-scale waves are present. If the coherence is larger than an arbitrary threshold (0.65) in a frequency range near or below the proton gyro-frequency, the phase angle between the two perpendicular magnetic components is also computed (fifth panel), allowing to determine the fluctuations handedness; as we are interested in circularly-polarized waves, we exclude the phases outside the intervals +90° ±45° (right-handed waves) and −90° ±45° (left-handed waves).

Using the above indicators, it is therefore possible to unambiguously identify regions with circularly-polarized wave activity both in time and frequency. The perpendicular magnetic power is finally integrated in the identified wave intervals and frequency band (with lower and higher frequencies f₁ and f₂ respectively to be identified according to the above criteria), namely within the wave patches clearly visible in the phase angle plot. This procedure provides the time series of one frequency-integrated local parameter $Q_w={\int }_{f_2}^{f_1}\delta B_{\perp }^2\text{d}f$, defined for each data point in the time series, indicating the total power associated to the wave-like fluctuations. We note that the parameter is not computed outside of the wave patches. The bottom panels of Fig. 2 shows two examples of wave parameter Q_w for samples 26 and 27. Interval 26 (left panel) with highly irregular and intermittent behaviour, capturing the corresponding wave patches observed in the scalograms. In the adjacent interval 27 (right panel), waves are nearly absent and, accordingly, the wave parameter values are negligible.

Finally, using the time series of frequency-integrated wave power, two similar global parameters can be computed to quantitatively assess the occurrence of waves within each interval. The first one is simply obtained as the time-integrated power ${\overline{Q}}_{\mathit{w}}=\int Q_{\mathit{w}}\mathrm{d}t$, with the integration intended for each interval. The second one is the average over the interval $⟨Q_{\mathit{w}}⟩={\overline{Q}}_{\mathit{w}}/\mathrm{\Delta }\mathcal{W}$, taking into account the density of waves within each interval. The obtained values are listed in Table A.3 for all intervals. In some occasions, when no waves were identified, the parameters were set to 0. The wave parameters ${\overline{Q}}_{\mathit{w}}$ and ⟨Q_w⟩ are used in Sect. 5 to determine possible correlations with the turbulence parameters.

5. Discussion

Once the turbulent properties of the fluctuations and the presence of ion-scale waves have been quantitatively assessed, it is possible to investigate correlations between the two phenomena. This may help understanding the dynamical processes of solar wind plasmas and, in particular, the cross-scale coupling between fluid and sub-ion processes. Using the results of the analysis for the 36 samples, correlation coefficients have been computed between pairs of parameters of solar wind (V_sw and θ_vb), turbulence (Hurst number, spectral exponent, kurtosis and intermittency), and waves (the two estimators presented in Sect. 4). Both linear (Pearson) and nonlinear (Spearman) coefficients have been computed. For each pair of parameters, the largest of the two has been considered. The complete list of coefficients for all pairs is presented in Table A.4. As expected, some of the parameters are trivially correlated with each other, such for example those related to different scaling exponents of the same field, or the two wave parameters. Others display known correlations, associated with the nature of the solar wind intervals. However, despite the high variability of the parameters and the experimental conditions, some non-trivially related pairs show moderate, non-negligible correlation. These are highlighted in bold in Table A.2. The most interesting correlation was found between the wave indicators ($\overline{Q_{\mathit{w}}}$ and ⟨Q_w⟩) and the intermittency parameters (μ and κ). For example, C(⟨Q_w⟩,μ) = − 0.5 indicates that intervals with substantial presence of ion-scale waves are likely to show reduced intermittency. This observation clearly highlights the link between the characteristics of the fluid-scale turbulent cascade to the excitation of waves at ion scales.

A more visual description is provided by the scatter plots of pairs of parameters listed in Tables 1 and A.1, shown in Fig. 9 for some pairs of parameters from Table A.4. In all panels, the samples are color-coded according to their group as determined in Sect. 3.1 (group 1: blue circles; group 2: red squares; group 3: green triangles). Whenever relevant, vertical or horizontal lines indicate typical value of the parameters for standard fluid turbulence. The top-left panel, plotting the wind speed and the density Hurst number, highlights the clear separation between group 1 (mostly large Hurst number) and groups 2 and 3 (smaller Hurst number). Additionally, it clearly shows that while intervals of group 1 belong to both fast and slow wind, nearly all intervals of groups 2 and 3 (with one single exception) belong to faster solar wind. It is worth noting that the flow speed does not necessarily act as an ordering parameter for Alfvénicity, though slow wind is generally less Alfvénic than fast wind (see for example D’Amicis et al. 2021, and references therein). Specifically, fast wind can exhibit different levels of Alfvénicity. In this respect, the top-left panel suggests that the lack of Alfvénicity in the slow wind assures a more developed turbulence, while possible enhancements in Alfvénic nature of the fluctuations in some (though not all) fast wind samples may prevent plasma from fully developing into a turbulent state. It turns out that in the samples studied here the fast wind can include both standard and reduced turbulence intervals (e.g., with shallower spectra), depending on the corresponding level of Alfvénic fluctuations. We recall here that while solar wind turbulence is likely strongly driven by Alfvénic fluctuations, these need to include both counterpropagating modes in order to effectively generate nonlinear interactions. On the other hand, if the fluctuations are unbalanced, with one mode prevailing over the other (typically resulting in a definite sign large cross-helicity), then the sweeping effect strongly reduces the nonlinear interactions, resulting in weaker turbulence and shallower spectra (Dobrowolny et al. 1980).

Fig. 9. Scatter plots of pairs of solar wind, turbulence and wave parameters (see Tables 1 and A.1). Blue, red and green points indicate intervals of group 1, 2, and 3 respectively. Blue open circles in the bottom central and left panels represent intervals with wave zero parameters, for clarity, their value is set to 0.1 in order to be represented in the logarithmic vertical axis. Vertical and horizontal lines indicate standard fluid turbulence reference vales. The correlation coefficient is indicated (Table A.2).

The top-center panel of Fig. 9 shows that the angle between the magnetic field and the radial direction (approximately corresponding to the velocity vector and, in turn, to the sampling direction, at such distances) is also relevant to the turbulence. In particular, for group 1 intervals the Kolmogorov spectrum (dashed line) is observed at all angles. On the contrary, groups 2 and 3, characterized by a shallower spectrum, only include intervals with nearly radial field. This suggests that during intervals belonging to groups 2 and 3, Solar Orbiter sampled parallel fluctuations (namely the slab component of turbulence), which are generally less evolved and likely more Alfvénic with respect to 2D turbulence. Interestingly, this result is in contrast with critical balance theory (Goldreich & Sridhar 1995; Telloni et al. 2019a), which predicts a steeper spectrum (with a scaling close to −2) for parallel fluctuations (an interesting discussion on the validity and relevance of critical balance in solar wind turbulence is provided in Oughton & Matthaeus 2020).

The top-right panel of Fig. 9 shows that strong intermittency (large kurtosis) is mostly observed in intervals with quasi-perpendicular field. Additionally, it highlights the good correlation existing between the angle and the intermittency exponent κ, demonstrating that θ_vb is a good ordering parameter for intermittency. As mentioned above, in the studied intervals the solar wind plasma is likely to be more Alfvénic at quasi-parallel angles, where the turbulence is only poorly developed. The stochastic nature of the Alfvénic fluctuations tends to reduce the intermittency, which is indeed lower at larger angles. On the other hand, at quasi-perpendicular angles, where the turbulence is more fully developed (possibly in association with reduced Alfvénic fluctuations), the mitigating effect of Alfvénicity is lower and the coherent structures advected by the wind tend to emerge, resulting in the observed increase in intermittency.

The bottom-left panel of Fig. 9 shows the strong overall correlation between the spectral exponents of electron density and magnetic field magnitude. For the intervals of group 1, spectral exponents of both fields are mostly consistent with the standard Kolmogorov value. On the contrary, the evident linear correlation of the more variable exponents of groups 2 and 3 strongly suggests the Alfvénic nature of the fluctuations, with well-correlated compressive magnetic magnitude and plasma density fluctuations. It is indeed worth recalling that density behaves as a passive scalar (it reproduces the magnetic field magnitude characteristics) only in the Alfvénic solar wind, where the contribution of compressive fluctuations is negligible. In this perspective, for intervals in groups 2 and 3 the plasma density can be considered as a proxy of the magnetic field for the turbulent properties (spectral scaling, intermittency, etc.).

In the bottom-center panel of Fig. 9, the correlation between the angle and the wave density parameter is shown. In this case, no clear separation between the three groups is observed. However, it is evident that for intervals with perpendicular field the wave density is always small (noting that five intervals of group 1 for which ⟨Q_w⟩ = 0 have been artificially represented on the logarithmic vertical axis by the open blue circles at 0.1). The four intervals of group 3 are also characterized by large presence of waves. This is in good agreement with the expectations. Indeed, the presence of KAWs (at quasi-perpendicular angles) and ICWs (at quasi-parallel angles) strictly depends on the presence of Alfvénic fluctuations at fluid scales. Larger Alfvénicity is associated with enhanced presence of waves at ion scales, as first shown by Bruno & Telloni (2015; see Fig. 2 therein) on a single case study, and then corroborated on a statistical data set by Telloni et al. (2019b; see Fig. 3c therein).

Finally, the bottom-right panel of Fig. 9 shows the correlation between the intermittency scaling exponent κ and the total wave power ⟨Q_w⟩. For group 1 intervals with no waves, the same representation as in the bottom-center panel has been adopted. The observed correlation is a very interesting result. Indeed, despite the scattered plot, a general trend is evident: stronger intermittency intervals have less wave activity. This is in striking agreement with a scenario in which for higher the Alfvénic fluctuations (which implies a lower intermittency), the presence of waves at ion scales (and, in turn, the related measured energy) is larger. This has been very recently validated by the statistical work by Telloni et al. (2019b).

The overall conclusion gained from the examination of Fig. 9 is that 2D fluctuations (fluctuations sampled at quasi-perpendicular angles with the magnetic field) are always characterized by strong Kolmogorov turbulence, strong intermittency, and absence of wave activity. A less Alfvénic content is also suggested. As discussed above, all these fluid and kinetic characteristics are strictly related to each other. On the other hand, slab fluctuations (fluctuations sampled at quasi-radial directions) are associated with less developed turbulence (smaller spectral exponent and intermittency parameters) and stronger ion-scale wave activity. A higher Alfvénic content (which acts to make less efficient the nonlinear interactions) can be inferred for these intervals. In this respect, the speed of the solar wind flows does not seem to be an ordering parameter, whereas Alfvénicity is shown to be a more suitable one.

6. Conclusions

In this work, we present the first analysis of turbulence and intermittency of the solar wind electron density measured by the RPW instrument on board the Solar Orbiter spacecraft. First, 36 intervals were selected during the month of June 2020, when Solar Orbiter was located in the inner heliosphere, approximately at 0.5 AU from the Sun. The study was performed using standard analysis techniques as well as empirical mode decomposition (EMD). It was found that the intervals could be separated in three groups, according to their agreement with standard turbulence parameters. The distinction was attributed to the different level of Alfvénic fluctuations, which in the absence of plasma measurements was inferred from the turbulence characteristics. Using the magnetic field measurements from the MAG instrument, the presence of ion-scale waves was detected through coherence analysis. These waves have been identified as being mostly ion cyclotron waves. A new parameter was introduced to determine quantitatively the energy associated to waves in a given frequency range. Making use of such parameter, two estimators were introduced to assess the overall wave activity within each interval. Comparing these wave parameters with the turbulence indicators, we found that intervals with enhanced presence of waves are also characterized by anomalous turbulence and weaker intermittency of the solar wind density. While the study of statistical correlations cannot determine the causality relation between the phenomena, it definitely provides constraints and, in this specific case, they are shown to help improve the understanding of the cross-scale connection between the fluid scales and the ion-scales. Here, we highlight the possible role of Alfvénic fluctuations in driving both the observed ion cyclotron waves and the reduced turbulence and intermittency. We note that the presence of ion cyclotron waves may play an important role in the excitation of electrostatic fluctuations at higher frequency (Valentini & Veltri 2009; Valentini et al. 2014).

The preliminary results described in this paper demonstrate the unprecedented high quality of the Solar Orbiter RPW electron density data estimated using the probe-to-spacecraft potential. The excellent performances of the EMD-based analysis allowed the accurate determination of turbulence parameters. The results described here represent the most detailed description of turbulence of solar wind density fluctuations thus far. While these results are generally in line with previous observations at 1 AU, future studies of Solar Orbiter measurements will finally allow us to determine the radial evolution of the properties of density turbulence.

Furthermore, the novel parameter providing quantitative assessment of the presence of waves helped identify an important relationship between fluid-scale turbulence and ion-scale phenomena in the solar wind plasma. The future study of measurements including magnetic field and velocity fluctuations, and more extended statistical analysis could facilitate a deeper understanding of such relationship.

Acknowledgments

Solar Orbiter is a space mission of international collaboration between ESA and NASA, operated by ESA. We thank the entire Solar Orbiter team and instrument PIs for data access and support. Solar Orbiter data are available at http://soar.esac.esa.int/soar. The RPW instrument has been designed and funded by CNES, CNRS, the Paris Observatory, The Swedish National Space Agency, ESA-PRODEX and all the participating institutes. Swedish co-authors are supported by the Swedish Research Council, grant 2016-05507, and Swedish National Space Agency grant 20/136. CNES and CDPP are acknowledged for the support to the French co-authors. Solar Orbiter magnetometer operations are funded by the UK Space Agency (grant ST/T001062/1). T.S.H. is supported by STFC grant ST/S000364/1. F.C. acknowledges the contribution received from EU-H2020 program ERA-PLANET through the project “Integrated Global Observing Systems for Persistent Pollutants” (iGOSP) (Grant Agreement: 689443), funded under H2020-SC5-15-2015 “Strengthening the European Research Area in the domain of Earth Observation”, from FET Proactive project “Towards new frontiers for distributed environmental monitoring based on an ecosystem of plant seed-like soft robots” (I-Seed), funded under Horizon 2020 research and innovation programme (Grant agreement: 101017940), and from EU-H2020 project “EuroGEO Showcases: Applications Powered by Europe” (e-shape) (Grant Agreement: 820852), funded under H2020-SC5-2018-2 “Strengthening the benefits for Europe of the Global Earth Observation System of Systems (GEOSS) – establishing EuroGEO”. L.S.V. was funded by the Swedish Contingency Agency grant 2016-2102 and by SNSA grant 86/20. CLV was partially supported by EPN project PIM-19-01.

References

Appendix A: Additional tables

All Tables

All Figures

	Fig. 1. Overview of the data. From top to bottom: solar wind electron density ne (RPW), deHoffmann-Teller solar wind velocity estimate Vsw (RPW) and interplanetary magnetic field components, Bi, and magnitude, |B| (MAG), measured by Solar Orbiter during the whole month of June 2020.
In the text

Fig. 2. Selected intervals. Two examples of adjacent intervals during day 20 June 2020. Left: interval 26; right panel: interval 27. From top to bottom: interplanetary magnetic field components, Bi, and magnitude, |B| (MAG); electron density ne (RPW); magnetic filed spectrogram; perpendicular magnetic field components coherence; perpendicular magnetic field components phase angle; wave parameter Qw(t) (see Sect. 4).

In the text

Fig. 3. Results of the EMD analysis. Left panel: example of IMFs ϕj(t) for sample 1 (black lines). The bottom plot (red line) represents the residual rn(t). For better readability, only odd IMFs have been plotted. Central panel: average timescale ⟨τj⟩ of each IMF of sample 1 as a function of the mode j. Error bars represent the 95% confidence bounds. The dashed line is a least square fit obtained from the relation ⟨τj⟩=αγj. Right panel: IMF variance Var(ϕj) as a function of the average period ⟨τj⟩ for sample 1, with the associated 95% confidence bounds. The dashed line represents the relation Var(ϕj)∼⟨τj⟩2ℋ, being ℋ the Hurst number (see the list for all intervals in Table A.1).

In the text

Fig. 4. Results of the Hilbert-Huang analysis. Left panel: second-order Hilbert spectrum ℒ2(f) (blue circles), the classical Fourier PSD (Welch 1967) E(f) (red line), and the second-order structure function S2(ℓt) (green squares, plotted as a function of the inverse timescale ℓt), for sample 1. Power scaling is present in the same frequency range for all methods. The dashed line represents the classical Kolmogorov scaling ℒ2(f)∼f−5/3, while the dotted line represent the expected scaling for the second-order structure function $S_2({\mathcal{l}}_t)\sim {\mathcal{l}}_t^{\zeta (2)}$ (ζ(2) = β2 − 1). Central panel: generalized Hilbert spectra ℒq(f) for q ∈ [1, 5], obtained for sample 1. The curves have been vertically shifted for clarity. The shaded area represents the frequency range of the bootstrapping least-square fit. Right panel: same as in the central panel but for sample 9. The power behavior is still present, but the power law exponents are considerably different.

In the text

	Fig. 5. Probability distribution function P(β2) of the scaling exponents β2, constructed via bootstrap resampling, for sample 1 (upper panel) and sample 9 (lower panel). In both panels, the vertical dotted line represents the median of the distribution (50th percentile), while the vertical dashed bars indicate the 95% confidence interval (enclosed between 2.5th and 97.5th percentile).
In the text

Fig. 6. Results of the equivalent structure function analysis. Top-left panel: HSA equivalent scaling exponents ξ(q) for group 1 (blue circles) and their average (black crosses). The exponents ζ(q) for velocity fluctuations in the inertial range of hydrodynamic turbulence (red squares) are shown for reference (Benzi et al. 1993). The dashed line represents the theoretical expectation ζ(q) = q/3 (Kolmogorov 1941). Top-center and top-right panels: equivalent scaling exponents for the magnetic field radial component Br and magnitude |B|, for group 1. Bottom panels: same as the top panels, but for the intervals of group 2.

In the text

	Fig. 7. Least-squares fit of the scaling exponents ξ(q) via the log-normal cascade model (5) (lines), for sample 10 (red circles) and for the average of all exponents (black stars). The fitting parameter μ is in agreement with that of classical hydrodynamic turbulence μ ≈ 0.02.
In the text

	Fig. 8. Scale-dependent kurtosis K(ℓt) for samples 26 (group 1) and 27 (group 2). Power-law fits $K({\mathcal{l}}_t)\sim {\mathcal{l}}_t^{-\kappa }$ give the indicated exponents κ. The Gaussian value K = 3 is also shown (horizontal dotted line).
In the text

Fig. 9. Scatter plots of pairs of solar wind, turbulence and wave parameters (see Tables 1 and A.1). Blue, red and green points indicate intervals of group 1, 2, and 3 respectively. Blue open circles in the bottom central and left panels represent intervals with wave zero parameters, for clarity, their value is set to 0.1 in order to be represented in the logarithmic vertical axis. Vertical and horizontal lines indicate standard fluid turbulence reference vales. The correlation coefficient is indicated (Table A.2).

In the text