Efficient kinetic Lattice Boltzmann simulation of three-dimensional Hall-MHD turbulence

Abstract

Simulating plasmas in the Hall-magnetohydrodynamics (Hall-MHD) regime represents a valuable approach for the investigation of complex nonlinear dynamics developing in astrophysical frameworks and fusion machines. The Hall electric field is computationally very challenging as it involves the integration of an additional term, proportional to $\boldsymbol {\nabla } \times ((\boldsymbol {\nabla }\times \boldsymbol {B})\times \boldsymbol {B})$, in 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 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 magnetohydrodynamics (MHD) model, plasma is treated as a single species quasi-neutral fluid with conductive properties sensitive to the action of the magnetic field (Galtier 2016). In the ideal MHD description, ions and electrons are tied to the magnetic field, moving with the same velocity. The Hall-MHD model relaxes the MHD prescriptions assuming ions disunite from the magnetic field due to their inertia, while electrons remain bound to it (Pandey & Wardle 2008). In this framework, the resistive Ohm's law is generalized through the introduction of the Hall electric field, proportional to $\boldsymbol {J} \times \boldsymbol {B}$, where $\boldsymbol {J}$ and $\boldsymbol {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 the 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, Commerçon & Chabrier 2018) to the solar atmosphere and the solar wind (Galtier & Buchlin 2007; González-Morales, Khomenko & Cally 2019), and it has been used also to investigate magnetic reconnection processes (Wang, Bhattacharjee & Ma 2001; Morales, Dasso & Gómez 2005; Ma et al. 2018) and the dynamo action (Mininni, Gómez & Mahajan 2002; Mininni, Gomez & Mahajan 2005; Gómez, Mininni & Dmitruk 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$. 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, Gomez & Mahajan 2003), finite-volume (Tóth, Ma & Gombosi 2008; Marchand 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). However, 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, Gross & Krook 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, Vergassola & Benzi 1991; Martínez, Chen & Matthaeus 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, Succi & Benzi 1989; Benzi, Succi & Vergassola 1992; d'Humieres 1994). Therefore, Pattison et al. (2008) and Riley, Richard & Girimaji (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 modelling to simulate Hall-MHD turbulence, an effort that has never been undertaken previously. In the present study, an 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 by De Rosis, Lévêque & Chahine (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. Section 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 graphics processing unit (GPU)-accelerated simulations of the three-dimensional Orszag–Tang (OT) 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:

(2.1)\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{U} = 0, \end{gather}

(2.2)\begin{gather} \partial_t \boldsymbol{U}+(\boldsymbol{U} \boldsymbol{\cdot} \boldsymbol{\nabla}) \boldsymbol{U}=\frac{1}{\rho_0}\boldsymbol{J}\times \boldsymbol{B} -\frac{1}{\rho_0}\boldsymbol{\nabla} {p} + \nu \boldsymbol{\nabla}^2 \boldsymbol{U}, \end{gather}

(2.3)\begin{gather} \partial_t \boldsymbol{B}=\boldsymbol{\nabla} \times (\boldsymbol{U}\times \boldsymbol{B} - \eta \boldsymbol{\nabla} \times \boldsymbol{B}), \end{gather}

(2.4)\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{B}=0, \end{gather}

where $t$ is the time, $\rho _0$ is the mass density of the fluid, $\nu$ is the kinematic viscosity and $\eta$ it the magnetic resistivity. The electric current density is expressed as $\boldsymbol {J}={1}/{\mu _0}\boldsymbol {\nabla } \times \boldsymbol {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

(2.5)\begin{gather} \rho[\partial_t \boldsymbol{U} + (\boldsymbol{U} \boldsymbol{\cdot} \boldsymbol{\nabla})\boldsymbol{U}] = en(\boldsymbol{E}+\boldsymbol{U}\times \boldsymbol{B})-\boldsymbol{\nabla} {p}_i + \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol \sigma + \boldsymbol{R}, \end{gather}

(2.6)\begin{gather} 0 =-en(\boldsymbol{E}+\boldsymbol{U_e}\times \boldsymbol{B})-\boldsymbol{\nabla} {p}_e - \boldsymbol{R}, \end{gather}

where $e$ is the unit electric charge, $\boldsymbol \sigma$ is the viscous stress tensor, $n$ is the particle density with $\rho =mn$, and $\boldsymbol {R}$ is the rate (per unit volume) of momentum exchange due to collisions between protons and electrons. The latter is given by $\boldsymbol {R}=-mn \,f_{ie}(\boldsymbol {U}-\boldsymbol {U_{e}})$, where $f_{ie}$ denotes the collision frequency and can be reformulated as $\boldsymbol {R}=-(m \,f_{ie}/e) \boldsymbol {J}$, with the density current $\boldsymbol {J}=en(\boldsymbol {U}-\boldsymbol {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

(2.7)\begin{equation} \partial_t \boldsymbol{U}+(\boldsymbol{U} \boldsymbol{\cdot} \boldsymbol{\nabla}) \boldsymbol{U}=\frac{1}{\rho}\boldsymbol{J}\times \boldsymbol{B}- \frac{1}{\rho}\boldsymbol{\nabla} {p} + \nu \boldsymbol{\nabla}^2 \boldsymbol{U}. \end{equation}

However, by replacing $\boldsymbol {U_e}$ by $\boldsymbol {U}-\boldsymbol {J}/ne$ and the expression for the rate of momentum exchange into (2.6), the Ohm's law becomes

(2.8)\begin{equation} \boldsymbol{E}=-\left(\boldsymbol{U}-\frac{1}{en}\boldsymbol{J}\right)\times \boldsymbol{B}+\frac{1}{en}\boldsymbol{\nabla} {p}_e + \frac{m \,f_{ie}}{e^2 n}\boldsymbol{J}. \end{equation}

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

(2.9)\begin{equation} \partial_t \boldsymbol{B}=\boldsymbol{\nabla} \times [(\boldsymbol{U}-\alpha_H\boldsymbol{J}) \times \boldsymbol{B}] + \eta \boldsymbol{\nabla}^2 \boldsymbol{B}, \end{equation}

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, $\boldsymbol {\nabla } \times (({1}/{en)}\boldsymbol {\nabla } {p}_e) = -({1}/{en^2})\boldsymbol {\nabla } n \times \boldsymbol {\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 = n T_e$, where $T_e$ is a constant plasma temperature. Therefore, $\boldsymbol {\nabla } n \times \boldsymbol {\nabla } p_e = 0$. The Hall-MHD equations mentioned earlier include a finite ion–electron collision frequency, responsible for the $\boldsymbol {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, however, 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 the electric current and electric field, and between the 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 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

(2.10)\begin{gather} \partial_t \rho + \boldsymbol{\nabla} \boldsymbol{\cdot} (\boldsymbol{\rho \boldsymbol{U}}) = 0, \end{gather}

(2.11)\begin{gather} \partial_t (\boldsymbol{\rho U}) + \boldsymbol{\nabla} \boldsymbol{\cdot} \left( \rho \boldsymbol{U} \otimes \boldsymbol{U} + {p} \mathbb{I} +\tfrac{1}{2}|\boldsymbol{B}|^2\mathbb{I}-\boldsymbol{B}\otimes \boldsymbol{B}\right) =\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\sigma}, \end{gather}

in which $\boldsymbol {\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 }=\tfrac {1}{2}|\boldsymbol {B}|^2\delta _{\alpha \beta }-B_{\alpha }B_{\beta }$¹ is adopted, and $\mu _0$ has been absorbed by replacing $\boldsymbol {B}$ with $\mu _0^{1/2}\boldsymbol {B}$. This (standard) normalization will be assumed hereafter, which allows for simplifying the Lorentz force as $(\boldsymbol {\nabla } \times \boldsymbol {B})\times \boldsymbol {B}$. The general form of the viscous stress is

(2.12)\begin{equation} \boldsymbol{\sigma} = \mu ( \boldsymbol{\nabla} \boldsymbol U + (\boldsymbol{\nabla} \boldsymbol U)^{\rm T} ) +\left(\zeta - \tfrac{2}{3} \mu\right) (\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol U) \mathbb{I}, \end{equation}

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

(2.13)\begin{equation} {p}=\rho c_s^2 \end{equation}

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

(2.14)\begin{equation} \partial_t \boldsymbol{B}+\boldsymbol{\nabla} \boldsymbol{\cdot} \left((\boldsymbol{U}-\alpha_H \boldsymbol{J}) \otimes \boldsymbol{B}-\boldsymbol{B} \otimes (\boldsymbol{U}-\alpha_H \boldsymbol{J}) \right) = \eta \boldsymbol{\nabla}^2 \boldsymbol{B}. \end{equation}

Let us remark that following the normalization of $\boldsymbol {B}$ by $\mu _0^{1/2}$, the Hall current $\alpha _H \boldsymbol {J}$ reads as ${\alpha _H}/{{\mu _0}^{1/2}}~\boldsymbol {\nabla } \times \boldsymbol {B}$. In the next sections, the developed LB scheme will conform to the set of (2.10), (2.11), (2.13) and (2.14). The divergence-free condition on $\boldsymbol {B}$ is preserved by (2.14), justifying that it is sufficient to impose $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {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 Appendix A. A central-moment collision operator (De Rosis et al. 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 behaviour 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 \boldsymbol {U}$ or stress tensor $\boldsymbol {\varSigma }$ are obtained as statistical moments of the particle distributions, i.e.

(3.1)\begin{gather} \rho = \sum_{i=0}^{N-1} f_i, \end{gather}

(3.2)\begin{gather}\rho \boldsymbol{U} = \sum_{i=0}^{N-1} f_i \boldsymbol{c}_i, \end{gather}

(3.3)\begin{gather}\boldsymbol{\varSigma} =\sum_{i=0}^{N-1} f_i \boldsymbol{c}_i \otimes \boldsymbol{c}_i, \end{gather}

by summing over the local mass densities $f_0,\ldots,\,f_{N-1}$ of particles moving with velocities $\boldsymbol {c}_0,\ldots,\boldsymbol {c}_{N-1}$, respectively. The sums replace here the integrals over $\boldsymbol {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}(\boldsymbol {x}, t)$ evolve according to a discrete-velocity analogue of the Boltzmann equation, which reads as

(3.4)\begin{equation} \partial_t f_i + \left( \boldsymbol c_i \boldsymbol{\cdot} \boldsymbol{\nabla}\right) f_i=- \frac{1}{\tau}\left(\,f_i -f_i^{(0)}(\rho, \boldsymbol{U})\right) \end{equation}

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, \boldsymbol {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 $\boldsymbol {c}_0,\ldots,\boldsymbol {c}_{N-1}$ chosen in a way such that particles travel from a lattice node to a neighbour lattice node in exactly one time step (see figure 1).

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

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

(3.5)\begin{equation} \bar{f}_i = f_i + {\Delta t}/{2\tau} (\,f_i-f_i^{(0)}) \end{equation}

originally introduced by He, Shan & Doolen (1998), as

(3.6)\begin{equation} \bar{f}_i(\boldsymbol{x}+\boldsymbol{c}_i \Delta t,t+\Delta t) = \bar{f}_i(\boldsymbol{x},t) - \omega\left(\bar{f}_i(\boldsymbol{x},t)-{f}^{(0)}_i(\rho,\boldsymbol{U})(\boldsymbol{x},t)\right), \end{equation}

where the discrete distribution functions $\bar {f}_i(\boldsymbol {x},t)$ depend on the three spatial coordinates $\boldsymbol {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

(3.7)\begin{equation} \frac{1}{\omega}=\left(\frac{\nu}{c_s^2 \Delta t}+\frac{1}{2}\right), \end{equation}

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

(3.8)\begin{align} \rho & = \sum_{i=0}^{N-1} \bar f_i \end{align}

(3.9)\begin{align} \mathrm{and}\quad \rho \boldsymbol{U} & = \sum_{i=0}^{N-1} \bar f_i \boldsymbol{c}_i. \end{align}

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

(3.10)\begin{gather} \bar{f}_i(\boldsymbol{x}+\boldsymbol{c}_i\Delta t,t+\Delta t) = \bar{f}_i^*(\boldsymbol{x},t), \end{gather}

(3.11)\begin{gather}\bar{f}^*_i(\boldsymbol{x},t) = \bar{f}_i(\boldsymbol{x},t) - \omega \left( \bar{f}_i(\boldsymbol{x},t)-{f}^{(0)}_i(\rho,\boldsymbol{U})(\boldsymbol{x},t) \right ). \end{gather}

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 $\boldsymbol {c_i}$, which reads as

(3.12)\begin{equation} f_i^{(0)}(\rho, \boldsymbol{U}) = w_i \rho \left( 1 +\frac{\boldsymbol{c}_i\boldsymbol{\cdot} \boldsymbol{U}}{c_s^2} +\frac{(\boldsymbol{c}_i\boldsymbol{\cdot} \boldsymbol{U})^2}{2c_s^4} - \frac{\boldsymbol{U}\boldsymbol{\cdot} \boldsymbol{U}}{2c_s^2} + \cdots \right) \end{equation}

with the weights $w_\mathrm {center}=8/27$, $w_\mathrm {face}=2/27$, $w_\mathrm {edge}=1/54$ and $w_\mathrm {corner}=1/216$ for the D3Q27 lattice (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; Coreixas, Chopard & Latt 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

(3.13)\begin{equation} {f}_i^{\mathrm{mhd}(0)}(\rho, \boldsymbol{U}, \boldsymbol{B}) = {f}_i^{(0)}(\rho, \boldsymbol{U}) + \frac{w_i}{2c_s^4} \left((\boldsymbol{B}\boldsymbol{\cdot} \boldsymbol B) (\boldsymbol{c}_i\boldsymbol{\cdot}\boldsymbol{c}_i)-(\boldsymbol{c}_i\boldsymbol{\cdot}\boldsymbol{B})^2\right) \end{equation}

so that the second-order moment becomes

(3.14)\begin{equation} \boldsymbol{\varSigma}^{\mathrm{mhd}(0)} = \sum_{i=0}^{N-1} f_i^{\mathrm{mhd}(0)}\boldsymbol{c}_i \otimes \boldsymbol{c}_i = \rho \boldsymbol{U} \otimes \boldsymbol{U} + {p} \mathbb{I} +\tfrac{1}{2}|\boldsymbol{B}|^2\mathbb{I}-\boldsymbol{B}\otimes \boldsymbol{B}. \end{equation}

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 artefacts 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, Greiner & Korvink 2006, 2007; Geier et al. 2015; De Rosis et al. 2018). This is the very framework adopted in lay-outing our code (details are given in 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 $\bar {\boldsymbol {c}}_i=\boldsymbol {c}_i-\boldsymbol {U}$ used for the CMs evaluation. Here, we consider the set of CMs as formally defined by

(3.15)\begin{equation} |{k}\rangle \equiv [{k}_0 \cdots{k}_{26}]^\top=\mathrm{T}^\top|{\bar f}\rangle, \end{equation}

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

(3.16)\begin{align} \left. \begin{gathered} |\mathrm{T}_0\rangle=|1\rangle, \\ |\mathrm{T}_1\rangle;|\mathrm{T}_2\rangle;|\mathrm{T}_3\rangle=[\bar{c}_{ix}]^\top;[\bar{c}_{iy}]^\top;[\bar{c}_{iz}]^\top,\\ |\mathrm{T}_4\rangle;|\mathrm{T}_5\rangle;|\mathrm{T}_6\rangle=[\bar{c}_{ix}\bar{c}_{iy}]^\top;[\bar{c}_{ix}\bar{c}_{iz}]^\top;[\bar{c}_{iy}\bar{c}_{iz}]^\top,\\ |\mathrm{T}_7\rangle;|\mathrm{T}_8\rangle;|\mathrm{T}_9\rangle=[\bar{c}_{ix}^2-\bar{c}_{iy}^2]^\top;[\bar{c}_{ix}^2-\bar{c}_{iz}^2]^\top;[\bar{c}_{ix}^2+\bar{c}_{iy}^2+\bar{c}_{iz}^2]^\top,\\ |\mathrm{T}_{10}\rangle;|\mathrm{T}_{11}\rangle;|\mathrm{T}_{12}\rangle=[\bar{c}_{ix}\bar{c}_{iy}^2+\bar{c}_{ix}\bar{c}_{iz}^2]^\top;[\bar{c}_{ix}\bar{c}_{iy}^2+\bar{c}_{iy}\bar{c}_{iz}^2]^\top;[\bar{c}_{ix}^2\bar{c}_{iy}+\bar{c}_{iy}^2\bar{c}_{iz}]^\top,\\ |\mathrm{T}_{13}\rangle;|\mathrm{T}_{14}\rangle;|\mathrm{T}_{15}\rangle=[\bar{c}_{ix}\bar{c}_{iy}^2-\bar{c}_{ix}\bar{c}_{iz}^2]^\top;[\bar{c}_{ix}\bar{c}_{iy}^2-\bar{c}_{iy}\bar{c}_{iz}^2]^\top;[\bar{c}_{ix}^2\bar{c}_{iy}-\bar{c}_{iy}^2\bar{c}_{iz}]^\top,\\ |\mathrm{T}_{16}\rangle=[\bar{c}_{ix}\bar{c}_{iy}\bar{c}_{iz}]^\top ,\\ |\mathrm{T}_{17}\rangle;|\mathrm{T}_{18}\rangle;|\mathrm{T}_{19}\rangle=[\bar{c}_{ix}^2\bar{c}_{iy}^2+\bar{c}_{ix}^2\bar{c}_{iz}^2+\bar{c}_{iy}^2\bar{c}_{iz}^2]^\top;[\bar{c}_{ix}^2\bar{c}_{iy}^2+\bar{c}_{ix}^2\bar{c}_{iz}^2-\bar{c}_{iy}^2\bar{c}_{iz}^2]^\top;[\bar{c}_{ix}^2\bar{c}_{iy}^2-\bar{c}_{ix}^2\bar{c}_{iz}^2]^\top,\\ |\mathrm{T}_{20}\rangle;|\mathrm{T}_{21}\rangle;|\mathrm{T}_{22}\rangle=[\bar{c}_{ix}^2\bar{c}_{iy}\bar{c}_{iz}]^\top;[\bar{c}_{ix}\bar{c}_{iy}^2\bar{c}_{iz}]^\top;[\bar{c}_{ix}\bar{c}_{iy}\bar{c}_{iz}^2]^\top,\\ |\mathrm{T}_{23}\rangle;|\mathrm{T}_{24}\rangle;|\mathrm{T}_{25}\rangle=[\bar{c}_{ix}\bar{c}_{iy}^2\bar{c}_{iz}^2]^\top;[\bar{c}_{ix}^2\bar{c}_{iy}\bar{c}_{iz}^2]^\top;[\bar{c}_{ix}^2\bar{c}_{iy}^2\bar{c}_{iz}]^\top,\\ |\mathrm{T}_{26}\rangle=[\bar{c}_{ix}^2\bar{c}_{iy}^2\bar{c}_{iz}^2]^\top. \end{gathered} \right\} \end{align}

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

(3.17)\begin{equation} |{k}^*\rangle =|{k}\rangle -S \left(|{k}\rangle-|{k}^{(0)}\rangle\right) \quad \mathrm{with}\ |{k}^{(0)}\rangle =\mathrm{T}^\top|{f^{\mathrm{mhd}(0)}}\rangle, \end{equation}

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

(3.18)\begin{equation} S=\mathrm{diag}[1,1,1,1,\omega,\omega,\omega,\omega,\omega,1,\ldots,1], \end{equation}

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

(3.19)\begin{equation} |\bar f^*\rangle = {\mathrm{T}^{-1}}^\top|{k}^*\rangle \end{equation}

before moving on to the streaming step (3.10).

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, Khanfir & Chanteu (1995) and Bouchut (1999), Dellar (2002) proposed a decomposition of the magnetic field as

(3.20)\begin{equation} \boldsymbol{B}(\boldsymbol{x},t)=\sum_{i=0}^{M-1} \boldsymbol{\bar g}_i(\boldsymbol{x},t), \end{equation}

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

The magnetic field $\boldsymbol {B}$ is here provided by the zeroth-order moment of $\boldsymbol {\bar 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 figure 1) shall prove to be satisfactory in reproducing the magnetic field of Hall-MHD turbulent plasmas. Analogously to the fluid case, an LB scheme can be derived to simulate the induction equation in the form

(3.21)\begin{equation} \boldsymbol{\bar g}_i(\boldsymbol{x}+\boldsymbol{\xi}_i\Delta t,t+\Delta t) = \boldsymbol{\bar g}_i(\boldsymbol{x},t) - \omega_B\left(\boldsymbol{\bar g}_i(\boldsymbol{x},t)-\boldsymbol{g}^{(0)}_i(\boldsymbol{U},\boldsymbol{B})(\boldsymbol{x},t)\right), \end{equation}

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

(3.22)\begin{equation} \frac{1}{\omega_B}=\left(\frac{\eta}{C^2\Delta t}+\frac{1}{2}\right) \end{equation}

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 $\boldsymbol {u}$, $\boldsymbol {B}$ or $\boldsymbol {J}$ may be exchanged between the two lattices without interpolation. This constraint imposes that

(3.23)\begin{equation} 2 C = \sqrt{3} c_s. \end{equation}

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

(3.24)\begin{equation} g_{i\alpha}^{(0)} (\boldsymbol{U}, \boldsymbol{B})=W_i\left(B_{\alpha}+\frac{1}{C^2}\xi_{i\beta}(U_{\beta}B_{\alpha}-B_{\beta}U_{\alpha})\right) \end{equation}

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

(3.25)\begin{equation} \sum_{i=0}^{M-1}\xi_{i\alpha} g^{(0)}_{i\beta}=\boldsymbol{U} \otimes \boldsymbol{B} - \boldsymbol{B}\otimes \boldsymbol{U} \end{equation}

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

(3.26)\begin{equation} \boldsymbol{{\varLambda}}^{ (0)}_{\alpha\beta}=\sum_{i=0}^{M-1}\xi_{i\alpha}g^{\mathrm{Hall} (0)}_{i\beta}=(\boldsymbol{U} - \alpha_H \boldsymbol{J}) \otimes \boldsymbol{B} - \boldsymbol{B} \otimes (\boldsymbol{U} - \alpha_H \boldsymbol{J}), \end{equation}

which is obviously possible by now considering

(3.27)\begin{equation} g_{i\alpha}^{\mathrm{Hall} (0)} (\boldsymbol{U}, \boldsymbol{B}, \boldsymbol{J}) = W_i\left(B_{\alpha}+\frac{1}{C^2}\xi_{i\beta}((U_{\beta} - \alpha_H J_\beta)B_{\alpha}-B_{\beta}(U_{\alpha}-\alpha_H J_\alpha))\right). \end{equation}

Nevertheless, $\boldsymbol {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 $\boldsymbol {U}$, they have only been extended up to the second order in $\boldsymbol {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 $\boldsymbol {J}$, are self-consistently obtained (within an $O({Ma}^3)$ error) from the first-order moment of the densities as

(3.28)\begin{equation} \mathrm{J}_{\gamma} = \varepsilon_{\alpha\beta\gamma} \frac{\partial B_{\alpha}}{\partial x_{\beta}}=-\varepsilon_{\alpha\beta\gamma} \frac{\omega_B}{C^2}\left(\boldsymbol{\boldsymbol{\varLambda}}_{\alpha\beta}-\boldsymbol{{\varLambda}}_{\alpha\beta}^{(0)}\right), \end{equation}

where $\varepsilon _{\alpha \beta \gamma }$ is the Levi–Civita tensor and $\boldsymbol {{\varLambda }}_{\alpha \beta }=\sum _{i=0}^{M-1}\xi _{i\alpha }\bar g_{i\beta }$ (Dellar 2002).

By replacing (3.26) in (3.28), we obtain a linear system readily solvable to obtain the current density $\boldsymbol {J}$, namely

(3.29)\begin{equation} \left(\mathbb{I}+\frac{2\alpha_H\omega_B}{C^2}\boldsymbol{M}\right)\boldsymbol{J}= \boldsymbol{J_0}, \end{equation}

where

(3.30a,b)\begin{equation} \boldsymbol{M}=\begin{bmatrix} 0 & B_z & -B_y \\ -B_z & 0 & B_x \\ B_y & -B_x & 0 \end{bmatrix}\quad \mathrm{and}\quad \boldsymbol{J_0}=\begin{bmatrix} \boldsymbol{{\varLambda}}_{yz}-\boldsymbol{{\varLambda}}_{zy} - 2\left(U_yB_z-U_zB_y\right) \\ \boldsymbol{{\varLambda}}_{zx}-\boldsymbol{{\varLambda}}_{xz} - 2\left(U_zB_x-U_xB_z\right) \\ \boldsymbol{{\varLambda}}_{xy}-\boldsymbol{{\varLambda}}_{yx} - 2\left(U_xB_y-U_yB_x\right) \end{bmatrix}. \end{equation}

Obviously, a solution exists only if it is possible to invert $\tilde {\boldsymbol {M}}=\mathbb {I}+(2\alpha _H\omega _B/C^2)\boldsymbol {M}$. It can be easily verified that $\mathrm {det}(\tilde {\boldsymbol {M}})\neq 0$, which proves this solution exists and is unique. The current density obtained by solving (3.29) can then be used to compute the equilibrium densities (3.27) 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.28) also provides a consistent approximation of the divergence of the magnetic field. Indeed by taking the trace of the magnetic tensor, one obtains

(3.31)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{B} \simeq-\frac{\omega_B}{C^2}\mathrm{Tr}\left(\boldsymbol{\varLambda}_{\alpha\beta}\right) \end{equation}

by noticing that $\mathrm {Tr}(\boldsymbol {\varLambda }_{\alpha \beta }^{(0)})=0$. Furthermore, the ${O}({Ma}^3)$ correction cancels out by taking the trace. Therefore, this correction is pushed to a higher order, so that the divergence-free $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {B}=0$ corresponds with high accuracy to the condition $\mathrm {Tr}(\boldsymbol {\varLambda }_{\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/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. To obtain a dimensionless induction equation, let us normalize 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

(3.32)\begin{equation} \left(\frac{U_0 B_0}{L_0}\right)\partial_t \boldsymbol{b}=\frac{1}{L_0}\boldsymbol{\nabla} \times \left[\left(U_0\boldsymbol{u}-\frac{\alpha_H B_0}{\sqrt{\mu_0} L_0}\boldsymbol{\nabla} \times \boldsymbol{b}\right)\times (B_0\boldsymbol{b})\right] + \frac{\eta B_0}{L_0^2} \boldsymbol{\nabla}^2 \boldsymbol{b}. \end{equation}

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

(3.33)\begin{equation} \partial_t \boldsymbol{b}=\boldsymbol{\nabla} \times \left[\left(\boldsymbol{u}-\epsilon_H\boldsymbol{\nabla} \times \boldsymbol{b}\right)\times \boldsymbol{b}\right] + \frac{1}{{Re}_m} \boldsymbol{\nabla}^2 \boldsymbol{b} \end{equation}

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

(3.34)\begin{equation} \epsilon_H = \frac{\alpha_H B_0}{\sqrt{\mu_0} L_0 U_0}. \end{equation}

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,

(3.35)\begin{equation} \rho \frac{U_0^2}{L_0}\left[\partial_t \boldsymbol{u}+\left(\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla}\right)\boldsymbol{u}\right] =-\frac{1}{L_0}\boldsymbol{\nabla} \rho c_s^2 + \rho\nu \frac{U_0}{L_0^2}\boldsymbol{\nabla}^2 \boldsymbol{u}+\frac{B_0^2}{L_0}\left(\boldsymbol{\nabla}\times \boldsymbol{b}\right)\times\boldsymbol{b}, \end{equation}

which gives

(3.36)\begin{equation} \partial_t \boldsymbol{u}+\left(\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla}\right)\boldsymbol{u} =-\frac{1}{{Ma}^2} \frac{1}{\rho} {\boldsymbol{\nabla} \rho} + \frac{1}{{Re}}\boldsymbol{\nabla}^2 \boldsymbol{u}+ \left( \frac{V_A}{U_0} \right)^2 \left (\boldsymbol{\nabla}\times \boldsymbol{b} \right )\times\boldsymbol{b}, \end{equation}

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

(3.37)\begin{equation} \epsilon_H^\mathrm{lbm} = \frac{\alpha_H}{\sqrt{\mu_0}} \frac{\left [ B_0/B^\mathrm{lbm}_0 \right ]}{\left [ L_0/L^\mathrm{lbm}_0 \right ] \left [U_0/U^\mathrm{lbm}_0 \right ]}=\epsilon_H \frac{U^\mathrm{lbm}_0 L_0^\mathrm{lbm}}{B_0^\mathrm{lbm}}. \end{equation}

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

(3.38)\begin{equation} \epsilon_H^\mathrm{lbm} = \epsilon_H L_0^\mathrm{lbm} = \epsilon_H N \end{equation}

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

(3.39)\begin{equation} \epsilon_H = \frac{V_A}{U_0L_0}\sqrt{\frac{m}{\mu_0ne^2}}=\frac{L_H}{L_0} \end{equation}

with

(3.40)\begin{equation} L_H = \frac{V_A}{U_0}\sqrt{\frac{m}{\mu_0ne^2}}. \end{equation}

In lattice units,

(3.41)\begin{equation} \epsilon_H^\mathrm{lbm} = \frac{L_H}{L_0} \frac{U^\mathrm{lbm}_0 L_0^\mathrm{lbm}}{B_0^\mathrm{lbm}} = \frac{L_H}{\Delta x} \frac{U^\mathrm{lbm}_0}{B_0^\mathrm{lbm}} = \frac{U^\mathrm{lbm}_0 L_H^\mathrm{lbm}}{B_0^\mathrm{lbm}}. \end{equation}

If $U_0=V_A$, $L_H$ is equal to the ion inertial length $d_i$ (or ion skin depth). In that case, $\epsilon _H^\mathrm {lbm} = L_H^\mathrm {lbm}$ and this corresponds to the number of lattice points per ion inertial length. It is assumed that the dynamics of an MHD plasma develops under the influence of the Hall effect at scales $\ell$ smaller that $L_H$.

3.4. CFL condition for Hall-MHD turbulence

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

(3.42)\begin{equation} \Delta t \le \Delta x/c_\mathrm{max}, \end{equation}

where $c_\mathrm {max}$ refers to the largest speed at which a signal propagates in the solution. In the context of Hall-MHD, $c_\mathrm {max}$ should be identified with the largest phase speed of 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)=k V_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 nm_i}$, $m_i$ being the mass of ions and the Hall parameter being $\alpha _H=1/ne$. Therefore, one obtains that the time step decreases quadratically with the grid spacing as

(3.43)\begin{equation} \Delta t \le \frac{\mu_0 (\Delta x)^2}{{\rm \pi} ~\alpha_H B}, \end{equation}

assuming the largest attainable wavenumber to be $k_\mathrm {max}= {\rm \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}$ as

(3.44)\begin{equation} \Delta t \le \frac{ (\Delta x)^2}{{\rm \pi}(\alpha_H/\sqrt{\mu_0}) B_0}, \end{equation}

which finally yields in lattice units to

(3.45)\begin{equation} 1 \le \frac{ 1}{{\rm \pi}~ \epsilon_H^\mathrm{lbm} L_0^\mathrm{lbm} U_0^\mathrm{lbm}} = \frac{1}{\rm \pi} \frac{B_0^\mathrm{lbm}}{\epsilon_H \left(N {U_0^\mathrm{lbm}}\right)^2}, \end{equation}

where (3.34) is used to retrieve $\alpha _H/\sqrt {\mu _0}$, and $N=L_0^\mathrm {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^\mathrm {lbm}}$, which is in turn written as

(3.46)\begin{equation} {Ma}\le \left(\frac{\sqrt{3}}{\rm \pi} \right)\frac{1}{\epsilon_H N^2}. \end{equation}

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 by 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 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 implementation of the LB scheme 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 central processing unit (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 centre (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 (Gómez et al. 2010; Meyrand & Galtier 2012; Ferrand et al. 2022; Yadav, Miura & Pandit 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

(4.1a,b)\begin{equation} \boldsymbol{u}(\boldsymbol{x},t) = \boldsymbol{u}^\prime(\boldsymbol{x},t) \quad \mathrm{and} \quad \boldsymbol{b}(\boldsymbol{x},t) = \hat{\boldsymbol{e}}_z + \boldsymbol{b}^\prime(\boldsymbol{x},t), \end{equation}

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

(4.2)\begin{align} \boldsymbol{u}^\prime(\boldsymbol{x},t) & = \left[ B (\hat{\boldsymbol{e}}_x + {\rm i}\hat{\boldsymbol{e}}_y) \exp({\rm i} {k}z - {\rm i}\omega_\pm t)\right.\nonumber\\ & \quad + C (\hat{\boldsymbol{e}}_y + {\rm i}\hat{\boldsymbol{e}}_z) \exp({\rm i} {k}x)\nonumber\\ & \left.\quad+ A (\hat{\boldsymbol{e}}_z + {\rm i}\hat{\boldsymbol{e}}_x) \exp({\rm i} {k}y) \right] {\rm e}^{-\nu {k}^2 t} \end{align}

and

(4.3)\begin{equation} \boldsymbol{b}^\prime(\boldsymbol{x},t) = \alpha_\pm \boldsymbol{u}^\prime(\boldsymbol{x},t) \end{equation}

in complex notation. 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 {\boldsymbol {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

(4.4)\begin{equation} \alpha_\pm=-\frac{1}{2} \epsilon_H {k} \pm \sqrt{\frac{\epsilon_H^2{k}^2}{4}+1}. \end{equation}

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.36). 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{\rm \pi}$. The evolution of the velocity field is deterministic from the initial condition

(4.5)\begin{equation} \left. \begin{gathered} {u}^\mathrm{lbm}_x(\boldsymbol{x},0) = {U_0^\mathrm{lbm}} \left(B\sin \left(\frac{4{\rm \pi} z_k}{N}\right)+A\cos \left(\frac{4{\rm \pi} y_j}{N}\right) \right),\\ {u}^\mathrm{lbm}_y(\boldsymbol{x},0) = {U_0^\mathrm{lbm}} \left(B\cos\left(\frac{4{\rm \pi} z_k}{N}\right) + C\sin \left(\frac{4{\rm \pi} x_i}{N}\right) \right),\\ {u}^\mathrm{lbm}_z(\boldsymbol{x},0) = {U_0^\mathrm{lbm}} \left(C\cos \left(\frac{4{\rm \pi} x_i}{N}\right) + A\sin \left(\frac{4{\rm \pi} y_j}{N}\right) \right), \end{gathered} \right\} \end{equation}

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

(4.6)\begin{equation} \left. \begin{gathered} \mathrm{b}^\mathrm{lbm}_x(\boldsymbol{x},0) = \alpha_+ {u}^\mathrm{lbm}_x(\boldsymbol{x},0),\\ \mathrm{b}^\mathrm{lbm}_y(\boldsymbol{x},0) = \alpha_+ {u}^\mathrm{lbm}_y(\boldsymbol{x},0),\\ \mathrm{b}^\mathrm{lbm}_z(\boldsymbol{x},0) = \alpha_+ {u}^\mathrm{lbm}_z(\boldsymbol{x},0) + U_0^\mathrm{lbm}, \end{gathered} \right\} \end{equation}

since $U_0=V_A$. For the 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^\mathrm {lbm}=N$ according to (3.38). This value ensures that the solution is affected by the Hall effect with $L_H = L_0$ from (3.39). A three-dimensional rendering of the initial conditions expressed in (4.5) and (4.6) is displayed in figure 2. With this initialization, the current density $\boldsymbol {J}=\boldsymbol {\nabla }\times \boldsymbol {B}$ is non-zero at $t=0$. The parameters used in the different simulations are reported in table 1. The Mach number is always small enough for the plasma to approach the incompressible limit and 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 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 ($\boldsymbol {b} = \alpha _+ \boldsymbol {u} + \hat {\boldsymbol {e}}_z$) are presented in the following.

Figure 2. Three-dimensional rendering of the initial condition as indicated in (4.5) and (4.6). The magnitudes of the (a) fluid velocity, (b) magnetic field and (c) density of electric current are here displayed for $N=128$.

Table 1. Parameters of LB simulations. The magnetic Prandtl number is kept unitary. The kinematic viscosity is given in dimensionless units, i.e. normalized by $U_0 L_0$, which means that the Reynolds number ${Re}=1/\nu$.

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 table 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 figure 3. The solution appears as a damped wave propagating in the direction of the ambient magnetic field. The amplitude and the phase of the solution are well captured in the LB simulation. The results obtained for different resolutions and viscosity values at Mach number ${Ma}=0.007$ are shown in figure 4 for ten 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 figure 5. The level of these relative fluctuations is typically of order $10^{-7}$–$10^{-6}$ 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.

Figure 3. Temporal evolution of the velocity magnitude $|\boldsymbol {u}|(\boldsymbol {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\times 10^{-4}$.

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

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

4.3. Dispersion and dissipation errors

The dispersion and dissipation errors of the LB scheme 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 figure 3). Therefore, tagging as $t^\mathrm {max}_{i}$ and $\bar t^\mathrm {max}_i$ 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

(4.7)\begin{equation} \varepsilon_{\phi}= 1-\frac{1}{M}\sum_{i=0}^{M-1}\frac{t^\mathrm{max}_{i+1}-t^\mathrm{max}_i}{\bar t^\mathrm{max}_{i+1}-\bar t^\mathrm{max}_i} \end{equation}

over $M$ oscillating periods. For practical purposes, we have used $M=10$. As expected, it can be observed in figure 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^2$. This confirms a second-order accuracy of the LB scheme. We also found that the dispersion error exhibits a rather constant behaviour when changing the Mach number, and does not seem to be affected by the value of the kinematic viscosity either. Let us remark that some results differ from the global trend, certainly due to the premise of (physical) instabilities at the lowest viscosity. After 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.

(4.8)\begin{equation} \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}. \end{equation}

The dissipation error provides a first measure of the numerical dissipation. Two different scaling behaviours are considered, namely the so-called acoustic and diffusive scaling (Krueger et al. 2016). The acoustic scaling consists in keeping the Mach number fixed while monitoring the convergence rate of the error $\varepsilon _r$ for different Reynolds numbers, as a function of the resolution (see figure 7). However, the diffusive scaling is obtained by keeping the lattice viscosity fixed (see figure 8). The behaviour of the numerical solution is consistent between the two regimes, showing a convergence of the dissipation error with respect to the grid resolution $\propto 1/N^2$, as expected for a second-order scheme.

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 an LB scheme.

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}$.

Figure 8. Diffusive scaling ($\nu$ constant) of the (relative) dissipation error for different values of the Mach number (${Ma}$) at fixed kinematic viscosity ($\nu$).

One of the advantages of dealing with a dissipative solution of the Hall-MHD equations is the possibility to identify an effective viscosity $\tilde \nu$ related to the damping $\propto \exp (-\tilde \nu k^2 t)$ of the numerical solution. By decomposing $\tilde \nu$ into the sum of a physical and a (spurious) numerical viscosity, $\tilde \nu = \nu + \nu _\mathrm {num}$, the ratio between these two contributions reads as

(4.9)\begin{equation} \varepsilon_{\nu} = \frac{\nu_\mathrm{num}}{\nu} = \frac{\tilde \nu-\nu}{\nu}. \end{equation}

The results obtained for the viscosity error $\varepsilon _{\nu }$ are shown in figure 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 of the physical viscosity and the Mach number, whereas it only depends on the lattice resolution.

Figure 9. Scaling of the ratio between the numerical and kinematic viscosities $\varepsilon _{\nu }=\nu ^\mathrm {num}/\nu$ with the grid resolution ($N$) for different Mach numbers (${Ma}$) and kinematic viscosity. The second-order accuracy of the LB algorithm is highlighted by the black lines, i.e. $\varepsilon _{\nu } \sim O(\Delta x^2)$.

Finally, despite Dellar (2002) stating that a D3Q7 lattice was sufficient to reliably account for the dynamics of each component of the magnetic field, 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 figure 1), thus confirming what was reported by 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.20)). 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 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 (Mininni et al. 2002, 2003; Mininni, Pouquet & Montgomery 2006; Gómez et al. 2010; Marino et al. 2013; Pouquet & Marino 2013; Marino et al. 2014; Marino, Pouquet & Rosenberg 2015a; Pouquet et al. 2019). It is a well-established community code available on https://github.com/pmininni/GHOST. The ghost code 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:

(4.10)\begin{gather} \boldsymbol{U(\boldsymbol{x},0})=U_0\left[-2\sin{y};2\sin{x};0\right] \end{gather}

(4.11)\begin{gather}\boldsymbol{B(\boldsymbol{x},0})=B_0\left[-2 \sin{2y}+\sin{z}; 2\sin{x}+\sin{z} ; \sin{x}+\sin{y}\right] \end{gather}

with $U_0=1$ and $B_0=0.8$ in a cubic box of size $2{\rm \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 |\boldsymbol {\nabla } \times \boldsymbol {U}|^2\rangle -\eta \langle |\boldsymbol {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 root mean square velocity and $L=2{\rm \pi} \int k^{-1}E_v(k)\,\textrm {d}k/\int E_v(k)\,\textrm {d}k$ 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 figure 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 $\boldsymbol {J}=\boldsymbol {\nabla } \times \boldsymbol {B}$. Let us recall at this stage that $\boldsymbol {J}$ is directly obtained from the magnetic densities in the LB simulation, and is not inferred by differentiating the magnetic field.

Figure 10. (a) Time evolution of the mean magnetic dissipation $\propto \langle |\boldsymbol {J}|^2 \rangle$ in freely evolving MHD turbulence for an LB simulation ($N=512$) and a pseudo-spectral simulation (N = $512$) performed with the ghost solver. (b) Time evolution of the mean kinetic ($E_v$) and magnetic $(E_B$) energies.

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 figure 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 with the standard BGK collision operator) but at the cost of an enhanced numerical dissipation (Coreixas et al. 2019). However, when increasing the spatial resolution to $N=768$, the numerical dissipation is reduced and the spectrum of the flame run gets very close to that of the pseudo-spectral solution. This observation is qualitatively consistent with the statement made by Shen et al. (2018) that LB needs approximately twice the spatial resolution of a pseudo-spectral simulation to achieve similar accuracy in turbulent flows.

Figure 11. (a) Kinetic and (b) magnetic energy spectra of MHD turbulence at the peak of magnetic dissipation. The spectra are normalized 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.

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 with the pseudo-spectral simulation). It should also be noted that, while the maximum wavenumber is $k_\mathrm {max}=N/3$ (due to the 2/3 rule for de-aliasing) in pseudo-spectral simulations, the range of resolved scales reaches the Nyquist cutoff $k_\mathrm {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

(4.12)\begin{equation} \mathcal{S}_\mathrm{MHD}(k) = \sum_{|\boldsymbol{k'}|< k} \mathrm{Re}{\left[\mathcal{F}({\boldsymbol{U}})^*\boldsymbol{\cdot} \left( \mathcal{F}({\boldsymbol{U}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{U}}) - \mathcal{F}({\boldsymbol{J}\times \boldsymbol{B}}) \right) - \mathcal{F}({\boldsymbol{B}})^*\boldsymbol{\cdot} \mathcal{F}(\boldsymbol{\nabla} \times ({\boldsymbol{U}\times \boldsymbol{B}}))\right]}, \end{equation}

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

(4.13)\begin{equation} \mathcal{D}(k) = \sum_{|\boldsymbol{k'}|< k} \nu k'^{2}|\mathcal{F}({\boldsymbol{U}})|^2+\eta k'^{2}|\mathcal{F}({\boldsymbol{B}})|^2, \end{equation}

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

(4.14)\begin{equation} \partial_t \sum_{|\boldsymbol{k}^\prime|< k} E(\boldsymbol{k}^\prime) =-\mathcal{S}_\mathrm{MHD}(k) -\mathcal{D}(k). \end{equation}

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 figure 12. A satisfactory agreement is observed in particular for the nonlinear 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 12. Various third-order energy fluxes across wavenumber $k$ as reported in (4.13). Comparison is made between (a) LB and (b) de-aliased pseudo-spectral ghost simulations. Let us notice that $\partial _t \sum _{|\boldsymbol {k}^\prime |< k} E(\boldsymbol {k}^\prime ) = -\mathcal {S}_\mathrm {MHD}(k) -\mathcal {D}(k) <0$ as expected for freely decaying MHD turbulence.

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

The flame solver 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{\rm \pi}$ with resolution $512^3$ and $768^3$, initialized with the OT vortex as described in the previous section, for different values of the Hall parameter (see table 2). The Mach number was adjusted to the Hall parameter to accommodate the CFL condition based on the time scale of the whistler waves (see (3.46)). The Reynolds number is estimated here at the peak of magnetic dissipation (indicated by the vertical dashed lines in figure 13), at which the plasma is assumed to have reached a fully developed turbulent state. For a $512^3$ lattice dimension, only three GPUs were used in parallel, resulting in a computational speed of approximately 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 for the three runs indicated in table 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 ($\boldsymbol {u}$, $\boldsymbol {B}$ and $\boldsymbol {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 $\boldsymbol {u}$ and $\boldsymbol {B}$ is shown in figure 14 for the simulation at the highest resolution (run IV in table 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 figure 15 for the electric current density $\boldsymbol {J}$ and the vorticity $\boldsymbol {\omega }=\boldsymbol {\nabla }\times \boldsymbol {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 behaviour captured in figure 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 figure 13) are plotted in figure 16 for each run at resolution $512^3$. As expected, increasing the value of the Hall parameter $\epsilon _H$ (indicated by the vertical dash-dotted line in figure 16) produces a shift of the Hall length-scale $L_H$ towards larger scales, and hence a shrinkage of the Kolmorogov'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 behaviour 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 by Kiyani, Osman & Chapman (2015).

Figure 13. Evolution of the magnetic dissipation over time for the three simulations performed with the OT initial condition (see table 2). The shaded areas around the peak of the current density (black dashed line) correspond to the range over which the energy spectra in figure 16 have been averaged.

Figure 14. Results of a 3-D Hall-MHD simulation at $N=768^3$ taken at the peak of dissipation. (a) A three-dimensional visualization of the normalized velocity field magnitude $v^\prime =(v-\bar {v})/\sigma _v$ with $v=|\boldsymbol {v}|$ and (b) the relative power density spectrum. The same in panels (c,d) for the normalized magnetic field $B^\prime$. The slopes $k^{-5/3}$ and $k^{-2.73}$ are given as a reference.

Figure 15. (a) Three-dimensional rendering of the normalized current density $J^\prime =(J-\bar {J})/\sigma _J$, with $J=|\boldsymbol {J}|$, and (b) vorticity $\omega ^\prime$ taken at the same time as in figure 14.

Figure 16. Kinetic (blue) and magnetic (orange) spectra corresponding to three different Hall-MHD turbulence regimes (see table 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{\rm \pi} /(\epsilon _H L_0)$ with $L_0=2{\rm \pi} \int E_v(k)k^{-1}\,\textrm {d}k/\int E_v(k)\,\textrm {d}k$.

Table 2. Parameters of Hall-MHD turbulence runs. Here, ${Re}$, ${Ma}$ and ${Pr}_m$ denote respectively the Reynolds number (at the peak of magnetic dissipation), the initial Mach number and the magnetic Prandtl number. The (dimensionless) Hall parameter is $\epsilon _H$. The number of grid points per dimension is $N$. The total duration of the run is $t_\mathrm {tot}/\tau _{0}$ and the time at which the peak of current density occurs is $t_\mathrm {peak}/\tau _{0}$ in units of the reference time scale $L_0/U_0$. The Mach number satisfies the CFL condition (3.46) imposed by whistler waves.

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 with 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 with 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 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 $\boldsymbol {U}$ with respect to the Hall current $\alpha _H \boldsymbol {J}$ in (3.27).

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, Cuesta & Matthaeus 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). Figure 17 shows the trace spectra of the magnetic field fluctuations measured by the PSP/FIELDS (Bale et al. 2016) and SO/MAG (Horbury et al. 2020) magnetometers on board 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 table 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, Lacombe & Mangeney 2008; Alexandrova et al. 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 figure 17, marking the transition from the MHD turbulent inertial range (where energy is adiabatically transferred to increasingly smaller scales), which 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 a matter for 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 that 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, Trenchi & Telloni 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 for 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 (figure 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 $3$-D 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 wind turbulence has been presented by Sorriso-Valvo 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.

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 coloured lines.

Table 3. Main solar wind parameters computed at the time where the solar wind samples used to produce the power spectra in figure 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}$.

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 used 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 time stepping to ensure stability with longer time steps than what is imposed by the CFL condition for explicit time stepping. In our (explicit) scheme, the time step was established by default to meet this latter condition according to (3.46). 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 undergo 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 behaviour of the magnetic energy spectrum at that 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, Dmitruk & Andrés 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.