# An in-depth numerical study of exact laws for compressible Hall
magnetohydrodynamic turbulence

R. Ferrand11affiliationmark:  , F. Sahraoui11affiliationmark:  , S.
Galtier11affiliationmark:  22affiliationmark:  , N. Andrés33affiliationmark:
44affiliationmark:  , P. Mininni33affiliationmark:  and P.
Dmitruk33affiliationmark<EMAIL_ADDRESS>1 Laboratoire
de Physique des Plasmas, CNRS, École polytechnique, Université Paris-Saclay,
Sorbonne Université, Observatoire de Paris-Meudon, F-91128 Palaiseau Cedex,
France  2 Institut Universitaire de France (IUF)  3 Departamento de Física,
UBA, Ciudad Universitaria, 1428, Buenos Aires, Argentina  4 Instituto de
Astronomía y Física del Espacio, CONICET-UBA, Ciudad Universitaria, 1428,
Buenos Aires, Argentina

###### Abstract

Various exact laws governing compressible magnetohydrodynamic (MHD) and Hall-
MHD (CHMHD) turbulence have been derived in recent years. Other than their
fundamental theoretical interest, these laws are generally used to estimate
the energy dissipation rate from spacecraft observations in order to address
diverse problems related, e.g., to heating of the solar wind (SW) and
magnetospheric plasmas. Here we use various $1024^{3}$ direct numerical
simulation (DNS) data of free-decay isothermal CHMHD turbulence obtained with
the GHOST code (Geophysical High-Order Suite for Turbulence) to analyze two of
the recently derived exact laws. The simulations reflect different intensities
of the initial Mach number and the background magnetic field. The analysis
demonstrates the equivalence of the two laws in the inertial range and relates
the strength of the Hall effect to the amplitude of the cascade rate at sub-
ion scales. When taken in their general form (i.e., not limited to the
inertial range) some subtleties regarding the validity of the stationarity
assumption or the absence of the forcing in the simulations are discussed. We
show that the free-decay nature of the turbulence induces a shift from a large
scale forcing towards the presence of a scale-dependent reservoir of energy
fueling the cascade or dissipation. The reduced form of the exact laws (valid
in the inertial range) ultimately holds even if the stationarity assumption is
not fully verified.

##  1 Introduction

Understanding the dynamics of turbulent magnetized flows has been a
longstanding problem in physics, and especially in astrophysics where
turbulence is thought to play a leading role in various physical processes.
Examples are the interstellar medium, the SW or planetary magnetospheres, in
which turbulence controls structures formation, energy whereabouts and
particle heating or acceleration (Matthaeus et al., 1999; Bruno & Carbone,
2013; Kritsuk et al., 2007; Arzoumanian et al., 2011; Sahraoui et al., 2020).
Due to the chaotic nature of turbulence, such media are often studied thanks
to the use of specific tools, which rely on statistical methods to uncover
trends in the behavior of turbulent flows. A prime example of such tools are
exact laws: these equations, which can be obtained through the sole hypothesis
of statistical homogeneity (and further refined by introducing time
stationarity and infinite Reynolds number) express the rate of energy flowing
towards the small/dissipative scales of a system as a function of two-point
structure functions, without requiring the use of closure models. Initiated by
the work of Kolmogorov and his four-fifths law for hydrodynamic turbulence
(Kolmogorov, 1941), the quest for exact laws has grown wider ever since with
more and more elaborate models being derived. The first steps into studying
plasma turbulence were taken by Politano & Pouquet (1998), who derived a law
for incompressible MHD turbulence. This result paved the way for more precise
studies of space plasmas (Sorriso-Valvo et al., 2007; MacBride et al., 2008;
Marino et al., 2008; Stawarz et al., 2009; Osman et al., 2011; Galtier, 2012).
More general exact laws have been derived subsequently, by considering the
influence of Hall physics in incompressible models (Galtier, 2008; Banerjee &
Galtier, 2017; Hellinger et al., 2018; Ferrand et al., 2019), or the
compression of the flow (Carbone et al., 2009; Galtier & Banerjee, 2011;
Banerjee & Galtier, 2013, 2014; Banerjee & Kritsuk, 2017; Andrés & Sahraoui,
2017; Andrés et al., 2018; Banerjee & Kritsuk, 2018; Lindborg, 2019; Ferrand
et al., 2021a; Simon & Sahraoui, 2021).

Thanks to those laws, several studies of astrophysical media have been made
possible, either through direct numerical simulations (DNSs) (Mininni &
Pouquet, 2009; Kritsuk et al., 2013; Verdini et al., 2015; Ferrand et al.,
2020) or in situ data analysis (Marino et al., 2011; Coburn et al., 2015;
Banerjee et al., 2016; Hadid et al., 2017, 2018; Andrés et al., 2019; Sorriso-
Valvo et al., 2019; Bandyopadhyay et al., 2020; Andrés et al., 2020). These
allowed testing the efficiency of the exact laws in practical situations in a
variety of different systems. Studies specifically designed to test the
validity of these exact laws were also conducted, allowing for a more in-depth
understanding of how the constituents of these equations relate to each other
(Carbone et al., 2009; Hellinger et al., 2018; Andrés et al., 2018). This kind
of work is especially important as the models are further refined, and exact
laws grow even more complex, such as those derived for compressible Hall MHD,
a model taking into account the small scale correction to the plasma dynamics
due to the Hall effect. For this last plasma model, where two different exact
laws have been derived, such a study has yet to be made. It is thus the aim of
this paper to fill this gap by proposing, through the use of an ensemble of
high-resolution DNSs for decaying compressible Hall MHD turbulence, a term-by-
term analysis of exact laws derived by Andrés et al. (2018) and Ferrand et al.
(2021a). The objective is not only to better understand the exact laws, but to
test their more general expressions (valid beyond the inertial range and
without time stationarity), which require the bare minimum of hypotheses, and
unveil the relations between the energy cascade, the dissipation, and the
dynamical variables that form these equations.

The structure of this paper is as follows: first we present the two models
from Andrés et al. (2018) and Ferrand et al. (2021a), which will be referred
to as A18 and F21 respectively, and provide compact and general expressions
requiring only the assumption of statistical homogeneity to be obtained. We
then describe the simulations and the numerical schemes used to compute all
the terms forming both exact laws, and then present and discuss the results
obtained. Finally, we give a conclusion on our work and on the behavior of the
exact laws.

##  2 Theoretical models

Prior to presenting the theory of the two exact laws we introduce some
notations and relations tied to the framework of both calculations. We
introduce the spatial increment $\boldsymbol{\ell}$, connecting two points r
and $\textbf{r}^{\prime}$ in the physical space as
$\textbf{r}^{\prime}=\textbf{r}+\boldsymbol{\ell}$, and define for any given
field $\xi$ : $\xi\equiv\xi(\textbf{r})$ and
$\xi^{\prime}\equiv\xi(\textbf{r}^{\prime})$. We also define the notations
$\delta\xi=\xi^{\prime}-\xi$, $\bar{\delta}\xi=\frac{1}{2}(\xi^{\prime}+\xi)$
and the differential operator in the direction $\boldsymbol{\ell}$ as
$\boldsymbol{\nabla}_{\boldsymbol{\ell}}$. This operator obeys the following
relation on ensemble averages $\left\langle\right\rangle$ :
$\langle\boldsymbol{\nabla}^{\prime}\cdot\rangle=-\langle\boldsymbol{\nabla}\cdot\rangle=\boldsymbol{\nabla}_{\boldsymbol{\ell}}\cdot\langle\rangle$.

###  2.1 F21 model

The isothermal CHMHD exact law derived by Ferrand et al. (2021a) is obtained
by considering the following three-dimensional compressible HMHD equations:

$\displaystyle\partial_{t}\rho+\boldsymbol{\nabla}\cdot(\rho\textbf{v})=$
$\displaystyle\,0\,,$ (1)
$\displaystyle\rho(\partial_{t}\textbf{v}+\textbf{v}\cdot\boldsymbol{\nabla}\textbf{v})=$
$\displaystyle-\boldsymbol{\nabla}P+\textbf{J}\times\textbf{B}+{\bf
d_{\nu}}+{\bf f}\,,$ (2) $\displaystyle\partial_{t}\textbf{B}=$
$\displaystyle\boldsymbol{\nabla}\times(\textbf{v}\times\textbf{B})-\lambda\boldsymbol{\nabla}\times(\textbf{J}_{\text{c}}\times\textbf{B})+{\bf
d_{\eta}}\,,$ (3) $\displaystyle\boldsymbol{\nabla}\cdot\textbf{B}=$
$\displaystyle\,0\,,$ (4)

where $\rho$ is the mass density, v the velocity, $P$ the pressure, B the
magnetic field, $\textbf{J}=(\boldsymbol{\nabla}\times\textbf{B})/\mu_{0}$ the
current density and $\textbf{J}_{\text{c}}=\textbf{J}/\rho$ the normalized
current. The Hall effect is introduced in the model through the presence of
the term
$\lambda\boldsymbol{\nabla}\times(\textbf{J}_{\text{c}}\times\textbf{B})$ in
the induction equation: $\lambda=m_{i}/q_{e}$, with $m_{i}$ the ion mass and
$q_{e}$ the magnitude of the electron charge, is the marker of the Hall effect
and is connected to the ion inertial length $d_{i}$ through the relation
$d_{i}=\lambda/\sqrt{\mu_{0}\rho}$, with $\mu_{0}$ the vacuum permeability.
The dissipation terms are:

$\displaystyle{\bf d_{\nu}}=$
$\displaystyle\,\nu\Delta\textbf{v}+\frac{\nu}{3}\boldsymbol{\nabla}\theta\,,$
(5) $\displaystyle{\bf d_{\eta}}=$ $\displaystyle\,\eta\Delta\textbf{B}\,,$
(6)

with $\theta=\boldsymbol{\nabla}\cdot\textbf{v}$ the dilatation, $\nu$ the
dynamic viscosity and $\eta$ the magnetic diffusivity. The term f represents a
large scale forcing. The isothermal closure writes $P=c_{s}^{2}\rho$ with
$c_{s}$ the constant speed of sound. These equations are used to derive a
dynamical equation for the modified second-order structure function

$\displaystyle\langle S\rangle\equiv$
$\displaystyle\left\langle\frac{1}{2}\bar{\delta}\rho|\delta\textbf{v}|^{2}+\frac{1}{2\mu_{0}}|\delta\textbf{B}|^{2}+\frac{1}{2}\delta\rho\delta
e\right\rangle.$ (7)

In the simulation code the internal energy follows the polytropic definition
$e=\frac{c_{s}^{2}}{\gamma(\gamma-1)}((\rho/\rho_{0})^{\gamma-1}-1)$ with a
polytropic index $\gamma=1.01$ close enough to unity so that the isothermal
approximation remains valid. Note that under this approximation entropy is
assumed to be constant, and thus the internal energy variation results only
from the work of the pressure force (Simon & Sahraoui, 2021).

Injecting equations (1)-(4) in the time derivative of (7) leads, after a hefty
amount of calculations, to the equation:

$\displaystyle\partial_{t}\langle S\rangle=$ $\displaystyle\partial_{t}\langle
E^{tot}\rangle-\frac{1}{2}\boldsymbol{\nabla}_{\boldsymbol{\ell}}\cdot\left\langle\bar{\delta}\rho|\delta\textbf{v}|^{2}\delta\textbf{v}\right\rangle+\frac{1}{4}\langle(\rho\theta^{\prime}+\rho^{\prime}\theta)|\delta\textbf{v}|^{2}\rangle-\langle\delta\rho\,\delta\textbf{v}\cdot\bar{\delta}(\textbf{J}_{\text{c}}\times\textbf{B})\rangle+{\lambda}\langle\delta(\textbf{J}\times\textbf{B})\cdot\delta\textbf{J}_{\text{c}}\rangle$
$\displaystyle-\frac{1}{2\mu_{0}}\boldsymbol{\nabla}_{\boldsymbol{\ell}}\cdot\left\langle|\delta\textbf{B}|^{2}\delta\textbf{v}-2(\delta\textbf{v}\cdot\delta\textbf{B})\delta\textbf{B}-\lambda|\delta\textbf{B}|^{2}\delta\textbf{J}_{\text{c}}+2\lambda(\delta\textbf{B}\cdot\delta\textbf{J}_{\text{c}})\delta\textbf{B}\right\rangle$
$\displaystyle-\frac{1}{2}\left\langle\left(1+\frac{\rho^{\prime}}{\rho}\right)\textbf{v}^{\prime}\cdot({\bf
d_{\nu}}+{\bf
f})+\left(1+\frac{\rho}{\rho^{\prime}}\right)\textbf{v}\cdot({\bf
d_{\nu}^{\prime}}+{\bf f^{\prime}})\right\rangle$
$\displaystyle+\frac{1}{2}\left\langle\frac{\rho}{\rho^{\prime}}\textbf{v}^{\prime}\cdot({\bf
d_{\nu}^{\prime}}+{\bf
f^{\prime}})+\frac{\rho^{\prime}}{\rho}\textbf{v}\cdot({\bf d_{\nu}}+{\bf
f})\right\rangle-\frac{1}{\mu_{0}}\langle\textbf{B}^{\prime}\cdot{\bf
d_{\eta}}+\textbf{B}\cdot{\bf d^{\prime}_{\eta}}-\textbf{B}\cdot{\bf
d_{\eta}}\rangle\,,$ (8)

with the total energy $E^{tot}=\rho u^{2}/2+B^{2}/(2\mu_{0})+\rho e$. Equation
(8), already reported in Ferrand et al. (2021a), is the most general equation
obtainable in this model under the sole assumption of statistical homogeneity.
With additional assumptions of time stationarity, forcing limited to the
largest scales and infinite (magnetic and kinetic) Reynolds numbers one can
then retrieve the exact law:

$\displaystyle-4\varepsilon_{F21}=$
$\displaystyle\boldsymbol{\nabla}_{\boldsymbol{\ell}}\cdot\left\langle\bar{\delta}\rho|\delta\textbf{v}|^{2}\delta\textbf{v}+\frac{1}{\mu_{0}}|\delta\textbf{B}|^{2}\delta\textbf{v}-\frac{2}{\mu_{0}}(\delta\textbf{v}\cdot\delta\textbf{B})\delta\textbf{B}\right\rangle-\frac{1}{2}\langle(\rho\theta^{\prime}+\rho^{\prime}\theta)|\delta\textbf{v}|^{2}\rangle+2\langle\delta\rho\,\delta\textbf{v}\cdot\bar{\delta}(\textbf{J}_{\text{c}}\times\textbf{B})\rangle$
$\displaystyle+$
$\displaystyle\frac{\lambda}{\mu_{0}}\boldsymbol{\nabla}_{\boldsymbol{\ell}}\cdot\left\langle
2(\delta\textbf{B}\cdot\delta\textbf{J}_{\text{c}})\delta\textbf{B}-|\delta\textbf{B}|^{2}\delta\textbf{J}_{\text{c}}\right\rangle-2{\lambda}\langle\delta(\textbf{J}\times\textbf{B})\cdot\delta\textbf{J}_{\text{c}}\rangle.$
(9)

Equation (9) (or similar ones obtained for other models) is generally the one
that is used in numerical simulations and spacecraft data to infer the cascade
rate (left-hand side) from measurable quantities (right-hand side) (Banerjee &
Galtier, 2016; Hadid et al., 2017, 2018; Andrés et al., 2018, 2019). However,
to be obtained, non-trivial approximations had to be made, about time
stationarity and on the forcing and dissipative terms (Ferrand et al., 2021a).
While these approximation are hard (if not impossible) to test on spacecraft
data, they can a priori be verified in DNSs, which is one of the goals of the
present study. The specific point about time stationarity and the absence of
driving is of a particular relevance in free-decay simulations as the ones we
are using in this study. Therefore, we will consider in the following the more
general equation (8), written in a more compact form:

$2\varepsilon_{F21}+\partial_{t}\left\langle
E^{tot}\right\rangle-\partial_{t}\left\langle
S\right\rangle+\mathcal{D}_{F21}+\mathcal{D}_{F21}^{loc}+\mathcal{F}_{F21}+\mathcal{F}_{F21}^{loc}=0$
(10)

where we define:

$\displaystyle\mathcal{D}_{F21}\equiv$
$\displaystyle-\frac{1}{2}\left\langle\left(1+\frac{\rho^{\prime}}{\rho}\right)\textbf{v}^{\prime}\cdot{\bf
d_{\nu}}+\left(1+\frac{\rho}{\rho^{\prime}}\right)\textbf{v}\cdot{\bf
d_{\nu}^{\prime}}\right\rangle-\frac{1}{\mu_{0}}\left\langle\textbf{B}^{\prime}\cdot{\bf
d_{\eta}}+\textbf{B}\cdot{\bf d^{\prime}_{\eta}}\right\rangle,$ (11)
$\displaystyle\mathcal{D}_{F21}^{loc}\equiv$
$\displaystyle~{}+\frac{1}{2}\left\langle\frac{\rho}{\rho^{\prime}}\textbf{v}^{\prime}\cdot{\bf
d_{\nu}^{\prime}}+\frac{\rho^{\prime}}{\rho}\textbf{v}\cdot{\bf
d_{\nu}}\right\rangle+\frac{1}{2\mu_{0}}\left\langle\textbf{B}\cdot{\bf
d_{\eta}}+\textbf{B}^{\prime}\cdot{\bf d^{\prime}_{\eta}}\right\rangle,$ (12)
$\displaystyle\mathcal{F}_{F21}\equiv$
$\displaystyle-\frac{1}{2}\left\langle\left(1+\frac{\rho^{\prime}}{\rho}\right)\textbf{v}^{\prime}\cdot{\bf
f}+\left(1+\frac{\rho}{\rho^{\prime}}\right)\textbf{v}\cdot{\bf
f^{\prime}}\right\rangle,$ (13) $\displaystyle\mathcal{F}_{F21}^{loc}\equiv$
$\displaystyle+\frac{1}{2}\left\langle\frac{\rho}{\rho^{\prime}}\textbf{v}^{\prime}\cdot{\bf
f^{\prime}}+\frac{\rho^{\prime}}{\rho}\textbf{v}\cdot{\bf f}\right\rangle.$
(14)

Note that in equation (10) $\varepsilon_{F21}$ is introduced by identifying it
in equation (8) to its expression given by relation (9). It does not result
from using the same assumptions that led to relation (9), and thus equation
(10) remains very general. The superscript loc is used for the forcing and
dissipation terms that result from the product of vectors taken locally (i.e.,
at the same position r or $\textbf{r}^{\prime}$), assuming that
$\rho\sim\rho^{\prime}$ for the velocity field terms. This contrasts with the
terms $\mathcal{D}_{F21}$ and $\mathcal{F}_{F21}$ that involve (distant) two-
point correlations. Assuming that $\rho\sim\rho^{\prime}$ in weakly
compressible simulations such as those of this study, the terms
$\mathcal{D}_{F21}$ and $\mathcal{F}_{F21}$ are likely to be scale-
independent. However, as we will show below they still have a significant
impact on the energy balance. The same remark can be made about the term
$\partial_{t}\left\langle E^{tot}\right\rangle$ in Eq. 10, which is clearly a
local term (unlike $\partial_{t}\left\langle S\right\rangle$ since $S$ is a
second order structure function).

###  2.2 A18 model

The other exact law derived by Andrés et al. (2018) relies on the same three-
dimensional isothermal HMHD equations, yet using the Alfvén speed
$\textbf{v}_{A}\equiv\textbf{B}/\sqrt{\mu_{0}\rho}$ instead of the magnetic
field:

$\displaystyle\partial_{t}\rho=$
$\displaystyle-\boldsymbol{\nabla}\cdot(\rho\textbf{v}),$ (15)
$\displaystyle\partial_{t}\textbf{v}=$
$\displaystyle-\textbf{v}\cdot\boldsymbol{\nabla}\textbf{v}+\textbf{v}_{\text{A}}\cdot\boldsymbol{\nabla}\textbf{v}_{\text{A}}-\frac{1}{\rho}\boldsymbol{\nabla}(P+P_{M})$
$\displaystyle-\textbf{v}_{\text{A}}\cdot(\boldsymbol{\nabla}\cdot\textbf{v}_{\text{A}})+{\bf
d_{\nu}}+{\bf f},$ (16) $\displaystyle\partial_{t}\textbf{v}_{\text{A}}=$
$\displaystyle-(\textbf{v}-\lambda\textbf{J}_{\text{c}})\cdot\boldsymbol{\nabla}\textbf{v}_{\text{A}}+\textbf{v}_{\text{A}}\cdot\boldsymbol{\nabla}(\textbf{v}-\lambda\textbf{J}_{\text{c}})$
$\displaystyle-\frac{\textbf{v}_{\text{A}}}{2}(\boldsymbol{\nabla}\cdot\textbf{v}-\lambda\boldsymbol{\nabla}\cdot\textbf{J}_{\text{c}})+\frac{{\bf
d_{\eta}}}{\sqrt{\rho}},$ (17)
$\displaystyle\textbf{v}_{\text{A}}\cdot\boldsymbol{\nabla}\rho=$
$\displaystyle-2\rho(\boldsymbol{\nabla}\cdot\textbf{v}_{\text{A}}),$ (18)
$\displaystyle\textbf{J}_{\text{c}}\cdot\boldsymbol{\nabla}\rho=$
$\displaystyle-\rho(\boldsymbol{\nabla}\cdot\textbf{J}_{\text{c}}),$ (19)

where $P_{M}\equiv\rho v_{\text{A}}^{2}/2$ is the magnetic pressure and $\rho$
is a time-space variable. The system is once again closed with an isothermal
closure. These equations are used to compute the time derivative of the two
point correlator
$R_{E}\equiv~{}\frac{\rho}{2}(\textbf{v}\cdot\textbf{v}^{\prime}+\textbf{v}_{\text{A}}\cdot\textbf{v}_{\text{A}}^{\prime})+\rho
e^{\prime}$, which ultimately leads to the exact law:

$\displaystyle-2\varepsilon_{A18}=$
$\displaystyle~{}\frac{1}{2}\boldsymbol{\nabla}_{\ell}\cdot\left\langle[(\delta(\rho\textbf{v})\cdot\delta\textbf{v}+\delta(\rho\textbf{v}_{\text{A}})\cdot\delta\textbf{v}_{\text{A}}+2\delta
e\delta\rho\big{]}\delta\textbf{v}-[\delta(\rho\textbf{v})\cdot\delta\textbf{v}_{\text{A}}+\delta\textbf{v}\cdot\delta(\rho\textbf{v}_{\text{A}})]\delta\textbf{v}_{\text{A}}\right\rangle$
$\displaystyle+$
$\displaystyle\left\langle[R_{E}^{\prime}-\frac{1}{2}(R_{B}^{\prime}+R_{B})+\frac{P_{M}^{\prime}-P^{\prime}}{2}-E^{tot^{\prime}}](\boldsymbol{\nabla}\cdot\textbf{v})\right\rangle+\left\langle[R_{E}-\frac{1}{2}(R_{B}+R_{B}^{\prime})+\frac{P_{M}-P}{2}-E^{tot}](\boldsymbol{\nabla}^{\prime}\cdot\textbf{v}^{\prime})\right\rangle$
$\displaystyle+$
$\displaystyle\left\langle[(R_{H}-R_{H}^{\prime})-\bar{\rho}(\textbf{v}^{\prime}\cdot\textbf{v}_{\text{A}})+H^{\prime}+\lambda\delta\rho\frac{\textbf{J}_{\text{c}}\cdot\textbf{v}_{\text{A}}^{\prime}}{2}](\boldsymbol{\nabla}\cdot\textbf{v}_{\text{A}})\right\rangle$
$\displaystyle+$
$\displaystyle\left\langle[(R_{H}^{\prime}-R_{H})-\bar{\rho}(\textbf{v}\cdot\textbf{v}_{\text{A}}^{\prime})+H-\lambda\delta\rho\frac{\textbf{J}_{\text{c}}^{\prime}\cdot\textbf{v}_{\text{A}}}{2}](\boldsymbol{\nabla}^{\prime}\cdot\textbf{v}_{\text{A}}^{\prime})\right\rangle$
$\displaystyle+$
$\displaystyle\frac{1}{2}\left\langle\big{(}e^{\prime}+\frac{v_{\text{A}}}{2}^{{}^{\prime}2}\big{)}\big{[}\boldsymbol{\nabla}\cdot(\rho\textbf{v})\big{]}+\big{(}e+\frac{v_{\text{A}}}{2}^{2}\big{)}\big{[}\boldsymbol{\nabla}^{\prime}\cdot(\rho^{\prime}\textbf{v}^{\prime})\big{]}\right\rangle-\frac{1}{2}\left\langle\beta^{-1^{\prime}}\boldsymbol{\nabla}^{\prime}\cdot(e^{\prime}\rho\textbf{v})+\beta^{-1}\boldsymbol{\nabla}\cdot(e\rho^{\prime}\textbf{v}^{\prime})\right\rangle$
$\displaystyle+$
$\displaystyle\frac{1}{2}\boldsymbol{\nabla}_{\ell}\cdot\left\langle
2\lambda[(\overline{\rho\textbf{J}_{\text{c}}\times\textbf{v}_{\text{A}}})\times\delta\textbf{v}_{\text{A}}-\delta(\textbf{J}_{\text{c}}\times\textbf{v}_{\text{A}})\times\overline{\rho\textbf{v}_{\text{A}}}]\right\rangle+\lambda\left\langle\frac{R_{B}-R_{B}^{\prime}}{2}(\boldsymbol{\nabla}\cdot\textbf{J}_{\text{c}})+\frac{R_{B}^{\prime}-R_{B}}{2}(\boldsymbol{\nabla}^{\prime}\cdot\textbf{J}_{\text{c}}^{\prime})\right\rangle,$
(20)

with the cross helicity
$H\equiv\rho(\textbf{v}\cdot\nobreak\textbf{v}_{\text{A}})$, its two-points
correlator
$R_{H}\equiv\rho(\textbf{v}\cdot\textbf{v}_{\text{A}}^{\prime}+\textbf{v}_{\text{A}}\cdot\textbf{v}^{\prime})/2$,
the correlator for magnetic energy
$R_{B}\equiv\rho\textbf{v}_{\text{A}}\cdot\textbf{v}_{\text{A}}^{\prime}/2$
and $\beta^{-1}\equiv v_{\text{A}}^{2}/2c_{s}^{2}$. Primed variables are
obtained by inverting the positions of the primes in the definitions. Just
like equation (9), this law is obtained under the assumptions of statistical
homogeneity, forcing limited to large scales, time stationarity and infinite
Reynolds number. Again, the assumptions on the forcing, dissipation and time
stationarity can be dropped by using the general expression obtained directly
from the dynamical equation of $\left\langle
R_{E}+R_{E}^{\prime}\right\rangle$:

$2\varepsilon_{A18}+\partial_{t}\left\langle
R_{E}+R_{E}^{\prime}\right\rangle+\mathcal{D}_{A18}+\mathcal{F}_{A18}=0.$ (21)

Here again the same caution as above applies regarding the introduction of
$\varepsilon_{A18}$, whose expression is given by (20). The dissipative and
forcing terms $\mathcal{D}_{A18}$ and $\mathcal{F}_{A18}$ are not explicitly
given in Andrés et al. (2018) but are easy to calculate. For the dissipation
term, the component stemming from the velocity field and Navier-Stokes
equation is the same as in (11), and the one originating from the magnetic
field and induction equation stems from:

$\partial_{t}\left\langle\frac{1}{2}(\rho+\rho^{\prime})\textbf{v}_{\text{A}}\cdot\textbf{v}_{\text{A}}^{\prime}\right\rangle,$
(22)

and reads:

$\displaystyle\frac{1}{2}\left\langle\rho\textbf{v}_{\text{A}}\cdot\frac{\textbf{d}_{\eta}^{\prime}}{\sqrt{\rho^{\prime}}}+\rho\textbf{v}_{\text{A}}^{\prime}\cdot\frac{\textbf{d}_{\eta}}{\sqrt{\rho}}+\rho^{\prime}\textbf{v}_{\text{A}}\cdot\frac{\textbf{d}_{\eta}^{\prime}}{\sqrt{\rho^{\prime}}}+\rho^{\prime}\textbf{v}_{\text{A}}^{\prime}\cdot\frac{\textbf{d}_{\eta}}{\sqrt{\rho}}\right\rangle$
$\displaystyle=$
$\displaystyle\frac{1}{2\mu_{0}}\left\langle(\frac{\sqrt{\rho}}{\sqrt{\rho^{\prime}}}+\frac{\sqrt{\rho^{\prime}}}{\sqrt{\rho}})(\textbf{B}\cdot\textbf{d}_{\eta}^{\prime}+\textbf{B}^{\prime}\cdot\textbf{d}_{\eta})\right\rangle,$
(23)

which ultimately leads to

$\displaystyle\mathcal{D}_{A18}=$
$\displaystyle~{}-\frac{1}{2}\left\langle\left(1+\frac{\rho^{\prime}}{\rho}\right)\textbf{v}^{\prime}\cdot{\bf
d_{\nu}}+\left(1+\frac{\rho}{\rho^{\prime}}\right)\textbf{v}\cdot{\bf
d_{\nu}^{\prime}}\right\rangle$ $\displaystyle-$
$\displaystyle~{}\frac{1}{2\mu_{0}}\left\langle(\frac{\sqrt{\rho}}{\sqrt{\rho^{\prime}}}+\frac{\sqrt{\rho^{\prime}}}{\sqrt{\rho}})(\textbf{B}\cdot\textbf{d}_{\eta}^{\prime}+\textbf{B}^{\prime}\cdot\textbf{d}_{\eta})\right\rangle.$
(24)

The forcing term is identical to the one of law F21, and reads

$\displaystyle\mathcal{F}_{A18}=$
$\displaystyle~{}-\frac{1}{2}\left\langle\left(1+\frac{\rho^{\prime}}{\rho}\right)\textbf{v}^{\prime}\cdot{\bf
f}+\left(1+\frac{\rho}{\rho^{\prime}}\right)\textbf{v}\cdot{\bf
f^{\prime}}\right\rangle.$ (25)

##  3 Numerical methods

###  3.1 Simulation data

The equations of CHMHD (1)-(4) are solved numerically using the pseudo-
spectral code GHOST (Gómez et al., 2005; Mininni et al., 2011) along with a
module for solving compressible HMHD flows with the eventual presence of a
background magnetic field. Three simulations were run in a cubic periodic box
of spatial resolution of $N=1024$ grid points and size $L_{0}=2\pi$ in all
three directions, and all simulations use dimensionless viscosity and magnetic
diffusivity of $\nu=\eta=3\times 10^{-4}$. Due to the normalization of the
variables, the corresponding Reynolds numbers are simply the inverse of the
viscosity and the diffusivity: $R_{e}=R_{m}\simeq 3\times 10^{3})$. The ion
inertial length $d_{i}$ is set to $d_{i}=0.02L_{0}$ for the three Runs. These
simulations do not feature any forcing: instead, the kinetic and magnetic
fields are set to an initial state built from a superposition of harmonic
modes with random phases whose energy in Fourier space is put in a sphere
between wave vectors $k_{down}$ and $k_{up}$, following the idea of Pouquet &
Patterson (1978). Then, the flow is left to evolve and decay. Varying the
intensity of background magnetic field $B_{0}$ and the initial Mach number
$M_{S}$ allows us to evaluate the influence of these parameters on the
dynamics of the system. All aforementioned parameters are reported in Table 1.

Run |  Resolution |  $B_{0}$ |  $M_{S}$ |  $d_{i}/L_{0}$ |  $\nu=\eta$ |  $k_{down}$ |  $k_{up}$
---|---|---|---|---|---|---|---
I |  $1024^{3}$ |  $2$ |  0.25 |  0.02 |  $3.0\times 10^{-4}$ |  1 |  3
II |  $1024^{3}$ |  $0$ |  0.25 |  0.02 |  $3.0\times 10^{-4}$ |  1 |  3
III |  $1024^{3}$ |  $0$ |  0.5 |  0.02 |  $3.0\times 10^{-4}$ |  1 |  3
Table 1: List of Runs and their relevant parameters.

As the simulations are free-decay, instead of looking for a stationary regime,
we wait for the simulations to reach a maximum of dissipation, indicating that
turbulence has had enough time to fully develop, and lead our study on times
selected around this moment. Around these times the sum of the kinetic and
magnetic energies of the system, $E_{k}+E_{m}=\rho u^{2}/2+B^{2}/(2\mu_{0})$
is expected to be decreasing steadily. Figures 1 and 2 show these dissipation
and energy for Runs I and II (Run III exhibits an almost identical behavior to
Run II). The steady decline of the energy is indeed observed for Run II, yet
Run I shows oscillations on top of the general behavior. These oscillations,
as shown in Fig. 1, match fluctuations of the internal energy and are thus
thought to be a consequence of exchanges between the kinetic plus magnetic
energy and internal energy, as was already reported in Yang et al. (2021), and
initiated by the presence of waves. Using a linear fit, one can estimate the
energy dissipation rate at the selected times: for Run I it is estimated at
$\sim-0.047$ and for Runs II and III at $\sim-0.087$. If the energy cascade
rate is indeed representative of the energy transferred to and dissipated at
small scales, its amplitude should match these estimates for the various
simulations. We will return to this point in the next section.

Figure 1: Top: Incompressible dissipation as a function of time for Run I. The
vertical dotted lines represent the 11 times selected to compute the exact
laws. Bottom: Fluctuations of kinetic, magnetic and internal energies (resp.
$E_{k}$, $E_{m}$ and $E_{int}$) as a function of time for Run I. The total
energy (green curve) shows a continuous exchange between kinetic plus magnetic
and internal energies. The narrow orange line is a linear fit of the times
studied in this paper, and its slope represents the rate of energy loss at
these times.

Figure 2: Top: Incompressible dissipation as a function of time for Run II.
The vertical dotted lines represent the 5 times selected to compute the exact
laws. Bottom: Fluctuations of kinetic plus magnetic energy as a function of
time for Run II. The narrow orange line is a linear fit of the times studied
in this paper, and its slope represents the rate of energy loss at these
times.

###  3.2 Methods of calculation

To compute all terms from equations (10) and (21) with large-enough statistics
we use two different numerical schemes depending on whether there is a
background magnetic field in the simulation or not. Both originate from the
discrete decomposition of space proposed by Taylor et al. (2003): increment
vectors $\boldsymbol{\ell}$ are selected along 73 directions defined by base
vectors connecting two points of the grid. Increments are taken as multiples
of these base vectors so that both r and
${\textbf{r}}^{\prime}={\textbf{r}}+\boldsymbol{\ell}$ lie on known grid
points, allowing for a well mapping of space without having to interpolate the
3D data, which would be a very time-consuming process. Such a decomposition
has already been successfully used to compute the two-point correlation
functions in simulation data of compressible MHD turbulence (Andrés et al.,
2018).

All source terms (i.e. terms that do not appear as divergences along the
increment vector) can be computed directly as long as the increment vector is
known. We only need to proceed to the ensemble average, which in our case is
taken on the full simulation domain:

$\left\langle\xi\right\rangle=\sum_{\textbf{r}}\frac{\xi({\textbf{r}})}{N^{3}},$
(26)

where $N=1024$ is the number of the grid points.

For flux terms (i.e. terms that appear as divergences along the increment
vector, of the form
$\boldsymbol{\nabla}_{\boldsymbol{\ell}}\cdot{\textbf{F}}$) we make the
assumption that our system is isotropic. In this situation we naturally use
spherical coordinates, in which $\boldsymbol{\ell}$ is defined as
$\boldsymbol{\ell}=(\ell,\phi,\theta)$ and the derivative operator
$\boldsymbol{\nabla}_{\boldsymbol{\ell}}$ reduces to
$\boldsymbol{\nabla}_{\boldsymbol{\ell}}\cdot\left\langle{\textbf{F}}\right\rangle=\frac{1}{\ell^{2}}\partial_{\ell}[\ell^{2}\left\langle
F_{\ell}\right\rangle(\ell)]$. Thus, for a given increment vector
$\boldsymbol{\ell}$ we only have to compute the projection of the vectorial
flux F on the direction of $\boldsymbol{\ell}$, which writes:

$\left\langle F_{\ell}\right\rangle(\ell,\phi,\theta)=\left\langle
cos(\phi)sin(\theta)F_{x}+sin(\phi)sin(\theta)F_{y}+cos(\theta)F_{z}\right\rangle.$
(27)

These projections are then averaged at fixed $\ell$:

$\left\langle F_{\ell}\right\rangle(\ell)=\sum_{\phi,\theta}\frac{\left\langle
F_{\ell}\right\rangle(\ell,\phi,\theta)}{n_{dir}},$ (28)

where $n_{dir}=73$ refers to the number of different directions taken for
$\boldsymbol{\ell}$.

Note that, the isotropy assumption stands for Runs II and III in which
$\textbf{B}_{0}=0$. A similar method based on the assumption of a symmetry of
revolution along the axis of $\textbf{B}_{0}$ was also used to study Run I.
While this method is a priori more suited for the study of simulations with
$\textbf{B}_{0}\neq 0$, the isotropic decomposition of the data ultimately
provided better results even on Run I, and is thus the only one used in this
paper. A more detailed discussion on this point is given in the Appendix.

###  3.3 Applied calculation

Using the method described above we compute the various terms of equations
(10) and (21). Note however that since the present simulations are free-decay
the forcing terms appearing in those equations are identically zero.
Therefore, the latter reduce to

$2\varepsilon_{F21}+\partial_{t}\left\langle
E^{tot}\right\rangle-\partial_{t}\left\langle
S\right\rangle+\mathcal{D}_{F21}+\mathcal{D}_{F21}^{loc}=0$ (29)

for the F21 model and to

$2\varepsilon_{A18}+\partial_{t}\left\langle
R_{E}+R_{E}^{\prime}\right\rangle+\mathcal{D}_{A18}=0.$ (30)

for the A18 model.

To lead the calculations we retained a number of snapshots for each Run: 11
snapshots evenly spaced in time for Run I (that present fluctuations on the
incompressible energy and thus require a time average) between 8 and 9
turnover times, and 5 evenly spaced snapshots for Runs II and III between 5.3
and 5.5 turnover times and 5.6 and 5.8 turnover times respectively. All time
derivatives are obtained by using a five-points finite differences method:

$\displaystyle
f^{\prime}(t)\approx\frac{f(t-2h)-8f(t-h)+8f(t+h)-f(t+2h)}{12h},$ (31)

where $h$ represents the time step between two selected snapshots. For Run I,
time derivatives (e.g. $\partial_{t}\left\langle S\right\rangle$) are
calculated on all possible subsets of 5 consecutive snapshots among the
initial 11, for a total of 7 calculations (one using times 8 to 8.4, one using
times 8.1 to 8.5 etc.), then the resulting derivatives are averaged over time.
For all other terms (e.g. $\varepsilon_{F21}$) we compute the time average
over all 11 snapshots. For Runs II and III, the derivatives are calculated on
the 5 snapshots retained using the same derivation method and all other terms
are calculated only on the central snapshot, respectively at times 5.4 and
5.7. Averaging these terms over the five snapshots was found to bring no
change to the results (not shown).

##  4 Calculation of the energy cascade rate in the inertial range

We first study the energy cascade rates obtained for laws F21 and A18
classically obtained under the full assumptions of space homogeneity, time
stationarity and infinite Reynolds number, i.e. equations (9) and (20)
respectively. For both laws, energy cascade rates are broken down into a Hall
component $\varepsilon^{Hall}$ and a MHD component $\varepsilon^{MHD}$, which
are made respectively of the terms in factor of $\lambda$ and of all the
remaining terms from equations (9) and (20). A comparison of the cascade rates
provided by the two exact laws is shown in Fig. 3 for Runs I and II. One can
observe that the two models yield closely similar components (MHD and Hall) of
the cascade rate nearly at all scales. This validates numerically the
equivalence of the two exact laws in the inertial range as anticipated in
Ferrand et al. (2021a). Note that in both Run I and Run II the value of the
cascade rate in the inertial range, centered around $\ell/d_{i}\simeq 2$,
roughly matches the energy dissipation rate estimated through the linear fit
on the energy (respectively $\sim-0.047$ and $\sim-0.087$, see Figs. 1 and 2),
suggesting that the energy cascade approximated in the inertial range is
representative of the dissipation in the system. The difference in the cascade
rate values between the two runs is due to differences in the initial/driving
amplitude of fluctuations, which is the lowest for Run I with $B_{0}=2$.

Figure 3: Comparison between the different components of the cascade rate
given by the two models of Ferrand et al. (2021a) (Eq. 9) and Andrés et al.
(2018) (Eq. 20) for Run I (top) and Run II (bottom). Plain lines represent
positive values whereas dashed lines represent negative ones.

Another question that can be addressed regarding the cascade rate in the
inertial range is its sensitivity (or not) to the turbulent sonic Mach number.
To do so we compare in Fig. 4 the results of Runs II and III that correspond
to the initial Mach numbers $M_{S}=0.25$ and $M_{S}=0.5$, respectively. We
observe that, overall, the cascade rate are very close to each other at all
scales, and appear to be similar among all three Runs, indicating that the
increase in the Mach number from $M_{S}=0.25$ to $M_{S}=0.5$ does not bring
significant changes to the total dynamics of the system. This result is in
agreement with the findings in Andrés et al. (2018) for compressible MHD
turbulence who already reported that, for Mach numbers up to $M_{S}=0.5$,
purely compressible components of the exact law remain negligible in
comparison to the flux terms, which only slightly deviate from their
incompressible counterparts. Note that such conclusions only hold a priori for
subsonic regimes: supersonic turbulent flows can develop a dominant
compressible (source-like) component of the energy cascade over the
traditional flux driven one (Ferrand et al., 2020).

Figure 4: Comparison between the cascade rate given by model Ferrand et al.
(2021a) (Eq. 9) for two different Mach numbers: $M_{S}=0.25$ (Run II) and
$M_{S}=0.5$ (Run III)

At this point, an important remark can be made: all the energy cascade rates
reported in Figs. 3 and 4 have a relatively low amplitude Hall component,
which never becomes dominant with regard to the MHD component in contrast with
results reported previously from 3D CGL simulations (Ferrand et al., 2021b)
(i.e., simulations using a closure with an anisotropic pressure tensor
introduced by Chew, Goldenberg and Low (Chew et al., 1956), hence the
acronym). An explanation can be given: the Hall effect remains too weak in our
simulations such that the dissipation inhibits the energy cascade before its
Hall component becomes dominant.

To test this hypothesis we ran an additional simulation (Run I-512) that is
akin to Run I but with a lower resolution $N=512$, a slightly higher
dissipation $\nu=\eta=8.0\times 10^{-4}$ and an ion inertial length
$d_{i}=0.05/L_{0}$. $L_{0}$ and $M_{S}$ are kept unchanged. Reducing the
resolution while increasing the value of $d_{i}$ (despite the slight increase
of the dissipation) allows for increasing the size of the sub-ion range, which
should in turn increase the importance of the Hall effect at the smallest
scales available. A first measure of this can be obtained by looking at the
power spectrum of the electric field, defined by the reduced Ohm’s law as
(omitting the resistive term):

${\bf
E}=-\textbf{v}\times\textbf{B}+\frac{1}{nq_{e}}\textbf{J}\times\textbf{B},$
(32)

where $n$ is the particle density. Fig. 5 shows the ratio of the power
spectrum density (PSD) of the Hall to the ideal components of the electric
field for Run I and the lower resolution Run I-512. As expected, we observe
both an increase in the amplitude of the Hall electric field and its extension
to smaller scales for Run I-512 in comparison with Run I.

Figure 5: Ratio of the power spectrum density of the Hall to the ideal
($\textbf{v}\times\textbf{B}$) components of the electric field for Run I and
the lower resolution Run I-512. Note that the window has been set to only show
the scales spanned by the spectrum of run I-512 ; the one for run I extends
further towards smaller values of k.

Then, we look at how the increased Hall effect witnessed in Run I-512 is
reflected on the other spectra and the energy cascade rate. Fig. 6 represents
the kinetic and magnetic power spectra of the simulations for both Run I and
Run I-512. The slope of the magnetic spectrum provides a means to pinpoint the
MHD inertial range (slope of -5/3) and the Hall range (slope of -7/3). It
appears that, in both runs, the MHD range corresponds to a wide range of
scales above the ion scale, as indicated by the relatively steady -5/3
scaling. Typically, for Run I, the -5/3 scaling roughly spans
$kd_{i}\in[0.1,0.7]$, which corresponds to $\ell/d_{i}\in[1.4,10]$ and is in
line with the results reported in Fig. 3 for the energy cascade rate. For the
Hall range however, it appears that the spectra quickly drop below the ion
scale, straying away from the -7/3 scaling. The drop is less noticeable in Run
I-512,that features a stronger Hall effect, still no clear -7/3 scaling can be
witnessed at sub-ion scales. However, a clearer knee with a change in slope
can be observed in the spectral scaling of this run at wave numbers slightly
smaller than the inverse of the ion scale. These observations suggest that, as
expected, the Hall regime enters in competition with the increasing
dissipation at small scales. Similarly, when looking at the energy cascade
rates that are reported in Fig. 7, we observe little enhancement in the
amplitude of the Hall component of Run I-512, but it still does not dominate
the cascade at sub-ion scales despite the Hall effect being stronger.

Figure 6: Kinetic and magnetic power spectra (in red and green respectively)
for run I (top) and run I-512 (bottom). Reference power laws of slopes -5/3
and -7/3 are represented by a black dashed line and a grey dash-dot line
respectively.

These results call for a further increase in the potency of the Hall effect
and/or the introduction of hyperviscosity and hyperdiffusivity in the GHOST
code, which would push the dissipation to the smallest possible scales (under
the assumption that the mechanisms responsible for energy dissipation in a
plasma act at scales much smaller than the ion inertial length, or are more
localized in spectral space than the usual viscosity and Ohmic dissipation).
The latter solution was indeed used in the CHMHD-CGL simulations where the
Hall cascade was found to dominate below the ion inertial length (Ferrand et
al., 2021b). The former however was not reasonably feasible in this study due
to computational constraints. Indeed, for the simulation code to run properly,
and as a result of the explicit time resolution of the dispersion of whistlers
in Hall MHD, one needs to set an adequate value for the time step of the
calculations. This time step has, among other constraints resulting from the
Courant-Friedrichs-Lewy condition, an upper bound that roughly behaves as
$\sim(dx)^{2}/(d_{i}B)$ where $dx$ is the grid resolution. Thus, increasing
$d_{i}$ too much while keeping the grid separation small drastically tightens
the constraint on the time step, up to a point where the computational cost
becomes prohibitive. Ultimately, at this point, hyperdissipation appears as
the best available solution.

Figure 7: Comparison between the cascade rate given by model Ferrand et al.
(2021a) (Eq. 9) for Run I and Run I-512.

##  5 Exact laws beyond the inertial range and time stationarity

As we already stated, switching from a forced turbulence model to a free-decay
model invalidates the stationarity hypothesis used to derive the final
expression of the compressible exact laws, namely equations (9)-(20). Thus, we
make use in this section of equations (29) and (30) describing the general
exact laws F21 and A18, free from these hypotheses and of the presence of an
external forcing. All the terms in those equations are computed for the three
Runs and are displayed in figures 8 and 9, along with their sums that,
according to the aforementioned equations, should amount to zero.

###  5.1 Full equation for law A18

Figure 8: Calculation of equation (30) (black line) and its components for all
three Runs. Plain lines represent positive values whereas dashed lines
represent negative ones.

Equation (30) gives a very similar behavior for all three Runs. The sum of all
terms, which is supposed to be zero, lies around 1 to 1.5 orders of magnitude
below all other terms, which is reasonable given the statistical and discrete
nature of the numerical calculations led here. Theoretical limit cases can be
easily evaluated here: at large and intermediate scales the dissipation term
should be negligible as it represents a mean of uncorrelated terms, resulting
in the equation:

$2\varepsilon_{A18}+\partial_{t}\left\langle
R_{E}+R_{E}^{\prime}\right\rangle=0.$ (33)

This relation is overall well verified by Run II and III at large scales as
reflected by the matching between the red and yellow curves in Fig. 8.
However, this does not seem to be the case in Run I and the black curve rises
at large scales as a result. The reason for this odd behavior is not fully
understood yet. At small scales the dissipation is expected to kick in and
take energy away from the cascade $\varepsilon_{A18}$. Also, at small scales
we have $\textbf{x}\to\textbf{x}^{\prime}$, and therefore
$\partial_{t}\left\langle R_{E}+R_{E}^{\prime}\right\rangle\to
2\partial_{t}\left\langle E^{tot}\right\rangle$, resulting in the equation:

$2\varepsilon_{A18}+2\partial_{t}\left\langle
E^{tot}\right\rangle+\mathcal{D}_{A18}=0.$ (34)

This equation too is verified in all three Runs, since the black curve in Fig.
8 that represents a measure of any departure from the perfect fulfillment of
the equation is at least 1 order of magnitude lower than the other components
of the equation.

The physical interpretation of Eq. (30) is rather simple and can be summed up
as follows: in free-decay simulations the term $\partial_{t}\left\langle
R_{E}+R_{E}^{\prime}\right\rangle$ in the A18 model plays the role of a
“forcing”, i.e., $-\partial_{t}\left\langle
R_{E}+R_{E}^{\prime}\right\rangle\equiv\mathcal{F}_{F21}$ that inputs energy
into the system at each time step. This reservoir of energy is then split into
a cascade component with a rate $\varepsilon_{A18}$ and a dissipation one with
a rate $\mathcal{D}_{A18}$. At large scale, since $\mathcal{D}_{A18}\to 0$ all
the energy is almost entirely injected in the cascade towards smaller scales.
The sum of all terms, which should be equal to zero, can be seen as an
estimation of the error induced in the calculations of time derivatives by the
fluctuations of energy.

###  5.2 Full equation for law F21

Figure 9: Calculation of equation (29) (black line) and its components for all
three Runs. Plain lines represent positive values whereas dashed lines
represent negative ones.

Fig. 9 shows that, again, the results of equation (29) do not change much
between the three Runs. The local terms (i.e., independent on the increment
$\boldsymbol{\ell}$) $\mathcal{D}_{F21}^{loc}$ and $\partial_{t}\left\langle
E_{tot}\right\rangle$ are almost equal, which can be intuited by the fact that
$\mathcal{D}_{F21}^{loc}\sim-\varepsilon_{F21}$ in the limit of
$\rho\sim\rho^{\prime}$. Note that in this case the model still perfectly
holds for Run I, whereas with equation (30) the term $\partial_{t}\left\langle
R_{E}+R_{E}^{\prime}\right\rangle$ showed an irregular increase at large
scale. This may suggest that the model depicted by equation (29) is more
robust than the previous one.

Interpreting the balance between the terms of the non-stationary F21 law is
more subtle than for A18 because of the presence of the local terms
$\partial_{t}\left\langle E^{tot}\right\rangle$ and $\mathcal{D}_{F21}^{loc}$.
Here it is the scale dependent term $-\partial_{t}\left\langle S\right\rangle$
that acts as a “forcing” to input energy into the system, but this energy is
diminished by the local (i.e., increment $\boldsymbol{\ell}=0$) energy
variation term $\partial_{t}\left\langle E^{tot}\right\rangle$ and the local
dissipation $\mathcal{D}_{F21}^{loc}$. Then, the resulting energy balance is
split, similarly to the A18 model, between the cascading part to small scales
and the dissipation (both being scale dependent quantities):

$2\varepsilon_{F21}+\Big{[}-\partial_{t}\left\langle
S\right\rangle+\partial_{t}\left\langle
E^{tot}\right\rangle+\mathcal{D}_{F21}^{loc}\Big{]}+\mathcal{D}_{F21}=0.$ (35)

##  6 Conclusion

Similarly to the case of IHMHD, several exact laws exist for estimating the
turbulence energy cascade rate in CHMHD. In this paper, making use of high
resolution free-decay simulations of CHMHD turbulence, we showed that the two
exact laws available for this model provide the same value of $\varepsilon$,
as was already proven for IHMHD laws (Ferrand et al., 2019). In the absence of
a direct mathematical proof of the equivalence between the two compressible
laws, this paper brings evidence that they indeed describe the same turbulent
cascade. The influence of the Hall effect on the energy cascade was also
investigated. It appears that the development of a Hall-driven energy cascade
in numerical simulations may be much more hindered by the action of
dissipation at near-ion scales than by the size of the sub-ionic range. This
underlines the potential importance of using a certain amount of
hyperviscosity when doing simulations of CHMHD as a model of the large- and
intermediate-scales in a plasma.

The question as to how the exact laws behave in absence of an external forcing
led us to investigate a more general form of compressible exact laws, dropping
the usual assumption of time stationarity and considering otherwise neglected
time derivatives. This study shows that a shift in the interpretation of both
laws occurs: instead of the continuous (in time) large scale forcing, the laws
point toward the existence of a scale dependent reservoir of energy, mainly
described by the time derivative terms, from which either the turbulent
cascade or the dissipation (depending on the considered scale) draw. This
reservoir overall coincides with the energy dissipation rate, which suggests
that the considered cascading energy is ultimately bound to be fully
dissipated.

These results confirm the non-trivial assertion that the final, well-known
form of exact laws, obtained with the assumption of stationarity, remain valid
within the inertial range even for decaying turbulent flows in which the
aforementioned assumption is not verified, as it is the case in some turbulent
space plasmas such as the SW taken far away from the sun.

## Acknowledgements

This work was granted access to the HPC resources of CINES under allocation
2021 A0090407714 made by GENCI. N.A. acknowledge financial support from
CNRS/CONICET Laboratoire International Associé (LIA) MAGNETO. N.A.
acknowledges financial support from the following grants: PICT 2018 1095 and
UBACyT 20020190200035BA.

## Appendix
axi-symmetric law calculation

Throughout this paper, all terms of the generalized exact laws are computed
using the isotropic decomposition described in section 3.2. However, for
simulations with $B_{0}\neq 0$ the isotropy assumption should not a priori be
used. Instead we thus propose an axi-symmetric scheme that may prove to be
more suited for the analysis of such simulation data.

We assume a symmetry of revolution around $\textbf{B}_{0}$, which is here
aligned with the z-axis. In this case we adopt cylindrical coordinates: the
increment vector is defined as
$\boldsymbol{\ell}=(\ell_{\perp},\phi,\ell_{z})$ and the derivative operator
becomes
$\boldsymbol{\nabla}_{\boldsymbol{\ell}}\cdot\left\langle{\textbf{F}}\right\rangle=\frac{1}{\ell_{\perp}}\partial_{\ell_{\perp}}[\ell_{\perp}\left\langle
F_{\ell_{\perp}}\right\rangle(\ell_{\perp},\ell_{z})]+\partial_{\ell_{z}}\left\langle
F_{\ell_{z}}\right\rangle(\ell_{\perp},\ell_{z})$. However, the discrete
decomposition adopted in this paper makes it impossible to effectively compute
$\partial_{\ell_{z}}\left\langle
F_{\ell_{z}}\right\rangle(\ell_{\perp},\ell_{z})$ at arbitrary values of
$(\ell_{\perp},\ell_{z})$ without resorting to multi-dimensional interpolation
on irregular grids, bringing lots of additional calculations and more
imprecision to the result. Consequently, we only consider the perpendicular
component of the flux that is averaged over the parallel increments
${\ell_{z}}$, thus the derivative operator :
$\boldsymbol{\nabla}_{\boldsymbol{\ell}}\cdot\left\langle{\textbf{F}}\right\rangle=\frac{1}{\ell_{\perp}}\partial_{\ell_{\perp}}[\ell_{\perp}\left\langle
F_{\ell_{\perp}}\right\rangle(\ell_{\perp})]$. To do so we first compute the
projection of F on the direction of $\boldsymbol{\ell}_{\perp}$:

$\left\langle
F_{\ell_{\perp}}\right\rangle(\ell_{\perp},\phi,\ell_{z})=\left\langle
cos(\phi)F_{x}+sin(\phi)F_{y}\right\rangle,$ (36)

and then take the average over all directions (i.e., over $\ell_{z}$ and
$\phi$):

$\left\langle
F_{\ell_{\perp}}\right\rangle(\ell_{\perp})=\sum_{\phi,\ell_{z}}\frac{\left\langle
F_{\ell_{\perp}}\right\rangle(\ell_{\perp},\phi,\ell_{z})}{57}.$ (37)

Here we only use 57 directions, corresponding to all directions of the
isotropic decomposition forming an angle of $45^{\circ}$ or more with the
parallel direction. This is done to obtain results pertaining to all values of
$\ell_{\perp}$ while avoiding redundant calculations due to the periodicity of
the data cubes.

This method presents two inconveniences: first, the statistics are weaker than
the isotropic decomposition as we probe a smaller number of directions for the
increment vector. This may lead to less precise calculations. Second,
disregarding parallel fluxes may lead to miss a small portion of the cascade.
Due to these limitations the isotropic decomposition may sometimes yield
slightly better results even in presence of a background magnetic field. A
good way to test the efficiency of this method is to compute the energy
cascade rate of both exact laws A18 and F21 and see if they match as one would
expect. This test was used on Run I and the results are reported in Fig. 10.
We immediately observe that the two laws, and especially the Hall components,
do not match as well as they are with the isotropic decomposition or for Run
II (see Fig. 3), which may be a consequence of the aforementioned limitations.
For this reason, we finally chose to keep using the isotropic decomposition to
study Run I in this paper.

Figure 10: Energy cascade rate calculated with laws A18 and F21 for Run I,
using the axi-symmetric method.

Still, note that Run I remains a weakly anisotropic simulation, with a
background magnetic field of only 2. Thus, the conclusion drawn in this
appendix may not be true for simulations featuring a stronger background
field, i.e. stronger anisotropies. Ultimately, checking which method is the
most suited one to lead the study should be done on a case-by-case basis.

## References

  * Andrés et al. (2018) Andrés, N., Galtier, S., & Sahraoui, F. 2018, Phys. Rev. E, 97, 013204
  * Andrés et al. (2020) Andrés, N., Romanelli, N., Hadid, L. Z., et al. 2020, The Astrophysical Journal, 902, 134
  * Andrés & Sahraoui (2017) Andrés, N., & Sahraoui, F. 2017, Phys. Rev. E, 96, 053205
  * Andrés et al. (2018) Andrés, N., Sahraoui, F., Galtier, S., et al. 2018, J. Plasma Physics, 84, 1
  * Andrés et al. (2019) Andrés, N., Sahraoui, F., Galtier, S., et al. 2019, Phys. Rev. Lett., 123, 245101
  * Arzoumanian et al. (2011) Arzoumanian, D., André, Ph., Didelon, P., et al. 2011, A&A, 529, L6
  * Bandyopadhyay et al. (2020) Bandyopadhyay, R., Sorriso-Valvo, L., Chasapis, A., et al. 2020, Phys. Rev. Lett., 124, 225101
  * Banerjee & Galtier (2013) Banerjee, S., & Galtier, S. 2013, Phys. Rev. E, 87, 013019
  * Banerjee & Galtier (2014) Banerjee, S., & Galtier, S. 2014, J. Fluid Mech., 742, 230
  * Banerjee & Galtier (2016) Banerjee, S., & Galtier, S. 2016, Phys. Rev. E, 93, 033120
  * Banerjee & Galtier (2017) —. 2017, J. Phys. A: Math. Theo., 50, 015501
  * Banerjee et al. (2016) Banerjee, S., Hadid, L., Sahraoui, F., & Galtier, S. 2016, Astrophys. J. Lett., 829, L27
  * Banerjee & Kritsuk (2017) Banerjee, S., & Kritsuk, A. 2017, Phys. Rev. E, 96, 053116
  * Banerjee & Kritsuk (2018) —. 2018, Phys. Rev. E, 97, 023107
  * Bruno & Carbone (2013) Bruno, R., & Carbone, V. 2013, Living Reviews in Solar Physics, 10, 2
  * Carbone et al. (2009) Carbone, V., Marino, R., Sorriso-Valvo, L., Noullez, A., & Bruno, R. 2009, Phys. Rev. Lett., 103, 061102
  * Carbone et al. (2009) Carbone, V., Sorriso-Valvo, L., & Marino, R. 2009, EPL (Europhysics Letters), 88, 25001
  * Chew et al. (1956) Chew, G. F., Goldberger, M. L., & Low, F. E. 1956, Proc. R. Soc. London. A, 236, 112
  * Coburn et al. (2015) Coburn, J., Forman, M., Smith, C., Vasquez, B., & Stawarz, J. 2015, Phil. Trans. R. Soc. A, 373, 20140150
  * Ferrand et al. (2021a) Ferrand, R., Galtier, S., & Sahraoui, F. 2021a, Journal of Plasma Physics, 87, 905870220
  * Ferrand et al. (2020) Ferrand, R., Galtier, S., Sahraoui, F., & Federrath, C. 2020, Astrophys. J., 904, 160
  * Ferrand et al. (2019) Ferrand, R., Galtier, S., Sahraoui, F., et al. 2019, Astrophys. J., 881, 50
  * Ferrand et al. (2021b) Ferrand, R., Sahraoui, F., Laveder, D., et al. 2021b, The Astrophysical Journal, 923, 122
  * Galtier (2008) Galtier, S. 2008, Phys. Rev. E, 77, 015302
  * Galtier (2012) —. 2012, Astrophys. J., 746, 184
  * Galtier & Banerjee (2011) Galtier, S., & Banerjee, S. 2011, Phys. Rev. Lett., 107, 134501
  * Gómez et al. (2005) Gómez, D. O., Mininni, P. D., & Dmitruk, P. 2005, Physica Scripta, 123
  * Hadid et al. (2017) Hadid, L., Sahraoui, F., & Galtier, S. 2017, Astrophys. J., 838, 11
  * Hadid et al. (2018) Hadid, L. Z., Sahraoui, F., Galtier, S., & Huang, S. Y. 2018, Phys. Rev. Lett., 120, 055102
  * Hellinger et al. (2018) Hellinger, P., Verdini, A., Landi, S., Franci, L., & Matteini, L. 2018, Astrophys. J., 857, L19
  * Kolmogorov (1941) Kolmogorov, A. 1941, Dokl. Akad. Nauk SSSR, 32, 16
  * Kritsuk et al. (2007) Kritsuk, A. G., Norman, M. L., Padoan, P., & Wagner, R. 2007, Astrophys. J., 665, 416
  * Kritsuk et al. (2013) Kritsuk, A. G., Wagner, R., & Norman, M. L. 2013, J. Fluid Mech., 729, 1
  * Lindborg (2019) Lindborg, E. 2019, JFM, 874, R2
  * MacBride et al. (2008) MacBride, B., Smith, C., & Forman, M. 2008, Astrophys. J., 679, 1644
  * Marino et al. (2008) Marino, R., Sorriso-Valvo, L., Carbone, V., et al. 2008, Astrophys. J., 677, L71
  * Marino et al. (2011) Marino, R., Sorriso-Valvo, L., Carbone, V., et al. 2011, Planetary and Space Science, 59, 592, cross-Scale Coupling in Plasmas
  * Matthaeus et al. (1999) Matthaeus, W. H., Zank, G. P., Smith, C. W., & Oughton, S. 1999, Phys. Rev. Lett., 82, 3444
  * Mininni & Pouquet (2009) Mininni, P. D., & Pouquet, A. 2009, Phys. Rev. E, 80, 025401
  * Mininni et al. (2011) Mininni, P. D., Rosenberg, D., Reddy, R., & Pouquet, A. 2011, Parallel Computing, 37, 316
  * Osman et al. (2011) Osman, K. T., Wan, M., Matthaeus, W. H., Weygand, J. M., & Dasso, S. 2011, Phys. Rev. Lett., 107, 165001
  * Politano & Pouquet (1998) Politano, H., & Pouquet, A. 1998, Phys. Rev. E, 57, 21
  * Pouquet & Patterson (1978) Pouquet, A., & Patterson, G. S. 1978, Journal of Fluid Mechanics, 85, 305–323
  * Sahraoui et al. (2020) Sahraoui, F., Zafer Hadid, L., & Huang, S. 2020, Reviews of Modern Plasma Physics, 4, doi:10.1007/s41614-020-0040-2
  * Simon & Sahraoui (2021) Simon, P., & Sahraoui, F. 2021, General exact law of compressible isentropic magnetohydrodynamic flows: theory and spacecraft observations in the solar wind, , , arXiv:2105.08011
  * Sorriso-Valvo et al. (2007) Sorriso-Valvo, L., Marino, R., Carbone, V., et al. 2007, Phys. Rev. Lett., 99, 115001
  * Sorriso-Valvo et al. (2019) Sorriso-Valvo, L., Catapano, F., Retinò, A., et al. 2019, Phys. Rev. Lett., 122, 035102
  * Stawarz et al. (2009) Stawarz, J. E., Smith, C. W., Vasquez, B. J., Forman, M. A., & MacBride, B. T. 2009, Astrophys. J., 697, 1119
  * Taylor et al. (2003) Taylor, M., Kurien, S., & Eyink, G. 2003, Phys. Rev. E, 68, 026310
  * Verdini et al. (2015) Verdini, A., Grappin, R., Hellinger, P., Landi, S., & Müller, W. C. 2015, The Astrophysical Journal, 804, 119
  * Yang et al. (2021) Yang, Y., Wan, M., Matthaeus, W. H., & Chen, S. 2021, Journal of Fluid Mechanics, 916, A4