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. Reference Bruno, Carbone, Vörös, D'Amicis, Bavassano, Cattaneo, Mura, Milillo, Orsini and Veltri2009; Chen et al. Reference Chen, Bale, Salem and Maruca2013; Kiyani, Osman & Chapman Reference Kiyani, Osman and Chapman2015) and in the Earth's magnetosheath (Stawarz et al. Reference Stawarz, Eriksson, Wilder, Ergun, Schwartz, Pouquet, Burch, Giles, Khotyaintsev and Le Contel2016; Chen & Boldyrev Reference Chen and Boldyrev2017) 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 Reference Iroshnikov1963; Kraichnan Reference Kraichnan1965). 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 Reference Verscharen, Klein and Maruca2019, 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. Reference Howes, Cowley, Dorland, Hammett, Quataert and Schekochihin2008; Schekochihin et al. Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Quataert and Tatsuno2009), eventually complemented by kinetic dissipative effects, such as Landau damping (Sulem et al. Reference Sulem, Passot, Laveder and Borgogno2016) and other effects associated with the deformation of the particles velocity distribution functions (Yang et al. Reference Yang, Matthaeus, Parashar, Haggerty, Roytershteyn, Daughton, Wan, Shi and Chen2017; Del Sarto & Pegoraro Reference Del Sarto and Pegoraro2018). 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. Reference Hellinger, Verdini, Landi, Franci and Matteini2018), 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 Reference Frisch1995; Marsch & Tu Reference Marsch and Tu1997), 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. Reference Camporeale, Sorriso-Valvo, Califano and Retinò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 Reference Boldyrev and Perez2012; Loureiro & Boldyrev Reference Loureiro and Boldyrev2017; Mallet, Schekochihin & Chandran Reference Mallet, Schekochihin and Chandran2017; Landi et al. Reference Landi, Franci, Papini, Verdini, Matteini and Hellinger2019).

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. Reference Howes, Tenbarge, Dorland, Quataert, Schekochihin, Numata and Tatsuno2011; Servidio et al. Reference Servidio, Valentini, Califano and Veltri2012; Wan et al. Reference Wan, Matthaeus, Roytershteyn, Karimabadi, Parashar, Wu and Shay2015; Haggerty et al. Reference Haggerty, Parashar, Matthaeus, Shay, Yang, Wan, Wu and Servidio2017). 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. Reference Franci, Landi, Matteini, Verdini and Hellinger2015a,Reference Franci, Verdini, Matteini, Landi and Hellingerb, Reference Franci, Hellinger, Guarrasi, Chen, Papini, Verdini, Matteini and Landi2018a,Reference Franci, Landi, Verdini, Matteini and Hellingerb, Reference Franci, Stawarz, Papini, Hellinger, Nakamura, Burgess, Landi, Verdini, Matteini and Ergun2020; Papini et al. Reference Papini, Franci, Landi, Verdini, Matteini and Hellinger2019b) and in the Earth magnetosheath (Franci et al. Reference Franci, Stawarz, Papini, Hellinger, Nakamura, Burgess, Landi, Verdini, Matteini and Ergun2020). 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. Reference Franci, Landi, Matteini, Verdini and Hellinger2016b; Cerri & Califano Reference Cerri and Califano2017; Franci et al. Reference Franci, Cerri, Califano, Landi, Papini, Verdini, Matteini, Jenko and Hellinger2017). Moreover, spectral properties at the reconnection exhausts consistent with a developed turbulent state were observed in a fully kinetic simulation (Pucci et al. Reference Pucci, Servidio, Sorriso-Valvo, Olshevsky, Matthaeus, Malara, Goldman, Newman and Lapenta2017). Finally, recent works (Franci et al. Reference Franci, Cerri, Califano, Landi, Papini, Verdini, Matteini, Jenko and Hellinger2017; Papini et al. Reference Papini, Franci, Landi, Hellinger, Verdini and Matteini2019a,Reference Papini, Franci, Landi, Verdini, Matteini and Hellingerb) 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. Reference Bruno, Carbone, Veltri, Pietropaolo and Bavassano2001; Chang, Tam & Wu Reference Chang, Tam and Wu2004; Consolini, Chang & Lui Reference Consolini, Chang and Lui2005; Horbury, Forman & Oughton Reference Horbury, Forman and Oughton2008; Lion, Alexandrova & Zaslavsky Reference Lion, Alexandrova and Zaslavsky2016). 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. (Reference Huang, Shen, Long, Wu, Shih, Zheng, Yen, Tung and Liu1998) 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. Reference Alberti, Consolini, Carbone, Yordanova, Marcucci and De Michelis2019). Unfortunately, EMD shown to be unstable in presence of noise, and the ensemble EMD (EEMD, Wu & Huang Reference Wu and Huang2009) 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 Reference Lin, Wang and Zhou2009; Cicone, Liu & Zhou Reference Cicone, Liu and Zhou2016; Cicone Reference Cicone2020). 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. Reference Bertello, Piersanti, Candidi, Diego and Ubertini2018; Piersanti et al. Reference Piersanti, Materassi, Cicone, Spogli, Zhou and Ezquer2018; Spogli et al. Reference Spogli, Piersanti, Cesaroni, Materassi, Cicone, Alfonsi, Romano and Ezquer2019). 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 Reference Cicone and Zhou2017, Reference Cicone and Zhou2020; Cicone Reference Cicone2020).

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. (Reference Papini, Franci, Landi, Verdini, Matteini and Hellinger2019b).

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. Reference Shay, Drake, Rogers and Denton2001). 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. Reference Landi, Del Zanna, Papini, Pucci and Velli2015; Papini, Landi & Del Zanna Reference Papini, Landi and Del Zanna2018, Reference Papini, Landi and Del Zanna2019c) and plasma turbulence (Papini et al. Reference Papini, Franci, Landi, Hellinger, Verdini and Matteini2019a,Reference Papini, Franci, Landi, Verdini, Matteini and Hellingerb). 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 Reference Wray1990).

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 Reference Mininni and Pouquet2009), 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. Reference Papini, Franci, Landi, Verdini, Matteini and Hellinger2019b) (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 (Reference Cicone and Zhou2017).

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. Reference Cicone, Liu and Zhou2016; Cicone & Zhou Reference Cicone and Zhou2017; Cicone & Dell'Acqua Reference Cicone and Dell'Acqua2020).

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. Reference Lin, Wang and Zhou2009; Cicone et al. Reference Cicone, Liu and Zhou2016). 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. Reference Cicone, Liu and Zhou2016). 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 (Reference Frisch1995), 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 (Reference Stallone, Cicone and Materassi2020).

Another quantity that measure the departure from a gaussian behaviour is the Kullback–Leibler (KL) divergence (e.g. Granero-Belinchón, Roux & Garnier Reference Granero-Belinchón, Roux and Garnier2018). 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. Reference Papini, Franci, Landi, Verdini, Matteini and Hellinger2019b), 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. Reference Papini, Franci, Landi, Verdini, Matteini and Hellinger2019b) 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.