# Impact of pressure anisotropy on the cascade rate of Hall magnetohydrodynamic turbulence with biadiabatic ions

Pauline A Simon, Fouad Sahraoui, Sébastien Galtier, Dimitri Laveder, Thierry Passot, Pierre-Louis Sulem

► To cite this version:

Pauline A Simon, Fouad Sahraoui, Sébastien Galtier, Dimitri Laveder, Thierry Passot, et al. Impact of pressure anisotropy on the cascade rate of Hall magnetohydrodynamic turbulence with biadiabatic ions. Physical Review E, 2025, 111 (1), pp.015210. 10.1103/PhysRevE.111.015210. hal-04920638

HAL Id: hal-04920638

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

Submitted on 30 Jan 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Impact of pressure anisotropy on the cascade rate of Hall magnetohydrodynamic turbulence with biadiabatic ions

Pauline A. Simon*

Department of Physics and Astronomy, Queen Mary University of London (QMUL), London E1 4NS, United Kingdom
and Laboratoire de Physique des Plasmas (LPP), CNRS, Observatoire de Paris, Sorbonne Université, Université Paris-Saclay,
École Polytechnique, Institut Polytechnique de Paris, 91120 Palaiseau, France

Fouad Sahraoui and Sébastien Galtier

Laboratoire de Physique des Plasmas (LPP), CNRS, Observatoire de Paris, Sorbonne Université, Université Paris-Saclay,
École Polytechnique, Institut Polytechnique de Paris, 91120 Palaiseau, France

Dimitri Laveder, Thierry Passot, and Pierre-Louis Sulem

Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire J.L. Lagrange,
Boulevard de l'Observatoire, CS 34229, F-06304 Nice Cedex 4, France

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

(Received 7 June 2024; accepted 9 December 2024; published 29 January 2025)

The impact of ion pressure anisotropy on the energy cascade rate of Hall-MHD turbulence with biadiabatic ions and isothermal electrons is evaluated in three-dimensional direct numerical simulations, using the exact (or third-order) law derived by Simon and Sahraoui in 2022. It is shown that pressure anisotropy can enhance or reduce the cascade rate, depending on the scales, in comparison with the prediction of the exact law with isotropic pressure, by an amount that correlates well with pressure anisotropy, $a_p = \frac{\ell_{\perp}}{p_{\parallel}} \neq 1$, that develops in simulations initialized with an isotropic pressure ($a_{p0} = 1$). A simulation with initial pressure anisotropy $a_{p0} = 4$ confirms this trend, exhibiting a stronger impact on the cascade rate, both in the inertial range and at larger scales, close to the forcing scales. Furthermore, a Fourier-based numerical method, to compute exact laws in numerical simulations in the full $(\ell_{\perp}, \ell_{\parallel})$ increment plane, is presented.

DOI: 10.1103/PhysRevE.111.015210

## I. INTRODUCTION

The turbulent state of fluids and plasmas (e.g., stars, nebula, stellar wind, interplanetary medium) is associated with the development of cascades of various quantities (such as energy or helicity) characterized by a constant transfer rate on a range of scales belonging to the inertial domain, where driving and dissipation are negligible. The total energy cascade rate can be used as a proxy of the energy dissipation and, thus, of the heating rate in a turbulent medium. It can be estimated using third-order laws (also often called exact laws) that relate the cascade rate to the increments of the fluid moments, which are more easily measurable experimentally. Such laws can be obtained by assuming only the statistical homogeneity of the fluid. A more specific set of hypotheses (separation of injection and dissipation scales, statistical stationarity, high Reynolds numbers) is usually applied to isolate the nonlinear terms of the third-order law associated with the cascade process in the inertial range [1]. Initiated by von Kármán

and Howarth [2] and Kolmogorov [3] in the framework of incompressible hydrodynamics turbulence, the formalism has been extended to incompressible magnetohydrodynamics (MHD) [4,5], Hall-MHD [6–9], and two-fluid flows [10,11], and then to various compressible regimes characterized by different equations of state [12–23]. Transfer rates of quantities different from the energy have also been considered in the literature [22,24–27].

In the solar wind and other accessible natural plasmas, the third-order laws have been used to estimate the heating rate using *in situ* observations [15,28–34]. Since these laws grow in complexity when including more physical effects (e.g., compressibility, nonideal terms in Ohm’s law, pressure anisotropy, etc.), one must then resort to numerical simulations to grasp the physics involved in each term, both for freely decaying [1,8,9,18,35–37] and driven turbulence [38,39].

Recently, Simon and Sahraoui [23,40] proposed a generic derivation of the third-order law for the total energy based on the internal-energy equation and applied it to the Hall-MHD model with either isotropic or gyrotropic pressures. This work paves the way to more realistic studies of compressible turbulence in weakly collisional plasmas, such as those of the near-Earth space where the pressure is all but isotropic [34,41,42]. In fluid descriptions of plasmas, pressure anisotropy was first included by Chew *et al.* [43] within an MHD model usually referred to as CGL, after the authors’ names, or biadiabatic because the heat fluxes are neglected. In

<footnote>*Contact author: pauline.simon@qmul.ac.uk</footnote>
the linear approximation, this model is known to permit the development of fire hose and mirror-type instabilities [44,45]. In this paper, we propose to evaluate the cascade rate in driven-turbulence simulations of the CGL-Hall-MHD system for a proton-electron plasma, obtained by adding the Hall effect to the CGL model, with the aim being to estimate quantitatively the effect of pressure anisotropy on the turbulent cascade, using the third-order law proposed by Simon and Sahraoui [40]. The study is performed on simulation data used in Ferrand et al. [38], complemented by a new simulation specifically designed for the present work. CGL-Hall-MHD is insufficient to provide a full description of the sub-ion range dynamics in space plasmas, where additional kinetic effects such as finite Larmor radius corrections, Landau damping, and cyclotron resonance may play an important role. However, it can be viewed as a relatively simple model to study the simultaneous influence of complex phenomena, such as compressibility, dispersion, and pressure anisotropy on the energy cascade in the MHD and the subionic ranges of scales.

This paper is organized as follows: Sec. II recalls the main formulas derived in Simon and Sahraoui [40] that are used in the rest of the paper. Section III provides a brief description of the numerical code and simulation setups. Section III B describes the Fourier-based method used to evaluate the third-order law in the numerical simulations. Section IV summarizes the main results of the study, which are discussed in Sec. V.

## II. THIRD-ORDER LAW FOR HALL-MHD TURBULENCE WITH BIADIABATIC IONS AND ISOTHERMAL ELECTRONS

The third-order laws for turbulent plasma described within the MHD and Hall-MHD approximations with a biadiabatic closure are derived in Ref. [40]. However, they should be extended to the regime of isothermal electrons considered in the simulations used in the present study.

The CGL-Hall-MHD simulations at hand are based on the following evolution equations for the plasma density $ρ$, the fluid velocity $v$, the magnetic field $B$, the parallel $p_‖$ and perpendicular $p_⊥$ components of the gyrotropic ion pressure relative to the local magnetic field, and the total (electron and ion) pressure tensor $\bar{P}$, assuming an isothermal electronic pressure $p_e$,

$$ \partial_t \rho + \nabla \cdot (\rho v) = -v_\rho \Delta_\alpha^4 \rho, $$

$$ \partial_t v + v \cdot \nabla v + \frac{1}{\rho} \nabla \cdot \bar{P} - \frac{j}{\rho} \times B = f - v_\rho \Delta_\alpha^4 v, $$

$$ \begin{aligned} \partial_t B - \nabla \times (v \times B) + \frac{1}{R_i} \nabla \times \left( \frac{j}{\rho} \times B \right) - \frac{1}{R_i} \nabla \times \left( \frac{1}{\rho} \nabla p_e \right) \\ = -v_B \Delta_\alpha^4 B, \end{aligned} $$

$$ \partial_t p_‖ + \nabla \cdot (p_‖ v) + 2p_‖ bb : \nabla v = -v_p \Delta_\alpha^4 p_‖, $$

$$ \partial_t p_⊥ + \nabla \cdot (p_⊥ v) + p_⊥ \nabla \cdot v - p_⊥ bb : \nabla v = -v_p \Delta_\alpha^4 p_⊥, $$

$$ \bar{P} = \frac{\beta_0}{2}(p_⊥ + p_e)\bar{I} + \frac{\beta_0}{2}(p_‖ - p_⊥)bb, $$

$$ p_e = \rho. \tag{1} $$

Here, $j = \nabla \times B$ is the current density, $\bar{I}$ is the identity tensor, and $b = \frac{B}{|B|}$ is the unit vector in the direction of the local magnetic field. The density, magnetic field, velocity, and pressures ($p_‖$, $p_⊥$ and $p_e$) are measured in units of the equilibrium density $\rho_0$, mean magnetic field $B_0$, Alfvén speed $c_A = B_0/\sqrt{4\pi\rho_0}$, and equilibrium ion parallel pressure $p_‖p_0$, respectively. Furthermore, $R_i = \frac{L}{d_i} = L\sqrt{4\pi e^2 \rho_0 / (m_p c)}$, where $L$ is the unit length, $d_i$ is the ion inertial length, $e$ is the elementary charge, $c$ is the light speed, and $m_p$ is the proton mass, measures the importance of the Hall term relative to the nonlinearity. The ratio of the parallel ion pressure to the magnetic pressure is given by $\beta_0/2 = p_‖p_0/(\rho_0 c_A^2)$. The time unit is $T = L/c_A$. The term $f$ describes the injection of kinetic energy in the system, while the hyperdissipative terms $v_X \Delta_\alpha^4 X$ induce small-scale dissipation of the generic field $X$. Note that the electronic pressure term, included in Ohm's law, cancels out for isentropic or isothermal electronic flows [46].

The third-order law in the present setting can easily be obtained from the one derived in [40] for the two-point correlation function of total energy:

$$ \mathcal{R} = \frac{1}{2} \langle \frac{1}{2}(\rho' + \rho)v \cdot v' + \frac{1}{2}(\rho' + \rho)v_A \cdot v'_A + \rho u' + \rho' u \rangle, \tag{2} $$

where $v_A = B/\sqrt{\rho}$ is the local Alfvén velocity, and $\rho u = \frac{\beta_0}{4}(2p_⊥ + p_‖) + \rho_e u_e$ is the density of internal energy, including the isothermal electronic contribution $\rho_e u_e = \frac{\beta_0}{2}\rho \ln \rho$. The operator $\langle \rangle$ indicates an ensemble average, which is generally replaced by a spatial average, assuming ergodicity holds. As usual, $\delta\varphi$ refers to the increment of a generic field $\varphi$ between the two positions $x$ and $x' = x + \ell$. Primed quantities are evaluated at the position $x'$. Based on the sole assumption of statistical homogeneity, the “von Karman-Howarth-Monin (KHM) equation” that gives the evolution of $\mathcal{R}$ will read (we will go through the details gradually afterward)

$$ \partial_t \mathcal{R} + \varepsilon_{\text{CGL}} - \varepsilon_F - \varepsilon_D = 0. \tag{3} $$

The nonlinear term $\varepsilon_{\text{CGL}}$ is the total cascade rate and is expected to be equal to the mean injection rate $\varepsilon_F$ in the inertial range of scales of stationary turbulence, where the dissipation $\varepsilon_D$ rate is negligible and the temporal derivative $\partial_t \mathcal{R}$ vanishes. Its expression is given by

$$ -4\varepsilon_{\text{CGL}}(\ell) = \nabla_{\ell} \cdot \mathcal{F}(\ell) + [\mathcal{S}(\ell) + \mathcal{S}(-\ell)] + \varepsilon_{\text{Hall}}(\ell), \tag{4} $$

with

$$ \mathcal{F}(\ell) = \langle [\delta(\rho v) \cdot \delta v + \delta(\rho v_A) \cdot \delta v_A + 2\delta\rho\delta u] \delta v \rangle + \left\langle -\delta\rho\delta\left(\frac{\bar{P}_*}{\rho}\right) \cdot \delta v \right\rangle - \langle [\delta(\rho v_A) \cdot \delta v + \delta(\rho v) \cdot \delta v_A] \delta v_A \rangle, $$
$$
\begin{align}
S(\ell) &= \left\langle \left( \rho \mathbf{v} \cdot \delta \mathbf{v} + \frac{1}{2} \rho \mathbf{v}_A \cdot \delta \mathbf{v}_A - \frac{1}{2} \mathbf{v}_A \cdot \delta(\rho \mathbf{v}_A) + 2 \rho \delta u \right) \nabla' \cdot \mathbf{v}' \right\rangle + \left\langle \left[-2 \rho \mathbf{v} \cdot \delta \mathbf{v}_A - \rho \mathbf{v}_A \cdot \delta \mathbf{v} + \delta(\rho \mathbf{v}) \cdot \mathbf{v}_A\right] \nabla' \cdot \mathbf{v}'_A \right\rangle \nonumber \\
&\quad + \left\langle -2 \rho \delta \left(\frac{\bar{P}}{\rho}\right) : \nabla' \mathbf{v}' \right\rangle + \left\langle \left( \delta \rho \frac{\bar{P}_*}{\rho} \cdot \mathbf{v} - \rho \delta \left(\frac{\bar{P}_*}{\rho}\right) \cdot \mathbf{v} \right) \cdot \frac{\nabla' \rho'}{\rho'} \right\rangle,
\end{align}
$$

$$
\epsilon_{\text{Hall}}(\ell) = 2 \frac{1}{R_i} \left\langle \overline{\mathbf{j} \times \mathbf{v}_A} \times \delta \mathbf{v}_A - \delta \left( \frac{\overline{\mathbf{j}}}{\rho} \times \mathbf{v}_A \right) \times \overline{\rho \mathbf{v}_A} \right\rangle + \frac{1}{R_i} \left\langle \frac{1}{2} (\delta \rho \mathbf{v}_A' \cdot \mathbf{v}_A) \nabla' \cdot \frac{\mathbf{j}'}{\rho'} - \left( \delta \rho \frac{\mathbf{j}'}{\rho'} \cdot \mathbf{v}_A \right) \nabla' \cdot \mathbf{v}_A' \right\rangle, \quad (5)
$$

where $\overline{P}_*$ is the total (thermal and magnetic) pressure tensor. The spatial derivative $\nabla'$ is performed relative to the variable $\mathbf{x}'$, while $\nabla_\ell$ is the spatial derivative operator in the scale space.

The contributions to the cascade rate [Eq. (5)] can be separated into the divergence of the flux density $\mathcal{F}$ constituted of structure functions, source $S$ and Hall ($\epsilon_{\text{Hall}}$) terms. To analyze them, we define the isotropic, $p = \frac{\beta_0}{6}(2p_\perp + p_\parallel) + \frac{\beta_0}{2}\rho$, and anisotropic, $\overline{\mathbf{\Pi}} = \overline{\mathbf{P}} - p\overline{\mathbf{I}}$, parts of the pressure. Therefore, the cascade rate (5) can be split into two parts, $\epsilon_{\text{iso}}$ and $\epsilon_{\overline{\mathbf{\Pi}}}$, the latter depending on the anisotropic part of the pressure tensor $\overline{\mathbf{\Pi}}$, while the former does not. One writes

$$ -4\epsilon_{\overline{\mathbf{\Pi}}} = \nabla_{\ell} \cdot \mathcal{F}_{\Pi} + S_{\Pi1} + S_{\Pi2} + S_{\Pi3}, \qquad (6) $$

with

$$ \mathcal{F}_{\Pi}(\ell) = -\left\langle \delta \rho \delta \left(\frac{\overline{\mathbf{\Pi}}}{\rho}\right) \cdot \delta \mathbf{v} \right\rangle, $$

$$ S_{\Pi1}(\ell) = [S(\ell) + S(-\ell)]_{\Pi1} = \left\langle 2\rho'\delta\left(\frac{\overline{\mathbf{\Pi}}}{\rho}\right) : \nabla\mathbf{v} - 2\rho\delta\left(\frac{\overline{\mathbf{\Pi}}}{\rho}\right) : \nabla'\mathbf{v}' \right\rangle, $$

$$ S_{\Pi2}(\ell) = [S(\ell) + S(-\ell)]_{\Pi3} = \left\langle \rho'\delta\left(\frac{\overline{\mathbf{\Pi}}}{\rho}\right) \cdot \mathbf{v}' \cdot \frac{\nabla\rho}{\rho} - \rho\delta\left(\frac{\overline{\mathbf{\Pi}}}{\rho}\right) \cdot \mathbf{v} \cdot \frac{\nabla'\rho'}{\rho'} \right\rangle, $$

$$ S_{\Pi3}(\ell) = [S(\ell) + S(-\ell)]_{\Pi2} = \left\langle \delta\rho \frac{\overline{\mathbf{\Pi}}}{\rho} \cdot \mathbf{v} \cdot \frac{\nabla'\rho'}{\rho'} - \delta\rho \frac{\overline{\mathbf{\Pi}}'}{\rho'} \cdot \mathbf{v}' \cdot \frac{\nabla\rho}{\rho} \right\rangle. \quad (7) $$

When computing third-order laws in numerical simulation, it is important to assess their level of accuracy. For this purpose, we will verify Eq. (3), by computing $\partial_t \mathcal{R}$, $\epsilon_F$, and $\epsilon_D$. In our simulations with the forcing and the hyperdissipative terms, $\epsilon_F$ and $\epsilon_D$ are given by

$$ \epsilon_F = \frac{1}{4} \left\langle (\rho + \rho')(\mathbf{v} \cdot \mathbf{f}' + \mathbf{v}' \cdot \mathbf{f}) \right\rangle, \qquad (8) $$

$$
\begin{align}
\epsilon_D = -\frac{\nu_v}{4} &\left\langle (\rho + \rho')(\mathbf{v} \cdot \Delta_{\alpha}^4 \mathbf{v}' + \mathbf{v}' \cdot \Delta_{\alpha}^4 \mathbf{v}) \right\rangle - \frac{\nu_{\rho}}{4} \left\langle (\Delta_{\alpha}^4 \rho + \Delta_{\alpha}^{4' \rho'}) \mathbf{v} \cdot \mathbf{v}' \right\rangle \nonumber \\
& - \frac{\nu_B}{4} \left\langle (\rho + \rho') \left( \frac{1}{\sqrt{\rho'}} \mathbf{v}_A \cdot \Delta_{\alpha}^{4'} (\sqrt{\rho'} \mathbf{v}'_A) + \frac{1}{\sqrt{\rho}} \mathbf{v}'_A \cdot \Delta_{\alpha}^4 (\sqrt{\rho} \mathbf{v}_A) \right) \right\rangle + \frac{\nu_{\rho}}{8} \left\langle (\rho - \rho') \left( \frac{\Delta_{\alpha}^4 \rho}{\rho} - \frac{\Delta_{\alpha}^{4' \rho'}}{\rho'} \right) \mathbf{v}_A \cdot \mathbf{v}'_A \right\rangle \nonumber \\
& - \frac{\nu_{\rho} \beta_0}{8} \left\langle \frac{\rho}{\rho'} \Delta_{\alpha}^{4'} (2p_{\perp} + p_{\parallel}) + \frac{\rho'}{\rho} \Delta_{\alpha}^4 (2p_{\perp} + p_{\parallel}) \right\rangle - \frac{\nu_{\rho} \beta_0}{8} \left\langle \left( \frac{\Delta_{\alpha}^4 \rho}{\rho} - \frac{\Delta_{\alpha}^{4' \rho'}}{\rho'} \right) \left( \frac{\rho}{\rho'} (2p_{\perp} + p_{\parallel}) - \frac{\rho'}{\rho} (2p_{\perp} + p_{\parallel}) \right) \right\rangle \nonumber \\
& - \frac{\nu_{\rho} \beta_0}{4} \left\langle \left( \frac{\rho'}{\rho} + \ln \rho' \right) \Delta_{\alpha}^4 \rho + \left( \frac{\rho}{\rho'} + \ln \rho \right) \Delta_{\alpha}^{4' \rho'} \right\rangle. \tag{9}
\end{align}
$$

In fact, when computed numerically, the left-hand side of Eq. (3) is not strictly zero but amounts to a value $\zeta$ that defines the numerical errors inherent to the computation of the various terms.

### III. NUMERICAL SETUP

#### A. Simulation model

The numerical code used here is based on Eq. (1). Assuming periodic boundary conditions, the Cartesian three-dimensional space variables $\{x, y, z\}$ are discretized via a Fourier pseudospectral method, where most of the aliasing is

suppressed by spectral truncation at 2/3 of the maximal wave number in each direction. The time stepping is performed with a third-order low-storage Runge-Kutta scheme with a prescribed time step $\delta t$.

In the driven turbulence simulations considered here, a Langevin forcing $f$ is added in the momentum equation [see Eq. (1)]. In Fourier space, this forcing appears as a Dirac distribution at the smallest wave vectors of the simulation domain and reflects the injection of balanced (in the sense that the same amount of energy is carried by forward- and backward-propagating waves) dispersive Alfvén waves with
<table><caption>TABLE I. Parameters of the various runs: $k_{\perp}^{(0)}d_i$, $\theta_i$, $E_{max}$, $A_f$ are forcing parameters (injection wave number, angle, threshold, and amplitude, respectively), while $v_v = v_B$, $v_\rho$, $v_p$ are hyperdissipation coefficients arising in the dynamical equations for the velocity, the magnetic field, the density, and the pressures, respectively. The anisotropy parameter $\alpha$ enters the definition of the anisotropic Laplacian used in the simulations. The parameters $\delta t$ and $t_I$ refer to the time step and to the time at which the analysis was performed, while $a_{p0}$ measures the initial pressure anisotropy.</caption><thead><tr><th>Name</th><th>CGL1</th><th>CGL2</th><th>CGL3</th><th>CGL3B</th><th>CGL5</th><th>CGL6</th></tr></thead><tbody><tr><td>Resolution</td><td>512<sup>3</sup></td><td>512<sup>3</sup></td><td>512<sup>2</sup> × 1024</td><td>512<sup>2</sup> × 1024</td><td>512<sup>2</sup> × 1024</td><td>512<sup>2</sup> × 1024</td></tr><tr><td>k<sub>⊥</sub><sup>(0)</sup>d<sub>i</sub></td><td>0.045</td><td>0.045</td><td>0.5</td><td>0.5</td><td>0.147</td><td>0.147</td></tr><tr><td>θ<sub>i</sub></td><td>83°</td><td>75°</td><td>75°</td><td>75°</td><td>75°</td><td>75°</td></tr><tr><td>E<sub>max</sub></td><td>1.6 × 10<sup>-2</sup></td><td>1.6 × 10<sup>-2</sup></td><td>4.5 × 10<sup>-2</sup></td><td>1.125 × 10<sup>-2</sup></td><td>1.6 × 10<sup>-2</sup></td><td>1.6 × 10<sup>-2</sup></td></tr><tr><td>A<sub>f</sub></td><td>1.0 × 10<sup>-3</sup></td><td>1.0 × 10<sup>-3</sup></td><td>8.0 × 10<sup>-3</sup></td><td>4.0 × 10<sup>-3</sup></td><td>3.0 × 10<sup>-3</sup></td><td>3.0 × 10<sup>-3</sup></td></tr><tr><td>v<sub>v</sub> = v<sub>B</sub></td><td>7.35 × 10<sup>-8</sup></td><td>7.35 × 10<sup>-8</sup></td><td>4.0 × 10<sup>-14</sup></td><td>1.0 × 10<sup>-14</sup></td><td>3.0 × 10<sup>-11</sup></td><td>3.0 × 10<sup>-11</sup></td></tr><tr><td>v<sub>ρ</sub></td><td>0</td><td>0</td><td>1.6 × 10<sup>-14</sup></td><td>1.0 × 10<sup>-14</sup></td><td>0</td><td>0</td></tr><tr><td>v<sub>p</sub></td><td>7.35 × 10<sup>-9</sup></td><td>7.35 × 10<sup>-9</sup></td><td>1.6 × 10<sup>-14</sup></td><td>1.0 × 10<sup>-14</sup></td><td>3.0 × 10<sup>-12</sup></td><td>3.0 × 10<sup>-12</sup></td></tr><tr><td>α</td><td>80</td><td>10</td><td>2.5</td><td>2.5</td><td>6</td><td>5</td></tr><tr><td>δt</td><td>6.25 × 10<sup>-2</sup></td><td>5 × 10<sup>-2</sup></td><td>2 × 10<sup>-4</sup></td><td>3 × 10<sup>-4</sup></td><td>5 × 10<sup>-3</sup></td><td>5 × 10<sup>-3</sup></td></tr><tr><td>t<sub>I</sub></td><td>6700</td><td>12900</td><td>361</td><td>410</td><td>12905</td><td>2730</td></tr><tr><td>a<sub>p0</sub></td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>4</td></tr></tbody></table>

amplitudes $A_f$, random phases, and a wave vector making an angle $\theta_i$ with the mean magnetic field. For two simulations (CGL1 and CGL2), forced in the MHD scales, the injected waves are close to MHD Alfvén waves. After some time, a quasistationary regime is established in which the sum of the perpendicular kinetic and magnetic energies fluctuates around a prescribed level. The forcing is turned on when the sum of the energies of the perpendicular kinetic and magnetic fluctuations reaches a floor level given by 93% of the prescribed maximal energy $E_{max}$ and is turned off when this latter level is reached. The internal energy $u$ can, however, continue to increase, with $p_\perp$ increasing and $p_\parallel$ decreasing (not shown), but on a typical time scale longer than the characteristic time of the turbulence dynamics. In some runs (e.g., CGL3), the averaged pressures reach a quasistationary level, while in others (e.g., CGL3B), they continue to evolve, but, due to the separation in the time scales, the analysis of the transfer rate can be carried out safely.

Hyperdissipation in the form of terms depending on $\Delta_\alpha^4$ and a coefficient is also added in the dynamical equations [Eq. (1)] to smooth the solution at the smallest scales. The coefficients are adjusted to optimize the scale separation between injection and dissipation. Since the simulations develop spatial anisotropies with respect to the mean magnetic field taken in the $z$ direction, an anisotropy parameter $\alpha$ is introduced in the Laplacian operator $\Delta_\alpha = \partial_x^2 + \partial_y^2 + \alpha \partial_z^2$.

A summary of the parameter setup for the reported simulations is given in Table I. Except for run CGL1 for which the driving angle is $\theta = 83^\circ$, this angle is taken equal to $75^\circ$ in all other simulations. In run CGL1 and CGL2, the driving acts at a transverse wave number $k_\perp^{(0)}$ located in the MHD range ($k_\perp^{(0)}d_i = 0.045$), while in runs CGL3 and CGL3B, it acts in the close Hall range ($k_\perp^{(0)}d_i = 0.5$). The difference between CGL3 and CGL3B originates from the prescribed energy level $E_{max}$ that is four times larger in the former. Run CGL5 is initialized as CGL6, except for the pressures, with the driving located in the MHD range close to the Hall

range ($k_\perp^{(0)}d_i = 0.145$). We can estimate a Reynolds number characterizing the prevalence of the nonlinear process over the dissipative one as $R_e \sim vL^7/v_v$ because of the use of a $(-\Delta)^4$ hyperdissipation, $L \sim 1/(k_\perp^{(0)}d_i)$ denoting the scale of energy injection. From the values of $v_v$ and $k_\perp^{(0)}d_i$, and knowing that the velocity amplitude is about $v \approx 0.1$ for all simulations, we get huge values from $R_e \approx 10^{14}$ to $R_e \approx 10^{24}$ according to the simulations, despite the moderate scale separation between injection and dissipation and thanks to the use of hyperdissipation.

The initial conditions for the different fields are $\rho_0 = 1$, $v_0 = [0, 0, 0]$, and $B_0 = [0, 0, 1]$. The ion beta $\beta_{i0}$ and the dimensionless ion inertial length $d_i$ are initialized to one. The pressure anisotropy $a_p = p_\perp/p_\parallel$ is initialized at $a_{p0} = 1$ for all the runs, except CGL6 where it is set equal to four. The runs CGL1, CGL2, and CGL3 were studied by Ferrand et al. [38] in a different framework. A three-dimensional (3D) rendering of the magnitude of the magnetic fluctuations, together with a few total magnetic field lines, is shown in Fig. 1.

The data set used for each simulation is extracted at the time $t_I$ given in Table I. At this time, the perpendicular magnetic energy spectrum has stabilized in a quasistationary regime and displays spectral exponent close to the expected values $-5/3$ and $-7/3$, signatures of an inertial range in the MHD and subion ranges of scales, respectively (see Fig. 2). These results indicate that the energy injected at the largest scales of the simulations by the forcing (that acts almost continuously) cascades through an inertial range extending on one to two decades, down to the smallest scales where the hyperdissipation can act efficiently.

## B. Pseudospectral calculation of the structure functions

As previously mentioned, the expression of a cascade rate can be written in the general form

$$-4\epsilon(\ell) = \nabla_{\ell} \cdot \mathcal{F}(\ell) + [\tilde{S}(\ell) + \tilde{S}(-\ell)],$$
<figure><img src="image_3.png" /><figcaption>FIG. 1. 3D rendering in the physical space of the magnetic field for CGL5. The color reflects the amplitude of the magnetic fluctuations $[B_x^2 + B_y^2 + (B_z - B_0)^2]^{1/2}$, while the white lines are the global field lines.</figcaption></figure>

where $\epsilon$ is the cascade rate and $\tilde{S}$ collects all the terms not included in the divergence of $\mathcal{F}$.

A straightforward method to compute it on the whole scale grid is to code a loop on all the possible separations $\ell$ whose components in each direction are multiple of the mesh size in this direction, up to half the corresponding linear extension of the domain in this direction, to compute spatial derivatives with a four-point difference scheme for good accuracy, and to approximate the ensemble average by a spatial average over all the collocation points of the separation domain. However, such a procedure is extremely time- and memory-consuming. Consequently, in previous works [1,18,38], the computation was performed on only a fraction of scales associated with wisely selected directions in the separation space. Here, we introduce a different method based on the equivalence between the two-point correlation function and a convolution product. Thanks to this remark, we can compute the whole scale

<figure><img src="image_4.png" /><figcaption>FIG. 2. Power spectra of the perpendicular magnetic field for the various simulations. The expected slope in the inertial (Hall) range is indicated in cyan (gray). The spectra of CGL1 and CGL2 are shifted down by a factor of 20 for clarity.</figcaption></figure>

grid at once. Indeed, for two fields A and B, the correlation product $\langle A' \cdot B \rangle = \int A(\mathbf{x} + \mathbf{l})B(\mathbf{x})d\mathbf{x}$ can be viewed as the convolution product $A(\mathbf{l}) * B(-\mathbf{l}) = \int A(\mathbf{x}) \cdot B[-(\mathbf{l} - \mathbf{x})]d\mathbf{x}$ that can easily be done in Fourier space as a simple product of the Fourier transform of A and of the complex conjugate of the Fourier transform of B. This protocol is possible for all terms in the source part $\tilde{S}$. For the terms entering $\mathcal{F}$, the nonlinearity of which is higher than quadratic, the procedure is to be performed iteratively. The outcome is a 3D array in the scale space.

Assuming that the scale space is axisymmetric with respect to the mean magnetic field and that all the considered statistical quantities are even in $\ell$, two-dimensional (2D) maps in function of $\ell_{\perp}$ and $\ell_{\parallel} > 0$ are obtained by angular averaging in each perpendicular plane. These full 2D plots provide information about the cascade anisotropy in turbulence experiments [47] and numerical simulations [48]. To describe the cascade in the perpendicular direction, one may consider the one-dimensional (1D) profiles, functions of $\ell_{\perp}$ only, which are obtained by restricting the considered quantity to the plane $\ell_{\parallel} = 0$. This procedure mirrors the integration in the Fourier space over all $k_{\parallel}$ to obtain the spectral distribution in $k_{\perp}$. Reciprocally, integration over all the parallel scales $\ell_{\parallel}$ is equivalent, in Fourier space, to consider the specific $k_{\parallel} = 0$ plane. Note that, for the sake of simplicity, we keep using the same notation $\epsilon$ for the 2D cascade rate $\epsilon(\ell_{\perp}, \ell_{\parallel})$ and the cascade rate in the perpendicular direction (in brief, perpendicular cascade rate) $\epsilon(\ell_{\perp}) = \epsilon(\ell_{\perp}, \ell_{\parallel} = 0)$ because the context prevents any ambiguity.

The time derivative needed to compute $\zeta$ is obtained using a two-point finite-difference scheme around the selected time $t_I$ (see Table I).

## IV. RESULTS

In this section, the results of the computation of the third-order law on the runs CGL1-CGL6 are analyzed. We recall that CGL1 and CGL2 cover mostly MHD scales, while CGL3 targets the smaller (Hall) scales. CGL3B is similar to CGL3 with a lower energy content, CGL5 was designed to cover the transition between the MHD and Hall ranges by prescribing a forcing at an intermediate scale, CGL6 is similar to CGL5 except for the initial pressure anisotropy. Combined, these simulations allow us to cover a range of scales extending on more than 2.5 decades, and spanning both the MHD and the subion scales. We start our analysis by considering the cascade in the perpendicular direction before turning to the full 2D maps.

### A. Global contributions to the KHM equation

The 1D perpendicular profiles ($\ell_{\parallel} = 0$) of the various terms entering Eq. (3) and of the numerical error $\zeta$ are shown for CGL1, CGL2, CGL3 and CGL6 in Fig. 3. In the increment space, the forcing term $\epsilon_F$ is nearly constant down to the smallest scales, consistent with a forcing in the form of a Dirac distribution in Fourier space. The term $\partial_t \mathcal{R}$ reflects the instantaneous time fluctuations that, as in decaying turbulence [1], do not prevent the establishment of an extended inertial range. The properties of $\epsilon_F$ and $\partial_t \mathcal{R}$ differ in the various simulations:
<figure><img src="image_5.png" /><figcaption>FIG. 3. 1D perpendicular profiles of the various terms of Eq. (3) for CGL1, CGL2, CGL3, and CGL6. Here and in the similar graphs presented in the following, solid lines refer to positive quantities and dashed lines to the absolute value of negative quantities. Note that the blue line is almost indistinguishable from the orange one, consistent with Eq. (3) and the low level of $\zeta$. The gray dot-dashed line corresponds to a -2 slope.</figcaption></figure>

for CGL1 and CGL2 (and CGL3B, not shown), both terms are positive at most of the scales, for CGL3 (and CGL5, not shown), they are mostly negative, while for CGL6, they have opposite signs. Furthermore, the instantaneous behavior of $\partial_t R$ seems to compensate that of the forcing in order to keep the inertial range quasisteady in time (not shown). Part of the injected energy is, in fact, directly converted to internal energy [38], which leads to enhancement and fluctuations of $\epsilon_{CGL}$ at the largest scales, as demonstrated in the following.

The dissipation term $\epsilon_D$ has an envelope consistent with a power law of exponent -2, at the smallest scales. This reflects a divergence of the two-point correlation function of a variable that has a power-law spectrum $k^{-3}$ or steeper, which is the case of the spectra shown in Fig. 2 that steepen significantly at large wave number because of hyperdissipation. This issue can be overcome by using higher-order correlation functions [49]. We will not further discuss the dissipation term in this paper considering that it represents about 1% of $\epsilon_{CGL}$ in the inertial range and thus does not significantly affect the cascade rate, but at the smallest scales.

We observe that $\epsilon_{CGL}$ (blue) coincides with $-\partial_t R + \epsilon_F + \epsilon_D$ (purple) with an error $\zeta$ (brown) between $10^{-7}$ and $10^{-8}$

for CGL1 and CGL2, between $10^{-5}$ and $10^{-6}$ for CGL3 and between $10^{-6}$ and $10^{-7}$ for CGL6. The uncertainty level is thus about 1% to 10% of $\epsilon_{CGL}$, which makes it possible to analyze the various components of $\epsilon_{CGL}$ in these simulations in a relatively reliable way.

Figure 4 shows the resulting 2D maps of the normalized total cascade rate $\epsilon_{CGL}$. The background colors indicate the variation of $\epsilon_{CGL}$ with the positive values in warm colors and the negative ones in cold colors. Superimposed are a few contour lines. As visible on the perpendicular cascade rate (Fig. 3), $\epsilon_{CGL}$ decreases slightly at small scales and varies significantly at the large ones, while displaying a positive plateau at the intermediate $(\ell_{\perp}, \ell_{\|})$ scales. This plateau defines the inertial range. We remark that in the $(\ell_{\perp}, \ell_{\|})$ plane, the inertial range covers a domain whose shape is close to a square because of the logarithmic axes. One also notices that the large-scale variations of $\epsilon_{CGL}$ are quite different in the various simulations: in CGL3B and CGL5 the sign changes, while in CGL3 and CGL6 $\epsilon_{CGL}$ is positive with an enhancement in the forcing range and, for CGL3, a slight decrease at the transition between the inertial and the forcing scales. The behavior of the cascade rate in CGL5 is similar to that in CGL3B in the perpendicular direction, including the sign variations. A
<figure><img src="image_6.png" /><figcaption>FIG. 4. 2D maps of $\epsilon_{CGL}/\max(|\epsilon_{CGL}|)$ for CGL3 (top left), CGL3B (top right), CGL5 (bottom left), and CGL6 (bottom right). An orange-red background indicates a positive nonlinear transfer and a green-blue background, a negative one. Contour lines are based on the linear scale displayed on the left of the color bar.</figcaption></figure>

possible interpretation of these features is presented in the following by analyzing the different components of $\epsilon_{CGL}$.

## B. Isotropic and anisotropic pressure-induced contributions, $\epsilon_{iso}$ and $\epsilon_{\bar{\Pi}}$, to the total cascade rate $\epsilon_{CGL}$

We now analyze, in the various simulations, the partition of the total cascade rate into $\epsilon_{\bar{\Pi}}$, originating from the pressure anisotropy, and the isotropic contribution $\epsilon_{iso}$. The results are shown in Fig. 5 as 1D perpendicular profiles in $\ell_{\perp}$. An estimate of the pressure-anisotropy contribution, through the signed ratio $\epsilon_{\bar{\Pi}}/|\epsilon_{CGL}|$, is given in Fig. 6, using ($\ell_{\perp}$, $\ell_{\|}$) 2D maps.

Figure 5 shows that the contribution of $\epsilon_{\bar{\Pi}}$ to the total cascade rate in the inertial range is smaller ($\epsilon_{\bar{\Pi}} \simeq 10\%$ of $\epsilon_{CGL}$) for CGL1 and CGL2 than for CGL3B (30%) and CGL6 (50%). However, at the largest scales, close to the forcing, this contribution reaches 30%–60% of $\epsilon_{CGL}$ for CGL1-2-3B, and 100% for CGL3-5-6. This increase is also visible in Fig. 6, where the maps become bluer or redder at large scales. The higher level of the full cascade rate for CGL6 compared to the other runs is due to the larger initial pressure anisotropy of this run, which impacts $\epsilon_{\bar{\Pi}}$. One also notices the similarities in the behavior of CGL3 and CGL5. The case of CGL3B will be discussed in Sec. V. As confirmed by CGL5, for which the

forcing scale lies between those of CGL1-2 and of CGL3B-3, the variations of the cascade rate occur mainly at the injection scales.

As seen in Fig. 5, for CGL1, CGL2, and CGL3B, $\epsilon_{\bar{\Pi}}$ is positive at the smallest scales ($\ell_{\perp}/d_i \lesssim 1$) and negative at larger ones. The change of sign occurs in the MHD range, close to $\ell_{\perp} = d_i$ in the perpendicular direction, and close to the forcing scales in the parallel direction (see Fig. 6). This behavior underlines the role of the forcing in the generation of sign fluctuations for $\epsilon_{\bar{\Pi}}$ (and possibly for other components of the cascade rate), rather than a physical effect occurring at the specific scales where the sign change is observed. These changes of sign induce a shallow decrease of the total cascade rate and limit the extension of the inertial range close to the forcing scale. While no change of sign of $\epsilon_{\bar{\Pi}}$ is visible for the perpendicular cascade of CGL3 and CGL5, the cascade rate behaves differently in the parallel direction for these simulations. CGL6 shows no sign change in either direction. For all the simulations, $\epsilon_{iso}$ and $\epsilon_{\bar{\Pi}}$ are significantly larger than the uncertainty level $\zeta$, indicating a negligible effect of the numerical error on the total cascade rate.

The above results highlight the contribution of pressure anisotropy to the cascade rate. The detailed analysis of $\epsilon_{iso}$ and $\epsilon_{\bar{\Pi}}$ is performed in the following subsections.
<figure><img src="image_7.png" /><figcaption>FIG. 5. 1D perpendicular-cascade $\epsilon_{CGL}$ (blue), $\epsilon_{iso}$ (orange), and $\epsilon_{\bar{\Pi}}$ (green) for the various simulations. The gray areas hold for the error level.</figcaption></figure>

### C. Contribution to $\epsilon_{iso}$ induced by density fluctuations

Here, we compare $\epsilon_{iso}$ that collects all the terms independent of the pressure anisotropy in $\epsilon_{CGL}$, as obtained in the present simulations, to the predictions of the incompressible cascade rate computed by Ferrand et al. [38], and with other results based on isotropic compressible models [18].

Figure 7 shows the ($\ell_{\perp}$, $\ell_{\|}$) 2D map of $\epsilon_{iso}$, obtained from CGL2 (left) and CGL3 (right). An inertial range is visible between $\ell_{\perp} \in [0.5, 7]$ and $\ell_{\|} \in [5, 90]$ for CGL2, and between $\ell_{\perp} \in [0.05, 1]$ and $\ell_{\|} \in [0.2, 10]$ for CGL3. At the smallest scales, $\epsilon_{iso}$ decreases due to dissipation and varies in sign and amplitude at the largest scales where the forcing acts. As expected, the inertial range (in the sense of a constant cascade rate) of CGL2 is evidenced in the MHD scales, while for CGL3 it is limited to the Hall scales. CGL1 (not shown) exhibits a broader inertial range in $\ell_{\perp}$ than CGL2, probably due to a more efficient generation of Alfvénic turbulence by a driving acting in a more quasiperpendicular direction ($\theta_i = 83^{\circ}$ for CGL1 and $75^{\circ}$ for CGL2).

These observations also hold for the 1D perpendicular cascade rate displayed in the middle graphs of Fig. 7 (thick blue line) for each simulation. The profile of $\epsilon_{iso}$ is close to that of the incompressible Hall-MHD contribution ($\epsilon_{inc-Hall-MHD}$, thin blue line), obtained with the third-order law used by Ferrand et al. [9]. The MHD ($\epsilon_{MHD}$, orange) and Hall ($\epsilon_{Hall}$, green) components are also close to their incompressible counterparts ($\epsilon_{inc-MHD}$ and $\epsilon_{inc-Hall}$, thin lines): the difference between compressible and incompressible models vary between 1% and 10%, and is slightly more important for CGL3. Such observations are expected since the present simulations are only weakly compressible (for CGL2: $\frac{\text{STD}(\rho)}{\langle\rho\rangle} = 2\%$ for CGL3: $\frac{\text{STD}(\rho)}{\langle\rho\rangle} = 8\%$, where STD refers to the standard deviation). The MHD and Hall contributions behave as expected: at scales $\ell_{\perp} \sim d_i$, the Hall contribution increases and overcomes the decreasing MHD one, thus ensuring a constant total cascade rate.

The various components of $\epsilon_{MHD}$ are displayed in the bottom panels of Fig. 7 for each simulation, with $\epsilon_{MHD}$ decomposed as $\epsilon_{MHD} = \nabla_{\ell} \cdot \mathcal{F}_{MHD} + S_{MHD}$, where

$$ \mathcal{F}_{MHD} = \mathcal{F}_{yaglom} + \mathcal{F}_{comp}, \qquad (10) $$
<figure><img src="image_8.png" /><figcaption>FIG. 6. 2D map of $\epsilon_{\bar{\pi}}/|\epsilon_{CGL}|$ for the different runs. Warm color indicates a positive contribution of pressure anisotropy varying from 0% (yellow) to 100% (red), and green-blue color, a negative contribution, while the purple domain circled by a red contour line corresponds to the range where the sign of $\epsilon_{CGL}$ fluctuates, because of the variations of $\epsilon_{iso}$.</figcaption></figure>

with

$$ \mathcal{F}_{\text{yaglom}} = -\frac{1}{4} \langle [\delta(\rho\mathbf{v}) \cdot \delta\mathbf{v} + \delta(\rho\mathbf{v}_A) \cdot \delta\mathbf{v}_A] \delta\mathbf{v} - [\delta(\rho\mathbf{v}_A) \cdot \delta\mathbf{v} + \delta(\rho\mathbf{v}) \cdot \delta\mathbf{v}_A] \delta\mathbf{v}_A \rangle, $$

$$ \mathcal{F}_{\text{comp}} = -\frac{1}{4} \langle \delta\rho \left[ 2\delta u - \delta \left( \frac{p+p_m}{\rho} \right) \right] \delta\mathbf{v} \rangle, $$
<figure><img src="image_9.png" /><figcaption>FIG. 7. (left) CGL2, (right) CGL3, (top) 2D map (with same graphic conventions as in Fig. 4) of $\epsilon_{iso}$ in the $(\ell_{\perp}, \ell_{\parallel})$ plane. (middle and bottom) 1D perpendicular profiles of the various contributions to $\epsilon_{iso}$. (iddle) comparison of the MHD ($\epsilon_{MHD}$, thick orange) and Hall ($\epsilon_{Hall}$, thick green) contributions to $\epsilon_{iso}$ (thick blue) with their incompressible counterparts (thin curves). (bottom) components of $\epsilon_{MHD}$ (orange), i.e., the flux part (thick blue) and source part (thin blue), and detail of the flux part, i.e., Yaglom part (thick green), and purely compressible part (thin green).</figcaption></figure>

and

$$S_{\text{MHD}} = [S(\ell) + S(-\ell)], \quad (11)$$

with

$$S(\ell) = \left\langle \left( \rho \mathbf{v} \cdot \delta \mathbf{v} + \frac{1}{2} \rho \mathbf{v}_A \cdot \delta \mathbf{v}_A - \frac{1}{2} \mathbf{v}_A \cdot \delta(\rho \mathbf{v}_A) + 2 \rho \delta u \right) \nabla' \cdot \mathbf{v}' \right\rangle + \left\langle \left[ -2 \rho \mathbf{v} \cdot \delta \mathbf{v}_A - \rho \mathbf{v}_A \cdot \delta \mathbf{v} + \delta(\rho \mathbf{v}) \cdot \mathbf{v}_A \right] \nabla' \cdot \mathbf{v}'_A \right\rangle \\ + \left\langle -2 \rho \delta \left( \frac{p}{\rho} \right) : \nabla' \mathbf{v}' \right\rangle + \left\langle \left[ \delta \rho \frac{p+p_m}{\rho} \cdot \mathbf{v} - \rho \delta \left( \frac{p+p_m}{\rho} \right) \cdot \mathbf{v} \right] \cdot \frac{\nabla' \rho'}{\rho'} \right\rangle$$

The quantity $\nabla \cdot \mathcal{F}_{yaglom}$ is the part of $\epsilon_{MHD}$ that survives when formally taking $\rho = \rho_0 = 1$, typical of an incompressible regime, yielding $\epsilon_{inc-MHD}$, while $\nabla \cdot \mathcal{F}_{comp}$ will vanish in this limit. Note that $\nabla \cdot \mathcal{F}_{comp}$ and $S_{MHD}$, referred to as
<figure><img src="image_10.png" /><figcaption>FIG. 8. 1D perpendicular profile of the different contributions to $\epsilon_{\bar{\Pi}}$ (blue) in run CGL5 and CGL6. The source term that does not explicitly depend on the density fluctuations is plotted in green. The gray area reflects the uncertainty level $\zeta$. Plain lines are for positive values and dashed lines are for negative ones.</figcaption></figure>

hybrid and source terms in Andrés et al. [18] and Simon and Sahraoui [23], are expressed here using the structure functions. Similarly to the results reported in isothermal turbulence simulations [18], their contribution to $\epsilon_{\text{MHD}}$ remains negligible (1% to 10% at the smallest scales), an effect possibly originating from the weak compressibility of the present simulations.

Similar behavior of $\epsilon_{\text{iso}}$ is obtained for the other simulations, in agreement with previous results published in the literature, thus validating the contribution to $\epsilon_{\text{CGL}}$ that is independent of the pressure anisotropy.

### D. Detailed analysis of $\epsilon_{\bar{\Pi}}$

All the simulations display similar results about the contributions of the different terms entering $\epsilon_{\bar{\Pi}}$. Since CGL6 is the run that shows the strongest effect of $\epsilon_{\bar{\Pi}}$, we concentrate here on this run and on CGL5 that is initialized with $a_{p0} = 1$, for comparison. The different terms contributing to $\epsilon_{\bar{\Pi}}$ are given in Eq. (7).

Figure 8 shows that the variations in amplitude and sign of $\epsilon_{\bar{\Pi}}$ (blue curve) is dominated by $S_{\Pi1}$ (green) and is only slightly impacted by the other terms whose contribution is less than 10% of $\epsilon_{\bar{\Pi}}$. We recall that $S_{\Pi1}$ is the term that does not explicitly depend on the density fluctuations, and thus survives in the incompressibility limit. Its impact is thus expected to be dominant in the present quasi-incompressible simulations. In order of importance, $\nabla \cdot \mathcal{F}_{\Pi}$ (orange) and $S_{\Pi2}$ (red) come next. These terms depend on the pressure fluctuation, and their relative importance depends on the runs: $|\nabla \cdot \mathcal{F}_{\Pi}| < |S_{\Pi2}|$ for CGL1 and CGL2, $|\nabla \cdot \mathcal{F}_{\Pi}| > |S_{\Pi2}|$ for CGL3 and CGL3B, and $|\nabla \cdot \mathcal{F}_{\Pi}| \sim |S_{\Pi2}|$ for CGL5 and CGL6. The last contribution $S_{\Pi3}$ (purple) that falls below the uncertainty level $\zeta$ nearly at all scales is the only source term that depends on the density fluctuations.

To compare the impact of the averaged pressure anisotropy to the one of its fluctuations, we decompose $S_{\Pi1}$ into

$$S_{\Pi1}^{(\Delta p)} = 2 \left\langle \frac{\Delta p}{\rho} \right\rangle \langle \delta(\mathbf{b}\mathbf{b}) : (\rho' \nabla \mathbf{v} - \rho \nabla' \mathbf{v}') \rangle$$

with

$$\Delta p = p_{\parallel} - p_{\perp}$$

and

$$S_{\Pi1} - S_{\Pi1}^{(\Delta p)}$$

The results are displayed in Fig. 9. For all simulations, the quantity depending on the fluctuations of $\Delta p = p_{\parallel} - p_{\perp}$ dominates the inertial range, while that depending on the mean value dominates the scales closest to the forcing. For CGL1, CGL2, and CGL3B, the change of sign is associated to a sign variation of the contribution depending on the fluctuations of $\Delta p$. Note that, similarly to the cascade rate (5), the pressure always appears divided by the density, e.g., $\Delta p/\rho$. Thus, the observed behavior could reflect the temperature rather than the pressure anisotropy, although they almost coincide in the present simulations, due to the weak density fluctuations.

## V. DISCUSSION

### A. Spectra of the pressure anisotropy

Section IV shows that the pressure-anisotropy contribution to the total cascade rate in the inertial range is led by the fluctuations of the pressure anisotropy. To better interpret this observation, a measure of the ratio of the spectral distributions of the perpendicular and parallel pressures for runs CGL3, CGL3B, CGL5, and CGL6 is presented in Fig. 10. These runs have been chosen because they show the largest contribution of the pressure-anisotropy terms to the cascade rate (CGL3B and CGL6 in the inertial range; CGL3, CGL5, and CGL6 at the largest scales).

We observe a change of pressure anisotropy from $|\hat{P}_{\perp}|(k_{\perp}) > |\hat{P}_{\parallel}|(k_{\perp})$ (see definitions in the figure caption) to the opposite ordering for runs CGL3 and CGL5 at $k_{\perp}d_i \approx 4$. Differently, CGL3B and CGL6 show a monotonic behavior $[|\hat{P}_{\perp}|(k_{\perp}) < |\hat{P}_{\parallel}|(k_{\perp})]$ at nearly all scales (excluding the dissipation range), with the strongest pressure anisotropy at the smallest scales. We acknowledge the difficulty of relating the behavior of this quantity to the complex nonlinear terms of the third-order law that involve pressure anisotropy. Nevertheless, we can speculate that the strong and uniform-in-scale pressure anisotropy observed for runs CGL3B and CGL6 could explain (or at least, is consistent with) the enhanced contribution of $\epsilon_{\bar{\Pi}}$ observed in the inertial range for those
<figure><img src="image_11.png" /><figcaption>FIG. 9. For each simulation, 1D perpendicular profile of the decomposition of $S_{\Pi 1}$ into terms depending or not on the (mean) quantity $\langle \frac{\Delta p}{\rho} \rangle = \langle \frac{p_{\parallel} - p_{\perp}}{\rho} \rangle$.</figcaption></figure>

runs in comparison to the others. On the other hand, the change from a perpendicular-pressure dominated regime to a parallel-pressure dominated one observed for runs CGL3 and CGL5, with potential fluctuations in sign of $\Delta p$, may result in contributions that would statistically cancel out, hence the weak contribution of the anisotropic terms in the cascade rate for those runs.

## B. Plasma-state distribution in the $(a_p, \beta_\parallel)$ plane

The distributions of the plasma representative points in the $(a_p, \beta_\parallel)$ diagram is displayed in Fig. 11, along with the thresholds of the linear fire hose and mirror instabilities (see Appendix). For all the runs, we observe that the representative points are localized in the stable region, except for a few points from CGL6 laying near the mirror-instability threshold. One can notice that the distributions associated with CGL1 to CGL5 are not centered around the initial values $a_{p0} = 1$ and $\beta_{\parallel 0} = 1$ but are shifted to a new regime of parameters.

The representative points from simulations CGL1, CGL2, and CGL3B stay close to $a_{p0} = 1$, while those corresponding to CGL3 and CGL5 shift up to an average value $\langle a_p \rangle = 2$, with, however, a larger standard deviation for the former run as shown in Table II. This migration is probably due to a direct conversion at large scales of part of the injected energy into internal energy, thus enhancing the average temperatures without contributing to $\epsilon_{iso}$, as shown by Ferrand et al. [38] with CGL1, CGL2, and CGL3 simulations. The difference in the migration appears to be consistent with the results reported above about the importance of the contribution of $\epsilon_{\bar{\Pi}}$ to the total cascade rate $\epsilon_{CGL}$ at large scales where the contribution of $\langle \frac{\Delta p}{\rho} \rangle$ is dominant. However, the proportion of points with $a_p < 1$, higher for CGL1, CGL2, and CGL3B, could explain the change of sign of $\epsilon_{\bar{\Pi}}$, as relatively more points with $\frac{\Delta p}{\rho} < 0$ contribute to $S_{\Pi 1} - S_{\Pi 1}^{(\frac{\Delta p}{\rho})}$. More detailed statistical studies using spacecraft data from the solar wind, for instance, will be necessary to disentangle the behavior of
<figure><img src="image_12.png" /><figcaption>FIG. 10. Ratio of the spectral distributions of the perpendicular and parallel pressures in the transverse Fourier plane, $|\hat{P}_{\perp}|/|\hat{P}_{\parallel}|$, where $|\hat{P}_{\perp}|(k_{\perp}) = [\sum_{k_{\parallel}} \sum_{k_{\perp} \le k'_{\perp} \le k_{\perp} + \Delta k_{\perp}} |\hat{p}_{\perp}(k'_{\perp}, k_{\parallel})|^2]^{1/2}$ and $|\hat{P}_{\parallel}|(k_{\perp}) = [\sum_{k_{\parallel}} \sum_{k_{\perp} \le k'_{\perp} \le k_{\perp} + \Delta k_{\perp}} |\hat{p}_{\parallel}(k'_{\perp}, k_{\parallel})|^2]^{1/2}$, plotted as a function of $k_{\perp}$, for runs CGL3 (yellow), CGL3B (green), CGL5 (red), and CGL6 (purple). Here, $\hat{P}_{\perp}$ and $\hat{P}_{\parallel}$ hold for the Fourier transform of the perpendicular and parallel pressures $p_{\perp}$ and $p_{\parallel}$, respectively. For comparison, the colored points on the left part of the figure indicate the ratio of the mean perpendicular and parallel pressures for the corresponding simulations.</figcaption></figure>

the cascade rate and the localization of the plasma state in the $(a_p, \beta_\parallel)$ plane and to validate the correlations reported here.

## VI. CONCLUSION

This paper summarizes the results of a study of the impact of pressure anisotropy on the turbulent cascade rate described by the third-order law obtained in Simon and Sahraoui [40], using six simulations of CGL-Hall-MHD turbulence. The contribution to the cascade rate of the terms independent of the anisotropic component of the pressure tensor has been found to agree with the results of previous studies using incompressible or weakly compressible models with isotropic pressure. In addition, the contribution of the anisotropic component of the pressure tensor reveals

(1) the relevance of this effect even in a quasi-incompressible case.

<figure><img src="image_13.png" /><figcaption>FIG. 11. Distribution of the plasma states in the $(a_p, \beta_\parallel)$ plane with their occurrence rate: occurrence points (edge of the distribution), 50% (inner ellipse) and 99% (outer curves points) in the various simulations. Dash-dotted lines indicate the fire hose (orange) and mirror (blue) linear instability thresholds (taking into account the isothermal electronic pressure).</figcaption></figure>

(2) an impact of the fluctuations of pressure anisotropy (visible both in real and Fourier space) on the cascade rate in the inertial range.

(3) a more important effect close to the forcing scales, linked to the anisotropy of the mean pressure.

(4) a sign variation of the contribution of $\epsilon_{\bar{\Pi}}$ to the cascade rate, potentially linked to the localization of the simulation in the $(\beta_\parallel, a_p)$ diagram.

Even if CGL6, which is initialized with $a_{p0} = 4$, shows a more prominent contribution of the pressure anisotropy term, the impact in the inertial range remains moderate. Extension of this study would require other simulations with different $a_p$ and possibly higher compressibility. Also, the third-order law of Simon and Sahraoui [40] is derived with a particular form of correlation function $\mathcal{R}$. This choice is not unique, and alternative quantities have been considered in the literature [27,35]. The sensitivity of the present conclusions to the retained correlation function is to be analyzed in further studies.

On the other hand, spacecraft data can also be used to conduct further studies, as they give access to a broad range of parameter space $(\beta_\parallel, a_p, \delta\rho/\rho)$. However, computing point-wise (spatial) derivatives involved in the third-order law using multi-spacecraft data (e.g., Cluster or MMS) requires some

<table><caption>TABLE II. Average and standard deviation of density $\rho$ (or electronic pressure), total ionic beta parameter $\beta_i = (2p_\perp + p_\parallel)/(3p_m)$ and parallel beta parameter $\beta_\parallel = p_\parallel/p_m$, pressure anisotropy $a_p$, percentage of data points where $a_p > 1$ or $a_p < 1$, and parameter $-\beta_\parallel^* = \beta_\parallel(a_p - 1)$, entering the terms of $\epsilon_{\bar{\Pi}}$ [see (7)].</caption><thead><tr><th>Name</th><th>CGL1</th><th>CGL2</th><th>CGL3</th><th>CGL3B</th><th>CGL5</th><th>CGL6</th></tr></thead><tbody><tr><td>⟨ρ⟩ ± σ<sub>ρ</sub></td><td>1.00 ± 0.02</td><td>1.00 ± 0.02</td><td>1.00 ± 0.08</td><td>1.00 ± 0.04</td><td>1.00 ± 0.02</td><td>1.00 ± 0.03</td></tr><tr><td>⟨β<sub>i</sub>⟩ ± σ<sub>β<sub>i</sub></sub></td><td>0.99 ± 0.07</td><td>0.99 ± 0.06</td><td>1.13 ± 0.39</td><td>1.03 ± 0.19</td><td>1.07 ± 0.12</td><td>2.98 ± 0.25</td></tr><tr><td>⟨β<sub>parallel</sub>⟩ ± σ<sub>β<sub>parallel</sub></sub></td><td>0.93 ± 0.09</td><td>0.93 ± 0.08</td><td>0.65 ± 0.25</td><td>0.87 ± 0.18</td><td>0.63 ± 0.07</td><td>1.01 ± 0.15</td></tr><tr><td>⟨a<sub>p</sub>⟩ ± σ<sub>a<sub>p</sub></sub></td><td>1.11 ± 0.08</td><td>1.11 ± 0.10</td><td>2.16 ± 0.49</td><td>1.30 ± 0.30</td><td>2.05 ± 0.12</td><td>3.97 ± 0.47</td></tr><tr><td>a<sub>p</sub> > 1(in %)</td><td>91</td><td>88</td><td>99.8</td><td>83</td><td>100</td><td>100</td></tr><tr><td>a<sub>p</sub> < 1(in %)</td><td>9</td><td>12</td><td>0.2</td><td>17</td><td>0</td><td>0</td></tr><tr><td>⟨-β<sub>parallel</sub><sup>*</sup>⟩ ± σ<sub>β<sub>parallel</sub><sup>*</sup></sub></td><td>0.09 ± 0.07</td><td>0.09 ± 0.08</td><td>0.71 ± 0.32</td><td>0.24 ± 0.24</td><td>0.66 ± 0.10</td><td>2.96 ± 0.25</td></tr></tbody></table>
caution because the uncertainties can be too large to allow reliable quantification of the small contributions due to density fluctuations or pressure anisotropy. The nine probes of the upcoming NASA/Helioswarm mission are likely to improve those estimates.

## ACKNOWLEDGMENTS

The work of P.S. was supported by DIM-ACAV + Doctoral fellowship of the Ile-de-France Region (at LPP, 2020-2023) and by UKRI Future Leaders Fellowship MR/W007657/1 (at QMUL, since 2023). P.S. thanks A. Jeandet for his help in the implementation of the data post-processing. Computations were performed at the mesocenter Cholesky of Ecole Polytechnique and mesocentre SIGAMM, hosted by Observatoire de la Côte d'Azur (OCA).

P.S. performed the present study. F.S. and S.G. supervised the work. The numerical simulation setups were designed and data were obtained in collaboration with D. Laveder, T. Passot, and P.L. Sulem who also provided Figs. 1, and 2,

and 10. All authors contributed to discussing the results and writing the paper. The work of P.S. was supported by DIM-ACAV + Doctoral fellowship of the Ile-de-France Region (at LPP, 2020-2023) and by UKRI Future Leaders Fellowship MR/W007657/1 (at QMUL, since 2023).

## APPENDIX: CGL FIRE HOSE AND MIRROR INSTABILITY THRESHOLDS IN THE PRESENCE OF ISOTHERMAL ELECTRONS

To obtain the CGL-MHD-instability thresholds with isothermal electrons for the fire hose and mirror instabilities, we linearize the dimensionless system (1), taking $d_i = 0$ (no Hall effect) and neglecting dissipative and forcing terms. The mean component of each quantity is isolated and indexed by the subscript 0. Furthermore, at equilibrium, $v_0 = 0$ and $\mathbf{B}_0 = [0, 0, B_0]_{x,y,z}$ are assumed. The electronic pressure is supposed isothermal ($\propto \rho$), consistent with the numerical simulations. Then, we apply the transformation to the Fourier space ($\partial_t \rightarrow -i\omega$, $\partial_x \rightarrow ik_\perp$, $\partial_y \rightarrow ik_y$, $\partial_z \rightarrow ik_\parallel$) with $k_y = 0$. The velocity equations become

$$ 
\begin{aligned}
0 &= -\frac{1}{v_{A0}^2}\omega v_x + \frac{\beta_\parallel}{2}a_p k_\perp \frac{p_\perp}{p_{\parallel 0}} + \left[\frac{\beta_\parallel}{2}(1-a_p) - 1\right]k_\parallel \frac{B_x}{B_0} + k_\perp \frac{B_z}{B_0} + \frac{\beta_e}{2}k_\perp \frac{p_e}{p_{e0}}, \\
0 &= -\frac{1}{v_{A0}^2}\omega v_y + \left[\frac{\beta_\parallel}{2}(1-a_p) - 1\right]k_\parallel \frac{B_y}{B_0}, \\
0 &= -\frac{1}{v_{A0}^2}\omega v_z + \frac{\beta_\parallel}{2}k_\parallel \frac{p_\parallel}{p_{\parallel 0}} + \frac{\beta_\parallel}{2}(1-a_p)k_\perp \frac{B_x}{B_0} + \frac{\beta_e}{2}k_\parallel \frac{p_e}{p_{e0}},
\end{aligned}
\quad (A1) $$

with $v_{A0}^2 = B_0^2/(4\pi\rho_0)$, $\beta_\parallel = \beta_0 4\pi p_{\parallel 0}/B_0^2$, $a_p = p_{\perp 0}/p_{\parallel 0}$, $D_i = d_i/\rho_0$, and $\beta_e = \beta_0 4\pi p_{e0}/B_0^2$.

The other equations read

$$ 
\begin{aligned}
0 &= -\omega\rho + \rho_0(k_\perp v_x + k_\parallel v_z), \\
0 &= -\omega B_x - B_0 k_\parallel v_x, \\
0 &= -\omega B_y - B_0 k_\parallel v_y, \\
0 &= -\omega B_z + B_0 k_\perp v_x, \\
0 &= -\omega p_\parallel + p_{\parallel 0}(k_\perp v_x + 3k_\parallel v_z), \\
0 &= -\omega p_\perp + p_{\perp 0}(2k_\perp v_x + k_\parallel v_z), \\
0 &= -\omega p_e + p_{e0}(k_\perp v_x + k_\parallel v_z).
\end{aligned}
\quad (A2) $$

Injecting (A2) into (A1) gives

$$ 
\begin{aligned}
0 &= \left[-\frac{1}{v_{A0}^2}\omega^2 + \left(2\frac{\beta_\parallel}{2}a_p + 1 + \frac{\beta_e}{2}\right)k_\perp^2 - \left(\frac{\beta_\parallel}{2}(1-a_p) - 1\right)k_\parallel^2\right]v_x + \left(\frac{\beta_\parallel}{2}a_p + \frac{\beta_e}{2}\right)k_\perp k_\parallel v_z, \\
0 &= \left[-\frac{1}{v_{A0}^2}\omega^2 - \left(\frac{\beta_\parallel}{2}(1-a_p) - 1\right)k_\parallel^2\right]v_y, \\
0 &= \left(\frac{\beta_\parallel}{2}a_p + \frac{\beta_e}{2}\right)k_\perp k_\parallel v_x + \left[-\frac{1}{v_{A0}^2}\omega^2 + \left(3\frac{\beta_\parallel}{2} + \frac{\beta_e}{2}\right)k_\parallel^2\right]v_z.
\end{aligned}
\quad (A3) $$

The determinant of the system (A3) involves two factors. The first one gives the dispersive relation of the incompressible Alfvén-wave, polarized such as $v_x = v_z = 0$ [i.e., $\omega = \pm k_\parallel v_{A0}\{1 - [\beta_\parallel(1 - a_p)/2]\}^{1/2}$]. The Alfvén wave is affected by the fire hose instability, which appears if

$$ 1 - \frac{\beta_\parallel}{2}(1 - a_p) < 0. $$
The second factor gives the dispersive relation of the magnetosonic waves polarized such as $v_y = 0$ and takes the form of a second-order polynomial $\omega^4 + A_1\omega^2 + A_0 = 0$ with

$$A_1 = -v_{A0}^2 \left[ \left( 2 \frac{\beta_{\parallel}}{2} a_p + 1 + \frac{\beta_e}{2} \right) k_{\perp}^2 + \left( \frac{\beta_{\parallel}}{2} (2 + a_p) + 1 + \frac{\beta_e}{2} \right) k_{\parallel}^2 \right], \\ A_0 = v_{A0}^2 k_{\parallel}^2 \left\{ \left[ \frac{\beta_{\parallel}}{2} \left( \frac{\beta_{\parallel}}{2} a_p (6 - a_p) + 3 \right) + \frac{\beta_e}{2} \left( 3 \frac{\beta_{\parallel}}{2} + 1 \right) \right] k_{\perp}^2 - \left( 3 \frac{\beta_{\parallel}}{2} + \frac{\beta_e}{2} \right) \left( \frac{\beta_{\parallel}}{2} (1 - a_p) - 1 \right) k_{\parallel}^2 \right\}.$$

Assuming that $k_{\perp} \gg k_{\parallel}$, $A_0$ gives the mirror threshold altered by the isothermal electronic pressure:

$$\frac{\beta_{\parallel}}{2} a_p (6 - a_p) + 3 + \frac{\beta_e}{\beta_{\parallel}} \left( 3 \frac{\beta_{\parallel}}{2} + 1 \right) = 0.$$

If $\beta_e = 0$, the threshold for the CGL-MHD [44] is recovered. In our simulations, $\beta_e = 1$.

[1] R. Ferrand, F. Sahraoui, S. Galtier, N. Andrés, P. Mininni, and P. Dmitruk, An in-depth numerical study of exact laws for compressible Hall magnetohydrodynamic turbulence, Astrophys. J. 927, 205 (2022).

[2] T. von Kármán and L. Howarth, On the statistical theory of isotropic turbulence, Proc. R. Soc. Lond. A Math. Phys. Sci. 164, 192 (1938).

[3] A. N. Kolmogorov, Dissipation of energy in locally isotropic turbulence, Proc. R. Soc. Lond. A 434, 15 (1991).

[4] H. Politano and A. Pouquet, Dynamical length scales for turbulent magnetized flows, Geophys. Res. Lett. 25, 273 (1998).

[5] H. Politano and A. Pouquet, von Kármán-Howarth equation for magnetohydrodynamics and its consequences on third-order longitudinal structure and correlation functions, Phys. Rev. E 57, R21 (1998).

[6] S. Galtier, von Kármán-Howarth equations for Hall magnetohydrodynamic flows, Phys. Rev. E 77, 015302(R) (2008).

[7] S. Banerjee and S. Galtier, An alternative formulation for exact scaling relations in hydrodynamic and magnetohydrodynamic turbulence, J. Phys. A: Math. Theor. 50, 015501 (2017).

[8] P. Hellinger, A. Verdini, S. Landi, L. Franci, and L. Matteini, von Kármán-Howarth equation for Hall magnetohydrodynamics: hybrid simulations, Astrophys. J. 857, L19 (2018).

[9] R. Ferrand, S. Galtier, F. Sahraoui, R. Meyrand, N. Andrés, and S. Banerjee, On exact laws in incompressible Hall magnetohydrodynamic turbulence, Astrophys. J. 881, 50 (2019).

[10] N. Andrés, S. Galtier, and F. Sahraoui, Exact scaling laws for helical three-dimensional two-fluid turbulent plasmas, Phys. Rev. E 94, 063206 (2016).

[11] N. Andrés, P. D. Mininni, P. Dmitruk, and D. O. Gómez, von Kármán-Howarth equation for three-dimensional two-fluid plasmas, Phys. Rev. E 93, 063202 (2016).

[12] S. Galtier and S. Banerjee, Exact relation for correlation functions in compressible isothermal turbulence, Phys. Rev. Lett. 107, 134501 (2011).

[13] S. Banerjee and S. Galtier, Exact relation with two-point correlation functions and phenomenological approach for compressible magnetohydrodynamic turbulence, Phys. Rev. E 87, 013019 (2013).

[14] S. Banerjee and S. Galtier, A Kolmogorov-like exact relation for compressible polytropic turbulence, J. Fluid Mech. 742, 230 (2014).

[15] S. Banerjee, L. Z. Hadid, F. Sahraoui, and S. Galtier, Scaling of compressible magnetohydrodynamic turbulence in the fast solar wind, Astrophys. J. Lett. 829, L27 (2016).

[16] N. Andrés and F. Sahraoui, Alternative derivation of exact law for compressible and isothermal magnetohydrodynamics turbulence, Phys. Rev. E 96, 053205 (2017).

[17] S. Banerjee and A. G. Kritsuk, Exact relations for energy transfer in self-gravitating isothermal turbulence, Phys. Rev. E 96, 053116 (2017).

[18] N. Andrés, F. Sahraoui, S. Galtier, L. Z. Hadid, P. Dmitruk, and P. D. Mininni, Energy cascade rate in isothermal compressible magnetohydrodynamic turbulence, J. Plasma Phys. 84, 21 (2018).

[19] N. Andrés, S. Galtier, and F. Sahraoui, Exact law for homogeneous compressible Hall magnetohydrodynamics turbulence, Phys. Rev. E 97, 013204 (2018).

[20] S. Banerjee and A. G. Kritsuk, Energy transfer in compressible magnetohydrodynamic turbulence for isothermal self-gravitating fluids, Phys. Rev. E 97, 023107 (2018).

[21] S. Banerjee and N. Andrés, Scale-to-scale energy transfer rate in compressible two-fluid plasma turbulence, Phys. Rev. E 101, 043212 (2020).

[22] P. Hellinger, A. Verdini, S. Landi, L. Franci, E. Papini, and L. Matteini, On cascade of kinetic energy in compressible hydrodynamic turbulence, arXiv:2004.02726 (2020).

[23] P. Simon and F. Sahraoui, General exact law of compressible isentropic magnetohydrodynamic flows: theory and spacecraft observations in the solar wind, Astrophys. J. 916, 49 (2021).

[24] O. Soulard and A. Briard, Cascade of circulicity in compressible turbulence, Phys. Rev. Fluids 7, 124604 (2022).

[25] P. Hellinger, E. Papini, A. Verdini, S. Landi, L. Franci, L. Matteini, and V. Montagud-Camps, Spectral transfer and Kármán-Howarth-Monin equations for compressible Hall magnetohydrodynamics, Astrophys. J. 917, 101 (2021).

[26] P. Hellinger, A. Verdini, S. Landi, E. Papini, L. Franci, and L. Matteini, Scale dependence and cross-scale transfer of kinetic energy in compressible hydrodynamic turbulence at moderate Reynolds numbers, Phys. Rev. Fluids 6, 044607 (2021).
[27] P. Hellinger, V. Montagud-Camps, L. Franci, L. Matteini, E. Papini, A. Verdini, and S. Landi, Ion-scale transition of plasma turbulence: pressure-strain effect, Astrophys. J. **930**, 48 (2022).

[28] L. Sorriso-Valvo, R. Marino, V. Carbone, A. Noullez, F. Lepreti, P. Veltri, R. Bruno, B. Bavassano, and E. Pietropaolo, Observation of inertial energy cascade in interplanetary space plasma, Phys. Rev. Lett. **99**, 115001 (2007).

[29] N. Andrés, F. Sahraoui, S. Galtier, L. Z. Hadid, R. Ferrand, and S. Y. Huang, Energy cascade rate measured in a collisionless space plasma with MMS data and compressible Hall magnetohydrodynamic turbulence theory, Phys. Rev. Lett. **123**, 245101 (2019).

[30] L. Sorriso-Valvo, F. Catapano, A. Retinò, O. Le Contel, D. Perrone, O. W. Roberts, J. T. Coburn, V. Panebianco, F. Valentini, S. Perri, A. Greco, F. Malara, V. Carbone, P. Veltri, O. Pezzi, F. Fraternale, F. Di Mare, R. Marino, B. Giles, T. E. Moore et al., Turbulence-driven ion beams in the magnetospheric Kelvin-Helmholtz instability, Phys. Rev. Lett. **122**, 035102 (2019).

[31] L. Z. Hadid, F. Sahraoui, and S. Galtier, Energy cascade rate in compressible fast and slow solar wind turbulence, Astrophys. J. **838**, 9 (2017).

[32] K.T. Osman, W. H. Matthaeus, K. H. Kiyani, B. Hnat, and S. C. Chapman, Proton kinetic effects and turbulent energy cascade rate in the solar wind, Phys. Rev. Lett. **111**, 201101 (2013).

[33] V. Carbone, R. Marino, L. Sorriso-Valvo, A. Noullez, and R. Bruno, Scaling laws of turbulence and heating of fast solar wind: the role of density fluctuations, Phys. Rev. Lett. **103**, 061102 (2009).

[34] L. Z. Hadid, F. Sahraoui, S. Galtier, and S. Y. Huang, Compressible magnetohydrodynamic turbulence in the Earth's magnetosheath: Estimation of the energy cascade rate using *in situ* spacecraft data, Phys. Rev. Lett. **120**, 055102 (2018).

[35] R. Ferrand, S. Galtier, F. Sahraoui, and C. Federrath, Compressible turbulence in the interstellar medium: new insights from a high-resolution supersonic turbulence simulation, Astrophys. J. **904**, 160 (2020).

[36] P. D. Mininni and A. Pouquet, Finite dissipation and intermittency in magnetohydrodynamics, Phys. Rev. E **80**, 025401(R) (2009).

[37] A. G. Kritsuk, R. Wagner, and M. L. Norman, Energy cascade and scaling in supersonic isothermal turbulence, J. Fluid Mech. **729**, R1 (2013).

[38] R. Ferrand, F. Sahraoui, D. Laveder, T. Passot, P. L. Sulem, and S. Galtier, Fluid energy cascade rate and kinetic damping: new insight from 3D Landau-fluid simulations, Astrophys. J. **923**, 122 (2021).

[39] B. Jiang, C. Li, Y. Yang, K. Zhou, W. H. Matthaeus, and M. Wan, Energy transfer and third-order law in forced anisotropic magnetohydrodynamic turbulence with hyper-viscosity, J. Fluid Mech. **974**, A20 (2023).

[40] P. Simon and F. Sahraoui, Exact law for compressible pressure-anisotropic magnetohydrodynamic turbulence: Toward linking energy cascade and instabilities, Phys. Rev. E **105**, 055111 (2022).

[41] S. D. Bale, J. C. Kasper, G. G. Howes, E. Quataert, C. Salem, and D. Sundkvist, Magnetic fluctuation power near proton temperature anisotropy instability thresholds in the solar wind, Phys. Rev. Lett. **103**, 211101 (2009).

[42] P. Hellinger, S. Landi, L. Matteini, A. Verdini, and L. Franci, Mirror instability in the turbulent solar wind, Astrophys. J. **838**, 158 (2017).

[43] G. F. Chew, M. Goldberger, and F. E. Low, The Boltzmann equation and the one-fluid hydromagnetic equations in the absence of particle collisions, Proc. R. Soc. Lond. A Math. Phys. Sci. **236**, 112 (1956).

[44] P. Hunana, A. Tenerani, G. P. Zank, E. Khomenko, M. L. Goldstein, G. M. Webb, P. S. Cally, M. Collados, M. Velli, and L. Adhikari, An introductory guide to fluid models with anisotropic temperatures. Part 1. CGL description and collisionless fluid hierarchy, J. Plasma Phys. **85**, 205850602 (2019).

[45] P. Hunana and G. P. Zank, On the parallel and oblique firehose instability in fluid models, Astrophys. J. **839**, 13 (2017).

[46] F. Sahraoui, G. Belmont, and L. Rezeau, Hamiltonian canonical formulation of Hall-magnetohydrodynamics: toward an application to weak turbulence theory, Phys. Plasmas **10**, 1325 (2003).

[47] C. Lamriben, P.-P. Cortet, and F. Moisy, Direct measurements of anisotropic energy transfers in a rotating turbulence experiment, Phys. Rev. Lett. **107**, 024503 (2011).

[48] D. Manzini, F. Sahraoui, F. Califano, and R. Ferrand, Local cascade and dissipation in incompressible Hall magnetohydrodynamic turbulence: The coarse-graining approach, Phys. Rev. E **106**, 035202 (2022).

[49] J. Cho and A. Lazarian, Simulations of electron magnetohydrodynamic turbulence, Astrophys. J. **701**, 236 (2009).