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

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).

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

$$\begin{eqnarray}&&{\partial }_{t}\rho +{\boldsymbol{\nabla }}\cdot (\rho {\boldsymbol{v}})=\,0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}\begin{array}{rcl}\rho \left({\partial }_{t}+{\boldsymbol{v}}\cdot {\boldsymbol{\nabla }}\right){\boldsymbol{v}} & = & -{\boldsymbol{\nabla }}P+({\boldsymbol{\nabla }}\times {\boldsymbol{B}})\times {\boldsymbol{B}}\\ & & +\ \nu \left[{\rm{\Delta }}{\boldsymbol{v}}+\displaystyle \frac{1}{3}{\boldsymbol{\nabla }}({\boldsymbol{\nabla }}\cdot {\boldsymbol{v}})\right],\end{array}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}\begin{array}{rcl}\left({\partial }_{t}+{\boldsymbol{v}}\cdot {\boldsymbol{\nabla }}\right)T & = & ({\rm{\Gamma }}-1)\left\{-({\boldsymbol{\nabla }}\cdot {\boldsymbol{v}})T+\eta \displaystyle \frac{| {\boldsymbol{\nabla }}\times {\boldsymbol{B}}{| }^{2}}{\rho }\right.\\ & & +\ \left.\displaystyle \frac{\nu }{\rho }\left[{\left({\boldsymbol{\nabla }}\times {\boldsymbol{v}}\right)}^{2}+\displaystyle \frac{4}{3}{\left({\boldsymbol{\nabla }}\cdot {\boldsymbol{v}}\right)}^{2}\right]\right\},\end{array}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}\begin{array}{rcl}{\partial }_{t}{\boldsymbol{B}} & = & {\rm{\nabla }}\times \left({\boldsymbol{v}}\times {\boldsymbol{B}}\right)+\eta {{\rm{\nabla }}}^{2}{\boldsymbol{B}}\\ & & -\ {\eta }_{H}{\rm{\nabla }}\times \displaystyle \frac{({\rm{\nabla }}\times {\boldsymbol{B}})\times {\boldsymbol{B}}}{\rho },\end{array}\end{eqnarray} \tag{ 4 }$$

where Γ = 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₀, a density ρ₀, an Alfvén speed ${c}_{{\rm{A}}}={B}_{0}/\sqrt{4\pi {\rho }_{0}}={{\rm{\Omega }}}_{i}{d}_{i}$, a pressure ${P}_{0}={\rho }_{0}{c}_{{\rm{A}}}^{2}$, and a temperature T₀ = (k_B/m_i) P₀/ρ₀. Here, the adimensional magnetic resistivity η is in units of d_ic_A and the Hall coefficient η_H = d_i/L is equal to 1. Ω_i = eB₀/m_ic 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 × L_y = 256 d_i × 256 d_i and a grid resolution of Δx = Δ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₀, 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 }}_{\mathrm{inj}}=2\pi {d}_{i}/{k}_{\perp }^{\mathrm{inj}}$, with ${k}_{\perp }^{\mathrm{inj}}{d}_{i}\simeq 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_rms = B_rms/B₀ ≃ 0.24. In the HPIC reference simulation, the plasma beta for ions and electrons and the magnetic resistivity are set to the values β_i = β_e = 1 and η = 5 × 10⁻⁴ respectively. In HMHD we set the plasma β = β_i + β_e = 2 accordingly. Instead, for the resistivity, we employ η = 10⁻³, 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 ≃ 40,000.

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 Ω_i⁻¹. 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, ${{\boldsymbol{J}}}_{z}$, at $t=45\,{{\rm{\Omega }}}_{i}^{-1}$, right after the 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.

Zoom In Zoom Out Reset image size

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

$$\begin{eqnarray}&&{\gamma }_{\mathrm{rec}}=\left|\displaystyle \frac{1}{{\rm{\Phi }}{| }_{{\rm{X}}}^{{\rm{O}}}}\displaystyle \frac{\partial {\rm{\Phi }}{| }_{{\rm{X}}}^{{\rm{O}}}}{\partial t}\right|,\end{eqnarray} \tag{ 5 }$$

where ${\rm{\Phi }}{| }_{{\rm{X}}}^{{\rm{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, ${\rm{\Phi }}{| }_{{\rm{X}}}^{{\rm{O}}}={A}_{z}^{{\rm{O}}}-{A}_{z}^{{\rm{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 $| {\boldsymbol{J}}|$ (c) and the rms of the current density ${\boldsymbol{J}}$ (d). For the HMHD run, the first reconnection events are detected at $t\simeq 30\,{{\rm{\Omega }}}_{i}^{-1}$ and concurrently with the first maximum of $| {\boldsymbol{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 }_{\mathrm{rec}}\rangle =0.08\,{{\rm{\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 }_{\mathrm{rec}}\rangle =1/\langle {\gamma }_{\mathrm{rec}}\rangle \simeq 12.2\,{{\rm{\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 }_{\mathrm{rec}}\lt 6\,{{\rm{\Omega }}}_{i}^{-1}$. After the maximum of turbulent activity is reached, i.e., after the maximum of rms($| {\boldsymbol{J}}|$) (Mininni & Pouquet 2009) at $t\simeq 165\,{{\rm{\Omega }}}_{i}^{-1}$, the number of reconnection events decreases.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Magnetic reconnection modifies the turbulent properties at kinetic scales soon after the first reconnection events are 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 }{{\boldsymbol{v}}}_{e}({k}_{\perp }))}^{-1}$, where ${{\boldsymbol{v}}}_{e}({k}_{\perp })$ is the electron velocity at the scale ℓ = 2π/k_⊥. 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 ≤ tΩ_i ≲ 30) in which the initial configuration relaxes and the first vortices and current sheets are formed, a transient phase (30 ≲ tΩ_i ≲ 50) triggered by the first reconnection events, and a third phase (50 ≲ tΩ_i) characterized by slowly evolving values of τ_nl(k_⊥) that lasts until the end of the simulation. In the first phase, at the beginning, the scale with the smallest value of τ_nl(k_⊥) is the injection scale ℓ_inj (marked by the horizontal dotted line), with ${\tau }_{{nl}}({k}_{\mathrm{inj}})\simeq 22\,{{\rm{\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\,{{\rm{\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 transfered to the smallest scales accessible (the 2/3 cutoff scale), where τ_nl (k_⊥) suddenly decreases. Then, during a transient phase that lasts between $t\simeq 30\,{{\rm{\Omega }}}_{i}^{-1}$ and $t\simeq 50\,{{\rm{\Omega }}}_{i}^{-1}$ (characterized by almost vertical isocontour lines of τ_nl(k_⊥); 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\,{{\rm{\Omega }}}_{i}^{-1}$. After this transient, the nonlinear time τ_nl(k_⊥) changes only slightly and then becomes roughly constant at $t\gtrsim 150\,{{\rm{\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.

Zoom In Zoom Out Reset image size

The evolution in the HPIC simulation is qualitatively the same, though the first reconnection events are triggered at a slightly later time ($t\simeq 40\,{{\rm{\Omega }}}_{i}^{-1}$) and the distribution of the reconnection rates is broader, with an average reconnection rate $\langle {\gamma }_{\mathrm{rec}}\rangle =0.22\,{{\rm{\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 ${{\boldsymbol{v}}}_{i,\parallel }$ 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 ${{\boldsymbol{v}}}_{i,\parallel }$ 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\,{{\rm{\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.

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 ${\boldsymbol{v}}$ in the HMHD description corresponds to the ion bulk velocity ${{\boldsymbol{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).

Zoom In Zoom Out Reset image size

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 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_⊥d_i ≃ 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 ${{\boldsymbol{B}}}_{\parallel }$ is compatible with a power law of spectral index −2.8, while the HMHD spectrum of ${{\boldsymbol{B}}}_{\parallel }$ is more compatible with a value of −3. At MHD scales the two spectra are different. There the HMHD spectrum of ${{\boldsymbol{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_⊥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 ≲ k_⊥d_i ≲ 1) we obtain a residual energy spectrum ${{\rm{E}}}_{R}={{\rm{P}}}_{{\boldsymbol{B}}}-{{\rm{P}}}_{{\boldsymbol{v}}}$ consistent with a slope of −2 (Figure 4(d)) and a total energy spectrum ${{\rm{E}}}_{T}={{\rm{P}}}_{{\boldsymbol{B}}}+{{\rm{P}}}_{{\boldsymbol{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:

$$\begin{eqnarray}&&{{\rm{E}}}_{R}/{{\rm{E}}}_{T}\approx A{\left({t}_{{\rm{A}}}/{t}_{{nl}}\right)}^{\alpha }\end{eqnarray} \tag{ 6 }$$

where the Alfvén time t_A = 1/(k_⊥b_rms) is built on the large-scale magnetic fluctuations, ${t}_{{nl}}=1/({k}_{\perp }\sqrt{{k}_{\perp }{{\rm{E}}}_{T}})$ is calculated using the total energy spectrum, the exponent on the r.h.s indicates a fast (α = 1) or slow (α = 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 (α = 1) and for the HMHD (α = 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 α = 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 ${{\rm{E}}}_{R}={{\rm{P}}}_{{{\boldsymbol{B}}}_{\perp }}-{{\rm{P}}}_{{{\boldsymbol{v}}}_{\perp }}$ and ${{\rm{E}}}_{T}\,={{\rm{P}}}_{{{\boldsymbol{B}}}_{\perp }}+{{\rm{P}}}_{{{\boldsymbol{v}}}_{\perp }}$). 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 α = 1,2 (thin solid and dashed lines, respectively). The thick and thin solid lines are parallel in the MHD range (0.3 ≲ k_⊥d_i ≲ 1) for both simulations, thus supporting a fast scenario (α = 1), while the slow scenario (α = 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 ≈ 3 and α = 1 are the same in both simulations.

Zoom In Zoom Out Reset image size

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 as a non-Gaussian behavior in the probability distribution functions (PDFs) of the increments of any field at a given scale ℓ. In Figure 7(a)–(c) we report the PDFs of the increments ΔB_y(x, y) = B_y(x + ℓ, y) − ΔB_y(x, y) of the y-component of the magnetic field fluctuations along the x-direction and at three different spatial separations ℓ/d_i = 5π, π, and π/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\,{{\rm{\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 (ℓ ≫ d_i), but as scales get smaller (e.g., already at ℓ = π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πd_i/ℓ ≈ 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πd_i/ℓ ≳ 0.3, well before the velocity spectra start to diverge (k_⊥d_i ≳ 2, see Figure 4 (b)).

Zoom In Zoom Out Reset image size

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 = ∂_tA_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₀ = 256d_i /2π, as done in the above works. Indeed, the PDF of $| {{\rm{E}}}_{{\rm{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.

Zoom In Zoom Out Reset image size

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 ≈−2. These indices are not consistent with 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 (α = 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 Cartesius based in the Netherlands at SURFsara through the DECI-13 (Distributed European Computing Initiative) call (project HybTurb3D), and CINECA for awarding access to HPC resources under the ISCRA initiative (grants HP10B2DRR4 and HP10C2EARF).

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, ∂_tA_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}_{\mathrm{rms}}^{2}$ at the X-points, normalized with respect to the mean-square magnetic field ${b}_{\mathrm{rms}}^{2}$. Instead, we use the following definition (see Equation (5))

$$\begin{eqnarray}&&{\gamma }_{\mathrm{rec}}=\left|\displaystyle \frac{1}{{\rm{\Phi }}{| }_{X}^{O}}\displaystyle \frac{\partial {\rm{\Phi }}{| }_{X}^{O}}{\partial t}\right|=\left|\displaystyle \frac{1}{{A}_{z}^{O}-{A}_{z}^{X}}{\partial }_{t}({A}_{z}^{O}-{A}_{z}^{X})\right|,\end{eqnarray} \tag{ 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

$$\begin{eqnarray}&&{\partial }_{x}{A}_{z}=-{B}_{y},\quad {\partial }_{y}{A}_{z}={B}_{x}.\end{eqnarray} \tag{ 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 $| {\rm{\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

$$\begin{eqnarray}{ \mathcal H }({x}_{i},{y}_{i})=\left|\begin{array}{ll}{\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{array}\right|\end{eqnarray} \tag{ 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, ∂_tA_z, is calculated by taking the z-component of the induction equation for the vector potential

$$\begin{eqnarray}&&{\partial }_{t}{A}_{z}={\left({\boldsymbol{v}}\times {\boldsymbol{B}}\right)}_{z}-\eta {\left({\rm{\nabla }}\times {\boldsymbol{B}}\right)}_{z}-{\eta }_{H}\displaystyle \frac{{\left(({\rm{\nabla }}\times {\boldsymbol{B}})\times {\boldsymbol{B}}\right)}_{z}}{\rho }.\end{eqnarray} \tag{ 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 ₀. The choice of ₀ is somewhat arbitrary. For our purposes, a good choice is to set ₀ equal to the maximum between 10% of the maximum amplitude of ${\boldsymbol{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 ∂_tA_z calculated in step 2.