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

Abstract

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 non-stationary 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 localized current structures to be detected and their contribution to the statistical and spectral properties of turbulence to be characterized. Overall, MIF arises as a very promising technique for the study of turbulent plasma environments.

1. Introduction: turbulence and intermittency in space plasmas

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, Osman & Chapman 2015) and in the Earth's magnetosheath (Stawarz et al. 2016; Chen & Boldyrev 2017) show that power spectra of magnetic fluctuations exhibit power-law behaviours 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 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 typically around $-3$), until dissipation scales are reached (for a comprehensive review, see Verscharen, Klein & Maruca 2019, and references therein).

Owing 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 with the deformation of the particles velocity distribution functions (Yang et al. 2017; Del Sarto & Pegoraro 2018). 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 is 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 that lead to a spatially inhomogeneous energy cascade and turbulent 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 generally 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, Schekochihin & Chandran 2017; Landi et al. 2019).

Nowadays, investigating plasma turbulence using direct approaches is becoming increasingly feasible. The increasing computational capabilities allow direct numerical simulations to be run that retain 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 (HPIC) numerical simulations (using a kinetic description for the ions and modelling 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. 2015a,b, 2018a,b, 2020; Papini et al. 2019b) and in the Earth magnetosheath (Franci et al. 2020). 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. 2016b; Cerri & Califano 2017; 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, HMHD 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 correct 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, Tam & Wu 2004; Consolini, Chang & Lui 2005; Horbury, Forman & Oughton 2008; Lion, Alexandrova & Zaslavsky 2016). Such methods, however, assume stationarity and/or linearity of the signal to be analysed. Yet, turbulence is intrinsically nonlinear and non-stationary.

To address these limitations, Huang et al. (1998) developed the Empirical Mode Decomposition (EMD), a technique designed specifically for decomposing non-stationary nonlinear one-dimensional signals into a set of intrinsic mode functions (IMFs), 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 developed recently (Lin, Wang & Zhou 2009; Cicone, Liu & Zhou 2016; Cicone 2020). Unlike EMD, they give a convergent solution for any square-integrable ($L^2$) signal, also in the presence of noise. IF methods have already been successfully employed in the analysis of time series from geomagnetic measurements (Bertello et al. 2018; Piersanti et al. 2018; Spogli et al. 2019). Multidimensional Iterative Filtering (MIF) generalizes IF to high-dimensional signals, and represents the fastest and most 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, 2020; 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.

2. Numerical simulations of plasma turbulence

The datasets used in this work were produced by two high-resolution numerical simulations of plasma turbulence, employing a fluid HMHD and a HPIC model, respectively. Their evolution into fully developed turbulence, as well as their statistical properties, have been thoroughly characterized in Papini et al. (2019b).

2.1. The HMHD model

The HMHD model takes into account two-fluid 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 $\boldsymbol {u}$ with the electron velocity $\boldsymbol {u}_e = \boldsymbol {u} - \boldsymbol {J}/ e n_e$, $\boldsymbol {J}$ being the current density. In their adimensionalized form, the HMHD equations take the form

(2.1)\begin{gather} \partial_t \rho + \boldsymbol{\nabla}\boldsymbol{\cdot}{(\rho \boldsymbol{u})} = 0, \end{gather}

(2.2)\begin{gather} \rho\left (\partial_{t} + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\right) \boldsymbol{u} = -\boldsymbol{\nabla} P + (\boldsymbol{\nabla}\times{\boldsymbol{B}})\times\boldsymbol{B} + \nu \left [ \nabla^2\boldsymbol{u} + \tfrac{1}{3}\boldsymbol{\nabla}(\boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{u}}) \right], \end{gather}

(2.3)\begin{gather} \left (\partial_{t} + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \right) T = (\varGamma - 1) \left \{ - (\boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{u}})T + \eta \frac{|\boldsymbol{\nabla}\times\boldsymbol{B}|^2}{\rho} + \frac{\nu}{\rho} \left [ (\boldsymbol{\nabla}\times{\boldsymbol{u}})^2 +\frac{4}{3} (\boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{u}})^2\right ] \right \} , \end{gather}

(2.4)\begin{gather} \partial_t{\boldsymbol{B}} = \boldsymbol{\nabla}\times\left ( \boldsymbol{u} \times \boldsymbol{B}\right ) + \eta \nabla^2 \boldsymbol{B} - \eta_\mathrm{H}\nabla\times \frac{(\nabla\times\boldsymbol{B})\times\boldsymbol{B}}{\rho}, \end{gather}

where $\varGamma =5/3$ is the adiabatic index and $\{\rho ,\boldsymbol {u},\boldsymbol {B},T,P\}$ are a function of time and space and denote the usual variables. The equation of state $P=\rho 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 $\rho _0$, an Alfvén velocity $c_A = B_0/\sqrt {4{\rm \pi} \rho _0} = \varOmega _i d_i$, a pressure $P_0=\rho _0 c_A^2$ and a plasma temperature $T_0 = (k_B/m_i) P_0/\rho _0$. Moreover, $\varOmega _i = e B_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 $\eta$ is in units of $d_i c_A$ and the Hall coefficient $\eta _H=d_i/L$ is equal to 1.

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, Landi & Del Zanna 2018, 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 $\{\rho ,\boldsymbol {u},\boldsymbol {B},T\}$ we use a third-order Runge–Kutta scheme (Wray 1990).

2.2. The HPIC 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 modelled as macroparticles that correspond to statistically representative 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 the position and velocity of the macroparticles 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. 2015a,b, 2016a,b, 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. 2020) (we refer the reader to Franci et al. (2018a) for further details and applications).

2.3. Numerical setup

Apart from one parameter, the HMHD and the HPIC simulations employ the same setup. We consider a two-dimensional box of size $L_x \times L_y = 256 d_i \times 256 d_i$ and a grid resolution of $\Delta x = \Delta y = d_i/8$, corresponding to $2048^2$ points. The system is initialized with a constant mean magnetic field $\boldsymbol {B}_0 = B_{0} \boldsymbol {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 (rms) amplitude $b_{\mathrm {rms}} = B_{\mathrm {rms}}/B_0 \simeq 0.24$, and wavenumbers spanning from the smallest non-zero value contained in the box up to the injection scale $\ell _{\mathrm {inj}}=2{\rm \pi} /k_{\perp }^{\mathrm {inj}}$, such that $k_{\perp }^{\mathrm {inj}} d_i \simeq 0.28$, with $k_{\perp } = \sqrt {k_x^2+k_y^2}$. In the HPIC simulation, we set the ion and electron plasma beta to $\beta _i=\beta _e=1$, whereas the magnetic resistivity has the value $\eta =5\times 10^{-4}$. Such value is chosen so to prevent the accumulation of magnetic energy at small scales caused by the use of a finite number of macroparticles (Franci et al. 2015a). The HPIC run employs 64 000 particles per cell (PPC). The HMHD simulation has a (total) plasma $\beta =\beta _i+\beta _e=2$, and a resistivity and a viscosity $\eta =\nu =10^{-3}$, such that the ion-kinetic scales and the Kolmogorov-dissipation scales are well separated.

2.4. 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 are disrupted by magnetic reconnection and release small-scale plasmoids that feed back to their turbulent surrounding. 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, corresponding to the maximum of the rms amplitude of the current density (Mininni & Pouquet 2009), is reached at $t=165 \tau _A$ and at $t=200 \tau _A$ for the HMHD and the HPIC run, respectively (see figure 1). At those times, the power spectrum of the total magnetic fluctuations (figure 1c) has a clear multiscale behaviour. At scales larger than the injection scale ($k_{\perp } < k_{\perp }^{\mathrm {inj}}$) we have the reservoir of energy that fuels the cascade. At fluid MHD scales ($k_{\perp }^{\mathrm {inj}} < k_{\perp } \lesssim 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_{\perp }^\mathrm {break} \sim 2/d_i$ (the so-called spectral break). Finally, at $k_{\perp } > k_{\perp }^{\mathrm {diss}} \simeq 12/d_i$ we reach the dissipation scales. The physical and statistical properties of these four regimes are quite different, owing to the diversity of the underlying physical mechanisms 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 turbulent dissipation is enhanced. In the following sections, we show how MIF methods can separate these regimes correctly.

Figure 1. Coloured contours of the amplitude of the total magnetic field fluctuations $|\boldsymbol {B}_{\textrm {tot}} - \boldsymbol {B}_0|$, at the maximum of the turbulent activity $(a)$ for the HMHD run at $t=165\tau _A$ and $(b)$ for the HPIC run at $t= 200 \tau _A$. The corresponding isotropized power spectra (Papini et al. 2019b) (solid red curve for the HMHD dataset and dashed blue curve for the HPIC run) are shown in $(c)$. Vertical dotted and dot-dashed lines denote the injection wavenumber $k_{\perp }^{\mathrm {inj}}$ and the Nyquist wavenumber, respectively.

3. Multidimensional Iterative Filtering

We now introduce the MIF technique. For more details about MIF methods, as well as applications and examples, the reader may refer to Cicone & Zhou (2017).

Given a (multidimensional) signal, $f(\boldsymbol {r})$ with $\boldsymbol {r} \epsilon \mathbb {R}^k$, MIF decomposes it into a finite number $N$ of (locally almost orthogonal) simple oscillating components $\,\hat {f}$ called IMFs

(3.1)\begin{equation} f(\boldsymbol {r}) = \sum_{j=1}^N \,\hat{f}_j(\boldsymbol {r}) + r_{f,N}(\boldsymbol {r}), \end{equation}

where $r_{f,N}$ is the residual of the decomposition (ideally, a trend signal as explained in the following). Each $\,\hat {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 $\lambda _j$, so as to isolate a fluctuating component whose average spatial frequency $\nu _j\simeq 1/\lambda _j$ is well behaved. For each IMF $\lambda _j$ is different and increasing with $j$. Therefore, IMFs with increasing $j$ will contain larger (smaller) scales (frequencies).

We first specify the low-pass filter operator

(3.2)\begin{equation} \mathcal{L}_j [s(\boldsymbol {r})] = \int_{\varOmega(\lambda_j)} s(\boldsymbol {r} + \boldsymbol {t}) w_j(\boldsymbol {t})\,\textrm{d}\boldsymbol {t} \end{equation}

that acts on a $L^2$ signal $s(\boldsymbol {r})$. Here $w_j(\boldsymbol {r}) \,\epsilon \, \varOmega (\lambda _j)$ is the kernel function associated to the filter, and $\varOmega (\lambda _j) \subset \mathbb {R}^k$ is the spherical support of $w_j$ with radius $\lambda _j$ (e.g. a circle in $\mathbb {R}^2$). In this work, we employ periodic boundary conditions and use a two-dimensional isotropic kernel function (see figure 2), whose radial profile is calculated by numerically solving a Fokker–Planck differential equation (Cicone et al. 2016; Cicone & Zhou 2017; Cicone & Dell'Acqua 2020).

Figure 2. Kernel function $w_j(x,y)$ for a MIF decomposition of a two-dimensional field. The red circle of radius $\lambda _j=5$ denotes the boundary of $\varOmega (\lambda _j)$.

Let us now define $S_{1,0}(\boldsymbol {r}) = f(\boldsymbol {r})$ and introduce the fluctuating function $S_{1,1}(\boldsymbol {r}) = S_{1,0} (\boldsymbol {r})- \mathcal {L}_1[S_{1,0}(\boldsymbol {r})]$. The scale $\lambda _1$ of $\mathcal {L}_1$ is chosen such that its length is comparable with the maximum spatial frequency contained in $S_{1,0}$ (for more details on the choice of $\lambda _1$, see Lin et al. 2009; Cicone et al. 2016). Through iteration, one can calculate $S_{1,n}(\boldsymbol {r}) = S_{1,n-1}(\boldsymbol {r}) - \mathcal {L}_1[S_{1,n-1}(\boldsymbol {r})]$. The first IMF is obtained in the limit of infinite iterations

(3.3)\begin{equation} \hat{f}_1(\boldsymbol {r}) = \lim_{n\rightarrow\infty} S_{1,n}(\boldsymbol {r}), \end{equation}

and the residual signal is

(3.4)\begin{equation} r_{f,1} (\boldsymbol {r})= f(\boldsymbol {r}) - \,\hat{f}_1(\boldsymbol {r}). \end{equation}

The iteration procedure that calculates the limit in (3.3) stops when the desired convergence accuracy is reached (Cicone et al. 2016). The second IMF $\,\hat {f}_2$ can be calculated by defining $S_{2,0}(\boldsymbol {r}) = r_{f,1}(\boldsymbol {r})$ and repeating the above procedure but using a kernel function with a larger radius $\lambda _2 > \lambda _1$, whose value is chosen based on $S_{2,0}$.

The $j$th IMF is given by

(3.5)\begin{equation} \hat{f}_j(\boldsymbol {r}) = \lim_{n\rightarrow\infty} S_{j,n}(\boldsymbol {r}), \end{equation}

with

(3.6a,b)\begin{equation} S_{j,n}(\boldsymbol {r}) = S_{j,n-1}(\boldsymbol {r}) - \mathcal{L}_j[S_{j,n-1}(\boldsymbol {r})], \quad S_{j,0}(\boldsymbol {r}) = r_{f,j}(\boldsymbol {r}) \end{equation}

and

(3.7)\begin{equation} \lambda_j > \lambda_{j-1} > \ldots > \lambda_1. \end{equation}

The decomposition ideally ends when the residual $r_{f,N}(\boldsymbol {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}(\boldsymbol {r})$ contains less than a given number of local extrema. We set this number to three, so that the residual will retain the largest wavelength contained in the domain. In the following, we drop the subscript $N$ to simplify notation and write the residual of a signal $f$ as $r_f$.

4. 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 § 2.4 and shown in figure 1. For each field, the computation of all the IMFs took approximately one minute on our workstation (equiped with an Intel Xeon Gold 5218 CPU and 192 GB of RAM). 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 $\{\hat {B}_{z,j}$ with $j=1,\ldots ,11\}$ have been extracted. The residual $r_{B_z}(x,y)$ is shown in panel (l). The dissipation and the ion kinetic scales (high spatial frequencies, $k_{\perp } d_i > 2$) are captured by the first four IMFs and are characterized by well-localized structures. Going to larger scales ($k_{\perp } d_i < 2$), these features become more homogeneously distributed. That happens, for instance, in the MHD inertial range, captured by $\hat {B}_{z,5},\hat {B}_{z,6}$, and $\hat {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.

Figure 3. Colored contours of the IMFs $\hat {B}_{z,j}$ resulting from the MIF decomposition of the out-of-plane magnetic field fluctuations $B_z$ of the HPIC dataset (from small to large scales going from (a,e,i) to (d,h,l) and from (a–d) to (i–l)). The residual $r_{B_z}$ $(h)$ contains the largest-scale field fluctuations and the mean field $B_0$. For each IMF $\hat {B}_{z,j}$, its average spatial wavenumber $k_{\perp }^{(j)} = 2{\rm \pi} \nu _j$ (see § 3) is reported.

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 $\{\,\hat {f}_i(x,y)\}$ is given by the (symmetric) orthogonality matrix $\boldsymbol {M}$, with elements

(4.1)\begin{equation} \mathrm{M}_{ij}=\langle \,\hat{f}_i,\,\hat{f}_j\rangle = \frac{1}{ \|\,\hat{f}_i\|\cdot\|\,\hat{f}_j\|} \left |\int_0^{L_x} \int_0^{L_y} \,\hat{f}_i (x,y) \,\hat{f}_j (x,y)\, \mathrm{d}x\,\mathrm{d}y \right |, \end{equation}

where

(4.2)\begin{equation} \|\,\hat{f}_i\| = \left ( \int_0^{L_x} \int_0^{L_y} \,\hat{f}_i (x,y)^2 \,\mathrm{d}x\,\mathrm{d}y \right )^{1/2}. \end{equation}

The set is orthogonal if $\mathrm {M}_{ij} = \delta _{ij}$ for each $i,j$.

Figure 4 shows the orthogonality matrix of the IMFs of the out-of-plane magnetic field fluctuations (see figure 3). As expected, the set is not orthogonal, because for neighbour 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 neighbours ($(i+2,i)$ pairs) the values drop to less than $0.08$ (except for one point).

Figure 4. IMF orthogonality matrix $\boldsymbol {M}$ of the out-of-plane magnetic field fluctuations, as given by (4.1). Only the lower triangle is shown.

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 $\{ \hat {\boldsymbol {B}}_{\mathrm {inj}}, \hat {\boldsymbol {B}}_\mathrm {MHD}, \hat {\boldsymbol {B}}_\mathrm {kin}, \hat {\boldsymbol {B}}_{\mathrm {diss}} \}$. Their components (e.g. $B_x$) are defined as

(4.3a–d)\begin{align} \hat{B}_{x,\mathrm{inj}} =\sum_{k_{\perp}^{(j)}=0}^{k_{\perp}^{\mathrm{inj}}} \hat{B}_{x,j} + r_{B_x},\quad \hat{B}_{x,\mathrm{MHD}} = \sum_{k_{\perp}^{(j)}>k_{\perp}^{\mathrm{inj}}}^{2/d_i} \hat{B}_{x,j}\nonumber\\ \hat{B}_{x,\mathrm{kin}} = \sum_{k_{\perp}^{(j)}>2/d_i}^{k_{\perp}^{\mathrm{diss}}} \hat{B}_{x,j},\quad \hat{B}_{x,\mathrm{diss}} = \sum_{k_{\perp}^{(j)}>k_{\perp}^{\mathrm{diss}}} \hat{B}_{x,j}, \end{align}

where $k_{\perp }^{(j)}$ is the average spatial wavenumber of $\hat {B}_{x,j}$ and, for each aggregated IMF, the sum is performed over the range of scales of interest.The amplitude of the aggregated IMFs (figure 5) reveals that the HPIC and the HMHD run are morphologically equivalent, characterized by homogeneously distributed features at large scales. As the scales decrease, such features become increasingly localized, self-organizing in a filamented network at the edge of the turbulent eddies, where ohmic dissipation is enhanced. The isotropized power spectra of the aggregated IMFs are shown in figure 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 successfully isolate the four different regimes, while retaining the full spatially local information of the fields.

Figure 5. From (a,e) to (d,h): amplitude of the aggregated IMFs (see (4.3a–d)) of magnetic fluctuations at the injection ($|\hat {\boldsymbol B}_{\mathrm {inj}}-\boldsymbol {B}_0|$), fluid ($|\hat {\boldsymbol B}_{\mathrm {MHD}}|$), ion kinetic ($|\hat {\boldsymbol B}_{\mathrm {kin}}|$), and dissipation ($|\hat {\boldsymbol B}_{\mathrm {diss}}|$) scales. (a–d) and (e–h) refer to the HMHD and to the HPIC run, respectively. The corresponding power spectra are shown in figure 6.

Figure 6. Isotropized one-dimensional power spectra of the total magnetic field fluctuations $\boldsymbol {B}$ (black dashes) and of the injection scales (blue), the MHD scales (orange), the ion kinetic scales (green) and the dissipation scales (red) aggregated IMFs, for $(a)$ the HMHD run and $(b)$ the HPIC run. Vertical dashed lines denote $k_{\perp }^{\mathrm {inj}}$, $k_{\perp } d_i=2$ and $k_{\perp } d_i =12$, i.e. the three wavenumbers that approximately separate the four regimes.

As a final remark, the aggregated IMFs have an increased orthogonality, as the maximum value out of the diagonal in their orthogonality matrix is approximately 0.10.

4.1. Current structures and intermittency

Intermittency in plasma turbulence is related to the dynamics of current sheets and localized coherent structures, because 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 displays the subregion $(x,y) \in [130,160]d_i \times [140,170] d_i$ of the HPIC 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, $\boldsymbol {J} = \boldsymbol {\nabla }\times \boldsymbol {B}$, is shown in panels (a) and (f) respectively. There, we can distinguish the small plasmoids in the magnetic amplitude, whereas the corresponding signal in the current density is almost swamped by the PPC noise and by dissipation. The amplitudes of the aggregated IMFs (a–e) and of their associated current densities (f–j) are shown in the other panels. Now the three plasmoids are clearly visible in $\hat {\boldsymbol {B}}_\mathrm {kin}$ and $\hat {\boldsymbol {J}}_\mathrm {kin}$, and their large-scale signature also appears at MHD scales.

Figure 7. (a–e) Amplitude of (from left to right) magnetic field fluctuations and corresponding aggregated IMFs of a subregion containing a current sheet undergoing plasmoid reconnection between two vortices, in the HPIC simulation and at the time of maximum turbulent activity. (f–j) Same as (a–e), but for the current density.

As expected, the aggregated IMFs also reveal that magnetic ohmic dissipation is mostly concentrated in strong current structures (high values of $\hat {\boldsymbol {J}}_{\mathrm {diss}}$ are found in areas of high $\hat {\boldsymbol {J}}_\mathrm {kin}$). Moreover, the highest dissipation (bright spots of $\hat {\boldsymbol {J}}_{\mathrm {diss}}$) takes place at the X-points between the plasmoids, a typical feature of magnetic reconnection events. All these morphological and physical multiscale properties, typically 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_{\perp }^{(j)}$, the kurtosis

(4.4)\begin{equation} K(\,f_j^{>}) = \frac{\langle (\,f_j^{>})^4 \rangle} {\langle (\,f_j^{>})^2 \rangle^2}, \end{equation}

where $f_j^{>}(x,y)$ is a high-pass filtered signal of $f(x,y)$, only containing frequencies $k_{\perp } \geq k_{\perp }^{(j)}$, and where $\langle ~\rangle$ 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

(4.5)\begin{equation} f_j^{>}(x,y) = \sum_{k_{\perp}^{(i)} \geq k_{\perp}^{(j)}} \,\hat{f}_i(x,y), \end{equation}

by summing over all the IMFs with an average frequency $k_{\perp }^{(i)} \geq k_{\perp }^{(j)}$. We point out that IF-based methods proved to be well suited to reconstructing the kurtosis, and more generally the multiscale statistical properties of non-stationary signals (including intermittent ones), as recently shown in Stallone, Cicone & Materassi (2020).

Another quantity that measure the departure from a gaussian behaviour is the Kullback–Leibler (KL) divergence (e.g. Granero-Belinchón, Roux & Garnier 2018). For a sample $X$ with probability density function (PDF) $p(x)$ of variance $\sigma _X^2$, the KL-divergence is defined as

(4.6)\begin{equation} \mathcal{K}_L (X) = H_G(X) - H(X), \end{equation}

where

(4.7)\begin{equation} H (X)= -\int_\mathbb{R} p (x) \log p(x)\,\mathrm{d}x \end{equation}

is the Shannon entropy of the sample, and $H_G = 0.5\log (2{\rm \pi} e \sigma _X^2)$ is the Shannon entropy that $X$ would have if $p(x)$ were a Gaussian. Here $\mathcal {K}_L (X)$ is always positive, being zero if the PDF of the sample is a Gaussian distribution. Finally, the KL divergence $\mathcal {K}_L (\,f_j^{>})$ of $f(x,y)$ at a given frequency $k_{\perp }^{(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$ (a zero excess denotes a Gaussian behaviour), together with the KL divergence of the magnetic field components of the HMHD (a,c) and the HPIC (b,d) datasets. The results, in qualitative agreement with our previous findings (Papini et al. 2019b), show a flat kurtosis at scales above the injection scale ($k_{\perp } d_i \lesssim 0.28$), denoting a self-similar behaviour 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 has been 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_{\perp } d_i = 1$, then abruptly increases with a behaviour similar 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$, whereas in the kinetic regime, dispersive effects couple the fluctuations with $B_z$.

Figure 8. Excess kurtosis (a,b) and KL divergence (c,d) of the magnetic field fluctuations as calculated from (4.4) and (4.6) respectively, by using the IMFs of the HMHD run (a,c) and of the HPIC run (b,d).

The simulations are initialized 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 figure 5 of Papini et al. 2019b) shows 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 (figure 8c,d) shows the same results. A departure from Gaussian behaviour ($\mathcal {K}_L >0$) is observed in the perpendicular fluctuations right below the injection scale 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_{\perp } d_i<1$), are not zero, which denote a non-Gaussian nature of the fluctuations. We finally note that at the smallest scales, where turbulent dissipation kicks in, the KL divergence decreases in both the datasets, as expected.

5. Discussion

In this work, we have investigated the magnetic multiscale properties of both fluid HMHD and hybrid ion-kinetic electron-fluid simulations of plasma turbulence by means of MIF, a novel technique for the decomposition of non-stationary multidimensional signals. By exploiting our large-scale high-resolution numerical datasets, we successfully 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.

5.1. Intermittency

Our results confirm that plasma turbulence is an intrinsic multiscale phenomenon. In the inertial range, the energy cascade consists of more or less homogenously distributed and weakly intermittent magnetic fluctuations (a slowly increasing kurtosis in the perpendicular fluctuations is observed). Ion kinetic scales are characterized by strongly localized coherent structures organized in a filamented network, which show a high degree of intermittency. Finally, when reaching the turbulent 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^>$ containing all spatial frequencies $k_{\perp } \geq k_{\perp }^{(j)}$ (see (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 employed extensively 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 $\Delta B_y^{\ell ,y} = B_y(x,y+\ell ) - B_y(x,y)$ (green curve) and $\Delta B_y^{\ell ,x}=B_y(x+\ell ,y)-B_y(x,y)$ (red curve), respectively parallel and perpendicular (in the plane) to the direction of the $B_y$ component (where $\ell = 2{\rm \pi} /k_{\perp }$). We also report the kurtosis obtained from the MIF decomposition of $B_y$ (already shown in figure 8b) and that obtained using fast Fourier transform (FFT) 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 $\Delta B_y^{\ell ,x}$ is almost constant in the inertial range, then linearly increases starting from $k_{\perp } d_i \simeq 1$, i.e. at the scales where Hall currents become important. Instead, $K(\Delta B_y^{\ell ,y})$ linearly increases already at the injection scales, but with a smaller slope than $K(\Delta B_y^{\ell ,x})$. Such a behaviour may be explained by considering the following argument. First, the increments $\Delta B_y^{\ell ,y}\propto \partial _y B_y$ and $\Delta B_y^{\ell ,x}\propto \partial _x B_y$ are proportional to the components of the magnetic field gradient parallel and perpendicular to the magnetic field direction, respectively. Second, localized and elongated magnetic structures, such as the thin current sheets between vortices that we observed, are aligned with the magnetic field. Therefore, $K(\Delta B_y^{\ell ,x})$ is sensitive to the thickness of such structures (which is of the order of $d_i$) whereas $K(\Delta B_y^{\ell ,y})$ probes their length (that is comparable with the size of the largest vortices). This may explain the observed increase in the parallel and perpendicular kurtosis at around those scales.

Figure 9. Excess kurtosis of the $y$-component of the magnetic field fluctuations from the HPIC dataset, as calculated from the increments along the $x$-direction ($\Delta B_y^{\ell ,x}$) and along the $y$-direction ($\Delta B_y^{\ell ,y}$). The excess kurtosis as given by (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 $B_y$, is also shown.

The kurtosis calculated with MIF, although it is roughly twice as large, reproduces the properties of both $K(\Delta B_y^{\ell ,y})$ and $K(\Delta B_y^{\ell ,x})$, following the former at the large scales down to the spectral break $k_{\perp } d_i \simeq 2$ and then steepening at the scales where $K(\Delta B_y^{\ell ,x}) > K(\Delta B_y^{\ell ,y})$. Such differences mainly arise from the fact that the kurtosis calculated from the increments at a given separation scale $\ell$ also depends on the fluctuations with wavenumber $k_{\perp } < 2{\rm \pi} /\ell$. Indeed, the kurtosis of the high-pass filtered increments of $\Delta B_y^{\ell ,x}$ (not shown here), almost overlaps with that obtained by using either the MIF or the FFT decomposition.

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_{\perp }$ increases). Analogous results are obtained for the $x$-component of the magnetic field. As already mentioned in § 4.1, the excess kurtosis (and the KL divergence) of $B_z$ from the HPIC dataset, although constant, is not zero at fluid scales ($k_{\perp } d_i<1$). The HMHD run, in contrast, shows an almost zero excess kurtosis for $B_z$ at those scales. At present, we do not have an explanation for such difference. However, it is reasonable to assume that the underlying physical mechanism is the same that causes the differences observed between the HMHD and the HPIC power spectra of $B_z$ (Papini et al. 2019b).

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 behaviour at dissipation scales. There is, however, an important difference. In the ion kinetic range, where AL19 observe a non-intermittent behaviour, we find high levels of intermittency. 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 $\sim -5/2$, different from that we report in our simulations ($-3$) and also from other solar wind conditions. Further investigation pursuing this path is currently underway, in order to measure the multifractal properties of our numerical simulations.

5.2. Reconnection and enhanced turbulent dissipation

Our multiscale analysis confirms that magnetic field turbulent 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 (Yang et al. 2017; Camporeale et al. 2018; Kuzzay, Alexandrova & Matteini 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, Kida & Ohkitani 1994), and consistent with the fact that areas of both high magnetic gradients and vorticities are correlated (Franci et al. 2016a; Kuzzay et al. 2019), especially at the reconnection sites, which produce high levels of vorticity (e.g. Widmer, Büchner & Yokoi 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 particle-per-cell noise and dissipation (see figure 7). In this context, MIF 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 a result of the presence of structures, as the results of this work seems to suggest. In a future work, we will exploit the multiscale capabilities of MIF methods to investigate the local properties of the plasma in the filamented network, so to further characterize turbulent dissipation.

Overall, MIF 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.