Statistical Analysis of Field-Aligned Alfvénic Turbulence and Intermittency in Fast Solar Wind

The statistical properties of fast Alfvénic solar wind turbulence have been analyzed by means of empirical mode decomposition and the associated Hilbert spectral analysis. The stringent criteria employed for the data selection in the Wind spacecraft database, has made possible to sample multiple k‖ field-aligned intervals of the three magnetic field components. The results suggest that the spectral anisotropy predicted by the critical balance theory is not observed in the selected database, whereas a Kolmogorov-like scaling (E(k‖) ∼ k−5/3) and a weak or absent level of intermittency are robust characteristics of the Alfvénic slab component of solar wind turbulence.


Introduction
The solar wind (SW) is the most studied environment for plasma turbulence [1]. Among the most interesting features of space plasmas, the existence of a strong ambient magnetic field leads to turbulence anisotropy [2] related to the preferential energy transfer in a perpendicular direction with respect to the mean magnetic field, as theoretically predicted [3][4][5]. Observations show different forms of such anisotropy. For example, elongation or deformation of eddies in the mean magnetic field direction, called wavevector anisotropy, produces k ⊥ > k ). On the other hand, different power levels in different directions relative to the local field, called power anisotropy, result in power spectra of different amplitude P(k ⊥ ) = P(k ). Finally, different scaling in different directions, measured through the power spectrum slope β ⊥ = β (where β represents the power-law exponent), is usually called spectral index anisotropy [1,6].
In the phenomenology of Alfvénic turbulence, the ratio between the linear Alfvén timescale and the nonlinear decay timescale regulates the balance of the energy flux across scales, thus determining the strength of turbulence. The linear Alfvén timescale, defined as τ A /υ A , represents the typical propagation time of an Alfvén wave along a structure of size , υ A being the Alfvén speed. The nonlinear decay time τ nl ⊥ /δu( ⊥ ) represents the perpendicular eddy turnover time, where δu( ⊥ ) is the fluctuation velocity across a structure of size ⊥ . When ratio of these two quantities χ = τ A /τ nl δu/ ⊥ υ A 1, oppositely propagating Alfvén wave packets may interact multiple times before their energy is transferred to smaller scales [7,8]. On the contrary, for χ 1 Alfvénic fluctuations undergo multiple decays before they interact. Both conditions naturally evolve towards comparable nonlinear and Alfvén timescales (χ 1), so that the nonlinear cascade occurs within one single interaction [3,9] (the reader is referred to Refs. [4,10] for more details on the MHD turbulence theory and on the evolution of the nonlinearity parameter χ).
In this case, it is possible to show that the scaling relation between perpendicular and parallel wave vectors is k ∝ k 2/3 ⊥ . Under this particular condition, usually called "critically balanced" state [4], wavevector anisotropy k ⊥ k would appear at small scales, as resulting from the anisotropy of the spectral index: E(k ⊥ ) ∝ k −5/3 ⊥ and E(k ) ∝ k −2 . Therefore, the critical balance condition can be tested by measuring the spectral index at various angles from SW measurements, a complicated issue using typical single-spacecraft SW measurements. Wavelet analysis has been often used to study the local mean magnetic field and to collect information on fluctuations at a particular local angle θ VB between the local mean magnetic field and the flow direction. The resulting spectral slope was close to −5/3 when θ VB ∼ 90 • and −2 when the angle is θ VB ∼ 0 • , providing an evidence for the critical balance theory [11][12][13][14]. Multi-point measurement techniques provided similar results [6,15]. The observed spectral anisotropy is also compatible with the presence of discontinuities in the SW, as for example current sheets [2,16]. Removal of intermittent events from data was accompanied by flatter parallel spectra [17].
However, both local mean-field calculation and removal of intermittent events might result in systematic errors in the reconstruction of the parallel spectrum. The observation of local alignment between magnetic field and wind direction is indeed prevented by pervasive large-amplitude Alfvénic fluctuations. Thus, very long time series are required to build the parallel spectrum [13] using randomly distributed data points, possibly uncorrelated, and artificially introducing the discontinuities that might generate the spectral anisotropy.
In order to study the SW turbulence in a direction parallel to the mean background magnetic field, a novel approach was recently proposed. A search of continuous Alfvénic, fast SW intervals was performed, producing a database mostly populated by k fluctuations. This allows to avoid the possible introduction of discontinuities and inhomogeneities. The selected intervals were studied using the Hilbert spectral analysis, so to remove the possible influence of small-scale noise and large-scale structures from spectral and intermittency properties in the inertial range. Preliminary results showed that spectral anisotropy is incompatible with the critical balance scenario. This work represents a follow-up of the results presented in [18]. Here, an in-depth statistical study based on the Hilbert spectral analysis and a different resampling method have been conducted on a subset of the intervals presented in the previous work, showing the robustness of the observations. In particular, despite of the power anisotropy, present in all samples, the distribution of the spectral exponents obtained from the resampling (and their relative scaling exponents) is characterized by a mean value fully compatible with the classical Kolmogorov spectrum (within the 95% of prediction interval), revealing the presence of short-range correlations.

Data Selection Criteria for k Magnetic Fluctuation
Alfvénic turbulence can be usually split in two components: a slab component, for fluctuations propagating along to background magnetic field B 0 , and a 2D component, whose magnetic fluctuations lie on a plane perpendicular to B 0 [19][20][21][22]. Therefore, wavevectors of the slab and 2D components are respectively parallel (k ) and perpendicular (k ⊥ ) to the global mean magnetic field.
Testing the critical balance model in SW fluctuations requires a reliable estimate of the spectral exponents of the slab (k ). However, sampling the SW in a purely parallel direction is rather challenging, in particular when Alfvénic turbulence is dominating. Indeed, even when the global mean background magnetic field is aligned to the sampling direction (which coincides with that of the bulk flow), the presence of high amplitude Alfvénic fluctuations can locally disrupt the alignment. As a result, the measured Alfvénic turbulence includes contributions from fluctuations along both k ⊥ and k directions.
In this work, a systematic search of SW intervals that provide k sampling has been performed using over 12 years of SW magnetic field observations, recorded by the Wind spacecraft between 2005 and 2016, with sampling frequency of ∆ω ≈ 11 Hz. In order to be selected as a k sample, a SW interval must satisfy the following requirements: 1. Magnetic field and velocity must be aligned within 15 • at all scales in the turbulent inertial range for the whole interval, i.e., both local and global mean magnetic field must be aligned with the bulk velocity; 2. The total duration cannot be less than 1 hour, with maximal percentage of missing data ≤20%; 3. The whole interval must be enclosed within a fast SW stream (V ≥ 550 km s −1 ); 4. It must have low magnetic compressibility (C B ≤ 25%), measured in the inertial range at scale of 20 min, where C B is defined as the ratio between the variance of the magnetic field intensity fluctuations and the total variance of the fluctuations, C B = σ 2 |B| / ∑ i=x,y,z σ 2 B i [23]; 5. The total pressure P m = 2nk B T + B 2 /(8π) (with n, T, and B are, respectively, the proton number density, the temperature, and magnetic field intensity) must satisfy the condition P m < 0.05 nPa, in order to avoid discontinuities or shocks [24].
If all of these conditions are fulfilled, the selected sample guarantees k fluctuations, due to undisturbed Alfvénic fast SW turbulence. Indeed, local sampling of parallel fluctuations requires angles between B and V vectors The more stringent criterion adopted here ensures thus local sampling of purely k fluctuations.
All selected samples are immersed in trailing edges of high-speed streams, in the strong turbulence regime. Some relevant parameters relative to the selected intervals are collected in Table 1: SW speed V , total pressure P m , magnetic compressibility evaluated at the fluid scale 20 min C B , and the mean angle θ VB between the velocity and magnetic field vectors. Table 1. Parameters of wind data considered for the analysis: the SW bulk speed, proton number density and temperature were measured by the Solar Wind Experiment (SWE) instrument [25] at ∆ω ≈ 0.01 Hz resolution; magnetic field at resolution ∆ω ≈ 11 Hz was measured by the Magnetic Field Investigation (MFI) magnetometer [26]. The magnetic compressibility C B has been evaluated at scale of 20 min, within in the inertial range.

Date
Start As customary, the link between the frequency space and the wavevector space is obtained by using the Taylor frozen hypothesis [1,27], which holds for all the intervals in this work. In It should be pointed out that in the selected intervals foreshock particles from the Earth's bow shock are unlikely but cannot be excluded [28,29]. However their possible role in the selected intervals is beyond the scope of this analysis, and it will be left out for a future work.

Empirical Mode Decomposition and Arbitrary Order Hilbert Spectra
In the Empirical Mode Decomposition (EMD) framework [30,31], the three components of the magnetic field B ≡ B x (t), B y (t), B z (t) can be decomposed each into a finite number n of oscillating basis functions, of increasing characteristic time scale τ, known as intrinsic mode functions (IMFs) (1) describes the mean trend. EMD is based on the local properties of the data [30], allowing to characterize non-stationary and nonlinear processes, and does not necessarily require a continuous time series, since it can be also applied to time series with non uniformly sampled data [32,33]. The decomposition process, called sifting process, consists principally of two stages: the local extrema of B i are identified and subsequently interpolated through cubic spline. Once interpolated, the envelopes of local maxima and local minima are defined; subsequently, the mean of the two envelope functions is subtracted from the original data. This difference is considered an IMF only if: (i) the number of local extrema and zero crossings does not differ by more than 1; (ii) at any point t i , the mean value of the extrema envelopes is zero. If the difference does not meet these criteria, the sifting procedure is repeated [30]. A general rule to stop the sifting is defined by a Cauchy-type convergence criterion, defining a standard deviation σ, evaluated from two consecutive steps of the algorithm. The iterative sifting stops when σ is smaller than a fixed threshold, usually σ thresh ≈ 0.2. An alternative method for stopping the sifting procedure is the so called 3-thresholds stoppage criterion, and was introduced in order to guarantee globally small fluctuations in the mean, while taking into account locally large amplitude variations [31]. In this work this second stoppage criterion has been adopted. The three threshold parameters used here are δ = 0.05, ξ 1 = 0.05, and ξ 2 = 10ξ 1 , in accordance with their typical values [31]. An example of the IMFs φ j (t) for the component B 1 ≡ B x (t) and the associated residual R φ (t) extracted from sample 2005/01/03, are shown in the left panel of Figure 2. Note that mode j = 1 (not shown) contains the fastest timescale, in many cases associated with the experimental noise embedded in the data-set. Once the IMFs have been obtained, an Hilbert transform is applied on each mode: where P represents the Cauchy principle value. This transformation allows to extract a time-dependent instantaneous frequency ω i,j (t) and a time-dependent amplitude modulation A i,j (t), from the analytical signal Z = X i,j (t) + iX i,j (t) ≡ A i,j (t)e iΦ i,j (t) [34]. Central panel of Figure 2 shows the characteristic periods τ 2,j (component B 2 (t), in log-linear plot), evaluated as the average inverse instantaneous frequency of the j-th mode: τ 2,j = ω 2,j −1 t ( · t representing a temporal average) [30], where ω i,j = 2π −1 dΦ i,j (t)/dt is the temporal variation of the instantaneous phase Φ i,j (t) [30]. Moreover, the Hilbert spectrum, defined as H(ω, t) = A 2 (ω, t), provides energy information in the time-frequency domain [35].
When the EMD is applied on turbulent processes [36][37][38] or multifractal processes [39,40], it intrinsically acts as a dyadic filter bank [41][42][43][44], where each IMF captures a narrow band in frequency space, and their characteristic period τ j follows an exponential law of the form: τ j = α × γ j , where γ = 2 for an exact dyadic decomposition. The dashed line in the central panel of Figure 2 represents the least square fit of the exponential relation obtained for τ j , with γ = 1.97 ± 0.01, showing very good agreement with the theoretical expectation. Moreover, the dyadic structure can be observed by analyzing the Fourier power spectral density (PSD) of the various IMFs, reported in Figure 2, right panel. The PSD of the data (black line) presents a power-law decay well reproduced by the superposition of the various IMFs Fourier spectra. The superposition behaves as a power-law of the form [35,45,46], with the same spectral exponent as the original PSD (Figure 2 right panel). In the classical, self-similar Kolmogorov theory [47,48], the spectral exponent β is linked to the scaling exponent of the second-order structure function (SF) via the relation A marginal integration performed on H(ω, t) provides the Hilbert marginal spectrum h(ω) = T −1 T 0 H(ω, t)dt, defined as the energy density at frequency ω [30,49]. Moreover, from the Hilbert spectrum, a joint probability density function P(ω, A) can be extracted, using the instantaneous frequency ω i,j and the amplitude A i,j of the j-th IMF as: h(ω) = ∞ to a second-order statistical moment [35]. This procedure can be generalized to the arbitrary order q ≥ 0 by defining the instantaneous ω-dependent q-th-order statistical moments: The second-order moment h(ω) ≡ L 2 is the analogous of the Fourier spectral energy density, and can be interpreted as the energy associated to the frequency ω. It should be pointed out that the definition of frequency in h(ω) is different from the definition in the Fourier framework. The interpretation of the Hilbert marginal spectrum should thus be given with more caution [30,35,50]. Figure 3 (left panel) shows an example of the second order moment L 2 (ω) obtained for the different components B i of the magnetic field of sample 2005/01/03. From the analysis, it is evident that the B x component is characterized by a lower power with respect to the two perpendicular components B y and B z . In Figure 3, the central panel shows the average isotropy ratio R(ω ) = L x 2 (ω )/L α 2 (ω ) (where α = z, y), obtained by averaging on all frequencies enclosed in the inertial range ω ∈ [10 −2 , 10 −1 ]. For comparison, the theoretical value for homogeneous and isotropic turbulence R(ω ) = 3/4, obtained from the Kolmogorov's second similarity hypothesis [51], is also indicated. in all cases, R(ω ) is always different from 3/4, indicating that the effects of the power anisotropy are always present despite of the short length of the samples. In addition, the the ratio x/z is always comparable with the ratio x/y, within the error. The PSD for sample 2005/01/03 (B x component) is plotted in Figure 2, right panel. The slope β 2 of the Hilbert spectra L 2 (ω) is fully compatible with the Fourier PSD, and clear power-law scaling. The right panel of Figure 3 displays a comparison between the Fourier and Hilbert spectra with the standard second-order structure function, S 2 ( t ), plotted against the inverse scale −1 t in order to match the frequency range. From the results, the agreement between the slope β 2 ≈ 1.61 ± 0.10 and the exponent ζ(2) = 0.57 ± 0.02 is evident, and both are in good agreement with the theoretical result ζ(2) ≡ β − 1 obtained from the classical Fourier PSD [48]. Similar values were found for all intervals selected for the analysis.

Statistical Analysis of HSA Slopes and Scaling Exponents for Alfvénic Turbulence
In order to construct a prediction interval for the various β q , a bootstrapping procedure has been adopted [52]. The bootstrap method is a general resampling procedure useful for the estimation of the distributions of statistics based on independent observations. The basic idea is that inference about a population can be modelled by resampling the data and performing inference from resampled data [52,53]. This allows to create multiple smaller data-sets (of equal size), by randomly extracting elements from the original data-set, calculating their statistic, and taking the average of the calculated statistics [54].
In order obtain a statistical distribution of the spectral slopes, two possibles approaches can be used. The first approach, known as residual resampling, consists of performing a least square fit on the original data, and the residuals are then resampled and added to the prediction of the original fit in order to obtain a new data-set, which will be fitted as a function the original frequency ω [52,55]. In the second approach, the original data L q (ω) and the corresponding frequencies ω are both resampled in order to construct the smaller subsets, and then least square fitted [56,57]. All those points have been selected in a range of frequency ω slightly larger than 1 decade in the inertial range. Both methods provided compatible results.
As a rule of thumb, for a simple statistical test the number N samp of random samples is usually chosen in the interval N samp ∈ [50, 100]. Indeed, while N samp = 50 provides a relatively good standard error estimate [55], for N samp > 100 the improvement in the estimation of standard errors becomes negligible [58]. On the other hand, when evaluating confidence or prediction interval a higher number of replicas are required. Here values N samp ∈ [600, 1000] have been selected, in accordance with recommended values [59].
The resampling has been performed on all intervals, and for all orders q. The PDFs of slopes , constructed for different replica number N samp , are shown in Figure 4. It is evident that all the curves collapse on the same distribution, indicating good statistical convergence. Moreover, the distributions exhibits a peak at a value β 2 ≈ 5/3, consistently with previous results from six-minute intervals, with similar physical properties, reported in the literature [60].   Figure 5 presents the PDFs of the slopes β q obtained from the Hilbert spectra L q (ω) for the three B i components. The maximum order q max = 3 was selected using the approximate rule q max ≈ log N − 1. For each component B i , distributions have been built by bootstrapping the ensemble of all the intervals listed in Table 1. The bin width ∆ was selected in accordance with the classical Freedman-Diaconis rule, ∆ = 2IQR[β q ]n −1/3 , where IQR is the interquartile range of β q , and n the number of samples.
Since no a priori assumption has been made on the distribution, this rule represents a suitable choice, as it is able to evidence the possible presence of heavy tails in the data. From the results of the second order Hilbert spectra, a slope β 2 ≈ 5/3 has been found, indicating that the parallel magnetic turbulence in fast, Alfvénic SW is characterized by the typical Kolmogorov spectrum. This was observed for all components B i [18,61]. Table 2 collects the slopes corresponding to various percentiles of the distribution of Figure 5, p 2.5 , p 50 , and p 97.5 , used to represent the 95% of prediction interval. Values are given for the three orders q of the various components B i . As discussed before, a relation between the PSD index and the scaling exponent of the second order SF ζ(2) exists. Such relationship can be extended to any arbitrary order q, allowing to define a family of generalized scaling exponents ξ(q) through the generalized Hilbert spectra L q [36,37,46] as ξ(q) ≡ β q − 1. Thus, the set of exponents ξ(q) is the Hilbert analogue of the standard set of scaling exponents of the structure functions or, if necessary, of their Extended Self-Similarity equivalent [62,63]. Table 2. Various percentile obtained from PDFs of Figure 5, for the three order Hilbert spectra L q : p 2.5 , p 50 , and p 97.5 . Relation 3, therefore, represents an alternative method of extracting the spectral slope β and the SF scaling exponents, with the advantage of constraining the effects of noise and large-scale structures. The dependence of the scaling exponents on the order, and in particular its deviation from a linear relation, are useful for a quantitative estimate of the intermittency level [35,45,48].

Component Order
The relation between ζ(q) and the Hurst exponent H is well known. In the classical Kolmogorov theory, the exponents are related to the qth-order moment of the fluctuations via the relation ζ(q) = qH, where the Hurst exponent describes the dynamical properties of the process: persistence (correlated fluctuations, long-term memory), or anti-persistence (anticorrelated fluctuations, short-term memory).
Examples of PDFs of ξ(q), with the values associated to percentiles p 2.5 , p 50 , and p 97.5 indicated by vertical bars, are shown in Figure 6 for the sample dated 2005/01/03, as obtained from the generalized Hilbert spectra. Finally, Figure 7 shows the order q dependence of the scaling exponents ξ(q), estimated through the HSA in the inertial range, for the three magnetic components (color coded) and for three of the samples used here (the three different panels). Markers and error bars represent the median and the uncertainty of the ensemble, obtained as described above as the p 2.5 , p 50 , and p 97.5 percentiles of the distributions. The dashed lines represent linear dependency of the exponents on the order, expected for homogeneous, non-intermittent scaling of the fluctuations. The Hurst exponents H are thus computed through a simple least-square fit of the linear relation, and the resulting values are included in each panel. From the analysis of the scaling exponents ξ(q), three different scenarios have emerged: (I) all the components of B(t) field have the same Hurst number H ( Figure 7); (II) the three H are substantially different for the components of B(t); (III) two components are similar and the third is different. Despite the small differences observed in the various B i components, the Hurst exponent never reach the value for uncorrelated processes. For all samples the Hurst exponent present a value H < 0.5 (Figure 7), revealing (consistently) the presence of an anti-persistent dynamic (short-range correlations), compatible with the Kolmogorov scaling observed in the tradiational fluid turbulence. In some cases, a small deviation from the pure linear scaling is observed as the moment q increase, signature of the presence of a weak intermittentcy level. The Kolmogorov-like scaling k −5/3 (nearly self-similar) and a weak intermittency, could be characteristics features of the slab component (k magnetic fluctuations) of the Alfvénic SW turbulence. Relatively weak levels of intremittency, or even its complete absence, have been observed in another experimental study, performed in the inverse cascade range of the two-dimensional turbulence, for an electromagnetically driven flow [64]. However, it should be pointed out that, contrary to what shown in Figure 3, in the 2D experiment no significant difference was found between parallel and perpendicular fluctuations [64]. Similar results have been reported also for others SW intervals [13,17] (in accordance with the selection criteria presented in this work). In particular it has been found that by removing the intermittency from the data, the PSD slope becomes close to k −5/3 losing the angle dependence between the local mean magnetic field and the flow direction.

Conclusions
In the attempt to describe the statistical properties of the k magnetic fluctuations of the fast Alfvénic solar wind, the EMD and the HSA have been applied to accurately selected Wind data intervals. By evaluating the generalized q-th order Hilbert spectra L q (ω), the analogous of the scaling exponents of the structure functions have been evaluated from the data. In order to obtain a robust indication of the scaling properties of the various samples, a bootstrap resampling procedure have been applied in order to construct a prediction interval for the spectral slopes β q , and scaling exponents ξ(q), measured via the Hilbert spectra L q , in the fast Alfvénic fast solar wind data. From the PDFs for the various slopes β 2 , estimated through the second order HSA analogous of the Fourier PSD, it has been found that β x 2 ≈ 1.65, β y 2 ≈ 1.70, and β z 2 ≈ 1.68, which robustly indicate that the parallel magnetic spectrum in fast, Alfvénic solar wind are characterized by the classical Kolmogorov β = 5/3 spectral index. Moreover, the scaling exponents ξ(q) shows an almost linear dependence on the order q, and despite the small differences among the three components B i , the Hurst number H ≡ ξ(1) never reaches the threshold value for uncorrelated processes. In other words, H < 0.5, revealing the anti-persistent behavior (short range correlations) of the k magnetic fluctuations. The relatively weak level of intermittency (or its absence) could also represents a key characteristics of k magnetic fluctuations. In light of this, since intermittency is a fundamental for the nonlinear turbulent cascade, it can be argued that the k turbulent properties could arise from a superposition of uncorrelated Alfvénic fluctuations, rather than from a fully developed nonlinear cascade. In this framework, the turbulent cascade could be activated among oppositely directed Alfvén modes, indicating that quasi uni-directional propagating Alfvén waves may produce a Kolmogorov-like turbulent spectrum [22].
The results presented in this work will be relevant to the Parker Solar Probe and Solar Orbiter missions, whose orbital characteristics are likely to sample solar wind intervals with the required characteristics to better confirm our observations and additionally study the radial evolution of such turbulence. Funding: F.C. was partially supported by the ERA-PLANET programme (www.era-planet.eu) (Contract.no. 689443) within the IGOSP project (www.igosp.eu).

Conflicts of Interest:
The authors declare no conflict of interest.