von Kármán–Howarth Equation for Hall Magnetohydrodynamics: Hybrid Simulations

Rarefied magnetized plasmas are ubiquitous in the astrophysical context. There, Coulomb collisions are so rare that nonlinear couplings driven by large-scale medium motions induce a turbulent cascade that reaches the kinetic scales of the plasma constituent (ions and electrons) before collisional effects start to dissipate its flux. Turbulence in such collisionless/weakly collisional plasmas is not well understood and remains one of the challenging problems of astrophysics.

The magnetized plasma flow of the solar origin, the solar wind, constitutes a natural laboratory for studying turbulence in collisionless plasmas (Bruno & Carbone 2013). Indeed, in situ observations show that the solar wind flow is strongly turbulent, and that the magnetic and velocity fields exhibiting power-law spectral properties as well as non-Gaussian statistical properties (Matthaeus et al. 2015). On relatively large scales the magnetic field power spectrum of the observed time series is close to f^−5/3 reminding of the Kolmogorov prediction for hydrodynamic (HD) turbulence. This spectrum is, however, anisotropic with respect to the ambient magnetic field; along the magnetic field direction the spectrum is steeper (close to f⁻²; Chen 2016), in agreement with theoretical expectations and numerical simulations based on the magnetohydrodynamic (MHD) approximation.

Around the ion characteristic scales (proton gyroradius ρ_i and proton inertial length d_i, that are typically close to each other in the solar wind) the spectrum steepens. This steepening is related to a change of the physical behavior; in the sub-ion range the MHD approximation breaks and a more accurate approximation, including Hall and kinetic effects, is needed. The physical phenomena responsible for the steepening are not yet clearly determined, and many processes related to the dispersive (Hall) and dissipative (collisionless damping) phenomena ranges appear at similar scales (Marsch 2006; Alexandrova et al. 2008). The position of the transition from the large, MHD scales to the sub-ion scales (so-called "ion spectral break") varies with the radial distance and depends on the ion temperature (or on the ratio between the ion and magnetic pressures β_i; Bruno & Trenchi 2014; Chen et al. 2014).

The understanding of turbulence is strongly facilitated by the existence of exact dynamical equations (which involve second- and third-order structures functions), obtained using a statistical approach assuming a homogeneity (and isotropy) of the medium. The classical incompressible HD results (de Karman & Howarth 1938; Kolmogorov 1941) have been extended to the incompressible MHD (Chandrasekhar 1951; Politano & Pouquet 1998a, 1998b; Carbone et al. 2009), the incompressible Hall-MHD (Galtier 2008), and recently to the compressible Hall-MHD (Andrés et al. 2018). In the context of the solar wind, effects of a homogeneous velocity shear (Wan et al. 2009) or an expansion (Hellinger et al. 2013) were also investigated. Exact scaling laws, such as the well-known 4/3 and 4/5 laws, are obtained from these dynamical equations for the inertial range once a stationary state and diverging Reynolds numbers are assumed to be still leaving at a finite dissipative rate. The scaling predicted by the incompressible MHD exact laws were measured in the solar wind and used to estimate the energy cascade rate (and the corresponding particle heating rate) in the solar wind (Sorriso-Valvo et al. 2007; MacBride et al. 2008; Marino et al. 2008; Stawarz et al. 2009; Coburn et al. 2015). In addition to the observations, numerical simulations constitute an important means for investigating the properties of turbulence. The dynamical equations were tested in HD (Gotoh et al. 2002; Ishihara et al. 2009) and MHD simulations (Sorriso-Valvo et al. 2002; Mininni & Pouquet 2009; Verdini et al. 2015), in which a relatively good agreement with theoretical predictions is observed. However, due to the limited resolution, the range of scales where the scaling laws are verified is reduced (the scaling laws are exact only in the infinite Reynolds number limit).

For the transition between large MHD and sub-ion scales, numerical simulations based on the hybrid approximation, where electrons are treated as a fluid whereas ions are described fully kinetically, turned out to be very useful (Parashar et al. 2009; Servidio et al. 2012, 2015; Valentini et al. 2014; Vasquez et al. 2014). Recent 2.5D hybrid simulations (Franci et al. 2015a, 2015b, 2016) exhibit a clear double power-law behavior of the magnetic field fluctuations in agreement with in situ observations in the solar wind. Such simulations are natural candidates for testing the exact laws. As we are also interested in the sub-ion range, we want to look at the incompressible Hall-MHD (e.g., Pezzi et al. 2017).

Assuming incompressibility, ${\boldsymbol{\nabla }}\cdot {\boldsymbol{u}}=0$, and a constant plasma density for simplicity, ρ = const., taking the magnetic field and the electric current in Alfvén units ${\boldsymbol{b}}={\boldsymbol{B}}/\sqrt{{\mu }_{0}\rho }$, and ${\boldsymbol{j}}={\boldsymbol{J}}/({en})$, the Hall-MHD equation will have this form:

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{u}}}{\partial t}+({\boldsymbol{u}}\cdot {\boldsymbol{\nabla }}){\boldsymbol{u}}-({\boldsymbol{b}}\cdot {\boldsymbol{\nabla }}){\boldsymbol{b}}=-{\boldsymbol{\nabla }}P/\rho +\nu {\rm{\Delta }}{\boldsymbol{u}}\\ &&\displaystyle \frac{\partial {\boldsymbol{b}}}{\partial t}+({\boldsymbol{u}}\cdot {\boldsymbol{\nabla }}){\boldsymbol{b}}-({\boldsymbol{b}}\cdot {\boldsymbol{\nabla }}){\boldsymbol{u}}\\ &&\quad =\,({\boldsymbol{j}}\cdot {\boldsymbol{\nabla }}){\boldsymbol{b}}-({\boldsymbol{b}}\cdot {\boldsymbol{\nabla }}){\boldsymbol{j}}+\eta {\rm{\Delta }}{\boldsymbol{b}};\end{eqnarray} \tag{ 1 }$$

here, P denotes the scalar total (particle and magnetic field) pressure, and ν and η denote the kinematic plasma viscosity and the electric resistivity, respectively. Following Carbone et al. (2009) we take the increment of the different quantities at ${\boldsymbol{x}}$ and ${{\boldsymbol{x}}}^{{\prime} }={\boldsymbol{x}}+{\boldsymbol{l}}$ (${\boldsymbol{l}}$ being the spatial lag), $\delta {\boldsymbol{u}}={\boldsymbol{u}}({{\boldsymbol{x}}}^{{\prime} })-{\boldsymbol{u}}({\boldsymbol{x}})$, ... and, assuming an isotropic turbulence averaging over ${\boldsymbol{x}}$ (the averaging is denoted by $\langle \ \rangle$), we get the following equation for the magnetic, velocity, and total second-order structure functions ${S}_{b}=\langle {| \delta {\boldsymbol{b}}| }^{2}\rangle$, ${S}_{u}=\langle {| \delta {\boldsymbol{u}}| }^{2}\rangle$, and $S={S}_{b}+{S}_{u}$ as functions of ${\boldsymbol{l}}$

$$\begin{eqnarray}&&\displaystyle \frac{\partial S}{\partial t}+{\boldsymbol{\nabla }}\cdot ({\boldsymbol{Y}}+{\boldsymbol{H}})+A=-4\epsilon +2\nu {\rm{\Delta }}{S}_{u}+2\eta {\rm{\Delta }}{S}_{b}\end{eqnarray} \tag{ 2 }$$

where the third-order structure function ${\boldsymbol{Y}}$ is the MHD turbulent cascade flux (Carbone et al. 2009; Verdini et al. 2015)

$$\begin{eqnarray}&&{\boldsymbol{Y}}=\langle \delta {\boldsymbol{u}}{| \delta {\boldsymbol{u}}| }^{2}+\delta {\boldsymbol{u}}{| \delta {\boldsymbol{b}}| }^{2}-2\delta {\boldsymbol{b}}(\delta {\boldsymbol{u}}\cdot \delta {\boldsymbol{b}})\rangle \end{eqnarray} \tag{ 3 }$$

and the third-order structure function ${\boldsymbol{H}}$ is its Hall correction (e.g., Galtier 2008 for a different form of the Hall contribution)

$$\begin{eqnarray}&&{\boldsymbol{H}}=\langle 2\delta {\boldsymbol{b}}(\delta {\boldsymbol{b}}\cdot \delta {\boldsymbol{j}})-\delta {\boldsymbol{j}}{| \delta {\boldsymbol{b}}| }^{2}\rangle .\end{eqnarray} \tag{ 4 }$$

The last term at the left-hand side of Equation (2), $A\,=\langle \delta {\boldsymbol{j}}\cdot \delta [({\boldsymbol{b}}\cdot {\boldsymbol{\nabla }}){\boldsymbol{b}}]\rangle$, is a correction that we expect to be negligible in homogeneous plasma turbulence (this term is small in the present numerical simulations), and is the dissipation rate $\epsilon =\nu \langle {\boldsymbol{\nabla }}{\boldsymbol{u}}:{\boldsymbol{\nabla }}{\boldsymbol{u}}\rangle +\eta \langle {\boldsymbol{\nabla }}{\boldsymbol{b}}:{\boldsymbol{\nabla }}{\boldsymbol{b}}\rangle$ (here ":" denotes the double contraction of two second-order tensors). Equation (2) is a dynamical equation that relates second- and third-order structure function generalizing the von Kármán–Howarth equation in the framework of the incompressible Hall-MHD equations. An exact scaling law for a formally infinite-extending inertial range can be obtained assuming a stationary turbulent state with diverging Reynolds numbers but retaining a finite dissipation rate, ${\boldsymbol{\nabla }}\cdot ({\boldsymbol{Y}}+{\boldsymbol{H}})=-4\epsilon$, and, assuming isotropy, one gets the scaling law (e.g., Carbone et al. 2009)

$$\begin{eqnarray}&&{Y}_{r}+{H}_{r}=-\displaystyle \frac{4}{3}\epsilon l\end{eqnarray} \tag{ 5 }$$

where $l=| {\boldsymbol{l}}|$, Y_r and H_r are the radial components (in the spherical coordinates corresponding to the lag space ${\boldsymbol{l}}$) of ${\boldsymbol{Y}}$ and ${\boldsymbol{H}}$, respectively. In a general anisotropic case of magnetized MHD the behavior is more complicated (e.g., Verdini et al. 2015).

Note that the Hall term ${\boldsymbol{H}}$, Equation (4), represents a correction to the mixed terms in the MHD term ${\boldsymbol{Y}}$, Equation (3), replacing the ion velocity by the electron one, ${{\boldsymbol{u}}}_{e}$; i.e., the magnetic field couples to the electron velocity field

$$\begin{eqnarray}&&{\boldsymbol{Y}}+{\boldsymbol{H}}=\langle \delta {\boldsymbol{u}}{| \delta {\boldsymbol{u}}| }^{2}+\delta {{\boldsymbol{u}}}_{e}{| \delta {\boldsymbol{b}}| }^{2}-2\delta {\boldsymbol{b}}(\delta {{\boldsymbol{u}}}_{e}\cdot \delta {\boldsymbol{b}})\rangle \end{eqnarray} \tag{ 6 }$$

as one may expect.

Now we can directly test the prediction of the dynamic law, Equation (2), in the hybrid simulation results where ions are described by a particle-in-cell model whereas electrons are a massless, charge-neutralizing fluid (Matthews 1994). The simulation setup is an extension of previous simulations (Franci et al. 2015a, 2015b, 2016) and uses the 2D version of the code Camelia (http://terezka.asu.cas.cz/helinger/camelia.html). The evolution of initially isotropic protons with different values of ${\beta }_{i}$ (1/16, 1/2, and 4) is numerically integrated in a 2D domain (x, y) of size $256{d}_{i}\times 256{d}_{i}$ and resolution Δx = Δy. In order to reduce the noise, a Gaussian smoothing on 3 × 3 points is used on the proton density and velocity in the code. A uniform ambient magnetic field ${{\boldsymbol{B}}}_{0}$, directed along z and perpendicular to the simulation domain, is present, whereas neutralizing electrons are assumed isotropic and isothermal with β_e = β_i. The system is perturbed with an isotropic 2D spectrum of modes with random phases, linear Alfvén polarization ($\delta {\boldsymbol{B}}\,\perp \,{{\boldsymbol{B}}}_{0}$), and vanishing correlation between magnetic field and velocity fluctuations. These modes are in the range $0.02\leqslant {{kd}}_{i}\leqslant 0.2$ and have a flat 1D power spectrum with rms fluctuations δB. The time step Δt for particle integration (the magnetic field is advanced with a smaller time step ${\rm{\Delta }}{t}_{B}={\rm{\Delta }}t/20$), the number of particle per cell N_ppc, and the resistivity η, used to avoid energy accumulation at the smallest scales, have been set specifically for each simulation; their values are reported in Table 1. We let the system evolve beyond the time when the fluctuations of the parallel current maximize. This indicates a presence of a well-developed turbulent cascade (Mininni & Pouquet 2009; Servidio et al. 2015); henceforth, we analyze properties of plasma turbulence at around these times (${t}_{d}=490{{\rm{\Omega }}}_{i}^{-1}$, $350{{\rm{\Omega }}}_{i}^{-1}$, and $290{{\rm{\Omega }}}_{i}^{-1}$) for each simulation.

Table 1.  List of Simulations and Their Relevant Parameters (η is Given in the Units of ${\mu }_{0}{v}_{A}^{2}/{{\rm{\Omega }}}_{i}$)

RUN	βi	δB/B0	${k}_{\mathrm{inj}}{d}_{i}$	Δx /di	η	Nppc	Δt Ωi
1	1/16	0.20	0.2	1/16	2 10−4	1024	0.0025
2	1/2	0.25	0.2	1/16	3 10−4	4096	0.01
3	4	0.25	0.2	1/8	5 10−4	32768	0.02

Download table as:  ASCIITypeset image

Figure 1 shows the power spectral density of the magnetic field for the three simulations. The simulated spectra exhibit two power laws with a smooth transition at ion scales (the ion spectral break), the shape and position of which depend on the plasma beta: its scale is close to d_i for small betas, whereas in high beta plasmas it is around ρ_i (Franci et al. 2016). The magnetic power spectral slopes on large scales is ∼ −5/3, whereas the velocity power spectrum is less steep, closer to −3/2 than to −5/3. In the sub-ion range, as already observed in (Franci et al. 2016), the magnetic power spectrum steepens, with slopes about −3.5, −3, −2.9 in the three simulations. The proton velocity fluctuating field decouples from the magnetic fluctuations around the proton gyroscales. The sub-ion velocity fluctuations have a limited scale range before reaching the noise level, so that it is difficult to distinguish between an exponential and a power-law dependence. Assuming the latter, as suggested by observations (Šafránková et al. 2016) and theoretical and numerical results (Meyrand et al. 2018), the power spectrum of the velocity fluctuations below the decoupling (above the noise level) is compatible with steep slopes and spectral indices between −5 and −6.

Zoom In Zoom Out Reset image size

Figure 1 also shows the power spectral density of the plasma density ρ. The initial, transverse fluctuations lead to the formation of density fluctuation; the rms relative density fluctuation $\delta \rho /{\rho }_{0}$ are 0.16, 0.081, and 0.017 for the three runs going from low to high beta. As one may expect, the density fluctuations are stronger in a low beta plasma.

The spectral features observed in Figure 1 are partly reflected in the second-order structure functions of the magnetic field, S_b, the proton velocity, S_u, and the proton density, S_ρ (see Figure 2). However, the properties of the inertial range, the spectral break, and the sub-ion range are less clear compared to the power spectra.

Zoom In Zoom Out Reset image size

For the test of the dynamical equation, Equation (2), we define the cascade rate * as

$$\begin{eqnarray}&&{\epsilon }^{* }=-\displaystyle \frac{1}{4}\displaystyle \frac{\partial S}{\partial t}-\displaystyle \frac{1}{4}{\boldsymbol{\nabla }}\cdot ({\boldsymbol{Y}}+{\boldsymbol{H}})+\displaystyle \frac{1}{2}\eta {\rm{\Delta }}{S}_{b}.\end{eqnarray} \tag{ 7 }$$

Figure 3 shows the test of the dynamic law, the cascade rate, *, normalized to the resistive heating rate, , with the different contributing terms to *. Note that in these simulations the correcting term A is not important, $| A| /4\lesssim 0.1\epsilon$ (not shown). The structure functions (and ) are calculated at two times (t_d + δt and t_d separated by $\delta t=10{{\rm{\Omega }}}_{i}^{-1}$) around the maximum turbulence activity and averaged; ∂S/∂t is estimated by $(S({t}_{d}+\delta t)-S({t}_{d}))/\delta t$.

Zoom In Zoom Out Reset image size

In all simulations we observe that at MHD scales ${\epsilon }^{* }\sim -\partial S/\partial t/4-{\boldsymbol{\nabla }}\cdot {\boldsymbol{Y}}$. The term $\partial S/\partial t$ is due to the energy decay at injection scales given by the (decaying) initial condition, whereas the region where ${\boldsymbol{\nabla }}\cdot {\boldsymbol{Y}}$ dominates covers the range of scales where a Kolomogorov-like spectrum is observed, (Figure 2) and represents the inertial MHD range. There, ${\boldsymbol{\nabla }}\cdot {\boldsymbol{Y}}$ varies only weakly and the radial component, Y_r (in the cylindrical coordinate system corresponding to ${\boldsymbol{l}}$), is roughly proportional to l (not shown here), representing an equivalent of the hydrodynamic exact scaling law.

Crossing l ∼ d_i, the Hall term flux ${\boldsymbol{\nabla }}\cdot {\boldsymbol{H}}$ starts to grow and becomes the dominant one at sub-ion scales (≲d_i), although this term is important only in a narrow range of scales and no linear scaling of H_r on l is observed. The sum of the two inertial contributions, ${\boldsymbol{\nabla }}\cdot ({\boldsymbol{Y}}+{\boldsymbol{H}})$ (the thin gray dashed lines in Figure 3), shows that the reduction in the MHD flux is partially compensated by the Hall flux and, where ${\boldsymbol{\nabla }}\cdot ({\boldsymbol{Y}}+{\boldsymbol{H}})$ is about constant, ${Y}_{r}+{H}_{r}$ is roughly proportional to l. This is especially true in the low β simulation, while a fraction of the flux becomes more and more lost as β increases. The position of the Larmor radius ρ_i for each simulation (the vertical dotted–dashed lines) suggests some kind of correlation between the decreasing inertial flux and ρ_i.

It is only at the smallest scales that the magnetic diffusion term, ηΔS_b/2, contributes significantly to the total cascade rate. In the low β simulations the sum of all contributions (black line) is quite constant at all scales. For the β_i = 1/2 a small reduction of the total cascade rate is observed near l ∼ d_i ∼ ρ_i, and such a reduction becomes significant in the high β simulation (around l ∼ ρ_i). In all cases, however, at the injection-MHD scales * is more than twice the Joule dissipation rate .

Even if we use a large number of particle per cell, a part of the smallest scales is influenced by the particle's noise, especially in the ion velocity field (the magnetic field is much less affected (Franci et al. 2015a)). However, ion velocity does not contribute to the Hall flux, ${\boldsymbol{H}}$, and comparison of two simulations with a different number of N_ppc (not shown here) confirms that the sub-ion cascade rate is only weakly affected by the noise level. The resistivity used in the hybrid simulations is a free parameter that is used to avoid accumulation of the cascading energy on small scales, and in some respects replaces the full electron physics. Increasing the resistivity typically leads to a reduction of the Hall cascade rate, because dissipative scales shift toward larger scales and overlap with the Hall term but * is not significantly modified.

In this work we have tested by means of hybrid simulations a dynamical vectorial law, Equation (2), which is the generalization of the von Kármán–Howarth equation for incompressible hydrodynamic turbulence in the framework of incompressible Hall-MHD equations. Simulations show that in the low β regime Equation (2) is reasonably well satisfied and an indication of the MHD inertial regime (where ${\boldsymbol{\nabla }}\cdot {\boldsymbol{Y}}\simeq \mathrm{const}$) is obtained. At sub-ion scales the decreasing of the MHD flux is compensated by an increase of the Hall term, ${\boldsymbol{\nabla }}\cdot {\boldsymbol{H}}$, which extends the inertial regime below the scale of the spectral break. Increasing the plasma β the Hall flux only partially compensates the reduction of the MHD one, meaning that in the intermediate and high β regimes Equation (2) breaks.

The break of Equation (2) at high β does not seem to be related to compressibility effects. Indeed, ratios between density and magnetic field power spectra reveal that compressible effects are less important, as β_i increases both at MHD and kinetic scales. However, the incompressible description of the Hall-MHD equations from which the dynamic law is derived does not take into account the other relevant quantity at ions scales, the ion-Larmor radius. For instance, this scale seems to be relevant in determining the break at ions scales in high beta plasmas (Chen et al. 2014; Franci et al. 2016) and characterizes the polarization properties of turbulence mediated by kinetic-Alfvén-wave-like fluctuations (e.g., Schekochihin et al. 2009; Boldyrev et al. 2013; Chen et al. 2013; Franci et al. 2015a). Alternatively, the reduction of the cascade rate, *, could be related to nonthermal features playing at the ion gyroradius also not retained in the Hall-MHD approximation with the scalar pressure. Indeed, the work of Del Sarto et al. (2016) indicates that pressure anisotropies/nongyrotropies appear at around the scale of ρ_i. In the low beta case ρ_i, is significantly smaller then d_i and at scale near the dissipative ones. As β increases, ρ_i becomes comparable to d_i and a larger portion of the ion velocity distribution function can interact with turbulent fluctuations. Consequently, the reduction of * moves toward larger scales and it is more pronounced.

Another difference between the kinetic simulations and the predicted cascade rate is its value. The cascade rate, *, is typically a factor two larger then resistive losses, . One possibility is that a non-negligible fraction of the dissipation is carried by numerical effects. The dissipation rate may be underestimated due to the second-order numerical scheme that introduces an effective dissipation ∝∇² in the magnetic field. If it is so, would be greater while * remains unchanged (except at the smallest scales), thus reducing the discrepancy between the two. The other possibility is that the energy cascading from MHD scales is partially transferred to ions via physical processes beyond the Hall-MHD approximation, possibly connected with nonthermal features as suggested by Yang et al. (2017). The ion heating rate in our simulations is indeed comparable to ${\epsilon }^{* }\mbox{--}\epsilon ;$ however, for a detailed study of the relation between * and the ion (and electron) heating rates one needs to include the full electron kinetics.

These results are very robust, as we observed similar behavior in all of our previous 2D hybrid simulations (Franci et al. 2016; Cerri et al. 2017); these simulations further indicate that the decrease of * in the sub-ion region evolves continuously from low to high beta plasma. The present work supports the validity of the estimates of the cascade rate in the solar wind (Sorriso-Valvo et al. 2007; MacBride et al. 2008; Stawarz et al. 2009; Marino et al. 2011; Coburn et al. 2015). However, it is necessary to include compressible effects (Banerjee & Galtier 2013; Hadid et al. 2017; Andrés et al. 2018) as well as the effects of the anisotropy of the cascade in the full 3D geometry (Verdini et al. 2015) in future work.

It is interesting to note that perpendicular spectral properties of 3D hybrid simulations are in many respects similar to their 2D counterparts (Franci et al. 2018); we expect an analogous behavior in the case of the exact law. The present work also needs to be extended to full particle simulations in order to compare the cascade rate determined from the exact law with the actual heating rates of the different particle species (Wu et al. 2013; Matthaeus et al. 2016). In conclusion, the present work suggests that the transition from the MHD to sub-ion scales is governed by a combination of an onset of Hall cascade (with the characteristic scale d_i) and a formation of ion nonthermal features related to dissipation (with the characteristic scale ρ_i).

The authors thank T. Tullio for continuous enriching support and acknowledge useful discussions with W. H. Matthaeus and L. Sorriso-Valvo. P.H. acknowledges grant 18-08861S of the Czech Science Foundation. The authors acknowledge CINECA for HPC resources under the ISCRA initiative (grants HP10BUUOJM, HP10BEANCY, HP10B2DRR4, HP10CGW8SW, HP10C04BTP). The (reduced) simulation data are available at the Virtual Mission Laboratory Portal (http://vilma.asu.cas.cz).