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

Raffaello Foldes, Emmanuel Lévêque, Raffaele Marino, Ermanno Pietropaolo, Alessandro de Rosis, Daniele Telloni, Fabio Feraco

► To cite this version:

Raffaello Foldes, Emmanuel Lévêque, Raffaele Marino, Ermanno Pietropaolo, Alessandro de Rosis, et al.. Efficient kinetic Lattice Boltzmann simulation of three-dimensional Hall-MHD Turbulence. Journal of Plasma Physics, 2023, 89 (4), pp.905890413. 10.1017/S0022377823000697. hal-04172983

HAL Id: hal-04172983

https://hal.science/hal-04172983

Submitted on 28 Jul 2023

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Efficient kinetic Lattice Boltzmann simulation
of three-dimensional Hall-MHD Turbulence

Raffaello Foldes¹,², Emmanuel Lévêque¹, , Raffaele Marino¹, , Ermanno Pietropaolo², Alessandro De Rosis³, Daniele Telloni⁴, and Fabio Feraco¹,⁵

¹Univ Lyon, CNRS, École Centrale de Lyon, INSA Lyon, Univ Claude Bernard Lyon I, LMFA
UMR 5509, F-69134 Ecully cedex, France

²Dipartimento di Scienze Fisiche e Chimiche, Università dell’Aquila, 67100 Coppito (AQ), Italy

³Department of Mechanical, Aerospace and Civil Engineering, The University of Manchester,
Manchester M13 9PL, United Kingdom

⁴National Institute for Astrophysics - Astrophysical Observatory of Torino, Via Osservatorio
20, I-10025 Pino Torinese, Italy

⁵Leibniz-Institute of Atmospheric Physics at the University of Rostock, Schlossstrasse 6,
Kühlungsborn, 18225, Germany

(Received xx; revised xx; accepted xx)

Simulating plasmas in the Hall-MagnetoHydroDynamics (Hall-MHD) regime 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 $\nabla \times ((\nabla \times \mathbf{B}) \times \mathbf{B})$ in the Faraday's induction law. The latter feeds back on the magnetic field 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). 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. 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. Spectral properties of the simulated plasma are finally compared with those obtained from numerical solutions from the well-established pseudo-spectral code GHOST. Furthermore, the LB simulations we present, varying the Hall parameter, highlight the transition from the MHD to the Hall-MHD regime, in excellent agreement with the magnetic field spectra measured in the solar wind.

# 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

† Email address for correspondence: raffaello.foldes@ec-lyon.fr
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 $J \times B$, where $J$ and $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 span 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 been 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 to resolve 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 account for 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 in terms of 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 (as mentioned before) when it comes to the integration of simulations in three dimensions and for many turnover times (Huba 2003). The 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). The first
three-dimensional MHD simulations based on the scheme proposed by Dellar have been
performed by Breyiannis & Valougeorgis (2004, 2006). 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. This
deficiency is not exclusive to MHD simulations but rather an inherent aspect of the BGK
collision operator itself. By utilizing a so-called Multi-Relaxation-Time (MRT) operator
defined in the space of moments, it becomes possible to explicitly dampen the non-
hydrodynamic modes and improve the stability (Higuera et al. 1989; Benzi et al. 1992;
d'Humieres 1994). Therefore, Pattison et al. (2008) and Riley et al. (2008) opted to use
MRT collision operators for the hydrodynamic parts of their lattice Boltzmann MHD
algorithms, whereas Dellar (2009) enhanced stability by considering MRT operators for
both the hydrodynamic and magnetic aspects. An entropic stabilization has also been
proposed by Flint & Vahala (2018), though leading to a more complicated scheme. These
advances encouraged us to pursue the LB modeling to simulate Hall-MHD turbulence,
an effort that has never been undertaken previously. In the present study, a MRT scheme
based on central-moments is considered for the evolution of the velocity field, while
dynamics of the magnetic field evolve under the action of a BGK collision operator,
following the scheme described in De Rosis et al. (2018). It is worth noting that Mendoza
& Muñoz (2008) had previously introduced a lattice Boltzmann algorithm for simulating
two charged species along with Maxwell's equations in the Hall MHD regime, as detailed
in the next section. Our approach, however, is more straightforward, neglecting the
electron inertia term. 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 ESA mission Solar Orbiter (Müller et al. 2020).

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 considered (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, we outline potential applications for the investigation
of space plasmas in §6, and draw conclusions in §7.

## 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:

$$ \nabla \cdot \mathbf{U} = 0 \qquad (2.1) $$

$$ \partial_t \mathbf{U} + (\mathbf{U} \cdot \nabla)\mathbf{U} = \frac{1}{\rho_0} \mathbf{J} \times \mathbf{B} - \frac{1}{\rho_0} \nabla \mathbf{p} + \nu \nabla^2 \mathbf{U} \qquad (2.2) $$

$$ \partial_t \mathbf{B} = \nabla \times (\mathbf{U} \times \mathbf{B} - \eta \nabla \times \mathbf{B}) \qquad (2.3) $$

$$ \nabla \cdot \mathbf{B} = 0 \qquad (2.4) $$

where $t$ is the time, $\rho_0$ is the mass density of the fluid, $\nu$ is the kinematic viscosity and $\eta$ is the magnetic resistivity. The electric current density is expressed as $\mathbf{J} = 1/\mu_0 \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 and electrons separately. For a fully ionized plasma in which the masses of ions (mainly 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 \nabla)\mathbf{U}] = en(\mathbf{E} + \mathbf{U} \times \mathbf{B}) - \nabla p_i + \nabla \cdot \sigma + \mathbf{R} \qquad (2.5) $$

$$ 0 = -en(\mathbf{E} + \mathbf{U}_e \times \mathbf{B}) - \nabla p_e - \mathbf{R} \qquad (2.6) $$

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

$$ \partial_t \mathbf{U} + (\mathbf{U} \cdot \nabla)\mathbf{U} = \frac{1}{\rho}\mathbf{J} \times \mathbf{B} - \frac{1}{\rho}\nabla p + \nu \nabla^2 \mathbf{U}. \qquad (2.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 (2.6), the Ohm's law becomes

$$ \mathbf{E} = -(\mathbf{U} - \frac{1}{en}\mathbf{J}) \times \mathbf{B} + \frac{1}{en}\nabla p_e + \frac{mf_{ie}}{e^2n}\mathbf{J}. \qquad (2.8) $$

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

$$ \partial_t \mathbf{B} = \nabla \times [(\mathbf{U} - \alpha_H \mathbf{J}) \times \mathbf{B}] + \eta \nabla^2 \mathbf{B} \qquad (2.9) $$

where $\alpha_H = 1/en$ is usually referred to as the Hall parameter and the magnetic resistivity $\eta = m f_{ie}/(e^2 n \mu_0)$. Let us note that, in general, $\nabla \times ((1/en)\nabla p_e) = -(1/en^2)\nabla n \times \nabla p_e$ (Kulsrud 2005). However, in the current context, we make the assumption that the electrons are isothermal, resulting in a dynamic pressure $p_e = nT_e$, where $T_e$ is a constant plasma temperature. Therefore, $\nabla n \times \nabla p_e = 0$. The Hall-MHD equations mentioned earlier include a finite ion-electron collision frequency, responsible for the $\mathbf{R}$-coupling term in (2.5) and (2.6). Additionally, they assume that the ion-ion collision frequency is large enough (much greater than the ion gyrofrequency) to permit the adoption of a standard Newtonian viscous stress in (2.7). The more comprehensive *Braginskii MHD* model, on the other hand, allows for the ion-ion collision frequency to be comparable to the ion gyrofrequency. Consequently, the Hall term emerges as just one component of the anisotropic relationship between electric current and electric field, and between stress and strain rate, with a preferred direction determined by the magnetic field. Dellar (2011)
provided a first LB approach to simulating the Braginskii MHD equations by modifying the hydrodynamics collision operator to depend on the magnetic field. Here, the main target of our simulations is represented by space plasmas providing a clear context for the use of Hall-MHD equations.

## 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 fluid 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, $Ma \equiv U_0/c_s$. Consequently, the incompressible equations should be replaced with the compressible formulation

$$ \partial_t \rho + \nabla \cdot (\rho \mathbf{U}) = 0 \qquad (2.10) $$

$$ \partial_t(\rho \mathbf{U}) + \nabla \cdot (\rho \mathbf{U} \otimes \mathbf{U} + \mathbf{p} \mathbb{I} + \frac{1}{2} |\mathbf{B}|^2 \mathbb{I} - \mathbf{B} \otimes \mathbf{B}) = \nabla \cdot \sigma \qquad (2.11) $$

in which $\sigma$ represents the viscous stress, the Lorentz force has been rewritten in a conservative form as the divergence of the Maxwell stress tensor $M_{\alpha\beta} = \frac{1}{2}|\mathbf{B}|^2\delta_{\alpha\beta} - B_{\alpha}B_{\beta}^{\dagger}$ 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, which allows for simplifying the Lorentz force as $(\nabla \times \mathbf{B}) \times \mathbf{B}$. The general form of the viscous stress is

$$ \sigma = \mu(\nabla \mathbf{U} + (\nabla \mathbf{U})^T) + (\zeta - \frac{2}{3}\mu)(\nabla \cdot \mathbf{U})\mathbb{I} \qquad (2.12) $$

where $\mu$ is the dynamic viscosity ($\nu = \mu/\rho$) and $\zeta$ is the bulk viscosity. 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

$$ p = \rho c_s^2 \qquad (2.13) $$

which is consistent with $O(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} + \nabla \cdot ((\mathbf{U} - \alpha_H \mathbf{J}) \otimes \mathbf{B} - \mathbf{B} \otimes (\mathbf{U} - \alpha_H \mathbf{J})) = \eta \nabla^2 \mathbf{B}. \qquad (2.14) $$

Let us remark 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} \nabla \times \mathbf{B}$. In the next sections, the developed LB scheme will conform to the set of equations (2.10), (2.11), (2.13) and (2.14). The divergence-free condition on $\mathbf{B}$ is preserved by (2.14), justifying that it is sufficient to impose $\nabla \cdot \mathbf{B} = 0$ initially. In the numerical modelling, 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.

<footnote>† The notation $(\mathbf{a} \otimes \mathbf{b})_{\alpha\beta} \equiv a_{\alpha}b_{\beta}$</footnote>
2018) and a high-connectivity D3Q27 lattice are used to integrate the dynamics of 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 satisfy sufficient isotropy, weakly-compressible Navier-Stokes dynamics can be reproduced with an $O(Ma^3)$ error. The macroscopic variables such as the fluid density $\rho$, momentum $\rho U$, or stress tensor $\Sigma$ are obtained as statistical moments of the particle distributions, i.e.

$$ \rho = \sum_{i=0}^{N-1} f_i \qquad (3.1) $$

$$ \rho U = \sum_{i=0}^{N-1} f_i c_i \qquad (3.2) $$

$$ \Sigma = \sum_{i=0}^{N-1} f_i c_i \otimes c_i \qquad (3.3) $$

by summing over the local mass densities $f_0, \dots, f_{N-1}$ of particles moving with velocities $c_0, \dots, c_{N-1}$, respectively. The sums replace here the integrals over $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 + (c_i \cdot \nabla) f_i = -\frac{1}{\tau} \left( f_i - f_i^{(0)}(\rho, U) \right) \qquad (3.4) $$

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, U)$, with a unique relaxation time $\tau = \nu/c_s^2$.

The Lattice keyword refers to the discretization in space and time of (3.4) with a set of microscopic velocity $c_0, \dots, 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).

The LB scheme can be expressed simply by using the change of variables

$$ \bar{f}_i = f_i + \Delta t / 2\tau (f_i - f_i^{(0)}) $$

originally introduced in (He et al. 1998), as

$$ \bar{f}_i(\mathbf{x} + c_i \Delta t, t + \Delta t) = \bar{f}_i(\mathbf{x}, t) - \omega \left( \bar{f}_i(\mathbf{x}, t) - f_i^{(0)}(\rho, U)(\mathbf{x}, t) \right) \qquad (3.5) $$
<figure><img src="image_2.png" /><figcaption>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 $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).</figcaption></figure>

where the discrete distribution functions $\bar{f}_i(\mathbf{x}, t)$ depend on the three spatial coordinates $\mathbf{x}$ and on time $t$. This change of variable comes from the trapezoidal rule used to approximate the integral of the collision operator (right-hand side of (3.4)) 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$ (Hénon 1987) so that

$$ \frac{1}{\omega} = \left( \frac{\nu}{c_s^2 \Delta t} + \frac{1}{2} \right) \qquad (3.6) $$

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

$$ \rho = \sum_{i=0}^{N-1} \bar{f}_i \qquad (3.7) $$

$$ \text{and } \rho \mathbf{U} = \sum_{i=0}^{N-1} \bar{f}_i \mathbf{c}_i. \qquad (3.8) $$

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

$$ \bar{f}_i(\mathbf{x} + \mathbf{c}_i \Delta t, t + \Delta t) = \bar{f}_i^*(\mathbf{x}, t) \qquad (3.9) $$

$$ \bar{f}_i^*(\mathbf{x}, t) = \bar{f}_i(\mathbf{x}, t) - \omega \left( \bar{f}_i(\mathbf{x}, t) - f_i^{(0)}(\rho, \mathbf{U})(\mathbf{x}, t) \right). \qquad (3.10) $$

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{\mathbf{c}_i \cdot \mathbf{U}}{c_s^2} + \frac{(\mathbf{c}_i \cdot \mathbf{U})^2}{2c_s^4} - \frac{\mathbf{U} \cdot \mathbf{U}}{2c_s^2} + \dots \right) \qquad (3.11) $$

with the weights $w_{\text{center}} = 8/27$, $w_{\text{face}} = 2/27$, $w_{\text{edge}} = 1/54$ and $w_{\text{corner}} = 1/216$ for the
D3Q27 lattice (He & Luo 1997). An expansion truncated at the second order is enough to recover the Navier-Stokes equations with an $O(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^{\text{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})(\mathbf{c}_i \cdot \mathbf{c}_i) - (\mathbf{c}_i \cdot \mathbf{B})^2 \right) \quad (3.12)$$

so that the second-order moment becomes,

$$\Sigma^{\text{mhd}(0)} = \sum_{i=0}^{N-1} f_i^{\text{mhd}(0)} \mathbf{c}_i \otimes \mathbf{c}_i = \rho \mathbf{U} \otimes \mathbf{U} + \text{pI} + \frac{1}{2} |\mathbf{B}|^2 \text{I} - \mathbf{B} \otimes \mathbf{B}. \quad (3.13)$$

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 (3.4) 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, due to the strongly nonlinear character of turbulent dynamics, spurious dissipative effects can occur as a result of the numerical integration of fluid-like equations over a very large number of grid points and of time-steps, as is the case for Hall-MHD turbulence.

A significant reduction of dissipation artifacts developing in turbulence simulations 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. 2006, 2007, 2015; De Rosis et al. 2018). This is the very framework adopted in lay-outing our code (details are given in the appendix A). A key ingredient of CM-LB schemes is the shift of the particle velocities by the local fluid velocity that defines a new set of local microscopic velocities $\tilde{\mathbf{c}}_i = \mathbf{c}_i - \mathbf{U}$ used for the CMs evaluation. Here, we consider the set of CMs as formally defined by

$$|\mathbf{k}\rangle \equiv [\mathbf{k}_0 \cdots \mathbf{k}_{26}]^\top = \mathbf{T}^\top |\bar{\mathbf{f}}\rangle \quad (3.14)$$

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

$$
\begin{align*}
|T_0\rangle &= |1\rangle \\
|T_1\rangle; |T_2\rangle; |T_3\rangle &= [\bar{c}_{ix}]^\top; [\bar{c}_{iy}]^\top; [\bar{c}_{iz}]^\top \\
|T_4\rangle; |T_5\rangle; |T_6\rangle &= [\bar{c}_{ix}\bar{c}_{iy}]^\top; [\bar{c}_{ix}\bar{c}_{iz}]^\top; [\bar{c}_{iy}\bar{c}_{iz}]^\top \\
|T_7\rangle; |T_8\rangle; |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 \\
|T_{10}\rangle; |T_{11}\rangle; |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}\bar{c}_{iz}^2]^\top \\
|T_{13}\rangle; |T_{14}\rangle; |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}\bar{c}_{iz}^2]^\top \\
|T_{16}\rangle &= [\bar{c}_{ix}\bar{c}_{iy}\bar{c}_{iz}]^\top \\
|T_{17}\rangle; |T_{18}\rangle; |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 \\
|T_{20}\rangle; |T_{21}\rangle; |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 \\
|T_{23}\rangle; |T_{24}\rangle; |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}^2]^\top \\
|T_{26}\rangle &= [\bar{c}_{ix}^2\bar{c}_{iy}^2\bar{c}_{iz}^2]^\top.
\end{align*}
$$

This set of vectors forms a simple relevant basis (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 (3.10) now generalizes as

$$
|k^*\rangle = |k\rangle - S(|k\rangle - |k^{(0)}\rangle) \quad \text{with} \quad |k^{(0)}\rangle = T^\top |f^{\text{mhd}(0)}\rangle \quad (3.15)
$$

where *S* is a diagonal matrix applied to each moment individually. Let us point out that the BGK collision is recovered by taking *S* = *ω*I. A proper choice for *S* is given by

$$
S = \operatorname{diag}[1, 1, 1, 1, \omega, \omega, \omega, \omega, \omega, \dots, 1] \quad (3.16)
$$

which ensures that mass and momentum are conserved by the collision operator and that kinematic viscosity is suitably taken into account. The bulk viscosity can be set separately from the shear viscosity and, here, it is implicitly defined by taking the trace of the second-order post-collision central-moment at equilibrium. Eventually, the post-collision distributions are obtained by returning to the space of the distributions through

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

before moving on to the streaming step (3.9).

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

We now present the LB scheme for the magnetic field introduced by Dellar (2002), here extended to encompass the Hall effect in simulating MHD turbulent plasmas. Following the works previously done by Croisille *et al.* (1995) and Bouchut (1999), Dellar (2002) proposed a decomposition of the magnetic field as:

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

where the sum spans a set of vector-valued densities $g_0, \dots, g_{M-1}$ associated with the
microscopic velocities $\xi_0, \dots, \xi_{M-1}$.

The magnetic field **B** is here provided by the zeroth-order moment of $\bar{\mathbf{g}}_i$ hinting that a
lattice with low connectivity should suffice to simulate its dynamics. In practice, a D3Q7
lattice with only seven velocities (see green arrows in Fig. 1) shall prove to be satisfactory in reproducing the magnetic field of Hall-MHD turbulent plasmas. Analogously to the fluid case, a LB scheme can be derived in order to simulate the induction equation in the form

$$ \bar{\mathbf{g}}_i(\mathbf{x} + \boldsymbol{\xi}_i \Delta t, t + \Delta t) = \bar{\mathbf{g}}_i(\mathbf{x}, t) - \omega_B \left( \bar{\mathbf{g}}_i(\mathbf{x}, t) - \mathbf{g}_i^{(0)}(\mathbf{U}, \mathbf{B})(\mathbf{x}, t) \right) \quad (3.19) $$

where the relaxation parameter $\omega_m$ is now related to the magnetic resistivity $\eta$ by

$$ \frac{1}{\omega_B} = \left( \frac{\eta}{C^2 \Delta t} + \frac{1}{2} \right) \quad (3.20) $$

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. \quad (3.21) $$

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) \quad (3.22) $$

with $W_{\text{center}} = 1/4$ and $W_{\text{face}} = 1/8$ for a D3Q7 lattice. By doing so, the first-order moment

$$ \sum_{i=0}^{M-1} \xi_{i\alpha} g_{i\beta}^{(0)} = \mathbf{U} \otimes \mathbf{B} - \mathbf{B} \otimes \mathbf{U} \quad (3.23) $$

would suitably reconstruct the transport term of the induction equation. Including the Hall correction in this equation is thus equivalent to upgrading the equilibrium densities, so that

$$ \Lambda_{\alpha\beta}^{(0)} = \sum_{i=0}^{M-1} \xi_{i\alpha} g_{i\beta}^{\text{Hall}(0)} = (\mathbf{U} - \alpha_H \mathbf{J}) \otimes \mathbf{B} - \mathbf{B} \otimes (\mathbf{U} - \alpha_H \mathbf{J}) \quad (3.24) $$

which is obviously possible by now considering

$$ g_{i\alpha}^{\text{Hall}(0)}(\mathbf{U}, \mathbf{B}, \mathbf{J}) = W_i \left( B_{\alpha} + \frac{1}{C^2} \xi_{i\beta} \left( (U_{\beta} - \alpha_H J_{\beta}) B_{\alpha} - B_{\beta} (U_{\alpha} - \alpha_H J_{\alpha}) \right) \right). \quad (3.25) $$

Nevertheless, $\mathbf{J}$ needs to be computed, possibly from the densities, in this expression. Note that although equilibrium distributions have been expanded up to the sixth order in $\mathbf{U}$, they have only been extended up to the second order in $\mathbf{B}$, which may sound contradictory. However, it is important to recognize that there is no continuous distribution for the magnetic field that is analogous to the Maxwell-Boltzmann distribution for the velocity. Thus, only the first two orders of the expansion can be reconstructed by matching the moments to the terms of the induction equation. Attempting to consider higher-order expansions would open a large variety of possibilities for defining non-physical moments, which is beyond the scope of the present work. An essential benefit of the LB framework is that the spatial derivatives of the magnetic field, thus $\mathbf{J}$, are self-consistently obtained (within an $O(Ma^3)$ error) from the first-order moment of the densities as

$$ J_{\gamma} = \varepsilon_{\alpha\beta\gamma} \frac{\partial B_{\alpha}}{\partial x_{\beta}} = -\varepsilon_{\alpha\beta\gamma} \frac{\omega_B}{C^2} \left( \Lambda_{\alpha\beta} - \Lambda_{\alpha\beta}^{(0)} \right) \quad (3.26) $$

where $\varepsilon_{\alpha\beta\gamma}$ is the Levi-Civita tensor and $\Lambda_{\alpha\beta} = \sum_{i=0}^{M-1} \xi_{i\alpha} \bar{g}_{i\beta}$ (Dellar 2002).
By replacing (3.24) in (3.26) we obtain a linear system readily solvable to obtain the current density **J**, namely

$$ \left( \mathbb{I} + \frac{2\alpha_H \omega_B}{C^2} \mathbf{M} \right) \mathbf{J} = \mathbf{J}_0 \qquad (3.27) $$

where

$$ \mathbf{M} = \begin{bmatrix} 0 & B_z & -B_y \\ -B_z & 0 & B_x \\ B_y & -B_x & 0 \end{bmatrix} \quad \text{and} \quad \mathbf{J}_0 = \begin{bmatrix} \mathbf{A}_{yz} - \mathbf{A}_{zy} - 2(U_y B_z - U_z B_y) \\ \mathbf{A}_{zx} - \mathbf{A}_{xz} - 2(U_z B_x - U_x B_z) \\ \mathbf{A}_{xy} - \mathbf{A}_{yx} - 2(U_x B_y - U_y B_x) \end{bmatrix}. \qquad (3.28) $$

Obviously, a solution exists only if it is possible to invert $\tilde{\mathbf{M}} = \mathbb{I} + (2\alpha_H \omega_B / C^2) \mathbf{M}$. It can be easily verified that $\det(\tilde{\mathbf{M}}) \neq 0$, which proves this solution exists and is unique. The current density obtained by solving (3.27) can then be used to compute the equilibrium densities (3.25) and proceed to the collision operation. It is fair to mention that in a similar vein, Dellar (2013) introduced a modification of the collision operator to incorporate MHD current-dependent resistivity, with the current being derived from the non-equilibrium components of the magnetic distribution functions. The expression of (3.26) also provides a consistent approximation of the divergence of the magnetic field. Indeed by taking the trace of the magnetic tensor, one obtains

$$ \nabla \cdot \mathbf{B} \approx -\frac{\omega_B}{C^2} \mathrm{Tr}(\mathbf{A}_{\alpha\beta}) \qquad (3.29) $$

by noticing that $\mathrm{Tr}(\mathbf{A}_{\alpha\beta}^{(0)}) = 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 $\nabla \cdot \mathbf{B} = 0$ corresponds with high accuracy to the condition $\mathrm{Tr}(\mathbf{A}_{\alpha\beta}) = 0$ in the LB framework (Dellar 2002). In practice, we have checked in our LB simulations that this condition was maintained throughout the runs, to machine round-off error.

### 3.3. Dimensionless formulation

In the following, the Hall-MHD equations are re-arranged in a dimensionless form in terms of the control parameter $\epsilon_H$, associated with the Hall parameter $\alpha_H = 1/\mathrm{ne}$. This control parameter is then recast in lattice units for practical LB purposes. Physical quantities in lattice units are hereafter indicated with the superscript $^{\mathrm{lbm}}$. In lattice units, the lattice spacing $\Delta x$ and the time-step $\Delta t$ 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$. Leveraging these characteristic quantities, (2.9) can be written in a dimensionless form as

$$ \left(\frac{U_0 B_0}{L_0}\right) \partial_t \mathbf{b} = \frac{1}{L_0} \nabla \times \left[ \left( U_0 \mathbf{u} - \frac{\alpha_H B_0}{\sqrt{\mu_0 L_0}} \nabla \times \mathbf{b} \right) \times (B_0 \mathbf{b}) \right] + \frac{\eta B_0}{L_0^2} \nabla^2 \mathbf{b} \qquad (3.30) $$

Dimensionless fields are here indicated by lowercase letters. This equation can be reduced to

$$ \partial_t \mathbf{b} = \nabla \times \left[ \left( \mathbf{u} - \epsilon_H \nabla \times \mathbf{b} \right) \times \mathbf{b} \right] + \frac{1}{\mathrm{Re}_m} \nabla^2 \mathbf{b} \qquad (3.31) $$

by defining the magnetic Reynolds number $\mathrm{Re}_m = U_0 L_0 / \eta$ and the dimensionless Hall parameter

$$ \epsilon_H = \frac{\alpha_H B_0}{\sqrt{\mu_0 L_0 U_0}}. \qquad (3.32) $$
We can treat in the same fashion the fluid momentum equation, where the reference scales are the same as those used to adimensionalize the induction equation. Therefore,

$$ \rho \frac{U_0^2}{L_0} [\partial_t \mathbf{u} + (\mathbf{u} \cdot \nabla) \mathbf{u}] = -\frac{1}{L_0} \nabla \rho c_s^2 + \rho \nu \frac{U_0}{L_0^2} \nabla^2 \mathbf{u} + \frac{B_0^2}{L_0} (\nabla \times \mathbf{b}) \times \mathbf{b} \quad (3.33) $$

which gives

$$ \partial_t \mathbf{u} + (\mathbf{u} \cdot \nabla) \mathbf{u} = -\frac{1}{\mathrm{Ma}^2} \frac{1}{\rho} \nabla \rho + \frac{1}{\mathrm{Re}} \nabla^2 \mathbf{u} + \left(\frac{V_A}{U_0}\right)^2 (\nabla \times \mathbf{b}) \times \mathbf{b} \quad (3.34) $$

where the control parameters are the Mach number Ma = $U_0/c_s$, the (fluid) Reynolds number Re = $U_0L_0/\nu$ and $V_A/U_0$, the Alfvén velocity being $V_A = B_0/\sqrt{\rho}$. The Hall number $\epsilon_H$ is given in lattice units by

$$ \epsilon_H^{\text{lbm}} = \frac{\alpha_H}{\sqrt{\mu_0}} \frac{[B_0/B_0^{\text{lbm}}]}{[L_0/L_0^{\text{lbm}}][U_0/U_0^{\text{lbm}}]} = \epsilon_H \frac{U_0^{\text{lbm}} L_0^{\text{lbm}}}{B_0^{\text{lbm}}}. \quad (3.35) $$

If one considers that the reference velocity $U_0$ corresponds to the Alfvén velocity ($U_0 = V_A$) and $\rho \approx 1$ for simplicity, one obtains that $U_0^{\text{lbm}} = B_0^{\text{lbm}}$ and

$$ \epsilon_H^{\text{lbm}} = \epsilon_H L_0^{\text{lbm}} = \epsilon_H N \quad (3.36) $$

with $N = L_0/\Delta x$ being the number of lattice points per reference length $L_0$. The Hall parameter $\epsilon_H$ can also be obtained as the ratio of two reference scales as

$$ \epsilon_H = \frac{V_A}{U_0 L_0} \sqrt{\frac{m}{\mu_0 n e^2}} = \frac{L_H}{L_0} \quad (3.37) $$

with

$$ L_H = \frac{V_A}{U_0} \sqrt{\frac{m}{\mu_0 n e^2}}. \quad (3.38) $$

In lattice units,

$$ \epsilon_H^{\text{lbm}} = \frac{L_H}{L_0} \frac{U_0^{\text{lbm}} L_0^{\text{lbm}}}{B_0^{\text{lbm}}} = \frac{L_H}{\Delta x} \frac{U_0^{\text{lbm}}}{B_0^{\text{lbm}}} = \frac{U_0^{\text{lbm}} L_H^{\text{lbm}}}{B_0^{\text{lbm}}}. \quad (3.39) $$

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^{\text{lbm}} = L_H^{\text{lbm}}$ and this corresponds to the number of lattice points per ion inertial length. It is assumed that the dynamics of a MHD plasma develops under the influence of the Hall effect at scales $l$ smaller than $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 \le \Delta x/c_{\max} \quad (3.40) $$

where $c_{\max}$ refers to the largest speed at which a signal propagates in the solution. In the context of Hall-MHD, $c_{\max}$ should be identified with the largest phase speed of the whistler waves. When the plasma dynamics in the direction of the magnetic field $B$ is dominant, the phase speed of the whistler waves varies as $c_w(k) = kV_A^2/\omega_{ci}$ with the wavenumber $k$; $V_A$ is the Alfvén velocity, while $\omega_{ci}$ is the ion cyclotron frequency. In physical units, $\omega_{ci} = eB/m_i$ and $V_A = B/\sqrt{\mu_0 n m_i}$, $m_i$ being the mass of ions and the Hall parameter being $\alpha_H = 1/\text{ne}$. Therefore, one obtains that the time-step decreases
<figure><img src="image_3.png" /><figcaption>FIGURE 2. Three-dimensional rendering of the initial condition as indicated in (4.5) and (4.6). The magnitudes of the fluid velocity (left), magnetic field (center) and density of electric current (right) are here displayed for $N = 128$.</figcaption></figure>

quadratically with the grid spacing, as

$$ \Delta t \le \frac{\mu_0 (\Delta x)^2}{\pi \alpha_H B} \qquad (3.41) $$

assuming the largest attainable wavenumber to be $k_{max} = \pi/\Delta x$ in the context of the Hall-MHD turbulence. This condition can be rewritten accounting for the rescaling of the magnetic field by $\mu_0^{1/2}$

$$ \Delta t \le \frac{(\Delta x)^2}{\pi (\alpha_H / \sqrt{\mu_0}) B_0} \qquad (3.42) $$

which finally yields in lattice units to

$$ 1 \le \frac{1}{\pi \epsilon_H^{\text{lbm}} L_0^{\text{lbm}} U_0^{\text{lbm}}} = \frac{B_0^{\text{lbm}}}{\pi \epsilon_H (N U_0^{\text{lbm}})^2} \qquad (3.43) $$

where (3.32) is used to retrieve $\alpha_H / \sqrt{\mu_0}$, and $N = L_0^{\text{lbm}}$. If $U_0 = V_A$, the CFL condition for whistler waves can be reformulated as a condition on the Mach number $Ma = \sqrt{3}U_0^{\text{lbm}}$, which is in turn written as

$$ Ma \le \left( \frac{\sqrt{3}}{\pi} \right) \frac{1}{\epsilon_H N^2}. \qquad (3.44) $$

This condition recalls 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.

## 4. Results

Our LB scheme and code FLAME is now validated against the analytical solution of the incompressible and dissipative Hall-MHD equations proposed in (Xia & Yang 2015). The latter is used as a benchmark to evaluate accuracy and convergence of the numerical solutions for different values of the control parameters (in a regime of low Reynolds numbers in which the aforementioned analytical solution holds). A further validation was done focusing on the MHD range of scales, this time in a regime of high Reynolds numbers. The solutions of the MHD dynamics produced by FLAME were
<table><caption>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_0L_0$, which means that the Reynolds number $Re = 1/ν$.</caption><thead><tr><th>Parameters</th><th>Values</th></tr></thead><tbody><tr><td>Resolution: <i>N</i></td><td>32, 64, 96, 128</td></tr><tr><td>Mach number: Ma [×10<sup>-2</sup>]</td><td>1.0, 0.7, 0.5, 0.3</td></tr><tr><td>Kinematic viscosity: ν [×10<sup>-3</sup>]</td><td>1.0, 0.5, 0.33, 0.25, 0.2, 0.17, 0.14, 0.13, 0.11</td></tr></tbody></table>

compared in this case with those obtained with a well-established pseudo-spectral solver, widely used for turbulent plasma simulations, namely the Geophysical High-Order Suite for Turbulence (GHOST, Mininni et al. (2011); Rosenberg et al. (2020)). Finally, the physical consistency of the output and the computational performance were evaluated when accounting explicitly for the Hall effect in the turbulent regime. This allowed us to assess the reliability of our code in simulating the multi-scale dynamics generated by turbulent flows at high Reynolds numbers, and in reproducing the transition from the MHD to the Hall-MHD spectral range (at sub-ion scales).

The FLAME code relies on a multi-GPU (Graphic Processing Unit) implementation of the LB scheme in order to reach high resolution that optimizes the computational times. Massive multi-threading is handled within the OpenCL (Open Computing Language) framework, allowing a high portability of the code. 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 the latter. Turbulence simulations were run on a cluster equipped with NVIDIA A100-40Gb GPU cards, hosted at the CINECA supercomputing center (Italy).

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

Due to their high computational cost, the availability in the literature of plasma simulations reproducing the Hall-MHD range of scales (in three dimensions) is much less than for the MHD case. Moreover, Hall-MHD simulations are in general performed using pseudo-spectral codes (Ferrand et al. 2022; Meyrand & Galtier 2012; Gómez et al. 2010; Yadav et al. 2022), which integrate of course the dynamical equations in the Fourier space. Interestingly, Mahajan & Krishan (2005) derived an analytical solution for the non-dissipative Hall-MHD equations, then extended by Xia & Yang (2015) with the inclusion of dissipative effects. This solution is used in the following to test the stability and convergence of FLAME. Encompassing dissipative effects Xia & Yang (2015), this analytical solution allowed us to quantify as well the numerical dissipation spuriously introduced by our scheme.

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

$$ \mathbf{u}(\mathbf{x}, t) = \mathbf{u}'(\mathbf{x}, t) \quad \text{and} \quad \mathbf{b}(\mathbf{x}, t) = \hat{\mathbf{e}}_z + \mathbf{b}'(\mathbf{x}, t) \qquad (4.1) $$

where the fluctuating velocity and magnetic fields are damped circular polarized waves given respectively by
$$ \mathbf{u}'(\mathbf{x}, t) = [B(\hat{\mathbf{e}}_x + i\hat{\mathbf{e}}_y) \exp(ikz - i\omega_{\pm}t) \\ + C(\hat{\mathbf{e}}_y + i\hat{\mathbf{e}}_z) \exp(ikx) \\ + A(\hat{\mathbf{e}}_z + i\hat{\mathbf{e}}_x) \exp(iky)]e^{-\nu k^2 t} \quad (4.2) $$

and

$$ \mathbf{b}'(\mathbf{x}, t) = \alpha_{\pm} \mathbf{u}'(\mathbf{x}, t) \quad (4.3) $$

in complex notations. The amplitudes A, B, and C are arbitrary real values. The ambient magnetic field here is 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 the real part is a solution. The pulsation $\omega_{\pm} = -\alpha_{\pm}k$, where $\alpha_{\pm}$ depends itself on the wavenumber k as

$$ \alpha_{\pm} = -\frac{1}{2}\epsilon_H k \pm \sqrt{\frac{\epsilon_H^2 k^2}{4} + 1} \quad (4.4) $$

The magnetic Prandtl number is assumed equal to unity in obtaining this solution and the reference velocity is assumed equal to the Alfvén velocity, i.e. $U_0 = V_A$ in (3.34). Finally, it is worth mentioning that this analytical solution holds in a strictly incompressible framework, which, given the intrinsically compressible nature of the LB scheme, prescribes that our simulations must be run at a (very) low Mach number so that relative density fluctuations generated by the code remain negligible. In our investigations, the Hall-MHD equations have been integrated in a cubic box of size $L_0 = 2\pi$. The evolution of the velocity field is deterministic from the initial condition

$$ \begin{aligned} u_x^{\text{lbm}}(\mathbf{x}, 0) &= U_0^{\text{lbm}} \left( B \sin \left( \frac{4\pi z_k}{N} \right) + A \cos \left( \frac{4\pi y_j}{N} \right) \right) \\ u_y^{\text{lbm}}(\mathbf{x}, 0) &= U_0^{\text{lbm}} \left( B \cos \left( \frac{4\pi z_k}{N} \right) + C \sin \left( \frac{4\pi x_i}{N} \right) \right) \\ u_z^{\text{lbm}}(\mathbf{x}, 0) &= U_0^{\text{lbm}} \left( C \cos \left( \frac{4\pi x_i}{N} \right) + A \sin \left( \frac{4\pi y_j}{N} \right) \right) \end{aligned} \quad (4.5) $$

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^{\text{lbm}}$ is related by construction to the Mach number through $U_0^{\text{lbm}} = Ma/\sqrt{3}$. The magnetic field is initially proportional to the fluid velocity with

$$ \begin{aligned} b_x^{\text{lbm}}(\mathbf{x}, 0) &= \alpha_+ u_x^{\text{lbm}}(\mathbf{x}, 0) \\ b_y^{\text{lbm}}(\mathbf{x}, 0) &= \alpha_+ u_y^{\text{lbm}}(\mathbf{x}, 0) \\ b_z^{\text{lbm}}(\mathbf{x}, 0) &= \alpha_+ u_z^{\text{lbm}}(\mathbf{x}, 0) + U_0^{\text{lbm}} \end{aligned} \quad (4.6) $$

since $U_0 = V_A$. For sake of simplicity, the initial density is set to one everywhere in the space. The (normalized) Hall parameter is fixed at $\epsilon_H = 1$, that is $\epsilon_H^{\text{lbm}} = N$ according to (3.36). This value ensures that the solution is affected by the Hall effect with $L_H = L_0$ from (3.37). A three-dimensional rendering of the initial conditions expressed in (4.5) and (4.6) is displayed in Fig. 2. With this initialization, the current density $\mathbf{J} = \nabla \times \mathbf{B}$ is non-zero at $t = 0$. The parameters used in the different simulations are reported in Tab. 1. The Mach number is always small enough for the plasma to approach the incompressible limit and in order to reduce the intrinsic discretization 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
<figure><img src="image_4.png" /><figcaption>FIGURE 3. Temporal evolution of the velocity magnitude $|\mathbf{u}|(0, t)$. Comparison between the analytical (dashed line) and the numerical solutions (symbols) at different lattice resolutions. The Mach number $Ma = 0.003$ and the kinematic viscosity $\nu = 3.3 \cdot 10^{-4}$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIGURE 4. Temporal evolution of the velocity magnitude $|\mathbf{u}|(0, t)$ for different values of the resolution ($N$) and viscosity ($\nu$) at fixed Mach number $Ma = 0.007$.</figcaption></figure>

in the latter case, requiring a significant reduction of the Mach number (with $\epsilon_H = 1$). Results obtained for $\alpha_+$ and the velocity field only ($\mathbf{b} = \alpha_+\mathbf{u} + \hat{\mathbf{e}}_z$) are presented in the following.

## 4.2. Stability and incompressibility

The stability of the scheme was tested exploring the parameter space defined through the Mach number, the lattice resolution and the kinematic viscosity (see Tab. 1). The analytical solution introduced by Xia & Yang (2015) is such that the nonlinear terms in the incompressible dissipative Hall MHD equations are strictly zero. In practice, physical instabilities triggered by numerical errors do naturally develop and grow in time in simulations whenever the viscosity is too small, eventually inducing the transition to a turbulent state. Therefore, the numerical stability and accuracy of FLAME were assessed in runs in which the viscosity was sufficiently high to prevent such transition to turbulence. The typical temporal evolution of the velocity at a fixed location in the simulation domain is shown in Fig. 3. The solution appears as a damped wave
<figure><img src="image_6.png" /><figcaption>FIGURE 5. Relative density fluctuations at different values of the Mach number (Ma) and kinematic viscosity (ν) as a function of the resolution (N). In our simulations, the reference density ρ₀ is fixed at unity.</figcaption></figure>

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 Ma = 0.007 are shown in Fig. 4 for 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⁻⁷ – 10⁻⁶, indicating a very good convergence towards the incompressible limit in all the simulations presented. Furthermore, the results confirm that the amplitude of density fluctuations decreases with the Mach number.

### 4.3. Dispersion and dissipation errors

The dispersion and dissipation errors of the LB scheme implemented in FLAME are now assessed. In this analysis, the dispersion error (or phase error) is computed by evaluating the shift in time between the local maxima of the numerical solution and the analytical wave solution (see Fig. 3). Therefore, tagging as $t_i^{\max}$ and $\bar{t}_i^{\max}$ the positions in time of the maxima of the numerical and analytical solution (at a fixed location) respectively, the average value of the relative dispersion error can be defined as

$$ \varepsilon_{\phi} = 1 - \frac{1}{M} \sum_{i=0}^{M-1} \frac{t_{i+1}^{\max} - t_i^{\max}}{\bar{t}_{i+1}^{\max} - \bar{t}_i^{\max}} \quad (4.7) $$

over M oscillating periods. For practical purposes we have used M = 10. As expected, it can be observed in Fig. 6 how the dispersion error is very small and decreases as the resolution N of the grid increases, showing a power-scaling law close to 1/N². This confirms a second-order accuracy of the LB scheme. We also found that the dispersion error exhibits a rather constant behavior 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 synchronizing the phases of numerical and analytical solutions, the (relative) dissipation error is evaluated by comparing the velocity magnitude of the two solutions, i.e.

$$ \varepsilon_r = \left[ \frac{1}{M} \sum_{i=0}^{M-1} \left( \frac{u(t_i) - \bar{u}(t_i)}{\bar{u}(t_i)} \right)^2 \right]^{1/2} \quad (4.8) $$
<figure><img src="image_7.png" /><figcaption>FIGURE 6. Relative dispersion error as defined by (4.7) for different values of the kinematic viscosity $\nu$ and the Mach number $Ma$. The error decreases as $N^{-2}$ as expected for a LB scheme.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIGURE 7. Acoustic scaling (Ma constant) of the relative dissipation error for different values of the kinematic viscosity ($\nu$) at fixed Mach number (Ma). The error decreases as $N^{-2}$.</figcaption></figure>

The dissipation error provides a first measure of the numerical dissipation. Two different scaling behaviors are considered, namely the so-called acoustic and diffusive scaling (Krueger et al. 2016). The acoustic scaling consists in keeping the Mach number fixed while monitoring the convergence rate of the error $\epsilon_r$ for different Reynolds numbers, 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 behavior 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 scheme.

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^2t)$ of the numerical solution. By decomposing $\tilde{\nu}$ into the sum of a physical and a (spurious) numerical viscosity, $\tilde{\nu} = \nu + \nu_{\text{num}}$, the ratio between these two contributions reads as

$$ \epsilon_{\nu} = \frac{\nu_{\text{num}}}{\nu} = \frac{\tilde{\nu} - \nu}{\nu}. \qquad (4.9) $$

The results obtained for the viscosity error $\epsilon_\nu$ are shown in Fig. 9. Here, we found that the numerical viscosity represents only a small percentage of the estimated total viscosity, and it 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 is independent from the physical viscosity and the Mach number, whereas it only depends on the lattice resolution.

Finally, despite Dellar (2002) showed that a D3Q7 lattice was sufficient to reliably
<figure><img src="image_9.png" /><figcaption>FIGURE 8. Diffusive scaling ($\nu$ constant) of the (relative) dissipation error for different values of the Mach number (Ma) at fixed kinematic viscosity ($\nu$).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIGURE 9. Scaling of the ratio between the numerical and kinematic viscosities $\epsilon_\nu = \nu^{\text{num}}/\nu$ with the grid resolution ($N$) for different Mach numbers (Ma) and kinematic viscosity. The second-order accuracy of the LB algorithm is highlighted by the black lines, i.e. $\epsilon_\nu \sim O(\Delta x^2)$.</figcaption></figure>

account for the dynamics of each component of the magnetic field, in order to check the validity of this statement, LB simulations with enhanced connectivity have been performed here to investigate whether a more isotropic representation of the magnetic densities would significantly 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), thus confirming what was reported in Dellar (2002). A plausible explanation of this lies on the fact that the magnetic field is represented as a zeroth order moment of densities for each component (see (3.18)). Therefore, a 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 between the dynamics of MHD plasmas simulated with FLAME and the outputs obtained with the GHOST pseudo-spectral solver for high-resolution simulations, when both codes perform the same decaying test run initialized with the classical Orszag-Tang (OT) vortex problem (Orszag & Tang 1979). Indeed, the OT solution is often considered as a prototypical flow to study freely evolving MHD turbulence. The GHOST solver has been widely used to tackle a variety of problems related to both geophysical fluids and space plasmas (Marino et al. 2013; Pouquet & Marino 2013; Marino et al. 2014, 2015a; Mininni et al. 2002, 2003, 2006; Gómez et al. 2010; Pouquet et al. 2019). It is a well-established community code available on https:
<figure><img src="image_11.png" /><figcaption>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.</figcaption></figure>

//github.com/pmininni/GHOST. GHOST is a hybrid MPI/OpenMP/CUDA-parallelized framework that hosts a variety of solvers having also GPU capability, delivering high performance, robust results and an optimal scaling up to hundreds of thousand computing cores. It relies on a second-order Runge-Kutta scheme for time integration and is de-aliased based on the classical two-third rule. As a pseudo-spectral de-aliased code, it provides very high accuracy in resolving the spatial scales (Patterson & Orszag 1971). The OT vortex problem prescribes the following initialization for the velocity and magnetic fields:

$$ \mathbf{U}(\mathbf{x}, 0) = U_0 [-2 \sin y; 2 \sin x; 0] $$

$$ \mathbf{B}(\mathbf{x}, 0) = B_0 [-2 \sin 2y + \sin z; 2 \sin x + \sin z; \sin x + \sin y] $$

with $U_0 = 1$ and $B_0 = 0.8$ in a cubic box of size $2\pi$.

In the simulation performed here, the Reynolds number attains values up to $Re = UL/\nu \simeq 1600$ when the flow reaches its peak of dissipation. The small-scale energy dissipation is 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$. In the definition of the Reynolds number, $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 $Ma = 0.025$. The number of grid points in each direction is $N = 512$.

The time evolution of the mean magnetic dissipation, as well as the kinetic and magnetic energies, are shown in Fig. 10 for two realizations of LB and pseudo-spectral simulations of the same OT problem. The simple visual inspection of the runs shows that the agreement between FLAME and GHOST is very satisfactory for the cases under study. Only a slight underestimation of the magnetic dissipation in the FLAME run can be observed for a few time steps after the peak of the current density $\mathbf{J} = \nabla \times \mathbf{B}$. Let us recall at this stage that $\mathbf{J}$ is directly obtained from the magnetic densities in the LB simulation, and is not inferred by differentiating the magnetic field.

A more detailed comparison is provided by looking at the Fourier decomposition of the fields obtained with the two codes. The kinetic and magnetic energy spectra are displayed in Fig.11 at the peak of the magnetic dissipation. The kinetic energy spectrum of the LB simulation seems over-damped at high wavenumbers. This is related to a known drawback of the moment-based collision operator, which ensures 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
<figure><img src="image_12.png" /><figcaption>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.</figcaption></figure>

numerical dissipation is reduced and the spectrum of the FLAME run gets very close to that of the pseudo-spectral solution. This observation is qualitatively consistent with the statement made in Shen et al. (2018) that LB needs about twice the spatial resolution of a pseudo-spectral simulation to achieve similar accuracy in turbulent flows.

Concerning the magnetic energy spectrum, the results from both simulations perfectly match, reflecting the fact that the 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 $k_{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_{max} = N/2$ in LB simulations. Particular attention is now paid to the wavenumber-by-wavenumber energy budget of the MHD equations. Starting from (2.11) and (2.14), the (total) energy flux across wavenumber k can be defined as

$$S_{\text{MHD}}(k) = \sum_{|\mathbf{k}'|<k} \Re[\mathcal{F}(\mathbf{U})^* \cdot (\mathcal{F}(\mathbf{U} \cdot \nabla \mathbf{U}) - \mathcal{F}(\mathbf{J} \times \mathbf{B})) - \mathcal{F}(\mathbf{B})^* \cdot \mathcal{F}(\nabla \times (\mathbf{U} \times \mathbf{B}))] \quad (4.10)$$

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

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

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}'|<k} E(\mathbf{k}') = -S_{\text{MHD}}(k) - \mathcal{D}(k). \quad (4.12)$$

We would like to mention that the contribution of the pressure term (not shown here) is negligible in the context of these simulations. The fluxes obtained for the LB and pseudo-spectral OT implementations (with N = 512) are displayed in detail in Fig. 12. A satisfactory agreement is observed in particular for the non-linear energy transfer terms, over the entire range of resolved wavenumbers. The slight over-dissipative nature of the LB scheme is again evidenced in the output of the dissipation term $\mathcal{D}(k)$ at very high wavenumbers.
<figure><img src="image_13.png" /><figcaption>FIGURE 12. Various third-order energy fluxes across wavenumber *k* as reported in (4.11). Comparison is made between LB (left) and de-aliased pseudo-spectral (right) GHOST simulations. Let us notice that $\partial_t \sum_{|k'|<k} E(k') = -S_{\text{MHD}}(k) - \mathcal{D}(k) < 0$ as expected for freely-decaying MHD turbulence.</figcaption></figure>

<table><caption>TABLE 2. Parameters of Hall-MHD turbulence runs. Re, Ma and Pr<sub>m</sub> 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_{tot}/\tau_0$ and the time at which the peak of current density occurs is $t_{peak}/\tau_0$ in units of the reference time scale $L_0/U_0$. The Mach number satisfies the CFL condition (3.4) imposed by whistler waves.</caption><thead><tr><th>Run</th><th>N</th><th>Ma [×10<sup>-4</sup>]</th><th>Re</th><th>Pr<sub>m</sub></th><th>ε<sub>H</sub></th><th>t<sub>tot</sub>/τ<sub>0</sub></th><th>t<sub>peak</sub>/τ<sub>0</sub></th></tr></thead><tbody><tr><td>I</td><td>512</td><td>7.0</td><td>4400</td><td>1</td><td>0.0025</td><td>48.2</td><td>31.8</td></tr><tr><td>II</td><td>512</td><td>1.0</td><td>5240</td><td>1</td><td>0.01</td><td>35.8</td><td>33.6</td></tr><tr><td>III</td><td>512</td><td>0.625</td><td>7150</td><td>1</td><td>0.025</td><td>29.6</td><td>26.5</td></tr><tr><td>IV</td><td>768</td><td>0.6</td><td>6000</td><td>1</td><td>0.015</td><td>36.6</td><td>32.0</td></tr></tbody></table>

## 5. High-Resolution simulations of 3D Hall-MHD plasmas

FLAME was used to simulate plasma dynamics in a regime in which the Hall-MHD term is non-negligible. In particular, the governing equations have been integrated in a triply periodic cubic lattice of size $L = 2\pi$ with resolution 512³ and 768³, initialized with the OT vortex as described in the previous section, for different values of the Hall parameter (see Tab. 2). The Mach number was adjusted to the Hall parameter in order to accommodate the CFL condition based on the time-scale of the whistler waves (see 3.4). The Reynolds number is estimated here at the peak of magnetic dissipation (indicated by the vertical dashed lines in Fig 13), at which the plasma is assumed to have reached a fully developed turbulent state. For a 512³ lattice dimension, only three 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). This led to a wall-clock computational time of 10, 55 and 69 hours respectively for the three runs indicated in Tab. 2 to pass the peak of magnetic-energy dissipation. The computational times reported above are comprehensive of the time required to transfer the three-dimensional vector fields (u, B and 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 u and B is shown in Fig. 14 for the
<figure><img src="image_14.png" /><figcaption>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. 16 have been averaged.</figcaption></figure>

simulation at the highest resolution (run IV in Tab. 2), again taken at the peak of the magnetic-energy dissipation. The three-dimensional visualization is displayed together with the kinetic and magnetic energy spectra, the latter showing two regimes above and below the ion inertial length $d_i$. At the same time, the small-scale activity visible in Fig. 15 for the electric current density **J** and the vorticity $\omega = \nabla \times u$, emphasizes the presence of current sheets, Kelvin-Helmholtz instabilities and vortices, emerging as the disordered structures characteristic of the Hall effect (Miura & Araki 2014). Furthermore, we have found that increasing the intensity of the Hall effect produces a faster development of turbulence in the plasmas under study due to the presence of both whistler and Hall-drift waves, propagating quicker than the Alfvén waves in the ideal MHD (Huba 2003). This is consistent with the behavior captured in Fig. 13 for the three runs, increasing the Hall parameter. The kinetic and magnetic energy spectra averaged over a time interval (around the peak of dissipation, as indicated by the shaded areas in Fig. 13) are plotted in Fig. 16 for each run at resolution 512³. As expected, increasing the value of the Hall parameter $\epsilon_H$ (indicated by the vertical dash-dotted line in Fig. 16) produces a shift of the Hall length-scale $L_H$ towards larger scales, hence a shrink of the Kolmogorov's $k^{-5/3}$ power law range in both kinetic and magnetic energy spectra. A very surprising and promising feature of these simulations is the behavior of the magnetic spectra in the Hall-MHD regime. In fact, at wavenumbers $k > k_H$, the spectrum develops (as $\epsilon_H$ increases) a power-law scaling that is in perfect agreement with the $k^{-2.73}$ scaling obtained from the spectral analysis of solar wind measurements at sub-ion scales, as reported in (Kiyani et al. 2015).

In the MHD regime, the time step of the (compressible) Lattice Boltzmann runs is constrained by the need for resolving sound waves. Therefore, the time-step of an LB simulation is typically much smaller compared to the time-step of equivalent (incompressible) pseudo-spectral simulation, the ratio between the two time-steps being typically the Mach number (Horstmann et al. 2022). Therefore, in the case of the MHD, the advantage for our LB scheme in terms of turn-around times is not that big compared to standard pseudo-spectral simulations. The situation is different when it comes to the simulation of plasmas in the Hall-MHD regime, where the time-steps of the two methods are identically constrained by the speed of whistler waves. In this case, the efficiency of the
<figure><img src="image_15.png" /><figcaption>FIGURE 14. Results of a 3D Hall-MHD simulation at $N = 768^3$ taken at the peak of dissipation. The upper panels show a three-dimensional visualization of the normalized velocity field magnitude $v' = (v - \bar{v})/\sigma_v$ with $v = |v|$ (upper left) and the relative power density spectrum (upper right). The same in the two bottom panels for the normalized magnetic field $B'$. The slopes $k^{-5/3}$ and $k^{-2.73}$ are given as a reference.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>FIGURE 15. Three-dimensional rendering of the normalized current density $J' = (J - \bar{J})/\sigma_J$, with $J = |J|$ (left), and vorticity $\omega'$ (right) taken at the same time as in Fig. 14.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>FIGURE 16. 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$.</figcaption></figure>

LB scheme (exploiting the computational power of GPU accelerators) becomes a major
advantage with respect to pseudo-spectral simulations, leading to wall-clock turn-around
times that are significantly smaller for LB schemes, and for FLAME in particular. Finally,
we would like to mention that an extension of FLAME allowing the simulation of the
electron MHD dynamics would simply consist in modifying the equilibrium distributions
for the magnetic field in the LB scheme, by neglecting the bulk velocity U with respect
to the Hall current $\alpha_H \mathbf{J}$ in (3.25).

**6. Hall-MHD simulations for space plasma turbulence investigations**

Space plasmas, whose dynamics involve the turbulent transport of the energy from
very large (Adhikari et al. 2015a,b) down to the very small scales (Cerri et al. 2016)
due to their large Reynolds numbers Matthaeus et al. (2008); Parashar et al. (2019),
do actually develop well-defined 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 two of the most recent space
missions: Solar Orbiter (SO; Müller et al. 2020) and Parker Solar Probe (PSP; Fox et al.
2016). Fig. 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 these state of the art 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 to a low-
speed solar wind stream measured on July 14 – 15, 2020, whose origin was identified
in a coronal streamer and pseudo-streamer configuration (D’Amicis et al. 2021). In
Tab. 3 we report the characteristic parameters of these solar wind samples. It is worth
recalling 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, the latter being in general significantly larger than the former. For values of
density $n$, temperature $T$ and $\beta$ typical for space plasmas, the relation $\rho_i \lesssim d_i$ is valid.
However, it has been remarked by several authors how, for $\beta \sim O(1)$, these 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 breaking point identifying the
transition between the end of the MHD range and the beginning of the range where plasma kinetic effects become relevant, 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 SO and the 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 the MHD turbulent inertial range (where energy is adiabatically transferred to smaller and smaller scales), that is characterized by a Kolmogorov-like spectrum Marino et al. (2011, 2012); Marino & Sorriso-Valvo (2023), to a range where the kinetic effects begin to dominate and in which the energy gets dissipated (at the bottom of such range), ultimately heating the solar wind plasma Marino et al. (2008). As is known from spacecraft observations, fluid and kinetic scales in the solar wind are characterized by different power-law spectral exponents. Features of these 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 controlling the evolution of turbulence in the interplanetary space are still matter of investigation. Nonlinear interactions (Bruno & Carbone 2013), expansion-driven 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 and dissipation 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 in the solar wind, are major open questions in the space plasma community that could be effectively targeted by means of numerical investigations produced with FLAME, which allows capturing the transition between MHD and Hall-MHD regimes (Fig.16), like the more standard pseudo-spectral codes. Another puzzle of solar and space plasma dynamics that can be tackled with our LB code is how magnetic switchbacks observed in the solar corona as well as in the solar wind do contribute to the local heating of the plasma. The switchbacks are intermittent magnetic-field polarity reversals widely observed in the heliosphere (Bale et al. 2019) and in the solar corona (Telloni et al. 2022b), that are thought to play a major role in the acceleration and heating of the solar wind. However, characterizing their contribution to the plasma energetics among other plasma processes is a challenging task, for which it is important to run highly accurate 3D Hall-MHD numerical simulations, able to resolve an extended dynamical range with the largest possible scale separation. A first implementation of FLAME aiming at demonstrating the decaying nature of solar
<table><caption>TABLE 3. Main solar wind parameters computed at time where the solar wind samples used to produce the power spectra in Fig. 17 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 (SO and PSP) from the Sun $D_{sun}$.</caption><thead><tr><th>Probe</th><th>ρ̄[cm<sup>3</sup>]</th><th>V̄[km/s]</th><th>T̄[MK]</th><th>B̄[nT]</th><th>β</th><th>d<sub>i</sub>[km]</th><th>D<sub>sun</sub>[AU]</th></tr></thead><tbody><tr><td>PSP</td><td>419</td><td>622</td><td>1.89</td><td>332</td><td>0.50</td><td>57.5</td><td>0.09</td></tr><tr><td>SO</td><td>16</td><td>429</td><td>1.50</td><td>6.76</td><td>3.59</td><td>11.1</td><td>0.64</td></tr></tbody></table>

<figure><img src="image_18.png" /><figcaption>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</figcaption></figure>

wind turbulence has been presented in (Sorriso-Valvo, L. et al. 2023), where a direct comparison between the simulated fields and the observations performed by the Helios 2 spacecraft is proposed, showing very good agreement.

## 7. Conclusions

The LB approach extends the horizon for the numerical investigation of plasma dynamics. Stability issues, which have long been a handicap for the implementation of the LB method to investigate turbulent flows, are now mostly solved thanks to the use of improved collision operators that do not compromise the accuracy of the numerical solutions. Furthermore, the computational efficiency of the LB schemes on many-core devices such as GPUs allows for advantageous turn-around times. A major advantage of dealing with a kinetic representation at the level of the numerical scheme, is that the derivatives of the magnetic field are directly embedded in the solution, allowing for an intrinsically accurate description of the current density since it does not require further implementations of a differentiation scheme. The study presented here shows that the LB approach provides a valuable and efficient numerical tool to simulate Hall-MHD plasma turbulence. Furthermore, the LB approach a priori allows us to add complexity
to the plasma, such as thermal effects, multi-species, complex geometries, etc. at the cost of new coupled lattice dynamics and boundary conditions, therefore preserving the computational performance. Extended MHD codes currently utilized for tokamak applications such as NIMROD (Sovinec & King 2010), BOUT++ (Dudson et al. 2015) or JOREK (Hoelzl et al. 2021) rely on implicit or semi-implicit timestepping to ensure stability with longer time-steps than what is imposed by the CFL condition for explicit timestepping. In our (explicit) scheme, the time-step was established by default to meet this later condition according to (3.4). It would be valuable to investigate the extent to which this constraint on $\Delta t$ could be relaxed while maintaining stability, due to the magnetic diffusivity (and to a lesser extent the numerical diffusivity) taming whistler waves at high frequencies. Preliminary tests indicate that there may indeed be room to increase the time-step in the context of Hall-MHD turbulence. The preliminary results provided by a simple benchmark based on the OT vortex problem anticipate that our code will be able to deliver excellent performances with the simulation of astrophysical and space plasmas in which the Hall term is expected to play a significant role in 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; Herbert et al. 2016). 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 (Pouquet & Marino 2013; Marino et al. 2015a). The computational efficiency of our LB model 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, though in a simplified configuration, that our LB model is able to capture the physical effects produced by the Hall term, such as faster dynamics due to the interplay of whistler waves and turbulence, the breakdown of the Kolmogorov spectrum at sub-ion scales and a behavior of the magnetic energy spectrum at those scales which has been already observed in solar wind measurements. All that provides FLAME with the potential to become a powerful tool for the investigation of magnetohydrodynamic plasmas in a variety of configurations of interest for heliospheric and magnetospheric studies. While an incompressible (or weakly-compressible) formulation is usually justified in the context of a space plasma (Andrés et al. 2022; Brodiano et al. 2023), it would be worth accounting explicitly for compressible effects to reach a more comprehensive representation of its dynamics. Adding compressibility effects would require resorting to an additional equation of state and a coupled lattice scheme to deal with the temperature or density space-and-time evolutions. The present analysis and tests were performed using a benchmark configuration (OT vortex problem) which is isotropic, hence does not embed the anisotropy introduced by the ambient magnetic field in which the solar wind develops its dynamics. However, LB simulations of Hall-MHD flows performed with FLAME will be suitable to investigate plasmas immersed in a background magnetic field, at scales that are nowadays within the reach of the high-resolution instruments on board of the latest solar and magnetospheric missions, such as Solar Orbiter or Parker Solar Probe.

## Acknowledgements

We gratefully thank Dr. Emmanuel Quemener for the technical support provided during the development of our multi-GPU code at the Centre Blaise Pascal computer testing platform of the École Normale de Lyon (France). Most of the simulations were ran on HPC facilities at the École Centrale de Lyon (PMCS2I), in Ecully (France), that is
supported by the Auvergne-Rhône-Alpes region through the GRANT CPRT07-13 CIRA and the national Equip@Meso grant (ANR-10-EQPX-29-01). We acknowledge as well CINECA (Italy) that provided cpu time under the ISCRA initiative (to perform high-resolution simulations of Hall-MHD turbulence) as well as support in the frame of the project LaB-HMHD - HP10C4HXCB. R.M., R.F. and F.F. acknowledge support from the project "EVENTFUL" (ANR-20-CE30-0011), funded by the French "Agence Nationale de la Recherche" - ANR through the program AAPG-2020. The collaboration of R.F. and R.M. was facilitated by support from the International Space Science Institute in ISSI Team 556. We kindly acknowledge the two anonymous referees for the relevant and interesting remarks which helped to significantly improve the presentation of our results.

## Data availability

The data that support the findings of this study are available from the corresponding author, R.F., upon reasonable request.

## Competing interests

The authors report no conflict of interests.

## 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 {$c_i$}$_{i=0,...,26}$ is defined in Cartesian components by

$$|c_x\rangle = [0, -1, 0, 0, -1, -1, -1, -1, -1, -1, -1, -1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1]^T$$

$$|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]^T$$

$$|c_z\rangle = [0, 0, 0, -1, 0, 0, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, 0, 0, 1, 0, 0, 1, -1, 1, -1, 1, -1, 1, -1]^T$$

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

$$f_i^{(0)}(\rho, \mathbf{U}) = w_i \rho \left\{ 1 + \frac{\mathbf{c}_i \cdot \mathbf{U}}{c_s^2} + \frac{1}{2c_s^4} \left[ H_{ixx}^{(2)} U_x^2 + H_{iyy}^{(2)} U_y^2 + H_{izz}^{(2)} U_z^2 + 2 \left( H_{ixy}^{(2)} U_x U_y + \right. \right. \right. \\
+ H_{ixz}^{(2)} U_x U_z + H_{iyz}^{(2)} U_y U_z \left. \left. \right) \right] + \frac{1}{2c_s^6} \left[ H_{ixxy}^{(3)} U_x^2 U_y + H_{ixxz}^{(3)} U_x^2 U_z + H_{ixyy}^{(3)} U_x U_y^2 + \right. \right. \\
+ H_{ixzz}^{(3)} U_x U_z^2 + H_{iyyz}^{(3)} U_y U_z^2 + H_{iyyz}^{(3)} U_y^2 U_z + 2 H_{ixyz}^{(3)} U_x U_y U_z \left. \right] + \frac{1}{4c_s^8} \left[ H_{ixxyy}^{(4)} U_x^2 U_y^2 + \right. \right. \\
+ H_{ixzzz}^{(4)} U_x^2 U_z^2 + H_{iyyzz}^{(4)} U_y^2 U_z^2 + 2 \left( H_{ixyzz}^{(4)} U_x U_y U_z^2 + H_{ixyyx}^{(4)} U_x U_y^2 U_z + \right. \right. \\
+ H_{ixxyz}^{(4)} U_x^2 U_y U_z \left. \right) + \frac{1}{4c_s^{10}} \left[ H_{ixxyzz}^{(5)} U_x^2 U_y U_z^2 + H_{ixxyyz}^{(5)} U_x^2 U_y^2 U_z + \right. \right. \\
+ H_{ixyyzz}^{(5)} U_x U_y^2 U_z^2 \left. \right] + \frac{1}{8c_s^{12}} H_{ixxyyzz}^{(6)} U_x^2 U_y^2 U_z^2 \right\}$$

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

$$f_i^{\text{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})(\mathbf{c}_i \cdot \mathbf{c}_i) - (\mathbf{c}_i \cdot \mathbf{B})^2 \right).$$

The set of central moments $k_i$ is computed by applying the (invertible) transformation matrix $T$ with the column vectors

$$ 
\begin{align*}
|T_0\rangle &= |1\rangle \\
|T_1\rangle; |T_2\rangle; |T_3\rangle &= [\bar{c}_{ix}]^\top; [\bar{c}_{iy}]^\top; [\bar{c}_{iz}]^\top \\
|T_4\rangle; |T_5\rangle; |T_6\rangle &= [\bar{c}_{ix}\bar{c}_{iy}]^\top; [\bar{c}_{ix}\bar{c}_{iz}]^\top; [\bar{c}_{iy}\bar{c}_{iz}]^\top \\
|T_7\rangle; |T_8\rangle; |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 \\
|T_{10}\rangle; |T_{11}\rangle; |T_{12}\rangle &= [\bar{c}_{ix}\bar{c}_{iy} + \bar{c}_{ix}\bar{c}_{iz}^2]^\top; [\bar{c}_{ix}\bar{c}_{iy} + \bar{c}_{iy}\bar{c}_{iz}^2]^\top; [\bar{c}_{ix}^2\bar{c}_{iy} + \bar{c}_{iy}^2\bar{c}_{iz}^2]^\top \\
|T_{13}\rangle; |T_{14}\rangle; |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}^2]^\top \\
|T_{16}\rangle &= [\bar{c}_{ix}\bar{c}_{iy}\bar{c}_{iz}]^\top \\
|T_{17}\rangle; |T_{18}\rangle; |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 \\
|T_{20}\rangle; |T_{21}\rangle; |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 \\
|T_{23}\rangle; |T_{24}\rangle; |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}^2]^\top \\
|T_{26}\rangle &= [\bar{c}_{ix}^2\bar{c}_{iy}^2\bar{c}_{iz}^2]^\top
\end{align*}
$$

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

$$S = \text{diag}[1, 1, 1, 1, \omega, \omega, \omega, \omega, 1, \dots, 1],$$

which leads to

$$ 
\begin{align*}
k_{0 \dots 3}^* &= \langle T_{0 \dots 3} | f^{\text{mhd}(0)} \rangle \\
k_{4 \dots 8}^* &= \omega \langle T_{4 \dots 8} | f^{\text{mhd}(0)} \rangle + (1 - \omega) \langle T_{4 \dots 8} | \bar{f} \rangle \\
k_{9 \dots 26}^* &= \langle T_{9 \dots 26} | f^{\text{mhd}(0)} \rangle.
\end{align*}
$$

## Appendix B. Calculation of the electric current density

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

$$ \tilde{\mathbf{M}} \begin{bmatrix} J_x \\ J_y \\ J_z \end{bmatrix} = - \frac{1}{2\alpha_H} \begin{bmatrix} \Lambda_{yz} - \Lambda_{zy} - 2(U_y B_z - U_z B_y) \\ \Lambda_{zx} - \Lambda_{xz} - 2(U_z B_x - U_x B_z) \\ \Lambda_{xy} - \Lambda_{yx} - 2(U_x B_y - U_y B_x) \end{bmatrix} $$

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

$$ \tilde{\mathbf{M}} = \begin{bmatrix} C^2/2\alpha_H\omega_B & B_z & -B_y \\ -B_z & C^2/2\alpha_H\omega_B & B_x \\ 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 (3.20) 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

$$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)$$

$$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_x^2 - 4\alpha_H^2\omega_B^3 \Lambda_{zx} B_x^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)$$

$$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_x B_z - 4\alpha_H^2\omega_B^3 \Lambda_{zx} B_x B_z - 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) \quad (B1)$$

where $D = C^6 + 4C^2\alpha_H^2\omega_B^2|\mathbf{B}|^2 > 0$.

## REFERENCES

* ADHIKARI, L., ZANK, G. P., BRUNO, R., TELLONI, D., HUNANA, P., DOSCH, A., MARINO, R. & Hu, Q. 2015a The transport of low-frequency turbulence in astrophysical flows. ii. solutions for the super-alfvénic solar wind. *The Astrophysical Journal* **805** (1), 63.

* ADHIKARI, L., ZANK, G P, BRUNO, R., TELLONI, D., HUNANA, P., DOSCH, A., MARINO, R & Hu, Q 2015b The transport of low-frequency turbulence in the super-Alfvénic solar wind. *Journal of Physics: Conference Series* **642** (1), 012001.

* ALEXANDROVA, O., LACOMBE, C. & MANGENEY, A. 2008 Spectra and anisotropy of magnetic fluctuations in the Earth's magnetosheath: Cluster observations. *Ann. Geophysicae* **26** (11), 3585–3596.

* 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.

* ANDRÉS, N., SAHRAOUI, F., HUANG, S., HADID, L. Z. & GALTIER, S. 2022 The incompressible energy cascade rate in anisotropic solar wind turbulence. *Astron. & Astroph.* **661**, A116, arXiv: 2112.13748.

* 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. *Nat.* **576** (7786), 237–242.

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.

BENZI, R., SUCCI, S. & VERGASSOLA, M. 1992 The lattice Boltzmann equation: Theory and applications. *Phys. Rep.* **222** (3), 145–197.

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, 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.

BREYIANNIS, G. & VALOUGEORGIS, D. 2004 Lattice kinetic simulations in three-dimensional magnetohydrodynamics. *Phys. Rev. E* **69**, 065702.

BREYIANNIS, GEORGE & VALOUGEORGIS, DIMITRIS 2006 Lattice kinetic simulations of 3-D MHD turbulence. *Comp. & Fl.* **35**, 920–924.

BRODIANO, M., MITRUK, P. & ANDRÉS, N. 2023 A statistical study of the compressible energy cascade rate in solar wind turbulence: Parker solar probe observations. *Physics of Plasmas* **30** (3), 032903.

BRUNO, ROBERTO & CARBONE, VINCENZO 2013 The solar wind as a turbulence laboratory. *Living Rev. Sol. Phys.* **10** (1), 2.

BRUNO, R. & TRENCHI, L. 2014 Radial dependence of the frequency break between fluid and kinetic scales in the solar wind fluctuations. *Astroph. J. Lett.* **787** (2), L24.

BRUNO, R., TRENCHI, L. & TELLONI, D. 2014 Spectral slope variation at proton scales from fast to slow solar wind. *Astroph. J. Lett.* **793** (1), L15.

CERRI, S. S., CALIFANO, F., JENKO, F., TOLD, D. & RINCON, F. 2016 Subproton-scale cascades in solar wind turbulence: Driven hybrid-kinetic simulations. *The Astrophysical Journal Letters* **822** (1), L12.

CHEN, SHIYI, CHEN, HUDONG, MARTNEZ, DANIEL & MATTHAEUS, WILLIAM 1991 Lattice Boltzmann model for simulation of magnetohydrodynamics. *Phys. Rev. Lett.* **67**, 3776–3779.

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, 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, 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, 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. & Astroph.* **656**, A21.

DE ROSIS, ALESSANDRO 2017 Nonorthogonal central-moments-based lattice Boltzmann scheme in three dimensions. *Phys. Rev. E* **95**, 013310.

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, ALESSANDRO, LÉVÈQUE, EMMANUEL & CHAHINE, ROBERT 2018 Advanced lattice Boltzmann scheme for high-Reynolds-number magneto-hydrodynamic flows. *J. of Turb.* **19** (6), 446-462.

DELLAR, PAUL J. 2002 Lattice kinetic schemes for magnetohydrodynamics. *J. of Comp. Phys.* **179** (1), 95-126.

DELLAR, P J 2009 Moment equations for magnetohydrodynamics. *J. of Stat. Mech.: Th. and Exp. 2009* (06), P06003.

DELLAR, PAUL J. 2011 Lattice boltzmann formulation for braginskii magnetohydrodynamics. *Comp. & Fl.* **46** (1), 201-205, 10th ICFD Conference Series on Numerical Methods for Fluid Dynamics (ICFD 2010).

DELLAR, PAUL J. 2013 Lattice Boltzmann magnetohydrodynamics with current-dependent resistivity. *J. of Comp. Phys.* **237**, 115-131.

DUDSON, B. D., ALLEN, A., BREYIANNIS, G., BRUGGER, E., BUCHANAN, J., EASY, L., FARLEY, S., JOSEPH, I., KIM, M., MCGANN, A. D. & ET AL. 2015 Bout: Recent and current developments. *Journal of Plasma Physics* **81** (1), 365810104.

D'HUMIERES, D 1994 Generalized lattice Boltzmann equations. *Prog. Aeronaut. Astronaut.* **159**, 450-458.

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. *Europhys. Lett.* **123**, 44002.

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. *Astroph. J.* **927** (2), 205.

FLINT, CHRISTOPHER & VAHALA, GEORGE 2018 A partial entropic lattice Boltzmann MHD simulation of the Orszag-Tang vortex. *Rad. Eff. and Def. in Sol.* **173** (1-2), 55-65.

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, SÉBASTIEN 2016 *Introduction to modern magnetohydrodynamics*. Cambridge University Press.

GALTIER, SEBASTIEN & BUCHLIN, ERIC 2007 Multiscale Hall-magnetohydrodynamic turbulence in the solar wind. *Astroph. J.* **656** (1), 560-566.

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, MARTIN, GREINER, ANDREAS & KORVINK, JAN G. 2006 Cascaded digital lattice Boltzmann automata for high Reynolds number flow. *Phys. Rev. E* **73**, 066705.

GEIER, MARTIN, SCHÖNHERR, MARTIN, PASQUALI, ANDREA & KRAFCZYK, MANFRED 2015 The cumulant lattice Boltzmann equation in three dimensions: Theory and validation. *Comp. & Math. with App.* **70** (4), 507-547.

GÓMEZ, DANIEL O., MININNI, PABLO D. & DMITRUK, PABLO 2010 Hall-magnetohydrodynamic small-scale dynamos. *Phys. Rev. E* **82**, 036406.

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. *Astroph. J.* **870** (2), 94.

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.
HE, XIAOYI, SHAN, XIAOWEN & DOOLEN, GARY D. 1998 Discrete Boltzmann equation model for nonideal gases. *Phys. Rev. E* **57**, R13-R16.

HERBERT, C., MARINO, R., ROSENBERG, D. & POUQUET, A. 2016 Waves and vortices in the inverse cascade regime of stratified turbulence with or without rotation. *J. of Fl. Mech.* **806**, 165-204.

HIGUERA, F. J., SUCCI, S. & BENZI, R. 1989 Lattice gas dynamics with enhanced collisions. *Europhysics Letters* **9** (4), 345.

HOELZL, M., HUIJSMANS, G.T.A., PAMELA, S.J.P., BÉCOULET, M., NARDON, E., ARTOLA, F.J., NKONGA, B., ATANASIU, C.V., BANDARU, V., BHOLE, A., BONFIGLIO, D., CATHEY, A., CZARNY, O., DVORNOVA, A., FEHÉR, T., FIL, A., FRANCK, E., FUTATANI, S., GRUCA, M., GUILLARD, H., HAVERKORT, J.W., HOLOD, I., HU, D., KIM, S.K., KORVING, S.Q., KOS, L., KREBS, I., KRIPNER, L., LATU, G., LIU, F., MERKEL, P., MESHCHERIAKOV, D., MITTERAUER, V., MOCHALSKYY, S., MORALES, J.A., NIES, R., NIKULSIN, N., ORAIN, F., PRATT, J., RAMASAMY, R., RAMET, P., REUX, C., SÄRKIMÄKI, K., SCHWARZ, N., VERMA, P. SINGH, SMITH, S.F., SOMMARIVA, C., STRUMBERGER, E., VAN VUGT, D.C., VERBEEK, M., WESTERHOF, E., WIESCHOLLEK, F. & ZIELINSKI, J. 2021 The jorek non-linear extended mhd code and applications to large-scale instabilities and their control in magnetically confined fusion plasmas. *Nuclear Fusion* **61** (6), 065001.

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. & Astroph.* **642**, A9.

HORSTMANN, JAN, TOUIL, HATEM, VIENNE, LUCIEN, RICOT, DENIS & LÉVÈQUE, EMMANUEL 2022 Consistent time-step optimization in the lattice Boltzmann method. *J. of Comp. Phys.* **462**, 111224.

HUBA, J. D. 2003 Hall Magnetohydrodynamics - A tutorial. In *Space plasma simulation*, vol. 615, pp. 166-192. Büchner, J. and Dum, C. and Scholer, M.

HÉNON, MICHEL 1987 Viscosity of a lattice gas. *Compl. Sys.* **462**.

IROSHNIKOV, P. S. 1963 Turbulence of a conducting fluid in a strong magnetic field. *Astron. Zh.* **40**, 742.

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. *Philos. Trans. of the R. Soc. A: Math., Phys. and Eng. Sc.* **373** (2041), 20140155.

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, 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, ROBERT H. 1965 Inertial-range spectrum of hydromagnetic turbulence. *Phys. of Fl.* **8** (7), 1385-1387.

KRUEGER, TIMM, KUSUMAATMAJA, HALIM, KUZMIN, ALEXANDR, SHARDT, OREST, SILVA, GONCALO & VIGGEN, ERLEND MAGNUS 2016 *The lattice Boltzmann method: Principles and practice*. Springer.

KULSRUD, RUSSELL M. 2005 *Plasma Physics for Astrophysics*. Princeton University Press.

LEWY, H., FRIEDRICHS, K. & COURANT, R. 1928 Über die partiellen Differenzengleichungen der mathematischen Physik. *Mathematische annalen* **100**, 32-74.

MA, YINGJUAN, RUSSELL, CHRISTOPHER T., TOOTH, 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. *J. of Geoph. Res.: Sp. Phys.* **123** (5), 3742–3763.

MAHAJAN, S. M. & KRISHAN, V. 2005 Exact solution of the incompressible Hall magnetohydrodynamics. *Mon. Not. of the R. Astron. Soc.: Lett.* **359** (1), L27–L29.

MALARA, F. & VELLI, M. 1996 Parametric instability of a large-amplitude nonmonochromatic Alfvén wave. *Phys. of Pl.* **3** (12), 4427–4433.

MALASPINAS, ORESTIS 2015 Increasing stability and accuracy of the lattice Boltzmann scheme: Recursivity and regularization. *arXiv:1505.06900 preprint*.

MARCHAND, P., COMMERÇON, B. & CHABRIER, G. 2018 Impact of the Hall effect in star formation and the issue of angular momentum conservation. *Astron. & Astroph.* **619**, A37.

MARINO, R., MININNI, P. D., ROSENBERG, D. & POUQUET, A. 2013 Inverse cascades in rotating stratified turbulence: Fast growth of large scales. *Europhys. Lett.* **102** (4), 44006.

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, R., POUQUET, A. & ROSENBERG, D. 2015a Resolving the paradox of oceanic large-scale balance and small-scale mixing. *Phys. Rev. Lett.* **114**, 114504.

MARINO, RAFFAELE, ROSENBERG, DUANE, HERBERT, CORENTIN & POUQUET, ANNICK 2015b Interplay of waves and eddies in rotating stratified turbulence and the link with kinetic-potential energy partition. *Europhys. Lett.* **112**.

MARINO, RAFFAELE & SORRISO-VALVO, LUCA 2023 Scaling laws for the energy transfer in space plasma turbulence. *Phys. Rep.* **1006**, 1–144.

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. J. Lett.* **677** (1), L71.

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. *Plan. and Sp. Sc.* **59** (7), 592–597.

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. J.* **750** (1), 41.

MARTÍNEZ, DANIEL O., CHEN, SHIYI & MATTHAEUS, WILLIAM H. 1994 Lattice Boltzmann magnetohydrodynamics. *Phys. of Pl.* **1** (6), 1850–1867.

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. J. Lett.* **678** (2), L141.

MENDOZA, M. & MUÑOZ, J. D. 2008 Three-dimensional lattice boltzmann model for magnetic reconnection. *Phys. Rev. E* **77**, 026713.

MEYRAND, ROMAIN & GALTIER, SÉBASTIEN 2012 Spontaneous chiral symmetry breaking of Hall magnetohydrodynamic turbulence. *Phys. Rev. Lett.* **109**, 194501.

MININNI, PABLO D., GÓMEZ, DANIEL O. & MAHAJAN, SWADESH M. 2002 Dynamo action in Hall magnetohydrodynamics. *Astroph. J.* **567** (1), L81–L83.

MININNI, PABLO D., GOMEZ, DANIEL O. & MAHAJAN, SWADESH M. 2003 Dynamo action in magnetohydrodynamics and Hall-magnetohydrodynamics. *Astroph. J.* **587** (1), 472–481.

MININNI, PABLO D., GOMEZ, DANIEL O. & MAHAJAN, SWADESH M. 2005 Direct simulations of helical Hall-MHD turbulence and dynamo action. *Astroph. J.* **619** (2), 1019–1027.

MININNI, P. D., POUQUET, A. G. & MONTGOMERY, D. C. 2006 Small-scale structures in three-dimensional magnetohydrodynamic turbulence. *Phys. Rev. Lett.* **97**, 244503.

MININNI, PABLO D., ROSENBERG, DUANE, REDDY, RAGHU & POUQUET, ANNICK 2011 A hybrid MPI–OpenMP scheme for scalable parallel pseudospectral computations for fluid turbulence. *Parall. Comp.* **37** (6), 316–326.

MIURA, H. & ARAKI, K. 2014 Structure transitions induced by the Hall term in homogeneous and isotropic magnetohydrodynamic turbulence. *Phys. of Pl.* **21** (7), 072313.

MONTGOMERY, DAVID & DOOLEN, GARY D. 1987 Magnetohydrodynamic cellular automata. *Phys. Lett. A* **120**, 229–231.
MORALES, LAURA, DASSO, S. & GÓMEZ, DANIEL 2005 Hall effect in incompressible magnetic reconnection. *J. of Geoph. Res.* **110**.

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. 2020 The Solar Orbiter mission. Science overview. *Astron. Astroph.* **642**, A1.

NORMAN, C. & HEYVAERTS, J. 1985 Anomalous magnetic field diffusion during star formation. *Astron. & Astroph.* **147** (2), 247–256.

ORSZAG, S. A. & TANG, C. M. 1979 Small-scale structure of two-dimensional magnetohydrodynamic turbulence. *J. of Fl. Mech.* **90**, 129–143.

PANDEY, B. P. & WARDLE, MARK 2008 Hall magnetohydrodynamics of partially ionized plasmas. *Mon. Not. of the R. Astron. Soc.* **385** (4), 2269–2278.

PAPINI, EMANUELE, FRANCI, LUCA, LANDI, SIMONE, VERDINI, ANDREA, MATTEINI, LORENZO & HELLINGER, PETR 2019 Can Hall magnetohydrodynamics explain plasma turbulence at sub-ion scales? *Astroph. J.* **870** (1), 52.

PARASHAR, T. N., CUESTA, M. & MATTHAEUS, W. H. 2019 Reynolds number and intermittency in the expanding solar wind: Predictions based on Voyager observations. *Astroph. J. Lett.* **884** (2), L57.

PATTERSON, G. S. & ORSZAG, STEVEN A. 1971 Spectral calculations of isotropic turbulence: Efficient removal of aliasing interactions. *Phys. of Fl.* **14** (11), 2538–2541.

PATTISON, M.J., PREMNATH, K.N., MORLEY, N.B. & ABDOU, M.A. 2008 Progress in lattice boltzmann methods for magnetohydrodynamic flows relevant to fusion applications. *Fusion Engineering and Design* **83** (4), 557–572.

POUQUET, A. & MARINO, R. 2013 Geophysical turbulence and the duality of the energy flow across scales. *Phys. Rev. Lett.* **111**, 234501.

POUQUET, A., ROSENBERG, D., STAWARZ, J. & MARINO, R. 2019 Helicity dynamics, inverse, and bidirectional cascades in fluid and magnetohydrodynamic turbulence: A brief review. *Earth Space Sci.* **6**, 1–19.

RILEY, B., RICHARD, J. & GIRIMAJI, S. S 2008 Progress in lattice boltzmann methods for magnetohydrodynamic schemes in turbulence and rectangular jets. *Int. J. Mod. Phys. C* **19**, 1211–1220.

ROSENBERG, D., MININNI, P.D., REDDY, R. & POUQUET, A. 2020 GPU parallelization of a hybrid pseudospectral geophysical turbulence framework using CUDA. *Atm.* **11**, 178.

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, XIAOWEN & HE, XIAOYI 1998 Discretization of the velocity space in the solution of the Boltzmann equation. *Phys. Rev. Lett.* **80**, 65–68.

SHEN, NAIJIAN, LI, YUAN, PULLIN, D. I., SAMTANEY, RAVI & WHEATLEY, VINCENT 2018 On the magnetohydrodynamic limits of the ideal two-fluid plasma equations. *Physics of Plasmas* **25** (12), 122113.

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. & Astroph.* **650**, A21.

SILVA, GONCALO & SEMIAO, VIRIATO 2014 Truncation errors and the rotational invariance of three-dimensional lattice models in the lattice Boltzmann method. *J. of Comp. Phys.* **269**, 259–279.

SMITH, CHARLES W., HAMILTON, KATHLEEN, VASQUEZ, BERNARD J. & LEAMON, ROBERT J. 2006 Dependence of the dissipation range spectrum of interplanetary magnetic fluctuations on the rate of energy cascade. *Astroph. J. Lett.* **645** (1), L85–L88.
SORRISO-VALVO, L., MARINO, R., FOLDES, R., LÉVÈQUE, E., D'AMICIS, R., BRUNO, R., TELLONI, D. & YORDANOVA, E. 2023 Helios 2 observations of solar wind turbulence decay in the inner heliosphere. *Astron. & Astroph.* **672**, A13.

SOVINEC, C.R. & KING, J.R. 2010 Analysis of a mixed semi-implicit/implicit algorithm for low-frequency two-fluid plasma modeling. *Journal of Computational Physics* **229** (16), 5803–5819.

SUCCI, S., VERGASSOLA, M. & BENZI, R. 1991 Lattice Boltzmann scheme for two-dimensional magnetohydrodynamics. *Phys. Rev. A* **43**, 4521–4524.

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. *Astroph. J. Lett.* **938** (2), L8.

TELLONI, DANIELE, CARBONE, FRANCESCO, BRUNO, ROBERTO, ZANK, GARY P., SORRISO-VALVO, LUCA & MANCUSO, SALVATORE 2019 Ion cyclotron waves in field-aligned solar wind turbulence. *Astroph. J. Lett.* **885** (1), L5.

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. *Astroph. J. Lett.* **912** (2), L21.

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. *Astroph. J. Lett.* **936** (2), L25.

TÓTH, GÁBOR, MA, YINGJUAN & GOMBOSI, TAMAS I. 2008 Hall magnetohydrodynamics on block-adaptive grids. *J. of Comp. Phys.* **227** (14), 6967–6984.

WANG, XIAOGANG, BHATTACHARJEE, A. & MA, Z. W. 2001 Scaling of collisionless forced reconnection. *Phys. Rev. Lett.* **87**, 265003.
XIA, ZHENWEI & YANG, WEIHONG 2015 Exact solutions of the incompressible dissipative Hall magnetohydrodynamics. *Phys. of Plas.* **22** (3), 032306.

YADAV, SHARAD K., MIURA, HIDEAKI & PANDIT, RAHUL 2022 Statistical properties of three-dimensional Hall magnetohydrodynamics turbulence. *Phys. of Fl.* **34** (9), 095135.

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