# Can Hall Magnetohydrodynamics Explain Plasma Turbulence at Sub-ion Scales?

Emanuele Papini, Luca Franci, Simone Landi, Andrea Verdini, Lorenzo Matteini, Petr Hellinger

► To cite this version:

Emanuele Papini, Luca Franci, Simone Landi, Andrea Verdini, Lorenzo Matteini, et al. Can Hall Magnetohydrodynamics Explain Plasma Turbulence at Sub-ion Scales?. The Astrophysical Journal, 2019, 870 (1), pp.52. 10.3847/1538-4357/aaf003. hal-03111701

HAL Id: hal-03111701

https://hal.science/hal-03111701v1

Submitted on 13 Aug 2025

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Can Hall Magnetohydrodynamics Explain Plasma Turbulence at Sub-ion Scales?

Emanuele Papini¹, Luca Franci²,³, Simone Landi¹,³, Andrea Verdini¹, Lorenzo Matteini⁴, and Petr Hellinger⁵

¹Dipartimento di Fisica e Astronomia, Università degli Studi di Firenze, via G. Sansone 1, I-50019 Sesto Fiorentino, Italy; papini@arcetri.inaf.it

²School of Physics and Astronomy, Queen Mary University of London, London E1 4NS, UK

³INAF, Osservatorio Astrofisico of Arcetri, Largo E. Fermi 5, 50125, I-50125 Firenze, Italy

⁴LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Univ. Paris Diderot, Sorbonne Paris Cité, place Jules Janssens 5, 92195 F-92195 Meudon, France

⁵Astronomical Institute, CAS, Bocni II/1401, CZ-14100 Prague, Czech Republic

Received 2018 September 6; revised 2018 November 5; accepted 2018 November 9; published 2019 January 7

## Abstract

We investigate the properties of plasma turbulence by means of two-dimensional Hall-magnetohydrodynamic (HMHD) and hybrid particle-in-cell (HPIC) numerical simulations. We find that the HMHD simulations exhibit spectral properties that are in most cases in agreement with the results of the HPIC simulations and with solar wind observations. The energy spectra of magnetic fluctuations exhibit a double power law with spectral index -5/3 at MHD scales and -3 at kinetic scales, while for velocity fluctuations the spectral index is -3/2 at MHD scales. The break between the MHD and the kinetic scales occurs at the same scale in both simulations. In the MHD range the slopes of the total energy and residual energy spectra satisfy a fast Alfvén-dynamo balance. The development of a turbulent cascade is concurrently characterized by magnetic reconnection events taking place in thin current sheets that form between large eddies. A statistical analysis reveals that reconnection is qualitatively the same and fast in both the HMHD and HPIC models, characterized by inverse reconnection rates much smaller than the characteristic large-eddy nonlinear time. The agreement extends to other statistical properties, such as the kurtosis of the magnetic field. Moreover, the observation of a direct energy transfer from the large vortices to the small sub-ion scales, triggered by magnetic reconnection, further supports the existence of a reconnection-mediated turbulent regime at kinetic scales. We conclude that the HMHD fluid description captures to a large extent the transition of the turbulent cascade between the large MHD scales and the sub-ion scales.

**Key words:** magnetic reconnection – plasmas – solar wind – turbulence

## 1. Introduction

Understanding the mechanisms responsible for transferring magnetic and kinetic energy from large to small scales in weakly collisional plasmas is of great interest in space, astrophysical, and laboratory plasmas, since such mechanisms are relevant in many fundamental processes: the formation of hot coronae, the heating and acceleration of stellar winds and their interaction with planetary magnetospheres, the acceleration of particles, and explosive phenomena such as solar flares and coronal mass ejections, to name a few. Indeed, in hot rarefied plasmas collisions between particles are not efficient in dissipating energy at scales above the particle's characteristic kinetic scale, namely the ion (electron) Larmor and demagnetization scales. This is the case of, e.g., the solar wind and the Earth's magnetosheath. In situ spacecraft observations with high spatial and temporal resolution (see, e.g., Chen et al. 2013b; Chen & Boldyrev 2017) have shown that the magnetic energy spectrum follows a Kolmogorov power-law cascade at scales larger than the ion kinetic length scales. As the cascade approaches these scales, namely the ion inertial length $d_i$ or the ion gyroradius, the one-fluid magnetohydrodynamic (MHD) description breaks and nonlinear interactions are strongly modified by the different propagation and polarization properties introduced by the specific species' dynamics. For instance, the heavy ions have a different gyromotion around the magnetic field lines than the light electrons, and at scales below $d_i$ the motion of the electrons decouples from that of the ions (Hall effect).

Indeed, observations have shown that around such scales the magnetic field spectrum steepens and the nature of the

fluctuations changes (see, e.g., Denskat et al. 1983; Goldstein et al. 1994). In particular, a power-law spectrum with a slope between -2.8 and -3 is routinely observed (see, e.g., Alexandrova et al. 2009; Chen & Boldyrev 2017). Concurrently, the electric to magnetic field power increases (Salem et al. 2012), together with compressive fluctuations (see, e.g., Chen et al. 2012, 2013a; Matteini et al. 2017), and the nature of the fluctuations seems to be compatible with strongly obliquely Alfvénic modes, the so-called kinetic Alfvén waves (see, e.g., Salem et al. 2012; Chen et al. 2013a).

Early attempts to describe the transition of the electromagnetic turbulence across ion scales use fluid extended MHD models, either in the low (Hall-MHD, or HMHD; Ghosh et al. 1996; Galtier 2006; Galtier & Buchlin 2007; Shaikh & Shukla 2009) or in the high (electron MHD; Biskamp et al. 1996, 1999; Cho & Lazarian 2004, 2009; Meyrand & Galtier 2013) frequency approximation. In these approaches, we observe a transition near the ion scales with a steepening in the magnetic field power and a simultaneous flattening in the electric field power fluctuations (Matthaeus et al. 2008). Steeper spectra, as observed in the solar wind and in magnetospheres, have been often reproduced by the use of both 2D (Franci et al. 2015a, 2015b, 2016b) and 3D (see, e.g., Howes et al. 2011; Franci et al. 2018b) kinetic models, as well as in fluid models that include kinetic dissipative effects such as ion and electron Landau damping (Sulem et al. 2016). It has been thus suggested that the steepening of the spectra is strictly correlated to deformations in the particle's velocity distribution function associated, for example, to wave-particle interactions such as instabilities and wave damping, to thermal anisotropies, and/or to non-gyrotropic terms in the full pressure tensor
(Del Sarto et al. 2016; Del Sarto & Pegoraro 2018) which can drive dissipative effects (Yang et al. 2017a). Recently, high-resolution hybrid particle-in-cell (HPIC) simulations (Franci et al. 2015a, 2015b), which use a kinetic description for the ions and a fluid one for the electrons, were able to reproduce the steep spectra observed at kinetic scales in the solar wind and to show that the transition between the MHD and the ion scales is correlated to the largest characteristic kinetic scale encountered by the nonlinear cascade (Franci et al. 2016b). It was also shown that the formation of the turbulent cascade in the sub-ion range is associated to the formation of small-scale reconnecting current sheets, which rapidly transfer the magnetofluid energy at sub-ion scales (Cerri & Califano 2017; Franci et al. 2017). Finally, analyses of HPIC simulations based on the von Karman-Howarth equations (Galtier 2008; Hellinger et al. 2018) indicate that a part of the MHD cascade continues at sub-ion scales via the Hall term. The two cascades are not, however, separable and overlap (the Hall part progressively dominating at smaller scales), both contributing to the total HMHD cascade.

The main goal of this work is to check whether a fluid model that allows for the formation of small-scale coherent structures, such as current sheets, can reproduce spectra steeper than expected by simpler phenomenological models, without the need to include purely kinetic effects. To that purpose, we performed a 2D high-resolution HMHD simulation for a warm polytropic plasma and compared the outcome with analogous results obtained from an HPIC simulation.

Our results suggest that the main processes for the formation of steeper spectra involve the formation of coherent small-scale structures and are consistent with recent models of reconnection-mediated turbulence, such as proposed by Boldyrev & Perez (2012), Loureiro & Boldyrev (2017), and Mallet et al. (2017).

## 2. HMHD Simulations: Numerical Setup

Our model integrates the viscous-resistive MHD equations but retaining the Hall term in the induction equation, that is, by substituting the fluid velocity $v$ with the electron velocity $v_e = v - J/\epsilon n_e$. In their adimensionalized form, the HMHD equations read

$$ \partial_t \rho + \nabla \cdot (\rho v) = 0, \quad (1) $$

$$ \rho(\partial_t + v \cdot \nabla)v = -\nabla P + (\nabla \times B) \times B + \nu \left[ \Delta v + \frac{1}{3} \nabla(\nabla \cdot v) \right], \quad (2) $$

$$ (\partial_t + v \cdot \nabla)T = (\Gamma - 1) \left\{ - (\nabla \cdot v) T + \eta \frac{|\nabla \times B|^2}{\rho} + \frac{\nu}{\rho} \left[ (\nabla \times v)^2 + \frac{4}{3} (\nabla \cdot v)^2 \right] \right\}, \quad (3) $$

$$ \partial_t \mathbf{B} = \nabla \times (v \times \mathbf{B}) + \eta \nabla^2 \mathbf{B} - \eta_H \nabla \times \frac{(\nabla \times \mathbf{B}) \times \mathbf{B}}{\rho}, \quad (4) $$

where $\Gamma = 5/3$ is the adiabatic index and the variables retain their usual meaning. All quantities are renormalized using $L = d_i$ as characteristic length and with respect to a field amplitude $B_0$, a density $\rho_0$, an Alfvén speed $c_A = B_0/\sqrt{4\pi\rho_0} = \Omega_i d_i$, a pressure

$P_0 = \rho_0 c_A^2$, and a temperature $T_0 = (k_B/m_i) P_0/\rho_0$. Here, 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. $\Omega_i = eB_0/m_i c$ is the ion-cyclotron angular frequency and $m_i$ is the mass of the ions.

Equations (1)–(4) are numerically solved using a code we employed for studies of magnetic reconnection (Landi et al. 2015; Papini et al. 2018), modified to include periodic boundaries in all directions. We consider a 2D $(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 (Orszag 1971), to avoid aliasing of the nonlinear quadratic terms. Aliasing of the cubic terms is mitigated by the presence of a finite dissipation (Ghosh et al. 1993). Time integration is performed via a third-order Runge-Kutta scheme. The other code employed in this work is the Lagrangian HPIC code CAMELIA (Matthews 1994; Franci et al. 2018a), which integrates the Vlasov-Maxwell equations by coupling fully kinetic ions to a charge-neutralizing and massless fluid of isothermal electrons. CAMELIA has been successfully used for numerical studies of plasma turbulence (Franci et al. 2015a, 2015b, 2016a, 2016b, 2017), and it reproduced many of the spectral properties observed in the solar wind (Franci et al. 2018a) and in the Earth’s magnetosheath by the Magnetospheric Multiscale (MMS) mission (L. Franci et al. 2019, in preparation).

The numerical setup used in this work is the same as in the HPIC simulation of Franci et al. (2017), in order to allow a straightforward comparison between HMHD and HPIC results. We consider a 2D 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² points. Out of the plane, along the z-direction (identified as the parallel direction), we set a background constant magnetic field $B_0$, which we refer to as the mean field. The initial state is populated by freely decaying large-amplitude Alfvénic-like fluctuations with zero mean cross helicity in the xy-plane perpendicular to the mean field (Franci et al. 2015a) and up to the injection scale $\ell_{\text{inj}} = 2\pi d_i / k_{\perp}^{\text{inj}}$, with $k_{\perp}^{\text{inj}} d_i \approx 0.2$, where $k_{\perp} = \sqrt{k_x^2 + k_y^2}$. The root-mean-square (rms) amplitude of these fluctuations is set to $b_{\text{rms}} = B_{\text{rms}}/B_0 \approx 0.24$. In the HPIC reference simulation, the plasma beta for ions and electrons and the magnetic resistivity are set to the values $\beta_i = \beta_e = 1$ and $\eta = 5 \times 10^{-4}$ respectively. In HMHD we set the plasma $\beta = \beta_i + \beta_e = 2$ accordingly. Instead, for the resistivity, we employ $\eta = 10^{-3}$, i.e., twice the value of the hybrid simulation. We also set the viscosity to the same value. At the global scales of the box, these values correspond to a viscous and magnetic Reynolds number of $R = R_m \approx 40,000$.

## 3. Role of Magnetic Reconnection in Developing Turbulence

The initial Alfvénic fluctuations quickly evolve to form vortices and localized current sheets, which then shrink down to a critical width of the order of $d_i$ in a time of about $30 \Omega_i^{-1}$. These current sheets quickly disrupt, due to the onset of fast magnetic reconnection processes, generating a variety of small-scale magnetic islands and fluctuations that are fed back into the turbulent surroundings (a similar dynamics is known to occur also in MHD turbulence; see, e.g., Matthaeus & Lamkin 1986; Biskamp & Welter 1989; Carbone et al. 1990). As an example, Figure 1 shows a contour plot of the out-of-plane current density, $J_z$, at $t = 45 \Omega_i^{-1}$, right after the
<figure><img src="image_3.png" /><figcaption>Figure 1. Colored contours of the out-of-plane current density $J_\parallel$ for the HMHD and the HPIC simulations at $t = 45 \Omega_i^{-1}$, in the whole grid (large panels on the left) and in a subgrid containing one reconnecting current sheet (small panels on the right). Color scales are saturated to ±0.5 for easy visualization of the structures.</figcaption></figure>

first reconnection events have been observed, for both the HPIC and the HMHD runs. In particular, a zoom on two current sheets (small panels) reveals a plasmoid-chain structure, a characteristic footprint of the tearing instability (see, e.g., Loureiro et al. 2007; Papini et al. 2018). As the system evolves, other reconnection events take place in newly formed current sheets until turbulence is fully developed.

To quantitatively support the above description, we analyzed in detail the reconnection events in both simulations. A reconnection event is identified by a magnetic X-point and its nearest O-point in a current sheet. For each event we calculated the reconnection rate

$$ \gamma_{\text{rec}} = \left| \frac{1}{\Phi|_{X}^{\text{O}}} \frac{\partial \Phi|_{X}^{\text{O}}}{\partial t} \right|, \qquad (5) $$

where $\Phi|_{X}^{\text{O}}$ is the reconnected magnetic flux density between the X-point and the O-point, i.e., the difference in the out-of-plane vector potential between those points, $\Phi|_{X}^{\text{O}} = A_{z}^{\text{O}} - A_{z}^{\text{X}}$ (for further details, see the Appendix). Figure 2 shows the temporal evolution of the distribution of reconnection rates for the HMHD run (a) and the HPIC run (b), together with the maximum of $|J|$ (c) and the rms of the current density $J$ (d). For the HMHD run, the first reconnection events are detected at $t \approx 30 \Omega_i^{-1}$ and concurrently with the first maximum of $|J|$ (marked by a vertical dashed line), which is a proxy for reconnection events (Franci et al. 2017). The reconnection rates follow a log-normal distribution, with a mean value of $\langle \gamma_{\text{rec}} \rangle = 0.08 \Omega_i$ (red solid line). Figure 3 shows the distribution of all the reconnection rates measured in the HMHD run at all times. As the system evolves, this value remains roughly constant and corresponds to an average reconnection time $\langle \tau_{\text{rec}} \rangle = 1/\langle \gamma_{\text{rec}} \rangle \approx 12.2 \Omega_i^{-1}$. We note, however, that a significant fraction (roughly 22%) of the reconnection events has a reconnection time smaller than the average, and 10% of the reconnection events are very fast, with $\tau_{\text{rec}} < 6 \Omega_i^{-1}$. After the maximum of turbulent activity is reached, i.e., after the maximum of rms($|J|$) (Mininni & Pouquet 2009) at $t \approx 165 \Omega_i^{-1}$, the number of reconnection events decreases.

Magnetic reconnection modifies the turbulent properties at kinetic scales soon after the first reconnection events are

<figure><img src="image_4.png" /><figcaption>Figure 2. Time evolution of the reconnection rates $\gamma_{\text{rec}}/\Omega_i$ for the HMHD (a) and the HPIC run (b): colored contours indicate the distribution of reconnection events, horizontal solid red lines denote the average reconnection rate, and the dashed red lines denote 10th and 90th percentile values. (c) Time evolution of the maximum of $|J|$, and (d) of the rms value of $J$. In all plots, the vertical dashed and dotted-dashed lines mark the first maximum in $|J|$ of the HMHD run and the HPIC run respectively.</figcaption></figure>

detected. This is made evident by looking at the nonlinear time associated to turbulence, that is the eddy turnover time. Such a characteristic time at a given scale is estimated as $\tau_{nl}(k_\perp) = (k_\perp v_e(k_\perp))^{-1}$, where $v_e(k_\perp)$ is the electron velocity at
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Histogram of the reconnection rates $\gamma_{\text{rec}}$ as measured in the HMHD run at all times. The blue curve is the log-normal fit to the histogram, corresponding to a Gaussian fit of $\log(\gamma_{\text{rec}})$. As in Figure 2(a), the red solid line denotes the mean reconnection rate, and the dashed red lines denote 10th and 90th percentile values. The black solid line denotes the mean of $\log(\gamma_{\text{rec}})$.</figcaption></figure>

the scale $\ell = 2\pi/k_{\perp}$. The use of this definition unifies the classic definition at MHD scales (where the electron velocity equals the fluid velocity) and the definition of the nonlinear time at kinetic scales. Figure 4 shows a contour plot of the nonlinear time as a function of time and wavenumber. Here one can identify three distinct phases: an early evolutionary phase ($0 \le t\Omega_i \lesssim 30$) in which the initial configuration relaxes and the first vortices and current sheets are formed, a transient phase ($30 \lesssim t\Omega_i \lesssim 50$) triggered by the first reconnection events, and a third phase ($50 \lesssim t\Omega_i$) characterized by slowly evolving values of $\tau_{nl}(k_{\perp})$ that lasts until the end of the simulation. In the first phase, at the beginning, the scale with the smallest value of $\tau_{nl}(k_{\perp})$ is the injection scale $\ell_{\text{inj}}$ (marked by the horizontal dotted line), with $\tau_{nl}(k_{\text{inj}}) \simeq 22 \Omega_i^{-1}$. As the system evolves and the energy gets redistributed, the nonlinear time at scales above the injection scale increases, while at the scales right below the injection scale the nonlinear time slightly decreases, due to the development of a direct cascade. At small scales and before the first reconnection events are detected ($t \lesssim 30 \Omega_i^{-1}$), the nonlinear time is very large, since there is not yet a relevant amount of energy at those scales. As soon as reconnection is triggered, the second phase begins. The energy is directly transferred to the smallest scales accessible (the $2/3$ cutoff scale), where $\tau_{nl}(k_{\perp})$ suddenly decreases. Then, during a transient phase that lasts between $t \simeq 30 \Omega_i^{-1}$ and $t \simeq 50 \Omega_i^{-1}$ (characterized by almost vertical isocontour lines of $\tau_{nl}(k_{\perp})$; see, e.g., the red contour line), the energy is fed to larger and larger scales. We interpret this behavior as the signature in Fourier space of the coalescence of magnetic islands (Finn & Kaw 1977), which we also observe in real space at the reconnection sites. During this transient, the number of reconnection events increases very rapidly and reaches a statistically stationary value at about $t = 50 \Omega_i^{-1}$. After this transient, the nonlinear time $\tau_{nl}(k_{\perp})$ changes only slightly and then becomes roughly constant at $t \gtrsim 150 \Omega_i^{-1}$, indicating that a stationary turbulent state has been reached (see Section 4). Figure 4 provides a direct quantitative evidence of the ability of magnetic reconnection to influence the dynamics of turbulence at kinetic scales.

The evolution in the HPIC simulation is qualitatively the same, though the first reconnection events are triggered at a

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Colored contours of the characteristic nonlinear time $\tau_{nl}(k_{\perp})$ as a function of wavenumber $k_{\perp}$ (scale $\ell$) and time $t$. The right axis denotes the corresponding scale $\ell = 2\pi/k_{\perp}$. The red contour line highlights the scale at which the nonlinear time equals the average reconnection time. The white area denotes the region where the nonlinear time exceeds $1000 \Omega_i^{-1}$. The horizontal dotted line indicates the injection scale.</figcaption></figure>

slightly later time ($t \simeq 40 \Omega_i^{-1}$) and the distribution of the reconnection rates is broader, with an average reconnection rate $\langle\gamma_{\text{rec}}\rangle = 0.22 \Omega_i$, roughly three times that of the HMHD run. We believe that nongyrotropic ion effects decrease the amplitude of the out-of-plane ion velocity $v_{i,\|}$ at the reconnection sites in the HPIC simulation, thus increasing the reconnection rate, as demonstrated by Yin et al. (2002). Indeed, in the current sheets of Figure 1 we measured an amplitude of $v_{i,\|}$ at the X-points in the HPIC run which is roughly 1/3 of that in the HMHD run.

Note that the maximum of turbulent activity in the HPIC run is reached later, at about $t \simeq 200 \Omega_i^{-1}$, due to the smaller value of the resistivity employed, which sets a dissipation scale smaller than that in the HMHD run, thus increasing the time needed to fully develop a turbulent cascade.

## 4. Intermittency and Spectral Properties in Fully Developed Turbulence

From the previous section we conclude that, despite the differences in the numerical approaches and theoretical models, the HMHD and HPIC simulations have remarkable similarities. The most interesting agreement is found in the spectral properties.

Figure 5 shows the omnidirectional power spectra of magnetic and velocity fluctuations of both the HPIC and the HMHD run, at the time when turbulence has fully developed and the rms of $J$ has reached its maximum. We recall that the fluid velocity $v$ in the HMHD description corresponds to the ion bulk velocity $v_i$ of the HPIC model. After the maximum, the spectra remain quite stable. At MHD scales, we observe a Kolmogorov-like cascade with slope $-5/3$ in the magnetic field fluctuations (see Figure 5(a)), while the power spectrum of the fluid (ion) velocity field fluctuations is flatter, with a spectral index of about $-3/2$ (see Figure 5(b)). This is in agreement with solar wind observations (Chen et al. 2011).

A transition occurs at ion kinetic scales, where the magnetic field steepens following a power law of $-3$. The HPIC and HMHD power spectra of total magnetic field fluctuations (Figure 5(a)) almost perfectly match at all scales, from the
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Isotropized power spectra $P(k_\perp)$ as a function of $k_\perp = \sqrt{k_x^2 + k_y^2}$ from the HMHD run (solid curves) at $t = 165 \Omega_i^{-1}$ and from the HPIC run (dashed curves) at $t = 200 \Omega_i^{-1}$, of total magnetic field (a) and perpendicular ion velocity (b) fluctuations, and of parallel magnetic field and parallel velocity fluctuations (c). The vertical dotted, dashed, and dotted-dashed lines denote the injection scale $k_\perp^{\text{inj}} d_i \simeq 0.2$, the 2/3 filter's cutoff of the HMHD model, and the Nyquist wavenumber respectively. Black dashed lines denote reference slopes. Panel (d) shows the residual energy spectrum $E_R = P_B - P_v$, compensated by $k_\perp^2$, from the HMHD simulation.</figcaption></figure>

injection scale (vertical dotted line) down to the scale where the HMHD spectrum has the filter's cutoff at 2/3 of the Nyquist frequency (vertical dashed line) and the HPIC spectrum rises due to numerical noise. The location of the spectral break at $k_\perp d_i \simeq 2$ also matches. The agreement extends to the spectra of the perpendicular magnetic fluctuations (not shown here).

For the spectra of parallel magnetic fluctuations (Figure 5(c)) we find some compelling differences. At kinetic scales, the parallel magnetic field spectra have the same quantitative and qualitative behavior, though the HPIC spectrum of $B_\parallel$ is

compatible with a power law of spectral index $-2.8$, while the HMHD spectrum of $B_\parallel$ is more compatible with a value of $-3$. At MHD scales the two spectra are different. There the HMHD spectrum of $B_\parallel$ is coupled to the spectrum of parallel velocity fluctuations (blue solid curve), thus suggesting that parallel Alfvénic fluctuations give the dominant contribution. In the HPIC case (dashed curves) such a coupling is not present. This picture is reversed in the spectra of the perpendicular velocity fluctuations (Figure 5(b)): at MHD scales and for both the HMHD and the HPIC simulations we recover a power-law cascade of spectral index $-3/2$. At kinetic scales ($k_\perp d_i > 1$) the two spectra diverge, due to the inability of HMHD to model ion kinetic effects.

In the MHD range of the HMHD run ($0.3 \lesssim k_\perp d_i \lesssim 1$) we obtain a residual energy spectrum $E_R = P_B - P_v$ consistent with a slope of $-2$ (Figure 4(d)) and a total energy spectrum $E_T = P_B + P_v \propto k_\perp^{-3/2}$ (Figure 4(e)). While the former agrees with the HPIC results (Franci et al. 2015b), the higher energy in the parallel component of velocity fluctuations found in the HMHD run causes its total energy spectrum to be flatter than $-5/3$. Müller & Grappin (2005) and Grappin et al. (2016) proposed that the residual energy spectrum and the total energy spectrum are related by a balance between a local dynamo effect and an Alfvénic coupling. Their relation reads:

$$E_R/E_T \approx A(t_A/t_{nl})^\alpha \quad (6)$$

where the Alfvén time $t_A = 1/(k_\perp b_{rms})$ is built on the large-scale magnetic fluctuations, $t_{nl} = 1/(k_\perp \sqrt{k_\perp E_T})$ is calculated using the total energy spectrum, the exponent on the r.h.s indicates a fast ($\alpha = 1$) or slow ($\alpha = 2$) Alfvénic coupling, and A is a constant of order 1. The above spectral slopes at the peak of the turbulent activity suggest a different scenario for the HPIC ($\alpha = 1$) and for the HMHD ($\alpha = 2$) runs. By comparing the l.h.s. and the r.h.s. of Equation (6) at different times, we found that in HPIC the exponent $\alpha = 1$ reproduces the slope of the ratio $E_R/E_T$ at all times past the peak of turbulent activity. In contrast, in HMHD, the exponent increases steadily from 2 to 3. This is due to the behavior of the parallel components in HMHD. In fact, the exponent matching the slope of the ratio $E_R/E_T$ is stable also in HMHD when only the perpendicular components of magnetic and velocity fluctuations are used, and in the following we will restrict our analysis to those components (i.e., we redefine $E_R = P_B - P_v$ and $E_T = P_B + P_v$). In Figure 6, the l.h.s. of Equation (6), averaged for about 1.5 large-eddy turnover times, is shown for the HPIC (top) and the HMHD (bottom) as thick solid lines, along with the average of the r.h.s. with $\alpha = 1,2$ (thin solid and dashed lines, respectively). The thick and thin solid lines are parallel in the MHD range ($0.3 \lesssim k_\perp d_i \lesssim 1$) for both simulations, thus supporting a fast scenario ($\alpha = 1$), while the slow scenario ($\alpha = 2$, dashed lines) has a steeper power law. Note that, although the ratio $E_R/E_T$ has a different power-law index (the ratio is compensated with different indexes in the two figures), the coefficients $A \approx 3$ and $\alpha = 1$ are the same in both simulations.

Finally, we focus on the intermittent properties of the two simulations. Localized current sheets and other coherent structures are intimately related to intermittency, that is the departure from Kolmogorov self-similarity law (Kolmogorov 1941). Among other properties, intermittency manifests itself
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Ratio of the residual energy to the total energy (thick solid lines) in the perpendicular (2D) components for the HPIC (top panel) and HMHD (bottom panel) simulations, averaged over 60 times after the peak of current density. The ratio of the characteristic Alfvén time and dynamo time, $(t_A/t_{nt})^\alpha$, is plotted as thin solid and dashed lines for the fast and slow scenario, respectively ($\alpha = 1$ and $\alpha = 2$; see Equation (6)).</figcaption></figure>

as a non-Gaussian behavior in the probability distribution functions (PDFs) of the increments of any field at a given scale $\ell$. In Figure 7(a)–(c) we report the PDFs of the increments $\Delta B_y(x, y) = B_y(x + \ell, y) - \Delta B_y(x, y)$ of the y-component of the magnetic field fluctuations along the x-direction and at three different spatial separations $\ell/d_i = 5\pi$, $\pi$, and $\pi/4$, i.e., in the inertial range (a), at the spectral break (b), and at kinetic scales (c) respectively. The PDFs are calculated at the same time $t = 200 \Omega_i^{-1}$ in both simulations. Indeed, the PDFs from the two models are comparable to each other at all scales. They are consistent with a Gaussian function at large scales in the inertial range ($\ell \gg d_i$), but as scales get smaller (e.g., already at $\ell = \pi d_i$ in Figure 7(b)), they depart from Gaussian behavior and develop fatter and fatter tails. To better quantify the level of intermittency, we calculated the scale-dependent excess kurtosis, shown in Figure 7 for the increments of the magnetic field (d) and of the velocity fields (e). The results for the magnetic field kurtosis are in remarkable agreement, being close to zero at large scales down to $2\pi d_i/\ell \approx 0.8$ and then increasing linearly down to the smallest scales where HPIC has a slightly larger kurtosis. The excess kurtosis of the velocity field increments are also similar at large scales, but the HPIC kurtosis becomes slightly larger already at intermediate scales, $2\pi d_i/\ell \gtrsim 0.3$, well before the velocity spectra start to diverge ($k_\perp d_i \gtrsim 2$, see Figure 4 (b)).

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Top panels: PDFs of the increments $\Delta B_y = B_y(x + \ell, y) - \Delta B_y(x, y)$ at three different scales $\ell$: in the inertial range (a), at the spectral break (b), and at sub-ion scales (c). The dotted lines drawing a triangle correspond to a Gaussian function. The PDFs are plotted against $\Delta B_y|\Delta B_y|/\sigma^2$, with $\sigma^2$ being the variance of each PDF. Bottom panels: excess kurtosis $K - 3$ for the magnetic field (d) and the velocity field (e), at $t = 200 \Omega_i^{-1}$ for both the HPIC run (dashed lines) and the HMHD run (solid lines). The horizontal dotted line denotes the zero excess kurtosis of a Gaussian function.</figcaption></figure>
## 5. Discussion

In this work, we provided numerical evidence that many of the statistical and spectral properties of plasma turbulence at sub-ion scales can be explained within the framework of HMHD. The results have been obtained by performing a study of turbulence generated by freely decaying Alfvénic-like fluctuations with zero mean cross helicity using both a full viscoresistive fluid HMHD model and a Vlasov-Maxwell HPIC model coupling fully kinetic ions to fluid and massless electrons, which was able to correctly reproduce in situ observations (Franci et al. 2018a). In the plasma regime we investigated, the HMHD and the hybrid model showed a remarkable agreement.

In the early evolution of our simulations and long before the formation of a power-law direct cascade at fluid scales, we observed, concurrently to the trigger of the first reconnection events, a transfer of energy from the large-scale vortices directly to small kinetic scales. After that, in a short transition phase, energy is fed into the whole kinetic range, which we interpret as the signature in Fourier space of the coalescence of magnetic islands and their expulsion from the first current sheets into the turbulent surroundings. Later, after a few eddy turnover times, once turbulence is almost developed, we are not able to separate the contribution of magnetic reconnection to the power spectrum. However, it is likely that such a reconnection-mediated energy transfer occurs whenever a current sheet reconnects, as suggested by Franci et al. (2017) and corroborated by recent studies that employ advanced techniques to measure the scale-to-scale energy transfer in localized coherent structures (e.g., Yang et al. 2017b; Camporeale et al. 2018). Our findings further strengthen the view that magnetic reconnection provides a direct channel to drive turbulence at kinetic scales (Franci et al. 2017), and thus supports theoretical models of reconnection-mediated turbulence (Loureiro & Boldyrev 2017; Mallet et al. 2017).

It is useful to compare our results on the reconnection rates with similar studies encompassing MHD (Servidio et al. 2010, hereafter SV10) and fully kinetic PIC (Haggerty et al. 2017, hereafter HG17) simulations of turbulence. In particular, our findings that the reconnection rates follow a log-normal distribution seem to contrast with the results of SV10 and HG17, who observed different distributions (see Figure 7 of SV10 and Figure 11 of HG17). To make a meaningful comparison with their work, we measured the values of the reconnecting out-of-plane electric field $E_X = \partial_x A_z$ at the X-points (that is the definition of reconnection rate used by both SV10 and HG17) and renormalized the values to the global Alfvén time, defined on the global scale $L_0 = 256d_i/2\pi$, as done in the above works. Indeed, the PDF of $|E_X|$ drawn in Figure 8 using a linear binning with roughly the same binsize of SV10 is both qualitatively and quantitatively similar to the PDFs of SV10 and HG17.

Once turbulence has fully developed, both the HMHD and the hybrid model reproduce a Kolmogorov-like cascade of spectral index $-5/3$ at fluid scales and a kinetic cascade of spectral index $-3$ in the total magnetic field spectra. The location of the spectral break is also recovered. A flatter power-law spectrum of slope $-3/2$ is observed in the velocity spectra at fluid scales, which is also responsible for the different slope in the total energy spectrum ($-5/3$ and $-3/2$ in HPIC and HMHD, respectively). We also obtain a residual energy spectrum of slope $\approx -2$. These indices are not consistent with

<figure><img src="image_10.png" /><figcaption>**Figure 8.** PDF of $|E_X|$ at the X-points, normalized with respect to the global Alfvén time $\tau_A = L_0/c_A$, with $L_0 = 256d_i/2\pi$. Solid red lines and filled circles denote a linear binning (constant $\Delta|E_X| = 0.01$) of the PDF. The vertical dashed line denotes the mean reconnection rate $\sim 0.023\tau_A^{-1}$ of the distribution.</figcaption></figure>

the same Alfvén-dynamo balance regulating the ratio of residual and total energy at inertial range scales (Müller & Grappin 2005; Grappin et al. 2016). However, when only perpendicular components are retained (i.e., 2D wavevectors and components are analyzed), both HPIC and HMHD display properties consistent with a fast Alfvén-dynamo balance scenario ($\alpha = 1$).

HPIC and HMHD simulations show a different behavior of the velocity field spectrum at ion scales, the latter being flatter and with a stronger Alfvénic coupling in the transition from MHD to ion scales. The overall agreement extends to other statistical properties, related to intermittency, and further confirms that the HMHD model may be accurate enough to describe the dynamics of magnetic fields fluctuations in a turbulent plasma.

Results of this work suggest that the Hall term is the dominant term shaping the turbulent properties and dynamics of the magnetic fields at sub-ion scales. That rules out other mechanisms not described by our HMHD model, for instance gyrotropic effects, temperature anisotropies, particle-wave interaction effects such as Landau damping, cyclotron resonance, and linear Vlasov instabilities. Some of these effects are likely responsible for the differences arising between the HMHD and HPIC models in the power spectra of the velocity fluctuations at kinetic scales. In particular, at large scales, the difference between the HMHD and HPIC power spectra of parallel magnetic field and velocity fluctuations may be due to nonzero off-diagonal terms in the ion pressure tensor generated by non-gyrotropic effects at the reconnection sites (Yin et al. 2002), which may also be responsible for the differences between the reconnection rates of the Hall and the hybrid model.

Recently, Gonzalez et al. investigated turbulence at kinetic scales by means of both a two-fluid incompressible HMHD model retaining electron inertia effects and a fully kinetic plasma model. Although the use of a proton-to-electron mass ratio of 25 in those models hinders a direct comparison with our results in the sub-ion range, they also find an agreement in spectral and turbulent properties between the two models, even at electron scales. Their results further confirm the potential use
of HMHD fluid models in the study of many properties of turbulence at kinetic scales.

The ability of HMHD to describe the plasma dynamics at sub-ion scales needs to be further investigated in other plasma regimes. A more detailed analysis highlighting the differences between HMHD and HPIC simulations may provide further constraints on additional phenomena beyond the Hall physics relevant for plasma turbulence at sub-ion scales. For instance, preliminary results from analyses based on the von Karman–Howarth equations show interesting similarities, but also differences, between the HMHD and HPIC simulations, thus complementing the present work. This will be the subject of future work.

The authors wish to thank William H. Matthaeus for useful discussions. 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. We acknowledge partial funding by Fondazione Cassa di Risparmio di Firenze under the project HIPERCRHEL. P.H. acknowledges grant 18-08861S of the Czech Science Foundation. 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 Cartesian 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).

# Appendix

## Reconnection Rates: Identification of Magnetic X/O-points

Characterizing the statistical properties of reconnection at a given time requires the measurements of several physical quantities at the magnetic X- and O-points inside all current sheets present in the simulation. In 2D models, these points are the critical points (i.e., minima, maxima, and saddles) of the out-of-plane vector potential $A_z$. To that purpose, we designed an algorithm that, for each output of the simulation, performs the following operations:

1. identification of critical X/O-points of $A_z$;

2. calculation of $A_z$, $\partial_t A_z$, and other ancillary quantities at the critical points;

3. identification of current sheet structures, if present;

4. selection and pairing of the X/O-points located inside each current sheet; and

5. calculation of the reconnection rate of each XO-pair.

The algorithm is implemented in IDL, and employs Fourier decomposition to calculate all derivatives. The first two steps are common to other algorithms used for the analysis of magnetic reconnection in turbulence (e.g., Servidio et al. 2010; Haggerty et al. 2017). In those studies, the reconnection rate is defined as the value of $\partial_t A_z / b_{rms}^2$ at the X-points, normalized with respect to the mean-square magnetic field $b_{rms}^2$. Instead, we use the following definition (see Equation (5))

$$ \gamma_{\text{rec}} = \left| \frac{1}{\Phi|_{X}^{O}} \frac{\partial \Phi|_{X}^{O}}{\partial t} \right| = \left| \frac{1}{A_{z}^{O} - A_{z}^{X}} \partial_{t}(A_{z}^{O} - A_{z}^{X}) \right|, \quad (7) $$

which self-consistently defines a quantity with the desired dimensions (i.e., the inverse of a time), without the need for any further normalization. This definition is also consistent with the classic definition of the growth rate adopted for the tearing instability, as shown in a previous study (Papini et al. 2018). The above steps are implemented as follows.

**Step 1:** We begin by taking the components of the gradient of $A_z$, which in our models correspond to the in-plane components of the magnetic field

$$ \partial_x A_z = -B_y, \quad \partial_y A_z = B_x. \quad (8) $$

We then compute the zero contour lines of $B_x$ and $B_y$, by using the IDL ISOCONTOUR procedure. The critical points are the intersection points of such lines, where $|\nabla A_z| = 0$. We found that this method is less prone to the detection of spurious critical points (Wan et al. 2010) than the use of other interpolation techniques. The remaining operations are analogous to Servidio et al. (2010). For each of the critical points, we compute the determinant and the eigenvalues of the Hessian matrix

$$ \mathcal{H}(x_i, y_i) = \begin{vmatrix} \partial_x^2 A_z(x_i, y_i) & \partial_x \partial_y A_z(x_i, y_i) \\ \partial_x \partial_y A_z(x_i, y_i) & \partial_y^2 A_z(x_i, y_i) \end{vmatrix} \quad (9) $$

at the coordinates $(x_i, y_i)$ of the critical point (which are usually located between the grid points), by means of bilinear interpolation of the derivatives of $A_z$ at $(x_i, y_i)$. Finally, positive and negative values of $\det(\mathcal{H})$ identify O-points and X-points respectively.

**Step 2:** $A_z$ is given as output by the HMHD code, while in CAMELIA is calculated via Fourier integration of the magnetic fields $B_x$ and $B_y$. The time derivative, $\partial_t A_z$, is calculated by taking the z-component of the induction equation for the vector potential

$$ \partial_t A_z = (v \times B)_z - \eta(\nabla \times B)_z - \eta_H \frac{((\nabla \times B) \times B)_z}{\rho}. \quad (10) $$

Finally, the values at the X/O-points are found via interpolation.

**Step 3:** Current sheets in the simulations are defined by the regions where the amplitude of the current density is larger than a given threshold $\epsilon_0$. The choice of $\epsilon_0$ is somewhat arbitrary. For our purposes, a good choice is to set $\epsilon_0$ equal to the maximum between 10% of the maximum amplitude of $J$ in the simulation box and 10 times its mean, at each different time.

**Step 4:** We select the X-points belonging to each current sheet separately and we pair them with the nearest O-point belonging to the same current sheet. We discard the pairs whose distance is less than four grid points.

**Step 5:** Using Equation (7), we obtain the reconnection rate of each XO-pair from the values of $A_z$ and $\partial_t A_z$ calculated in step 2.

## ORCID iDs

Emanuele Papini https://orcid.org/0000-0002-7969-7415

Luca Franci https://orcid.org/0000-0002-7419-0527

Simone Landi https://orcid.org/0000-0002-1322-8712

Andrea Verdini https://orcid.org/0000-0003-4380-4837

Lorenzo Matteini https://orcid.org/0000-0002-6276-7771

Petr Hellinger https://orcid.org/0000-0002-5608-0834
References

Alexandrova, O., Saur, J., Lacombe, C., et al. 2009, PhRvL, 103, 165003

Biskamp, D., Schwarz, E., & Drake, J. F. 1996, PhRvL, 76, 1264

Biskamp, D., Schwarz, E., Zeiler, A., Celani, A., & Drake, J. F. 1999, PhPl, 6, 751

Biskamp, D., & Welter, H. 1989, PhFlB, 1, 1964

Boldyrev, S., & Perez, J. C. 2012, ApJL, 758, L44

Camporeale, E., Sorriso-Valvo, L., Califano, F., & Retinò, A. 2018, PhRvL, 120, 125101

Carbone, V., Veltri, P., & Mangeney, A. 1990, PhFlA, 2, 1487

Cerri, S. S., & Califano, F. 2017, NJPh, 19, 025007

Chen, C. H. K., Bale, S. D., Salem, C., & Mozer, F. S. 2011, ApJL, 737, L41

Chen, C. H. K., Bale, S. D., Salem, C. S., & Maruca, B. A. 2013a, ApJ, 770, 125

Chen, C. H. K., & Boldyrev, S. 2017, ApJ, 842, 122

Chen, C. H. K., Boldyrev, S., Xia, Q., & Perez, J. C. 2013b, PhRvL, 110, 259901

Chen, C. H. K., Salem, C. S., Bonnell, J. W., Mozer, F. S., & Bale, S. D. 2012, PhRvL, 109, 035001

Cho, J., & Lazarian, A. 2004, ApJL, 615, L41

Cho, J., & Lazarian, A. 2009, ApJ, 701, 236

Del Sarto, D., & Pegoraro, F. 2018, MNRAS, 475, 181

Del Sarto, D., Pegoraro, F., & Califano, F. 2016, PhRvE, 93, 053203

Denskat, K. U., Beinroth, H. J., & Neubauer, F. M. 1983, JGZG, 54, 60

Finn, J. M., & Kaw, P. K. 1977, PhFl, 20, 72

Franci, L., Hellinger, P., Guarrasi, M., et al. 2018a, J. Phys.: Conf. Series, 1031, 012002

Franci, L., Hellinger, P., Matteini, L., Verdini, A., & Landi, S. 2016a, in AIP Conf. Ser. 1720, Solar Wind 14, ed. L. Wang et al. (Melville, NY: AIP), 040003

Franci, L., Landi, S., Matteini, L., Verdini, A., & Hellinger, P. 2015a, ApJ, 812, 21

Franci, L., Landi, S., Matteini, L., Verdini, A., & Hellinger, P. 2016b, ApJ, 833, 91

Franci, L., Landi, S., Verdini, A., Matteini, L., & Hellinger, P. 2018b, ApJ, 853, 26

Franci, L., Verdini, A., Matteini, L., Landi, S., & Hellinger, P. 2015b, ApJL, 804, L39

Franci, L., Cerri, S. S., Califano, F., et al. 2017, ApJL, 850, L16

Galtier, S. 2006, JPIPh, 72, 721

Galtier, S. 2008, PhRvE, 77, 015302

Galtier, S., & Buchlin, E. 2007, ApJ, 656, 560

Ghosh, S., Hossain, M., & Matthaeus, W. H. 1993, CoPhC, 74, 18

Ghosh, S., Siregar, E., Roberts, D. A., & Goldstein, M. L. 1996, JGR, 101, 2493

Goldstein, M. L., Roberts, D. A., & Fitch, C. A. 1994, JGR, 99, 11519

Grappin, R., Müller, W.-C., & Verdini, A. 2016, A&A, 589, A131

Haggerty, C. C., Parashar, T. N., Matthaeus, W. H., et al. 2017, PhPl, 24, 102308

Hellinger, P., Verdini, A., Landi, S., Franci, L., & Matteini, L. 2018, ApJL, 857, L19

Howes, G. G., Tenbarge, J. M., Dorland, W., et al. 2011, PhRvL, 107, 035004

Kolmogorov, A. 1941, DoSSR, 30, 301

Landi, S., Del Zanna, L., Papini, E., Pucci, F., & Velli, M. 2015, ApJ, 806, 131

Loureiro, N. F., & Boldyrev, S. 2017, ApJ, 850, 182

Loureiro, N. F., Schekochihin, A. A., & Cowley, S. C. 2007, PhPl, 14, 100703

Mallet, A., Schekochihin, A. A., & Chandran, B. D. G. 2017, JPIPh, 83, 905830609

Matteini, L., Alexandrova, O., Chen, C. H. K., & Lacombe, C. 2017, MNRAS, 466, 945

Matthaeus, W. H., & Lamkin, S. L. 1986, PhFl, 29, 2513

Matthaeus, W. H., Servidio, S., & Dmitruk, P. 2008, PhRvL, 101, 149501

Matthews, A. P. 1994, JCoPh, 112, 102

Meyrand, R., & Galtier, S. 2013, PhRvL, 111, 264501

Mininni, P. D., & Pouquet, A. 2009, PhRvE, 80, 025401

Müller, W.-C., & Grappin, R. 2005, PhRvL, 95, 114502

Orszag, S. A. 1971, JAtS, 28, 1074

Papini, E., Landi, S., & Zanna, L. D. 2018, J. Phys.: Conf. Series, 1031, 012020

Salem, C. S., Howes, G. G., Sundkvist, D., et al. 2012, ApJL, 745, L9

Servidio, S., Matthaeus, W. H., Shay, M. A., et al. 2010, PhPl, 17, 032315

Shaikh, D., & Shukla, P. K. 2009, PhRvL, 102, 045004

Sulem, P. L., Passot, T., Laveder, D., & Borgogno, D. 2016, ApJ, 818, 66

Wan, M., Oughton, S., Servidio, S., & Matthaeus, W. H. 2010, PhPl, 17, 082308

Yang, Y., Matthaeus, W. H., Parashar, T. N., et al. 2017a, PhRvE, 95, 061201

Yang, Y., Matthaeus, W. H., Parashar, T. N., et al. 2017b, PhPl, 24, 072306

Yin, L., Winske, D., Gary, S. P., & Birn, J. 2002, PhPl, 9, 2575