2. Numerical models

The dynamical behaviour of a plasma strongly depends on its frequencies. Here, we provide a very brief description of the three models used in the present work, which can properly focus on different ranges of frequency. At the lowest frequencies, where ions and electrons are locked together, and the plasma behaves like an electrically conducting fluid, a MHD model is good enough to describe the system. In this study, we use a HMHD simulation, which retains dispersive effects that become relevant at scales comparable to the proton skin depth. However, when collisions can be neglected, pressure anisotropy can develop in the plasma. In this case, and at somewhat higher frequency than the previous regime, a two-fluid description is needed, since electrons and ions can move relatively to each other. Therefore, we use a LF model which is able to take into account both pressure anisotropy and low-frequency kinetic effects, such as Landau damping. At higher frequencies, a kinetic model is required. We consider the HVM approach which describes the evolution of the proton distribution function, while electrons are treated as a fluid. All three models include the electron inertial effects with a proton-to-electron mass ratio $m_p /m_e =100$.

The HMHD and LF codes use a two-dimensional Fourier pseudo-spectral method to compute the spatial derivatives, while the time advance is performed with a third-order Runge–Kutta scheme. In these fluid codes, aliasing errors in the evaluation of nonlinear terms are treated by a $2/3$ truncation in the spectral space. In the HVM code, time advance is based on the time-splitting scheme (Cheng & Knorr Reference Cheng and Knorr1976), combined with the current advance method (Matthews Reference Matthews1994). In the Vlasov equation, advection in space and velocity is solved through an explicit upwind scheme (Van Leer Reference Van Leer1977). The HVM code uses fast Fourier transform operations only when computing the fields, and truncation is not needed since possible dealiasing errors are damped by the intrinsic dissipation of the numerical scheme.

The vector fields will evolve in a double periodic domain $D = [L,L]$, in the $(x,y)$ plane. $L = 2{\rm \pi} \times 25 d_p$ is the length of each direction $x$ and $y$. This domain is discretized with $N_x = N_y = 1024$ grid points. Additionally, the HVM simulation is run in a three-dimensional velocity domain $- 5 v_{\textrm {th},p} \leq \boldsymbol {v} \leq 5 v_{\textrm {th},p}$, discretized with $51^3$ points, and where $v_{\textrm {th},p} = \sqrt {T_p / m_p}$ is the thermal velocity of the protons. At $t=0$, the equilibrium was set by a homogeneous plasma embedded in a uniform background magnetic field, $\boldsymbol {B}_0 = (0,0,B_0)$, oriented in the $z$ direction. We set the proton plasma beta as $\beta = 2 v_{th,p} / v_A = 1$, where $v_A = B_0 / \sqrt {4 {\rm \pi} n_0 m_p}$ is the Alfvén speed. The initial perturbation of this equilibrium is performed by velocity (and magnetic) field fluctuations with random phases in the interval $2 \leq m \leq 6$, for the wavenumbers $k = 2 {\rm \pi} m /L$. Although the amplitudes of the magnetic perturbations are $b / B_0 \approx 0.3$, neither density fluctuations nor parallel variances are imposed. In fluid models, a magnetic diffusivity term with a coefficient $\eta = 2 \times 10^{-2}$ has been added.

2.1. Hall magnetohydrodynamics

The set of HMHD equations, in a non-dimensional form, reads

(2.1)\begin{gather} \frac{\partial n}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} (n \boldsymbol{u}) = 0, \end{gather}

(2.2)\begin{gather} {\frac{\partial \boldsymbol{u}}{\partial t} + (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla}) \boldsymbol{u} ={-} \frac{\beta}{2 n} \boldsymbol{\nabla} P + \frac{1}{n} \left[ (\boldsymbol{\nabla} \times \boldsymbol{B}) \times \boldsymbol{B} \right] - \eta_4 \nabla^4 \boldsymbol{u} }, \end{gather}

(2.3)\begin{gather} {\frac{\partial \boldsymbol{B}}{\partial t} ={-}\boldsymbol{\nabla} \times \boldsymbol{E} - \eta_4 \nabla^4 \boldsymbol{B}}, \end{gather}

(2.4)\begin{gather} {(1 - d_e^2 \nabla^2)} \boldsymbol{E} + \boldsymbol{u} \times \boldsymbol{B} -\frac{1}{n} \left[ (\boldsymbol{\nabla} \times \boldsymbol{B}) \times \boldsymbol{B} \right] {- \eta (\boldsymbol{\nabla} \times \boldsymbol{B})} = 0, \end{gather}

(2.5)\begin{gather} P = n^\gamma, \end{gather}

where $n$ is the plasma density, $\boldsymbol {u}$ is the hydrodynamic velocity, $\boldsymbol {E}$ and $\boldsymbol {B}$ are the electric and magnetic fields, respectively, $P$ is the isotropic total pressure, $\eta = 2 \times 10^{-2}$ is the magnetic diffusivity and $\gamma = 5/3$ is the adiabatic index. This non-dimensional set of equations is obtained normalizing $n$ to $n_0 m_p$, $\boldsymbol {u}$ to the Alfvén speed $v_A$, $t$ to the inverse proton cyclotron frequency $\varOmega _{\textrm {cp}}^{-1}$ and the length units to the proton skin depth $d_p = v_A / \varOmega _{\textrm {cp}}$. In Ohm's law, the electron inertia is retained in the form used in Valentini et al. (Reference Valentini, Trávníček, Califano, Hellinger and Mangeney2007), with a non-dimensional electron skin depth $d_e^2 = m_e/m_i$ taken equal to $1/100$. For regularization purposes, and to reproduce the same-time behaviour for the maximum value of the integrated parallel current density $\langle\,j_z^2\rangle$, bi-Laplacian hyperviscosity and hyperdiffusivity, with coefficients equal to $\nu _4=\eta _4=5\times 10^{-4}$, have been added to the right-hand side of velocity and magnetic field equations. Note that hyperviscosity and hyperdiffusivity strongly mitigate the effect of the electron inertia, and therefore the $d_e^2$ term is negligible in the present run.

2.2. Landau fluid

The LF system of dynamical equations solves for the evolution of the magnetic field, density, velocity, parallel and perpendicular pressures and heat fluxes in a way consistent with the low-frequency limit of the linear kinetic theory, up to a scale of the order of the ion gyroradius. The electric field is computed using a generalized Ohm's law, which includes the Hall term, together with electron inertia with the same simplified form and the same artificial mass ratio as in the HMHD case, and the same magnetic diffusivity. The electron pressure gradient, which is formally present in Ohm's law and in the momentum equation, is here neglected to closely match the conditions of the Vlasov simulations, in which electrons are cold. The ion velocity equation reads

(2.6)\begin{equation} \frac{\partial \boldsymbol{u}}{\partial t} + (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla}) \boldsymbol{u} ={-} \frac{\beta}{2 n} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{P} + \frac{1}{n} \left[ (\boldsymbol{\nabla} \times \boldsymbol{B}) \times \boldsymbol{B} \right], \end{equation}

while the density and the magnetic fields obey (2.1) and (2.3) of HMHD. The ion pressure tensor is conveniently written as $\boldsymbol {P} = p_{\perp }\boldsymbol {n} + p_{\|} \boldsymbol {\tau } + \boldsymbol {\varPi }$, with $\boldsymbol {\tau } = \hat {b}\otimes \hat {b}$, $\boldsymbol {n} = \boldsymbol {I} - \boldsymbol {\tau }$ and $\hat {b}=\boldsymbol {B}/|\boldsymbol {B}|$ is the unit vector along the local magnetic field. The parallel and perpendicular ion pressures obey

(2.7)\begin{gather} \partial_t p_{\|} + \boldsymbol{\nabla} \boldsymbol{\cdot} (\boldsymbol{u} \, p_{\|}) + 2 p_{\|}\, \hat{b}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}\, \boldsymbol{\cdot}\hat{b} + \hat{b}\boldsymbol{\cdot} (\boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{Q}) \boldsymbol{\cdot} \hat{b} \nonumber\\ +\left(\left(\boldsymbol{\varPi} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}\right)^S : \boldsymbol{\tau} - \boldsymbol{\varPi} :\frac{\textrm{d}\boldsymbol{\tau}}{\textrm{d}t}\right) =0, \\ \partial_t p_{{\perp}} + \boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{u} \, p_{{\perp}})+ p_{{\perp}} \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u} - p_{{\perp}}\,\hat{b}\boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u} \boldsymbol{\cdot} \hat{b} \nonumber\\ + \frac{1}{2}\left(\textrm{tr}\,\boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{Q} - \hat{b}\boldsymbol{\cdot} \left(\boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{Q}\right) \boldsymbol{\cdot} \hat{b} \right)\nonumber\end{gather}

(2.8)\begin{gather}+\frac{1}{2}\left(\textrm{tr}\, \left(\boldsymbol{\varPi}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}\right)^S -\left(\boldsymbol{\varPi} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}\right)^S : \boldsymbol{\tau} + \boldsymbol{\varPi} :\frac{\textrm{d}\boldsymbol{\tau} }{\textrm{d}t} \right) =0 , \end{gather}

thus involving the ion heat flux tensor $\boldsymbol {Q}$ along with the non-gyrotropic ion pressure tensor $\boldsymbol {\varPi }$. Here, and in (2.6), the latter will be neglected in the energy transfer budget since we are discarding finite Larmor radius effects, even though they may be relevant for energy transfers at the scale of the present simulations, and the trace term explicitly contributes to the energy conservation (see e.g. Yang et al. (Reference Yang, Matthaeus, Parashar, Wu, Wan, Shi, Chen, Roytershteyn and Daughton2017) for full kinetic simulations and Cerri & Camporeale (Reference Cerri and Camporeale2020) for explicit formulation in hybrid kinetic models). The non-gyrotropic pressure closure used in the present LF model (Passot et al. Reference Passot, Henri, Laveder and Sulem2014) is indeed valid only in a weakly nonlinear regime, while as we numerically checked it may lead to non-physical behaviour in the present context of a too strong nonlinearity, especially when electrons are cold. So here we focused on kinetic damping only. This point should also be kept in mind when interpreting the possible discrepancies between ion-scale features obtained from HVM simulation and from the fluid simulations (see e.g. § 5).

On the same footing, the ion heat flux tensor $\boldsymbol {Q}$ is taken in its gyrotropic form, which in (2.7) and (2.8) implies $\hat {b}\boldsymbol {\cdot } (\boldsymbol {\nabla } \boldsymbol {\cdot }\boldsymbol {Q}) \boldsymbol {\cdot } \hat {b} = \boldsymbol {\nabla }\boldsymbol {\cdot } (q_{\|} \hat {b} ) - 2q_{\perp r} \boldsymbol {\nabla } \boldsymbol {\cdot } \hat {b}$ and $( \textrm {tr} ( \boldsymbol {\nabla } \boldsymbol {\cdot }\boldsymbol {Q} ) - \hat {b}\boldsymbol {\cdot } (\boldsymbol {\nabla } \boldsymbol {\cdot }\boldsymbol {Q}) \boldsymbol {\cdot } \hat {b} )/2 = \boldsymbol {\nabla } \boldsymbol {\cdot } (q_{\perp } \hat {b}) + q_{\perp } \boldsymbol {\nabla } \boldsymbol {\cdot } \hat {b}$, where the gyrotropic heat flux components obey (still neglecting higher-order non-gyrotropic contributions)

(2.9)\begin{gather} \partial_t q_{\|} + \boldsymbol{\nabla}\boldsymbol{\cdot} (q_{\|} \boldsymbol{u}) + 3 q_{\|} \hat{b}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} \boldsymbol{\cdot} \hat{b} + 3 p_{\|} \left(\hat{b}\boldsymbol{\cdot} \boldsymbol{\nabla}\right)\left(\frac{p_{\|}}{\rho} \right) \nonumber\\ + \boldsymbol{\nabla}\boldsymbol{\cdot} ({\tilde{r}_{\|\|}} \hat{b})- 3 {\tilde{r}_{\|\perp}} \boldsymbol{\nabla}\boldsymbol{\cdot} \hat{b}=0,\\ \partial_t q_{{\perp}}+ \boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{u} q_{{\perp}}) + q_{{\perp}}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u} + p_{\|}\left(\hat{b}\boldsymbol{\cdot}\boldsymbol{\nabla}\right)\left(\frac{p_{{\perp}}}{\rho}\right) + \boldsymbol{\nabla} \boldsymbol{\cdot} \left({\tilde{r}_{\|\perp}} \hat{b}\right) \nonumber\end{gather}

(2.10)\begin{gather}+ \left(\left(p_{\|} -p_{{\perp}}\right)\frac{p_{{\perp}}}{\rho} -{\tilde{r}_{{\perp}{\perp}}}+ {\tilde{r}_{\|\perp}}\right) \left(\boldsymbol{\nabla}\boldsymbol{\cdot} \hat{b}\right) =0. \end{gather}

The gyrotropic fourth-rank cumulants ${\tilde {r}_{\|\|r} }$, ${\tilde {r}_{\|\perp r}}$ and ${\tilde {r}_{\perp \perp r}}$ needed to close the hierarchy are given by (5.1), (5.2) and (5.4) of Sulem & Passot (Reference Sulem and Passot2015), where ion Landau damping is retained only in a quasilinear form, in the approximation discussed in Passot et al. (Reference Passot, Henri, Laveder and Sulem2014).

Together with the same hyperviscosity and hyperdiffusivity included in HMHD equations, additional bi-Laplacian dissipative terms have been supplemented in the equations for the density and the pressures (with coefficient equal to $2.5 \times 10^{-3}$) and in the equations for the ion heat flux (with coefficient $10^{-4}$) in order to deal with the high level of compressibility in the simulation (the Mach number reaching values up to $0.4$). In fact, we note that compressibility of these fluid models will be important when comparing our results.

2.3. Hybrid Vlasov–Maxwell

The HVM code (Valentini et al. Reference Valentini, Trávníček, Califano, Hellinger and Mangeney2007) solves the Vlasov equation for the proton distribution function $f = f (\boldsymbol {r}, \boldsymbol {v}, t)$, while the electron response is taken into account through a generalized Ohm's law for the electric field. The Vlasov equation, in normalized units,

(2.11)\begin{equation} \frac{\partial f}{\partial t} + \boldsymbol{v} \boldsymbol{\cdot} \frac{\partial f}{\partial \boldsymbol{r}} + (\boldsymbol{E} + \boldsymbol{u} \times \boldsymbol{B}) \boldsymbol{\cdot} \frac{\partial f}{\partial \boldsymbol{v}} = 0, \end{equation}

is solved in a 2D-3V phase-space domain (two dimensions in physical space and three dimensions in velocity space), coupling it with (2.3) and (2.4). Velocities $\boldsymbol {v}$ and $\boldsymbol {u}$ (proton bulk velocity) are scaled to $v_A$. Quasi-neutrality ($n = n_e = n_p$), and cold electrons, is considered. The time step is chosen in order to satisfy the Courant–Friedrichs–Lewy condition for the numerical stability. The proton distribution function is initialized with a homogeneous-density Maxwellian function. In this procedure, displacement current is neglected in the Ampère law, making the assumption of low frequencies.

3. Local energy transfer analysis

The analysis developed in this work is performed in a period of maximal turbulence activity. As described by Perrone et al. (Reference Perrone, Passot, Laveder, Valentini, Sulem, Zouganelis, Veltri and Servidio2018), this state is reached at $t^* = 60 \varOmega _{\textrm {cp}}^{-1}$, for the same set of simulations considered here. They established $t^*$ after following $\langle\,j_z^2 \rangle$ in time, for HMHD, LF and HVM descriptions, where the same $\beta$ was held in order to ensure a similar level of density fluctuations. In figure 1, we plot the out-of-plane component of the current density, $j_z(x,y)$, at $t = t^*$. The structure locations (and intensities) present on this component are consistent for the three simulations: HMHD (figure 1a), LF (figure 1b) and HVM (figure 1c).

Figure 1.

Contour plot of the out-of-plane current density $j_z(x,y)$ for HMHD (a), LF (b) and HVM (c) simulations at $t = t^{\ast }$ (see text).

Figure 2 shows the spectra of $\boldsymbol {B}$ and $\boldsymbol {u}$ for the HMHD (figure 2a,d,g,j), LF (figure 2b,e,h,k) and HVM (figure 2c,f,i,l) simulations. Magnetic field spectra (bidimensional $|B(k_x,k_y)|^2$, and reduced integrated $|B_{ki}|^2$) are presented in figure 2(a–f), while the velocity field spectra (bidimensional $|u(k_x,k_y)|^2$, and reduced integrated $|u_{ki}|^2$) can be seen in figure 2(g–l). The integrated spectra are plotted with respect to $k_i = k_x$ with a red-dashed line, and with respect to $k_i = k_y$ with a black-solid line. Bidimensional spectra maps of HMHD and LF simulations show the $2/3$ truncation, which was not performed in the HVM simulation.

Figure 2.

Magnetic field (a–f) and velocity field (g–l) spectra. The integrated spectra are presented with respect to $k_i$, $i = x$ (red-dashed line), $y$ (black-solid line). Standard phenomenological power laws ($k_{i}^{-5/3}$ in the MHD range (see e.g. Frisch Reference Frisch1995) and $k_{i}^{-8/3}$ in the HMHD range (Meyrand & Galtier Reference Meyrand and Galtier2013)) are reported as reference.

Considering that, in the small-scale limit, HMHD physics was observed to match hybrid kinetic model results (Vásconez et al. Reference Vásconez, Pucci, Valentini, Servidio, Matthaeus and Malara2015; Pucci et al. Reference Pucci, Vásconez, Pezzi, Servidio, Valentini, Matthaeus and Malara2016), Hellinger et al. (Reference Hellinger, Verdini, Landi, Franci and Matteini2018) presented a generalization of the von Kármán–Howarth equation for incompressible hydrodynamic turbulence, in the framework of incompressible HMHD equations. This generalization was revisited by Ferrand et al. (Reference Ferrand, Galtier, Sahraoui, Meyrand, Andrés and Banerjee2019), who also considered homogeneous turbulence, and worked with two-point correlation tensors depending only on the relative displacement $\boldsymbol {\ell }$, and not on the absolute positions, as $\delta \boldsymbol {g} = \boldsymbol {g} (\boldsymbol {r} + \boldsymbol {\ell }) - \boldsymbol {g} (\boldsymbol {r})$. In this context, the mean rate of total energy injection $\epsilon$ is expressed as the combination of the third-order structure functions corresponding to the MHD turbulent cascade flux (Carbone, Sorriso-Valvo & Marino Reference Carbone, Sorriso-Valvo and Marino2009; Verdini et al. Reference Verdini, Grappin, Hellinger, Landi and Müller2015) and to Hall corrections (Galtier Reference Galtier2008):

(3.1)\begin{gather} {-}4 \epsilon = \boldsymbol{\nabla} \boldsymbol{\cdot} \langle (\delta \boldsymbol{u} \boldsymbol{\cdot} \delta \boldsymbol{u}) \delta \boldsymbol{u} + (\delta \boldsymbol{B} \boldsymbol{\cdot} \delta \boldsymbol{B}) \delta \boldsymbol{u} - 2 (\delta \boldsymbol{u} \boldsymbol{\cdot} \delta \boldsymbol{B}) \delta \boldsymbol{B} - \frac{d_p}{2} (\delta \boldsymbol{B} \boldsymbol{\cdot} \delta \boldsymbol{B}) \delta \boldsymbol{j}\nonumber\\ + d_p (\delta \boldsymbol{B} \boldsymbol{\cdot} \delta \boldsymbol{j}) \delta \boldsymbol{B} \rangle. \end{gather}

In fact, we can identify the single contribution of each component that corresponds to the Yaglom effect contributions,

(3.2)\begin{gather} \boldsymbol{Y}_1 (x,y) = (\delta \boldsymbol{u} \boldsymbol{\cdot} \delta \boldsymbol{u}) \delta \boldsymbol{u}, \end{gather}

(3.3)\begin{gather}\boldsymbol{Y}_2 (x,y) = (\delta \boldsymbol{B} \boldsymbol{\cdot} \delta \boldsymbol{B}) \delta \boldsymbol{u}, \end{gather}

(3.4)\begin{gather}\boldsymbol{Y}_3 (x,y) ={-}2 (\delta \boldsymbol{u} \boldsymbol{\cdot} \delta \boldsymbol{B}) \delta \boldsymbol{B}, \end{gather}

and to the Hall effect contributions,

(3.5)\begin{gather} \boldsymbol{H}_1 (x,y) = d_p (\delta \boldsymbol{B} \boldsymbol{\cdot} \delta \boldsymbol{j}) \delta \boldsymbol{B}, \end{gather}

(3.6)\begin{gather}\boldsymbol{H}_2 (x,y) ={-} \frac{d_p}{2} (\delta \boldsymbol{B} \boldsymbol{\cdot} \delta \boldsymbol{B}) \delta \boldsymbol{j}. \end{gather}

In this way, (3.1) would be written as $-4 \epsilon = \boldsymbol {\nabla } \boldsymbol {\cdot } \langle \boldsymbol {Y}_1 + \boldsymbol {Y}_2 + \boldsymbol {Y}_3 + \boldsymbol {H}_1 + \boldsymbol {H}_2 \rangle$. Compressive fluctuations (i.e. $\delta B$ and/or $\delta n$) naturally develop as the cascade proceeds towards smaller scales, becoming important at ion scales and below, even when these types of fluctuations are not directly injected at large (MHD) scales (which is the case for the simulations presented here). This effect was explicitly shown for two-dimensional simulations, for instance, in Cerri et al. (Reference Cerri, Franci, Califano, Landi and Hellinger2017). Direct numerical simulations of compressible MHD turbulence (Andrés et al. Reference Andrés, Sahraoui, Galtier, Hadid, Dmitruk and Mininni2018) and measurements in a collisionless space plasma using Magnetospheric Multiscale data in comparison with simulations of compressible HMHD turbulence (Andrés et al. Reference Andrés, Sahraoui, Galtier, Hadid, Ferrand and Huang2019) indicate that neglecting the contribution of compressibility to the exact laws results in an underestimation of the total energy flux and dissipation rate, especially near the ion scale. However, in Andrés et al. (Reference Andrés, Sahraoui, Galtier, Hadid, Dmitruk and Mininni2018) it was also found that except in cases in which the largest scales already display strong density fluctuations, the compressible Yaglom-like flux gives the dominant contribution to the energy flux at the MHD scales, and remains close to the incompressible Yaglom (or Politano–Pouquet) flux, while all other contributions of the compressibility to the energy flux remain small.

On the other hand, a heuristic proxy has been recently introduced. It focuses on the local turbulent energy transfer rate towards the smallest resolved scale. This proxy, named as LET, was constructed in order to extend the Yaglom law to MHD turbulence, and matches the latter law when small density fluctuations (at the scale $\ell$) can be neglected (Carbone et al. Reference Carbone, Marino, Sorriso-Valvo, Noullez and Bruno2009). It is estimated, for two-dimensional incompressible turbulence, through the combined third-order fluctuations of velocity, magnetic field and current density (Sorriso-Valvo et al. Reference Sorriso-Valvo, Perrone, Pezzi, Valentini, Servidio, Zouganelis and Veltri2018), and neglecting unity-order multiplicative factors, as

(3.7)\begin{gather} - 2 \ell \epsilon_\ell \equiv (\delta \boldsymbol{u} \boldsymbol{\cdot} \delta \boldsymbol{u}) \delta u_\ell + (\delta \boldsymbol{B} \boldsymbol{\cdot} \delta \boldsymbol{B}) \delta u_\ell - 2 (\delta \boldsymbol{u} \boldsymbol{\cdot} \delta \boldsymbol{B}) \delta B_\ell - \frac{d_p}{2} (\delta \boldsymbol{B} \boldsymbol{\cdot} \delta \boldsymbol{B}) \delta j_\ell\nonumber\\ + d_p (\delta \boldsymbol{B} \boldsymbol{\cdot} \delta \boldsymbol{j}) \delta B_\ell. \end{gather}

Consistently with our notation, the latter equation could be rewritten as

(3.8)\begin{equation} -2 \ell \epsilon_\ell = \epsilon_Y + \epsilon_H, \end{equation}

where $\epsilon _Y = \epsilon _{Y_1} + \epsilon _{Y_2} + \epsilon _{Y_3}$ and $\epsilon _H = \epsilon _{H_1} + \epsilon _{H_2}$. For this, $\epsilon _{Y_1} = (\delta \boldsymbol {u} \boldsymbol {\cdot } \delta \boldsymbol {u}) \delta u_\ell$ measures the kinetic energy available to be transported by the longitudinal component of $\delta \boldsymbol {u}$, $\epsilon _{Y_2} = (\delta \boldsymbol {B} \boldsymbol {\cdot } \delta \boldsymbol {B}) \delta u_\ell$ quantifies the magnetic energy that will be transported by $\delta u_\ell$, $\epsilon _{Y_3} = - 2 (\delta \boldsymbol {u} \boldsymbol {\cdot } \delta \boldsymbol {B}) \delta B_\ell$ is related to the velocity–magnetic field correlations coupled to the longitudinal magnetic field fluctuations, $\epsilon _{H_1} = -d_p (\delta \boldsymbol {B} \boldsymbol {\cdot } \delta \boldsymbol {B}) \delta j_\ell /2$ indicates the quantity of magnetic energy advected by longitudinal current density field fluctuations and $\epsilon _{H_2} = d_p (\delta \boldsymbol {B} \boldsymbol {\cdot } \delta \boldsymbol {j}) \delta B_\ell$ is related to the magnetic–current density field correlations coupled to the longitudinal magnetic field fluctuations. In (3.7), the two latter terms vanish when $d_p \to 0$, which recovers the classic Yaglom law. It should be noted that the unaveraged third-order quantity described in (3.7) cannot be rigorously derived, and may therefore display convergence issues. In other words, it does not describe a turbulent cascade, but just the local contribution of the turbulent fluctuations to the energy transfer. On the other hand, its use as a proxy is partly supported by the existence of local versions of the third-order scaling law, whose rigorous derivation is based on the use of scale- and space-localized filters (Duchon & Robert Reference Duchon and Robert2000; Eyink Reference Eyink2003; Camporeale et al. Reference Camporeale, Sorriso-Valvo, Califano and Retinò2018; Coburn & Sorriso-Valvo Reference Coburn and Sorriso-Valvo2019).

Since we want to explore the energy flux distribution resulting from the turbulent nonlinear transfer, we first focus on the scale $\sim d_p$, i.e. at the bottom of the inertial range, where the intermittent structures carrying most of the energy have been formed. As we noted in § 2, $N_x = N_y = 1024$, then ${\rm \Delta} x = {\rm \Delta} y \approx 0.15$, resulting in $\ell \approx d_{p}$. At this scale, we will identify the Yaglom $\boldsymbol {Y}_m(x,y) = (\boldsymbol {Y}_{mx} + \boldsymbol {Y}_{my})/2$ and Hall $\boldsymbol {H}_n(x,y) = (\boldsymbol {H}_{n x} + \boldsymbol {H}_{n y})/2$ contributions, according to their spatial distribution, where $m = 1, 2, 3$ and $n = 1, 2$. These two-dimensional maps are presented in figure 3 (HMHD), figure 4 (LF) and figure 5 (HVM). They clearly highlight the intermittent nature of the cross-scale energy transfer. This is evidenced by the presence of strong contributions to the scaling law in small-scale localized, sparse structures. Most of the energy transfer, and consequently of the energy dissipation, is taking place at those locations, which was also reported in Yang et al. (Reference Yang, Matthaeus, Parashar, Wu, Wan, Shi, Chen, Roytershteyn and Daughton2017) and in Camporeale et al. (Reference Camporeale, Sorriso-Valvo, Califano and Retinò2018). This feature may be consistent with the idea that reconnection of current sheets, and the formation of coherent structures, could provide a viable mechanism for (likely cross-scale) energy transfer across ion-kinetic scales (Cerri & Califano Reference Cerri and Califano2017; Franci et al. Reference Franci, Cerri, Califano, Landi, Papini, Verdini, Matteini, Jenko and Hellinger2017). It is interesting to note that different terms of the HMHD law may present similar structures at the same locations. However, in several instances the structures of different contributions are not collocated, suggesting that the nature of the fluctuations producing the energy transfer may change with position, possibly enabling different dissipation mechanisms (Sorriso-Valvo et al. Reference Sorriso-Valvo, Catapano, Retinò, Le Contel, Perrone, Roberts, Coburn, Panebianco, Valentini and Perri2019). It is also apparent that most of the more intense structures are positive, suggesting that the energy is being transferred towards the small scales in a ‘direct’ cascade. However, several negative structures are also present, where the contribution from the various terms is acting so as to remove energy from the smaller scales. We do not imply here that the sign of the local contributions is related to the cascade direction, and the latest statements should be taken as a qualitative indication. Note that the overall energy flux will result from the average over the whole domain, and is mostly positive, as expected and as we discuss in the next section. Comparing the different simulations reveals that the various terms have similar behaviour. In particular, the $\boldsymbol {Y}_3$ term seems to be dominating in the simulations. Interestingly, in the HMHD simulation the Hall terms appear weaker than in the two other simulations. This observation will be discussed later.

Figure 3.

The HMHD simulation. Bidimensional maps of Yaglom components $\boldsymbol {Y}_1$ (equation (3.2)), $\boldsymbol {Y}_2$ (equation (3.3)) and $\boldsymbol {Y}_3$ (equation (3.4)); and Hall components $\boldsymbol {H}_1$ (equation (3.5)) and $\boldsymbol {H}_2$ (equation (3.6)), of the mean rate of total energy injection, estimated at the scale $\ell \simeq d_p$.

Figure 4.

The LF simulation. Bidimensional maps of Yaglom components $\boldsymbol {Y}_1$ (equation (3.2)), $\boldsymbol {Y}_2$ (equation (3.3)) and $\boldsymbol {Y}_3$ (equation (3.4)); and Hall components $\boldsymbol {H}_1$ (equation (3.5)) and $\boldsymbol {H}_2$ (equation (3.6)), of the mean rate of total energy injection, estimated at the scale $\ell \simeq d_p$.

Figure 5.

The HVM simulation. Two-dimensional maps of Yaglom components $\boldsymbol {Y}_1$ (equation (3.2)), $\boldsymbol {Y}_2$ (equation (3.3)) and $\boldsymbol {Y}_3$ (equation (3.4)); and Hall components $\boldsymbol {H}_1$ (equation (3.5)) and $\boldsymbol {H}_2$ (equation (3.6)), of the mean rate of total energy injection, estimated at the scale $\ell \simeq d_p$.

Finally, the maps of $\boldsymbol {Y} + \boldsymbol {H}$ are separated by components in figure 6 for the HVM simulation (similar behaviour was found in HMHD and LF simulations – not shown here). While the general behaviour is similar, some differences between the two components of the energy transfer exist, possibly an effect of the finite size of the ensemble.

Figure 6.

The $x$ (a) and $y$ (b) components of the term $\boldsymbol {Y} + \boldsymbol {H}$ (c), computed from the HVM simulation, estimated at the scale $\ell \simeq d_p$.

4. Global energy transfer analysis

We continue our analysis computing the local turbulent energy transfer rate $- 2 \epsilon _\ell$ across the scales $\ell$. Figure 7 presents the averaged (3.7) over the whole domain, corresponding to the HMHD (figure 7a), LF (figure 7b) and HVM (figure 7c) simulations. We separate the Yaglom (black-solid line) and Hall (red-solid line) contributions in order to highlight their relative amplitude across the scales $\ell$. The sum of both contributions, giving the Hall–Yaglom law, is plotted with asterisks. As already observed in the previous section, in the HMHD and LF simulations, the Hall contribution $\epsilon _{H(\ell )}$ is always less than the Yaglom one. However, both contributions become closer for scales $\ell \lesssim d_p$. In the case of the HVM simulation, that occurs in the range $\ell \lesssim d_p$, $\epsilon _{H(\ell )} > \epsilon _{Y(\ell )}$. Then, at $\ell \approx d_p$, $\epsilon _{Y(\ell )}$ becomes more relevant and remains so for larger scales. This behaviour is consistent with the settings of our simulations, where the mean field is out-of-plane, and there is strong whistler/magnetosonic activity. Linear theory studied by Camporeale & Burgess (Reference Camporeale and Burgess2017) supports this interpretation, which is also compatible with earlier results from hybrid and full kinetic two-dimensional simulations reported by Cerri et al. (Reference Cerri, Califano, Jenko, Told and Rincon2016) and by Grošelj et al. (Reference Grošelj, Cerri, Navarro, Willmott, Told, Loureiro, Califano and Jenko2017). If we recall that the Hall term $\epsilon _H$ manifests itself through compressible activity, the cascade rates are expected to be affected in the fluid simulations; while in the HVM simulation, this compressible activity is suppressed through in-plane Landau damping and ion-cyclotron resonances, which re-inject in the system incompressible Alfvénic-like fluctuations. In this two-dimensional setting, the LF model, in which Landau dissipation is introduced at a quasilinear level, partially quenches compressible fluctuations, but not as efficiently as the HVM one. In addition, we can confirm that in the interval $2 d_p \lesssim \ell \lesssim 10 d_p$ (roughly corresponding to the MHD turbulence range), the scaling of $- 2 \epsilon _\ell$ is compatible with a linear scaling law (black-dashed line) for our three cases of study, as previously reported by Sorriso-Valvo et al. (Reference Sorriso-Valvo, Perrone, Pezzi, Valentini, Servidio, Zouganelis and Veltri2018).

Figure 7.

Cascade rate $-2 \epsilon _\ell$ law as a function of the scale $\ell$, conducted in the HMHD (a), LF (b) and HVM (c) simulations. The contributing terms $\epsilon _Y$ (black-solid line) and $\epsilon _H$ (red-solid line) are presented together with their direct sum (asterisks). The dashed line represents the linear scaling.

Moreover, in figure 8 we can see how Yaglom $-\epsilon _{Y_m}$ and Hall $-\epsilon _{H_n}$ terms are supporting the latter global behaviour. We use open diamonds to picture the positive terms, while negative quantities are plotted with red-filled diamonds. Figure 8(a,d,g,j,m) shows the HMHD simulation results, figure 8(b,e,h,k,n) corresponds to the LF simulation, while the results from the HVM simulation are presented in figure 8(c,f,i,l,o). If we focus on the amplitude and sign of these quantities, we cannot note substantial differences between the individual Hall contributions among the simulations. However, the sign of $-\epsilon _{Y_2}$ (associated with the magnetic energy transported by the velocity field) is passing from negative to positive at $\sim 2 d_p$ in the HVM simulation. This is not seen in the HMHD and LF counterparts.

Figure 8.

Terms $-\epsilon _{Y_1}$ (a–c), $-\epsilon _{Y_2}$ (d–f), $-\epsilon _{Y_3}$ (g–i), $-\epsilon _{H_1}$ (j–l) and $-\epsilon _{H_2}$ (m–o). The components are computed from the HMHD (a,d,g,j,m), LF (b,e,h,k,n) and HVM (c,f,i,l,o) simulations. Open diamonds correspond to the positive sign, while red-filled diamonds show negative values. Linear scaling is plotted with a black-dashed line.

Finally, we use our results to compare the amplitude of each of the Yaglom and Hall terms across the scales. In figure 9, we present these comparisons for the Yaglom terms: $\epsilon _{Y_1} / \epsilon _{Y_3}$ (figure 9a) and $\epsilon _{Y_2} / \epsilon _{Y_3}$ (figure 9b). Comparison of Hall terms $\epsilon _{H_1} / \epsilon _{H_2}$ is shown in figure 9(c). Figure 9(d) shows $\epsilon _{Y} / \epsilon _{H}$. Black-solid line represents the HMHD simulation, blue-dash-dotted line represents the LF results and red-dashed line is for the HVM simulation. From a global point of view, we note that the comparisons computed from the LF and HVM simulations are quite similar, especially for scales $\ell \lesssim 10 d_p$. After this range, the sign of $\epsilon _{Y_3}$ fluctuates with respect to that of $\epsilon _{Y_1}$ and $\epsilon _{Y_2}$. On the other hand, the sign of $\epsilon _{H_1}$, as compared with $\epsilon _{H_2}$, remains negative for all of the scales. Once again, the LF and HVM simulations seem quite similar when comparing the Hall terms. The comparison $\epsilon _{Y} / \epsilon _{H}$ shows that the HVM simulation is useful for study the transition from fluid to kinetic scales ($\sim d_p$) as we can see when comparing this ratio.

Figure 9.

(a,b) Comparison of the Yaglom terms $\epsilon _{Y_1}$ and $\epsilon _{Y_2}$, with respect to $\epsilon _{Y_3}$. (c) Comparison between the Hall terms ($\epsilon _{H_1}$ and $\epsilon _{H_2}$). (d) Comparison between $\epsilon _{Y}$ and $\epsilon _{H}$, where the horizontal grey-dashed line marks the barrier of $\epsilon _{Y}/\epsilon _{H} = 1$. Black-solid line corresponds to the HMHD simulation, blue-dash-dotted line corresponds to the LF simulation and the HVM simulation is represented by the red-dashed line.