Spectral Transfer and Kármán–Howarth–Monin Equations for Compressible Hall Magnetohydrodynamics

Turbulence is ubiquitous in astrophysical plasma environments (Matthaeus & Velli 2011). The solar wind constitutes a natural laboratory for studying turbulence in weakly collisional plasmas (Bruno & Carbone 2013). As it expands from the solar corona, it exhibits a strong nonadiabatic behavior, and to sustain its thermal ion energetic properties, it needs to be heated (Vasquez et al. 2007; Cranmer et al. 2009; Hellinger et al. 2011, 2013). The situation is less clear for electrons that carry a strong heat flux (Štverák et al. 2015). The solar wind is a strongly turbulent flow with large-amplitude fluctuations of the magnetic field, particle velocity field, and other quantities as well. These turbulent fluctuations are the usual suspect for the observed particle heating.

Besides phenomenological approaches, the Kármán–Howarth–Monin (KHM) equation (de Kármán & Howarth 1938; Kolmogorov 1941; Monin & Yaglom 1975; Frisch 1995) represents a way to determine the cascade/dissipation rate in turbulence. This equation was originally derived for incompressible (constant-density) hydrodynamic (HD) turbulence and further extended to incompressible magnetohydrodynamic (MHD) turbulence (Politano & Pouquet 1998) and to incompressible Hall MHD turbulence (Galtier 2008; Hellinger et al. 2018; Ferrand et al. 2019). Starting from the KHM equations and under further assumptions (e.g., isotropy), the so-called exact (scaling) laws can be derived to connect the third-order structure functions with the energy cascade/dissipation rate. In magnetized plasmas, the assumption of isotropy is, however, questionable owing to the anisotropy introduced by the ambient magnetic field (Shebalin et al. 1983; Oughton et al. 1994). Once relaxed, the cascade rate is given by the divergence of third-order structure functions.

In situ spacecraft observations of turbulence are based on 1D time series, and it is not fully clear to what extent these characterize the inherently 3D turbulent fluctuations and how well the KHM equation can be used to estimate the cascade rate (Podesta et al. 2009; Smith et al. 2018). First applications of the KHM equation to in situ observations showed linear scaling of the third-order structure functions (Sorriso-Valvo et al. 2007; Marino et al. 2008), and isotropy was assumed to estimate the corresponding cascade rate. To account for the expected anisotropy, a hybrid model that combines a 2D perpendicular (with respect to the ambient magnetic field) cascade with a parallel 1D one was developed (MacBride et al. 2008; Stawarz et al. 2009). The turbulent energy cascade and its anisotropy can be partly constrained by multispacecraft observations (Osman et al. 2011), but observational works need to be complemented by numerical simulations. Verdini et al. (2015) analyzed the effect of using different turbulence models for estimating the cascade rates using direct MHD simulations, showing that the isotropic approximation may lead to large errors in the estimation of the cascade rate. The hybrid model (MacBride et al. 2008; Stawarz et al. 2009) gives, on the other hand, relatively good energy cascade estimates. In this respect, Franci et al. (2020) compared in situ observations of turbulence with simulation results based on observation-driven plasma parameters, finding good quantitative agreement between spectral properties of the observed and simulated turbulent fluctuations, but also a good agreement between the observed and simulated third-order structure functions and the resulting cascade rates.

To date, there are not many studies estimating the energy cascade rate in the solar wind, and estimates based on the KHM equation are typically sufficient to explain the observed temperature radial profiles (Coburn et al. 2015; Smith & Vasquez 2021); however, estimated cascade rates exhibit a large variability and may even reach negative values in fast solar wind streams with a large cross-helicity (Smith et al. 2009). Most of the work is done at 1 au (or farther away from the Sun. Recently, Bandyopadhyay et al. (2020a) used data from the Parker Solar Probe at 0.17 and 0.25 au and showed (assuming isotropy) that the estimated cascade rate is consistent with the radial profile of the proton temperature. However, a systematic study of the radial dependence of the cascade rate is missing. On the other hand, it is not clear whether solar wind turbulence is strong enough to keep heating the plasma. While the numerical study of Montagud-Camps et al. (2018) shows that turbulence is able to sustain a radial profile of the temperature similar to the observed one in the inner heliosphere (at least within the MHD approximation with an isotropic temperature), phenomenological transport models of solar wind turbulence (Zhou & Matthaeus 1990; Oughton et al. 2011) indicate that, eventually at larger radial distances from the Sun, turbulent fluctuations would be exhausted and need to be, in turn, sustained by some energy injection.

The KHM equations can also be used to understand at which scales the dissipation becomes important. Solar wind turbulence exhibits a transition in the form of a spectral break, at ion scales (Chen et al. 2014), and a similar transition is also observed in direct numerical simulations (e.g., Franci et al. 2016; Papini et al. 2019a, and references therein). Analyses of in situ observations and numerical simulations based on the incompressible Hall MHD KHM equation suggest that this transition is a combination of the onset of Hall physics and dissipation (Hellinger et al. 2018; Bandyopadhyay et al. 2020b). These results have some limitations: the simulations are 2D, and the observations are based on 1D time series and assume isotropy.

The solar wind exhibits weak density fluctuations, δ n/n ∼ 0.1, so that the incompressible approximation (or the nearly incompressible one; cf. Zank et al. 2017) is likely applicable. On the other hand, closer to the Sun the compressibility is expected to be larger (cf. Krupar et al. 2020). Furthermore, at ion characteristic scales the level of compressibility (density fluctuations) increases (Chen et al. 2013; Pitňa et al. 2019), so that the incompressible approximation may not hold. Also in high cross-helicity flows, the incompressible nonlinearity may be reduced and the compressible coupling may become important. Banerjee et al. (2016) and Hadid et al. (2017) show that the compressible KHM equation gives larger cascade rates. Andrés et al. (2019) show that the compressibility may be important at sub-ion scales, where the Hall physics becomes important. However, these results (Banerjee et al. 2016; Hadid et al. 2017; Andrés et al. 2019, and references therein) have a major limitation: they include the internal energy in the isothermal closure (consequently, they estimate the cascade rate of the total energy including the redistribution of the internal energy).

In order to avoid issues in the calculation of the divergence of the third-order structure function, Banerjee & Galtier (2017) proposed an alternative form of the KHM equation (originally formulated for the incompressible Hall MHD, and extended to the compressible flow by Banerjee & Kritsuk 2018) where the cascade rate is expressed using second-order structure functions. In this paper we analyze three different approaches that may be used to estimate the cascade rates of the kinetic + magnetic energy in compressible Hall MHD turbulence. We derive a new form of the compressible (standard) KHM equation (Hellinger et al. 2018; Ferrand et al. 2019), and we also generalize the alternative formulation of Banerjee & Galtier (2017) to the compressible Hall MHD. We compare these KHM approaches with an isotropic spectral transfer equation (Hellinger et al. 2021). We use these methods to analyze results of a 2D weakly compressible Hall MHD simulation. We also look at reduced 1D results in the context of in situ observations. The paper is organized as follows. In Section 2, we derive the standard (Section 2.1) and the alternative KHM (Section 2.2) equations. In Section 3, we analyze the results of a 2D Hall MHD simulation using the KHM equations. In Section 4, we present and apply the ST analysis to the Hall MHD simulation. Results of the ST and KHM analyses are then compared. In Section 5, we analyze the energy cascade rates from 1D reduced forms of the two KHM equation and of the ST one. In Section 6 we discuss the obtained results.

We investigate a system governed by the following Hall MHD equations for the plasma density ρ, the plasma mean velocity u , and the magnetic field B :

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+({\boldsymbol{u}}\cdot {\boldsymbol{\nabla }})\rho =-\rho {\boldsymbol{\nabla }}\cdot {\boldsymbol{u}},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\rho \displaystyle \frac{\partial {\boldsymbol{u}}}{\partial t}+\rho ({\boldsymbol{u}}\cdot {\boldsymbol{\nabla }}){\boldsymbol{u}}=({\boldsymbol{\nabla }}\times {\boldsymbol{B}})\times {\boldsymbol{B}}-{\boldsymbol{\nabla }}p+{\boldsymbol{\nabla }}\cdot {\boldsymbol{\tau }},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}={\boldsymbol{\nabla }}\times \left[({\boldsymbol{u}}-{\boldsymbol{j}})\times {\boldsymbol{B}}\right]+\eta {\rm{\Delta }}{\boldsymbol{B}},\end{eqnarray} \tag{ 3 }$$

where p is the plasma pressure, τ is the viscous stress tensor (given by ${\tau }_{\mathrm{ij}}=\mu \left(\partial {u}_{i}/\partial {x}_{j}+\partial {u}_{j}/\partial {x}_{i}-2/3{\delta }_{\mathrm{ij}}\partial {u}_{k}/\partial {x}_{k}\right)$, where the dynamic viscosity μ is assumed to be constant), η is the electric resistivity, j is the electric current density in velocity units, and j = J /ρ_c = u − u _e (ρ_c and u _e being the charge density and the electron velocity, respectively). Here we assume SI units except for the magnetic permeability μ₀, which is set to 1 (SI results can be obtained by the rescaling ${\boldsymbol{B}}\to {\boldsymbol{B}}{\mu }_{0}^{-1/2}$).

For the formulation of the KHM equation in terms of structure functions, we use the density-weighted velocity field w = ρ^1/2 u (Kida & Orszag 1990; Praturi & Girimaji 2019; Schmidt & Grete 2019) to take into account a variable density. In order to represent the kinetic and magnetic energy, we use the following structure functions:

$$\begin{eqnarray*}&&{{ \mathcal S }}_{w}=\left\langle | \delta {\boldsymbol{w}}{| }^{2}\right\rangle \,\mathrm{and}\ \ {{ \mathcal S }}_{B}=\left\langle | \delta {\boldsymbol{B}}{| }^{2}\right\rangle ,\end{eqnarray*}$$

respectively. Here $\delta {\boldsymbol{w}}={\boldsymbol{w}}({{\boldsymbol{x}}}^{{\prime} })-{\boldsymbol{w}}({\boldsymbol{x}})$, where ${{\boldsymbol{x}}}^{{\prime} }={\boldsymbol{x}}+{\boldsymbol{l}}$, and 〈•〉 denotes spatial averaging (over x ); the same definition holds for δ B and other quantities. Henceforth, we denote the value of any quantity a at x and ${{\boldsymbol{x}}}^{{\prime} }$ as a = a( x ) and ${a}^{{\prime} }=a({{\boldsymbol{x}}}^{{\prime} })$, respectively.

2.1. Standard KHM

For the total structure function, ${ \mathcal S }={{ \mathcal S }}_{w}+{{ \mathcal S }}_{B}$, one can get a dynamic KHM equation, taking Equations (2)–(3) at two positions, ${{\boldsymbol{x}}}^{{\prime} }$ and x . Subtracting those, one gets equations for ∂δ w /∂t and ∂δ B /∂t. Assuming a statistically homogeneous system, one can get then the following form of the KHM equation (for the detailed derivation see Appendix):

$$\begin{eqnarray}\,\begin{array}{rcl}\displaystyle \frac{\partial { \mathcal S }}{\partial t}+{{\boldsymbol{\nabla }}}_{{\boldsymbol{l}}}\cdot \left({ \mathcal Y }+{ \mathcal H }\right)+{ \mathcal R } & = & { \mathcal C }+2\left\langle \delta \theta \delta p\right\rangle -2\left\langle \delta {\boldsymbol{\Sigma }}:\delta {\boldsymbol{\tau }}\right\rangle \\ & & -\,4{Q}_{\eta }+2\eta {{\rm{\Delta }}}_{{\boldsymbol{l}}}{{ \mathcal S }}_{B},\end{array}\,\end{eqnarray} \tag{ 4 }$$

where θ = ∇ · u is the dilatation field, Σ = ∇ u is the stress tensor, and Δ _l denotes the Laplace operator (with respect to the separation distance l ). Equation (4) combines the second-order structure functions with the third-order ones related to the energy transfer/cascade rate

$$\begin{eqnarray*}\begin{array}{rcl}{ \mathcal Y } & = & \left\langle \delta {\boldsymbol{u}}\left(| \delta {\boldsymbol{w}}{| }^{2}+| \delta {\boldsymbol{B}}{| }^{2}\right)-2\delta {\boldsymbol{B}}\left(\delta {\boldsymbol{u}}\cdot \delta {\boldsymbol{B}}\right)\right\rangle \\ { \mathcal H } & = & \left\langle \delta {\boldsymbol{B}}\left(\delta {\boldsymbol{j}}\cdot \delta {\boldsymbol{B}}\right)-\displaystyle \frac{1}{2}\delta {\boldsymbol{j}}| \delta {\boldsymbol{B}}{| }^{2}\right\rangle .\end{array}\end{eqnarray*}$$

Equation (4) also contains the explicitly compressible contribution to the cascade rate,

$$\begin{eqnarray*}&&{ \mathcal R }=\left\langle \delta {\boldsymbol{w}}\cdot \left({\theta }^{{\prime} }{\boldsymbol{w}}-\theta {{\boldsymbol{w}}}^{{\prime} }\right)\right\rangle ,\end{eqnarray*}$$

of hydrodynamic origin (cf. Hellinger et al. 2021). The dissipation part includes the resistive incompressible-like terms ${Q}_{\eta }-\eta {{\rm{\Delta }}}_{{\boldsymbol{l}}}{{ \mathcal S }}_{B}/2$ and the generally compressible viscous term $\left\langle \delta {\boldsymbol{\tau }}:\delta {\boldsymbol{\Sigma }}\right\rangle /2$, where ${Q}_{\eta }=\eta \left\langle {\boldsymbol{\nabla }}{\boldsymbol{B}}:{\boldsymbol{\nabla }}{\boldsymbol{B}}\right\rangle =\eta \left\langle | {\boldsymbol{J}}{| }^{2}\right\rangle$ is the Joule dissipation rate per unit volume ( J = ∇ × B being the electric current density). The correction term ${ \mathcal C }$ can be given in the form

$$\begin{eqnarray*}\begin{array}{rcl}{ \mathcal C } & = & 2{C}_{\sqrt{\rho }}\left[{\boldsymbol{u}},{\boldsymbol{\nabla }}p\right]-2{C}_{\sqrt{\rho }}\left[{\boldsymbol{u}},{\boldsymbol{\nabla }}\cdot {\boldsymbol{\tau }}\right]\\ & & +2{C}_{\sqrt{\rho }}\left[{\boldsymbol{u}},{\boldsymbol{B}}\times {\boldsymbol{J}}\right]+{C}_{\rho }\left[{\boldsymbol{B}}\times {\boldsymbol{j}},{\boldsymbol{J}}\right],\end{array}\end{eqnarray*}$$

where

$$\begin{eqnarray*}&&{C}_{\rho }\left[{\boldsymbol{a}},{\boldsymbol{b}}\right]=\left\langle \left(\displaystyle \frac{{\rho }^{{\prime} }}{\rho }-1\right){{\boldsymbol{a}}}^{{\prime} }\cdot {\boldsymbol{b}}+\left(\displaystyle \frac{\rho }{{\rho }^{{\prime} }}-1\right){\boldsymbol{a}}\cdot {{\boldsymbol{b}}}^{{\prime} }\right\rangle .\end{eqnarray*}$$

Note that this term explicitly depends on the density variation and disappears in the constant-density approximation.

Equation (4) can be cast in the following form:

$$\begin{eqnarray}&&\displaystyle \frac{1}{4}\displaystyle \frac{{\rm{\partial }}{ \mathcal S }}{{\rm{\partial }}t}-{{ \mathcal K }}_{{\rm{H}}}-{{ \mathcal K }}_{{\rm{MHD}}}={\rm{\Psi }}-{ \mathcal D },\end{eqnarray} \tag{ 5 }$$

where we have combined some terms as

$$\begin{eqnarray}\begin{array}{c}\begin{array}{c}\,\begin{array}{rcl}{{ \mathcal K }}_{{\rm{MHD}}} & = & -\displaystyle \frac{1}{4}{\boldsymbol{\nabla }}\cdot { \mathcal Y }-\displaystyle \frac{1}{4}{ \mathcal R }+\displaystyle \frac{1}{2}{C}_{\sqrt{\rho }}\left[{\boldsymbol{u}},{\boldsymbol{B}}\times {\boldsymbol{J}}\right]\\ {{ \mathcal K }}_{{\rm{H}}} & = & -\displaystyle \frac{1}{4}{\boldsymbol{\nabla }}\cdot { \mathcal H }+\displaystyle \frac{1}{4}{C}_{\rho }\left[{\boldsymbol{B}}\times {\boldsymbol{j}},{\boldsymbol{J}}\right]\\ {\rm{\Psi }} & = & \displaystyle \frac{1}{2} \langle \delta p\delta \theta \rangle +\displaystyle \frac{1}{2}{C}_{\sqrt{\rho }}\left[{\boldsymbol{u}},{\boldsymbol{\nabla }}p\right]\\ { \mathcal D } & = & {Q}_{\eta }-\displaystyle \frac{\eta }{2}{\rm{\Delta }}{{ \mathcal S }}_{B}+\displaystyle \frac{1}{2} \langle \delta {\boldsymbol{\tau }}:\delta {\boldsymbol{\Sigma }} \rangle +\displaystyle \frac{1}{2}{C}_{\sqrt{\rho }}\left[{\boldsymbol{u}},{\boldsymbol{\nabla }}\cdot {\boldsymbol{\tau }}\right],\end{array}\,\end{array}\end{array}\end{eqnarray} \tag{ 6 }$$

and where we dropped the l index for ∇ and Δ. Here, ${{ \mathcal K }}_{\mathrm{MHD}}$ and ${{ \mathcal K }}_{{\rm{H}}}$ are the MHD and Hall cascade rates, respectively, Ψ represents the pressure dilatation effect, whereas ${ \mathcal D }$ accounts for the effects of dissipation and heating.

As noted by Hellinger et al. (2021), the KHM structure functions in hydrodynamic (HD) turbulence have a cumulative behavior that is complementary to the isotropic ST equation. We will see in Section 4.1 that this is also true in Hall MHD. Here we just note that the viscous term $\left\langle \delta {\boldsymbol{\tau }}:\delta {\boldsymbol{\Sigma }}\right\rangle$ is the generalization to the compressible case of the two incompressible dissipative terms $2\nu \left\langle {\boldsymbol{\nabla }}{\boldsymbol{u}}:{\boldsymbol{\nabla }}{\boldsymbol{u}}\right\rangle -\nu {\rm{\Delta }}\left\langle \delta {{\boldsymbol{u}}}^{2}\right\rangle$. At large separation distances, ∣ l ∣ → ∞ , where the correlations $\left\langle {\boldsymbol{\tau }}({{\boldsymbol{x}}}^{{\prime} }):{\boldsymbol{\Sigma }}({\boldsymbol{x}})\right\rangle \to 0$, the viscous term becomes twice the viscous heating rate Q_μ , $\left\langle \delta {\boldsymbol{\tau }}:\delta {\boldsymbol{\Sigma }}\right\rangle \to 2\left\langle {\boldsymbol{\tau }}:{\boldsymbol{\Sigma }}\right\rangle =2{Q}_{\mu }$, and the dissipation term reaches ${ \mathcal D }\to {Q}_{\mathrm{tot}}$, where we denote the total heating rate as Q_tot =`Q_μ + Q_η .

2.2. Alternative Formulation

The cross-scale energy transfer/cascade rates in Equation (4) contain terms with a divergence of third-order structure functions that are difficult to evaluate from 1D time series. Banerjee & Galtier (2017) presented an alternative formulation of the KHM equation for incompressible Hall MHD turbulence where these terms are replaced by second-order structure functions. This approach can be easily generalized to the compressible case, and the resulting equation can be given the following form:

$$\begin{eqnarray}\,\begin{array}{rcl}\displaystyle \frac{\partial { \mathcal S }}{\partial t}+\tilde{{ \mathcal Y }}+\tilde{{ \mathcal H }}+\tilde{{ \mathcal R }} & = & \tilde{{ \mathcal C }}+2\left\langle \delta \theta \delta p\right\rangle -2\left\langle \delta {\boldsymbol{\Sigma }}:\delta {\boldsymbol{\tau }}\right\rangle \\ & & -\,4{Q}_{\eta }+2\eta {\rm{\Delta }}{{ \mathcal S }}_{B},\end{array}\,\end{eqnarray} \tag{ 7 }$$

where

$$\begin{eqnarray}\,\begin{array}{rcl}\tilde{{ \mathcal Y }} & = & -\,2\left\langle \delta {\boldsymbol{w}}\cdot \delta \left({\boldsymbol{u}}\times {\boldsymbol{\varpi }}\right)\right\rangle -2\left\langle \delta {\boldsymbol{w}}\cdot \delta \left(\displaystyle \frac{{\boldsymbol{J}}\times {\boldsymbol{B}}}{\sqrt{\rho }}\right)\right\rangle \\ & & -\,2\left\langle \delta {\boldsymbol{J}}\times \delta \left({\boldsymbol{u}}\times {\boldsymbol{B}}\right)\right\rangle \end{array}\,\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\tilde{{ \mathcal H }}=2\left\langle \delta {\boldsymbol{J}}\times \delta \left({\boldsymbol{j}}\times {\boldsymbol{B}}\right)\right\rangle \end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&\tilde{{ \mathcal R }}=\left\langle \delta {\boldsymbol{w}}\cdot \delta \left({\boldsymbol{w}}\theta \right)\right\rangle +2\left\langle \delta {\boldsymbol{w}}\cdot \delta \left[({\boldsymbol{\nabla }}{\boldsymbol{w}})\cdot {\boldsymbol{u}}\right]\right\rangle \end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&\tilde{{ \mathcal C }}=2{C}_{\sqrt{\rho }}\left[{\boldsymbol{u}},{\boldsymbol{\nabla }}p\right]-2{C}_{\sqrt{\rho }}\left[{\boldsymbol{u}},{\boldsymbol{\nabla }}\cdot {\boldsymbol{\tau }}\right]\end{eqnarray} \tag{ 11 }$$

and ϖ = ∇ × w . The other terms are defined in the previous section.

From Equation (7) we get the MHD and Hall cascade rates as

$$\begin{eqnarray}&&{{ \mathcal K }}_{{\rm{M}}{\rm{H}}{\rm{D}}}=-\displaystyle \frac{1}{4}\tilde{{ \mathcal Y }}-\displaystyle \frac{1}{4}\tilde{{ \mathcal R }},\,\,{{ \mathcal K }}_{{\rm{H}}}=-\displaystyle \frac{1}{4}\tilde{{ \mathcal H }}.\end{eqnarray} \tag{ 12 }$$

Note that only the first term in $\tilde{{ \mathcal R }}$ disappears in the incompressible, constant-density limit. The second term, in this limit, is proportional to 〈δ u · δ(∇∣ u ∣²)〉 (cf. Banerjee & Galtier 2017) and is generally nonnegligible.

Now we test the predictions of the compressible KHM equations using a Hall MHD simulation. The Hall MHD model consists of the nonlinear, compressible viscous-resistive MHD equations, modified only by the presence of the Hall term in the induction equation. The system is described by Equations (1)–(3), complemented by the equation for the plasma temperature

$$\begin{eqnarray}&&\displaystyle \frac{\partial T}{\partial t}+({\boldsymbol{u}}\cdot {\boldsymbol{\nabla }})T=(\gamma -1)\left(-T{\boldsymbol{\nabla }}\cdot {\boldsymbol{u}}+\displaystyle \frac{\eta | {\boldsymbol{J}}{| }^{2}}{\rho }+\displaystyle \frac{{\boldsymbol{\tau }}:{\boldsymbol{\Sigma }}}{\rho }\right),\end{eqnarray} \tag{ 13 }$$

where γ = 5/3. These equations are solved using the pseudo-spectral approach and a third-order Runge–Kutta scheme (Papini et al. 2019a, 2019b). We consider a 2D (x, y) periodic domain and use Fourier decomposition to calculate the spatial derivatives. In the Fourier space we also filter according to the 2/3 Orszag rule (Orszag 1971), to avoid aliasing of the quadratic nonlinear terms. Aliasing of the cubic terms is mitigated by the presence of a finite dissipation (Ghosh et al. 1993). We consider a 2D box of size 256 d_i × 256 d_i and a grid resolution of Δx = Δy = d_i /8, corresponding to 2048² points, where d_i is the ion inertial length. We set a constant background magnetic field B ₀ along the z (out-of-plane) direction. The initial state is populated by large-amplitude Alfvénic fluctuations in the xy-plane with $k={({k}_{x}^{2}+{k}_{y}^{2})}^{1/2}\leqslant {k}^{\mathrm{inj}}$, where the injection scale is ${k}^{\mathrm{inj}}=0.2{d}_{i}^{-1}$. The relative rms amplitude of these fluctuations is set to 0.17 in units of B₀ = ∣ B ₀∣. Initially, the system has no cross-helicity, 〈 u · B 〉 ≃ 0 (here 〈•〉 denotes spatial averaging over the simulation box), and the density and the temperature are uniform with the plasma beta β = 0.5. No forcing is present, so the simulation resolves the evolution of freely decaying turbulence, and we assume the viscosity and resistivity μ = η = 10⁻³.

In the simulation the total energy E_tot = E_kin + E_int + E_mag is well conserved. Here E_kin = 〈ρ u²〉/2 is the kinetic energy, E_int = 〈ρ T〉/(γ − 1) is the internal one, and E_mag = 〈B²〉/2. Figure 1(a) displays the evolution of the relative changes in these energies, ΔE(t) = [E(t) − E(0)]/E_tot(0), where E =E_{kin,tot,mag,int}. The relative change of the total energy is negligible, ΔE_tot(t = 8) ∼ 3 × 10⁻⁵. Figure 1(b) shows the evolution of the rms values of the out-of-plane components of the current J_z and the vorticity ω_z . 〈∣J_z ∣²〉 reaches a maximum at $t\simeq 255{{\rm{\Omega }}}_{i}^{-1}$, which is a signature of a fully developed turbulent cascade (Mininni & Pouquet 2009). The z component of the vorticity exhibits a similar evolution. Figure 1(c) shows the maximum value (over the simulation box) of J_z . The sharp increase of $\max | {J}_{z}|$ for 40 ≲ tΩ_i ≲ 60 indicates the formation of thin current sheets, and the saturation at tΩ_i ∼ 60 is due to the onset of magnetic reconnection (Franci et al. 2017; Papini et al. 2019a). Figure 1(d) displays the rms value of the density, 〈(ρ − 〈ρ〉)²〉^1/2, which steadily increases until at $t\gtrsim 210\,{{\rm{\Omega }}}_{i}^{-1}$ there is an indication of a saturation. Figure 1(e) shows the dissipation and the pressure dilatation rates. The dissipation rate increases with time and becomes quasi-stationary at later times. The pressure dilatation rate becomes initially negative, but at later times, $t\gtrsim 210\,{{\rm{\Omega }}}_{i}^{-1}$, it oscillates around zero.

Zoom In Zoom Out Reset image size

To check the prediction of Equation (5) in the simulation, we define the validity test as

$$\begin{eqnarray}&&{ \mathcal O }=-\displaystyle \frac{1}{4}\displaystyle \frac{{\rm{\partial }}{ \mathcal S }}{{\rm{\partial }}t}+{{ \mathcal K }}_{{\rm{H}}}+{{ \mathcal K }}_{{\rm{MHD}}}+{\rm{\Psi }}-{ \mathcal D }.\end{eqnarray} \tag{ 14 }$$

The top panel of Figure 2 shows this validity check ${ \mathcal O }$ in the 2D simulation at $t=255\,{{\rm{\Omega }}}_{i}^{-1}$, as well as the contributing terms. The compressible KHM equation is very well conserved: $| { \mathcal O }| /Q\lt 0.3 \%$, where Q is the heating/dissipation rate at $t=255\,{{\rm{\Omega }}}_{i}^{-1}$, $Q={Q}_{\mathrm{tot}}(255\,{{\rm{\Omega }}}_{i}^{-1});$ henceforth we use Q to normalize all the relevant quantities. At large separation distances, ${\rm{\partial }}{ \mathcal S }/{\rm{\partial }}t/4-{\Psi }\sim -{ \mathcal D }\sim -Q$, which represents the energy conservation. The pressure dilatation term Ψ is small. At intermediate separations the MHD cascade dominates, ${{ \mathcal K }}_{\mathrm{MHD}}\,\sim 0.85\,Q$. This term is mostly compensated by the dissipation term ${ \mathcal D }$. At small separation distances the Hall cascade sets in and becomes comparable to the MHD one. On these scales the viscous and resistive dissipation also starts to act.

Zoom In Zoom Out Reset image size

The pressure dilatation effect is weak but not negligible. Since the energy exchanges due to this effect oscillate around zero, it is interesting to look at the average behavior of the KHM equation over one period of these oscillations (see Figure 1; cf. Hellinger et al. 2021). The bottom panel of Figure 2 displays with solid lines the different terms contributing to the KHM equation averaged over such a period (240 ≤ tΩ_i ≤ 264; henceforth 〈•〉_t denotes the average over this time); the dotted lines show the corresponding minimum and maximum values. The pressure dilatation term varies within a few percent of Q over this period of time, and, on average, the effect of the pressure dilatation effect is negligible. The terms $-\partial { \mathcal S }/\partial t/4$ and ${{ \mathcal K }}_{\mathrm{MHD}}$ exhibit similar temporal variability (they have similar ranges of values during this period; see the dotted lines in the bottom panel of Figure 2), whereas ${{ \mathcal K }}_{{\rm{H}}}$ and ${ \mathcal D }$ are about constant in time.

3.1. Evolution

The compressible KHM equation is valid during the whole simulation, and the homogeneity condition, $\langle ({{\boldsymbol{\nabla }}}_{{{\boldsymbol{x}}}^{{\prime} }}+{{\boldsymbol{\nabla }}}_{{\boldsymbol{x}}})\bullet \rangle \,=0$, is automatically satisfied for the periodic boundary conditions. Figure 3 displays the validity test ${ \mathcal O }$, as well as the different contributing terms (normalized to Q), as functions of time and l. Figure 3(a) shows that the compressible KHM equation is well conserved, $\max | { \mathcal O }| /Q\lesssim 1 \%$ for $t\gt 20\,{{\rm{\Omega }}}_{i}^{-1}$. During early times, the maximum relative error in $| { \mathcal O }| /Q$ is about 2%. At later times, the $| { \mathcal O }| /Q$ is below 1% and is fluctuating with time and roughly constant at all scales at a given time. This is caused by numerical errors; in particular, the time derivative $\partial { \mathcal S }/\partial t$ is estimated by the finite difference with ${\rm{\Delta }}t={{\rm{\Omega }}}_{i}^{-1}$. This is likely not sufficient for the transition from the initial superposition of large-scale Alfvénic fluctuations. The interpretation is supported by the fact that ${ \mathcal O }$ increases when we increase Δt.

Zoom In Zoom Out Reset image size

Figures 3(b) and (c) show that initially $\partial { \mathcal S }/\partial t$ and ${{ \mathcal K }}_{\mathrm{MHD}}$ compensate each other, $\partial { \mathcal S }/\partial t+4{{ \mathcal K }}_{\mathrm{MHD}}\sim 0$. $\partial { \mathcal S }/\partial t$ evolves and, at later times, dominates at large separation distances, indicating a decay of the kinetic + magnetic energy at large scales. ${{ \mathcal K }}_{\mathrm{MHD}}$ becomes dominant at intermediate separation distances. Figure 3(d) shows that the Hall term becomes important at 40 ≲ tΩ_i ≲ 60 for l ≲ 3 d_i ; this is about the time when thin current sheets form and start to reconnect (Papini et al. 2019a). Interestingly, there is an indication of a negative Hall cascade rate, starting earlier 40 ≲ tΩ_i for 3 d_i ≲ l ≲ 20 d_i , indicating an energy transfer from small to large scales, which may be part of the onset of turbulence (cf. Franci et al. 2017). The Hall term settles to the asymptotic value after tΩ_i ≃ 150, when the energy at large scales has had time to participate in the cascade.

Figure 3(e) demonstrates that the dissipation gradually develops as the (MHD and Hall) cascade evolves and brings energy to small scales. Figure 3(f) displays the pressure dilatation term Ψ. After the initial transition, this term is important on relatively large scales for a long time; its negative value indicates a compressible plasma heating. Then, it becomes weakly oscillating around zero for $t\gtrsim 210\,{{\rm{\Omega }}}_{i}^{-1}$.

Note that the alternative KHM equation, Equation (7), gives, unsurprisingly, almost the same results as Equation (4). The cascade rates ${{ \mathcal K }}_{\mathrm{MHD}}$ and ${{ \mathcal K }}_{{\rm{H}}}$ differ within 10⁻³ Q between the two methods. On the other hand, for the weakly compressible simulation, the incompressible approximation (Hellinger et al. 2018; Ferrand et al. 2019) exhibits an error that is mostly related to the neglected pressure dilatation term. However, the incompressible MHD and Hall cascade rates are close to their compressible counterparts (not shown).

Another way to analyze the scale dependence of turbulence and its processes is the spectral (Fourier) decomposition. To characterize the kinetic energy in the compressible flow, one can use the density-weighted velocity field w as in the KHM approach. The evolution for the energy in a given Fourier mode

$$\begin{eqnarray}&&\hat{{\boldsymbol{w}}}({\boldsymbol{k}})=\displaystyle \sum _{{\boldsymbol{x}}}{\boldsymbol{w}}({\boldsymbol{x}})\exp (i{\boldsymbol{k}}\cdot {\boldsymbol{x}}),\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&\hat{{\boldsymbol{B}}}({\boldsymbol{k}})=\displaystyle \sum _{{\boldsymbol{x}}}{\boldsymbol{B}}({\boldsymbol{x}})\exp (i{\boldsymbol{k}}\cdot {\boldsymbol{x}})\end{eqnarray} \tag{ 16 }$$

follows from Equations (1)–(3) (cf. Mininni et al. 2007; Grete et al. 2017),

$$\begin{eqnarray}\begin{array}{c}\,\begin{array}{rcl}\displaystyle \frac{1}{2}\displaystyle \frac{{\rm{\partial }}\left(| \hat{{\boldsymbol{w}}}{| }^{2}+| \hat{{\boldsymbol{B}}}{| }^{2}\right)}{{\rm{\partial }}t} & = & -{T}_{{\rm{MHD}}}-{T}_{{\rm{H}}}-{\mathfrak{R}}{\hat{{\boldsymbol{w}}}}^{\ast }\cdot \hat{\displaystyle \frac{{\boldsymbol{\nabla }}p}{\sqrt{\rho }}}\\ & & +\,{\mathfrak{R}}{\hat{{\boldsymbol{w}}}}^{\ast }\cdot \hat{\displaystyle \frac{{\boldsymbol{\nabla }}\cdot {\boldsymbol{\tau }}}{\sqrt{\rho }}}-\eta {k}^{2}| \hat{{\boldsymbol{B}}}{| }^{2},\end{array}\,\end{array}\end{eqnarray} \tag{ 17 }$$

where the MHD and Hall transfer terms are

$$\begin{eqnarray*}\begin{array}{c}\begin{array}{rcl}{T}_{{\rm{MHD}}} & = & {\mathfrak{R}}\left[{\hat{{\boldsymbol{w}}}}^{\ast }\cdot \hat{({\boldsymbol{u}}\cdot {\boldsymbol{\nabla }}){\boldsymbol{w}}}+\displaystyle \frac{1}{2}{\hat{{\boldsymbol{w}}}}^{\ast }\cdot \hat{{\boldsymbol{w}}\theta }\right.\\ & & \left.-\,{\hat{{\boldsymbol{w}}}}^{\ast }\cdot \hat{\displaystyle \frac{({\boldsymbol{\nabla }}\times {\boldsymbol{B}})\times {\boldsymbol{B}}}{\sqrt{\rho }}}-{\hat{{\boldsymbol{B}}}}^{\ast }\cdot \hat{{\boldsymbol{\nabla }}\times \left({\boldsymbol{u}}\times {\boldsymbol{B}}\right)}\right],\\ {T}_{{\rm{H}}} & = & {\mathfrak{R}}\left[{\hat{{\boldsymbol{B}}}}^{\ast }\cdot \hat{{\boldsymbol{\nabla }}\times \left({\boldsymbol{j}}\times {\boldsymbol{B}}\right)}\right],\end{array}\end{array}\end{eqnarray*}$$

respectively; here the wide hat denotes the Fourier transform, the star denotes the complex conjugate, and ${\mathfrak{R}}$ denotes the real part.

In the inertial range, one expects that the time derivative terms in Equation (17) are zero, and the same is expected for the dissipation and pressure dilatation terms. The transfer term T_MHD + T_H is also expected to be zero, as all the energy that arrives from larger scales proceeds to smaller scales. To characterize (isotropic) turbulence, we use a low-pass-filtered kinetic + magnetic energy (i.e., the energy in modes with wavevector magnitudes smaller than or equal to k; cf. Hellinger et al. 2021)

$$\begin{eqnarray}&&{E}_{k}=\displaystyle \frac{1}{2}\sum _{| {{\boldsymbol{k}}}^{{\rm{{\prime} }}}| \leqslant k}\left(| \hat{{\boldsymbol{w}}}{| }^{2}+| \hat{{\boldsymbol{B}}}{| }^{2}\right).\end{eqnarray} \tag{ 18 }$$

For E_k one gets the following dynamic equation:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\partial }}{E}_{k}}{{\rm{\partial }}t}+{S}_{{\rm{MHD}}k}+{S}_{{\rm{H}}k}={{\rm{\Psi }}}_{k}-{D}_{k},\end{eqnarray} \tag{ 19 }$$

where

$$\begin{eqnarray}\,\begin{array}{c}\begin{array}{c}\,\begin{array}{rcl}{S}_{{\rm{MHD}}k} & = & \displaystyle \sum _{\left|{{\boldsymbol{k}}}^{{\rm{{\prime} }}}| \leqslant k\right.}{T}_{{\rm{MHD}}}({{\boldsymbol{k}}}^{{\rm{{\prime} }}}),\\ {S}_{{\rm{H}}k} & = & \displaystyle \sum _{\left|{{\boldsymbol{k}}}^{{\rm{{\prime} }}}| \leqslant k\right.}{T}_{{\rm{H}}}({{\boldsymbol{k}}}^{{\rm{{\prime} }}}),\end{array}\,\end{array}\end{array}\,\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{k}=-\displaystyle \sum _{| {{\boldsymbol{k}}}^{{\rm{{\prime} }}}| \leqslant k}{\hat{{\boldsymbol{w}}}}^{\ast }\cdot \hat{{\rho }^{-1/2}{\boldsymbol{\nabla }}p},\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&{D}_{k}=\eta \displaystyle \sum _{| {{\boldsymbol{k}}}^{{\rm{{\prime} }}}| \leqslant k}{{k}^{{\rm{{\prime} }}}}^{2}| \hat{{\boldsymbol{B}}}{| }^{2}-\displaystyle \sum _{| {{\boldsymbol{k}}}^{{\rm{{\prime} }}}| \leqslant k}{\hat{{\boldsymbol{w}}}}^{\ast }\cdot \hat{{\rho }^{-1/2}{\boldsymbol{\nabla }}\cdot {\boldsymbol{\tau }}}.\end{eqnarray} \tag{ 22 }$$

Here S_MHDk and S_Hk represent the MHD and Hall energy transfer rates, respectively, Ψ_k describes the pressure dilatation effect, and D_k is the (viscous and resistive) dissipation rate for modes with wavevector magnitude smaller than or equal to k. When S_MHDk or S_Hk are constant, the corresponding energy transfer rate is constant and can be identified with the cascade rate. For large wavevectors, one gets the unfiltered values

$$\begin{eqnarray}&&{E}_{k}\to {E}_{\mathrm{kin}}+{E}_{\mathrm{mag}},\ \ {{\rm{\Psi }}}_{k}\to \langle p\theta \rangle ,\ \mathrm{and}\ \ {D}_{k}\to Q.\end{eqnarray} \tag{ 23 }$$

To check the spectral transfer of energy given by Equation (19), we define the validity test (similarly to Equation (14)) as

$$\begin{eqnarray}&&{O}_{k}=\displaystyle \frac{{\rm{\partial }}{E}_{k}}{{\rm{\partial }}t}+{S}_{{\rm{MHD}}k}+{S}_{{\rm{H}}k}-{{\rm{\Psi }}}_{k}+{D}_{k}.\end{eqnarray} \tag{ 24 }$$

Figure 4 (top panel) shows the validity test and the contributing terms as functions of k at $t=255\,{{\rm{\Omega }}}_{i}^{-1}$. Equation (19) is well satisfied, ∣O_k ∣/Q < 0.3%. At large scales, small positive values of ∂E_k /∂t compensated by S_MHDk indicate an inverse cascade/energy transfer from small to large scales. At medium scales, the MHD cascade term dominates, reaching a maximum value of about 0.85 Q; the MHD cascade term is compensated mostly by ∂E_k /∂t. At small scales, the Hall cascade sets in, but, at the same time, the dissipation gets important and a weak pressure dilatation effect appears. For large k one recovers the energy conservation

$$\begin{eqnarray}&&\displaystyle \frac{\partial ({E}_{\mathrm{kin}}+{E}_{\mathrm{mag}})}{\partial t}=\langle p\theta \rangle -Q.\end{eqnarray} \tag{ 25 }$$

Zoom In Zoom Out Reset image size

As in the KHM case, it is interesting to look at the properties of the ST equation averaged over one period of the pressure dilatation oscillation. Figure 4 (bottom panel) displays the terms contributing to the ST equation (Equation (19)) averaged over the time interval 240 ≤ tΩ_i ≤ 264, along with their minimum and maximum values. The pressure dilatation term Ψ_k exhibits a variation of a few percent of Q, and, on average, it is negligible. The decay ∂E_k /∂t and the MHD cascade S_MHDk also have a similar variability, whereas the Hall cascade S_Hk and the dissipation D_k tend to be about constant during the oscillation period.

The ST equation is valid during the whole simulation owing to the periodic boundary conditions as in the case of the KHM equation. Figure 5 displays the evolution of the ST equation in the simulation in a format similar to that of Figure 3.

Zoom In Zoom Out Reset image size

Figure 5 shows that the ST and KHM equations give similar results. Note, moreover, that a good quantitative agreement between the two approaches appears at later times when turbulence is well developed.

4.1. Comparison between KHM and ST Approaches

The KHM and ST equations give similar results. They are empirically related through the inverse proportionality between the wavevector k and the spatial separation distance l as ${kl}\simeq \sqrt{2}$. The top panel of Figure 6 shows the power spectra of the velocity (blue), magnetic field (red), and the total power spectrum as functions of k (at $t=255{{\rm{\Omega }}}_{i}^{-1}$) as a reference. The dashed line denotes a k^−5/3 power law for comparison.

Zoom In Zoom Out Reset image size

The bottom panel of Figure 6 displays a direct comparison between the ST and KHM contributing terms as functions of k (through the inverse proportionality ${kl}=\sqrt{2}$). Figure 6 shows that the MHD and Hall cascade rates in both approaches are comparable: $\langle {S}_{\mathrm{MHD}k}{\rangle }_{t}\simeq \langle {{ \mathcal K }}_{\mathrm{MHD}}{\rangle }_{t}$, $\langle {S}_{{\rm{H}}k}{ \rangle }_{t}\simeq \langle {{ \mathcal K }}_{{\rm{H}}}{ \rangle }_{t}$. Note that the MHD cascade term dominates at scales where the total power spectrum is close to a k^−5/3 power law. The average pressure dilatation effect is negligible. Similarly to the hydrodynamic case (Hellinger et al. 2021), the ST and KHM approaches are complementary:

$$\begin{eqnarray}&&\langle \partial {E}_{k}/\partial t{\rangle }_{t}+\langle \partial { \mathcal S }/\partial t{\rangle }_{t}/4\simeq -Q\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&\langle {D}_{k}{\rangle }_{t}+\langle { \mathcal D }{\rangle }_{t}\simeq Q.\end{eqnarray} \tag{ 27 }$$

∂E_k /∂t represents the rate of change of the kinetic energy at scales with wavevector magnitudes smaller than or equal to k, whereas $\partial { \mathcal S }/\partial t/4$ represents approximately the rate for wavevector magnitudes larger than k. Similarly, D_k is the dissipation rate on the scales ≤ k, whereas ${ \mathcal D }$ represents the dissipation rate on the scales > k.

In the present simulation, due to the reduced dimensionality, we cannot address the question of the spectral anisotropy (Verdini et al. 2015). We can, however, test what estimates of the cascade rate we can get from 1D cuts representing observed 1D time series. We take 1D cuts (in both the x- and y-direction) to calculate the ST predictions ${S}_{\mathrm{MHD}}^{(1{\rm{D}})}$ and ${S}_{{\rm{H}}}^{(1{\rm{D}})}$ and the KHM predictions ${{ \mathcal K }}_{\mathrm{MHD}}^{(1{\rm{D}})}$ and ${{ \mathcal K }}_{{\rm{H}}}^{(1{\rm{D}})}$.

For the reduced 1D ST equation we take the 1D Fourier transform in Equation (19) and retain all the spatial derivatives, and, similarly, we use all the spatial derivatives (in the real space) for the alternative KHM equation. On the other hand, for the standard KHM equation, Equation (4), only the derivative along the 1D direction in the separation space l is used to estimate ${\boldsymbol{\nabla }}\cdot ({ \mathcal Y }+{ \mathcal H })$. Results from these estimates are shown in Figure 7. This figure shows the MHD and Hall cascade rates obtained from the reduced 1D analyses as functions of k for the ST approach and k = 1/l for the KHM approaches. The calculation is done for the time $t=255\,{{\rm{\Omega }}}_{i}^{-1}$; see top panels of Figures 4 and 2. The reduced ST equation gives cascade rates relatively similar to those obtained from the whole simulation domain. On the other hand, the reduced 1D results based on the standard KHM equation, Equation (4), are quite different from the full results: the maximum value of the MHD cascade rate is notably about half the expected value. This can be partly amended by assuming isotropy, ${\boldsymbol{\nabla }}\cdot { \mathcal Y }\sim 2\partial {{ \mathcal Y }}_{l}/\partial l$. However, even after this correction, the agreement between the 1D ST and standard KHM method is poor. Interestingly, the reduced 1D results based on the alternative KHM equation, Equation (7), are in a good agreement with the full results, as well as with the results of the reduced 1D ST equation. This is likely related to the fact that for the alternative, divergence-less formulation we use higher-dimensionality derivatives that correspond to the divergence in the standard KHM equation. Finally, we note that the agreement between the 1D ST and KHM results appears for kl = α, where α ∼ 1.

Zoom In Zoom Out Reset image size

In this paper we derived two new forms of the KHM equation for compressible Hall MHD turbulence. We tested these equations, along with an isotropic ST equation, on the results of a 2D Hall MHD simulation of weakly compressible turbulence with a moderate Reynolds number. The KHM and ST equations are well satisfied in the simulation cross-validating these equations and simulation results. The two KHM equations give the same results, and they are equivalent and complementary to the ST equation via the inverse proportionality kl = α with $\alpha \simeq \sqrt{2}$. Note that a similar correspondence is observed in 3D HD turbulence (Hellinger et al. 2021) for $\alpha \simeq \sqrt{3}$, and preliminary results of isotropic (unmagnetized) 3D Hall MHD turbulence suggest the same relationship. In the reduced, 1D analyses we observe α ≃ 1. This indicates that the relationship between the ST and KHM equation kl = α depends on the dimensionality d of the analyzed space, $\alpha \simeq \sqrt{d}$. This simple relationship is useful to interpret the KMH results in the context of spectral analyses. Equations analogous to the ST and KHM ones can be obtained using the low-pass filtering/coarse graining of the energy conservation laws (Eyink & Aluie 2009; Aluie 2011, 2013; Yang et al. 2016, 2017; Camporeale et al. 2018).

The KHM and ST equations are valid during the whole simulation, owing to the periodic boundary conditions, and they can be used to analyze the onset of turbulence. The agreement between the two approaches is at least qualitative during the whole simulation, but only at later times, when turbulence is well developed, is there a good quantitative agreement; this indicates that the KHM and ST approaches become equivalent when turbulence is well developed and the energy transfer is expected to be local so that one can talk about the energy cascade. The simulation results suggest that the onset of the Hall cascade is related to the formation of reconnecting current sheets (Papini et al. 2019a) and that, at least during the onset, the Hall term leads to an energy transfer from small to large scales at a range of scales above d_i (cf. Franci et al. 2017); at this stage, this transfer is likely not local, and one cannot speak about an energy cascade. The locality can be tested using a (spectral) shell-to-shell transfer analysis (Alexakis et al. 2005). For the weakly compressible simulation, the incompressible approximation is not correct to a large extent owing to the neglected pressure dilatation effect. This error is important during the turbulence onset. When turbulence is well developed, the time-averaged pressure dilatation becomes negligible, and the compressible and incompressible KHM predictions for MHD and Hall cascade rates are close to each other.

It is interesting to note that the standard compressible KHM equation, Equation (4), does not depend (except for the correction term ${ \mathcal C }$) on the background magnetic field B ₀ = 〈 B 〉, similarly to the incompressible case (Oughton et al. 2013). On the other hand, there is in principle a contribution from the mean fluid velocity u ₀ = 〈 u 〉 (Hadid et al. 2017). In the incompressible, constant-density limit one recovers the incompressible results (Banerjee & Galtier 2017; Hellinger et al. 2018; Ferrand et al. 2019).

There are a couple of differences between our standard KHM equation (Equation (4)) and that of Andrés et al. (2018): We describe the kinetic energy by the structure function 〈∣δ w ∣²〉 that guarantees positive values in contrast with 〈δ(ρ u ) · δ u 〉. Similarly, we represent the magnetic energy contribution by the structure function S_B = 〈∣δ B ∣²〉 instead of 〈δ(ρ^1/2 B ) · δ( B /ρ^1/2)〉. Similarly, the alternative KHM equation (Equation (7)) is based on 〈∣δ w ∣²〉, whereas the work of Banerjee & Kritsuk (2018) corresponds to 〈δ(ρ u ) · δ u 〉. This difference is very likely not substantial. In the case of compressible HD turbulence, the corresponding two approaches give very similar results (Hellinger et al. 2021). However, the important difference between our results and those of Andrés et al. (2018) and of Banerjee & Kritsuk (2018) is that we formulate the KHM equations for the kinetic +magnetic energy. The simulation results exhibit no net energy exchange between the kinetic + magnetic energy and the internal one for later times when turbulence is well developed, and the pressure dilatation effect does not lead to cross-scale energy transfer (cf. Aluie et al. 2012). The inclusion of the internal energy in the KHM equation (Andrés et al. 2018; Banerjee & Kritsuk 2018, and references therein) is therefore not needed. Moreover, the HD results of Hellinger et al. (2021) indicate that it is hard to represent the kinetic energy and the internal one by compatible structure functions. Furthermore, Andrés et al. (2018) and Banerjee & Kritsuk (2018) use the isothermal closure to modify the description of the pressure dilatation effect; this closure is not generally applicable. Finally, estimates of the energy cascade rate in the solar wind based on such approaches (Banerjee et al. 2016; Hadid et al. 2017; Andrés et al. 2019) include the scale redistribution of the internal energy and cannot be simply related to the plasma heating rates (cf. MacBride et al. 2008).

For our study, we used one 2D Hall MHD simulation of weakly compressible turbulence for a moderated Reynolds number and zero cross-helicity. We observe a qualitatively similar behavior in other simulations. Quantitatively, the results depend on the spatial resolution, the box size, and other parameters. A larger simulation box (with respect to d_i ) allows for a wider MHD range (depending on the injection scale), and a better spatial resolution (with respect to d_i ) allows for a wider Hall range (depending on the viscosity and resistivity). The present results are sufficient for testing and comparing the KHM and ST equations but need to be extended to more compressible and/or larger systems that extend further to the Hall range of scales, and/or to systems with larger cross (canonical) helicity, etc. Our results need to be extended to 3D in order to investigate the anisotropy of the turbulent cascade. Both the KHM equations can be naturally used (cf. Verdini et al. 2015). However, it is not clear how to extend the isotropic ST approach to an anisotropic situation (cf. Verma 2017); the low-pass filtering/coarse graining approaches also usually assume isotropy.

The Hall cascade rate in the simulation is a fraction of the MHD one because of the dissipation. In both kinetic simulations and in situ observations the situation is similar (Bandyopadhyay et al. 2020b). The decrease of the cascade rate from the MHD to the Hall range is likely due to some sort of dissipation/particle energization. The present work assumes a scalar pressure p that is relevant for collision-dominated plasmas. In weakly collisional/collisionless systems, such as in the solar wind, it is necessary to employ the full pressure tensor P and to replace the pressure dilatation coupling, p θ, by the pressure-strain one, P : Σ (Yang et al. 2019; Matthaeus et al. 2020). The KHM and ST approaches can be easily generalized to account for the pressure tensor. For instance, 〈δ p δ θ〉 becomes simply 〈δ P : δ Σ〉 (see the viscous dissipation term).

Our reduced 1D analysis suggests that the ST equation and the alternative KHM equation give better estimations of the cascade rate from in situ observed time series compared to the usually used standard KHM. However, this is likely owing to the usage of multidimensional spatial derivatives that may possibly be only estimated using multipoint data of Cluster or MMS missions. Furthermore, the ST equation is formulated for an isotropic situation. Nevertheless, we believe that the ST and alternative KHM approaches are worth pursuing, as they provide other methods for measuring the energy cascade rates from in situ observations.

P.H. acknowledges grant 18-08861S of the Czech Science Foundation. This research was conducted with high-performance computing (HPC) resources provided by the Cineca ISCRA initiative (grant HP10C2EARF). L.F. is supported by the UK Science and Technology Facilities Council (STFC) grant ST/T00018X/1.

We start by reformulating the compressible Hall MHD equation in terms of w = ρ^1/2 u :

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{w}}}{\partial t}+({\boldsymbol{u}}\cdot {\boldsymbol{\nabla }}){\boldsymbol{w}}+\displaystyle \frac{{\boldsymbol{w}}\theta }{2}=\displaystyle \frac{{\boldsymbol{f}}}{\sqrt{\rho }}\end{eqnarray} \tag{ A1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}+({\boldsymbol{u}}\cdot {\boldsymbol{\nabla }}){\boldsymbol{B}}=({\boldsymbol{B}}\cdot {\boldsymbol{\nabla }}){\boldsymbol{u}}-{\boldsymbol{B}}\theta +\eta {\rm{\Delta }}{\boldsymbol{B}}+{\boldsymbol{h}},\end{eqnarray} \tag{ A2 }$$

where we denoted

$$\begin{eqnarray}&&{\boldsymbol{f}}=({\boldsymbol{\nabla }}\times {\boldsymbol{B}})\times {\boldsymbol{B}}-{\boldsymbol{\nabla }}p+{\boldsymbol{\nabla }}\cdot {\boldsymbol{\tau }}\end{eqnarray} \tag{ A3 }$$

$$\begin{eqnarray}&&{\boldsymbol{h}}=-({\boldsymbol{B}}\cdot {\boldsymbol{\nabla }}){\boldsymbol{j}}+({\boldsymbol{j}}\cdot {\boldsymbol{\nabla }}){\boldsymbol{B}}+{\boldsymbol{B}}\left({\boldsymbol{\nabla }}\cdot {\boldsymbol{j}}\right).\end{eqnarray} \tag{ A4 }$$

We take the Hall MHD equations at two different points, ${{\boldsymbol{x}}}^{{\prime} }={\boldsymbol{x}}+{\boldsymbol{l}}$ and x , and subtract them to get

$$\begin{eqnarray}&&\,\begin{array}{l}\displaystyle \frac{\partial \delta {\boldsymbol{w}}}{\partial t}+(\delta {\boldsymbol{u}}\cdot {\boldsymbol{\nabla }}^{\prime} )\delta {\boldsymbol{w}}+\left[{\boldsymbol{u}}\cdot \left({\boldsymbol{\nabla }}+{\boldsymbol{\nabla }}^{\prime} \right)\right]\delta {\boldsymbol{w}}\\ \,+\,\displaystyle \frac{1}{2}({{\boldsymbol{w}}}^{{\prime} }{\theta }^{{\prime} }-{\boldsymbol{w}}\theta )=\delta \left(\displaystyle \frac{{\boldsymbol{f}}}{\sqrt{\rho }}\right),\end{array}\,\end{eqnarray} \tag{ A5 }$$

where we denote a = a( x ), ${a}^{{\prime} }=a({{\boldsymbol{x}}}^{{\prime} })$, and $\delta a={a}^{{\prime} }-a$ for any variable a; analogically, ∇ = ∇ _x , ${{\boldsymbol{\nabla }}}^{{\prime} }={{\boldsymbol{\nabla }}}_{{{\boldsymbol{x}}}^{{\prime} }}$ and Δ = Δ _x , ${{\rm{\Delta }}}^{{\prime} }={{\rm{\Delta }}}_{{{\boldsymbol{x}}}^{{\prime} }}$, Δ being the Laplace operator. Here and henceforth this relationship is repeatedly applied (Carbone et al. 2009),

$$\begin{eqnarray}\,\begin{array}{rcl}\delta \left[({\boldsymbol{a}}\cdot {\boldsymbol{\nabla }}){\boldsymbol{b}}\right] & = & ({{\boldsymbol{a}}}^{{\prime} }\cdot {\boldsymbol{\nabla }}^{\prime} ){{\boldsymbol{b}}}^{{\prime} }-({\boldsymbol{a}}\cdot {\boldsymbol{\nabla }}){\boldsymbol{b}}\\ & = & ({{\boldsymbol{a}}}^{{\prime} }\cdot {\boldsymbol{\nabla }}^{\prime} )\delta {\boldsymbol{b}}+({\boldsymbol{a}}\cdot {\boldsymbol{\nabla }})\delta {\boldsymbol{b}}\\ & = & (\delta {\boldsymbol{a}}\cdot {\boldsymbol{\nabla }}^{\prime} )\delta {\boldsymbol{b}}+[{\boldsymbol{a}}\cdot ({\boldsymbol{\nabla }}+{\boldsymbol{\nabla }}^{\prime} )]\delta {\boldsymbol{b}},\end{array}\,\end{eqnarray} \tag{ A6 }$$

and x and ${{\boldsymbol{x}}}^{{\prime} }$ are assumed to be independent, i.e., ${\boldsymbol{\nabla }}^{\prime} {\boldsymbol{a}}\,={\boldsymbol{\nabla }}{{\boldsymbol{a}}}^{{\prime} }=0$ for any a ,

$$\begin{eqnarray}&&\displaystyle \frac{\partial | \delta {\boldsymbol{w}}{| }^{2}}{\partial t}+{\boldsymbol{\nabla }}^{\prime} \cdot \left(\delta {\boldsymbol{u}}| \delta {\boldsymbol{w}}{| }^{2}\right)+\left({\boldsymbol{\nabla }}+{\boldsymbol{\nabla }}^{\prime} \right)\cdot \left({\boldsymbol{u}}| \delta {\boldsymbol{w}}{| }^{2}\right)\end{eqnarray} \tag{ A7 }$$

$$\begin{eqnarray}&&+\delta {\boldsymbol{w}}\cdot \left({\theta }^{{\prime} }{\boldsymbol{w}}-\theta {{\boldsymbol{w}}}^{{\prime} }\right)=2\delta {\boldsymbol{w}}\cdot \delta \left(\displaystyle \frac{{\boldsymbol{f}}}{\sqrt{\rho }}\right).\end{eqnarray} \tag{ A8 }$$

Taking a spatial average,

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{ \mathcal S }}_{w}}{\partial t}+{{\boldsymbol{\nabla }}}_{{\boldsymbol{l}}}\cdot \left\langle \delta {\boldsymbol{u}}| \delta {\boldsymbol{w}}{| }^{2}\right\rangle +{ \mathcal R }=2\left\langle \delta {\boldsymbol{w}}\cdot \delta \left(\displaystyle \frac{{\boldsymbol{f}}}{\sqrt{\rho }}\right)\right\rangle ,\end{eqnarray} \tag{ A9 }$$

where ${{ \mathcal S }}_{w}=\langle | \delta {\boldsymbol{w}}{| }^{2}\rangle$ and ${ \mathcal R }=\langle \delta {\boldsymbol{w}}\cdot \left({\theta }^{{\prime} }{\boldsymbol{w}}-\theta {{\boldsymbol{w}}}^{{\prime} }\right)\rangle$. Here we assume that the system is homogeneous, i.e.,

$$\begin{eqnarray}&&\langle \left((,{\boldsymbol{\nabla }}+{\boldsymbol{\nabla }}^{\prime} \right)\cdot {\boldsymbol{a}}\rangle =0\end{eqnarray} \tag{ A10 }$$

for any quantity a (Frisch 1995).

Similarly for the magnetic field, we have

$$\begin{eqnarray}&&\,\begin{array}{l}\displaystyle \frac{\partial | \delta {\boldsymbol{B}}{| }^{2}}{\partial t}+{\boldsymbol{\nabla }}^{\prime} \cdot \left(\delta {\boldsymbol{u}}| \delta {\boldsymbol{B}}{| }^{2}\right)+\left({\boldsymbol{\nabla }}+{\boldsymbol{\nabla }}^{\prime} \right)\cdot \left({\boldsymbol{u}}| \delta {\boldsymbol{B}}{| }^{2}\right)\\ \quad =\,2\delta {\boldsymbol{B}}\cdot (\delta {\boldsymbol{B}}\cdot {\boldsymbol{\nabla }}^{\prime} )\delta {\boldsymbol{u}}+\left({B}^{{\prime} 2}-{B}^{2}\right)\left(\theta -{\theta }^{{\prime} }\right)\\ \qquad +\,2\delta {\boldsymbol{B}}\cdot \left[{\boldsymbol{B}}\cdot \left({\boldsymbol{\nabla }}+{\boldsymbol{\nabla }}^{\prime} \right)\right])\delta {\boldsymbol{u}}\\ \qquad +\,2\eta \delta {\boldsymbol{B}}\cdot \left({{\rm{\Delta }}}^{{\prime} }{{\boldsymbol{B}}}^{{\prime} }-{\rm{\Delta }}{\boldsymbol{B}}\right)+2\delta {\boldsymbol{B}}\cdot \delta {\boldsymbol{h}}.\end{array}\,\end{eqnarray} \tag{ A11 }$$

Taking the spatial average,

$$\begin{eqnarray}&&\,\begin{array}{l}\displaystyle \frac{\partial {{ \mathcal S }}_{B}}{\partial t}+{{\boldsymbol{\nabla }}}_{{\boldsymbol{l}}}\cdot \left\langle \delta {\boldsymbol{u}}| \delta {\boldsymbol{B}}{| }^{2}\right\rangle -2\left\langle \delta {\boldsymbol{B}}\cdot (\delta {\boldsymbol{B}}\cdot {\boldsymbol{\nabla }}^{\prime} )\delta {\boldsymbol{u}}\right\rangle \\ \quad =\,-\left\langle \delta \left({B}^{2}\right)\delta \theta \right\rangle +2\left\langle \delta {\boldsymbol{B}}\cdot \left[{\boldsymbol{B}}\cdot \left({\boldsymbol{\nabla }}+{\boldsymbol{\nabla }}^{\prime} \right)\right])\delta {\boldsymbol{u}}\right\rangle \\ \qquad -\,4{Q}_{\eta }+2\eta {{\rm{\Delta }}}_{{\boldsymbol{l}}}{S}_{B}+2\left\langle \delta {\boldsymbol{B}}\cdot \delta {\boldsymbol{h}}\right\rangle ,\end{array}\,\end{eqnarray} \tag{ A12 }$$

where S_B = 〈∣δ B ∣²〉, Q_η = η〈∇ B : ∇ B 〉 = η〈J²〉. In order to derive an equation for the structure function ${ \mathcal S }={{ \mathcal S }}_{w}+{{ \mathcal S }}_{B}$ and to obtain a simple equation, we introduce a correction term ${{ \mathcal C }}_{\mathrm{MHD}}=2\langle \delta {\boldsymbol{w}}\cdot \delta \left({\rho }^{-1/2}{\boldsymbol{f}}\right)\rangle -2\langle \delta {\boldsymbol{u}}\cdot \delta {\boldsymbol{f}}\rangle .$ The last term can be expressed as

$$\begin{eqnarray}\,\begin{array}{rcl}2\left\langle \delta {\boldsymbol{u}}\cdot \delta {\boldsymbol{f}}\right\rangle & = & 2\left\langle \delta {\boldsymbol{u}}\cdot \left({\boldsymbol{B}}\cdot \left({\boldsymbol{\nabla }}+{\boldsymbol{\nabla }}^{\prime} \right)\right)\delta {\boldsymbol{B}}\right\rangle \\ & & +\,2\left\langle \delta {\boldsymbol{u}}\cdot \left(\delta {\boldsymbol{B}}\cdot {\boldsymbol{\nabla }}^{\prime} \right)\delta {\boldsymbol{B}}\right\rangle -\left\langle \delta {\boldsymbol{u}}\cdot \delta \left({\boldsymbol{\nabla }}{B}^{2}\right)\right\rangle \\ & & +\,2\left\langle \delta \theta \delta p\right\rangle -2\left\langle \delta {\boldsymbol{\Sigma }}:\delta {\boldsymbol{\tau }}\right\rangle .\end{array}\,\end{eqnarray} \tag{ A13 }$$

Combining the previous results, we get for ${ \mathcal S }$

$$\begin{eqnarray}\,\begin{array}{rcl}\displaystyle \frac{\partial { \mathcal S }}{\partial t}+{{\boldsymbol{\nabla }}}_{{\boldsymbol{l}}}\cdot { \mathcal Y }+{ \mathcal R } & = & {{ \mathcal C }}_{\mathrm{MHD}}+2\left\langle \delta \theta \delta p\right\rangle -2\left\langle \delta {\boldsymbol{\Sigma }}:\delta {\boldsymbol{\tau }}\right\rangle \\ & & -\,4{Q}_{\eta }+2\eta {{\rm{\Delta }}}_{{\boldsymbol{l}}}{{ \mathcal S }}_{B}+2\langle \delta {\boldsymbol{B}}\cdot \delta {\boldsymbol{h}}\rangle ,\end{array}\,\end{eqnarray} \tag{ A14 }$$

where ${ \mathcal Y }=\left\langle \delta {\boldsymbol{u}}\left(| \delta {\boldsymbol{w}}{| }^{2}+| \delta {\boldsymbol{B}}{| }^{2}\right)-2\delta {\boldsymbol{B}}\left(\delta {\boldsymbol{B}}\cdot \delta {\boldsymbol{u}}\right)\right\rangle$.

To finish the calculation, we need to evaluate the Hall contribution 2〈δ B · δ h 〉. First, for δ h it is easy to get

$$\begin{eqnarray}\,\begin{array}{rcl}\delta {\boldsymbol{h}} & = & -(\delta {\boldsymbol{B}}\cdot {\boldsymbol{\nabla }}^{\prime} )\delta {\boldsymbol{j}}-({\boldsymbol{B}}\cdot \left({\boldsymbol{\nabla }}+{\boldsymbol{\nabla }}^{\prime} \right))\delta {\boldsymbol{j}}\\ & & +\,(\delta {\boldsymbol{j}}\cdot {\boldsymbol{\nabla }}^{\prime} )\delta {\boldsymbol{B}}+({\boldsymbol{j}}\cdot \left({\boldsymbol{\nabla }}+{\boldsymbol{\nabla }}^{\prime} \right))\delta {\boldsymbol{B}}\\ & & +\,\delta {\boldsymbol{B}}\left({\boldsymbol{\nabla }}^{\prime} \cdot \delta {\boldsymbol{j}}\right)+{\boldsymbol{B}}\left(\left[{\boldsymbol{\nabla }}+{\boldsymbol{\nabla }}^{\prime} \right]\cdot \delta {\boldsymbol{j}}\right).\end{array}\,\end{eqnarray} \tag{ A15 }$$

The term 2〈δ B · δ h 〉 may then be given in the form (Hellinger et al. 2018; Ferrand et al. 2019)

$$\begin{eqnarray}&&2\left\langle \delta {\boldsymbol{B}}\cdot \delta {\boldsymbol{h}}\right\rangle =-2{{\boldsymbol{\nabla }}}_{{\boldsymbol{l}}}\cdot { \mathcal H }+2{ \mathcal A },\end{eqnarray} \tag{ A16 }$$

where

$$\begin{eqnarray}\,\begin{array}{rcl}{ \mathcal H } & = & \left\langle \delta {\boldsymbol{B}}\left(\delta {\boldsymbol{B}}\cdot \delta {\boldsymbol{j}}\right)-\delta {\boldsymbol{j}}{\left|\delta {\boldsymbol{B}}\right|}^{2}/2\right\rangle \\ & = & \left\langle \delta {\boldsymbol{B}}\times \left(\delta {\boldsymbol{B}}\times \delta {\boldsymbol{j}}\right)+\delta {\boldsymbol{j}}{\left|\delta {\boldsymbol{B}}\right|}^{2}/2\right\rangle \end{array}\,\end{eqnarray} \tag{ A17 }$$

and

$$\begin{eqnarray}&&{ \mathcal A }=\left\langle {\boldsymbol{j}}\cdot \left({{\boldsymbol{B}}}^{{\prime} }\times {{\boldsymbol{J}}}^{{\prime} }\right)+{{\boldsymbol{j}}}^{{\prime} }\cdot \left({\boldsymbol{B}}\times {\boldsymbol{J}}\right)\right\rangle .\end{eqnarray} \tag{ A18 }$$

The term ${ \mathcal A }$ can be expressed as

$$\begin{eqnarray}\begin{array}{rcl}{ \mathcal A } & = & {{ \mathcal C }}_{{\rm{H}}}+\left\langle {\boldsymbol{J}}\cdot \left({{\boldsymbol{B}}}^{{\prime} }\times {{\boldsymbol{j}}}^{{\prime} }\right)+{{\boldsymbol{J}}}^{{\prime} }\cdot \left({\boldsymbol{B}}\times {\boldsymbol{j}}\right)\right\rangle \\ & = & {{ \mathcal C }}_{{\rm{H}}}+{\boldsymbol{\nabla }}^{\prime} \cdot \left[\left\langle \delta {\boldsymbol{B}}\times \left({{\boldsymbol{B}}}^{{\prime} }\times {{\boldsymbol{j}}}^{{\prime} }\right)+\delta {\boldsymbol{B}}\times \left({\boldsymbol{B}}\times {\boldsymbol{j}}\right)\right\rangle \right],\end{array}\,\end{eqnarray} \tag{ A19 }$$

where

$$\begin{eqnarray}\,\begin{array}{rcl}{{ \mathcal C }}_{{\rm{H}}} & = & \left\langle {\boldsymbol{j}}\cdot \left({{\boldsymbol{B}}}^{{\prime} }\times {{\boldsymbol{J}}}^{{\prime} }\right)+{{\boldsymbol{j}}}^{{\prime} }\cdot \left({\boldsymbol{B}}\times {\boldsymbol{J}}\right)\right\rangle \\ & & -\left\langle {\boldsymbol{J}}\cdot \left({{\boldsymbol{B}}}^{{\prime} }\times {{\boldsymbol{j}}}^{{\prime} }\right)+{{\boldsymbol{J}}}^{{\prime} }\cdot \left({\boldsymbol{B}}\times {\boldsymbol{j}}\right)\right\rangle .\end{array}\,\end{eqnarray} \tag{ A20 }$$

It can also be transformed into this form:

$$\begin{eqnarray}\,\begin{array}{rcl}{ \mathcal A } & = & -\left\langle {\boldsymbol{J}}\cdot \left({{\boldsymbol{B}}}^{{\prime} }\times {{\boldsymbol{j}}}^{{\prime} }\right)+{{\boldsymbol{J}}}^{{\prime} }\cdot \left({\boldsymbol{B}}\times {\boldsymbol{j}}\right)\right\rangle \\ & = & -{\boldsymbol{\nabla }}^{\prime} \cdot \left\langle \delta {\boldsymbol{B}}\times \left({{\boldsymbol{B}}}^{{\prime} }\times {\boldsymbol{j}}\right)\right\rangle -{\boldsymbol{\nabla }}^{\prime} \cdot \left\langle \delta {\boldsymbol{B}}\times \left({\boldsymbol{B}}\times {{\boldsymbol{j}}}^{{\prime} }\right)\right\rangle \\ & & -\left\langle \left({{\boldsymbol{j}}}^{{\prime} }\cdot {\boldsymbol{\nabla }}\right)| \delta {\boldsymbol{B}}{| }^{2}/2\right\rangle -\left\langle \left({\boldsymbol{j}}\cdot {\boldsymbol{\nabla }}^{\prime} \right)| \delta {\boldsymbol{B}}{| }^{2}/2\right\rangle .\end{array}\,\end{eqnarray} \tag{ A21 }$$

Combining Equation (A19) and Equation (A21) gives

$$\begin{eqnarray}&&2{ \mathcal A }={{ \mathcal C }}_{H}+{{\boldsymbol{\nabla }}}_{{\boldsymbol{l}}}\cdot { \mathcal H }.\end{eqnarray} \tag{ A22 }$$

The final version of the KHM equation in compressible Hall MHD gets this form:

$$\begin{eqnarray}\,\begin{array}{rcl}\displaystyle \frac{\partial { \mathcal S }}{\partial t}+{{\boldsymbol{\nabla }}}_{{\boldsymbol{l}}}\cdot \left({ \mathcal Y }+{ \mathcal H }\right)+{ \mathcal R } & = & { \mathcal C }+2\left\langle \delta \theta \delta p\right\rangle -2\left\langle \delta {\boldsymbol{\Sigma }}:\delta {\boldsymbol{\tau }}\right\rangle \\ & & -\,4{Q}_{\eta }+2\eta {{\rm{\Delta }}}_{{\boldsymbol{l}}}{{ \mathcal S }}_{B},\end{array}\,\end{eqnarray} \tag{ A23 }$$

with ${ \mathcal C }={{ \mathcal C }}_{\mathrm{MHD}}+{{ \mathcal C }}_{H}$.