Multidimensional Iterative Filtering: a new approach for investigating plasma turbulence in numerical simulations

Turbulent space and astrophysical plasmas exhibit a complex dynamics, which involves nonlinear coupling across different temporal and spatial scales. There is growing evidence that impulsive events, such as magnetic reconnection instabilities, lead to a spatially localized enhancement of energy dissipation, thus speeding up the energy transfer at small scales. Capturing such a diverse dynamics is challenging. Here, we employ the Multidimensional Iterative Filtering (MIF) method, a novel technique for the analysis of nonstationary multidimensional signals. Unlike other traditional methods (e.g., based on Fourier or wavelet decomposition), MIF does not require any previous assumption on the functional form of the signal to be identified. Using MIF, we carry out a multiscale analysis of Hall-magnetohydrodynamic (HMHD) and hybrid particle-in-cell (HPIC) numerical simulations of decaying plasma turbulence. The results assess the ability of MIF to spatially identify and separate the different scales (the MHD inertial range, the sub-ion kinetic, and the dissipation scales) of the plasma dynamics. Furthermore, MIF decomposition allows to detect localized current structures and to characterize their contribution to the statistical and spectral properties of turbulence. Overall, MIF arises as a very promising technique for the study of turbulent plasma environments.

(Received xx; revised xx; accepted xx) Turbulent space and astrophysical plasmas exhibit a complex dynamics, which involves nonlinear coupling across different temporal and spatial scales. There is growing evidence that impulsive events, such as magnetic reconnection instabilities, lead to a spatially localized enhancement of energy dissipation, thus speeding up the energy transfer at small scales. Capturing such a diverse dynamics is challenging. Here, we employ the Multidimensional Iterative Filtering (MIF) method, a novel technique for the analysis of nonstationary multidimensional signals. Unlike other traditional methods (e.g., based on Fourier or wavelet decomposition), MIF does not require any previous assumption on the functional form of the signal to be identified. Using MIF, we carry out a multiscale analysis of Hall-magnetohydrodynamic (HMHD) and hybrid particle-in-cell (HPIC) numerical simulations of decaying plasma turbulence. The results assess the ability of MIF to spatially identify and separate the different scales (the MHD inertial range, the sub-ion kinetic, and the dissipation scales) of the plasma dynamics. Furthermore, MIF decomposition allows to detect localized current structures and to characterize their contribution to the statistical and spectral properties of turbulence. Overall, MIF arises as a very promising technique for the study of turbulent plasma environments.
Space and astrophysical plasmas are often found in a turbulent state, characterized by a disordered and chaotic dynamics encompassing many different spatiotemporal scales. A key aspect of turbulence studies concerns unraveling the physical mechanisms responsible for the transfer and dissipation of energy across such scales. In-situ spacecraft observations of plasma turbulence in the solar wind (Bruno et al. 2009;Chen et al. 2013;Kiyani et al. 2015) and in the Earths magnetosheath (Stawarz et al. 2016;Chen & Boldyrev 2017), show that power spectra of magnetic fluctuations exhibit power-law behaviors encompassing several orders of magnitude in frequency. At large scales, where the plasma can be described as a fluid within the framework of Magnetohydrodynamics (MHD), magnetic spectra follow a Kolmogorov-like power-law, which denotes the existence of an inertial range where the scale-to-scale energy transfer takes place, without losses, via interactions between the turbulent eddies (Iroshnikov 1963;Kraichnan 1965). As the ion characteristic scales (i.e., the ion inertial length d i and/or the ion gyroradius) are reached, multifluid effects (such as, e.g., the appearance of Hall currents) and ion-kinetic effects become important. Below such scales, we observe a transition to a magnetic power spectrum with a steeper slope (with values ranging from -2 to -4, but tipically around -3), until dissipation scales are reached (for a comprehensive review, see Verscharen et al. 2019, and references therein).
Due to the complexity of the kinetic plasma dynamics, we still lack a definite explanation for the existence of a turbulent cascade beyond the inertial range. Several attempts invoke nonlinear wavelike interactions of dispersive modes (described in terms of, e.g., kinetic alfvén and/or whistler waves, Howes et al. 2008;Schekochihin et al. 2009), eventually complemented by kinetic dissipative effects, such as Landau damping (Sulem et al. 2016) and other effects associated to the deformation of the particles velocity distribution functions (Del Sarto & Pegoraro 2018;Yang et al. 2017). From the other side, numerical evidence hints that the MHD inertial range extends beyond ion scales, provided that one includes the Hall term in the generalized Ohm's law (Hellinger et al. 2018), so that the steeper slope is not caused by any dissipative process.
This scenario gets further complicated by the existence of coherent structures, such as current sheets, which naturally form in turbulent environments. These are related to a phenomenon known as intermittency (Frisch 1995;Marsch & Tu 1997), that is, the occurrence of sudden changes in the magnetic fluctuations which lead to a spatially inhomogeneous energy cascade and dissipation. Indeed, there is growing evidence that plasma instabilities, such as magnetic reconnection triggered in spatially localized current sheets, enhance magnetic energy dissipation (Camporeale et al. 2018). This casts some doubts on the interpretation of the turbulent cascade in terms of wavelike modes only, and more in general on models that do not include intermittent effects from coherent structures. In this context, theoretical models introducing the concept of reconnection-mediated turbulence, have been proposed (among others, Boldyrev & Perez 2012;Loureiro & Boldyrev 2017;Mallet et al. 2017;Landi et al. 2019).
Nowadays, investigating plasma turbulence using direct approaches is becoming more and more feasible. The increasing computational capabilities allow to run direct numerical simulations, retaining the main physics ingredients at microscales (e.g., Howes et al. 2011;Servidio et al. 2012;Wan et al. 2015;Haggerty et al. 2017). In particular, large high-resolution hybrid Particle-In-Cell (PIC) numerical simulations (using a kinetic description for the ions and modeling the electrons as a fluid), later complemented by Hall-magnetohydrodynamic (HMHD) simulations, were able to reproduce most of the turbulent properties observed in the solar wind (Franci et al. 2015b,a;Franci et al. 2018a,b;Franci et al. 2019a;Papini et al. 2019b) and in the Earth magnetosheath . Such simulations also showed that the development of a turbulent cascade at sub-ion scales is concurrent to the onset of reconnection events in ion-scale current sheets (Franci et al. 2016c;Franci et al. 2017). Moreover, spectral properties at the reconnection exhausts consistent with a developed turbulent state were observed in a fully kinetic simulation (Pucci et al. 2017). Finally, recent works (Franci et al. 2017;Papini et al. 2019a,b) have quantitatively shown that current sheets undergoing reconnection in developing turbulence trigger an energy transfer directly from large to small scales, and can initiate a turbulent cascade that later establishes a proper inertial range, regardless of the model (MHD, Hall-MHD, or ion-kinetic) employed.
The ability of numerical experiments to reproduce the turbulent plasma properties is encouraging, as it confirms that the physical models employed are the right ones to explain spacecraft observations (in that range of scales and/or in those specific plasma conditions). Moreover, unlike in situ observations, numerical simulations provide the full spatiotemporal information needed to understand the plasma dynamics. Nevertheless, extracting such information is challenging. Traditional analysis techniques, based on Fourier or wavelet decomposition, have been successful in describing some statistical properties of turbulence (e.g., Bruno et al. 2001;Chang et al. 2004;Consolini et al. 2005;Horbury et al. 2008;Lion et al. 2016). Such methods, however, assume stationarity and/or linearity of the signal to be analyzed. Yet, turbulence is intrinsically nonlinear and nonstationary.
To address these limitations, Huang et al. (1998) developed the Empirical Mode Decomposition (EMD), a technique specifically designed for decomposing nonstationary nonlinear one-dimensional signals into a set of Intrinsic Mode Functions (IMF), that oscillate around zero but with varying frequency and amplitude. Such decomposition is adaptive, based on the local characteristic scales of the signal, and does not require any assumption on the shape of the signal to be extracted. EMD has proven to be a very powerful tool in many research areas and has recently been used to measure the multifractal properties of the solar wind (Alberti et al. 2019). Unfortunately, EMD shown to be unstable in presence of noise, and the Ensemble EMD (EEMD, Wu & Huang 2009) and similar alternative methods, which address this issue, greatly increase the computational costs and lack a rigourous mathematical theory behind them.
As an alternative to (E)EMD and equivalent techniques, algorithms based on Iterative Filtering (IF) have been recently developed (Lin et al. 2009;Cicone et al. 2016;Cicone 2020). Unlike EMD, they give a convergent solution for any square-integrable (L 2 ) signal, also in presence of noise. IF methods have already been successfully employed in the analysis of time-series from geomagnetic measurements Bertello et al. 2018;Spogli et al. 2019). Multidimensional Iterative Filtering (MIF) generalizes IF to high-dimensional signals, and represents the fastest and more robust adaptive multidimensional decomposition technique currently available in the literature. It outperforms other methods in terms of computational costs and, at the same time, it retains all the convergence properties of the one-dimensional IF algorithms (for more details, see Cicone & Zhou 2017;Cicone 2020).
In this work, we carry out the first multiscale analysis of numerical simulations of plasma turbulence by means of MIF decomposition. We focus on two numerical datasets, obtained from one HMHD and one HPIC simulation respectively.
Our results demonstrate the ability of MIF to: (i) separate the different turbulent regimes (the Energy injection scales, the MHD inertial range, the sub-ion kinetic regime, and the dissipation scales) while retaining the information about the magnetic field spatial configuration, (ii) disentangle the morphological and physical features of magnetically reconnecting current sheets, and (iii) quantify the statistical properties of turbulence.

Numerical simulations of plasma turbulence
The datasets used in this work were produced by two high-resolution numerical simulations of plasma turbulence, thoroughly caracterized in Papini et al. (2019b). a Hall-MHD and a Hybrid-PIC simulation.

The Hall-MHD model
The HMHD model takes into account two-fluids effects that describe the separate dynamics of ions and electrons at sub-ion scales. Different HMHD models can introduce several levels of complexity, depending on whether they retain a description for the pressure tensors and/or for electron inertia effects (see, e.g., Shay et al. 2001). Here we use a model that consists of the nonlinear viscous-resistive MHD equations, modified only by the presence of the Hall term in the induction equation. This is done by substituting the fluid velocity u with the electron velocity u e = u − J /en e . In their adimensionalized form, the HMHD equations take the form where Γ = 5/3 is the adiabatic index and {ρ, u, B, T, P } are a function of time and space and denote the usual variables. The equation of state P = ρT relates the gas pressure to the other two thermodynamic variables. All quantities are renormalized with respect to a characteristic length L = d i , a magnetic field amplitude B 0 , a plasma density ρ 0 , an Alfvén velocity Moreover Ω i = eB 0 /m i c is the ion-cyclotron angular frequency and m i is the mass of the ions. With this normalization, the (adimensional) magnetic resistivity η is in units of d i c A and the Hall coefficient The equations (2.1-2.4) are numerically solved by using a pseudospectral code we developed, already employed for studies of magnetic reconnection (Landi et al. 2015;Papini et al. 2018;Papini et al. 2019c) and plasma turbulence (Papini et al. 2019a,b). We consider a two-dimensional (x, y) periodic domain and use Fourier decomposition to calculate the spatial derivatives. In Fourier space we also filter according to the 2/3 Orszag rule, to avoid aliasing of the nonlinear terms. For the temporal evolution of {ρ, u, B, T } we use a 3rd-order Runge-Kutta scheme.

The Hybrid-PIC model
The second dataset was produced by using the Lagrangian HPIC code CAMELIA (Current Advance Method Et cycLIc leApfrog, Matthews 1994;Franci et al. 2018a). In CAMELIA, the ions are modeled as macroparticles that correspond to statisticallyrepresentative portions of the distribution function in the phase space. The plasma charge is neutralized by a massless and isothermal electron fluid. The system is governed by the Vlasov-Maxwell equations. Electron inertia effects and the displacement current in the Maxwell's equations are neglected. Therefore, only macroparticle's position and velocity inside each grid cell, as well as magnetic fields defined at the cell nodes, need to be evolved in time. All other quantities and moments are functions of the above quantities, including the electric field (Matthews 1994).
Among many applications, CAMELIA has been employed for numerical studies of plasma turbulence (e.g. Franci et al. 2015bFranci et al. ,a, 2016cFranci et al. ,a, 2017. It reproduced many of the spectral properties observed in the solar wind (Franci et al. 2018a) and in the Earth magnetosheath (Franci et al. 2019a) (we refer the reader to Franci et al. (2018a) for further details and applications).

Numerical setup
Apart from few parameters, the HMHD and the HPIC simulations employ the same setup. We consider a 2D box of size L x × L y = 256 d i × 256 d i and a grid resolution of ∆x = ∆y = d i /8, corresponding to 2048 2 points. The system is initialized with a constant mean magnetic field B 0 = B 0 e z out of the plane, along the z direction (that we will refer to as the parallel direction). The xy-plane (i.e. the perpendicular plane) is filled with freely-decaying random Alfvénic-like sinusoidal fluctuations. These are characterized by a root-mean-square amplitude b rms = B rms /B 0 0.24, and wavenumbers spanning from the smallest nonzero value contained in the box up to the injection scale inj = 2π/k inj ⊥ , such that k inj ⊥ d i 0.28, with k ⊥ = k 2 x + k 2 y . In the HPIC simulation, we set the ion and electron plasma beta to β i = β e = 1, while the magnetic resistivity has the value η = 5 × 10 −4 . The HMHD simulation has a (total) plasma β = β i + β e = 2, and a resistivity and a viscosity η = ν = 10 −3 .

Datasets of fully developed turbulence
In both simulations, the initial Alfvénic fluctuations quickly evolve to form coherent structures, namely vortices and, in between them, current sheets. The latter get disrupted by magnetic reconnection and release small-scale plasmoids which feed back to their turbulent sorrounding. The subsequent evolution is characterized by the formation and disruption of many other current sheets. Concurrently, a turbulent cascade develops at large scales, until a quasi-stationary state, shown in Figure 1, is reached at t = 165 τ A and at t = 200 τ A for the HMHD and the HPIC run respectively. At those times, the power spectrum of the total magnetic fluctuations (right panel of Fig. 1) has a clear multiscale behavior. At scales larger than the injection scale (k ⊥ < k inj ⊥ ) we have the reservoir of energy that fuels the cascade. At fluid MHD scales (k inj ⊥ < k ⊥ 2/d i ) a Kolmogorov-like power law of spectral index −5/3 is present, which then transitions to a slope of −3 at sub-ion kinetic scales at about k break ⊥ ∼ 2/d i (the so called spectral break). Finally, at k ⊥ > k diss ⊥ 12/d i we reach the dissipation scales. The physical and statistical properties of these four regimes are quite different, due to the diversity of the underlying physical mechanims acting at those scales. For instance, sub-ion scales are characterized by increasing levels of intermittency generated by the presence of thin localized current structures where dissipation is enhanced. We refer the reader to In the next sections, we will show how MIF methods can correctly separate these regimes.

Multidimensional Iterative Filtering
We now introduce the Multidimensional Iterative Filtering technique. For more details about MIF methods, as well as applications and examples, we remind the reader to Cicone & Zhou (2017).
Given a (multidimensional) signal, f (r) with r R k , MIF decomposes it into a finite number N of (locally almost orthogonal) simple oscillating components f called Intrinsic Mode Functions (IMF) where r f,N is the residual of the decomposition (ideally, a trend signal). Each f j is the result of an iterative procedure that uses a low-pass filter to extract the moving average of the signal at a given scale λ j , so to isolate a fluctuating component whose average frequency ν j 1/λ j is well behaved. λ j is different for each IMF and increasing with j. Therefore, IMFs with increasing j will contain larger (smaller) scales (frequencies). We first specify the low-pass filter operator that acts on a L 2 signal s(r). Here w j (r) Ω(λ j ) is the kernel function associated to the filter, and Ω(λ j ) ⊂ R k is the spherical support of w j with radius λ j (e.g., a circle in R 2 ). In this work we use a two-dimensional isotropic kernel function (see Figure 2) Let us now define S 1,0 (r) = f (r) and introduce the fluctuating function S 1,1 (r) = S 1,0 (r) − L 1 [S 1,0 (r)]. The scale λ 1 of L 1 is chosen such that its lenght is comparable to the maximum frequency contained in S 1,0 (for more details on the choice of λ 1 , see Lin et al. 2009;Cicone et al. 2016). Through iteration, one can calculate S 1,n (r) = S 1,n−1 (r) − L 1 [S 1,n−1 (r)]. The first IMF is obtained in the limit of infinite iterations f 1 (r) = lim n→∞ S 1,n (r), (3.3) and the residual signal is The second f can be calculated by defining S 2,0 (r) = r 1 (r) and repeating the above procedure but using a kernel function with a larger radius λ 2 > λ 1 , whose value is chosen based on S 2,0 . The j-th IMF is given by and λ j > λ j−1 > . . . > λ 1 .
( 3.7) The decomposition ideally ends when the residual r f,N (r) is a trend signal, that is, r f,N contains no local extrema. In practice this is achieved by requiring that r f,N (r) contains less than a given number of local extrema. In the following, we drop the N subscript to simplify notation and write the residual of a signal f as r f .

Multiscale analysis of fully developed Turbulence
We now describe the multiscale analysis performed on the turbulent magnetic fields from the two numerical datasets described in Section 2.4 and shown in Fig. 1. In Figure  3, we report the MIF decomposition of the out-of-plane component of the magnetic field, B z , from the HPIC run. Eleven IMFs { B z,j with j = 1, ..., 11} have been extracted. The residual r Bz (x, y) is shown on the bottom right panel. The dissipation and the ion kinetic scales (high spatial frequencies, k ⊥ d i > 2) are captured by the first four IMFs and are characterized by well localized structures. Going to larger scales (k ⊥ d i < 2), these features become more homogenously distributed. That happens, for instance, in the MHD inertial range, captured by B z,5 , B z,6 , and B z,7 . Finally, the largest scales (above the injection) are contained in the last IMFs and in the residual, which also retain the mean field B 0 . Similar results (not shown here) are obtained for B x and B y , as well as for the MIF decomposition of the HMHD simulation. Unlike wavelets and Fourier modes, the IMFs are only locally almost orthogonal. Therefore, it may be useful to assess the degree of orthogonality. The orthogonality of a set { f i (x, y)} is given by the (symmetric) orthogonality matrix M, with elements where The set is orthogonal if M ij = δ ij for each i, j. Figure 4 shows the orthogonality matrix of the IMFs of the out-of-plane magnetic field fluctuations (see Fig. 3). As espected, the  set is not orthogonal, since for neighbor IMFS the lower diagonal (indices (i + 1, i)) of the orthogonality matrix reaches a maximum value of 0.77 and a mean of 0.34. However, for second neighbors ((i + 2, i) pairs) the values drop to less than 0.08 (except for one point).
The components of the magnetic fluctuations at large injection, inertial range/MHD, ion kinetic, and at dissipation scales are further separated by regrouping the IMFs in four aggregated (vector) IMFs { B inj , B MHD , B kin , B diss }. Their components (for, e.g.,  Fig. 6. ⊥ is the average spatial wavenumber of B x,j and, for each aggregated IMF, the sum is performed over the range of scales of interest. The amplitude of the agregated IMFs (Fig. 5) reveals that the HPIC and the HMHD run are morfologically equivalent, characterized by homogenously distributed features at large scales. As the scales decrease, such features become more and more localized, self-organizing in a filamented network at the edge of the turbulent eddies, where magnetic dissipation is enhanced. The isotropized power spectra of the agregated IMFs are shown in Fig. 6. The MHD inertial range, the kinetic range, and the dissipation scales, as well as the injection scales range, are well separated also in Fourier space. This further confirms the ability of the MIF decomposition to succesfully isolate the four different regimes, while retaining the full spatially local information of the fields.
As a final remark, the agregated IMFs have an increased orthogonality, as the maximum value out of the diagonal in their ortoghonality matrix is about 0.10.

Current structures and intermittency
Intermittency in plasma turbulence is related to the dynamics of current sheets and localized coherent structures, since they break the self-similarity of the system. Usually, such structures are found in a sort of multifractal configuration, such as the filamented network we observed in Figure 5. With MIF we can easily isolate these features. As an example, Figure 7   . Isotropized 1D power spectra of the total magnetic field fluctuations B (black dashes) and of the injection scales (blue), the MHD scales (orange), the ion kinetic scales (green), and the dissipation scales (red) agregated IMFs, for the HMHD run (left) and for the HPIC run (right). Vertical dashed lines denote k inj ⊥ , k ⊥ di = 2, and k ⊥ di = 12, i.e., the three wavenumbers that approximately separate the four regimes. simulation box, which hosts a chain of three plasmoids originated from the disruption of a reconnected current sheet. The amplitude of the original magnetic field fluctuations and of its current density, J = ∇ × B, is shown on the leftmost top and bottom panels respectively. There, we can distinguish the small plasmoids in the magnetic amplitudes, while the corresponding signal in the current density is almost swamped by the particleper-cell (PPC) noise and by dissipation. The amplitudes of the aggregated IMFs (top) and of their associated current densities (bottom) are shown in the other panels. Now the three plasmoids are clearly visible in B kin and J kin , and their large-scale signature also appears at MHD scales. As espected, the aggregated IMFs also reveal that magnetic dissipation is mostly concentrated in strong current structures (high values of J diss are found in areas of high J kin ). Moreover, the highest dissipation (bright spots of J diss ) takes place at the X-points between the plasmoids, a typical feature of magnetic reconnection events. All these morphological and physical multiscale properties, tipically difficult to isolate, are nicely and directly disentangled using MIF decomposition and in a straightforward manner.
Following Frisch (1995), a quantitative measure of the intermittency of a signal f (x, y) can be provided by measuring, at different frequencies k (j) ⊥ , the kurtosis where f > j (x, y) is a high-pass filtered signal of f (x, y), only containing frequencies k ⊥ k (j) ⊥ , and where denote a spatial average. The signal f (x, y) is intermittent if its kurtosis grows without bounds with frequency. Using MIF decomposition, we can write the filtered signal as by summing over all the IMF with an average frequency k ⊥ . We point out that Iterative Filtering based methods proved to be well suited to reconstruct the kurtosis, and more in general the multiscale statistical properties of nonstationary signals (including intermittent ones), as recently shown in Stallone et al. (2020).
Another quantity that measure the departure from a gaussian behavior is the Kullback-Leibler (KL) divergence (e.g., Granero-Belinchón et al. 2018). For a sample X with probability density function (PDF) p(x) of variance σ 2 X , the KL divergence is defined as where is the Shannon entropy of the sample, and H G = 0.5 log(2πeσ 2 X ) is the Shannon entropy that X would have if p(x) were a Gaussian. K L (X) is always positive, being zero if the PDF of the sample is a Gaussian distribution. Finally, the KL divergence K L (f > j ) of f (x, y) at a given frequency k (j) ⊥ is obtained by calculating the PDF and the variance of the values of f > j (x, y). Figure 8 shows the excess kurtosis, K(f > j ) − 3, together with the KL divergence of the magnetic field components of the HMHD (left) and the HPIC (right) datasets. The results, in qualitative agreement with our previous findings , show a flat kurtosis at scales above the injection scale (k ⊥ d i 0.28), denoting a selfsimilar behavior at those scales. The kurtosis of the perpendicular fluctuations start increasing in the MHD inertial range, then it steepens abruptly at kinetic scales, where Hall-current effects become important. This is in agreement with what observed in the aggregated IMFs: the presence of filamented networked structures at kinetic scales implies strong intermittency. The kurtosis of the parallel fluctuations B z , instead, remains constant in the inertial MHD range down to k ⊥ d i = 1, then abrutly increases with a behavior seimilar to the perpendicular fluctuations. We interpret such difference between parallel and perpendicular fluctuations at MHD scales as due to the particular setup we choose. Alternatively, this could be the signature of the different behaviour of the two regimes: at the MHD scales, turbulence is leaded by Alfvénic-like fluctuations, which are predominately polarized in the direction perpendicular to B 0 , while in the kinetic regime, dispersive effects couples the fluctuations with B z . The simulations are inizialized with Alfvénic fluctuations in the perpendicular plane. In the inertial range, this causes the formation of a turbulent cascade in the perpendicular magnetic fluctuations, which then transition to kinetic scales. Instead, the magnetic power spectrum of parallel fluctuations (see Fig. 5 of Papini et al. 2019b) show a cascade only at kinetic scales. Consequently, no intermittency develops until the disruption scales of current sheet (at around the ion inertial length d i ) are reached.
The KL divergence (bottom panels of Figure 8) shows the same results. A departure from gaussian behavior (K L > 0) is observed in the perpendicular fluctuations right below the injection scales and in the parallel fluctuations at kinetic scales. Interestingly, unlike the HMHD dataset, both the kurtosis and the KL divergence of B z in the HPIC dataset, although constant (i.e. not intermittent) at fluid scales (k ⊥ d i < 1), are not zero, which denote a non-gaussian nature of the fluctuations. We finally note that at the smallest scales, where dissipation kicks in, the KL divergence decreases in both the datasets, as expected.

Discussion
In this work, we investigated the magnetic multiscale properties of both fluid Hall-MHD and hybrid ion-kinetic electron-fluid simulations of plasma turbulence by means of Multidimensional Iterative Filtering (MIF), a novel technique for the decomposition of nonstationary multidimensional signals. By exploiting our large-scale high-resolution numerical datasets, we succesfully separated the four ranges of scales relevant to turbulence, namely the large injection scales, the inertial-range MHD scales, the sub-ion kinetic scales, and the dissipation scales. Moreover, we were able to reproduce the spectral and statistical properties of such regimes, while preserving the spatial information about morphology and localization of features and coherent structures, such as current sheets and vortices.

Intermittency
Our results confirm that plasma turbulence is an intrinsic multiscale phenomenon. In the inertial range, the energy cascade consists of more or less homegenously distributed and weakly intermittent magnetic fluctuations (a slowly increasing kurtosis in the perpendicular fluctuations is observed). Ion kinetic scales are characterized by stronglylocalized coherent structures organized in a filamented network, which show a high degree of intermittency. Finally, when reaching the dissipation scales, the kurtosis tends to flatten and the KL-divergence decreases, suggesting that intermittency is switching off. These results have been obtained by measuring the scale-dependent kurtosis and the KL-divergence of the magnetic field fluctuations, by means of a statistical analysis that exploits the MIF decomposition to calculate the high-pass filtered field f > j contaning all spatial frequencies k ⊥ k (j) ⊥ (see Eq. 4.5). It is instructive to compare these results with those obtained by using (i) Fourier transform in place of MIF decomposition to compute f > j (that is the exact definition of Frisch 1995) and (ii) the magnetic field increments instead of f > j . The latter method has been applied by Papini et al. (2019b) to the same simulation dataset used in this work, and it is extensively employed both in solar wind and magnetosheath observations (see, e.g., Koga et al. 2007;Chian & Miranda 2009;Wu et al. 2013;Bandyopadhyay et al. 2018), as well as in numerical simulations (e.g., Wan et al. 2012;Franci et al. 2015a;Haggerty et al. 2017).
In Figure 9 we report, for the HPIC run, the kurtosis of the increments ∆B ,y y = B y (x, y + ) − B y (x, y) (green curve) and ∆B ,x y = B y (x + , y) − B y (x, y) (red curve), respectively parallel and perpendicular (in the plane) to the direction of the B y component (where = 2π/k ⊥ ). We also report the kurtosis obtained from the MIF decomposition of B y (already shown in the top-right panel of Figure 8) and the one obtained using Fourier transform to calculate the high-pass filtered field f > j from B y . Overall, the results are in qualitative agreement. There are, however, some noticeable differences. The kurtosis of the increments ∆B ,x y is almost constant in the inertial range, then linearly increases starting from k ⊥ d i 1, i.e., at the scales where Hall currents become important. Instead, K(∆B ,y y ) linearly increases already at the injection scales, but with a smaller slope than K(∆B ,x y ). The kurtosis calculated with MIF, although it is roughly twice as larger, reproduces the properties of both K(∆B ,y y ) and K(∆B ,x y ), following the former at the large scales down to the spectral break k ⊥ d i 2 and then steepening at the scales where K(∆B ,x y ) > K(∆B ,y y ). Such differences may be explained by considering the following argument. Firstly, the increments ∆B ,y y ∝ ∂ y B y and ∆B ,x y ∝ ∂ x B y are proportional to the components of the magnetic field gradient parallel and perpendicular to the magnetic field direction, respectively. Secondly, localized and elongated magnetic structures, such as the thin current sheets between vortices that we observed, are aligned with the magnetic field. Therefore, K(∆B ,x y ) is sensitive to the thickness of such structures (which is of the order of d i ) while K(∆B ,y y ) probes their length (that is comparable to the size of the largest vortices). This may explain the observed increase in the parallel and perpendicular . Excess kurtosis of the y-component of the magnetic field fluctuations from the HPIC dataset, as calculated from the increments along the x-direction (∆B ,x y ) and along the y-direction (∆B ,y y ). The excess kurtosis as given by Eq. 4.4, by using MIF (solid orange curve and circles) and FFT (blue dashed curve) methods to calculate the high-pass filtered function f > j of By, is also shown.
kurtosis at around those scales. We finally remark that the kurtosis calculated using the Fourier and the MIF decomposition are in remarkable agreement, although the former shows smaller values at high frequencies and even bring to a sudden decrease at the highest frequencies (because of the increasing inability of the Fourier high-pass filtered field to localize structures, as k ⊥ increases). Analogous results are obtained for the xcomponent of the magnetic field. Overall, the multiscale statistical properties we recovered are consistent with the findings of Alberti et al. (2019) (hereafter AL19). They measured the multifractal nature of solar wind turbulence by using CLUSTER data and found increasing levels of intermittency in the MHD/inertial range, with a tendency toward a non-intermittent/monofractal behavior at dissipation scales. There is, however, an important difference. In the ion kinetic range, where AL19 observe a monofractal behavior, we find high levels of intermittency, which denote a multifractal nature of the fluctuations. We are not certain whether such difference is due to the particular dataset chosen by AL19 (a fast solar wind stream) or by our HMHD and HPIC models. We note, however, that at ion kinetic scales AL19 measure magnetic power spectra with a slope of ∼ −5/2, different from the one we report in our simulations (−3) and also from other solar wind conditions. Further investigation pursuing this path is currently underway, in order to assess the multifractal properties of our numerical simulations.

Reconnection and enhanced dissipation
Our multiscale analysis confirms that magnetic field dissipation is mostly concentrated in the filamented magnetic network, especially at the X-point of reconnecting current sheets. This is in agreement with previous studies of plasma turbulence that measured the scale-to-scale energy transfer by means of scale-filtering approaches Camporeale et al. 2018;Kuzzay et al. 2019). The filamented magnetic network observed at kinetic scales is somewhat reminiscent of the vortex filaments in hydrodynamic turbulence (e.g., Kida & Ohkitani 1992;Moffatt et al. 1994), and consistent with the fact that areas of both high magnetic gradients and vorticities are correlated (Franci et al. 2016b;Kuzzay et al. 2019), especially at the reconnection sites, which produce high levels of vorticity (e.g., Widmer et al. 2016). Furthermore, the physical and geometrical features typical of magnetic reconnection were easily disentangled and identified, even when the signal of the structure was swamped by PPC noise and dissipation (see Figure 7). In this context, Multidimensional Iterative Filtering also stands as a powerful tool for automatically identifying and removing noise from the physical signal of interest.
Currently, we are conducting a time-frequency analysis of our simulations, by employing IF. The aim is to confirm whether the turbulent dynamics at kinetic scales is wavelike (e.g. mediated by kinetic alfvén wavelike interactions) or it is due to the presence of structures, as the results of this work seems to suggest. Overall, Multidimensional Iterative Filtering is a promising technique to support the study of plasma turbulence and its properties, with many potential applications in both numerical simulations and spacecraft observations. E. Papini and S. Landi thank T. Alberti and G. Consolini for useful discussion. M. Piersanti thanks the Italian Space Agency for the financial support under the contract ASI LIMADOU scienza n • 2016-16-H0." This research was partially supported by the UK Science and Technology Facilities Council (STFC) grant ST/P000622/. This work was supported by the Programme National PNST of CNRS/INSU co-funded by CNES P. Hellinger acknowledges grant 18-08861S of the Czech Science Foundation. We acknowledge partial funding by "Fondazione Cassa di Risparmio di Firenze" under the project HIPERCRHEL. The authors acknowledge the "Accordo Quadro INAF-CINECA (2017)", for the availability of high performance computing resources and support, PRACE for awarding access to the resource Cartesius based in the Netherlands at SURFsara through the DECI-13 (Distributed European Computing Initiative) call (project HybTurb3D), and CINECA for awarding access to HPC resources under the ISCRA initiative (grants HP10B2DRR4 and HP10C2EARF).