# Efficient kinetic Lattice Boltzmann simulation of three-dimensional Hall-MHD
Turbulence

Raffaello Foldes1,2 Emmanuel Lévêque2<EMAIL_ADDRESS>Raffaele
Marino2  Ermanno Pietropaolo1  Alessandro De Rosis3  Daniele Telloni4  Fabio
Feraco2  1Dipartimento di Scienze Fisiche e Chimiche, Università dell’Aquila,
67100 Coppito (AQ), Italy 2Univ Lyon, CNRS, École Centrale de Lyon, INSA Lyon,
Univ Claude Bernard Lyon I, LMFA UMR 5509, F-69134 Ecully cedex, France
3Department of Mechanical, Aerospace and Civil Engineering, The University of
Manchester, Manchester M13 9PL, United Kingdom 4National Institute for
Astrophysics - Astrophysical Observatory of Torino, Via Osservatorio 20,
I-10025 Pino Torinese, Italy

###### Abstract

Simulating plasmas in the Hall-MagnetoHydroDynamics (Hall-MHD) regime, taking
advantage of modern high-performance computing facilities, represents a
valuable approach for the investigation of complex non-linear dynamics
developing in astrophysical frameworks and fusion machines. Taking into
account the Hall electric field is computationally very challenging as it
involves the integration of an additional term, proportional to
$\boldsymbol{\nabla}\times((\boldsymbol{\nabla}\times\mathbf{B})\times\mathbf{B})$
in the Faraday’s induction law. The latter feeds back on the magnetic field
$\mathbf{B}$ at small scales (between the ion and electron inertial scales),
requiring very high resolutions in both space and time in order to properly
describe its dynamics. The computational advantage provided by the kinetic
Lattice Boltzmann (LB) approach is exploited here to develop a new code, the
Fast Lattice-Boltzmann Algorithm for Mhd Experiments (flame). A major
advantage of dealing with a kinetic representation is that the derivatives of
the magnetic field are directly encompassed in the solution, thus allowing for
an immediate and intrinsically accurate implementation of the Hall electric
field. The flame code integrates the plasma dynamics in lattice units coupling
two kinetic schemes, one for the fluid protons (including the Lorentz force),
the other to solve the induction equation describing the evolution of the
magnetic field. The incompressibility condition is attained in the low-Mach-
number limit, whereas the divergence-free condition for the magnetic field is
preserved within machine precision. Here, the newly developed algorithm is
tested against an analytical wave-solution of the dissipative Hall-MHD
equations, pointing out its stability and second-order convergence, over a
wide range of the control parameters. Then, flame was used to reproduce
turbulent dynamics for different values of the Hall parameter. Spectral
properties of the simulated plasma are finally compared with those obtained
from numerical solutions from the well-established pseudo-spectral code,
ghost, in the MHD regime. Furthermore, the LB simulations we present, varying
the Hall parameter, highlight the transition from the MHD to the Hall-MHD
regime at large wavenumbers, exhibiting scaling exponents of the magnetic
energy spectrum at sub-ion scales in excellent agreement with that measured in
solar wind plasmas. Finally, we emphasized how our code scales efficiently
over many-core Graphic Processing Units accelerators, demonstrating it is
suited for the investigation of Hall-MHD turbulence in astrophysical and
laboratory plasmas, requiring very high resolution in order to reproduce their
dynamics.

    For the purpose of Open Access, a CC-BY public copyright licence has been
    applied by the authors to the present document and will be applied to all
    subsequent versions up to the Author Accepted Manuscript arising from this
    submission

##  1 Introduction

In the frame of the MHD model plasma is treated as a single species quasi-
neutral fluid with conductive properties sensitive to the action of the
magnetic field (Galtier, 2016). In the ideal MHD description, ions and
electrons are tied to the magnetic field, moving with the same velocity. The
Hall-MHD model relaxes the MHD prescriptions assuming ions disunite from the
magnetic field due to their inertia, while electrons remain bound to it
(Pandey & Wardle, 2008). In this framework the resistive Ohm’s law is
generalized through the introduction of the Hall electric field, proportional
to $\mathbf{J}\times\mathbf{B}$, where $\mathbf{J}$ and $\mathbf{B}$ denote
the current density and the magnetic field, respectively. The Hall electric
field has an effect on the magnetic field at length scales shorter than the
ion inertial length $d_{i}=c/\omega_{pi}$ ($\omega_{pi}$ being the ion plasma
frequency, $c$ the speed of light) as well as at time scales shorter than the
ion cyclotron period $1/\omega_{ci}$ (Huba, 2003). The scale $d_{i}$
corresponds to the scale at which ions and electrons decouple, and the
magnetic field becomes frozen into the electron fluid rather than in the bulk
plasma. Hall-MHD has been already adopted in literature to describe a variety
of astrophysical, space and laboratory environments, and to provide a detailed
description of plasma dynamics. Its applications spans from the star formation
(Norman & Heyvaerts, 1985; Marchand, P. et al., 2018) to the solar atmosphere
and the solar wind (Galtier & Buchlin, 2007; González-Morales et al., 2019),
and it has bee used also to investigate magnetic reconnection processes (Wang
et al., 2001; Morales et al., 2005; Ma et al., 2018) and the dynamo action
(Mininni et al., 2002, 2005; Gómez et al., 2010). A major difficulty in
simulating Hall-MHD is related to the need for resolving whistler waves,
evolving on fast dynamics with a phase speed $c_{w}(k)\propto k$ increasing
linearly with the wavenumber $k$. In order to properly capture the propagation
of the perturbations caused by the Hall effect, it is therefore necessary to
capture those plasma waves with $\max(c_{w})\propto 1/\Delta x$, at the
smallest resolved wavelength $\Delta x$. The Courant-Friedrichs-Lewy (CFL)
condition then yields ${\Delta t}\propto\Delta x^{2}$. This scaling implies a
rapid decrease of the time-step as the spatial resolution increases, which
poses severe limitations implying a high computational cost. Nevertheless,
Hall-MHD simulations have been proposed over the years in numerous studies,
through the integration of the equations with pseudo-spectral (Mininni et al.,
2003), finite-volume (Tóth et al., 2008; Marchand, P. et al., 2018) or hybrid
particle-in-cell codes (Ma et al., 2018; Papini et al., 2019). When dealing
with turbulent flows, pseudo-spectral methods are usually recognized as the
best option that allows for an equally-accurate representation of the fields
at the resolved dynamical scales (Patterson & Orszag, 1971). On the other
hand, their computational cost can be prohibitive when it comes to the
integration of simulations in three dimensions and for may turnover times
(Huba, 2003). Main purpose of the novel code that we developed here, flame
(Fast Lattice-Boltzmann Algorithm for MHD Experiments), is to overcome this
issue. Indeed, the Lattice Boltzmann (LB) implementation provides an
alternative to achieve a convenient trade-off between accuracy and
computational efficiency. Unlike more traditional methods that solve the
dynamics of flows at the macroscopic level, LB methods operate at an
underlying mesoscopic kinetic level. The flow complexity emerges from re-
iterating simple rules of collision and streaming of populations of particles
moving along the links of a regular cubic lattice (Krueger et al., 2016). The
connection between such an idealized representation and the macroscopic
dynamics is by now well-established and accepted, placing the method on a
solid theoretical and mathematical ground (Shan & He, 1998). Furthermore, due
to its intrinsically discrete nature and its focus on the local dynamics, it
is also computationally extremely efficient (Körner et al., 2006). A decisive
contribution to make possible the simulation of ideal MHD plasmas by means of
LB methods was made by Dellar (2002), who showed that the native LB framework
based on the Bhatnagar-Gross-Krook (BGK) collision (Bhatnagar et al., 1954)
could be consistently extended to encompass both the fluid dynamics driven by
the Lorentz force and the induction equation for the magnetic field. The
scheme introduced by Dellar overcomes the major limitations of previous
efforts (Montgomery & Doolen, 1987; Chen et al., 1991; Succi et al., 1991;
Martínez et al., 1994) and fully complies with the macroscopic MHD equations
in a weakly-compressible formulation (see §3). Nevertheless, it is prone to
develop numerical instabilities when strong gradients emerge in the flow, thus
delaying in the community its implementation for the simulation of turbulent
fluid frameworks. More recently, De Rosis et al. (2018) have shown that it is
possible to prevent the flow to become unstable by replacing the original BGK
collision with a collision operator defined in the space of moments, allowing
the explicit damping of the _non-hydrodynamic modes_ responsible for the onset
of spurious instabilities. An entropic stabilization has also been proposed by
Flint & Vahala (2018), though leading to a more complicated scheme. These
advances encouraged us to pursue with the LB modeling to simulate Hall-MHD
turbulence, an effort that has never been undertaken previously. The
development of flame was also strongly motivated by the need of the community
for innovative numerical tools for the study of space plasma turbulent
dynamics at scales that are by now within the reach of high-resolution
instruments on board spacecrafts, such as the the ESA mission Solar Orbiter
(Müller et al., 2020a). The paper is organized as follows. In §2, the Hall-MHD
equations are presented in a form that is relevant for LB developments. The LB
scheme implemented in flame is introduced and discussed in §3. The coupling
between the fluid and the magnetic lattices is explained, as well as the
inclusion of the Hall effect in the collision operator. The conversion from
physical to lattice units is discussed in great detail. §4 is devoted to the
validation of the code against an analytical solution of the dissipative Hall-
MHD equations (Xia & Yang, 2015). This section provides an assessment of the
numerical stability and a quantitative estimation of the dispersion and
dissipation errors. The computational efficiency is discussed in §5, where
GPU-accelerated simulations of the three-dimensional Orszag-Tang vortex
problem are discussed (Orszag & Tang, 1979). In a regime of high Reynolds
numbers, we show that LB simulations are able to reproduce the break in the
magnetic energy spectrum at sub-ion scales, in perfect agreement with solar-
wind measurements. Finally, in §7 we draw conclusions and outline potential
applications for the investigation of geophysical fluids and space plasmas.

##  2 The Hall-MHD equations

In this section, the Hall-MHD equations are introduced in the standard
incompressible approximation and in a weakly-compressible formulation,
suitable for LB developments.

###  2.1 Incompressible formulation

In this context, when we refer to the macroscopic description of the plasma
what we mean is the description of the prognostic fields appearing in the
model equations. Thus, at the macroscopic level, the incompressible resistive
MHD equations for an electrically conductive quasi-neutral fluid consist of
the incompressible Navier-Stokes equations, with the addition of the Lorentz
force, coupled with the resistive induction equation for the magnetic field,

$\boldsymbol{\nabla}\cdot\mathbf{U}=0$ (1)
$\partial_{t}\mathbf{U}+(\mathbf{U}\cdot\boldsymbol{\nabla})\mathbf{U}=\frac{1}{\rho_{0}}\mathbf{J}\times\mathbf{B}-\frac{1}{\rho_{0}}\boldsymbol{\nabla}\mathrm{p}+\nu\boldsymbol{\nabla}^{2}\mathbf{U}$
(2)
$\partial_{t}\mathbf{B}=\boldsymbol{\nabla}\times(\mathbf{U}\times\mathbf{B}-\eta\boldsymbol{\nabla}\times\mathbf{B})$
(3) $\boldsymbol{\nabla}\cdot\mathbf{B}=0$ (4)

where $\rho_{0}$ is the mass density of the fluid, $\nu$ is the kinematic
viscosity and $\eta$ it the magnetic resistivity. The electric current density
is expressed as $\mathbf{J}={1}/{\mu_{0}}\boldsymbol{\nabla}\times\mathbf{B}$,
where $\mu_{0}$ is the magnetic permeability in the vacuum. To account for the
Hall effect, it is necessary to take a step back in the mathematical
developments and resort to a two-fluid description that includes the fluid
equations for both ions (_i_) and electrons (_e_) separately. For a fully
ionized plasma in which the masses of protons and electrons (hereafter $i$ and
$e$) are $m_{e}\ll m_{i}\approx m$, the momentum equations read as:

$\rho[\partial_{t}\mathbf{U}+(\mathbf{U}\cdot\boldsymbol{\nabla})\mathbf{U}]=en(\mathbf{E}+\mathbf{U}\times\mathbf{B})-\boldsymbol{\nabla}\mathrm{p}_{i}+\boldsymbol{\nabla}\cdot\mathbf{\sigma}+\mathbf{R}$
(5)
$0=-en(\mathbf{E}+\mathbf{U_{e}}\times\mathbf{B})-\boldsymbol{\nabla}\mathrm{p}_{e}-\mathbf{R}$
(6)

where $\mathbf{\sigma}$ is the viscous stress tensor, $n$ the particle density
($\rho=mn$) and $\mathbf{R}$ the rate (per unit volume) of momentum exchange
due to collisions between protons and electrons. The latter is given by
$\mathbf{R}=-mnf_{ie}(\mathbf{U}-\mathbf{U_{e}})$ where $f_{ie}$ denotes the
collision frequency and can be reformulated as $\mathbf{R}=-(mf_{ie}/e)\
\mathbf{J}$, with the density current
$\mathbf{J}=en(\mathbf{U}-\mathbf{U_{e}})$. By summing (5) and (6), assuming
$\sigma_{ij}=\rho\nu(\partial_{i}U_{j}+\partial_{j}U_{i})$, one obtains

$\partial_{t}\mathbf{U}+(\mathbf{U}\cdot\boldsymbol{\nabla})\mathbf{U}=\frac{1}{\rho}\mathbf{J}\times\mathbf{B}-\frac{1}{\rho}\boldsymbol{\nabla}\mathrm{p}+\nu\boldsymbol{\nabla}^{2}\mathbf{U}.$
(7)

On the other hand, by replacing $\mathbf{U_{e}}$ by $\mathbf{U}-\mathbf{J}/ne$
and the expression for the rate of momentum exchange into (6), the Ohm’s law
becomes

$\mathbf{E}=-(\mathbf{U}-\frac{1}{en}\mathbf{J})\times\mathbf{B}+\frac{1}{en}\boldsymbol{\nabla}\mathrm{p}_{e}+\frac{mf_{ie}}{e^{2}n}\mathbf{J}.$
(8)

Taking the curl of this equation gives in the end an induction equation with
the Hall’s current correction in standard physical units as

$\partial_{t}\mathbf{B}=\boldsymbol{\nabla}\times[(\mathbf{U}-\alpha_{H}\mathbf{J})\times\mathbf{B}]+\eta\boldsymbol{\nabla}^{2}\mathbf{B}$
(9)

where $\alpha_{H}=1/en$ is usually referred to as the Hall parameter and the
magnetic resistivity $\eta=mf_{ie}/(e^{2}n\mu_{0})$.

###  2.2 Weakly-compressible formulation

Incompressibility is an assumption made at the macroscopic level and cannot be
implemented in the mesoscopic representation as this would imply that the
particles move at infinite speed, in order to adapt instantaneously the
pressure. Incompressibility can nevertheless be approached in the so-called
_weakly-compressible_ limit, in which the speed of sound waves $c_{s}$ becomes
much larger than the typical fluid velocity $U_{0}$, or equivalently, the
pressure field adapts in a time shorter than the time over which the flow
evolves. This regime is attained for vanishing Mach number, $\mathrm{Ma}\equiv
U_{0}/c_{s}$. Consequently, the incompressible equations should be replaced
with the compressible formulation

$\partial_{t}\rho+\boldsymbol{\nabla}\cdot(\mathbf{\rho\mathbf{U}})=0$ (10)
$\partial_{t}(\mathbf{\rho
U})+\boldsymbol{\nabla}\cdot\bigg{(}\rho\mathbf{U}\otimes\mathbf{U}+\mathrm{p}\mathbb{I}+\frac{1}{2}|\mathbf{B}|^{2}\mathbb{I}-\mathbf{B}\otimes\mathbf{B}\bigg{)}=\rho\nu\boldsymbol{\nabla}^{2}\mathbf{U}$
(11)

in which the Lorentz force has been rewritten in a conservative form, since
the divergence of the Maxwell stress tensor
$M_{\alpha\beta}=\frac{1}{2}|\mathbf{B}|^{2}\delta_{\alpha\beta}-B_{\alpha}B_{\beta}$111The
notation $\mathbf{a}\otimes\mathbf{b}\equiv a_{\alpha}b_{\beta}$ is adopted
and $\mu_{0}$ has been absorbed by replacing $\mathbf{B}$ with
$\mu_{0}^{1/2}\mathbf{B}$. This (standard) normalization will be assumed
hereafter, that allows to simplify the Lorentz force as
$(\boldsymbol{\nabla}\times\mathbf{B})\times\mathbf{B}$. Compressibility
requires resorting to an equation of state linking pressure, mass density and
temperature. Here, the low-Mach limit justifies the use of a simple isothermal
relation

$\mathrm{p}=\rho c_{s}^{2}$ (12)

which is consistent with $O(\mathrm{Ma}^{2})$ mass-density fluctuations. The
induction equation describing the evolution of the magnetic field can be
rewritten in the same fashion as

$\partial_{t}\mathbf{B}+\boldsymbol{\nabla}\cdot\left((\mathbf{U}-\alpha_{H}\mathbf{J})\otimes\mathbf{B}-\mathbf{B}\otimes(\mathbf{U}-\alpha_{H}\mathbf{J})\right)=\eta\boldsymbol{\nabla}^{2}\mathbf{B}.$
(13)

Let us remark here that following the normalization of $\mathbf{B}$ by
$\mu_{0}^{1/2}$, the Hall current $\alpha_{H}\mathbf{J}$ reads as
${\alpha_{H}}/{{\mu_{0}}^{1/2}}~{}\boldsymbol{\nabla}\times\mathbf{B}$. In the
next sections, the developed LB scheme will conform to the set of equations
(10), (11), (12) and (13). The divergence-free condition on $\mathbf{B}$ is
preserved by (13), justifying that it is sufficient to impose
$\boldsymbol{\nabla}\cdot\mathbf{B}=0$ initially. In the numerical modelling,
a particular attention will be paid to verify that this condition is indeed
preserved with accuracy.

##  3 Hall-MHD Lattice Boltzmann scheme

In this section, the standard LB method for classical fluid dynamics is
briefly introduced, focusing on key steps, then it is extended to encompass
Hall-MHD. Further details are provided in the appendix A. A central-moment
collision operator (De Rosis et al., 2018) and a high-connectivity D3Q27
lattice are used to integrate the fluid protons, while the evolution of the
magnetic field is accounted by a Bhatnagar-Gross-Krook (BGK) collision
operator (Bhatnagar et al., 1954) and a low-connectivity D3Q7 lattice. Our
original contribution to these developments is the self-consistent integration
of the Hall term in the LB scheme by suitably redefining the equilibrium state
for the magnetic field.

###  3.1 Lattice Boltzmann scheme for the fluid dynamics

####  3.1.1 Standard BGK Lattice Boltzmann scheme

The LB method (Krueger et al., 2016) is based on the idea that fluid motions
can be represented by the collective behavior of fictitious (introduced in the
frame of the LB integration strategy) particle populations evolving along the
links of a cubic lattice. When the lattice connectivity, which accounts for
the discrete directions of propagation of the particles, is high enough to
satisfies sufficient isotropy, weakly-compressible Navier-Stokes dynamics can
be simulated with an $O(\mathrm{Ma}^{3})$ error. The macroscopic variables
such as the fluid density, momentum, or stress tensor are obtained as
statistical moments of the particle distributions, _i.e._

$\displaystyle\rho$ $\displaystyle=$ $\displaystyle\sum_{i=0}^{N-1}f_{i}$ (14)
$\displaystyle\rho\mathbf{U}$ $\displaystyle=$
$\displaystyle\sum_{i=0}^{N-1}f_{i}\mathbf{c}_{i}$ (15)
$\displaystyle\boldsymbol{\Sigma}$ $\displaystyle=$
$\displaystyle\sum_{i=0}^{N-1}f_{i}\mathbf{c}_{i}\otimes\mathbf{c}_{i}$ (16)

by summing over the local mass densities $f_{0},\cdots,f_{N-1}$ of particles
moving with velocities $\boldsymbol{c}_{0},\cdots,\boldsymbol{c}_{N-1}$,
respectively. The sums replace here the integrals over $\mathbf{c}$ of the
classical kinetic theory as the result of a drastic decimation in velocity of
the phase space. From a theoretical viewpoint, the LB method is derived by
expanding the solution of the continuum Boltzmann equation onto a finite basis
of Hermite polynomials in velocity, and by resorting to a Gaussian quadrature
formula to express the statistical moments (He & Luo, 1997). As a consequence,
the particle densities ${f_{i}}(\mathbf{x},t)$ evolve according to a discrete-
velocity analogue of the Boltzmann equation, which reads as

$\partial_{t}f_{i}+\left(\mathbf{c}_{i}\cdot\boldsymbol{\nabla}\right)f_{i}=-\frac{1}{\tau}\bigg{(}f_{i}-f_{i}^{(0)}(\rho,\mathbf{U})\bigg{)}$
(17)

under the BGK approximation (Bhatnagar et al., 1954). The latter assumes that
collisions are responsible for the relaxation of the particle densities
towards their equilibrium state $f_{i}^{(0)}(\rho,\mathbf{U})$, with a unique
relaxation time $\tau=\nu/c_{s}^{2}$.

The _Lattice_ keyword refers to the discretization in space and time of (17)
with a set of microscopic velocity
$\boldsymbol{c}_{0},~{}\cdots,~{}\boldsymbol{c}_{N-1}$ chosen in a way such
that particles travel from a lattice node to a neighbour lattice node in
exactly one time-step (see Fig. 1).

Figure 1: Typical cubic lattices with 15, 19 and 27 velocities, respectively.
At each lattice node, the microscopic velocities point towards the centre
(black), the 6 centres of faces (green), the 12 centres of the edges (red) or
the 8 corners (blue) of a cube. The arrows represent the local displacements
$\boldsymbol{c}_{i}\Delta t$ of particles from a lattice node to a
neighbouring node during exactly one time-step. In the present study, a D3Q27
lattice that is more appropriate to simulate strongly non-linear fluid
dynamics is considered (Silva & Semiao, 2014).

The LB scheme then expresses simply as

$\bar{f}_{i}(\mathbf{x}+\mathbf{c}_{i}\Delta t,t+\Delta
t)=\bar{f}_{i}(\mathbf{x},t)-\omega\bigg{(}\bar{f}_{i}(\mathbf{x},t)-{f}^{(0)}_{i}(\rho,\mathbf{U})(\mathbf{x},t)\bigg{)}$
(18)

by considering the change of variables $\bar{f}_{i}=f_{i}+{\Delta
t}/{2\tau}(f_{i}-f_{i}^{(0)}).$ This change of variable comes from the
trapezoidal rule used to approximate the integral of the collision operator
(right-hand side of (17)) between $t$ and $t+\Delta t$ (Krueger et al., 2016).
It also calls for a redefinition of the relaxation time as $\tau+\Delta t/2$
so that

$\frac{1}{\omega}=\bigg{(}\frac{\nu}{c_{s}^{2}\Delta t}+\frac{1}{2}\bigg{)}$
(19)

where the speed of sound $c_{s}$ is linked to the lattice spacing by $\Delta
x/\Delta t=\sqrt{3}c_{s}$ for the D3Q27 lattice. The expressions of the mass
density and fluid momentum as statistical moments remain unchanged with

$\displaystyle\rho$ $\displaystyle=$
$\displaystyle\sum_{i=0}^{N-1}\bar{f}_{i}$ (20)
$\displaystyle\mathrm{and}\quad\rho\mathbf{U}$ $\displaystyle=$
$\displaystyle\sum_{i=0}^{N-1}\bar{f}_{i}\mathbf{c}_{i}.$ (21)

In practice, (18) is divided into a two-step algorithm with a streaming step
consecutive to a local collision operation, _i.e._

$\displaystyle\bar{f}_{i}(\mathbf{x}+\mathbf{c}_{i}\Delta t,t+\Delta t)$
$\displaystyle=$ $\displaystyle\bar{f}_{i}^{*}(\mathbf{x},t)$ (22)
$\displaystyle\bar{f}^{*}_{i}(\mathbf{x},t)$ $\displaystyle=$
$\displaystyle\bar{f}_{i}(\mathbf{x},t)-\omega\left(\bar{f}_{i}(\mathbf{x},t)-{f}^{(0)}_{i}(\rho,\mathbf{U})(\mathbf{x},t)\right).$
(23)

To complete the algorithm, the particle densities at the equilibrium
${f}_{i}^{(0)}$ need to be specified. By construction, ${f}_{i}^{(0)}$ is
defined as a truncated Hermite expansion of the continuous Maxwell-Boltzmann
distribution evaluated in $\mathbf{c_{i}}$, which reads as

$f_{i}^{(0)}(\rho,\mathbf{U})=w_{i}\rho\left(1+\frac{\boldsymbol{c}_{i}\cdot\mathbf{U}}{c_{s}^{2}}+\frac{(\boldsymbol{c}_{i}\cdot\mathbf{U})^{2}}{2c_{s}^{4}}-\frac{\mathbf{U}\cdot\mathbf{U}}{2c_{s}^{2}}+\cdots\right)$
(24)

with the weights $w_{\mathrm{center}}=8/27$, $w_{\mathrm{face}}=2/27$,
$w_{\mathrm{edge}}=1/54$ and $w_{\mathrm{corner}}=1/216$ for the D3Q27
lattice. An expansion truncated at the second order is enough to recover the
Navier-Stokes equations with an $O(\mathrm{Ma}^{3})$ error. However, several
groups (Malaspinas, 2015; Coreixas et al., 2017, 2019; De Rosis & Luo, 2019)
have recently shown that accounting for high-order terms results in a gain in
accuracy and stability. In our code, ${f}_{i}^{(0)}$ has been developed up to
the sixth order. The extension of the standard LB algorithm to encompass the
Lorentz force is straightforward and relies on the fundamental property that
the second-order statistical moment at equilibrium gives the conservative part
of the stress tensor. Therefore, incorporating the Lorentz force in the
equation describing the fluid dynamics, or equivalently, the Maxwell tensor in
the stress tensor amounts to upgrading the equilibrium state as

${f}_{i}^{\mathrm{mhd}(0)}(\rho,\mathbf{U},\mathbf{B})={f}_{i}^{(0)}(\rho,\mathbf{U})+\frac{w_{i}}{2c_{s}^{4}}\left((\mathbf{B}\cdot\mathbf{B})(\boldsymbol{c}_{i}\cdot\boldsymbol{c}_{i})-(\mathbf{c}_{i}\cdot\mathbf{B})^{2}\right)$
(25)

in order to obtain

$\boldsymbol{\Sigma}^{\mathrm{mhd}(0)}=\sum_{i=0}^{N-1}f_{i}^{\mathrm{mhd}(0)}\mathbf{c}_{i}\otimes\mathbf{c}_{i}=\rho\mathbf{U}\otimes\mathbf{U}+\mathrm{p}\mathbb{I}+\frac{1}{2}|\mathbf{B}|^{2}\mathbb{I}-\mathbf{B}\otimes\mathbf{B}.$
(26)

This concludes the introduction of the standard BGK-LB algorithm for MHD.

####  3.1.2 Central-moment Lattice Boltzmann scheme

Despite its simplicity, effectiveness and large popularity, the BGK-LB scheme
is known to suffer from numerical instability when large velocity gradients
develop in the flow. This issue made it necessary to adapt either the
numerical discretization of (17) or the collision operator (Krueger et al.,
2016). If the former leads to more stable schemes, accuracy is also
considerably degraded. This drawback motivated the remarkable efforts made
towards developing collision operators with improved stability, as recently
reviewed by Coreixas et al. (2019). Moment-based collision operators rely on
relaxing statistical moments rather than distributions, in addition, different
relaxation times can be chosen to individually over-damp non-hydrodynamic
moments (mainly responsible for instabilities) while ensuring the correct
relaxation of hydrodynamic moments, _e.g._ density, velocity or stress tensor.
By doing so stability can be considerably enhanced while preserving physical
consistency. Nevertheless, fluid dynamics being intrinsically non-linear,
coupling artefacts arise, typically, in the form of spurious dissipative
effects.
A significant reduction of these artefacts can be obtained by considering
statistical moments expressed in the reference frame of the moving fluid
rather than in the laboratory inertial frame, referring to a class of so-
called central-moment (CM) collision operators (Geier et al., 2007, 2015; De
Rosis et al., 2018). This is the framework adopted in our modelling and now
outlined. Further details are given in the appendix A.

A key ingredient of CM-LB scheme is the shift of the particle velocities by
the local fluid velocity, thus defining a new set of local microscopic
velocities $\bar{\mathbf{c}}_{i}=\mathbf{c}_{i}-\mathbf{U}$ against which CMs
are evaluated. Here, we consider the set of CMs defined formally by

$|\mathrm{k}\rangle\equiv[\mathrm{k}_{0}\cdots\mathrm{k}_{26}]^{\top}=\mathrm{T}^{\top}|{\bar{f}}\rangle$
(27)

where the transformation matrix $\mathrm{T}$ applies to the set of
distributions $|\bar{f}\rangle\equiv[\bar{f}_{0}\cdots\bar{f}_{26}]^{\top}$
and is explicitly defined by the column vectors

$\displaystyle|\mathrm{T}_{0}\rangle=|1\rangle$
$\displaystyle|\mathrm{T}_{1}\rangle;|\mathrm{T}_{2}\rangle;|\mathrm{T}_{3}\rangle=[\bar{c}_{ix}]^{\top};[\bar{c}_{iy}]^{\top};[\bar{c}_{iz}]^{\top}$
$\displaystyle|\mathrm{T}_{4}\rangle;|\mathrm{T}_{5}\rangle;|\mathrm{T}_{6}\rangle=[\bar{c}_{ix}\bar{c}_{iy}]^{\top};[\bar{c}_{ix}\bar{c}_{iz}]^{\top};[\bar{c}_{iy}\bar{c}_{iz}]^{\top}$
$\displaystyle|\mathrm{T}_{7}\rangle;|\mathrm{T}_{8}\rangle;|\mathrm{T}_{9}\rangle=[\bar{c}_{ix}^{2}-\bar{c}_{iy}^{2}]^{\top};[\bar{c}_{ix}^{2}-\bar{c}_{iz}^{2}]^{\top};[\bar{c}_{ix}^{2}+\bar{c}_{iy}^{2}+\bar{c}_{iz}^{2}]^{\top}$
$\displaystyle|\mathrm{T}_{10}\rangle;|\mathrm{T}_{11}\rangle;|\mathrm{T}_{12}\rangle=[\bar{c}_{ix}\bar{c}_{iy}^{2}+\bar{c}_{ix}\bar{c}_{iz}^{2}]^{\top};[\bar{c}_{ix}\bar{c}_{iy}^{2}+\bar{c}_{iy}\bar{c}_{iz}^{2}]^{\top};[\bar{c}_{ix}^{2}\bar{c}_{iy}+\bar{c}_{iy}^{2}\bar{c}_{iz}]^{\top}$
$\displaystyle|\mathrm{T}_{13}\rangle;|\mathrm{T}_{14}\rangle;|\mathrm{T}_{15}\rangle=[\bar{c}_{ix}\bar{c}_{iy}^{2}-\bar{c}_{ix}\bar{c}_{iz}^{2}]^{\top};[\bar{c}_{ix}\bar{c}_{iy}^{2}-\bar{c}_{iy}\bar{c}_{iz}^{2}]^{\top};[\bar{c}_{ix}^{2}\bar{c}_{iy}-\bar{c}_{iy}^{2}\bar{c}_{iz}]^{\top}$
$\displaystyle|\mathrm{T}_{16}\rangle=[\bar{c}_{ix}\bar{c}_{iy}\bar{c}_{iz}]^{\top}$
$\displaystyle|\mathrm{T}_{17}\rangle;|\mathrm{T}_{18}\rangle;|\mathrm{T}_{19}\rangle=[\bar{c}_{ix}^{2}\bar{c}_{iy}^{2}+\bar{c}_{ix}^{2}\bar{c}_{iz}^{2}+\bar{c}_{iy}^{2}\bar{c}_{iz}^{2}]^{\top};[\bar{c}_{ix}^{2}\bar{c}_{iy}^{2}+\bar{c}_{ix}^{2}\bar{c}_{iz}^{2}-\bar{c}_{iy}^{2}\bar{c}_{iz}^{2}]^{\top};[\bar{c}_{ix}^{2}\bar{c}_{iy}^{2}-\bar{c}_{ix}^{2}\bar{c}_{iz}^{2}]^{\top}$
$\displaystyle|\mathrm{T}_{20}\rangle;|\mathrm{T}_{21}\rangle;|\mathrm{T}_{22}\rangle=[\bar{c}_{ix}^{2}\bar{c}_{iy}\bar{c}_{iz}]^{\top};[\bar{c}_{ix}\bar{c}_{iy}^{2}\bar{c}_{iz}]^{\top};[\bar{c}_{ix}\bar{c}_{iy}\bar{c}_{iz}^{2}]^{\top}$
$\displaystyle|\mathrm{T}_{23}\rangle;|\mathrm{T}_{24}\rangle;|\mathrm{T}_{25}\rangle=[\bar{c}_{ix}\bar{c}_{iy}^{2}\bar{c}_{iz}^{2}]^{\top};[\bar{c}_{ix}^{2}\bar{c}_{iy}\bar{c}_{iz}^{2}]^{\top};[\bar{c}_{ix}^{2}\bar{c}_{iy}^{2}\bar{c}_{iz}]^{\top}$
$\displaystyle|\mathrm{T}_{26}\rangle=[\bar{c}_{ix}^{2}\bar{c}_{iy}^{2}\bar{c}_{iz}^{2}]^{\top}.$

This set of vectors forms a simple relevant basis ($\mathrm{T}$ is reversible)
allowing for a suitable separation between hydrodynamic and non-hydrodynamic
moments (De Rosis, 2017). In the space of CMs, the collision step (23) now
generalizes as

$|\mathrm{k}^{*}\rangle=|\mathrm{k}\rangle-S\bigg{(}|\mathrm{k}\rangle-|\mathrm{k}^{(0)}\rangle\bigg{)}\quad\mathrm{with}~{}~{}|\mathrm{k}^{(0)}\rangle=\mathrm{T}^{\top}|{f^{\mathrm{mhd}(0)}}\rangle$
(28)

where $S$ is a diagonal matrix applied to each moment individually. Let us
point out that the BGK collision is recovered by taking $S=\omega\mathbb{Id}$.
A proper choice for $S$ is given by

$S=\mathrm{diag}[1,1,1,1,\omega,\omega,\omega,\omega,\omega,1,...,1]$ (29)

which ensures that mass and momentum are conserved by the collision operator
and that viscosity is suitably taken into account. On the other hand, higher-
order moments are relaxed to equilibrium. The post-collision distributions are
finally obtained by returning in the space of distributions with

$|\bar{f}^{*}\rangle={\mathrm{T}^{-1}}^{\top}|\mathrm{k}^{*}\rangle$ (30)

before moving on to the streaming step (22).

###  3.2 Vector-valued Lattice Boltzmann scheme for the magnetic field

We now proceed to the introduction of the LB scheme for the magnetic field.
The scheme introduced by Dellar (2002) is extended to encompass the Hall
effect. Following previous works by Croisille et al. (1995) and Bouchut
(1999), Dellar (2002) proposed to decompose the magnetic field as

$\mathbf{B}(\mathbf{x},t)=\sum_{i=0}^{M-1}\mathbf{\bar{g}}_{i}(\mathbf{x},t)$
(31)

where the sum spans a set of vector-valued densities
$\mathbf{g}_{0},\cdots,\mathbf{g}_{M-1}$ associated with the microscopic
velocities $\boldsymbol{\xi}_{0},\cdots,\boldsymbol{\xi}_{M-1}$.

The magnetic field $\mathbf{B}$ is here given by the zeroth-order moment of
$\mathbf{\bar{g}}_{i}$ suggesting that a lattice with low connectivity should
be sufficient. In practice, a D3Q7 lattice with only seven velocities (see
green arrows in Fig. 1) will prove satisfactory for Hall-MHD turbulence.
Analogously to fluid dynamics, a LB scheme can be derived to simulate the
induction equation in the form

$\mathbf{\bar{g}}_{i}(\mathbf{x}+\boldsymbol{\xi}_{i}\Delta t,t+\Delta
t)=\mathbf{\bar{g}}_{i}(\mathbf{x},t)-\omega_{B}\bigg{(}\mathbf{\bar{g}}_{i}(\mathbf{x},t)-\mathbf{g}^{(0)}_{i}(\mathbf{U},\mathbf{B})(\mathbf{x},t)\bigg{)}$
(32)

where the relaxation parameter $\omega_{\mathrm{m}}$ is now related to the
magnetic resistivity $\eta$ by

$\frac{1}{\omega_{B}}=\bigg{(}\frac{\eta}{C^{2}\Delta t}+\frac{1}{2}\bigg{)}$
(33)

with $\Delta x/\Delta t=2C$ for the D3Q7 lattice. In practice, it is desirable
that the nodes of the D3Q7 and D3Q27 lattices coincide so that the macroscopic
quantities such as $\mathbf{u}$, $\mathbf{B}$ or $\mathbf{J}$ may be exchanged
between the two lattices without interpolation. This constraint imposes that

$2C=\sqrt{3}c_{s}.$ (34)

In the context of ideal MHD, the densities at equilibrium are given by

$g_{i\alpha}^{(0)}(\mathbf{U},\mathbf{B})=W_{i}\left(B_{\alpha}+\frac{1}{C^{2}}\xi_{i\beta}(U_{\beta}B_{\alpha}-B_{\beta}U_{\alpha})\right)$
(35)

with $W_{\mathrm{center}}=1/4$ and $W_{\mathrm{face}}=1/8$ for a D3Q7 lattice.
In this way, the first-order moment

$\sum_{i=0}^{M-1}\xi_{i\alpha}g^{(0)}_{i\beta}=\mathbf{U}\otimes\mathbf{B}-\mathbf{B}\otimes\mathbf{U}$
(36)

suitably reconstruct the transport term of the induction equation. Including
the Hall correction in the magnetic induction equation therefore amounts to
upgrading the equilibrium densities so that

$\boldsymbol{\mathrm{\Lambda}}^{(0)}_{\alpha\beta}=\sum_{i=0}^{M-1}\xi_{i\alpha}g^{\mathrm{Hall}(0)}_{i\beta}=(\mathbf{U}-\alpha_{H}\mathbf{J})\otimes\mathbf{B}-\mathbf{B}\otimes(\mathbf{U}-\alpha_{H}\mathbf{J})$
(37)

which is obviously possible by now considering

$g_{i\alpha}^{\mathrm{Hall}(0)}(\mathbf{U},\mathbf{B},\mathbf{J})=W_{i}\left(B_{\alpha}+\frac{1}{C^{2}}\xi_{i\beta}((U_{\beta}-\alpha_{H}J_{\beta})B_{\alpha}-B_{\beta}(U_{\alpha}-\alpha_{H}J_{\alpha}))\right).$
(38)

Nevertheless, $\mathbf{J}$ needs to be computed, possibly from the densities,
in this expression.

An essential benefit of the LB framework is that the spatial derivatives of
the magnetic field, and thus $\mathbf{J}$, are self-consistently obtained
(within an $O(\mathrm{Ma}^{3})$ error) from the first-order moment of the
densities as

$\mathrm{J}_{\gamma}=\epsilon_{\alpha\beta\gamma}\frac{\partial
B_{\alpha}}{\partial
x_{\beta}}=-\epsilon_{\alpha\beta\gamma}\frac{\omega_{B}}{C^{2}}\bigg{(}\mathbf{\boldsymbol{\Lambda}}_{\alpha\beta}-\boldsymbol{\mathrm{\Lambda}}_{\alpha\beta}^{(0)}\bigg{)}$
(39)

where $\epsilon_{\alpha\beta\gamma}$ is the Levi-Civita tensor and
$\boldsymbol{\mathrm{\Lambda}}_{\alpha\beta}=\sum_{i=0}^{M-1}\xi_{i\alpha}\bar{g}_{i\beta}$
(Dellar, 2002).

By replacing (37) in (39) we obtain a linear system readily solvable to obtain
the current density $\mathbf{J}$, namely

$\left(\mathbb{I}+\frac{2\alpha_{H}\omega_{B}}{C^{2}}\mathbf{M}\right)\mathbf{J}=\mathbf{J_{0}}$
(40)

where

$\mathbf{M}=\begin{bmatrix}0&B_{z}&-B_{y}\\\\[3.00003pt]
-B_{z}&0&B_{x}\\\\[3.00003pt]
B_{y}&-B_{x}&0\end{bmatrix}\;\mathrm{and}\;\;\mathbf{J_{0}}=\begin{bmatrix}\boldsymbol{\mathrm{\Lambda}}_{yz}-\boldsymbol{\mathrm{\Lambda}}_{zy}-2\left(U_{y}B_{z}-U_{z}B_{y}\right)\\\\[3.00003pt]
\boldsymbol{\mathrm{\Lambda}}_{zx}-\boldsymbol{\mathrm{\Lambda}}_{xz}-2\left(U_{z}B_{x}-U_{x}B_{z}\right)\\\\[3.00003pt]
\boldsymbol{\mathrm{\Lambda}}_{xy}-\boldsymbol{\mathrm{\Lambda}}_{yx}-2\left(U_{x}B_{y}-U_{y}B_{x}\right)\end{bmatrix}.$
(41)

Obviously, a solution exists if it is possible to invert
$\widetilde{\mathbf{M}}=\mathbb{I}+\left(2\alpha_{H}\omega_{B}/C^{2}\right)\mathbf{M}$.
It can be easily verified that $\mathrm{det}(\widetilde{\mathbf{M}})\neq 0$,
which proves that this solution exists and is unique. The current density
obtained by solving (40) can then be used to compute the equilibrium densities
(38) and proceed to the collision operation.

The expression of (39) also provides a consistent approximation of the
divergence of the magnetic field. Indeed by taking the trace of the magnetic
tensor, one obtains

$\boldsymbol{\nabla}\cdot\mathbf{B}\simeq-\frac{\omega_{B}}{C^{2}}\mathrm{Tr}\left(\boldsymbol{\Lambda}_{\alpha\beta}\right)$
(42)

by noticing that
$\mathrm{Tr}\left(\boldsymbol{\Lambda}_{\alpha\beta}^{(0)}\right)=0$.
Furthermore, the ${O}(\mathrm{Ma}^{3})$ correction cancels out by taking the
trace. Therefore, this correction is pushed to a higher order, so that the
divergence-free $\boldsymbol{\nabla}\cdot\mathbf{B}=0$ corresponds with high
accuracy to the condition
$\mathrm{Tr}\left(\boldsymbol{\Lambda}_{\alpha\beta}\right)=0$ in the LB
framework (Dellar, 2002). In practice, we have checked in our LB simulation
that this last condition was maintained to machine round-off error.

###  3.3 Dimensionless formulation

In the following, the dynamical equations are first rewritten in a
dimensionless form in order to highlight a control parameter $\epsilon_{H}$
associated with the Hall parameter $\alpha_{H}=1/ne$. This parameter is then
recast in lattice units for practical LB simulations. Physical quantities in
lattice units are hereafter indicated with the superscript lbm. In lattice
units, the lattice spacing and the time-step of the scheme define the units of
length and time, respectively. In order to obtain a dimensionless induction
equation, let us normalise the magnetic field with a reference value, $B_{0}$,
the fluid velocity with $U_{0}$, the current density with $B_{0}/L_{0}$, the
length with $L_{0}$ and the time with $L_{0}/U_{0}$. In this way (9) can be
written in dimensionless form as

$\bigg{(}\frac{U_{0}B_{0}}{L_{0}}\bigg{)}\partial_{t}\mathbf{b}=\frac{1}{L_{0}}\boldsymbol{\nabla}\times\bigg{[}\bigg{(}U_{0}\mathbf{u}-\frac{\alpha_{H}B_{0}}{\sqrt{\mu_{0}}L_{0}}\boldsymbol{\nabla}\times\mathbf{b}\bigg{)}\times(B_{0}\mathbf{b})\bigg{]}+\frac{\eta
B_{0}}{L_{0}^{2}}\boldsymbol{\nabla}^{2}\mathbf{b}$ (43)

Dimensionless fields are here indicated by lowercase letters. This equation
simplifies as

$\partial_{t}\mathbf{b}=\boldsymbol{\nabla}\times\bigg{[}\bigg{(}\mathbf{u}-\epsilon_{H}\boldsymbol{\nabla}\times\mathbf{b}\bigg{)}\times\mathbf{b}\bigg{]}+\frac{1}{\mathrm{Re}_{m}}\boldsymbol{\nabla}^{2}\mathbf{b}$
(44)

by defining the magnetic Reynolds number $\mathrm{Re}_{m}=U_{0}L_{0}/\eta$ and
the dimensionless Hall control parameter

$\epsilon_{H}=\frac{\alpha_{H}B_{0}}{\sqrt{\mu_{0}}L_{0}U_{0}}.$ (45)

We can proceed in the same way for the fluid momentum equation, where the
reference scales are the same as those adopted for the induction equation.
Therefore,

$\rho\frac{U_{0}^{2}}{L_{0}}\left[\partial_{t}\mathbf{u}+\left(\mathbf{u}\cdot\boldsymbol{\nabla}\right)\mathbf{u}\right]=-\frac{1}{L_{0}}\boldsymbol{\nabla}\rho
c_{s}^{2}+\rho\nu\frac{U_{0}}{L_{0}^{2}}\boldsymbol{\nabla}^{2}\mathbf{u}+\frac{B_{0}^{2}}{L_{0}}\left(\boldsymbol{\nabla}\times\mathbf{b}\right)\times\mathbf{b}$
(46)

which gives

$\partial_{t}\mathbf{u}+\left(\mathbf{u}\cdot\boldsymbol{\nabla}\right)\mathbf{u}=-\frac{1}{\mathrm{Ma}^{2}}\frac{1}{\rho}{\boldsymbol{\nabla}\rho}+\frac{1}{\mathrm{Re}}\boldsymbol{\nabla}^{2}\mathbf{u}+\left(\frac{V_{A}}{U_{0}}\right)^{2}\left(\boldsymbol{\nabla}\times\mathbf{b}\right)\times\mathbf{b}$
(47)

where the control parameters are the Mach number $\mathrm{Ma}=U_{0}/c_{s}$,
the (fluid) Reynolds number $\mathrm{Re}=U_{0}L_{0}/\nu$ and $V_{A}/U_{0}$
with the Alfvèn velocity $V_{A}=B_{0}/\sqrt{\rho}$. The Hall number
$\epsilon_{H}$ is given in lattice units by

$\epsilon_{H}^{\mathrm{lbm}}=\frac{\alpha_{H}}{\sqrt{\mu_{0}}}\frac{\left[B_{0}/B^{\mathrm{lbm}}_{0}\right]}{\left[L_{0}/L^{\mathrm{lbm}}_{0}\right]\left[U_{0}/U^{\mathrm{lbm}}_{0}\right]}=\epsilon_{H}\frac{U^{\mathrm{lbm}}_{0}L_{0}^{\mathrm{lbm}}}{B_{0}^{\mathrm{lbm}}}.$
(48)

If one considers that the reference velocity $U_{0}$ corresponds to the Alfvèn
velocity ($U_{0}=V_{A}$) and $\rho\simeq 1$ for simplicity, one obtains that
$U_{0}^{\mathrm{lbm}}=B_{0}^{\mathrm{lbm}}$ and

$\epsilon_{H}^{\mathrm{lbm}}=\epsilon_{H}L_{0}^{\mathrm{lbm}}=\epsilon_{H}N$
(49)

with $N=L_{0}/\Delta x$ being the number of lattice points per reference
length $L_{0}$. The Hall number $\epsilon_{H}$ can also be viewed as the ratio
of two reference scales as

$\epsilon_{H}=\frac{V_{A}}{U_{0}L_{0}}\sqrt{\frac{m}{\mu_{0}ne^{2}}}=\frac{L_{H}}{L_{0}}$
(50)

with

$L_{H}=\frac{V_{A}}{U_{0}}\sqrt{\frac{m}{\mu_{0}ne^{2}}}.$ (51)

In lattice units,

$\epsilon_{H}^{\mathrm{lbm}}=\frac{L_{H}}{L_{0}}~{}\frac{U^{\mathrm{lbm}}_{0}L_{0}^{\mathrm{lbm}}}{B_{0}^{\mathrm{lbm}}}=\frac{L_{H}}{\Delta
x}~{}\frac{U^{\mathrm{lbm}}_{0}}{B_{0}^{\mathrm{lbm}}}=\frac{U^{\mathrm{lbm}}_{0}L_{H}^{\mathrm{lbm}}}{B_{0}^{\mathrm{lbm}}}.$
(52)

If $U_{0}=V_{A}$, $L_{H}$ is equal to the ion inertial length $d_{i}$ (or ion
skin depth). In that case, $\epsilon_{H}^{\mathrm{lbm}}=L_{H}^{\mathrm{lbm}}$
and corresponds to the number of lattice points per ion inertial length. The
dynamics is affected by the Hall effect at scales $\ell$ smaller that $L_{H}$.

###  3.4 CFL condition for Hall-MHD turbulence

The Courant–Friedrichs–Lewy (CFL) condition (Lewy et al., 1928) determines,
for an explicit time-marching scheme, the maximum time-step for convergence,
as

$\Delta t\leq\Delta x/c_{\mathrm{max}}$ (53)

where $c_{\mathrm{max}}$ refers to the largest speed at which a signal
propagates in the solution. In the context of Hall-MHD, $c_{\mathrm{max}}$
should be identified with the largest phase speed of whistler waves. When the
plasma dynamics along the magnetic field $B$ is dominant, the phase speed of
whistler waves varies as $c_{w}(k)=kV_{A}^{2}/\omega_{ci}$ with the wavenumber
$k$. $V_{A}$ is the Alfvèn velocity and $\omega_{ci}$ is the ion cyclotron
frequency. In physical units, $\omega_{ci}=eB/m_{i}$ and
$V_{A}=B/\sqrt{\mu_{0}nm_{i}}$ with $m_{i}$ being the mass of ions, the Hall
parameter $\alpha_{H}=1/ne$. Therefore, one obtains that the time-step
decreases quadratically with the grid spacing as

$\Delta t\leq\frac{\mu_{0}(\Delta x)^{2}}{\pi~{}\alpha_{H}B}$ (54)

by considering the maximum attainable wavenumber $k_{\mathrm{max}}=\pi/\Delta
x$ in a context of Hall-MHD turbulence. This condition rewrites by accounting
for the rescaling of the magnetic field by ${\mu_{0}}^{1/2}$

$\Delta t\leq\frac{(\Delta x)^{2}}{\pi~{}(\alpha_{H}/\sqrt{\mu_{0}})~{}B_{0}}$
(55)

which finally yields in lattice units

$1\leq\frac{1}{\pi~{}\epsilon_{H}^{\mathrm{lbm}}L_{0}^{\mathrm{lbm}}U_{0}^{\mathrm{lbm}}}=\frac{1}{\pi}~{}\frac{B_{0}^{\mathrm{lbm}}}{\epsilon_{H}\left(N{U_{0}^{\mathrm{lbm}}}\right)^{2}}$
(56)

where (45) is used to express $\alpha_{H}/\sqrt{\mu_{0}}$ and
$N=L_{0}^{\mathrm{lbm}}$. If $U_{0}=V_{A}$, the CFL condition for whistler
waves reformulates as a condition on the Mach number
${\mathrm{Ma}=\sqrt{3}U_{0}^{\mathrm{lbm}}}$, which is written

$\mathrm{Ma}\leq\left(\frac{\sqrt{3}}{\pi}\right)\frac{1}{\epsilon_{H}N^{2}}.$
(57)

This condition resembles the quadratic dependence of the time step on the
resolution obtained with conventional CFD methods (Gómez et al., 2010). It
also confirms that Hall-MHD turbulence is computationally very demanding due
to the presence of whistler waves.

Figure 2: Three-dimensional rendering of the initial condition as indicated in
(62) and (63). The magnitudes of the fluid velocity (left), magnetic field
(center) and density of electric current (right) are here displayed for
$N=128$.

##  4 Results

The LB scheme introduced in the previous section was validated against an
exact analytical solution of the incompressible and dissipative Hall-MHD
equations (Xia & Yang, 2015). This solution provides a valuable benchmark for
evaluating the accuracy and convergence for different values of the control
parameters. A second validation was performed in the MHD regime at a high
Reynolds number with a well-established pseudo-spectral solver widely used for
turbulent plasma simulations, namely the Geophysical High-Order Suite for
Turbulence (ghost, Rosenberg et al. (2020)). Finally, the performance of the
implementation was evaluated with a typical problem for turbulent plasma
flows, namely the Orszag-Tang vortex (Orszag & Tang, 1979). This allowed us to
observe the reliability of our code in simulating turbulent flows at high
Reynolds numbers, and in reproducing the transition from MHD to Hall-MHD at
sub-ion scales.

Parameters | Values
---|---
Resolution: $N$ | $32,\ 64,\ 96,\ 128$
Mach number: $\mathrm{Ma}\hphantom{0}[\times 10^{-2}]$ | $1.0,\ 0.7,\ 0.5,\ 0.3$
Kinematic viscosity: $\nu\hphantom{0}[\times 10^{-3}]$000 | $1.0,\ 0.5,\ 0.33,\ 0.25,\ 0.2,\ 0.17,\ 0.14,\ 0.13,\ 0.11$
Table 1: Parameters of LB simulations. The magnetic Prandtl number is kept
unitary. The kinematic viscosity is given in dimensionless units, _i.e._
normalised by $U_{0}L_{0}$, which means that the Reynolds number
$\mathrm{Re}=1/\nu$.

###  4.1 Exact solution of the dissipative Hall-MHD

MHD in the Hall regime do not offer many reference test-cases to evaluate the
numerical efficiency and accuracy of our numerical modelling. In the three-
dimensional turbulent regime, existing simulations are very few due to the
very high computational cost, and are generally based pseudo-spectral codes
(Ferrand et al., 2022; Meyrand & Galtier, 2012; Gómez et al., 2010; Yadav et
al., 2022) which integrate the dynamical equations in the Fourier space rather
than in the physical space. Nevertheless, Mahajan & Krishan (2005) derived
analytically an exact solution for the non-dissipative Hall-MHD equations,
which has then been extended by Xia & Yang (2015) to include dissipative
effects. This solution will be used in the following to test the stability and
convergence of our LB scheme. Since this analytical solution encompasses
dissipative effects, it will also give us the opportunity to evaluate
quantitatively the artificial increase of the dissipation due to the numerical
modelling.

The solution provided by Xia & Yang (2015) may be rewritten in a dimensionless
form (see §3.3) as

$\mathbf{u}(\mathbf{x},t)=\mathbf{u}^{\prime}(\mathbf{x},t)\quad\mathrm{and}\quad\mathbf{b}(\mathbf{x},t)=\hat{\mathbf{e}}_{z}+\mathbf{b}^{\prime}(\mathbf{x},t)$
(58)

where the fluctuating velocity and magnetic fields are damped circular
polarised waves given respectively by

$\displaystyle\mathbf{u}^{\prime}(\mathbf{x},t)$
$\displaystyle=\big{[}B(\hat{\mathbf{e}}_{x}+i\hat{\mathbf{e}}_{y})\exp(i\mathrm{k}z-i\omega_{\pm}t)$
(59)
$\displaystyle+C(\hat{\mathbf{e}}_{y}+i\hat{\mathbf{e}}_{z})\exp(i\mathrm{k}x)$
$\displaystyle+A(\hat{\mathbf{e}}_{z}+i\hat{\mathbf{e}}_{x})\exp(i\mathrm{k}y)\big{]}e^{-\nu\mathrm{k}^{2}t}$

and

$\mathbf{b}^{\prime}(\mathbf{x},t)=\alpha_{\pm}\mathbf{u}^{\prime}(\mathbf{x},t)$
(60)

in complex notations. The amplitudes $A$, $B$ and $C$ are arbitrary real
values. The ambient magnetic field is here assumed to be oriented along the
unit vector $\hat{\mathbf{e}}_{z}$. Since the dynamical equations only consist
of real variables, either the imaginary part or real part is a solution. The
pulsation $\omega_{\pm}=-\alpha_{\pm}\mathrm{k}$, where $\alpha_{\pm}$ depends
itself on the wavenumber $\mathrm{k}$ as

$\alpha_{\pm}=-\frac{1}{2}\epsilon_{H}\mathrm{k}\pm\sqrt{\frac{\epsilon_{H}^{2}\mathrm{k}^{2}}{4}+1}.$
(61)

The magnetic Prandtl number is unitary in this solution, and it is supposed
that the reference velocity corresponds to the Alfvèn velocity, _i.e._
$U_{0}=V_{A}$ in (47). Finally, it should be mentioned that this exact
solution is established in a strictly incompressible framework, which imposes
on our simulations a (very) low Mach number so that relative density
fluctuations are negligible. In our investigations, the Hall-MHD equations
were integrated in a cubic box of size $L_{0}=2\pi$. The evolution of the
velocity field is deterministic from the initial condition

$\displaystyle\mathrm{u}^{\mathrm{lbm}}_{x}(\mathbf{x},0)={U_{0}^{\mathrm{lbm}}}\left(B\sin\left(\frac{4\pi
z_{k}}{N}\right)+A\cos\left(\frac{4\pi y_{j}}{N}\right)\right)$ (62)
$\displaystyle\mathrm{u}^{\mathrm{lbm}}_{y}(\mathbf{x},0)={U_{0}^{\mathrm{lbm}}}\left(B\cos\left(\frac{4\pi
z_{k}}{N}\right)+C\sin\left(\frac{4\pi x_{i}}{N}\right)\right)$
$\displaystyle\mathrm{u}^{\mathrm{lbm}}_{z}(\mathbf{x},0)={U_{0}^{\mathrm{lbm}}}\left(C\cos\left(\frac{4\pi
x_{i}}{N}\right)+A\sin\left(\frac{4\pi y_{j}}{N}\right)\right)$

expressed in lattice units with $A=0.3$, $B=0.2$, $C=0.1$ and $N$ being the
number of lattice nodes per reference length $L_{0}$. The reference velocity
$U_{0}^{\mathrm{lbm}}$ is related by construction to the Mach number through
$U_{0}^{\mathrm{lbm}}=\mathrm{Ma}/\sqrt{3}$. The magnetic field is initially
proportional to the fluid velocity with

$\displaystyle\mathrm{b}^{\mathrm{lbm}}_{x}(\mathbf{x},0)=\alpha_{+}\mathrm{u}^{\mathrm{lbm}}_{x}(\mathbf{x},0)$
(63)
$\displaystyle\mathrm{b}^{\mathrm{lbm}}_{y}(\mathbf{x},0)=\alpha_{+}\mathrm{u}^{\mathrm{lbm}}_{y}(\mathbf{x},0)$
$\displaystyle\mathrm{b}^{\mathrm{lbm}}_{z}(\mathbf{x},0)=\alpha_{+}\mathrm{u}^{\mathrm{lbm}}_{z}(\mathbf{x},0)+U_{0}^{\mathrm{lbm}}$

since $U_{0}=V_{A}$. For simplicity, the initial density is uniformly equal to
one. The (normalised) Hall parameter is fixed at $\epsilon_{H}=1$, that is
$\epsilon_{H}^{\mathrm{lbm}}=N$ according to (49). This value ensures that the
solution is affected by the Hall effect with $L_{H}=L_{0}$ from (50). A three-
dimensional rendering of the initial conditions expressed in (62) and (63) is
represented in Fig. 2. With this initial condition, the current density
$\mathbf{J}=\boldsymbol{\nabla}\times\mathbf{B}$ is non-zero at $t=0$. The
parameters of the different simulations are reported in Tab. 1. The Mach
number is small enough to approach the incompressible limit and reduce the
discretisation error of the LB method. The CFL condition imposed by this
solution is also satisfied. Finally, let us mention that analogous simulations
were performed with the phase speed $\alpha_{-}$ yielding very similar results
on accuracy and stability. However, the phase speed is much larger in the
latter case, therefore requiring a significant reduction of the Mach number
(with $\epsilon_{H}=1$). Only results for $\alpha_{+}$ and the velocity field
($\mathbf{b}=\alpha_{+}\mathbf{u}+\hat{\mathbf{e}}_{z}$) are presented in the
following.

Figure 3: Temporal evolution of the velocity magnitude
$|\mathbf{u}|(\mathbf{0},t)$. Comparison between the analytical (dashed line)
and the numerical solutions (symbols) at different lattice resolutions. The
Mach number $\mathrm{Ma}=0.003$ and the kinematic viscosity $\nu=3.3\cdot
10^{-4}$.

###  4.2 Stability and incompressibility

The stability was tested by exploring the parameter space defined by the Mach
number, the lattice resolution and the kinematic viscosity (see Tab. 1). The
analytical solution introduced by Xia & Yang (2015) renders strictly zero the
nonlinear terms of the incompressible dissipative Hall MHD equations. In
practice, physical instabilities triggered by numerical errors naturally
develop and grow in time in simulations when the viscosity is too small,
eventually resulting in a transition to a turbulent state. The numerical
stability and accuracy were measured in situations where the viscosity was
sufficiently high to prevent this transition to turbulence.

Figure 4: Temporal evolution of the velocity magnitude
$|\mathbf{u}|(\mathbf{0},t)$ for different values of the resolution ($N$) and
viscosity ($\nu$) at fixed Mach number $\mathrm{Ma}=0.007$.

A typical temporal evolution of the velocity at a fixed location is shown in
Fig. 3. The solution appears as a damped wave propagating in the direction of
the ambient magnetic field. The amplitude and the phase of the solution are
well captured in the LB simulation. The results obtained for different
resolutions and viscosity values at Mach number $\mathrm{Ma}=0.007$ are shown
in Fig. 4 over 10 wave periods. All simulations remained numerically stable in
the explored range of parameters. The temporal averages of relative density
fluctuations at different values of Mach number and kinematic viscosity are
displayed in Fig. 5. The level of these relative fluctuations is typically of
order $10^{-7}$ – $10^{-6}$ indicating that the incompressible limit is
approached very satisfactorily in these simulations. Furthermore, the results
confirm that the amplitude of density fluctuations decreases with the Mach
number.

Figure 5: Relative density fluctuations at different values of the Mach number
($\mathrm{Ma}$) and kinematic viscosity ($\nu$) as a function of the
resolution ($N$). In our simulations, the reference density $\rho_{0}$ is
fixed at unity.

###  4.3 Dispersion and dissipation errors

The dispersion and dissipation errors of the LB scheme are now evaluated. In
this analysis, the dispersion error (or phase error) is computed by evaluating
the shift in time between the maxima of the numerical solution and the maxima
of the analytical wave solution (see Fig. 3). Therefore, by noting
$t^{\mathrm{max}}_{i}$ and $\bar{t}^{\mathrm{max}}_{i}$ the respective
positions in time of the maxima of the numerical and analytical solution (at a
fixed location), the average value of the relative dispersion error is defined
as

$\varepsilon_{\phi}=1-\frac{1}{M}\sum_{i=0}^{M-1}\frac{t^{\mathrm{max}}_{i+1}-t^{\mathrm{max}}_{i}}{\bar{t}^{\mathrm{max}}_{i+1}-\bar{t}^{\mathrm{max}}_{i}}$
(64)

over $M$ oscillating periods. In practice, we have used $M=10$. As expected,
it is observed in Fig. 6 that the dispersion error is very small and decreases
when increasing the resolution $N$ of the grid, with a power-scaling law close
to $1/N^{2}$. This confirms a second-order accuracy of the LB scheme. It is
also found that the dispersion error manifests a constant behaviour when
changing the Mach number, and does not seem to be affected by the value of the
kinematic viscosity either. Let us remark that some results differ from the
global trend, certainly due to the premise of (physical) instabilities at the
lowest viscosity. After realigning in phase the numerical and analytical
solution, the (relative) dissipation error is evaluated by comparing the
velocity magnitude of the numerical and analytical solutions, _i.e._

$\varepsilon_{r}=\sum_{i=0}^{M-1}\frac{\left|\mathrm{u}(t_{i})-\mathrm{\bar{u}}(t_{i})\right|}{\mathrm{\bar{u}}(t_{i})}.$
(65)

The dissipation error gives a first measure of the numerical dissipation. Two
different scaling behaviours are considered, namely the so-called acoustic and
diffusive scaling (Krueger et al., 2016). The acoustic scaling consists of
keeping the Mach number fixed and observing the convergence rate of the error
$\varepsilon_{r}$ for different Reynolds numbers, and as a function of the
resolution (see Fig. 7). On the other hand, the diffusive scaling is obtained
by keeping the lattice viscosity fixed (see Fig. 8). The behaviour of the
numerical solution is consistent between the two regimes, showing a
convergence of the dissipation error with respect to the grid resolution
$\propto 1/N^{2}$ as expected for a second-order accurate scheme.

Figure 6: Relative dispersion error as defined by (64) for different values of
the kinematic viscosity $\nu$ and the Mach number $\mathrm{Ma}$. The error
decreases as $N^{-2}$ as expected for a LB scheme. Figure 7: Acoustic scaling
($\mathrm{Ma}$ constant) of the relative dissipation error for different
values of the kinematic viscosity ($\nu$) at fixed Mach number
($\mathrm{Ma}$). The error decreases as $N^{-2}$.

One of the advantages of dealing with a _dissipative_ solution of the Hall-MHD
equations is the possibility to identify an effective viscosity $\tilde{\nu}$
related to the damping $\propto\exp(-\tilde{\nu}k^{2}t)$ of the numerical
solution, and compare it to the physical viscosity $\nu$. By decomposing
$\tilde{\nu}$ into the sum of a physical and a numerical viscosity,
$\tilde{\nu}=\nu+\nu_{\mathrm{num}}$, the ratio between the numerical and
physical viscosities writes

$\varepsilon_{\nu}=\frac{\nu_{\mathrm{num}}}{\nu}=\frac{\tilde{\nu}-\nu}{\nu}.$
(66)

The results obtained for the viscosity error $\varepsilon_{\nu}$ are shown in
Fig. 9. It is found that the numerical viscosity represents only a small
percentage of the physical viscosity and decreases as $1/N^{2}$ with the
resolution, which is once again consistent with a second-order accuracy of the
LB scheme. Interestingly, it is observed that the (relative) viscosity error
remains independent of the physical viscosity and the Mach number, but only
depends on the lattice resolution.

Figure 8: Diffusive scaling ($\nu$ constant) of the (relative) dissipation
error for different values of the Mach number ($\mathrm{Ma}$) at fixed
kinematic viscosity ($\nu$). Figure 9: Scaling of the ratio between the
numerical and kinematic viscosities $\varepsilon_{\nu}=\nu^{\mathrm{num}}/\nu$
with the grid resolution ($N$) for different Mach numbers ($\mathrm{Ma}$) and
kinematic viscosity. The second-order accuracy of the LB algorithm is
highlighted by the black lines, _i.e._ $\varepsilon_{\nu}\sim O(\Delta
x^{2})$.

Finally, despite Dellar (2002) showed that a D3Q7 lattice was sufficient to
reliably account for the dynamics of each component of the magnetic field, LB
simulations with an enhanced connectivity have been performed to investigate
if a more isotropic representation of the magnetic densities would notably
improve the level of accuracy of the algorithm (Silva & Semiao, 2014).
Interestingly, our results showed no significant improvement when upgrading
the magnetic lattice to D3Q15 or D3Q27 lattices (see Fig. 1). A plausible
explanation lies in the fact that the magnetic field is represented as a zero-
th order moment of densities for each component (see (31)). Therefore, few
degrees of freedom are certainly sufficient to accurately reconstruct the
moments and describe the magnetic field dynamics.

###  4.4 Comparison with pseudo-spectral simulations of MHD turbulence

In this section, comparisons are made with the ghost pseudo-spectral solver
for high-resolution MHD simulations within the framework of the Orszag-Tang
(OT) vortex problem (Orszag & Tang, 1979). The latter is an initial value
problem, therefore allowing for a direct comparison between different
numerical simulations. The OT solution is also often considered as a
prototypical flow to study freely-evolving MHD turbulence. The ghost solver is
widely used for various geophysical and astrophysical studies (Feraco et al.,
2021; Marino et al., 2014; Mininni et al., 2002, 2003, 2006; Gómez et al.,
2010). It is a well-established open-source community code available on
https://github.com/pmininni/GHOST. ghost is a hybrid MPI/OpenMP/CUDA-
parallelized code, which proved high performance, robust results and an
optimal scaling up to hundreds of thousand processors. It relies on a second-
order Runge-Kutta scheme for the time integration. As a pseudo-spectral de-
aliased code, it provides infinite accuracy in resolved spatial scales
(Patterson & Orszag, 1971). The dynamics of the OT vortex problem are
initialised for the velocity and magnetic fields by

$\displaystyle\mathbf{U(\mathbf{x},0})$ $\displaystyle=$ $\displaystyle
U_{0}\left[-2\sin{y}~{};~{}2\sin{x}~{};~{}0\right]$
$\displaystyle\mathbf{B(\mathbf{x},0})$ $\displaystyle=$ $\displaystyle
B_{0}\left[-2\sin{2y}+\sin{z}~{};~{}2\sin{x}+\sin{z}~{};~{}\sin{x}+\sin{y}\right]$

with $U_{0}=1$ and $B_{0}=0.8$ in a cubic box of size $2\pi$.

In the performed simulation, the Reynolds number attains
$\mathrm{Re}=UL/\nu\simeq 1600$ when the flow reaches its peak of dissipation.
The dissipation of energy is here defined as
$\epsilon=-\nu\langle|\nabla\times\mathbf{U}|^{2}\rangle-\eta\langle|\mathbf{J}|^{2}\rangle$
and encompasses both the kinetic and magnetic dissipation with $\nu/\eta=1$.
$U$ refers to the r.m.s velocity and $L=2\pi\int k^{-1}E_{v}(k)dk/\int
E_{v}(k)dk$ is the integral length scale, where $E_{v}(k)$ is the energy
spectrum of the velocity field. The Mach number is fixed at
$\mathrm{Ma}=0.025$. The number of grid points in each direction is $N=512$.

The time evolution of the mean magnetic dissipation, and kinetic and magnetic
energies are shown in Fig. 10 from the two LB and pseudo-spectral simulations
of the same problem. The agreement is very satisfactory. Only a slight
underestimation of the magnetic dissipation is temporarily observed after the
peak. Let us recall at this stage that the current density
$\mathbf{J}=\nabla\times\mathbf{B}$ is directly accessible from the magnetic
densities in the LB simulation, and is not obtained by differentiating the
magnetic field.

Figure 10: (Left) Time evolution of the mean magnetic dissipation
$\propto\langle|\mathbf{J}|^{2}\rangle$ in freely-evolving MHD turbulence for
a LB simulation ($N=512$) and a pseudo-spectral simulation (N=$512$) performed
with the ghost solver. (Right) Time evolution of the mean kinetic ($E_{v}$)
and magnetic $(E_{B}$) energies. Figure 11: Kinetic (left) and magnetic
(right) energy spectra of MHD turbulence at the peak of magnetic dissipation.
The spectra are normalised by the total kinetic and magnetic energies
respectively. Comparison between LB simulations at two different resolutions
($N=512$, $N=768$) and a de-aliased pseudo-spectral simulation ($N=512$)
performed with the ghost solver.

A more detailed comparison is provided by a scale-by-scale decomposition. The
kinetic and magnetic energy spectra are displayed in Fig.11 at the peak of
magnetic dissipation. The kinetic energy spectrum of the LB simulation is
obviously over-damped at high wavenumbers. This is related to a known drawback
of the moment-based collision operator, which ensures a higher stability
(compared to the standard BGK collision operator) but at the cost of an
enhanced numerical dissipation (Coreixas et al., 2019). However, when
increasing the spatial resolution to $N=768$, the numerical dissipation is
reduced and the spectrum comes very close to the pseudo-spectral solution.
Concerning the magnetic energy spectrum, it is found that the results from
both simulations perfectly match, reflecting the fact that BGK collision
operator adopted for the magnetic scheme does not add numerical dissipation
(as compared to the pseudo-spectral simulation). It should also be noted that,
while the maximum wave-number is limited to $k_{\mathrm{max}}=N/3$ (due to the
2/3 rule for de-aliasing) in pseudo-spectral simulations, the range of
resolved scales reaches the Nyquist cut-off $k_{\mathrm{max}}=N/2$ in LB
simulations. A particular attention is now paid to the wavenumber-by-
wavenumber energy budget of the MHD equations. Starting from (11) and (13),
the (total) energy flux across wavenumber $k$ can be identified as

$\mathcal{S}_{\mathrm{MHD}}(k)=\sum_{|\mathbf{k^{\prime}}|<k}\Re{\left[\mathcal{F}({\mathbf{U}})^{*}\cdot\left(\mathcal{F}({\mathbf{U}\cdot\nabla\mathbf{U}})-\mathcal{F}({\mathbf{J}\times\mathbf{B}})\right)-\mathcal{F}({\mathbf{B}})^{*}\cdot\mathcal{F}(\nabla\times({\mathbf{U}\times\mathbf{B}}))\right]}$
(67)

whereas the (total) dissipation in the range $[0,k[$ is given by

$\mathcal{D}(k)=\sum_{|\mathbf{k^{\prime}}|<k}\nu k^{\prime
2}|\mathcal{F}({\mathbf{U}})|^{2}+\eta k^{\prime
2}|\mathcal{F}({\mathbf{B}})|^{2}$ (68)

where $\mathcal{F}(\cdot)$ means the Fourier transform and ∗ is the complex
conjugate. The wavenumber-by-wavenumber energy budget then writes

$\partial_{t}\sum_{|\mathbf{k}^{\prime}|<k}E(\mathbf{k}^{\prime})=-\mathcal{S}_{\mathrm{MHD}}(k)-\mathcal{D}(k).$
(69)

Let us mention that the contribution of the pressure term (not shown here) is
negligible in the present context. The results obtained for the LB and pseudo-
spectral simulations (with $N=512$) are detailed in Fig. 12. A very
satisfactory agreement is obtained in particular for the non-linear energy
transfer terms in the entire range of resolved wavenumbers. The slight over-
dissipation of the LB scheme is again evidenced for the dissipation term
$\mathcal{D}(\mathrm{k})$ at very high wavenumbers.

Figure 12: Various third-order energy fluxes across wavenumber $k$ as reported
in (68). Comparison is made between the LB (left) and the de-aliased pseudo-
spectral (right) ghost simulations. Let us notice that
$\partial_{t}\sum_{|\mathbf{k}^{\prime}|<k}E(\mathbf{k}^{\prime})=-\mathcal{S}_{\mathrm{MHD}}(k)-\mathcal{D}(k)<0$
as expected for freely-decaying MHD turbulence.

##  5 High-Resolution simulation of 3D Hall-MHD

Run | $N$ | $\mathrm{Ma}\,[\times 10^{-4}]$ | $\mathrm{Re}$ | $\mathrm{Pr}_{m}$ | $\epsilon_{H}$ | $t_{\mathrm{tot}}/\tau_{0}$ | $t_{\mathrm{peak}}/\tau_{0}$
---|---|---|---|---|---|---|---
I | 512 | 7.0 | 4400 | 1 | 0.0025 | 48.2 | 31.8
II | 512 | 1.0 | 5240 | 1 | 0.01 | 35.8 | 33.6
III | 512 | 0.625 | 7150 | 1 | 0.025 | 29.6 | 26.5
IV | 768 | 0.6 |  | 1 | 0.015 | 36.6 | 32.0
Table 2: Parameters of Hall-MHD turbulence runs. $\mathrm{Re}$, $\mathrm{Ma}$
and $\mathrm{Pr}_{m}$ denote respectively the Reynolds number (at the peak of
magnetic dissipation), the initial Mach number and the magnetic Prandtl
number. The (dimensionless) Hall parameter is $\epsilon_{H}$. The number of
grid points per dimension is $N$. The total duration of the run is
$t_{\mathrm{tot}}/\tau_{0}$ and the time at which the peak of current density
occurs is $t_{\mathrm{peak}}/\tau_{0}$ in units of the reference time scale
$L_{0}/U_{0}$. The Mach number satisfies the CFL condition (57) imposed by
whistler waves. Figure 13: Evolution of the magnetic dissipation over time for
the three simulations performed with the OT initial condition (see Tab.2). The
shaded areas around the peak of the current density (black dashed line)
correspond to the range over which the energy spectra in Fig. 14 have been
averaged.

In the Hall-MHD regime, three runs were performed in a triply periodic cubic
lattice of size $L=2\pi$ with resolution $512^{3}$ starting from the same
initial condition of the OT vortex problem (as described in the previous
section). The parameters of the runs are kept fixed except for the Hall
parameter, which is different for each simulation (see runs I–III in Tab. 2),
and for the Mach number, which is adjusted according to the Hall parameter to
respect the CFL condition related to whistler waves (see 57). Let us note that
the Reynolds number is computed at the peak of magnetic dissipation (indicated
by the vertical dashed lines in Fig 13), at which the OT solution is expected
to reach a developed turbulent state. Our flame code relies on a multi-GPU
(Graphic Processing Unit) implementation of the LB scheme in order to reach
high resolution with affordable computational times. Massive multi-threading
is handled by using the OpenCL (Open Computing Language) framework allowing
for the portability of a unique code across different targeted devices. The
spatial domain is split along a single direction and each GPU is assigned a
sub-domain. A one-to-one mapping operates between the host CPU processes and
the GPUs. Therefore, the exchange of boundary nodes between the GPUs is
handled through memory transfers with the CPU processes and a message-passing
interface (MPI) between these latter. Simulations were run on a workstation
equipped with NVIDIA A100-40Gb GPUs at the CINECA computing center (Italy).
For a $512^{3}$ lattice dimension, 3 GPUs were used in parallel resulting in a
computational speed of about 20 iterations per second, or equivalently, in
$2.7$ BLUPS (Billions of Lattice-node Updates Per Second). In practice, this
led to a wall-clock computational time of 10, 55 and 69 hours respectively for
the three runs to pass over the peak of magnetic-energy dissipation (see Fig.
13). The computational times reported above are comprehensive of the time
required to transfer the three-dimensional vector fields ($\mathbf{u}$,
$\mathbf{B}$ and $\mathbf{J}$) between the CPUs and GPUs, and perform post-
processing operations such as the tracking of the mean kinetic and magnetic
energies, and mean energy dissipation rates. All computations were performed
in double precision. A rendering of the large-scale fields $\mathbf{u}$ and
$\mathbf{B}$ is depicted in Fig. 15 for a simulation at higher resolution
$N=768$ (run IV in Tab. 2), always taken at the peak of magnetic-energy
dissipation. The three-dimensional visualization are reported together with
the kinetic and magnetic spectra, the latter clearly showing two different
regimes above and below the ion inertial length $d_{i}$. At the same time in
Fig. 16 the small-scale fields, represented by the electric current density
$\mathbf{J}$ and the vorticity
$\boldsymbol{\omega}=\boldsymbol{\nabla}\times\mathbf{u}$, highlight coherent
structures such as current sheets, Kelvin-Helmholtz instabilities and
vortices, with disordered small-scale structures in the electric current
density characteristic of the Hall effect (Miura & Araki, 2014). Furthermore,
it is established that increasing the intensity of the Hall effect produces a
faster development of turbulence in plasmas due to the presence of whistler
and Hall-drift waves, propagating faster than classical Alfvèn waves in ideal
MHD (Huba, 2003). This is consistent with the behaviour captured in Fig. 13
for the three runs with increasing Hall parameters. The kinetic and magnetic
energy spectra are plotted in Fig. 14 for each run. The spectra have been
averaged over a time interval (around the peak of dissipation) indicated by
the shaded areas in Fig. 13. As expected, the increase of the Hall parameter
$\epsilon_{H}$ is responsible for the displacement of the Hall length-scale
$L_{H}$ towards large scales (indicated by the vertical dash-dotted line in
Fig. 14) and hence, for a cut-off of the Kolmorogov’s $k^{-5/3}$ scaling law
at larger and larger scales for both the kinetic and magnetic energy spectra.
The most surprising and promising result of these simulations is provided by
the behaviour of the magnetic spectra in the Hall-MHD regime. At wavenumbers
$k>k_{H}$, the spectrum develops (as $\epsilon_{H}$ increases) a power-law
scaling that is found in perfect agreement with the $k^{-2.73}$ scaling
obtained from solar wind measurements at sub-ion scales (Kiyani et al., 2015).

Figure 14: Kinetic (blue) and magnetic (orange) spectra corresponding to three
different Hall-MHD turbulence regimes (see Tab. 2). The black dashed and
dotted lines indicate respectively the $-5/3$ and $-2.73$ slopes in reference
to solar wind measurements (Kiyani et al., 2015). All the spectra are taken at
the peak of magnetic dissipation. The (black) vertical dash-dotted lines
indicate the cross-over wavenumbers between MHD and Hall-MHD regimes given by
$k_{H}=2\pi/(\epsilon_{H}L_{0})$ with $L_{0}=2\pi\int E_{v}(k)k^{-1}dk/\int
E_{v}(k)dk$. Figure 15: Results of a 3D Hall-MHD simulation at $N=768^{3}$
taken at the peak of dissipation. Upper panels show a three-dimensional
visualization of the normalized velocity field magnitude
$v^{\prime}=(v-\bar{v})/\sigma_{v}$ with $v=|\mathbf{v}|$ (upper left) and the
relative power density spectrum (upper right). The same in the two bottom
panels for the normalized magnetic field $B^{\prime}$. The slopes $k^{-5/3}$
and $k^{-2.73}$ are given as a reference. Figure 16: Three-dimensional
rendering of the normalized current density
$J^{\prime}=(J-\bar{J})/\sigma_{J}$, with $J=|\mathbf{J}|$ (left), and
vorticity $\omega^{\prime}$ (right) taken at the same time as in Fig. 15.

In the MHD regime, pseudo-spectral simulation relies on the incompressibility
condition that allows for large time-steps, whereas the time-stepping of the
(compressible) Lattice Boltzmann scheme is constrained to resolve sound waves.
Therefore, the time-step of a LB simulation is very small compared to the
time-step of a pseudo-spectral simulation, the ratio between the two time-
steps being typically the Mach number (Horstmann et al., 2022). The
superiority of LB simulations in terms of turn-around times is thus not so big
as compared to pseudo-spectral simulations. The things are different in the
Hall regime, where the time-steps of the two methods are constrained
identically by the speed of whistler waves. In that situation, the
computational efficiency of the LB scheme (exploiting the computational power
of GPU accelerators) is a major advantage and leads to wall-clock turn-around
times that are significantly smaller. Finally, let us mention that extension
to electron MHD is here straightforward and amounts to modifying the
equilibrium distributions for the magnetic field in the LB scheme, by
neglecting the bulk velocity $\mathbf{U}$ with respect to the Hall current
$\alpha_{H}\mathbf{J}$ in (38).

##  6 Hall-MHD simulations for space plasma turbulence investigations

Space plasmas, whose dynamics involve a vast range of scales due to their high
Reynolds numbers Matthaeus et al. (2008); Parashar et al. (2019), do actually
develop distinct MHD and Hall-MHD power-law spectral ranges, with a distinct
transition between them. This clearly emerges from the observations performed
with plasma and magnetic field instruments on board of state of the art
missions such as Solar Orbiter (SO; Müller et al., 2020b) and Parker Solar
Probe (PSP; Fox et al., 2016). Figure 17 shows the trace spectra of the
magnetic field fluctuations measured by the PSP/FIELDS (Bale et al., 2016) and
SO/MAG (Horbury et al., 2020) magnetometers on board the two spacecrafts. In
particular, the PSP (red) magnetic field sample, measured on November $20$,
$2021$, is relative to the fast solar wind plasma stream coming from an
equatorial coronal hole, while the SO (blue) sample is relative a low-speed
solar wind stream measured on July $14-15$, $2020$, whose origin was
identified in a coronal streamer and pseudostreamer configuration (D’Amicis et
al., 2021). In table 3 we report characteristic parameters relative to these
solar wind samples. Here it is worth to recall that the ion gyroradius
$\rho_{i}=v_{T,i}/\omega_{ci}$ (with $v_{T,i}$ ion thermal speed) and inertial
length $d_{i}=c/\omega_{pi}$ are defined in terms of ion cyclotron
$\omega_{ci}$ and plasma frequency $\omega_{pi}$, respectively (as indicated
in their formulas), the latter being in general significantly larger than the
former. For values of density $n$, temperature $T$ and $\beta$ typical of
space plasmas, the relation $\rho_{i}\lesssim d_{i}$ is valid. However, it was
observed by several authors how, for $\beta\sim\mathcal{O}(1)$, these two
characteristic length scales are comparable $\rho_{i}\simeq d_{i}$
(Alexandrova et al., 2008, 2009; Sahraoui et al., 2009; Kiyani et al., 2015).
Thus, in the solar wind the break point identifying the transition between the
MHD range and the kinetic range in the magnetic field spectrum, at the sub-ion
scales, is often referred as occurring either at the ion gyroradius or at the
inertial length scale, when $\beta\sim 1$. In spite of the different speeds,
both the SOLO and PSP solar wind samples we considered here are Alfvénic,
i.e., they are characterized by a high correlation between velocity and
magnetic field fluctuations (see Bruno & Carbone, 2013, and references
therein, for a comprehensive review on the solar wind turbulence). A clear
frequency break is observed at $k\rho_{i}\sim 1$ separating fluid and kinetic
scales, as shown in Fig.17, marking the transition from a MHD turbulent
inertial range (where energy is adiabatically transferred to smaller and
smaller scales), characterized Kolmogorov-like spectrum Marino et al. (2011);
Marino & Sorriso-Valvo (2023), to a range where kinetic effect begins to
dominate in which energy gets eventually dissipated, heating the solar wind
plasma Marino et al. (2008). As it is known, fluid and kinetic scales in the
solar wind are characterized by different power-law exponents in the spectra
based on spacecraft measurements. Features of their spectral ranges mostly
depend on the distance from the Sun at which observations are made, i.e., on
the observed stage of evolution of the solar wind turbulence (see, e.g.,
Telloni et al., 2021, 2022a). The physical phenomena as well as the governing
parameters that control the evolution of turbulence in the interplanetary
space are still matter of investigation. Nonlinear interactions (Bruno &
Carbone, 2013), expansion-driven velocity and magnetic (Shi et al., 2021) and
velocity shears Marino et al. (2012), as well as the parametric decay of
Alfvén waves (Malara & Velli, 1996), all certainly play some role. However, to
date, there is not a clear consensus on how turbulence evolves from a spectrum
resembling the one predicted by the Iroshnikov-Kraichnan phenomenology
(Iroshnikov, 1963; Kraichnan, 1965) to a Kolmogorov-like spectrum (Kolmogorov,
1941) as the solar wind expands from regions within the solar corona, or very
close to it, to the outer heliosphere. Moreover, the slope of the magnetic-
field spectrum beyond the ion skin depth (or ion inertial length) is highly
variable, with power-law exponents ranging from $\sim-4$ to $\sim-2$ (Smith et
al., 2006; Bruno et al., 2014), being also affected by the redistribution of
the magnetic field energy at the (larger) fluid scales: in general, the larger
is overall the power spectral density (PSD) within the MHD inertial range, the
steeper is the spectrum at the kinetic scales. A number of dissipative wave-
particle mechanisms are supposedly involved in the energy transfer process at
the very small scales. Among these, cyclotron-resonant dissipation certainly
plays an important role (see, e.g., Bruno & Trenchi, 2014; Telloni et al.,
2019), though the way energy is first brought to the small scales then
dissipated in the collision-less solar wind plasma is still a matter of
debate. Both the evolution of turbulence in the heliosphere and how energy is
dissipated within the solar wind, are two major open questions in the space
plasma community that could be effectively addressed with the numerical
simulation produced by flame, which allows to capture the transition between
MHD and Hall-MHD regimes. Another puzzle of solar and space plasma dynamics
that could be tackled using this new LB code is how magnetic switchbacks
observed in the solar corona as well as in the solar wind do contribute to the
heating of the plasma in the space. The switchbacks are intermittent magnetic-
field polarity reversals widely observed in the heliosphere (Bale et al.,
2019) and also in the solar corona (Telloni et al., 2022b). They are thought
to play a key role in the acceleration and heating of the solar wind. However,
disentangling their contribution to the energetics of the plasma from other
plasma processes is a challenging task, for which it is important to have
better quality (and highly resolved) $3$D Hall-MHD numerical simulations, able
to fully resolve a large dynamical range with largest possible scale
separation.

Probe | $\bar{\rho}\,[cm^{3}]$ | $\bar{V}\,[km/s]$ | $\bar{T}\,[MK]$ | $\bar{B}\,[n\mathrm{T}]$ | $\beta$ | $d_{i}\,[km]$ | $D_{sun}\,[\mathrm{AU}]$
---|---|---|---|---|---|---|---
PSP | 419 | 622 | 1.89 | 332 | 0.50 | 57.5 | 0.09
SOlO | 16 | 429 | 1.50 | 6.76 | 3.59 | 11.1 | 0.64
Table 3: Main solar wind parameters computed at time where the solar wind
samples used to produce the power spectra in Fig.017 have been collected. Here
we report solar wind density $\bar{\rho}$, solar wind bulk velocity $\bar{V}$,
proton temperature $\bar{T}$, average magnetic field $\bar{B}$, ion inertial
length $d_{i}$ and the distance of the spacecraft (SOLO and PSP) from the Sun
$D_{sun}$. Figure 17: Power Spectra Density (PSD) of the magnetic field
fluctuations observed by PSP (red) and SO (blue) on November $20$, $2021$ and
July $14-15$, $2020$, respectively. The $k^{-1.67}$, $k^{-1.5}$ and
$k^{-2.73}$ scalings are shown for reference as colored lines, as reported in
the legend.

##  7 Conclusions

The LB approach opens a new horizon for the numerical investigation of plasma
dynamics. Stability issues, which has long been a handicap for the application
of the LB method to turbulent flows, are now solved thanks to the use of
improved collision operators without notably sacrificing accuracy.
Furthermore, the computational efficiency of the LB schemes on many-core
devices such as GPUs allows for advantageous turn-around times. Another
significant advantage of dealing with a kinetic representation is that the
derivatives of the magnetic field are directly embedded in the solution, thus
allowing for an intrinsically accurate description of the current density
without resorting to any differentiation scheme. Altogether, our study shows
that the LB approach provides a very valuable and efficient means to simulate
Hall-MHD turbulence. Furthermore, the kinetic LB approach allows us to add
physical ingredients such as thermal effects, multi-species, complex
geometries and boundary conditions etc. at the cost of new coupled lattice
dynamics, therefore preserving the computational performance. The preliminary
results provided by a simple benchmark based on the OT vortex problem show
that our code is promising to simulate astrophysical and space plasmas where
the Hall term may play a significant role for the dynamics of the system.
Indeed, in plasmas as well as in anisotropic fluids, turbulence has to compete
with waves in transferring energy across the scales (Marino et al., 2015b).
The interplay of waves and turbulence is responsible for the emergence of new
characteristic length scales and the existence of different regimes (Marino et
al., 2013; Feraco et al., 2018) in which various forms of energy can cascade
to small or to large scales (Marino et al., 2014), or undergoing a dual energy
cascade (Marino et al., 2015a). The computational efficiency of our LB model
for the investigation of dynamics and energetics of anisotropic turbulent
flows will allow us to run simulations of fluids and plasmas separating
regimes (in terms of spatial and temporal scales) where different physical
phenomena dominate. We proved, even with a simple configuration, that our
kinetic model is able to capture various physical effects introduced by the
the Hall term, such as a faster dynamics due to the interplay of whistler
waves and turbulence, the breakdown of the Kolmogorov spectrum at sub-ion
scales and a behaviour of the magnetic energy spectrum at that scales which
has been already observed in solar wind measurements. All that provides our LB
approach with the potential to become a major tool of investigation for
magnetohydrodynamic plasmas in a variety of configurations of interest for
heliospheric and magnetospheric studies. The present results are produced by
using a benchmark configuration (OT vortex problem) which is isotropic, and
hence, does not reflect the strong anisotropy related to the ambient magnetic
field observed in the solar wind. However, we are confident that, with further
developments of our code, LB simulations of Hall-MHD flows will be suitable to
investigate rarefied heliospheric and magnetospheric plasmas at scales which
are nowadays within the reach of the high resolution instruments on board of
the latest solar wind and magnetospheric missions, such as Solar Orbiter or
Parker Solar Probe.

## Acknowledgements

We gratefully thank Dr. Emmanuel Quemener for his technical support in the
development of our multi-GPU code within the Centre Blaise Pascal’s computer
testing platform at ENS de Lyon (France). Most simulations were run on the
local HPC facilities at École Centrale de Lyon (PMCS2I) supported by the
Auvergne-Rhône-Alpes region (GRANT CPRT07-13 CIRA) and the national Equip@Meso
grant (ANR-10-EQPX-29-01). Finally, we acknowledge the CINECA computing center
under the ISCRA initiative, for the availability of HPC resources necessary
for high-resolution simulations of Hall-MHD turbulence, and support within the
project LaB-HMHD - HP10C4HXCB.

##  Appendix A Central-moment-based LB scheme for fluid dynamics

For the fluid, the discretization (in velocity) of the phase space refers to
the D3Q27 lattice. The set of adopted microscopic velocities
$\\{\mathbf{c}_{i}\\}_{i=0,..,26}$ is defined in Cartesian components by

$\begin{split}|c_{x}\rangle=[0,-1,0,0,-1,-1,-1,-1,0,0,-1,-1,-1,-1,1,0,0,1,1,1,1,0,0,1,1,1,1]^{\top}\\\
|c_{y}\rangle=[0,0,-1,0,-1,1,0,0,-1,-1,-1,-1,1,1,0,1,0,1,-1,0,0,1,1,1,1,-1,-1]^{\top}\\\
|c_{z}\rangle=[0,0,0,-1,0,0,-1,1,-1,1,-1,1,-1,1,0,0,1,0,0,1,-1,1,-1,1,-1,1,-1]^{\top}\end{split}$

The equilibrium densities (without accounting for the Lorentz force) are
developed up to the sixth-order as

$\begin{split}{f}_{i}^{(0)}(\rho,\mathbf{U})={}&w_{i}\rho\bigg{\\{}1+\frac{\mathbf{c}_{i}\cdot\mathbf{U}}{c_{s}^{2}}+\frac{1}{2c_{s}^{4}}\bigg{[}H^{(2)}_{ixx}U^{2}_{x}+H^{(2)}_{iyy}u^{2}_{y}+H^{(2)}_{izz}U^{2}_{z}+2\bigg{(}H^{(2)}_{ixy}U_{x}U_{y}+\\\
&+H^{(2)}_{ixz}U_{x}U_{z}+H^{(2)}_{iyz}U_{y}U_{z}\bigg{)}\bigg{]}+\frac{1}{2c_{s}^{6}}\bigg{[}H^{(3)}_{ixxy}U^{2}_{x}U_{y}+H^{(3)}_{ixxz}U^{2}_{x}U_{z}+H^{(3)}_{ixyy}U_{x}U^{2}_{y}+\\\
&+H^{(3)}_{ixzz}U_{x}U^{2}_{z}+H^{(3)}_{iyzz}U_{y}U^{2}_{z}+H^{(3)}_{iyyz}U^{2}_{y}U_{z}+2H^{(3)}_{ixyz}U_{x}U_{y}U_{z}\bigg{]}+\frac{1}{4c_{s}^{8}}\bigg{[}H^{(4)}_{ixxyy}U^{2}_{x}U^{2}_{y}+\\\
&+H^{(4)}_{ixxzz}U^{2}_{x}U^{2}_{z}+H^{(4)}_{iyyzz}U^{2}_{y}U^{2}_{z}+2\bigg{(}H^{(4)}_{ixyzz}U_{x}U_{y}U^{2}_{z}+H^{(4)}_{ixyyx}U_{x}U^{2}_{y}U_{z}+\\\
&+H^{(4)}_{ixxyz}U^{2}_{x}U_{y}U_{y}\bigg{)}\bigg{]}+\frac{1}{4c_{s}^{10}}\bigg{[}H^{(5)}_{ixxyzz}U^{2}_{x}U_{y}U^{2}_{z}+H^{(5)}_{ixxyyz}U^{2}_{x}U^{2}_{y}U_{z}+\\\
&+H^{(5)}_{ixyyzz}U_{x}U^{2}_{y}U^{2}_{z}\bigg{]}+\frac{1}{8c_{s}^{12}}H^{(6)}_{ixxyyzz}U^{2}_{x}U^{2}_{y}U^{2}_{z}\bigg{\\}}\end{split}$

where the weights $w_{i}$ are related to the lattice connectivity with
$w_{\mathrm{center}}=8/27$, $w_{\mathrm{face}}=2/27$, $w_{\mathrm{edge}}=8/27$
and $w_{\mathrm{corner}}=1/216$ for the D3Q27 lattice (see Fig. 1), and
$H^{(n)}_{i}$ refers to the nth-order Hermite polynomial tensor in velocity
$\mathbf{c}_{i}$. The Lorentz force is eventually taken into account by
upgrading the densities as

${f}_{i}^{\mathrm{mhd}(0)}(\rho,\mathbf{U},\mathbf{B})={f}_{i}^{(0)}(\rho,\mathbf{U})+\frac{w_{i}}{2c_{s}^{4}}\left((\mathbf{B}\cdot\mathbf{B})(\boldsymbol{c}_{i}\cdot\boldsymbol{c}_{i})-(\mathbf{c}_{i}\cdot\mathbf{B})^{2}\right).$

The set of central moments $\mathrm{k_{i}}$ is computed by applying the
(invertible) transformation matrix $\mathrm{T}$ with the column vectors

$\displaystyle|\mathrm{T}_{0}\rangle=|1\rangle$
$\displaystyle|\mathrm{T}_{1}\rangle;|\mathrm{T}_{2}\rangle;|\mathrm{T}_{3}\rangle=[\bar{c}_{ix}]^{\top};[\bar{c}_{iy}]^{\top};[\bar{c}_{iz}]^{\top}$
$\displaystyle|\mathrm{T}_{4}\rangle;|\mathrm{T}_{5}\rangle;|\mathrm{T}_{6}\rangle=[\bar{c}_{ix}\bar{c}_{iy}]^{\top};[\bar{c}_{ix}\bar{c}_{iz}]^{\top};[\bar{c}_{iy}\bar{c}_{iz}]^{\top}$
$\displaystyle|\mathrm{T}_{7}\rangle;|\mathrm{T}_{8}\rangle;|\mathrm{T}_{9}\rangle=[\bar{c}_{ix}^{2}-\bar{c}_{iy}^{2}]^{\top};[\bar{c}_{ix}^{2}-\bar{c}_{iz}^{2}]^{\top};[\bar{c}_{ix}^{2}+\bar{c}_{iy}^{2}+\bar{c}_{iz}^{2}]^{\top}$
$\displaystyle|\mathrm{T}_{10}\rangle;|\mathrm{T}_{11}\rangle;|\mathrm{T}_{12}\rangle=[\bar{c}_{ix}\bar{c}_{iy}^{2}+\bar{c}_{ix}\bar{c}_{iz}^{2}]^{\top};[\bar{c}_{ix}\bar{c}_{iy}^{2}+\bar{c}_{iy}\bar{c}_{iz}^{2}]^{\top};[\bar{c}_{ix}^{2}\bar{c}_{iy}+\bar{c}_{iy}^{2}\bar{c}_{iz}]^{\top}$
$\displaystyle|\mathrm{T}_{13}\rangle;|\mathrm{T}_{14}\rangle;|\mathrm{T}_{15}\rangle=[\bar{c}_{ix}\bar{c}_{iy}^{2}-\bar{c}_{ix}\bar{c}_{iz}^{2}]^{\top};[\bar{c}_{ix}\bar{c}_{iy}^{2}-\bar{c}_{iy}\bar{c}_{iz}^{2}]^{\top};[\bar{c}_{ix}^{2}\bar{c}_{iy}-\bar{c}_{iy}^{2}\bar{c}_{iz}]^{\top}$
$\displaystyle|\mathrm{T}_{16}\rangle=[\bar{c}_{ix}\bar{c}_{iy}\bar{c}_{iz}]^{\top}$
$\displaystyle|\mathrm{T}_{17}\rangle;|\mathrm{T}_{18}\rangle;|\mathrm{T}_{19}\rangle=[\bar{c}_{ix}^{2}\bar{c}_{iy}^{2}+\bar{c}_{ix}^{2}\bar{c}_{iz}^{2}+\bar{c}_{iy}^{2}\bar{c}_{iz}^{2}]^{\top};[\bar{c}_{ix}^{2}\bar{c}_{iy}^{2}+\bar{c}_{ix}^{2}\bar{c}_{iz}^{2}-\bar{c}_{iy}^{2}\bar{c}_{iz}^{2}]^{\top};[\bar{c}_{ix}^{2}\bar{c}_{iy}^{2}-\bar{c}_{ix}^{2}\bar{c}_{iz}^{2}]^{\top}$
$\displaystyle|\mathrm{T}_{20}\rangle;|\mathrm{T}_{21}\rangle;|\mathrm{T}_{22}\rangle=[\bar{c}_{ix}^{2}\bar{c}_{iy}\bar{c}_{iz}]^{\top};[\bar{c}_{ix}\bar{c}_{iy}^{2}\bar{c}_{iz}]^{\top};[\bar{c}_{ix}\bar{c}_{iy}\bar{c}_{iz}^{2}]^{\top}$
$\displaystyle|\mathrm{T}_{23}\rangle;|\mathrm{T}_{24}\rangle;|\mathrm{T}_{25}\rangle=[\bar{c}_{ix}\bar{c}_{iy}^{2}\bar{c}_{iz}^{2}]^{\top};[\bar{c}_{ix}^{2}\bar{c}_{iy}\bar{c}_{iz}^{2}]^{\top};[\bar{c}_{ix}^{2}\bar{c}_{iy}^{2}\bar{c}_{iz}]^{\top}$
$\displaystyle|\mathrm{T}_{26}\rangle=[\bar{c}_{ix}^{2}\bar{c}_{iy}^{2}\bar{c}_{iz}^{2}]^{\top}$

where $\bar{\mathbf{c}}_{i}=\mathbf{c}_{i}-\mathbf{U}$ is the set of
microscopic velocities obtained by the shift of particle velocities by the
local fluid velocity. The $27\times 27$ collision matrix $S$ for the central
moments is a diagonal matrix with the respective relaxation rates

$S=\mathrm{diag}[1,1,1,1,\omega,\omega,\omega,\omega,\omega,1,...,1],$

which leads to

$\displaystyle\mathrm{k}_{0\cdots 3}^{*}=\langle\mathrm{T}_{0\cdots
3}|{f^{\mathrm{mhd}(0)}}\rangle$ $\displaystyle\mathrm{k}_{4\cdots
8}^{*}=\omega\langle\mathrm{T}_{4\cdots
8}|{f^{\mathrm{mhd}(0)}}\rangle+(1-\omega)\langle\mathrm{T}_{4\cdots
8}|{\bar{f}}\rangle$ $\displaystyle\mathrm{k}_{9\cdots
26}^{*}=\langle\mathrm{T}_{9\cdots 26}|{f^{\mathrm{mhd}(0)}}\rangle.$

##  Appendix B Calculation of the electric current density

The electric current is obtained by solving the linear system (40). By using
(41), this system can be re-expressed as

$\widetilde{\mathbf{M}}\begin{bmatrix}J_{x}\\\\[3.00003pt]
J_{y}\\\\[3.00003pt]
J_{z}\end{bmatrix}=-\frac{1}{2\alpha_{H}}\begin{bmatrix}{\mathrm{\Lambda}}_{yz}-{\mathrm{\Lambda}}_{zy}-2\left(U_{y}B_{z}-U_{z}B_{y}\right)\\\\[3.00003pt]
{\mathrm{\Lambda}}_{zx}-{\mathrm{\Lambda}}_{xz}-2\left(U_{z}B_{x}-U_{x}B_{z}\right)\\\\[3.00003pt]
{\mathrm{\Lambda}}_{xy}-{\mathrm{\Lambda}}_{yx}-2\left(U_{x}B_{y}-U_{y}B_{x}\right)\end{bmatrix}$

where
${\mathrm{\Lambda}}_{\alpha\beta}=\sum_{i=0}^{M-1}\xi_{i\alpha}\bar{g}_{i\beta}$
and $\widetilde{\mathbf{M}}$ is the invertible matrix

$\widetilde{\mathbf{M}}=\begin{bmatrix}C^{2}/2\alpha_{H}\omega_{B}&B_{z}&-B_{y}\\\\[3.00003pt]
-B_{z}&C^{2}/2\alpha_{H}\omega_{B}&B_{x}\\\\[3.00003pt]
B_{y}&-B_{x}&C^{2}/2\alpha_{H}\omega_{B}\end{bmatrix}$

where $C$ represents the characteristic speed of _magnetic_ particles on the
D3Q7 lattice and $\omega_{B}$ is the relaxation pulsation (33) associated with
the BGK collision operator for the magnetic field. The expression for the
three components of the electric current density obtained by solving the
previous linear system reads as

$\begin{split}J_{x}={}&\frac{1}{D}(-2C^{4}\omega_{B}B_{y}U_{z}+2C^{4}\omega_{B}B_{z}U_{y}-C^{4}\omega_{B}\Lambda_{yz}+C^{4}\omega_{B}\Lambda_{zy}-4C^{2}\alpha_{H}\omega_{B}^{2}B_{x}B_{y}U_{y}-\\\
&-4C^{2}\alpha_{H}\omega_{B}^{2}B_{x}B_{z}U_{z}+4C^{2}\alpha_{H}\omega_{B}^{2}B_{y}^{2}U_{x}-2C^{2}\alpha_{H}\omega_{B}^{2}\Lambda_{xy}B_{y}+2C^{2}\alpha_{H}\omega_{B}^{2}\Lambda_{yx}B_{y}+\\\
&+4C^{2}\alpha_{H}\omega_{B}^{2}B_{z}^{2}U_{x}-2C^{2}\alpha_{H}\omega_{B}^{2}\Lambda_{xz}B_{z}+2C^{2}\alpha_{H}\omega_{B}^{2}\Lambda_{zx}B_{z}-4\alpha_{H}^{2}\omega_{B}^{3}\Lambda_{yz}B_{x}^{2}+\\\
&+4\alpha_{H}^{2}\omega_{B}^{3}\Lambda_{zy}B_{x}^{2}+4\alpha_{H}^{2}\omega_{B}^{3}\Lambda_{xz}B_{x}B_{y}-4\alpha_{H}^{2}\omega_{B}^{3}\Lambda_{zx}B_{x}B_{y}-4\alpha_{H}^{2}\omega_{B}^{3}\Lambda_{xy}B_{x}B_{z}+\\\
&+4\alpha_{H}^{2}\omega_{B}^{3}\Lambda_{yx}B_{x}B_{z})\end{split}$
$\begin{split}J_{y}={}&\frac{1}{D}(2C^{4}\omega_{B}B_{x}U_{z}+2C^{4}\omega_{B}B_{z}U_{x}+C^{4}\omega_{B}\Lambda_{xz}-C^{4}\omega_{B}\Lambda_{zx}+4C^{2}\alpha_{H}\omega_{B}^{2}B_{x}^{2}U_{y}-\\\
&-4C^{2}\alpha_{H}\omega_{B}^{2}B_{x}B_{y}U_{x}+2C^{2}\alpha_{H}\omega_{B}^{2}\Lambda_{xy}B_{x}-2C^{2}\alpha_{H}\omega_{B}^{2}\Lambda_{yx}B_{x}-4C^{2}\alpha_{H}\omega_{B}^{2}B_{y}B_{z}U_{z}+\\\
&+4C^{2}\alpha_{H}\omega_{B}^{2}B_{z}^{2}U_{y}-2C^{2}\alpha_{H}\omega_{B}^{2}\Lambda_{yz}B_{z}+2C^{2}\alpha_{H}\omega_{B}^{2}\Lambda_{zy}B_{z}-4\alpha_{H}^{2}\omega_{B}^{3}\Lambda_{yz}B_{x}B_{y}+\\\
&+4\alpha_{H}^{2}\omega_{B}^{3}\Lambda_{zy}B_{x}B_{y}+4\alpha_{H}^{2}\omega_{B}^{3}\Lambda_{xz}B_{y}^{2}-4\alpha_{H}^{2}\omega_{B}^{3}\Lambda_{zx}B_{y}^{2}-4\alpha_{H}^{2}\omega_{B}^{3}\Lambda_{xy}B_{y}B_{z}+\\\
&+4\alpha_{H}^{2}\omega_{B}^{3}\Lambda_{yx}B_{y}B_{z})\end{split}$
$\begin{split}J_{z}={}&\frac{1}{D}(-2C^{4}\omega_{B}B_{x}U_{y}+2C^{4}\omega_{B}B_{y}U_{x}-C^{4}\omega_{B}\Lambda_{xy}+C^{4}\omega_{B}\Lambda_{yx}+4C^{2}\alpha_{H}\omega_{B}^{2}B_{x}^{2}U_{z}-\\\
&-4C^{2}\alpha_{H}\omega_{B}^{2}B_{x}B_{z}U_{x}+2C^{2}\alpha_{H}\omega_{B}^{2}\Lambda_{xz}B_{x}-2C^{2}\alpha_{H}\omega_{B}^{2}\Lambda_{zx}B_{x}+4C^{2}\alpha_{H}\omega_{B}^{2}B_{y}^{2}U_{z}-\\\
&-4C^{2}\alpha_{H}\omega_{B}^{2}B_{y}B_{z}U_{y}+2C^{2}\alpha_{H}\omega_{B}^{2}\Lambda_{yz}B_{y}-2C^{2}\alpha_{H}\omega_{B}^{2}\Lambda_{zy}B_{y}-4\alpha_{H}^{2}\omega_{B}^{3}\Lambda_{yz}B_{x}B_{z}+\\\
&+4\alpha_{H}^{2}\omega_{B}^{3}\Lambda_{zy}B_{x}B_{y}+4\alpha_{H}^{2}\omega_{B}^{3}\Lambda_{xz}B_{y}B_{z}-4\alpha_{H}^{2}\omega_{B}^{3}\Lambda_{zx}B_{y}B_{z}-4\alpha_{H}^{2}\omega_{B}^{3}\Lambda_{xy}B_{z}^{2}+\\\
&+4\alpha_{H}^{2}\omega_{B}^{3}\Lambda_{yx}B_{z}^{2})\end{split}$ (70)

where $D=C^{6}+4C^{2}\alpha_{H}^{2}\omega_{B}^{2}|\mathbf{B}|^{2}>0$.

## References

  * Alexandrova et al. (2008) Alexandrova, O., Lacombe, C. & Mangeney, A. 2008 Spectra and anisotropy of magnetic fluctuations in the earth’s magnetosheath: Cluster observations. Annales Geophysicae 26 (11), 3585–3596.
  * Alexandrova et al. (2009) Alexandrova, O., Saur, J., Lacombe, C., Mangeney, A., Mitchell, J., Schwartz, S. J. & Robert, P. 2009 Universality of solar-wind turbulent spectrum from mhd to electron scales. Phys. Rev. Lett. 103, 165003.
  * Bale et al. (2019) Bale, S. D., Badman, S. T., Bonnell, J. W., Bowen, T. A., Burgess, D., Case, A. W., Cattell, C. A., Chandran, B. D. G., Chaston, C. C., Chen, C. H. K., Drake, J. F., de Wit, T. Dudok, Eastwood, J. P., Ergun, R. E., Farrell, W. M., Fong, C., Goetz, K., Goldstein, M., Goodrich, K. A., Harvey, P. R., Horbury, T. S., Howes, G. G., Kasper, J. C., Kellogg, P. J., Klimchuk, J. A., Korreck, K. E., Krasnoselskikh, V. V., Krucker, S., Laker, R., Larson, D. E., MacDowall, R. J., Maksimovic, M., Malaspina, D. M., Martinez-Oliveros, J., McComas, D. J., Meyer-Vernet, N., Moncuquet, M., Mozer, F. S., Phan, T. D., Pulupa, M., Raouafi, N. E., Salem, C., Stansby, D., Stevens, M., Szabo, A., Velli, M., Woolley, T. & Wygant, J. R. 2019 Highly structured slow solar wind emerging from an equatorial coronal hole. Nature 576 (7786), 237–242.
  * Bale et al. (2016) Bale, S. D., Goetz, K., Harvey, P. R., Turin, P., Bonnell, J. W., Dudok de Wit, T., Ergun, R. E., MacDowall, R. J., Pulupa, M., Andre, M., Bolton, M., Bougeret, J. L., Bowen, T. A., Burgess, D., Cattell, C. A., Chandran, B. D. G., Chaston, C. C., Chen, C. H. K., Choi, M. K., Connerney, J. E., Cranmer, S., Diaz-Aguado, M., Donakowski, W., Drake, J. F., Farrell, W. M., Fergeau, P., Fermin, J., Fischer, J., Fox, N., Glaser, D., Goldstein, M., Gordon, D., Hanson, E., Harris, S. E., Hayes, L. M., Hinze, J. J., Hollweg, J. V., Horbury, T. S., Howard, R. A., Hoxie, V., Jannet, G., Karlsson, M., Kasper, J. C., Kellogg, P. J., Kien, M., Klimchuk, J. A., Krasnoselskikh, V. V., Krucker, S., Lynch, J. J., Maksimovic, M., Malaspina, D. M., Marker, S., Martin, P., Martinez-Oliveros, J., McCauley, J., McComas, D. J., McDonald, T., Meyer-Vernet, N., Moncuquet, M., Monson, S. J., Mozer, F. S., Murphy, S. D., Odom, J., Oliverson, R., Olson, J., Parker, E. N., Pankow, D., Phan, T., Quataert, E., Quinn, T., Ruplin, S. W., Salem, C., Seitz, D., Sheppard, D. A., Siy, A., Stevens, K., Summers, D., Szabo, A., Timofeeva, M., Vaivads, A., Velli, M., Yehle, A., Werthimer, D. & Wygant, J. R. 2016 The FIELDS Instrument Suite for Solar Probe Plus. Measuring the Coronal Plasma and Magnetic Field, Plasma Waves and Turbulence, and Radio Signatures of Solar Transients. Space Sci. Rev. 204 (1-4), 49–82.
  * Bhatnagar et al. (1954) Bhatnagar, P. L., Gross, E. P. & Krook, M. 1954 A model for collision processes in gases. i. small amplitude processes in charged and neutral one-component systems. Phys. Rev. 94, 511–525.
  * Bouchut (1999) Bouchut, François 1999 Construction of bgk models with a family of kinetic entropies for a given system of conservation laws. J. Stat. Phys. 95, 113–170.
  * Bruno & Carbone (2013) Bruno, Roberto & Carbone, Vincenzo 2013 The Solar Wind as a Turbulence Laboratory. Living Rev. Sol. Phys. 10 (1), 2.
  * Bruno & Trenchi (2014) Bruno, R. & Trenchi, L. 2014 Radial Dependence of the Frequency Break between Fluid and Kinetic Scales in the Solar Wind Fluctuations. Astrophys. J., Lett. 787 (2), L24.
  * Bruno et al. (2014) Bruno, R., Trenchi, L. & Telloni, D. 2014 Spectral Slope Variation at Proton Scales from Fast to Slow Solar Wind. Astrophys. J., Lett. 793 (1), L15.
  * Chen et al. (1991) Chen, Shiyi, Chen, Hudong, Martnez, Daniel & Matthaeus, William 1991 Lattice boltzmann model for simulation of magnetohydrodynamics. Phys. Rev. Lett. 67, 3776–3779.
  * Coreixas et al. (2019) Coreixas, Christophe, Chopard, Bastien & Latt, Jonas 2019 Comprehensive comparison of collision models in the lattice boltzmann framework: Theoretical investigations. Phys. Rev. E 100, 033305.
  * Coreixas et al. (2017) Coreixas, Christophe, Wissocq, Gauthier, Puigt, Guillaume, Boussuge, Jean-François & Sagaut, Pierre 2017 Recursive regularization step for high-order lattice boltzmann methods. Phys. Rev. E 96 (3), 033306.
  * Croisille et al. (1995) Croisille, JP., Khanfir, R. & Chanteu, G. 1995 Numerical simulation of the mhd equations by a kinetic-type method. J. Sci. Comput. 10, 81–92.
  * D’Amicis et al. (2021) D’Amicis, R., Bruno, R., Panasenco, O., Telloni, D., Perrone, D., Marcucci, M. F., Woodham, L., Velli, M., De Marco, R., Jagarlamudi, V., Coco, I., Owen, C., Louarn, P., Livi, S., Horbury, T., André, N., Angelini, V., Evans, V., Fedorov, A., Genot, V., Lavraud, B., Matteini, L., Müller, D., O’Brien, H., Pezzi, O., Rouillard, A. P., Sorriso-Valvo, L., Tenerani, A., Verscharen, D. & Zouganelis, I. 2021 First Solar Orbiter observation of the Alfvénic slow wind and identification of its solar source. Astron. Astrophys. 656, A21.
  * De Rosis (2017) De Rosis, Alessandro 2017 Nonorthogonal central-moments-based lattice boltzmann scheme in three dimensions. Phys. Rev. E 95, 013310.
  * De Rosis & Luo (2019) De Rosis, Alessandro & Luo, Kai H. 2019 Role of higher-order hermite polynomials in the central-moments-based lattice boltzmann framework. Phys. Rev. E 99 (1), 013301.
  * De Rosis et al. (2018) De Rosis, Alessandro, Lévêque, Emmanuel & Chahine, Robert 2018 Advanced lattice boltzmann scheme for high-reynolds-number magneto-hydrodynamic flows. J. Turb. 19 (6), 446–462.
  * Dellar (2002) Dellar, Paul J. 2002 Lattice kinetic schemes for magnetohydrodynamics. Journal of Computational Physics 179 (1), 95–126.
  * Feraco et al. (2021) Feraco, F., Marino, R., Primavera, L., Pumir, A., Mininni, P. D., Rosenberg, D., Pouquet, A., Foldes, R., Lévêque, E., Camporeale, E., Cerri, S. S., Charuvil Asokan, H., Chau, J. L., Bertoglio, J. P., Salizzoni, P. & Marro, M. 2021 Connecting large-scale velocity and temperature bursts with small-scale intermittency in stratified turbulence. EPL (Europhysics Letters) 135 (1), 14001.
  * Feraco et al. (2018) Feraco, Fabio, Marino, Raffaele, Pumir, Alain, Primavera, Leonardo, Mininni, Pablo, Pouquet, Annick & Rosenberg, Duane 2018 Vertical drafts and mixing in stratified turbulence: Sharp transition with froude number. EPL (Europhysics Letters) 123, 44002.
  * Ferrand et al. (2022) Ferrand, R., Sahraoui, F., Galtier, S., Andrés, N., Mininni, P. & Dmitruk, P. 2022 An in-depth numerical study of exact laws for compressible hall magnetohydrodynamic turbulence. The Astrophysical Journal 927 (2), 205.
  * Flint & Vahala (2018) Flint, Christopher & Vahala, George 2018 A partial entropic lattice boltzmann mhd simulation of the orszag–tang vortex. Radiation Effects and Defects in Solids 173 (1-2), 55–65, arXiv: https://doi.org/10.1080/10420150.2018.1442460.
  * Fox et al. (2016) Fox, N. J., Velli, M. C., Bale, S. D., Decker, R., Driesman, A., Howard, R. A., Kasper, J. C., Kinnison, J., Kusterer, M., Lario, D., Lockwood, M. K., McComas, D. J., Raouafi, N. E. & Szabo, A. 2016 The Solar Probe Plus Mission: Humanity’s First Visit to Our Star. Space Sci. Rev. 204 (1-4), 7–48.
  * Galtier (2016) Galtier, Sébastien 2016 Introduction to Modern Magnetohydrodynamics. Cambridge University Press.
  * Galtier & Buchlin (2007) Galtier, Sebastien & Buchlin, Eric 2007 Multiscale hall-magnetohydrodynamic turbulence in the solar wind. Astrophysical Journal 656 (1), 560–566.
  * Geier et al. (2007) Geier, M., Greiner, A. & Korvink, J.G. 2007 Properties of the cascaded lattice Boltzmann automaton. Int. J. Mod. Phys. C 18 (04), 455–462.
  * Geier et al. (2015) Geier, Martin, Schönherr, Martin, Pasquali, Andrea & Krafczyk, Manfred 2015 The cumulant lattice boltzmann equation in three dimensions: Theory and validation. Computers & Mathematics with Applications 70 (4), 507–547.
  * Gómez et al. (2010) Gómez, Daniel O., Mininni, Pablo D. & Dmitruk, Pablo 2010 Hall-magnetohydrodynamic small-scale dynamos. Phys. Rev. E 82, 036406.
  * González-Morales et al. (2019) González-Morales, P. A., Khomenko, E. & Cally, P. S. 2019 Fast-to-alfvén mode conversion mediated by hall current. II. application to the solar atmosphere. Astrophysical Journal 870 (2), 94.
  * He & Luo (1997) He, Xiaoyi & Luo, Li-Shi 1997 Theory of the lattice boltzmann method: From the boltzmann equation to the lattice boltzmann equation. Phys. Rev. E 56, 6811–6817.
  * Horbury et al. (2020) Horbury, T. S., O’Brien, H., Carrasco Blazquez, I., Bendyk, M., Brown, P., Hudson, R., Evans, V., Oddy, T. M., Carr, C. M., Beek, T. J., Cupido, E., Bhattacharya, S., Dominguez, J. A., Matthews, L., Myklebust, V. R., Whiteside, B., Bale, S. D., Baumjohann, W., Burgess, D., Carbone, V., Cargill, P., Eastwood, J., Erdös, G., Fletcher, L., Forsyth, R., Giacalone, J., Glassmeier, K. H., Goldstein, M. L., Hoeksema, T., Lockwood, M., Magnes, W., Maksimovic, M., Marsch, E., Matthaeus, W. H., Murphy, N., Nakariakov, V. M., Owen, C. J., Owens, M., Rodriguez-Pacheco, J., Richter, I., Riley, P., Russell, C. T., Schwartz, S., Vainio, R., Velli, M., Vennerstrom, S., Walsh, R., Wimmer-Schweingruber, R. F., Zank, G., Müller, D., Zouganelis, I. & Walsh, A. P. 2020 The Solar Orbiter magnetometer. Astron. Astrophys. 642, A9.
  * Horstmann et al. (2022) Horstmann, Jan, Touil, Hatem, Vienne, Lucien, Ricot, Denis & Lévêque, Emmanuel 2022 Consistent time-step optimization in the lattice boltzmann method. Journal of Computational Physics 462, 111224.
  * Huba (2003) Huba, J. D. 2003 Hall Magnetohydrodynamics - A Tutorial. In Space Plasma Simulation (ed. J. Büchner, C. Dum & M. Scholer), , vol. 615, pp. 166–192.
  * Iroshnikov (1963) Iroshnikov, P. S. 1963 Turbulence of a Conducting Fluid in a Strong Magnetic Field. Astron. Zh. 40, 742.
  * Kiyani et al. (2015) Kiyani, Khurom H., Osman, Kareem T. & Chapman, Sandra C. 2015 Dissipation and heating in solar wind turbulence: from the macro to the micro and back again. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 373 (2041), 20140155.
  * Kolmogorov (1941) Kolmogorov, A. 1941 The Local Structure of Turbulence in Incompressible Viscous Fluid for Very Large Reynolds’ Numbers. Dokl. Akad. Nauk SSSR 30, 301–305.
  * Körner et al. (2006) Körner, Carolin, Pohl, Thomas, Rüde, Ulrich, Thürey, Nils & Zeiser, Thomas 2006 Parallel lattice boltzmann methods for cfd applications. In Numerical Solution of Partial Differential Equations on Parallel Computers (ed. Are Magnus Bruaset & Aslak Tveito), pp. 439–466. Berlin, Heidelberg: Springer Berlin Heidelberg.
  * Kraichnan (1965) Kraichnan, Robert H. 1965 Inertial-Range Spectrum of Hydromagnetic Turbulence. Phys. Fluids 8 (7), 1385–1387.
  * Krueger et al. (2016) Krueger, Timm, Kusumaatmaja, Halim, Kuzmin, Alexandr, Shardt, Orest, Silva, Goncalo & Viggen, Erlend Magnus 2016 The Lattice Boltzmann Method: Principles and Practice. Springer.
  * Lewy et al. (1928) Lewy, H, Friedrichs, K & Courant, R 1928 Über die partiellen differenzengleichungen der mathematischen physik. Mathematische annalen 100, 32–74.
  * Ma et al. (2018) Ma, Yingjuan, Russell, Christopher T., Toth, Gabor, Chen, Yuxi, Nagy, Andrew F., Harada, Yuki, McFadden, James, Halekas, Jasper S., Lillis, Rob, Connerney, John E. P., Espley, Jared, DiBraccio, Gina A., Markidis, Stefano, Peng, Ivy Bo, Fang, Xiaohua & Jakosky, Bruce M. 2018 Reconnection in the martian magnetotail: Hall-mhd with embedded particle-in-cell simulations. Journal of Geophysical Research: Space Physics 123 (5), 3742–3763.
  * Mahajan & Krishan (2005) Mahajan, S. M. & Krishan, V. 2005 Exact solution of the incompressible Hall magnetohydrodynamics. Monthly Notices of the Royal Astronomical Society: Letters 359 (1), L27–L29.
  * Malara & Velli (1996) Malara, F. & Velli, M. 1996 Parametric instability of a large-amplitude nonmonochromatic Alfvén wave. Phys. Plasmas 3 (12), 4427–4433.
  * Malaspinas (2015) Malaspinas, Orestis 2015 Increasing stability and accuracy of the lattice boltzmann scheme: recursivity and regularization. arXiv preprint arXiv:1505.06900 .
  * Marchand, P. et al. (2018) Marchand, P., Commerçon, B. & Chabrier, G. 2018 Impact of the hall effect in star formation and the issue of angular momentum conservation. A&A 619, A37.
  * Marino et al. (2013) Marino, R., Mininni, P. D., Rosenberg, D. & Pouquet, A. 2013 Inverse cascades in rotating stratified turbulence: Fast growth of large scales. EPL (Europhysics Letters) 102 (4), 44006.
  * Marino et al. (2014) Marino, R., Mininni, P. D., Rosenberg, D. L. & Pouquet, A. 2014 Large-scale anisotropy in stably stratified rotating flows. Phys. Rev. E 90, 023018.
  * Marino et al. (2015a) Marino, R., Pouquet, A. & Rosenberg, D. 2015a Resolving the paradox of oceanic large-scale balance and small-scale mixing. Phys. Rev. Lett. 114, 114504.
  * Marino et al. (2015b) Marino, Raffaele, Rosenberg, Duane, Herbert, Corentin & Pouquet, Annick 2015b Interplay of waves and eddies in rotating stratified turbulence and the link with kinetic-potential energyartition. EPL (Europhysics Letters) 112.
  * Marino & Sorriso-Valvo (2023) Marino, R. & Sorriso-Valvo, L. 2023 Scaling laws for the energy transfer in space plasma turbulence. Physics Reports (IN PRESS) .
  * Marino et al. (2008) Marino, R., Sorriso-Valvo, L., Carbone, V., Noullez, A., Bruno, R. & Bavassano, B. 2008 Heating the Solar Wind by a Magnetohydrodynamic Turbulent Energy Cascade. Astroph. Journal Lett. 677 (1), L71.
  * Marino et al. (2011) Marino, Raffaele, Sorriso-Valvo, Luca, Carbone, Vincenzo, Veltri, Pierluigi, Noullez, Alain & Bruno, Roberto 2011 The magnetohydrodynamic turbulent cascade in the ecliptic solar wind: Study of Ulysses data. Planetary and Space Science 59 (7), 592–597.
  * Marino et al. (2012) Marino, R., Sorriso-Valvo, L., D’Amicis, R., Carbone, V., Bruno, R. & Veltri, P. 2012 On the Occurrence of the Third-order Scaling in High Latitude Solar Wind. Astroph. Journal 750 (1), 41.
  * Martínez et al. (1994) Martínez, Daniel O., Chen, Shiyi & Matthaeus, William H. 1994 Lattice boltzmann magnetohydrodynamics. Physics of Plasmas 1 (6), 1850–1867, arXiv: https://doi.org/10.1063/1.870640.
  * Matthaeus et al. (2008) Matthaeus, W. H., Weygand, J. M., Chuychai, P., Dasso, S., Smith, C. W. & Kivelson, M. G. 2008 Interplanetary Magnetic Taylor Microscale and Implications for Plasma Dissipation. Astroph. Journal Lett. 678 (2), L141.
  * Meyrand & Galtier (2012) Meyrand, Romain & Galtier, Sébastien 2012 Spontaneous chiral symmetry breaking of hall magnetohydrodynamic turbulence. Phys. Rev. Lett. 109, 194501.
  * Mininni et al. (2002) Mininni, Pablo D., Gómez, Daniel O. & Mahajan, Swadesh M. 2002 Dynamo action in hall magnetohydrodynamics. Astrophysical Journal 567 (1), L81–L83.
  * Mininni et al. (2003) Mininni, Pablo D., Gomez, Daniel O. & Mahajan, Swadesh M. 2003 Dynamo action in magnetohydrodynamics and hall-magnetohydrodynamics. Astrophysical Journal 587 (1), 472–481.
  * Mininni et al. (2005) Mininni, Pablo D., Gomez, Daniel O. & Mahajan, Swadesh M. 2005 Direct simulations of helical hall-MHD turbulence and dynamo action. Astrophysical Journal 619 (2), 1019–1027.
  * Mininni et al. (2006) Mininni, P. D., Pouquet, A. G. & Montgomery, D. C. 2006 Small-scale structures in three-dimensional magnetohydrodynamic turbulence. Phys. Rev. Lett. 97, 244503\.
  * Miura & Araki (2014) Miura, H. & Araki, K. 2014 Structure transitions induced by the hall term in homogeneous and isotropic magnetohydrodynamic turbulence. Physics of Plasmas 21 (7), 072313.
  * Montgomery & Doolen (1987) Montgomery, David & Doolen, Gary D. 1987 Magnetohydrodynamic cellular automata. Physics Letters A 120, 229–231.
  * Morales et al. (2005) Morales, Laura, Dasso, S. & Gómez, Daniel 2005 Hall effect in incompressible magnetic reconnection. Journal of Geophysical Research 110.
  * Müller et al. (2020a) Müller, D., St. Cyr, O. C., Zouganelis, I., Gilbert, H. R., Marsden, R., Nieves-Chinchilla, T., Antonucci, E., Auchère, F., Berghmans, D., Horbury, T. S., Howard, R. A., Krucker, S., Maksimovic, M., Owen, C. J., Rochus, P., Rodriguez-Pacheco, J., Romoli, M., Solanki, S. K., Bruno, R., Carlsson, M., Fludra, A., Harra, L., Hassler, D. M., Livi, S., Louarn, P., Peter, H., Schühle, U., Teriaca, L., del Toro Iniesta, J. C., Wimmer-Schweingruber, R. F., Marsch, E., Velli, M., De Groof, A., Walsh, A. & Williams, D. 2020a The Solar Orbiter mission. Science overview. Astronomy & Astrophysics 642, A1.
  * Müller et al. (2020b) Müller, D., St. Cyr, O. C., Zouganelis, I., Gilbert, H. R., Marsden, R., Nieves-Chinchilla, T., Antonucci, E., Auchère, F., Berghmans, D., Horbury, T. S., Howard, R. A., Krucker, S., Maksimovic, M., Owen, C. J., Rochus, P., Rodriguez-Pacheco, J., Romoli, M., Solanki, S. K., Bruno, R., Carlsson, M., Fludra, A., Harra, L., Hassler, D. M., Livi, S., Louarn, P., Peter, H., Schühle, U., Teriaca, L., del Toro Iniesta, J. C., Wimmer-Schweingruber, R. F., Marsch, E., Velli, M., De Groof, A., Walsh, A. & Williams, D. 2020b The Solar Orbiter mission. Science overview. Astron. Astrophys. 642, A1.
  * Norman & Heyvaerts (1985) Norman, C. & Heyvaerts, J. 1985 Anomalous magnetic field diffusion during star formation. Astronomy & Astrophysics 147 (2), 247–256.
  * Orszag & Tang (1979) Orszag, S. A. & Tang, C. M. 1979 Small-scale structure of two-dimensional magnetohydrodynamic turbulence. Journal of Fluid Mechanics 90, 129–143.
  * Pandey & Wardle (2008) Pandey, B. P. & Wardle, Mark 2008 Hall magnetohydrodynamics of partially ionized plasmas. Monthly Notices of the Royal Astronomical Society 385 (4), 2269–2278, arXiv: https://onlinelibrary.wiley.com/doi/pdf/10.1111/j.1365-2966.2008.12998.x.
  * Papini et al. (2019) Papini, Emanuele, Franci, Luca, Landi, Simone, Verdini, Andrea, Matteini, Lorenzo & Hellinger, Petr 2019 Can hall magnetohydrodynamics explain plasma turbulence at sub-ion scales? Astrophysical Journal 870 (1), 52.
  * Parashar et al. (2019) Parashar, T. N., Cuesta, M. & Matthaeus, W. H. 2019 Reynolds Number and Intermittency in the Expanding Solar Wind: Predictions Based on Voyager Observations. Astroph. Journal Lett. 884 (2), L57.
  * Patterson & Orszag (1971) Patterson, G. S. & Orszag, Steven A. 1971 Spectral calculations of isotropic turbulence: Efficient removal of aliasing interactions. The Physics of Fluids 14 (11), 2538–2541.
  * Rosenberg et al. (2020) Rosenberg, Duane, Mininni, Pablo D., Reddy, Raghu & Pouquet, Annick 2020 Gpu parallelization of a hybrid pseudospectral geophysical turbulence framework using cuda. Atmosphere 11 (2).
  * Sahraoui et al. (2009) Sahraoui, F., Goldstein, M. L., Robert, P. & Khotyaintsev, Yu. V. 2009 Evidence of a cascade and dissipation of solar-wind turbulence at the electron gyroscale. Phys. Rev. Lett. 102, 231102.
  * Shan & He (1998) Shan, Xiaowen & He, Xiaoyi 1998 Discretization of the velocity space in the solution of the boltzmann equation. Phys. Rev. Lett. 80, 65–68.
  * Shi et al. (2021) Shi, C., Velli, M., Panasenco, O., Tenerani, A., Réville, V., Bale, S. D., Kasper, J., Korreck, K., Bonnell, J. W., Dudok de Wit, T., Malaspina, D. M., Goetz, K., Harvey, P. R., MacDowall, R. J., Pulupa, M., Case, A. W., Larson, D., Verniero, J. L., Livi, R., Stevens, M., Whittlesey, P., Maksimovic, M. & Moncuquet, M. 2021 Alfvénic versus non-Alfvénic turbulence in the inner heliosphere as observed by Parker Solar Probe. Astron. Astrophys. 650, A21.
  * Silva & Semiao (2014) Silva, Goncalo & Semiao, Viriato 2014 Truncation errors and the rotational invariance of three-dimensional lattice models in the lattice boltzmann method. Journal of Computational Physics 269, 259–279.
  * Smith et al. (2006) Smith, Charles W., Hamilton, Kathleen, Vasquez, Bernard J. & Leamon, Robert J. 2006 Dependence of the Dissipation Range Spectrum of Interplanetary Magnetic Fluctuationson the Rate of Energy Cascade. Astrophys. J., Lett. 645 (1), L85–L88.
  * Succi et al. (1991) Succi, S., Vergassola, M. & Benzi, R. 1991 Lattice boltzmann scheme for two-dimensional magnetohydrodynamics. Phys. Rev. A 43, 4521–4524.
  * Telloni et al. (2022a) Telloni, Daniele, Adhikari, Laxman, Zank, Gary P., Hadid, Lina Z., Sánchez-Cano, Beatriz, Sorriso-Valvo, Luca, Zhao, Lingling, Panasenco, Olga, Shi, Chen, Velli, Marco, Susino, Roberto, Verscharen, Daniel, Milillo, Anna, Alberti, Tommaso, Narita, Yasuhito, Verdini, Andrea, Grimani, Catia, Bruno, Roberto, D’Amicis, Raffaella, Perrone, Denise, Marino, Raffaele, Carbone, Francesco, Califano, Francesco, Malara, Francesco, Stawarz, Julia E., Laker, Ronan, Liberatore, Alessandro, Bale, Stuart D., Kasper, Justin C., Heyner, Daniel, de Wit, Thierry Dudok, Goetz, Keith, Harvey, Peter R., MacDowall, Robert J., Malaspina, David M., Pulupa, Marc, Case, Anthony W., Korreck, Kelly E., Larson, Davin, Livi, Roberto, Stevens, Michael L., Whittlesey, Phyllis, Auster, Hans-Ulrich & Richter, Ingo 2022a Observation and Modeling of the Solar Wind Turbulence Evolution in the Sub-Mercury Inner Heliosphere. Astrophys. J., Lett. 938 (2), L8.
  * Telloni et al. (2019) Telloni, Daniele, Carbone, Francesco, Bruno, Roberto, Zank, Gary P., Sorriso-Valvo, Luca & Mancuso, Salvatore 2019 Ion Cyclotron Waves in Field-aligned Solar Wind Turbulence. Astrophys. J., Lett. 885 (1), L5.
  * Telloni et al. (2021) Telloni, Daniele, Sorriso-Valvo, Luca, Woodham, Lloyd D., Panasenco, Olga, Velli, Marco, Carbone, Francesco, Zank, Gary P., Bruno, Roberto, Perrone, Denise, Nakanotani, Masaru, Shi, Chen, D’Amicis, Raffaella, De Marco, Rossana, Jagarlamudi, Vamsee K., Steinvall, Konrad, Marino, Raffaele, Adhikari, Laxman, Zhao, Lingling, Liang, Haoming, Tenerani, Anna, Laker, Ronan, Horbury, Timothy S., Bale, Stuart D., Pulupa, Marc, Malaspina, David M., MacDowall, Robert J., Goetz, Keith, de Wit, Thierry Dudok, Harvey, Peter R., Kasper, Justin C., Korreck, Kelly E., Larson, Davin, Case, Anthony W., Stevens, Michael L., Whittlesey, Phyllis, Livi, Roberto, Owen, Christopher J., Livi, Stefano, Louarn, Philippe, Antonucci, Ester, Romoli, Marco, O’Brien, Helen, Evans, Vincent & Angelini, Virginia 2021 Evolution of Solar Wind Turbulence from 0.1 to 1 au during the First Parker Solar Probe-Solar Orbiter Radial Alignment. Astrophys. J., Lett. 912 (2), L21.
  * Telloni et al. (2022b) Telloni, Daniele, Zank, Gary P., Stangalini, Marco, Downs, Cooper, Liang, Haoming, Nakanotani, Masaru, Andretta, Vincenzo, Antonucci, Ester, Sorriso-Valvo, Luca, Adhikari, Laxman, Zhao, Lingling, Marino, Raffaele, Susino, Roberto, Grimani, Catia, Fabi, Michele, D’Amicis, Raffaella, Perrone, Denise, Bruno, Roberto, Carbone, Francesco, Mancuso, Salvatore, Romoli, Marco, Deppo, Vania Da, Fineschi, Silvano, Heinzel, Petr, Moses, John D., Naletto, Giampiero, Nicolini, Gianalfredo, Spadaro, Daniele, Teriaca, Luca, Frassati, Federica, Jerse, Giovanna, Landini, Federico, Pancrazzi, Maurizio, Russano, Giuliana, Sasso, Clementina, Biondo, Ruggero, Burtovoi, Aleksandr, Capuano, Giuseppe E., Casini, Chiara, Casti, Marta, Chioetto, Paolo, Leo, Yara De, Giarrusso, Marina, Liberatore, Alessandro, Berghmans, David, Auchère, Frédéric, Cuadrado, Regina Aznar, Chitta, Lakshmi P., Harra, Louise, Kraaikamp, Emil, Long, David M., Mandal, Sudip, Parenti, Susanna, Pelouze, Gabriel, Peter, Hardi, Rodriguez, Luciano, Schühle, Udo, Schwanitz, Conrad, Smith, Phil J., Verbeeck, Cis & Zhukov, Andrei N. 2022b Observation of a Magnetic Switchback in the Solar Corona. Astrophys. J., Lett. 936 (2), L25.
  * Tóth et al. (2008) Tóth, Gábor, Ma, Yingjuan & Gombosi, Tamas I. 2008 Hall magnetohydrodynamics on block-adaptive grids. Journal of Computational Physics 227 (14), 6967–6984.
  * Wang et al. (2001) Wang, Xiaogang, Bhattacharjee, A. & Ma, Z. W. 2001 Scaling of collisionless forced reconnection. Phys. Rev. Lett. 87, 265003.
  * Xia & Yang (2015) Xia, Zhenwei & Yang, Weihong 2015 Exact solutions of the incompressible dissipative hall magnetohydrodynamics. Physics of Plasmas 22 (3), 032306, arXiv: https://doi.org/10.1063/1.4914931.
  * Yadav et al. (2022) Yadav, Sharad K., Miura, Hideaki & Pandit, Rahul 2022 Statistical properties of three-dimensional Hall magnetohydrodynamics turbulence. Physics of Fluids 34 (9), 095135\.