Nonlinear fourth order Taylor expansion of lattice Boltzmann schemes

TL;DR

This paper develops a nonlinear fourth-order Taylor expansion for lattice Boltzmann schemes, providing explicit formulas and validation with the D2Q9 scheme to improve understanding of scheme accuracy.

Contribution

It introduces a nonlinear fourth-order expansion framework for lattice Boltzmann schemes, coupling asymptotic corrections with explicit formulas for conserved and nonconserved moments.

Findings

Validated algebraic expressions with previous results

Derived explicit formulas for nonlinear fourth-order accuracy

Applied framework to isothermal D2Q9 lattice Boltzmann scheme

Abstract

We propose a formal expansion of multiple relaxation times lattice Boltzmann schemes in terms of a single infinitesimal numerical variable. The result is a system of partial differential equations for the conserved moments of the lattice Boltzmann scheme. The expansion is presented in the nonlinear case up to fourth order accuracy. The asymptotic corrections of the nonconserved moments are developed in terms of equilibrium values and partial differentials of the conserved moments. Both expansions are coupled and conduct to explicit compact formulas. The new algebraic expressions are validated with previous results obtained with this approach. The example of isothermal D2Q9 lattice Boltzmann scheme illustrates the theoretical framework.

Equations

$$f_j=f_j(x,t),x\in {\mathbb{R}}^d,t\ge 0,0\le j<q.\text{vspace}\ast -6pt$$

$$\frac{\partial f_j}{\partial t}+\sum _{1\le \alpha \le d}{v}_j^{\alpha }\frac{\partial f_j}{\partial x_{\alpha }}=\frac{1}{\epsilon }{Q}_j(f),0\le j<q.\text{vspace}\ast -7pt$$

$$m_k\equiv \sum _{0\le j<q}{M}_{kj}{f}_j.\text{vspace}\ast -7pt$$

$$\sum _{0\le j<q}{M}_{kj}{Q}_j(f)=0,\forall f\in {\mathbb{R}}^q,0\le k<N.\text{vspace}\ast -7pt$$

$$m\equiv \left(\begin{array}{c}W\\ Y\end{array}\right).\text{vspace}\ast -7pt$$

$$\frac{\partial W_k}{\partial t}+\sum _{1\le \alpha \le d}\sum _{0\le j<q}{M}_{kj}{v}_j^{\alpha }\frac{\partial f_j}{\partial x_{\alpha }}=0,0\le k<N.\text{vspace}\ast -6pt$$

$$Q(f^{\mathrm{eq}})=0\text{vspace}\ast -6pt$$

$$f^{\mathrm{eq}}=M^{-1}\left(\begin{array}{c}W\\ \Phi (W)\end{array}\right).\text{vspace}\ast -7pt$$

$$\displaystyle Z={\rm diag}\,\Big{(}{{1}\over{\tau_{1}}},\,\dots,\,{{1}\over{\tau_{q-N}}}\Big{)}\vspace*{-6pt}$$

$$M\mathrm{d}Q(f^{\mathrm{eq}})M^{-1}=-\left(\begin{array}{cc}0& 0\\ 0& Z\end{array}\right).\text{vspace}\ast -7pt$$

$$\displaystyle\Lambda=M\,\,{\rm diag}\,\Big{(}\sum_{1\leq\alpha\leq d}v^{\alpha}\,\partial_{\alpha}\Big{)}\,\,M^{-1}\,.\vspace*{-7pt}$$

$$\Lambda \equiv \left(\begin{array}{cc}A& B\\ C& D\end{array}\right).\text{vspace}\ast -7pt$$

$$Q(f)\simeq Q(f^{\mathrm{eq}})+\mathrm{d}Q(f^{\mathrm{eq}}).(f-f^{\mathrm{eq}}).\text{vspace}\ast -6pt$$

$$\frac{\partial f_j}{\partial t}+\sum _{1\le \alpha \le d}{v}_j^{\alpha }\frac{\partial f_j}{\partial x_{\alpha }}=\frac{1}{\epsilon }\mathrm{d}{Q}_j(f^{\mathrm{eq}}).(f-f^{\mathrm{eq}}),0\le j<q.\text{vspace}\ast -7pt$$

$$\frac{\partial W}{\partial t}+AW+BY=0,\frac{\partial Y}{\partial t}+CW+DY=-\frac{1}{\epsilon }Z(Y-\Phi (W)).\text{vspace}\ast -7pt$$

$$\frac{\partial W}{\partial t}+{\Gamma }_1(W)+\epsilon {\Gamma }_2(W)+{\epsilon }^2{\Gamma }_3(W)+{\epsilon }^3{\Gamma }_4(W)=\mathrm{O}({\epsilon }^4).\text{vspace}\ast -7pt$$

$$Y=\Phi (W)+\epsilon Z^{-1}{\Psi }_1(W)+{\epsilon }^2{Z}^{-1}{\Psi }_2(W)+{\epsilon }^3{Z}^{-1}{\Psi }_2(W)+\mathrm{O}({\epsilon }^4)\text{vspace}\ast -7pt$$

$$\displaystyle\left\{\begin{array}[]{rl}\Gamma_{1}(W)&\!\!\!=A\,W+B\,\Phi(W)\\ \Psi_{1}(W)&\!\!\!={\rm d}\Phi(W).\Gamma_{1}(W)-\big{(}C\,W+D\,\Phi(W)\big{)}\\ \Gamma_{2}(W)&\!\!\!=B\,Z^{-1}\,\Psi_{1}(W)\\ \Psi_{2}(W)&\!\!\!=Z^{-1}\,{\rm d}\Psi_{1}(W).\Gamma_{1}(W)+{\rm d}\Phi(W).\Gamma_{2}(W)-D\,Z^{-1}\,\Psi_{1}(W)\\ \Gamma_{3}(W)&\!\!\!=B\,Z^{-1}\,\Psi_{2}(W)\\ \Psi_{3}(W)&\!\!\!=Z^{-1}\,{\rm d}\Psi_{1}(W).\Gamma_{2}(W)+{\rm d}\Phi(W).\Gamma_{3}(W)-D\,Z^{-1}\,\Psi_{2}(W)\\ &+Z^{-1}\,{\rm d}\Psi_{2}(W).\Gamma_{1}(W)\\ \Gamma_{4}(W)&\!\!\!=B\,Z^{-1}\,\Psi_{3}(W)\,.\end{array}\right.\vspace*{-7pt}$$

$$\displaystyle{{\partial W}\over{\partial t}}+A\,W+B\,\big{(}\Phi(W)+\varepsilon\,Z^{-1}\,\Psi_{1}(W)+\varepsilon^{2}\,Z^{-1}\,\Psi_{2}(W)+\varepsilon^{3}\,Z^{-1}\,\Psi_{2}(W)\big{)}={\rm O}(\varepsilon^{4})\,\vspace*{-6pt}$$

$$\frac{\partial W}{\partial t}+AW+B\Phi (W)+\epsilon B{Z}^{-1}{\Psi }_1(W)+{\epsilon }^{2}B{Z}^{-1}{\Psi }_2(W)+{\epsilon }^{3}B{Z}^{-1}{\Psi }_3(W)=\mathrm{O}({\epsilon }^4).\text{vspace}\ast -6pt$$

$$\displaystyle\left\{\begin{array}[]{rl}\Gamma_{1}(W)&\!\!\!=A\,W+B\,\Phi(W)\\ \Gamma_{2}(W)&\!\!\!=B\,Z^{-1}\,\Psi_{1}(W)\\ \Gamma_{3}(W)&\!\!\!=B\,Z^{-1}\,\Psi_{2}(W)\\ \Gamma_{4}(W)&\!\!\!=B\,Z^{-1}\,\Psi_{3}(W)\end{array}\right.\vspace*{-6pt}$$

$$Y=\Phi -\epsilon Z^{-1}({\partial }_{t}Y+CW+DY).\text{vspace}\ast -6pt$$

$$\displaystyle\left\{\begin{array}[]{l}\Psi_{1}(W)={\rm d}\Phi(W)\,.\,\Gamma_{1}-(C\,W+D\,\Phi)\\ \Psi_{2}(W)=Z^{-1}\,{\rm d}\Psi_{1}(W)\,.\,\Gamma_{1}+{\rm d}\Phi(W)\,.\,\Gamma_{2}(W)-D\,Z^{-1}\,\Psi_{1}(W)\\ \Psi_{3}(W)=Z^{-1}\,{\rm d}\Psi_{1}(W)\,.\,\Gamma_{2}+{\rm d}\Phi(W)\,.\,\Gamma_{3}-D\,Z^{-1}\,\Psi_{2}+Z^{-1}\,{\rm d}\Psi_{2}(W)\,.\,\Gamma_{1}\end{array}\right.\vspace*{-6pt}$$

$$f_j=f_j(x,t),x\in \mathcal{L},t=n\Delta t,n\in \mathbb{N},v_j\in \mathcal{V}.\text{vspace}\ast -6pt$$

$$\frac{\partial f}{\partial t}=\frac{1}{\epsilon }\mathrm{d}Q(f^{\mathrm{eq}}).(f-f^{\mathrm{eq}})\text{vspace}\ast -6pt$$

$$\frac{\partial f_j}{\partial t}+\sum _{1\le \alpha \le d}{v}_j^{\alpha }\frac{\partial f_j}{\partial x_{\alpha }}=0,0\le j<q.\text{vspace}\ast -6pt$$

$$\frac{\partial W}{\partial t}=0,\frac{\partial Y}{\partial t}=-\frac{1}{\epsilon }Z(Y-\Phi (W)).\text{vspace}\ast -6pt$$

$$m^{\ast }\equiv \left(\begin{array}{c}{W}^{\ast }\\ Y^{\ast }\end{array}\right).\text{vspace}\ast -7pt$$

$$W^{\ast }=W.\text{vspace}\ast -7pt$$

$$\frac{Y^{\ast }-Y}{\Delta t}=-\frac{1}{\epsilon }Z(Y-\Phi (W)),\text{vspace}\ast -6pt$$

Full text

Nonlinear fourth order Taylor expansion

of lattice Boltzmann schemes

François Duboisab

a* Laboratoire de Mathématiques d’Orsay, Faculté des Sciences d’Orsay,*

Université Paris-Saclay, France.

b* Conservatoire National des Arts et Métiers, LMSSC laboratory, Paris, France.*

19 June 2020 ***  This contribution to appear in Asymptotic Analysis has been first presented at the Institut Henri Poincaré (Paris, France) the 10 January 2018. Edition 14 January 2021.

Keywords: partial differential equations, asymptotic analysis

AMS classification: 65Q05, 76N15, 82C20.

PACS numbers: 02.70.Ns, 05.20.Dd, 47.10.+g

**Abstract **

We propose a formal expansion of multiple relaxation times lattice Boltzmann schemes in terms of a single infinitesimal numerical variable. The result is a system of partial differential equations for the conserved moments of the lattice Boltzmann scheme. The expansion is presented in the nonlinear case up to fourth order accuracy. The asymptotic corrections of the nonconserved moments are developed in terms of equilibrium values and partial differentials of the conserved moments. Both expansions are coupled and conduct to explicit compact formulas. The new algebraic expressions are validated with previous results obtained with this framework. The example of isothermal D2Q9 lattice Boltzmann scheme illustrates the theoretical framework.

**1)  Introduction **

The lattice Boltzmann schemes in their modern form have been developed with the contributions of d’Humières, Lallemand, Succi [39, 51, 38, 44] and many others. An underlying lattice Boltzmann equation is discretised on a cartesian grid with a finite set of velocities chosen in such a way that during one time step, an exact transport is done between two vertices of the mesh. A lattice Boltzmann scheme is composed with two steps: a nonlinear local relaxation step, followed by a linear advection scheme coupling a given vertex with a given family of neighbours. The relaxation step follows in general the “BGK” approximation introduced by Bhatnagar, Gross and Krook [5]. This numerical method allows the simulation of an important number of physical phenomena as isothermal flows, compressible flows with heat transfer, non-ideal fluids, multiphase and multi-component flows, microscale gas flows, soft-matter flows,… up to quantum mechanics. For the actual status of the method and the various applications, we refer e.g. to the books of Guo and Shu [33], Krüger et al. [36] and to the prospective article of Succi [57]. The lattice Boltzmann method is founded on a mesoscopic Boltzmann model but is not able in general to solve the Boltzmann equation or associated kinetic models. It is admitted that essentially macroscopic models are approximated with the lattice Boltzmann schemes.

The link between mesoscopic and macroscopic models is never straightforward and an analysis is necessary. The usual method of formal analysis is founded on the Chapman-Enskog [9] expansion in the spirit of the book of Huang [40]. The classical approach for the Boltzmann equation and other mesoscopic models posed in a continuum space-time. The Chapman-Enskog method has been adapted to lattice Boltzmann schemes with discrete space and time equations by Chen and Doolen [10] and by Qian and Zhou [52]. It allows the derivation the equivalent partial differential equations founded on a multiscale time analysis. This classical Chapman Enskog expansion is popular in the lattive Boltzmann community, as reported in the contributions of Lallemand and Luo [44], Philippi and Hegele [49] and Shan et al. [54] among others. A main drawback of this approach is the fact that multiscale expansions are used without a clear mathematical signification of the various variables and associated functions.

Independently of this framework, we have proposed in [15, 16] the Taylor expansion method for isothermal fluid problems to obtain formally equivalent partial differential equations at second order accuracy. This work was inspired by the method of equivalent equations for finite difference schemes, in a spirit proposed by Lerat and Peyret [46] and Warming and Hyett [58]. In this approach, the infinitesimal variable is simply the time step (proportional to the space step with the acoustic scaling) and there is only one infinitesimal scale for the analysis. This approach has been experimentaly validated with Pierre Lallemand in various contributions [22, 23, 25, 43]. It was extended to third order accuracy for thermal and fluid problems in [17], automated at an arbitrary order in the linear case in [4], extended with an external drift in [26] and to relative velocities in [19, 20]. Our point of view is a priori not to consider the third and fourth order terms as approximations of Burnett of super-Burnett equations (see e.g. [2, 56]) but simply as errors of the lattice Boltzmann scheme when solving partial differential equations of second order type. Nevertheless, progress in the approximation of fourth order partial differential equations have been obtained in [48]. Here, we have no particular fluid hypothesis. It can be applied to scalar or to vectorial [18, 30] lattice Bolzmann schemes. This Taylor expansion method has been also revisited with the Maxwell iterations of Yong et al. [59].

In this contribution, we first recall in Section 2 classic results related to the Boltzmann equation with discrete velocities. The link with multiple relaxation times lattice Boltzmann schemes as initially proposed by d’Humières [38] is presented in Section 3. Then we present the Taylor expansion method in a very general way and the result is explicited with compact formulas. Our main result (Section 4) concerns the derivation of asymptotic partial differential equations in the nonlinear case, up to fourth order accuracy. The proof for first and second orders is presented in Section 5, with a validation for the isothermal fluid case in the specific case of the D2Q9 scheme. In Section 6, the expansion up to third and fourth orders is detailed, with a confrontation with our previous work [17]. The proofs of these results need specific algebraic developments, presented in this contribution. The linear case with constant coefficients is described in Section 7.

**2)  Boltzmann equation with discrete velocities **

In the space $\,\mathbb{R}^{d}\,$ of dimension $\,d$, we consider a finite set of $\,q\,$ discrete velocities $\,v_{j}\in{\cal V}\,$ with components $\,v_{j}^{\alpha}\,$ for $\,1\leq\alpha\leq d$. The unknowns of the Boltzmann equation are the particle densities $\,f_{j}$. They are functions of space $\,x$, time $\,t\,$ and discrete velocities $\,v_{j}\,$:

[TABLE]

The vector $\,f(x,\,t)\in\mathbb{R}^{q}\,$ is constructed with the numbers $\,f_{j}(x,\,t)\,$ for $\,0\leq j<q$. We introduce a collision vector $\,\mathbb{R}^{q}\ni f\longmapsto Q(f)\in\mathbb{R}^{q}\,$ of components $\,Q_{j}(f)\,$ for $\,0\leq j<q\,$ with given regular functions $\,Q_{j}$. We introduce also a small parameter $\,\varepsilon>0\,$ that can be interpreted as a Knudsen number in the context of gas dynamics (see e.g. [9]).

The Boltzmann model with discrete velocities takes the form

[TABLE]

It has been proposed by Carleman [8], Gross [32], Broadwell [7] and intensively developed by Gatignol [28] and her co-workers.

$\bullet$  **Moments **

Our framework concerns multi relaxation times: we introduce a constant invertible matrix $\,M\,$ called “d’Humières matrix” [38] in this contribution. This matrix defines the vector of moments $\,m\in\mathbb{R}^{q}\,$ by a simple product:

[TABLE]

We introduce the number of conservation laws $\,N\,$ ($1\leq N<q$) such that the $\,N\,$ first moments of the collision kernel are equal to zero:

[TABLE]

Then it is natural to divide the vector of moments into two families:

[TABLE]

The conserved moments or macroscopic variables $\,W\,$ constitute a linear space of dimension $\,N\,$ and if the Boltzmann model with discrete velocities (1) is satisfied, we have $\,N\,$ conservation laws due to the cancellation (3):

[TABLE]

Observe that the nonconserved moments or microscopic variables $\,Y\,$ generate a linear space of dimension $\,q-N$.

$\bullet$  **Equilibrium states **

We suppose that the equilibrim states $\,f^{\rm eq}\,$ defined by the conditions

[TABLE]

are characterized with the help of a regular nonlinear vector field $\,\Phi\,:\mathbb{R}^{N}\longrightarrow\mathbb{R}^{q-N}\,$ such that

[TABLE]

In other words, the vector field $\,W\longmapsto Y^{\rm eq}\equiv\Phi(W)\,$ defines the set of equilibrium states.

We suppose moreover that the jacobian matrix $\,{\rm d}Q(f^{\rm eq})\,$ at equilibrium is diagonalizable with real eigenvalues and real eigenvectors. More precisely, taking into account the hypothesis (3), we suppose that there exists a diagonal matrix

[TABLE]

of order $\,q-N\,$ with strictly positive coefficients $\,\tau_{\ell}\,$ such that

[TABLE]

In other words, the jacobian operator $\,{\rm d}Q(f^{\rm eq})\,$ admits the matrix $\,M^{-1}\,$ as a matrix of right eigenvectors and the associated eigenvalues are all non positive.

$\bullet$  Momentum-velocity operator

For a linear space of dimension $\,d$, we introduce the momentum-velocity operator matrix $\,\Lambda\,$ defined by the relation

[TABLE]

It is a $\,q\times q\,$ operator matrix composed by first-order space differential operators. It is obtained by conjugation of the first order advection operator $\,v.\nabla\,$ by the d’Humières matrix $\,M$. The operator matrix $\,\Lambda\,$ is nothing else than the advection operator seen in the basis of moments. We introduce a block decomposition of the momentum-velocity operator matrix associated to the decomposition (4) of the moments. We define a $\,N\times N\,$ operator matrix $\,A$, a $\,N\times(q-N)\,$ operator matrix $\,B$, a $\,(q-N)\times N\,$ operator matrix $\,C\,$ and a $\,(q-N)\times(q-N)\,$ operator matrix $\,D\,$ according to

[TABLE]

Remember that in the following, the matrices $\,A$, $\,B$, $\,C\,$ and $\,D\,$ are matrices composed with first order space operators.

$\bullet$  **Boltmann-BGK system with discrete velocities **

We approximate now the previous discrete model with a BGK [5] type hypothesis: the state $\,f\,$ is close to equilibrium and we approach the collision kernel $\,Q\,$ by its first order expansion around the equilibrium state:

[TABLE]

Due to the condition $\,Q(f^{\rm eq})=0$, we obtain with this approximation that we call “BGK” in this contribution, even it is not exactly the hypothesis done in the original article [5], a new system of partial differential equations:

[TABLE]

Observe that $\,f^{\rm eq}\,$ is a function of the particle distribution $\,f\,$ through the relations (2), (4) and (5). Taking into account the hypothesis (6), the definition (7) and the notation (8), we can write the Boltmann-BGK system with discrete velocities (9) under the form of a system of two coupled partial differential equations:

[TABLE]

Proposition 1.  Fourth order Chapman-Enskog expansion of the Boltmann-BGK system with discrete velocities

As $\,\varepsilon\,$ tends to zero, the conserved moments $\,W\,$ of the Boltmann-BGK system with discrete velocities (10) satisfy formally the following asymptotic system of partial differential equations:

[TABLE]

The nonlinear operators are related to the expansion of the microscopic moments:

[TABLE]

and we have the interlaced relations

[TABLE]

In applications to fluid dynamics, the first order term $\,\Gamma_{1}(W)\,$ is associated to the Euler equations of gas dynamics and the second order term $\,\Gamma_{2}(W)\,$ to the viscous terms of the Navier-Stokes equations. The operators $\,\Gamma_{3}(W)\,$ and $\,\Gamma_{4}(W)\,$ represent the Burnett and super-Burnett terms respectively. With this Proposition 1, we essentially reformulate the classic Chapman Enskog expansion [9] in the discrete velocity case, as done e.g. in the book of Gatignol [28].

$\bullet$  Proof of Proposition 1.

First we inject the representation (12) inside the first equation of (10). It comes

[TABLE]

then after rearranging the terms,

[TABLE]

The confrontation of this partial differential equation with the Ansatz (11) shows that

[TABLE]

as proposed in (13).

Secondly, we write the second equation of (10) under the form

[TABLE]

We observe that

$-\partial_{t}Y={\rm d}Y\,.\,(-\partial_{t}W)$

\,\,\,\,={\rm d}\big{(}\Phi(W)+\varepsilon\,Z^{-1}\,\Psi_{1}(W)+\varepsilon^{2}\,Z^{-1}\,\Psi_{2}(W)\big{)}\,.\,\big{(}\Gamma_{1}+\varepsilon\,\Gamma_{2}(W)+\varepsilon^{2}\,\Gamma_{3}(W)\big{)}+{\rm O}(\varepsilon^{3})

\,\,\,\,={\rm d}\Phi(W)\,.\,\Gamma_{1}+\varepsilon\,\big{(}Z^{-1}\,{\rm d}\Psi_{1}(W)\,.\,\Gamma_{1}+{\rm d}\Phi(W)\,.\,\Gamma_{2}\big{)}

+\varepsilon^{2}\,\big{(}Z^{-1}\,{\rm d}\Psi_{2}(W)\,.\,\Gamma_{1}+Z^{-1}\,{\rm d}\Psi_{1}(W)\,.\,\Gamma_{2}+{\rm d}\Phi(W)\,.\,\Gamma_{3}\big{)}+{\rm O}(\varepsilon^{3})\,.

We report this expression inside the value of $\,Y\,$:

$Y=\Phi+\varepsilon\,Z^{-1}\,(-\partial_{t}Y-C\,W-D\,Y)$

\,=\Phi+\varepsilon\,Z^{-1}\,\big{[}{\rm d}\Phi(W)\,.\,\Gamma_{1}+\varepsilon\,(Z^{-1}\,{\rm d}\Psi_{1}(W)\,.\,\Gamma_{1}+{\rm d}\Phi(W)\,.\,\Gamma_{2})

$+\varepsilon^{2}\,(Z^{-1}\,{\rm d}\Psi_{2}(W)\,.\,\Gamma_{1}+Z^{-1}\,{\rm d}\Psi_{1}(W)\,.\,\Gamma_{2}+{\rm d}\Phi(W)\,.\,\Gamma_{3})$

-\,C\,W-D\,\big{(}\Phi(W)+\varepsilon\,Z^{-1}\,\Psi_{1}(W)+\varepsilon^{2}\,Z^{-1}\,\Psi_{2}(W)\big{)}\big{]}+{\rm O}(\varepsilon^{4})

$\,=\Phi+\varepsilon\,Z^{-1}\,({\rm d}\Phi(W)\,.\,\Gamma_{1}-C\,W-D\,\Phi)+\varepsilon^{2}\,Z^{-1}\,(Z^{-1}\,{\rm d}\Psi_{1}(W)\,.\,\Gamma_{1}+{\rm d}\Phi(W)\,.\,\Gamma_{2}-D\,Z^{-1}\,\Psi_{1})$

$+\varepsilon^{3}\,Z^{-1}\,(Z^{-1}\,{\rm d}\Psi_{2}(W)\,.\,\Gamma_{1}+Z^{-1}\,{\rm d}\Psi_{1}(W)\,.\,\Gamma_{2}+{\rm d}\Phi(W)\,.\,\Gamma_{3}-D\,Z^{-1}\,\Psi_{2})+{\rm O}(\varepsilon^{4})\,.$

Comparatively to the relation (12), we obtain

[TABLE]

and the last relations of (13) are explicited. The proof is completed. $\square$

**3)  Multiple Relaxation Times Lattice Boltzmann schemes **

We introduce now a regular cartesian lattice $\,{\cal L}\,$ composed by vertices $\,x\,$ separated by distances that are simple expressions of the space step $\,\Delta x$. A discrete time $\,t\,$ is supposed to be an integer multiple of a time step $\,\Delta t>0$. The unknowns of a lattice Boltzmann scheme with $\,q\,$ discrete velocities are the particle densities $\,f_{j}$. They are now functions of discrete space $\,x$, discrete time $\,t\,$ and discrete velocities $\,v_{j}\,$ for $\,0\leq j<q$:

[TABLE]

The discrete velocity set $\,{\cal V}\,$ does not depend on space or time.

Multiple relaxation times lattice Boltzmann scheme considered in this contribution consists in a simple Lie-Trotter splitting applied to the Boltzmann-BGK model with discrete velocities (9). First a relaxation step $\,f(t)\longrightarrow R(\Delta t)\,f(t)\equiv f^{*}\,$ for the approximation of the system of ordinary differential equations

[TABLE]

and secondly a transport step $\,f^{*}\longrightarrow T(\Delta t)\,f^{*}\equiv f(t+\Delta t)\,$ for the free advection

[TABLE]

$\bullet$  **Relaxation step **

Due to the formulation (10), the relaxation scheme is simply described in the moment representation and we have

[TABLE]

After relaxation, the vector of moments $\,m\,$ is transformed into a new vector $\,m^{*}\,$:

[TABLE]

The conserved moments are not modified during this relaxation step:

[TABLE]

As remarked in [16], we use an explicit first order Euler scheme for the microscopic moments:

[TABLE]

*id est *

[TABLE]

with

[TABLE]

Usually, as pointed by Lallemand and Luo [44], the relaxation matrix $\,S\,$ is a diagonal matrix, in coherence with the hypothesis concerning the diagonalization of the operator $\,{\rm d}Q(f^{\rm eq})\,$ presented in the previous section:

[TABLE]

Moreover, the matrix $\,S\,$ is a fixed invertible square matrix of dimension $\,(q-N)\times(q-N)$. In other words, the time step $\,\Delta t\,$ is if the order of the Knudsen number $\,\varepsilon\,$ in our asymptotic study. Finally the particle distribution $\,f^{*}\,$ after relaxation is given according to

[TABLE]

$\bullet$  **Necessary conditions for stability **

We deduce from the relation (16)

[TABLE]

and a natural stability condition takes the form $\,||\,{\rm I}-S\,||\leq 1$, id est

[TABLE]

The case of over-relaxation $\,1<s_{j}<2\,$ is widely used when implementing lattice Boltzmann schemes to take into account very law viscosities in fluid mechanics. We will emphasize this point in the next section.

$\bullet$  **Advection step **

In the second step of advection, the method of characteristics is applied in the very particular case when it is exact. In one time step each particle distribution $\,f_{j}\,$ is exactly transported from the note $\,x-v_{j}\,\Delta t\,$ of the lattice $\,{\cal L}\,$ to the vertex $\,x\in{\cal L}$. In other words, the advection step occurs with a Courant-Friedrichs-Lewy number $\,cfl\equiv 1\,$ for all discrete velocities $\,v_{j}$.

$\bullet$  **Lattice Boltzmann scheme as a Lie-Trotter splitting scheme **

The iteration of the lattice Boltzmann scheme is finally given by the relation

[TABLE]

and it is usefull to explicit it for each component:

[TABLE]

Observe that this framework is very general. The so-called “BGK version” of the lattice Boltzmann scheme [6, 39, 51] corresponds to the choice $\,M={\rm I}\,$ and $\,Z={{1}\over{\tau}}\,{\rm I}$. Usual multi-relaxation lattice Boltzmann schemes [21, 22, 38, 44] suppose in general some orthogonality property for the d’Humières matrix. This paradigm permits also the introduction of two particle distributions [3, 21, 34, 41] or eventually more with the vectorial schemes developed by Graille [30] or a variant proposed in [18]. In those last cases, the numbering $\,j\longmapsto v_{j}\,$ of the velocities is simply non injective.

It is naturally possible to think to a second order accurate splitting scheme. A natural idea is the Strang splitting [55]:

[TABLE]

In this case, as put in evidence by Dellar [12], one can write \,\widetilde{f}(t)\equiv R\big{(}{{\Delta t}\over{2}}\big{)}^{-1}\,f(t)\, and the previous Strang scheme can be written as

[TABLE]

This scheme, up to the shift of one half relaxation step, is identical to the Lie-Trotter splitting (19), except that the product \,R\big{(}{{\Delta t}\over{2}}\big{)}\,R\big{(}{{\Delta t}\over{2}}\big{)}\, plays the role of $\,R(\Delta t)$. A first possibility is to divide the time step $\,\Delta t\,$ by a factor $\,2\,$ in the relations (16)(17):

[TABLE]

Then

[TABLE]

and over-relaxation is no more possible because $\,s_{j}-{1\over 4}\,s_{j}^{2}\leq 1\,$ for all $\,s_{j}\in\mathbb{R}$. So such a numerical scheme is not used by the engineers at our knowledge.

A second possibility is to enforce the relation \,R\big{(}{{\Delta t}\over{2}}\big{)}=\sqrt{R(\Delta t)}, id est

[TABLE]

But this relation has a sense only if $\,s_{j}\leq 1\,$ and over-relaxation is again excluded.

An other idea is to exchange the roles of advection and relaxation in the Strang scheme:

[TABLE]

or to use more elaborated splittings as proposed by Drui et al. [14]. In this case, the transport step is nomore at a Courant number equal to unity and an interpolation procedure, costly in terms of numerical viscosity or complex to implement as the semi-Lagrangian method [27], is not usually chosen by experts developping lattice Boltzmann schemes.

In conclusion of this sub-section, the naive Lie-Trotter splitting scheme (19) gives very good numerical results. It is important to put in evidence precise qualities and defects of this scheme through the following formal asymptotic analysis.

Proposition 2.  Exponential form of lattice Boltzmann schemes

Consider a lattice Boltzmann scheme as defined at the previous section by the relations (2), (4), (5), (7), (8), (14), (15), (16), (18) and (20). Then we have an exponential form of the discrete iteration (20):

[TABLE]

$\bullet$  Proof of Proposition 2.

We have simply a long sequence of identities:

$\displaystyle\,m_{k}(x,t+\Delta t)=\sum_{j}M_{kj}\,f_{j}(x,t+\Delta t)$ due to (2)

$\,\displaystyle=\sum_{j}M_{kj}\,f_{j}^{*}(x-v_{j}\,\Delta t,\,t)$ due to (20)

$\,\displaystyle=\sum_{j\,\ell}M_{kj}\,(M^{-1})_{{}_{\scriptstyle\!j\ell}}\,m_{\ell}^{*}(x-v_{j}\,\Delta t,\,t)$ due to (18)

\,\displaystyle=\sum_{j\,\ell}M_{kj}\,(M^{-1})_{{}_{\scriptstyle\!j\ell}}\,\sum_{n=0}^{\infty}{{1}\over{n!}}\big{(}-\Delta t\sum_{\alpha}v_{j}^{\alpha}\,\partial_{\alpha}\big{)}^{n}\,m_{\ell}^{*}(x,\,t) Taylor expansion

\,\displaystyle=\sum_{\ell}\sum_{n=0}^{\infty}{{1}\over{n!}}\,\sum_{j}M_{kj}\,\Big{(}\!-\Delta t\,\sum_{\alpha}v_{j}^{\alpha}\,\partial_{\alpha}\Big{)}^{n}\,(M^{-1})_{{}_{\scriptstyle\!j\ell}}\,m_{\ell}^{*}(x,\,t) elementary algebra

\,\displaystyle=\sum_{\ell}\sum_{n=0}^{\infty}{{1}\over{n!}}\,(-\Delta t)^{n}\,\bigg{(}\Big{(}M\,\big{(}\sum_{\alpha}v^{\alpha}\,\partial_{\alpha}\big{)}\,M^{-1}\Big{)}^{\!\!n}\bigg{)}_{\scriptstyle\!\!k\ell}\,\,m_{\ell}^{*}(x,\,t)

because $\,M\,$ is a constant matrix

\,\displaystyle=\sum_{\ell}\sum_{n=0}^{\infty}{{1}\over{n!}}\,(-\Delta t)^{n}\,\big{(}\Lambda^{n}\big{)}_{\scriptstyle\!k\ell}\,\,m_{\ell}^{*}(x,\,t) due to (7)

\,\displaystyle=\sum_{\ell}\bigg{[}\,\sum_{n=0}^{\infty}{{1}\over{n!}}\,\big{(}-\Delta t\,\Lambda\big{)}^{n}_{k\ell}\,\bigg{]}\,\,m_{\ell}^{*}(x,\,t) elementary algebra

$\,\displaystyle=\sum_{\ell}\,{\rm exp}(-\Delta t\,\Lambda)_{k\ell}\,\,m_{\ell}^{*}(x,\,t)$ elementary algebra

\,\displaystyle=\Big{(}{\rm exp}(-\Delta t\,\Lambda)\,\,m^{*}(x,\,t)\Big{)}_{\!\!k} elementary algebra

and the proof is completed. $\square$

We can see the relation (21) as a discrete form of the Duhamel formula for ordinary differential equations. The idea of introducing a formal exponential expansion is also present in the work of Boghosian and Coveney in the BGK case [6].

$\bullet$  **Momentum-velocity operator for a D2Q9 lattice Boltzmann scheme **

We can concretize the matrix $\,\Lambda\,$ introduced in (7) with the popular two-dimensional “D2Q9” scheme. A set of nine velocities is defined by the Figure 1.

To fix the ideas, we use the moment matrix presented in Lallemand and Luo [44]:

[TABLE]

Observe that other choices are possible with the Hermite moments developed in [49, 54]. The moments are associated with the lines of the matrix (22). They are denominated (in this order) by $\,\rho,\,J_{x},\,J_{y},\,\varepsilon,\,x_{x},\,x_{y},\,q_{x},\,q_{y},\,h$. The lines of this invertible matrix are chosen orthogonal. For isothermal flows, the conserved moments

[TABLE]

correspond to the three first lines of the d’Humières matrix (22). The non-conserved moments complete the family:

[TABLE]

They are linked to the six last lines of the matrix defined in (22). With the velocities defined in Figure 1 and the moment matrix (22), we construct without difficulty the operator matrix $\,\Lambda\,$ for the thermal D2Q9 scheme:

[TABLE]

In the relation (25), we have put in evidence the block decomposition (8) for the D2Q9 scheme with the conserved moments precised in (23). Even if the space differential operators commute, the matrices that compose the momentum-velocity operator matrix does not commute! We have the very important but elementary structure

[TABLE]

and to fix the ideas

[TABLE]

4)  Taylor expansion method

In this section, we revisit the Taylor expansion method: we specify the hypotheses, express our general result and begin the proof with the order zero.

$\bullet$  **Hypotheses for a formal expansion **

The formal Taylor expansion supposes some precise hypotheses. First, we adopt the acoustic scaling: the ratio $\,\lambda\equiv{{\Delta x}\over{\Delta t}}\,$ is supposed to be constant in all this work. Moreover, the ratio $\,{{\Delta t}\over{\varepsilon}}\,$ remains constant in our asymptotic analysis. Then the relaxation matris $\,S\,$ is fixed and invertible; it is also the case for the matrix $\,S^{-1}\,$ proportional to the matrix $\,Z^{-1}\,$ who plays an important role in the relations (13) of the Chapman-Enskog expansion (11)(12)(13) for the continuous in time lattice Boltzmann-BGK system.

A unintuitive result has been discovered by Hénon [37]. Due to the fact that the fully discretized in time lattice Boltzmannn scheme (20) is substantially different from the continuous in time Boltzmann-BGK system with discrete velocities (9), the asymptotic expansions have similarities and differences. In particular, the lattice Boltzmannn scheme put in evidence what we call the Hénon matrix $\,\Sigma\,$ in this contribution. Roughtly speaking, it plays a role analogous to the matrix $\,Z^{-1}\,$ in the relations (13). But it is defined by

[TABLE]

This matrix emerges from the very classic second order analysis presented thereafter at Proposition 5 and detailed in the relations (40)(41). For applications to fluid dynamics, this matrix is closely related to viscosities, as explained in Section 5 for the D2Q9 example in relations (50) and (51). Then some values of the Hénon matrix $\,\Sigma\,$ are chosen as small as possible in order to simulate flows with high Reynolds number. In consequence, over-relaxation is a mandatory practice for lattice Boltzmann schemes applied to high Reynolds number flows. Recall that this matrix remains fixed in our analysis.

Secondly, we suppose also that we can differentiate all the expansions relative to space and/or time. This hypothesis is mathematically absolutly non trivial and expresses that all the formal expansions should take place in very regular functional spaces. Thirdly, we use the notation $\,{\rm d}\zeta(W).\xi\,$ for the action of the differential of some regular function $\,W\longmapsto\zeta(W)\,$ againts a test vector $\,\xi$. In particular,

[TABLE]

For second order derivatives, we introduce

[TABLE]

The asymptotic analysis occurs for a time step $\,\Delta t\,$ (or a space step $\,\Delta x$) tending to zero. Emerging partial differential equations for the conserved variables are denoted as

[TABLE]

The vector $\,\Gamma_{j}(W)\,$ belongs in $\,\mathbb{R}^{N}$. It is obtained after $j$ space derivations of the conserved moments $\,W\,$ and of the equilibrium vector $\,\Phi(W)$. For this reason, we will speak in the following of *e.g. * $\,\Gamma_{1}(W)\,$ as “first order term”, even if it is a term of order zero relative to the infinitesimal $\,\Delta t$. For non-conserved moments (or microscopic variables), we suppose

[TABLE]

where $\,{\rm I}\,$ is the identity matrix of dimension $\,q-N$. Observe that the two lines of (30) are equivalent thanks to (16) and (26). The vectors \,\Phi_{j}(W)\ are analogous to $\,\Gamma_{j}(W)\,$ but not with the same dimension: $\,\Gamma_{j}(W)\in\mathbb{R}^{N}\,$ and $\,\Phi_{j}(W)\in\mathbb{R}^{q-N}$.

$\bullet$  Main result. Recurrence formulas for fourth order expansion

We prove in this contribution that we have the following recurrence formulas for the explicitation of the vectors $\,\Gamma_{j}(W)\,$ and $\,\Phi_{j}(W)\,$ up to fourth order accuracy:

[TABLE]

Observe that the order of the relations in (31) is mandatory. Each step has to be explicited before the evaluation of the next expression. At second order accuracy, the three first lines of the relations (31) summarizes what is needed for a lot of applications with only first and second order partial differential equations. Observe that the compact form $\,\Gamma_{2}=B\,\Sigma\,\Psi_{1}\,$ for the second order term put in evidence the important choice of a Hénon matrix “as small as possible”. Observe that the relations (31) have similarities with the ones obtained in (13) for the Chapman-Enskog expansion. But essentially the matrix $\,\Sigma\,$ has replaced the matrix $\,Z^{-1}\,$ and new terms appear due to the post-processing with Taylor expansions.

**Proposition 3.  Equilibrium state and zero order expansion **

When $\,\Delta t\,$ tends to zero, the microscopic moments are close to their equilibrium value:

[TABLE]

$\bullet$  Proof of Proposition 3.

We expand one iteration of the scheme (21) at order zero:

[TABLE]

Then we deduce from (33):

[TABLE]

Then due to the iteration (16) and the fact that the matrix $\,S\,$ is supposed fixed, we have

[TABLE]

The first relation of (32) is satisfied. It is then immediate to deduce from the scheme iteration (16) the second relation of (32). The proposition is established. $\square$

In the following, we detail all the ingredients that conduct to the main result (29) to (31).

5)  Taylor expansion method at first and second order accuracy

In this section, we prove the two first orders of the expansion (29)(30)(31). We make the link with the Taylor expansion as presented in [15, 16] and we apply the result for the isothermal D2Q9 scheme.

**Proposition 4.  First order expansion **

When $\,\Delta t\,$ tends to zero, the macroscopic moments $\,W\,$ satisfy asymptotically the following system of $\,N\,$ first order equations:

[TABLE]

Moreover, the vector $\,\Gamma_{1}(W)\,$ for the first order dynamics satisfy

[TABLE]

$\bullet$  Proof of Proposition 4.

We expand one iteration of the scheme (21) at first order:

[TABLE]

We can replace the vector $\,m\,$ by its two components $\,W\,$ and $\,Y$. We have for the first line

[TABLE]

We simplify by the constant $\,W\,$ term, divide by $\,\Delta t\,$ and take into account the expansion (32):

[TABLE]

This relation is exactly the expansion (35) with the first order dynamics $\,\Gamma_{1}(W)\,$ evaluated according to the relation (36).

$\bullet$  Example of the D2Q9 fluid scheme

The results presented here are essentially a reformulation of the classic work of Lallemand and Luo [44]. With the moments introduced in (22), (23) and (24), we denote by $\,\Phi_{\varepsilon}$, $\,\Phi_{xx}\,$ and $\,\Phi_{xy}\,$ the equilibium values of the three first non conserved moments (24). The dynamics at first order is given by (36). We have in this D2Q9 case,

[TABLE]

These relations can be easily fitted with the Euler equations of gas dynamics in two space dimensions

[TABLE]

We identify the two expressions with $\,J_{x}\equiv\rho\,u\,$ and $\,J_{y}\equiv\rho\,v$:

[TABLE]

and we deduce the expression for the three first noncenserved moments:

[TABLE]

The other nonequilibrium moments $\,q_{x}$, $\,q_{y}\,$ and $\,h\,$ does not play any role in the first order partial equivalent equations.

**Proposition 5.  Taylor expansion method at second order accuracy **

An essential result for the applications is the second order asymptotic analysis for a time step $\,\Delta t\,$ (or a space step $\,\Delta x$) tending to zero. The expansion of the microscopic variables takes the form

[TABLE]

and the vector $\,\Psi_{1}(W)\,$ satisfies the relation

[TABLE]

A set of second order partial differential equations emerges

[TABLE]

The vector $\,\Gamma_{1}(W)\,$ has been precised in (36) the vector $\,\Gamma_{2}(W)\,$ is obtained after two derivations of the conserved moments $\,W\,$ and the equilibrium vector $\,\Phi(W)$:

[TABLE]

The expression of the second order term puts in evidence the importance of taking over-relaxation in the applications where $\,\Sigma\,$ should also be as small as possible.

$\bullet$  Proof of Proposition 5.

We first consider the first order expansion of microscopic moments. From the expansion (37) of one iteration of the scheme at first order, we can extract relations for the second component. We have for the microscopic moments

[TABLE]

Then

[TABLE]

The explicitation of $\,\partial_{t}Y\,$ is given by the chain rule from the expansion (34):

\displaystyle\partial_{t}Y=\partial_{t}\big{(}\Phi(W)+{\rm O}(\Delta t)\big{)}

$\qquad\displaystyle={\rm d}\Phi(W)\,.\,\partial_{t}W+{\rm O}(\Delta t)$

$\qquad\displaystyle={\rm d}\Phi(W).(-\Gamma_{1})+{\rm O}(\Delta t)\,.$

Then

[TABLE]

The relation (38) is established, and the expression of $\,\Psi_{1}(W)\,$ is obtained by the relation (39). Before using the expansion (38), we must expand the microscopic moments $\,Y^{*}\,$ after relaxation up to first order accuracy. We have the following calculus

$\displaystyle Y^{*}=Y+S\,(\Phi(W)-Y)$

\,\,\,\displaystyle=Y-\big{(}\Delta t\,\Psi_{1}+{\rm O}(\Delta t^{2})\big{)} due to (38)

\,\,\,\displaystyle=\Phi(W)+\Big{(}\Sigma+{{1}\over{2}}\,{\rm I}\Big{)}\,\big{(}\Delta t\,\Psi_{1}(W)+{\rm O}(\Delta t^{2})\big{)}-\big{(}\Delta t\,\Psi_{1}+{\rm O}(\Delta t^{2})\big{)} due to (26) and (38)

\,\,\,\displaystyle=\Phi(W)+\Big{(}\Sigma-{{1}\over{2}}\,{\rm I}\Big{)}\,\Delta t\,\Psi_{1}(W)+{\rm O}(\Delta t^{2})

and we have the relation

[TABLE]

The two relations (38) and (42) are nothing else that the first order terms of the general expansion (30).

$\bullet$  Second order partial differential equations

We expand now one iteration of the scheme (21) at second order accuracy:

[TABLE]

We replace the vector $\,m\,$ by its two components $\,W\,$ and $\,Y$. We use the decomposition (8) of the momentum-velocity operator matrix $\,\Lambda$. We have for the first component

[TABLE]

We explicit some terms present in the relation (44). We have

[TABLE]

Then, using the formal rule stated at the beginning of Section 4, we keep the order of accuracy after derivation:

\displaystyle\partial^{2}_{t}W=-\partial_{t}\big{(}\Gamma_{1}+{\rm O}(\Delta t)\big{)}

$\qquad\,\,\displaystyle=-{\rm d}\Gamma_{1}\,.\,\partial_{t}W+{\rm O}(\Delta t)$

$\qquad\,\,\displaystyle={\rm d}\Gamma_{1}\,.\,\Gamma_{1}+{\rm O}(\Delta t)$

$\qquad\,\,\displaystyle=A\,\Gamma_{1}+B\,{\rm d}\Phi\,.\,\Gamma_{1}+{\rm O}(\Delta t)\,.$

We deduce from (44):

[TABLE]

We order the various terms by powers of $\,\Delta t$:

\displaystyle\partial_{t}W=-A\,W-B\,\Big{(}\Phi+(\Sigma-{1\over 2}{\rm I})\,\Delta t\,\Psi_{1}\Big{)}-{1\over 2}\,\Delta t\,(A\,\Gamma_{1}+B\,{\rm d}\Phi.\Gamma_{1})

\qquad\qquad\qquad\displaystyle+{1\over 2}\,\Delta t\,\big{(}(A^{2}+B\,C)\,W+(A\,B+B\,D)\,\Phi\big{)}+{\rm O}(\Delta t^{2})

\qquad\,\displaystyle=-A\,W-B\,\Phi+\Delta t\,\Big{[}-B\,\Sigma\,\Psi_{1}+{1\over 2}B\,\big{(}{\rm d}\Phi\,.\,\Gamma_{1}-C\,W-D\,\Phi\big{)}-{1\over 2}A\,\,(A\,W+B\,\Phi)

\qquad\qquad\qquad\displaystyle-{1\over 2}\,B\,{\rm d}\Phi.\Gamma_{1}+{1\over 2}\,(A^{2}+B\,C)\,W+{1\over 2}\,(A\,B+B\,D)\,\Phi\Big{]}+{\rm O}(\Delta t^{2})

$\qquad\,\displaystyle=-\Gamma_{1}(W)-\Delta t\,B\,\,\Sigma\,\,\Psi_{1}(W)+{\rm O}(\Delta t^{2})$

and the relation (40) is proven with $\,\Gamma_{2}(W)=B\,\,\Sigma\,\,\Psi_{1}\,$ as proposed in (41). $\square$

$\bullet$  **Link with the original Taylor expansion method **

In [15, 16], the moments at equilibrium are denoted $\,m_{\ell}^{\rm eq}$:

[TABLE]

the coefficients $\,\widetilde{\Lambda}_{i\ell}^{\beta}\,$ are given according to

[TABLE]

are closely related to the operator matrix $\,\Lambda\,$ introduced in (7). The coefficients of the Hénon’s matrix (26) have a shifted numbering

[TABLE]

and the defect of conservation $\,\widetilde{\theta}_{k}\,$ is defined according to

[TABLE]

**Proposition 6.  Equivalence at second order accuracy **

The second order expansion (40), (38), (36), (39) and (41) is equivalent to the expansion proposed in [15, 16]:

[TABLE]

$\bullet$  Proof of Proposition 6.

Remark first that

[TABLE]

Then on one hand, the component number $\,i\,$ of the first order term \,\big{(}A\,W+B\,\Phi(W)\big{)}\, can be expanded in the following way:

\displaystyle\big{(}A\,W+B\,\Phi(W)\big{)}_{i}=\sum_{k<N}\widetilde{\Lambda}_{ik}^{\beta}\,\,\partial_{\beta}W_{k}+\sum_{k\geq N}\widetilde{\Lambda}_{ik}^{\beta}\,\,\partial_{\beta}\Phi_{k}

$\displaystyle\qquad\qquad\qquad\qquad\,=\sum_{k<N}\widetilde{\Lambda}_{ik}^{\beta}\,\,\partial_{\beta}m_{k}^{\rm eq}+\sum_{k\geq N}\widetilde{\Lambda}_{ik}^{\beta}\,\,\partial_{\beta}m_{k}^{\rm eq}$

$\displaystyle\qquad\qquad\qquad\qquad\,=\sum_{0\leq k<q}\widetilde{\Lambda}_{ik}^{\beta}\,\,\partial_{\beta}m_{k}^{\rm eq}\,$

and the first order terms of (40) and (48) are identical. Secondly, due to (47), we have for $\,k\geq N\,$

$\displaystyle\widetilde{\theta}_{k}=\partial_{t}m_{k}^{\rm eq}+\sum_{\ell\beta}\widetilde{\Lambda}_{k\ell}^{\beta}\,\,\partial_{\beta}m_{\ell}^{\rm eq}$

\displaystyle\quad\,=\partial_{t}\Phi_{k-N}+\big{(}\Lambda.m^{\rm eq}\big{)}_{k}

\displaystyle\quad\,=\big{(}{\rm d}\Phi\,.\,\partial_{t}W\big{)}_{k-N}+\big{(}\Lambda.m^{\rm eq}\big{)}_{k}

\displaystyle\quad\,=\big{(}{\rm d}\Phi\,.\,(\Gamma_{1}+{\rm O}(\Delta t))\big{)}_{k-N}-\big{(}C\,W+D\,\Phi\big{)}_{k-N}

and due to (39),

[TABLE]

We deduce that for second order terms with $\,0\leq i<N\,$:

\displaystyle\big{(}B\,\,\widetilde{\sigma}\,\,\Psi_{1}\big{)}_{i}=\sum_{0\leq\ell<q-N}B_{i\ell}\,\widetilde{\sigma}_{\ell+N}\,\,\big{(}\Psi_{1}(W)\big{)}_{\ell}

\displaystyle\qquad\qquad\,\,\,\,=\sum_{k\,\geq N}\,\sum_{0\leq\beta\leq d}\widetilde{\Lambda}_{ik}^{\beta}\,\,\partial_{\beta}\Big{(}\sigma_{k}\,\big{(}-\widetilde{\theta}_{k}+{\rm O}(\Delta t)\big{)}\Big{)} with $\,k=\ell+N$

\displaystyle\qquad\qquad\,\,\,\,=-\sum_{k\,\geq N}\,\sum_{0\leq\beta\leq d}\widetilde{\Lambda}_{ik}^{\beta}\,\widetilde{\sigma}_{k}\,\big{(}\partial_{\beta}\widetilde{\theta}_{k}\big{)}+{\rm O}(\Delta t)

because the coefficients $\,\widetilde{\sigma}_{k}\,$ are supposed constant. Then the relation (48) is established. $\square$

$\bullet$  Fluid D2Q9 diffusive tensor

If we compute the second order term $\,\Gamma_{2}\,$ for the D2Q9 scheme introduced previously, the holy grail would be to recover the viscous terms of the compressible Navier Stokes equations in two space dimensions. More precisely, with a given shear viscosity $\,\mu\,$ and a given bulk viscosity $\,\zeta$, we write the viscous terms of the lattice Boltzmann expansion as

[TABLE]

with the notation $\,\Psi_{x}$, $\,\Psi_{y}\,$ for the error for the dissipation of momentum.

The problem to enforce $\,\Psi_{x}=\Psi_{y}=0\,$ is ill-posed, as recognized by Dellar [11, 13], Házi and Kávrán [35] and Prasianakis at al. [50] among others. Nevertheless, it can be proven that the unphysical terms involving second order space derivatives involving the density can be removed with a pressure law of the type

[TABLE]

id est a sound velocity equal to $\,{{1}\over{\sqrt{3}}}\,$ in the velocity unit of the lattice Boltzmann scheme. Nevertheless, with the notation $\,\Phi_{qx}\,$ and $\,\Phi_{qy}\,$ for the equilibrium fonctions of the moments $\,q_{x}\,$ and $\,q_{y}\,$ and the coefficients of the Hénon matrix $\,\Sigma\,$ fixed as

[TABLE]

the choice

[TABLE]

conducts to the relation (49) with associated viscosities fixed in a very classical way:

[TABLE]

The equilibrium value $\,\Phi_{h}\,$ of the last fourth order moment $\,h\,$ (see the relation (24)) has no incidence on the value of the diffusive term $\,\Gamma_{2}$. The error is then reduced to third order terms relative to the velocity and we have precisely

[TABLE]

6)  Taylor expansion method at third and fourth order accuracy

The two following orders of the Taylor expansion does not set any theoretical difficulty, except the care of algebra with not so short formal expressions…

**Proposition 7.  Nonlinear third-order expansion **

When the time step $\,\Delta t\,$ has an infinitesimal value, the expansion of the microscopic variables takes the form

[TABLE]

The vector $\,\Psi_{1}(W)\,$ is evaluated at the relation (39) and $\,\Psi_{2}(W)\,$ satisfies the relation

[TABLE]

It is possible to precise a system of third order partial differential equations

[TABLE]

The vectors $\,\Gamma_{1}(W)\,$ and $\,\Gamma_{2}(W)\,$ have been precised at the relations (36) and (41) respectively. We have moreover

[TABLE]

$\bullet$  Proof of Proposition 7.

We first establish the relations (52) and (53). We start by the expansion of the scheme (43) at order 2. For the second component of nonconserved moments, we have

$Y+\Delta t\,\partial_{t}Y+{1\over 2}\,\Delta t^{2}\,\partial_{t}^{2}Y+{\rm O}(\Delta t^{3})=Y^{*}-\Delta t\,(C\,W+D\,Y^{*})+{1\over 2}\,\Delta t^{2}\,(C_{2}\,W+D_{2}\,Y^{*})+{\rm O}(\Delta t^{3})$

and, due to the relaxation process $\,Y-Y^{*}\equiv S\,(Y-\Phi(W))$, we have

[TABLE]

We can precise the two first time derivatives of the non-conserved moments. The relation (38) can also be written as

[TABLE]

We differentiate this previous relation relative to time:

\,\partial_{t}Y={\rm d}\Phi.\,\partial_{t}W+\big{(}\Sigma+{{1}\over{2}}\,{\rm I}\big{)}\,\,\Delta t\,\,{\rm d}\Psi_{1}.\,\partial_{t}W+{\rm O}(\Delta t^{2})

=-{\rm d}\Phi.\,\big{(}\Gamma_{1}(W)+\Delta t\,\Gamma_{2}(W)+{\rm O}(\Delta t^{2})\big{)}+\Delta t\,\big{(}\Sigma+{{1}\over{2}}\,{\rm I}\big{)}\,{\rm d}\Psi_{1}.\,\big{(}-\Gamma_{1}+{\rm O}(\Delta t)\big{)}+{\rm O}(\Delta t^{2})

and

[TABLE]

Because we need the second time derivative $\,\partial_{t}^{2}Y\,$ in the right hand side of (56), we differentiate the relation (57) relative to time. Using the notation (27), we deduce

[TABLE]

We precise now the expansions of $\,C\,W+D\,Y^{*}\,$ and $\,C_{2}\,W+D_{2}\,Y^{*}$. We have

C\,W+D\,Y^{*}=C\,W+D\,\big{[}\Phi(W)+\big{(}\Sigma-{{1}\over{2}}\,{\rm I}\big{)}\,\Delta t\,\Psi_{1}\big{]}+{\rm O}(\Delta t^{2})

and

[TABLE]

With a new set of operator matrices,

C_{2}\,W+D_{2}\,Y^{*}=C_{2}\,W+D_{2}\big{(}\Phi(W)+{\rm O}(\Delta t)\big{)}

$\,\,\,=(C\,A+D\,C)\,W+(C\,B+D^{2})\,\Phi(W)+{\rm O}(\Delta t)$

$\,\,\,=C\,\Gamma_{1}+D\,({\rm d}\Phi.\Gamma_{1}-\Psi_{1})+{\rm O}(\Delta t)$

and because $\,\,\,{\rm d}\Phi.\Gamma_{1}={\rm d}\gamma_{1}.\Gamma_{1}-(C\,\Gamma_{1}+D\,{\rm d}\Phi.\Gamma_{1})$,

[TABLE]

Then due to (56), (57), (58), (59) and (60), we have

$S\,(Y-\Phi(W))=-\Delta t\,(C\,W+D\,Y^{*}+\partial_{t}Y)+{1\over 2}\,\Delta t^{2}\,(C_{2}\,W+D_{2}\,Y^{*}-\partial_{t}^{2}Y)+{\rm O}(\Delta t^{3})$

=-\Delta t\,\big{(}\gamma_{1}-\Psi_{1}+\Delta t\,\,D\,\big{(}\Sigma-{{1}\over{2}}\,{\rm I}\big{)}\,\Psi_{1}\big{)}-\,\Delta t\,\big{(}-{\rm d}\Phi.\Gamma_{1}-\Delta t\,\big{[}{\rm d}\Phi.\Gamma_{2}+\big{(}\Sigma+{{1}\over{2}}\,{\rm I}\big{)}\,\,{\rm d}\Psi_{1}.\Gamma_{1}\big{]}\big{)}

+{1\over 2}\,\Delta t^{2}\,\big{(}{\rm d}(\gamma_{1}-\Psi_{1}).\Gamma_{1}-D\,\Psi_{1}\big{)}-{1\over 2}\,\Delta t^{2}\,\,{\rm d}\gamma_{1}.\,\Gamma_{1}+{\rm O}(\Delta t^{3})

=\Delta t\,\Psi_{1}+\Delta t^{2}\,\big{(}-D\,\Sigma\,\Psi_{1}+{1\over 2}\,D\,\Psi_{1}+{\rm d}\Phi.\Gamma_{2}+\Sigma\,{\rm d}\Psi_{1}.\Gamma_{1}-{1\over 2}\,D\,\Psi_{1}\big{)}+{\rm O}(\Delta t^{3})

=\Delta t\,\Psi_{1}+\Delta t^{2}\,\big{(}-D\,\Sigma\,\Psi_{1}+{\rm d}\Phi.\Gamma_{2}+\Sigma\,{\rm d}\Psi_{1}.\Gamma_{1}\big{)}+{\rm O}(\Delta t^{3})

and the relations (52) and (53) are established.

$\bullet$  Third order partial differential equations

We consider now the Taylor expansion of the relation (21) at third order:

[TABLE]

We consider the first conserved components of the moments:

[TABLE]

We simplify the constant term $\,W\,$ and divide by $\,\Delta t$. We deduce

[TABLE]

We explicit the partial derivatives $\,\partial_{t}^{2}W\,$ and $\,\partial_{t}^{3}W\,$ at orders one and zero respectively. From the relation

[TABLE]

we have

\partial_{t}^{2}W=\partial_{t}\,\big{(}-\Gamma_{1}-\Delta t\,\Gamma_{2}+{\rm O}(\Delta t^{2})\big{)}

\,\,={\rm d}\big{(}-\Gamma_{1}-\Delta t\,\Gamma_{2}+{\rm O}(\Delta t^{2})\big{)}.(\partial_{t}W)

\,\,={\rm d}\big{(}-\Gamma_{1}-\Delta t\,\Gamma_{2}+{\rm O}(\Delta t^{2})\big{)}\,.\,\big{(}-\Gamma_{1}-\Delta t\,\Gamma_{2}+{\rm O}(\Delta t^{2})

and

[TABLE]

Then

$\,\partial_{t}^{3}W=\partial_{t}(\partial_{t}^{2}W)$

\,\,\,=\partial_{t}\big{(}{\rm d}\Gamma_{1}.\,\Gamma_{1}\big{)}+{\rm O}(\Delta t)

\,\,\,={\rm d}\big{(}{\rm d}\Gamma_{1}.\,\Gamma_{1}\big{)}.\partial_{t}W+{\rm O}(\Delta t)

\,\,\,=-{\rm d}\big{(}{\rm d}\Gamma_{1}.\,\Gamma_{1}\big{)}.\Gamma_{1}+{\rm O}(\Delta t)

[TABLE]

In order to compute the coefficient $\,\Gamma_{3}$, we precise now the expansions of the quantities $\,A_{j}\,W+B_{j}\,Y^{*}\,$ for $\,j=1$, $\,2\,$ and $\,3$. We have

A\,W+B\,Y^{*}=A\,W+B\,\big{(}\Phi(W)+\big{(}\Sigma-{{1}\over{2}}\,{\rm I}\big{)}\,\big{(}\Delta t\,\Psi_{1}+\Delta t^{2}\,\Psi_{2}\big{)}+{\rm O}(\Delta t^{3})

\,\,\,\,=A\,W+B\,\Phi+\Delta t\,B\,\big{(}\Sigma-{{1}\over{2}}\,{\rm I}\big{)}\,\Psi_{1}+\Delta t^{2}\,B\,\big{(}\Sigma-{{1}\over{2}}\,{\rm I}\big{)}\,\Psi_{2}+{\rm O}(\Delta t^{3})

A_{2}\,W+B_{2}\,Y^{*}=A_{2}\,W+B_{2}\,\big{(}\Phi(W)+\big{(}\Sigma-{{1}\over{2}}\,{\rm I}\big{)}\,\Delta t\,\Psi_{1}\big{)}+{\rm O}(\Delta t^{2})

\,\,\,=(A^{2}+B\,C)\,W+(A\,B+B\,D)\,\Phi+\Delta t\,B_{2}\,\big{(}\Sigma-{{1}\over{2}}\,{\rm I}\big{)}\,\Psi_{1}+{\rm O}(\Delta t^{2})

\,\,\,=A\,\Gamma_{1}+B\,(\gamma_{1}-\Psi_{1})+\Delta t\,B_{2}\,\big{(}\Sigma-{{1}\over{2}}\,{\rm I}\big{)}\,\Psi_{1}+{\rm O}(\Delta t^{2})

because $\,\,C\,W+D\,\Psi_{1}=\gamma_{1}-\Psi_{1}$

\,\,\,={\rm d}\Gamma_{1}.\Gamma_{1}-B\,\Psi_{1}+\Delta t\,B_{2}\,\big{(}\Sigma-{{1}\over{2}}\,{\rm I}\big{)}\,\Psi_{1}+{\rm O}(\Delta t^{2})

because $\,\,{\rm d}\Gamma_{1}.\Gamma_{1}=A\,\Gamma_{1}+B\,\gamma_{1}$,

$A_{3}\,W+B_{3}\,Y^{*}=A_{3}\,W+B_{3}\,\Phi(W)+{\rm O}(\Delta t)$

$\,\,\,=(A_{2}\,A+B_{2}\,C)\,W+(A_{2}\,B+B_{2}\,D)\,\Phi+{\rm O}(\Delta t)$

$\,\,\,=A_{2}\,\Gamma_{1}+B_{2}\,(\gamma_{1}-\Psi_{1})+{\rm O}(\Delta t)$

$\,\,\,=(A^{2}+B\,C)\,\Gamma_{1}+(A\,B+B\,D)\,\gamma_{1}-B_{2}\,\Psi_{1}+{\rm O}(\Delta t)$

$\,\,\,=A\,(A\,\Gamma_{1}+B\,\gamma_{1})+B\,(C\,\Gamma_{1}+D\,\gamma_{1})-B_{2}\,\Psi_{1}+{\rm O}(\Delta t)$

$\,\,\,=A\,{\rm d}\Gamma_{1}.\Gamma_{1}+B\,(C\,\Gamma_{1}+D\,\gamma_{1})-B_{2}\,\Psi_{1}+{\rm O}(\Delta t)$

because $\,\,{\rm d}\Gamma_{1}.\Gamma_{1}=A\,\Gamma_{1}+B\,\gamma_{1}$

$\,\,\,=A\,{\rm d}\Gamma_{1}.\Gamma_{1}+B\,{\rm d}(\gamma_{1}-\Psi_{1}).\Gamma_{1}-B_{2}\,\Psi_{1}+{\rm O}(\Delta t)$

because $\,\,C\,\Gamma_{1}+D\,\gamma_{1}={\rm d}(\gamma_{1}-\Psi_{1}).\Gamma_{1}$

$\,\,\,=\partial^{2}\Gamma_{1}.\Gamma_{1}-B\,{\rm d}\Psi_{1}.\Gamma_{1}-B_{2}\,\Psi_{1}+{\rm O}(\Delta t)$

because $\,\,\partial^{2}\Gamma_{1}.\Gamma_{1}\equiv{\rm d}({\rm d}\Gamma_{1}.\Gamma_{1}).\Gamma_{1}=A\,{\rm d}\Gamma_{1}.\Gamma_{1}+B\,{\rm d}\gamma_{1}.\Gamma_{1}$.

We deduce now from (62), (63), (64) and the previous expressions:

$\partial_{t}W=-A\,W-B\,Y^{*}+{1\over 2}\,\Delta t\,(A_{2}\,W+B_{2}\,Y^{*}-\partial_{t}^{2}W)-{1\over 6}\,\Delta t^{2}\,(A_{3}\,W+B_{3}\,Y^{*}+\partial_{t}^{3}W)+{\rm O}(\Delta t^{3})$

\,=-\big{(}A\,W+B\,\Phi+\Delta t\,B\,\big{(}\Sigma-{{1}\over{2}}\,{\rm I}\big{)}\,\Psi_{1}+\Delta t^{2}\,B\,\big{(}\Sigma-{{1}\over{2}}\,{\rm I}\big{)}\,\Psi_{2}\big{)}

\qquad+{1\over 2}\,\Delta t\,\big{(}{\rm d}\Gamma_{1}.\Gamma_{1}-B\,\Psi_{1}+\Delta t\,B_{2}\,\big{(}\Sigma-{{1}\over{2}}\,{\rm I}\big{)}\,\Psi_{1}-{\rm d}\Gamma_{1}.\,\Gamma_{1}-\Delta t\,({\rm d}\Gamma_{1}.\,\Gamma_{2}+{\rm d}\Gamma_{2}.\,\Gamma_{1})\big{)}

\qquad-{1\over 6}\,\Delta t^{2}\,\big{(}\partial^{2}\Gamma_{1}.\Gamma_{1}-B\,{\rm d}\Psi_{1}.\Gamma_{1}-B_{2}\,\Psi_{1}-\partial^{2}\Gamma_{1}.\,\Gamma_{1}\big{)}+{\rm O}(\Delta t^{3})

\,=-\Gamma_{1}-\Delta t\,B\,\Sigma\,\Psi_{1}+\Delta t^{2}\,\big{(}-B\,\Sigma\,\Psi_{2}+{1\over 2}\,B\,\Psi_{2}+{1\over 2}\,B_{2}\,\Sigma\,\Psi_{1}-{1\over 4}\,B_{2}\,\Psi_{1}-{1\over 2}\,({\rm d}\Gamma_{1}.\,\Gamma_{2}+{\rm d}\Gamma_{2}.\,\Gamma_{1})

\qquad+{1\over 6}\,B\,{\rm d}\Psi_{1}.\Gamma_{1}+{1\over 6}\,B_{2}\,\Psi_{1}\big{)}+{\rm O}(\Delta t^{3}) because $\,\,\partial^{2}\Gamma_{1}.\Gamma_{1}=A\,{\rm d}\Gamma_{1}.\Gamma_{1}+B\,{\rm d}\gamma_{1}.\Gamma_{1}$

\,=-\Gamma_{1}-\Delta t\,B\,\Sigma\,\Psi_{1}+\Delta t^{2}\,\big{(}-B\,\Sigma\,\Psi_{2}+{1\over 2}\,B\,(\Sigma\,{\rm d}\Psi_{1}.\Gamma_{1}+\gamma_{2}-D\,\Sigma\,\Psi_{1})+{1\over 2}\,(A\,B+B\,D)\,\Sigma\,\Psi_{1}

\qquad-{1\over 12}\,B_{2}\,\Psi_{1}-{1\over 2}\,(A\,\Gamma_{2}+B\,\gamma_{2})-{1\over 2}\,B\,\Sigma\,{\rm d}\Psi_{1}.\Gamma_{1}+{1\over 6}\,B\,{\rm d}\Psi_{1}.\Gamma_{1}\big{)}+{\rm O}(\Delta t^{3})

\,=-\Gamma_{1}-\Delta t\,\Gamma_{2}-\Delta t^{2}\,\big{(}B\,\Sigma\,\Psi_{2}+{1\over 12}\,B_{2}\,\Psi_{1}-{1\over 6}\,B\,{\rm d}\Psi_{1}.\Gamma_{1}\big{)}+{\rm O}(\Delta t^{3})

and the relations (54) and (55) are established. $\square$

In [17], we have proposed a third order expansion of nonlinear lattice Boltzmann schemes for two specific applications: scalar advection diffusion ($N=1$) and fluid system ($N=d+1$) with mass and momentum conservation. In the following proposition, we show that our previous result at third order accuracy can be reformulated with the compact formulation (54) (55) of Proposition 6.

**Proposition 8.  Third-order formal expansion for thermics and fluids **

With the notations (45), (46) and (47) introduced at Proposition 5, we suppose moreover

[TABLE]

The first moments are denoted by $\,\rho\equiv\sum_{j}f_{j}\,$ and $\,J_{\alpha}\equiv\sum_{j}v_{j}^{\alpha}\,f_{j}$. For the scalar advection diffusion ($N=1$), the third order equivalent partial differential equation (54) (55) can be written

[TABLE]

with $\,J_{\alpha}^{\rm eq}\equiv\Phi_{\alpha}(\rho)\,$ and the defects of conservation $\,\widetilde{\theta}_{k}(W)\,$ defined in (47). For the fluid system ($N=d+1$), we have

[TABLE]

The relations (65) and (66) are axactly the ones proposed with numbers (35), (40) and (41) in the reference [17].

$\bullet$  Proof of Proposition 8.

We first consider the conservation defect $\,\widetilde{\theta}_{k}(W)\equiv\partial_{t}m_{k}^{\rm eq}+\sum_{\ell\beta}\widetilde{\Lambda}_{k\ell}^{\beta}\,\partial_{\beta}m_{\ell}^{\rm eq}\,$ introduced in (47). We define a vector conservation defect $\,\theta\,$ of dimension $\,q-N\,$ just obtained by a shift of component numbering:

[TABLE]

We have for $\,0\leq k<q-N$:

$\theta_{k}=({\rm d}\Phi.\,\partial_{t}W)_{k}+(\widetilde{\Lambda}.m^{\rm eq})_{k+N}$

$\quad\,=({\rm d}\Phi.\,\partial_{t}W)_{k}+(C\,W+D\,\Phi(W))_{k}$

$\quad\,=({\rm d}\Phi.\,(-\Gamma_{1}-\Delta t\,\Gamma_{2}))_{k}+(C\,W+D\,\Phi(W))_{k}+{\rm O}(\Delta t^{2})$ due to the expansion (29)

\quad\,=-\big{(}\Psi_{1}+\Delta t\,{\rm d}\Phi.\,\Gamma_{2}\big{)}_{k}+{\rm O}(\Delta t^{2})\, due to the relation (31)

and

[TABLE]

In particular, $\,\partial_{t}\theta=-\partial_{t}\Psi_{1}+{\rm O}(\Delta t)=-{\rm d}\Psi_{1}.\,(-\Gamma_{1})+{\rm O}(\Delta t)$ and

[TABLE]

Due to (41) and (55), we have the following calculus:

$\,\Gamma_{2}+\Delta t\,\Gamma_{3}=B\,\Sigma\,\Psi_{1}+\Delta t\,(B\,\Sigma\,\Psi_{2}-{1\over 12}\,B_{2}\,\Psi_{1}+{{1}\over{6}}\,B\,{\rm d}\Psi_{1}.\,\Gamma_{1})$

\qquad\qquad\quad=B\,\Sigma\,(-\theta+\Delta t\,\,{\rm d}\Phi.\,\Gamma_{2}+{\rm O}(\Delta t^{2}))+\Delta t\,B\,\Sigma\,\big{(}D\,\Sigma\,\Psi_{1}-{\rm d}\Phi.\,\Gamma_{2}-\Sigma\,{\rm d}\Psi_{1}.\,\Gamma_{1}\big{)}

$-{1\over 12}\,\Delta t\,B_{2}\,\Psi_{1}+{{1}\over{6}}\,\Delta t\,B\,{\rm d}\Psi_{1}.\,\Gamma_{1}$

\qquad\qquad\quad=-B\,\Sigma\,\theta+\Delta t\,\big{[}(B\,\Sigma\,D\,\Sigma\,-{1\over 12}\,B_{2})\,\Psi_{1}+B\,({{1}\over{6}}-\Sigma^{2})\,{\rm d}\Psi_{1}.\,\Gamma_{1}\big{]}+{\rm O}(\Delta t^{2})

\qquad\qquad\quad=-B\,\Sigma\,\theta+\Delta t\,\big{[}\,-(B\,\Sigma\,D\,\Sigma\,-{1\over 12}\,B_{2})\,\theta+B\,(\Sigma^{2}-{{1}\over{6}})\,\partial_{t}\theta\,\big{]}+{\rm O}(\Delta t^{2})\,.

The third-order partial equivalent equations (54) can be written under the form

[TABLE]

$\bullet$  We explicit now the following cartesian components:

$(B\,\Sigma\,D\,\Sigma\,\theta)_{{}_{\scriptstyle i}}=-\widetilde{\Lambda}_{ik}^{\beta}\,\,\widetilde{\sigma}_{k}\,(-\widetilde{\Lambda}_{k\ell}^{\gamma})\,\widetilde{\sigma}_{\ell}\,\,\partial_{\beta}\partial_{\gamma}\widetilde{\theta}_{\ell}=\widetilde{\Lambda}_{ik}^{\beta}\,\,\widetilde{\sigma}_{k}\,\,\widetilde{\Lambda}_{k\ell}^{\gamma}\,\widetilde{\sigma}_{\ell}\,\,\partial_{\beta}\partial_{\gamma}\widetilde{\theta}_{\ell}\,,$

$(B_{2}\,\theta)_{{}_{\scriptstyle i}}=(\widetilde{\Lambda}_{ik}^{\beta}\,\partial_{\beta})\,(\widetilde{\Lambda}_{k\ell}^{\gamma}\,\partial_{\gamma})\,\widetilde{\theta}_{\ell}=\widetilde{\Lambda}_{ik}^{\beta}\,\widetilde{\Lambda}_{k\ell}^{\gamma}\,\,\partial_{\beta}\partial_{\gamma}\widetilde{\theta}_{\ell}\,,$

((B\,\Sigma\,D\,\Sigma\,-{1\over 12}\,B_{2})\,\widetilde{\theta})_{{}_{\scriptstyle i}}=\widetilde{\Lambda}_{ik}^{\beta}\,\widetilde{\Lambda}_{k\ell}^{\gamma}\,\big{(}\widetilde{\sigma}_{k}\,\widetilde{\sigma}_{\ell}-{1\over 12}\big{)}\,\partial_{\beta}\partial_{\gamma}\widetilde{\theta}_{\ell}

and \,\,-(B\,(\Sigma^{2}-{{1}\over{6}})\,\partial_{t}\widetilde{\theta})_{{}_{\scriptstyle i}}=\widetilde{\Lambda}_{ik}^{\beta}\,\big{(}\widetilde{\sigma}_{k}^{2}\,-{1\over 6}\big{)}\,\partial_{\beta}\partial_{t}\widetilde{\theta}_{k}.

We deduce a new form of the third-order partial equivalent equations for $\,0\leq i<N$:

[TABLE]

$\bullet$  We focus on the mass conservation ($i=0$). We have

$\widetilde{\Lambda}_{0\ell}^{\alpha}=M_{0j}\,v_{j}^{\alpha}\,(M^{-1})_{{}_{\scriptstyle j\ell}}=M_{\alpha j}\,(M^{-1})_{{}_{\scriptstyle j\ell}}=\delta_{\alpha\ell}$.

Then

$\widetilde{\Lambda}_{0k}^{\beta}\,\widetilde{\sigma}_{k}\,\partial_{\beta}\,\widetilde{\theta}_{k}=\delta_{\beta k}\,\widetilde{\sigma}_{k}\,\partial_{\beta}\,\widetilde{\theta}_{k}=\widetilde{\sigma}_{\beta}\,\partial_{\beta}\,\widetilde{\theta}_{\beta}$.

In the thermal case, we have only one conservation law and $\,\widetilde{\sigma}_{\beta}\not=0$. Thus we have enlightened the second order term of the left hand side of the relation (65). In the fluid case, the momentum (id est moments numbered from $\,1\,$ to $\,d$) are conserved. In consequence, $\,\widetilde{\sigma}_{\beta}=0\,$ and we have no second order term in the mass conservation (66).

Now for the first term relative to third order term of (67), we have

\widetilde{\Lambda}_{0k}^{\beta}\,\widetilde{\Lambda}_{k\ell}^{\gamma}\,\big{(}\widetilde{\sigma}_{k}\,\widetilde{\sigma}_{\ell}-{1\over 12}\big{)}\,\partial_{\beta}\partial_{\gamma}\theta_{\ell}=\delta_{\beta k}\,\widetilde{\Lambda}_{k\ell}^{\gamma}\,\big{(}\widetilde{\sigma}_{k}\,\widetilde{\sigma}_{\ell}-{1\over 12}\big{)}\,\partial_{\beta}\partial_{\gamma}\theta_{\ell}=\widetilde{\Lambda}_{\beta\ell}^{\gamma}\,\big{(}\widetilde{\sigma}_{\beta}\,\widetilde{\sigma}_{\ell}-{1\over 12}\big{)}\,\partial_{\beta}\partial_{\gamma}\theta_{\ell}.

In the termal case, $\,\widetilde{\sigma}_{\beta}\not=0\,$ and this term is exactly the first of the two third order terms of the relation (65). In the fluid case, $\,\widetilde{\sigma}_{\beta}=0\,$ and this term has a more compact expression. Finally,

[TABLE]

For the second term relative to third order, we have

\widetilde{\Lambda}_{0k}^{\beta}\,\big{(}\widetilde{\sigma}_{k}^{2}\,-{1\over 6}\big{)}\,\partial_{\beta}\partial_{t}\theta_{k}=\big{(}\widetilde{\sigma}_{\beta}^{2}\,-{1\over 6}\big{)}\,\partial_{t}(\partial_{\beta}\theta_{\beta})

In the thermal case, $\,\widetilde{\sigma}_{\beta}\not=0\,$ and \widetilde{\Lambda}_{0k}^{\beta}\,\big{(}\widetilde{\sigma}_{k}^{2}\,-{1\over 6}\big{)}\,\partial_{\beta}\partial_{t}\theta_{k}=\big{(}\widetilde{\sigma}_{\beta}^{2}\,-{1\over 6}\big{)}\,\partial_{t}(\partial_{\beta}\theta_{\beta}). In the tluid case, the momentum is conserved and $\,\theta_{\beta}={\rm O}(\Delta t)\,$ and \,\widetilde{\Lambda}_{0k}^{\beta}\,\big{(}\widetilde{\sigma}_{k}^{2}\,-{1\over 6}\big{)}\,\partial_{\beta}\partial_{t}\theta_{k}={\rm O}(\Delta t). In a synthetic way,

[TABLE]

We deduce from (68) and (69) the expansion (65) in the thermal case and the first relation of (66) in the fluid case.

$\bullet$  For the momentum equation of the fluid case, $\,1\leq i=\alpha\leq d\,$ and the relation (67) can be rewritten as

[TABLE]

This relation is exactly the second relation of (66). $\square$

Proposition 9.  Nonlinear expansion at fourth order

When the time step $\,\Delta t\,$ has an infinitesimal value, the expansion of the microscopic variables takes the form

[TABLE]

The vectors $\,\Psi_{1}(W)\,$ and $\,\Psi_{2}(W)\,$ have been evaluated at the relations (39) and (53). We have

[TABLE]

A system of fourth order partial differential equations is asymptotically satisfied by the lattice Boltzmann scheme:

[TABLE]

The vectors $\,\Gamma_{1}(W)$, $\,\Gamma_{2}(W)\,$ and $\,\Gamma_{3}(W)\,$ have been precised at the relations (36), (41) and (55) respectively. We have now

[TABLE]

$\bullet$  Proof of Proposition 9.

We first establish the relations (70) and (71). We start by the expansion of the scheme (61) at order 3. For the second component of nonconserved moments, we have

\left\{\begin{array}[]{l}Y+\Delta t\,\partial_{t}Y+{1\over 2}\,\Delta t^{2}\,\partial_{t}^{2}Y+{1\over 6}\,\Delta t^{3}\,\partial_{t}^{3}Y+{\rm O}(\Delta t^{4})=Y^{*}-\Delta t\,(C\,W+D\,Y^{*})\\ \qquad\qquad+{1\over 2}\,\Delta t^{2}\,(C_{2}\,W+D_{2}\,Y^{*})-{1\over 6}\,\Delta t^{3}\,(C_{3}\,W+D_{3}\,Y^{*})+{\rm O}(\Delta t^{4})\,.\end{array}\right.

The relaxation process can be written $\,Y-Y^{*}\equiv S\,(Y-\Phi(W))$. We deduce

[TABLE]

We have to precise the three first time derivatives of the non-conserved moments. The relation (70) can also be written at second order accuracy:

[TABLE]

We differentiate this previous relation relative to time:

\,\partial_{t}Y={\rm d}\Phi.\,\partial_{t}W+\big{(}\Sigma+{{1}\over{2}}\,{\rm I}\big{)}\,\,\Delta t\,\,{\rm d}\Psi_{1}.\,\partial_{t}W+\big{(}\Sigma+{{1}\over{2}}\,{\rm I}\big{)}\,\,\Delta t^{2}\,\,{\rm d}\Psi_{2}.\,\partial_{t}W+{\rm O}(\Delta t^{3})

=-{\rm d}\Phi.\,\big{(}\Gamma_{1}(W)+\Delta t\,\Gamma_{2}(W)+\Delta t^{2}\,\Gamma_{3}(W)+{\rm O}(\Delta t^{3})\big{)}

-\Delta t\,\big{(}\Sigma+{{1}\over{2}}\,{\rm I}\big{)}\,{\rm d}\Psi_{1}.\,\big{(}\Gamma_{1}+\Delta t\,\Gamma_{2}+{\rm O}(\Delta t^{2})\big{)}+\,\Delta t^{2}\,\big{(}\Sigma+{{1}\over{2}}\,{\rm I}\big{)}\,{\rm d}\Psi_{2}.\,\big{(}-\Gamma_{1}+{\rm O}(\Delta t)\big{)}+{\rm O}(\Delta t^{3})

and

[TABLE]

We need the second time derivative $\,\partial_{t}^{2}Y\,$ and third order time derivative $\,\partial_{t}^{3}Y\,$ in the right hand side of (74). We differentiate the relation (75) relative to time. Using the notation (27), we deduce

\,\partial_{t}^{2}Y=-{\rm d}\gamma_{1}.(-\Gamma_{1}-\Delta t\,\Gamma_{2})-\Delta t\,\big{[}-{\rm d}\gamma_{2}.\Gamma_{1}-\big{(}\Sigma+{{1}\over{2}}\,{\rm I}\big{)}\,\,\partial^{2}\Phi(W).\Gamma_{1}(W)\big{]}+{\rm O}(\Delta t^{2})

and

[TABLE]

Finally,

[TABLE]

We precise now the expansions of $\,C\,W+D\,Y^{*}$, $\,C_{2}\,W+D_{2}\,Y^{*}\,$ and $\,C_{3}\,W+D_{3}\,Y^{*}$. We have

C\,W+D\,Y^{*}=C\,W+D\,\big{[}\Phi(W)+\big{(}\Sigma-{{1}\over{2}}\,{\rm I}\big{)}\,\Delta t\,\Psi_{1}+\big{(}\Sigma-{{1}\over{2}}\,{\rm I}\big{)}\,\Delta t^{2}\,\Psi_{2}\big{]}+{\rm O}(\Delta t^{2})

and

[TABLE]

For the second order term,

C_{2}\,W+D_{2}\,Y^{*}=C_{2}\,W+D_{2}\,\big{[}\Phi(W)+\Delta t\,\big{(}\Sigma-{{1}\over{2}}\,{\rm I}\big{)}\,\Psi_{1}\big{]}+{\rm O}(\Delta t^{2})

and due to (60),

[TABLE]

We have finally for the third order terms

$C_{3}\,W+D_{3}\,Y^{*}=C_{3}\,W+D_{3}\,\Phi(W)+{\rm O}(\Delta t)$

$\,\,\,=(C_{2}\,A+D_{2}\,C)\,W+(C_{2}\,B+D_{2}\,D)\,\Phi(W)+{\rm O}(\Delta t)$

$\,\,\,=C_{2}\,\Gamma_{1}+D_{2}\,(\gamma_{1}-\Psi_{1})+{\rm O}(\Delta t)$

$\,\,\,=(C\,A+D\,C)\,\Gamma_{1}+(C\,B+D^{2})\,\gamma_{1}-D_{2}\,\Psi_{1}+{\rm O}(\Delta t)$

$\,\,\,=C\,{\rm d}\gamma_{1}.\,\Gamma_{1}+D\,{\rm d}(\gamma_{1}-\Psi_{1}).\Gamma_{1}-D_{2}\,\Psi_{1}+{\rm O}(\Delta t)$

$\,\,\,=C\,{\rm d}\gamma_{1}.\,\Gamma_{1}+D\,{\rm d}\gamma_{1}.\Gamma_{1}-D\,{\rm d}\Psi_{1}.\Gamma_{1}-D_{2}\,\Psi_{1}+{\rm O}(\Delta t)$

and

[TABLE]

Then due to (74), (75), (76), (77), (78), (79) and (80), we have

$S\,(Y-\Phi(W))=-\Delta t\,(C\,W+D\,Y^{*}+\partial_{t}Y)+{1\over 2}\,\Delta t^{2}\,(C_{2}\,W+D_{2}\,Y^{*}-\partial_{t}^{2}Y)$

$-{1\over 6}\,\Delta t^{3}\,(C_{3}\,W+D_{3}\,Y^{*}+\partial_{t}^{3}Y)+{\rm O}(\Delta t^{4})$

=-\Delta t\,\big{(}\gamma_{1}-\Psi_{1}+\Delta t\,\,D\,\big{(}\Sigma-{{1}\over{2}}\,{\rm I}\big{)}\,\Psi_{1}+\,\Delta t^{2}\,\,D\,\big{(}\Sigma-{{1}\over{2}}\,{\rm I}\big{)}\,\Psi_{2}\big{)}

-\Delta t\,\big{[}-{\rm d}\Phi.\Gamma_{1}-\Delta t\,\big{(}{\rm d}\Phi.\Gamma_{2}+\big{(}\Sigma+{{1}\over{2}}\,{\rm I}\big{)}\,\,{\rm d}\Psi_{1}.\Gamma_{1}\big{)}

-\Delta t^{2}\,\big{(}{\rm d}\Phi.\Gamma_{3}+\big{(}\Sigma+{{1}\over{2}}\,{\rm I}\big{)}\,\,{\rm d}\Psi_{1}.\Gamma_{2}+\big{(}\Sigma+{{1}\over{2}}\,{\rm I}\big{)}\,\,{\rm d}\Psi_{2}.\Gamma_{1}\big{)}\big{]}

+{1\over 2}\,\Delta t^{2}\,\big{(}{\rm d}(\gamma_{1}-\Psi_{1}).\Gamma_{1}-D\,\Psi_{1}+\Delta t\,D_{2}\,\big{(}\Sigma-{{1}\over{2}}\,{\rm I}\big{)}\,\Psi_{1}\big{)}

-{1\over 2}\,\Delta t^{2}\,\big{(}{\rm d}\gamma_{1}.\,\Gamma_{1}+\,\Delta t\,\big{[}{\rm d}\gamma_{1}.\Gamma_{2}+{\rm d}\gamma_{2}.\Gamma_{1}+\big{(}\Sigma+{{1}\over{2}}\,{\rm I}\big{)}\,\,\partial^{2}\Psi_{1}.\Gamma_{1}\big{]}\big{)}

-{1\over 6}\,\Delta t^{3}\,\big{(}\partial^{2}(\gamma_{1}-\Psi_{1}).\Gamma_{1}-D\,{\rm d}\Psi_{1}.\Gamma_{1}-D_{2}\,\Psi_{1}\big{)}-{1\over 6}\,\Delta t^{3}\,\big{(}-\partial^{2}\gamma_{1}.\Gamma_{1}\big{)}+{\rm O}(\Delta t^{4})

=\Delta t\,\Psi_{1}+\Delta t^{2}\,\Psi_{2}+\Delta t^{3}\,\big{[}-D\,(\Sigma-{1\over 2}\,{\rm I})\,\Psi_{2}+{\rm d}\Phi.\Gamma_{3}+(\Sigma+{1\over 2}\,{\rm I})\,{\rm d}\Psi_{1}.\Gamma_{2}+(\Sigma+{1\over 2}\,{\rm I})\,{\rm d}\Psi_{2}.\Gamma_{1}

$+{1\over 2}\,D_{2}\,(\Sigma-{1\over 2}\,{\rm I})\,\Psi_{1}-{1\over 2}\,{\rm d}\gamma_{1}.\Gamma_{2}-{1\over 2}\,{\rm d}\gamma_{2}.\Gamma_{1}-{1\over 2}\,(\Sigma+{1\over 2}\,{\rm I})\,\partial^{2}\Psi_{1}.\Gamma_{1}+{1\over 6}\,\partial^{2}\Psi_{1}.\Gamma_{1}+{1\over 6}\,D\,{\rm d}\Psi_{1}.\Gamma_{1}$

+{1\over 6}\,D_{2}\,\Psi_{1}\big{]}+{\rm O}(\Delta t^{4})

=\Delta t\,\Psi_{1}+\Delta t^{2}\,\Psi_{2}+\Delta t^{3}\,\big{[}-D\,\Sigma\,\Psi_{2}+{1\over 2}\,D\,(\Sigma\,{\rm d}\Psi_{1}.\Gamma_{1}+\gamma_{2}-D\,\Sigma\,\Psi_{1})+{\rm d}\Phi.\Gamma_{3}+\Sigma\,{\rm d}\Psi_{1}.\Gamma_{2}

$+{1\over 2}\,{\rm d}\Psi_{1}.\Gamma_{2}+\Sigma\,{\rm d}\Psi_{2}.\Gamma_{1}+{1\over 2}\,(\Sigma\,\,\partial^{2}\Psi_{1}.\Gamma_{1}+{\rm d}\gamma_{2}.\Gamma_{1}-D\,\Sigma\,{\rm d}\Psi_{1}.\Gamma_{1})+{1\over 2}\,(C\,B+D^{2})\,\Sigma\,\Psi_{1}$

-{1\over 4}\,D_{2}\,\Psi_{1}-{1\over 2}\,{\rm d}\gamma_{1}.\Gamma_{2}-{1\over 2}\,{\rm d}\gamma_{2}.\Gamma_{1}-{1\over 2}\,\Sigma\,\,\partial^{2}\Psi_{1}.\Gamma_{1}-{1\over 12}\,\partial^{2}\Psi_{1}.\Gamma_{1}+{1\over 6}\,D\,{\rm d}\Psi_{1}.\Gamma_{1}+{1\over 6}\,D_{2}\,\Psi_{1}\big{]}

$+{\rm O}(\Delta t^{4})$

=\Delta t\,\Psi_{1}+\Delta t^{2}\,\Psi_{2}+\Delta t^{3}\,\big{[}{\rm d}\Phi.\Gamma_{3}-D\,\Sigma\,\Psi_{2}+{1\over 2}\,D\,\gamma_{2}+\Sigma\,{\rm d}\Psi_{1}.\Gamma_{2}+{1\over 2}\,{\rm d}\Psi_{1}.\Gamma_{2}+\Sigma\,\,{\rm d}\Psi_{2}.\Gamma_{1}

+\,{1\over 2}\,C\,B\,\Sigma\,\,\Psi_{1}-{1\over 12}\,D_{2}\,\Psi_{1}-{1\over 2}\,{\rm d}\gamma_{1}.\Gamma_{2}-{1\over 12}\,\partial^{2}\Psi_{1}.\Gamma_{1}+{1\over 6}\,D\,{\rm d}\Psi_{1}.\Gamma_{1}\big{]}+{\rm O}(\Delta t^{4})

=\Delta t\,\Psi_{1}+\Delta t^{2}\,\Psi_{2}+\Delta t^{3}\,\big{[}{\rm d}\Phi.\Gamma_{3}-D\,\Sigma\,\Psi_{2}+{1\over 2}\,D\,\gamma_{2}+\Sigma\,{\rm d}\Psi_{1}.\Gamma_{2}+{1\over 2}\,(-C\,\Gamma_{2}-D\,\gamma_{2}+{\rm d}\gamma_{1}.\Gamma_{2})

+\Sigma\,{\rm d}\Psi_{2}.\Gamma_{1}+\,{1\over 2}\,C\,B\,\Sigma\,\,\Psi_{1}-{1\over 12}\,D_{2}\,\Psi_{1}-{1\over 2}\,{\rm d}\gamma_{1}.\Gamma_{2}-{1\over 12}\,\partial^{2}\Psi_{1}.\Gamma_{1}+{1\over 6}\,D\,{\rm d}\Psi_{1}.\Gamma_{1}\big{]}+{\rm O}(\Delta t^{4})

=\Delta t\,\Psi_{1}+\Delta t^{2}\,\Psi_{2}+\Delta t^{3}\,\big{[}{\rm d}\Phi.\Gamma_{3}-D\,\Sigma\,\Psi_{2}+\Sigma\,{\rm d}\Psi_{1}.\Gamma_{2}+\Sigma\,{\rm d}\Psi_{2}.\Gamma_{1}-{1\over 12}\,D_{2}\,\Psi_{1}-{1\over 12}\,\partial^{2}\Psi_{1}.\Gamma_{1}

+{1\over 6}\,D\,{\rm d}\Psi_{1}.\Gamma_{1}\big{]}+{\rm O}(\Delta t^{4})\,.

The relations (70) and (71) are established.

$\bullet$  Fourth order partial differential equations

We consider now the Taylor expansion of the relation (21) at fourth order:

[TABLE]

For the first components (conserved moments):

[TABLE]

We simplify the constant term $\,W\,$ and divide by $\,\Delta t$. We deduce

[TABLE]

We explicit the partial derivatives $\,\partial_{t}^{2}W$, $\,\partial_{t}^{3}W\,$ and $\,\partial_{t}^{4}W\,$ at orders two, one and zero respectively. We start with the relation

[TABLE]

where $\,\Gamma_{1}$, $\,\Gamma_{2}\,$ and $\,\Gamma_{3}\,$ have been evaluated previously. We have

\partial_{t}^{2}W=\partial_{t}\big{(}-\Gamma_{1}-\Delta t\,\Gamma_{2}-\Delta t^{2}\,\Gamma_{3}+{\rm O}(\Delta t^{3})\big{)}

\,\,={\rm d}\big{(}-\Gamma_{1}-\Delta t\,\Gamma_{2}-\Delta t^{2}\,\Gamma_{3}+{\rm O}(\Delta t^{3})\big{)}.(\partial_{t}W)

\,\,={\rm d}\big{(}-\Gamma_{1}-\Delta t\,\Gamma_{2}-\Delta t^{2}\,\Gamma_{3}+{\rm O}(\Delta t^{2})\big{)}\,.\,\big{(}-\Gamma_{1}-\Delta t\,\Gamma_{2}-\Delta t^{2}\,\Gamma_{3}\big{)}+{\rm O}(\Delta t^{3})

and

[TABLE]

Then

$\,\partial_{t}^{3}W=\partial_{t}(\partial_{t}^{2}W)$

\,\,\,=\partial_{t}\big{(}A\,\Gamma_{1}+B\,{\rm d}\Phi.\Gamma_{1}+\Delta t\,\partial_{t}({\rm d}\Gamma_{2}.\,\Gamma_{1}+{\rm d}\Gamma_{1}.\,\Gamma_{2})\big{)}+{\rm O}(\Delta t^{2})

\,\,\,=(A\,{\rm d}\Gamma_{1}+B\,{\rm d}({\rm d}\Phi.\Gamma_{1})).\,\big{(}-\Gamma_{1}-\Delta t\,\Gamma_{2}+{\rm O}(\Delta t^{2})\big{)}

$+\Delta t\,\partial_{t}\,(B\,\Sigma\,{\rm d}\Psi_{1}.\,\Gamma_{1}+A\,\Gamma_{2}+B\,{\rm d}\Phi.\,\Gamma_{2})+{\rm O}(\Delta t^{2})$

\,\,\,=\partial_{t}\big{(}{\rm d}\Gamma_{1}.\,\Gamma_{1}+\Delta t\,({\rm d}\Gamma_{1}.\,\Gamma_{2}+{\rm d}\Gamma_{2}.\,\Gamma_{1})\big{)}+{\rm O}(\Delta t^{2})

\,\,\,={\rm d}\big{(}{\rm d}\Gamma_{1}.\,\Gamma_{1}+\Delta t\,({\rm d}\Gamma_{1}.\,\Gamma_{2}+{\rm d}\Gamma_{2}.\,\Gamma_{1})\big{)}.\,\partial_{t}W+{\rm O}(\Delta t^{2})

\,\,\,={\rm d}\big{(}{\rm d}\Gamma_{1}.\,\Gamma_{1}+\Delta t\,({\rm d}\Gamma_{1}.\,\Gamma_{2}+{\rm d}\Gamma_{2}.\,\Gamma_{1})\big{)}.\big{(}-\Gamma_{1}-\Delta t\,\Gamma_{2})+{\rm O}(\Delta t^{2})

\,\,\,=-{\rm d}\big{(}{\rm d}\Gamma_{1}.\,\Gamma_{1}\big{)}.\Gamma_{1}-\Delta t\,\big{[}{\rm d}(A\,\Gamma_{2}+B\,\gamma_{2}+B\,\Sigma\,{\rm d}\Psi_{1}.\,\Gamma_{1}).\,\Gamma_{1}+{\rm d}(A\,\Gamma_{1}+B\,\gamma_{1}).\,\Gamma_{2}\big{]}+{\rm O}(\Delta t^{2})

\,\,\,=-\partial^{2}\Gamma_{1}.\,\Gamma_{1}-\Delta t\,\big{[}A\,{\rm d}\Gamma_{2}.\,\Gamma_{1}+B\,{\rm d}\gamma_{2}.\,\Gamma_{1}+B\,\Sigma\,\partial^{2}\Psi_{1}.\,\Gamma_{1}+A\,{\rm d}\Gamma_{1}.\,\Gamma_{2}+B\,{\rm d}\gamma_{1}.\,\Gamma_{2}\big{]}+{\rm O}(\Delta t^{2})

and we have at first order for this third order derivative

[TABLE]

Last but not least, we have

\partial_{t}^{4}W=\partial_{t}\big{(}-\partial^{2}\Gamma_{1}.\,\Gamma_{1}+{\rm O}(\Delta t)\big{)}

\,\,=-\partial_{t}\big{(}A\,{\rm d}\Gamma_{1}.\,\Gamma_{1}+B\,{\rm d}\gamma_{1}.\,\Gamma_{1})+{\rm O}(\Delta t)\big{)}

\,\,={\rm d}\big{(}A\,{\rm d}\Gamma_{1}.\,\Gamma_{1}+B\,{\rm d}\gamma_{1}.\,\Gamma_{1}\big{)}.\,(\Gamma_{1}+{\rm O}(\Delta t))\,

and finally

[TABLE]

For the determination of the coefficient $\,\Gamma_{4}$, we precise the expansions of the quantities $\,A_{j}\,W+B_{j}\,Y^{*}\,$ for $\,j=1\,$ to $\,j=4$. We use the previous evaluations done in the proof of the previous propositions:

A\,W+B\,Y^{*}=A\,W+B\,\big{[}\Phi(W)+\big{(}\Sigma-{{1}\over{2}}\,{\rm I}\big{)}\,\big{(}\Delta t\,\Psi_{1}+\Delta t^{2}\,\Psi_{2}+\Delta t^{3}\,\Psi_{3}\big{)}\big{]}+{\rm O}(\Delta t^{4})

\,\,\,\,=A\,W+B\,\Phi+\Delta t\,B\,\big{(}\Sigma-{{1}\over{2}}\,{\rm I}\big{)}\,\Psi_{1}+\Delta t^{2}\,B\,\big{(}\Sigma-{{1}\over{2}}\,{\rm I}\big{)}\,\Psi_{2}+\Delta t^{3}\,B\,\big{(}\Sigma-{{1}\over{2}}\,{\rm I}\big{)}\,\Psi_{3}+{\rm O}(\Delta t^{4})

A_{2}\,W+B_{2}\,Y^{*}=A_{2}\,W+B_{2}\,\big{[}\Phi(W)+\big{(}\Sigma-{{1}\over{2}}\,{\rm I}\big{)}\,(\Delta t\,\Psi_{1}+\Delta t^{2}\,\Psi_{2})\big{]}+{\rm O}(\Delta t^{3})

=(A^{2}+B\,C)\,W+(A\,B+B\,D)\,\Phi+\Delta t\,B_{2}\,\big{(}\Sigma-{{1}\over{2}}\,{\rm I}\big{)}\,\Psi_{1}+\Delta t^{2}\,B_{2}\,\big{(}\Sigma-{{1}\over{2}}\,{\rm I}\big{)}\,\Psi_{2}+{\rm O}(\Delta t^{3})

={\rm d}\Gamma_{1}.\Gamma_{1}-B\,\Psi_{1}+\Delta t\,B_{2}\,\big{(}\Sigma-{{1}\over{2}}\,{\rm I}\big{)}\,\Psi_{1}+\Delta t^{2}\,B_{2}\,\big{(}\Sigma-{{1}\over{2}}\,{\rm I}\big{)}\,\Psi_{2}+{\rm O}(\Delta t^{3})

A_{3}\,W+B_{3}\,Y^{*}=A_{3}\,W+B_{3}\,[\Phi(W)+\Delta t\,\big{(}\Sigma-{{1}\over{2}}\,{\rm I}\big{)}\,\Psi_{1}]+{\rm O}(\Delta t^{2})

\,\,\,=(A_{2}\,A+B_{2}\,C)\,W+(A_{2}\,B+B_{2}\,D)\,\Phi+\Delta t\,B_{3}\,\big{(}\Sigma-{{1}\over{2}}\,{\rm I}\big{)}\,\Psi_{1}+{\rm O}(\Delta t^{2})

\,\,\,=\partial^{2}\Gamma_{1}.\Gamma_{1}-B\,{\rm d}\Psi_{1}.\Gamma_{1}-B_{2}\,\Psi_{1}+\Delta t\,B_{3}\,\big{(}\Sigma-{{1}\over{2}}\,{\rm I}\big{)}\,\Psi_{1}+{\rm O}(\Delta t^{2})

$A_{4}\,W+B_{4}\,Y^{*}=A_{4}\,W+B_{4}\,\Phi(W)+{\rm O}(\Delta t)$

$\,\,\,=(A_{3}\,A+B_{3}\,C)\,W+(A_{3}\,B+B_{3}\,D)\,\Phi+{\rm O}(\Delta t)$

$\,\,\,=A_{3}\,\Gamma_{1}+B_{3}\,(\gamma_{1}-\Psi_{1})+{\rm O}(\Delta t)$

$\,\,\,=(A_{2}\,A+B_{2}\,C)\,\Gamma_{1}+(A_{2}\,B+B_{2}\,D)\,\gamma_{1}-B_{3}\,\Psi_{1}+{\rm O}(\Delta t)$

$\,\,\,=A_{2}\,{\rm d}\Gamma_{1}.\,\Gamma_{1}+B_{2}\,{\rm d}(\gamma_{1}-\Psi_{1}).\,\Gamma_{1}-B_{3}\,\Psi_{1}+{\rm O}(\Delta t)$

$\,\,\,=(A^{2}+B\,C)\,{\rm d}\Gamma_{1}.\,\Gamma_{1}+(A\,B+B\,D)\,{\rm d}\gamma_{1}.\,\Gamma_{1}-B_{2}\,{\rm d}\Psi_{1}.\,\Gamma_{1}-B_{3}\,\Psi_{1}+{\rm O}(\Delta t)$

$\,\,\,=A\,\partial^{2}\Gamma_{1}.\,\Gamma_{1}+B\,\partial^{2}(\gamma_{1}-\Psi_{1}).\,\Gamma_{1}-B_{2}\,{\rm d}\Psi_{1}.\,\Gamma_{1}-B_{3}\,\Psi_{1}+{\rm O}(\Delta t)$

$\,\,\,=A\,\partial^{2}\Gamma_{1}.\,\Gamma_{1}+B\,\partial^{2}\gamma_{1}.\,\Gamma_{1}-B\,\partial^{2}\Psi_{1}.\,\Gamma_{1}-B_{2}\,{\rm d}\Psi_{1}.\,\Gamma_{1}-B_{3}\,\Psi_{1}+{\rm O}(\Delta t)\,.$

We deduce now from (81), (82), (83), (84) and the previous expressions:

$\partial_{t}W=-(A\,W+B\,Y^{*})+{1\over 2}\,\Delta t\,(A_{2}\,W+B_{2}\,Y^{*}-\partial_{t}^{2}W)-{1\over 6}\,\Delta t^{2}\,(A_{3}\,W+B_{3}\,Y^{*}+\partial_{t}^{3}W)$

$+{1\over 24}\,\Delta t^{3}\,(A_{4}\,W+B_{4}\,Y^{*}-\partial_{t}^{4}W)+{\rm O}(\Delta t^{4})$

\,=-\big{(}A\,W+B\,\Phi+\Delta t\,B\,\big{(}\Sigma-{{1}\over{2}}\,{\rm I}\big{)}\,\Psi_{1}+\Delta t^{2}\,B\,\big{(}\Sigma-{{1}\over{2}}\,{\rm I}\big{)}\,\Psi_{2}+\Delta t^{3}\,B\,\big{(}\Sigma-{{1}\over{2}}\,{\rm I}\big{)}\,\Psi_{3}\big{)}

\qquad+{1\over 2}\,\Delta t\,\big{[}{\rm d}\Gamma_{1}.\Gamma_{1}-B\,\Psi_{1}+\Delta t\,B_{2}\,\big{(}\Sigma-{{1}\over{2}}\,{\rm I}\big{)}\,\Psi_{1}+\Delta t^{2}\,B_{2}\,\big{(}\Sigma-{{1}\over{2}}\,{\rm I}\big{)}\,\Psi_{2}

\qquad\quad-\big{(}{\rm d}\Gamma_{1}.\,\Gamma_{1}+\Delta t\,({\rm d}\Gamma_{1}.\,\Gamma_{2}+{\rm d}\Gamma_{2}.\,\Gamma_{1})+\Delta t^{2}\,({\rm d}\Gamma_{1}.\,\Gamma_{3}+{\rm d}\Gamma_{2}.\,\Gamma_{2}+{\rm d}\Gamma_{3}.\,\Gamma_{1})\big{)}\big{]}

\qquad-{1\over 6}\,\Delta t^{2}\,\big{[}\partial^{2}\Gamma_{1}.\Gamma_{1}-B\,{\rm d}\Psi_{1}.\Gamma_{1}-B_{2}\,\Psi_{1}+\Delta t\,B_{3}\,\big{(}\Sigma-{{1}\over{2}}\,{\rm I}\big{)}\,\Psi_{1}

\qquad\quad-\partial^{2}\Gamma_{1}.\,\Gamma_{1}-\Delta t\,\big{(}A\,({\rm d}\Gamma_{2}.\,\Gamma_{1}+{\rm d}\Gamma_{1}.\,\Gamma_{2})+B\,({\rm d}\gamma_{2}.\,\Gamma_{1}+{\rm d}\gamma_{1}.\,\Gamma_{2}+\Sigma\,\partial^{2}\Psi_{1}.\,\Gamma_{1})\big{)}\big{]}

\qquad+{1\over 24}\,\Delta t^{3}\,\big{[}A\,\partial^{2}\Gamma_{1}.\,\Gamma_{1}+B\,\partial^{2}\gamma_{1}.\,\Gamma_{1}-B\,\partial^{2}\Psi_{1}.\,\Gamma_{1}-B_{2}\,{\rm d}\Psi_{1}.\,\Gamma_{1}-B_{3}\,\Psi_{1}

\qquad\quad-\big{(}A\,\partial^{2}\Gamma_{1}.\,\Gamma_{1}+B\,\partial^{2}\gamma_{1}.\,\Gamma_{1}\big{)}\big{]}+{\rm O}(\Delta t^{4})

\,=-\Gamma_{1}-\Delta t\,\Gamma_{2}-\Delta t^{2}\,\Gamma_{3}+\Delta t^{3}\,\big{[}-B\,\big{(}\Sigma-{{1}\over{2}}\,{\rm I}\big{)}\,\Psi_{3}+{1\over 2}\,B_{2}\,\big{(}\Sigma-{{1}\over{2}}\,{\rm I}\big{)}\,\Psi_{2}

\qquad-{1\over 2}\,({\rm d}\Gamma_{1}.\,\Gamma_{3}+{\rm d}\Gamma_{2}.\,\Gamma_{2}+{\rm d}\Gamma_{3}.\,\Gamma_{1})+{1\over 6}\,\big{(}-B_{3}\,\big{(}\Sigma-{{1}\over{2}}\,{\rm I}\big{)}\,\Psi_{1}+A\,({\rm d}\Gamma_{2}.\,\Gamma_{1}+{\rm d}\Gamma_{1}.\,\Gamma_{2})

\qquad+B\,({\rm d}\gamma_{2}.\,\Gamma_{1}+{\rm d}\gamma_{1}.\,\Gamma_{2}+\Sigma\,\partial^{2}\Psi_{1}.\,\Gamma_{1})\big{)}+{1\over 24}\,\big{(}A\,\partial^{2}\Gamma_{1}.\,\Gamma_{1}+B\,\partial^{2}\gamma_{1}.\,\Gamma_{1}

\qquad-B\,\partial^{2}\Psi_{1}.\,\Gamma_{1}-B_{2}\,{\rm d}\Psi_{1}.\,\Gamma_{1}-B_{3}\,\Psi_{1}-A\,\partial^{2}\Gamma_{1}.\,\Gamma_{1}-B\,\partial^{2}\gamma_{1}.\,\Gamma_{1}\big{)}\big{]}+{\rm O}(\Delta t^{4})

\,=-\Gamma_{1}-\Delta t\,\Gamma_{2}-\Delta t^{2}\,\Gamma_{3}-\Delta t^{3}\,\big{[}B\,\Sigma\,\Psi_{3}-{1\over 2}\,B\,\big{(}\gamma_{3}+\Sigma\,{\rm d}\Psi_{1}.\Gamma_{2}+\Sigma\,{\rm d}\Psi_{2}.\Gamma_{1}-D\,\Sigma\,\Psi_{2}

\qquad+{1\over 6}\,D\,{\rm d}\Psi_{1}.\Gamma_{1}-{1\over 12}\,D_{2}\,\Psi_{1}-{1\over 12}\,\partial^{2}\Psi_{1}.\Gamma_{1}\big{)}-{1\over 2}\,(A\,B+B\,D)\,\Sigma\,\Psi_{2}+{1\over 4}\,B_{2}\,\Psi_{2}

\qquad+{1\over 2}\,\big{(}A\,\Gamma_{3}+B\,\gamma_{3}+B\,\Sigma\,{\rm d}\Psi_{1}.\Gamma_{1}+{\rm d}\Gamma_{3}.\,\Gamma_{1}\big{)}+{1\over 6}\,B_{3}\,\Sigma\,\Psi_{1}-{1\over 12}\,B_{3}\,\Psi_{1}-{1\over 6}\,A\,B\,\Sigma\,{\rm d}\Psi_{1}.\Gamma_{1}

\qquad-{1\over 6}\,A\,(A\,\Gamma_{2}+B\,\gamma_{2})-{1\over 6}\,B\,\big{(}{\rm d}\gamma_{2}.\,\Gamma_{1}+{\rm d}\gamma_{1}.\,\Gamma_{2}+\Sigma\,\partial^{2}\Psi_{1}.\,\Gamma_{1}\big{)}+{1\over 24}\,B\,\partial^{2}\Psi_{1}.\,\Gamma_{1}+{1\over 24}\,B_{2}\,{\rm d}\Psi_{1}.\,\Gamma_{1}

\qquad+{1\over 24}\,B_{3}\,\Psi_{1}\big{]}+{\rm O}(\Delta t^{4})

\,=-\Gamma_{1}-\Delta t\,\Gamma_{2}-\Delta t^{2}\,\Gamma_{3}-\Delta t^{3}\,\big{[}B\,\Sigma\,\Psi_{3}-{1\over 2}\,B\,\Sigma\,{\rm d}\Psi_{2}.\Gamma_{1}-{1\over 12}\,B\,D\,{\rm d}\Psi_{1}.\Gamma_{1}+{1\over 24}\,B\,D_{2}\,\Psi_{1}

\qquad+{1\over 24}\,B\,\partial^{2}\Psi_{1}.\Gamma_{1}-{1\over 2}\,A\,B\,\Sigma\,\Psi_{2}+{1\over 4}\,B_{2}\,\Psi_{2}+{1\over 2}\,A\,\big{(}B\,\Sigma\,\Psi_{2}-{1\over 6}\,B\,{\rm d}\Psi_{1}.\Gamma_{1}+{1\over 12}\,B_{2}\,\Psi_{1}\big{)}

\qquad+{1\over 2}\,\big{(}B\,\Sigma\,{\rm d}\Psi_{2}.\Gamma_{1}-{1\over 6}\,B\,\partial^{2}\Psi_{1}.\Gamma_{1}+{1\over 12}\,B_{2}\,{\rm d}\Psi_{1}.\Gamma_{1}\big{)}+{1\over 6}\,B_{3}\,\Sigma\,\Psi_{1}-{1\over 24}\,B_{3}\,\Psi_{1}-{1\over 6}\,A\,B\,\Sigma\,{\rm d}\Psi_{1}.\Gamma_{1}

\qquad-{1\over 6}\,A^{2}\,\Gamma_{2}-{1\over 6}\,A\,B\,\gamma_{2}-{1\over 6}\,B\,\big{(}{\rm d}\gamma_{2}.\,\Gamma_{1}+{\rm d}\gamma_{1}.\,\Gamma_{2}+\Sigma\,\partial^{2}\Psi_{1}.\,\Gamma_{1}\big{)}+{1\over 24}\,B\,\partial^{2}\Psi_{1}.\,\Gamma_{1}

\qquad+{1\over 24}\,B_{2}\,{\rm d}\Psi_{1}.\,\Gamma_{1}\big{]}+{\rm O}(\Delta t^{4})

\,=-\Gamma_{1}-\Delta t\,\Gamma_{2}-\Delta t^{2}\,\Gamma_{3}-\Delta t^{3}\,\big{[}B\,\Sigma\,\Psi_{3}-{1\over 12}\,B_{2}\,{\rm d}\Psi_{1}.\,\Gamma_{1}+{1\over 24}\,B\,D_{2}\,\Psi_{1}+{1\over 4}\,B_{2}\,\Psi_{2}+{1\over 24}\,A\,B_{2}\,\Psi_{1}

$\qquad+{1\over 6}\,B_{3}\,\Sigma\,\Psi_{1}-{1\over 24}\,B_{3}\,\Psi_{1}-{1\over 6}\,A\,B\,\Sigma\,{\rm d}\Psi_{1}.\Gamma_{1}-{1\over 6}\,A^{2}\,B\,\Sigma\,\Psi_{1}-{1\over 6}\,A\,B\,\gamma_{2}$

\qquad-{1\over 6}\,B\,\big{(}{\rm d}\gamma_{2}.\,\Gamma_{1}+{\rm d}\gamma_{1}.\,\Gamma_{2}+\Sigma\,\partial^{2}\Psi_{1}.\,\Gamma_{1}\big{)}+{1\over 12}\,B_{2}\,{\rm d}\Psi_{1}.\,\Gamma_{1}\big{]}+{\rm O}(\Delta t^{4})

\,=-\Gamma_{1}-\Delta t\,\Gamma_{2}-\Delta t^{2}\,\Gamma_{3}-\Delta t^{3}\,\big{[}B\,\Sigma\,\Psi_{3}+{1\over 4}\,B_{2}\,\Psi_{2}+{1\over 6}\,B_{3}\,\Sigma\,\Psi_{1}-{1\over 6}\,A\,B\,\Sigma\,{\rm d}\Psi_{1}.\Gamma_{1}

\qquad-{1\over 6}\,A^{2}\,B\,\Sigma\,\Psi_{1}-{1\over 6}\,A\,B\,\big{(}\Psi_{2}-\Sigma\,{\rm d}\Psi_{1}.\,\Gamma_{1}+D\,\Sigma\,\Psi_{1}\big{)}-{1\over 6}\,B\,\big{(}{\rm d}\gamma_{2}.\,\Gamma_{1}+{\rm d}\gamma_{1}.\,\Gamma_{2}+\Sigma\,\partial^{2}\Psi_{1}.\,\Gamma_{1}\big{)}\big{]}

$\qquad+{\rm O}(\Delta t^{4})$ because $\,A\,B_{2}+B\,D_{2}=B_{3}$

\,=-\Gamma_{1}-\Delta t\,\Gamma_{2}-\Delta t^{2}\,\Gamma_{3}-\Delta t^{3}\,\big{[}B\,\Sigma\,\Psi_{3}+{1\over 4}\,B_{2}\,\Psi_{2}+{1\over 6}\,B\,D_{2}\,\Sigma\,\Psi_{1}-{1\over 6}\,A\,B\,\Psi_{2}

\qquad-{1\over 6}\,B\,\big{(}{\rm d}\gamma_{2}.\,\Gamma_{1}+{\rm d}\gamma_{1}.\,\Gamma_{2}+\Sigma\,\partial^{2}\Psi_{1}.\,\Gamma_{1}\big{)}\big{]}+{\rm O}(\Delta t^{4})

because $\,B_{3}-A\,(A\,B+B\,D)=B\,D_{2}$. The relations (72) and (73) are established and the proposition is proven. $\square$

**7)  Revisiting the “Berlin algorithm” in the linear case **

In [4], we have presented the Berlin explicit algorithm in order to determine, with the help of formal calculus, the equivalent partial differential equations of a linear lattice Boltzmann scheme. We do not enter here into the details of our previous contribution because the linearized version of the present algorithm is much more simple.

**Proposition 10.  Linearized general expansion at fourth order **

We suppose that the equilibrium value of the nonconserved moments is a linear function:

[TABLE]

with a given fixed $\,(q-N)\times N\,$ rectangle matrix $\,E$. Then the nonequilibium moments can be expanded as

[TABLE]

The $\,(q-N)\times N\,$ matrices $\,\beta_{j}\,$ are linear operators of order $\,j$. The equivalent partial equivalent equations are linear and can be written as

[TABLE]

and $\,\alpha_{j}\,$ is a $\,N\times N\,$ differential matrix of order $\,j$. In other terms, $\,\Phi_{j}(W)=\beta_{j}\,W\,$ and $\,\Gamma_{j}(W)=\alpha_{j}\,W\,$ are linear functions of the conserved moments. The coefficent matrices $\,\alpha_{j}\,$ and $\,\beta_{j}\,$ are computed according to a linearized version of (31):

[TABLE]

$\bullet$  Proof of Proposition 10.

First observe that $\,{\rm d}\Phi.\,\xi=E\,\xi\,$ for any test vector $\,\xi$. Due to the relation (36), we have

$\Gamma_{1}=A\,W+B\,\Phi(W)=A\,W+B\,E\,W=(A+B\,E)\,W\,$ and $\,\alpha_{1}=A+B\,E$. Due to (39), we have moreover

\Psi_{1}={\rm d}\Phi(W).\,\Gamma_{1}(W)-\big{(}C\,W+D\,\Phi(W)\big{)}

$\,\,=E\,\alpha_{1}\,W-C\,W-D\,E\,W$

$\,\,=(E\,\alpha_{1}-C-D\,E)\,W$

and the second relation of (85) is established. Due to (41),

$\Gamma_{2}(W)=B\,\Sigma\,\,\Psi_{1}(W)$

$\quad\,=B\,\Sigma\,\,\beta_{1}\,W\,$.

Then $\,\alpha_{2}=B\,\Sigma\,\,\beta_{1}\,$ and the third line of (85) is proven. We have now from (53),

$\Psi_{2}(W)=\Sigma\,{\rm d}\Psi_{1}(W).\Gamma_{1}(W)+{\rm d}\Phi(W).\Gamma_{2}(W)-D\,\Sigma\,\Psi_{1}(W)$

$\quad\,\,=\Sigma\,\beta_{1}\,\alpha_{1}\,W+E\,\alpha_{2}\,W-D\,\Sigma\,\beta_{1}\,W$

\quad\,\,=\big{(}\Sigma\,\beta_{1}\,\alpha_{1}+E\,\alpha_{2}-D\,\Sigma\,\beta_{1}\big{)}\,W

and $\,\beta_{2}=\Sigma\,\beta_{1}\,\alpha_{1}+E\,\alpha_{2}-D\,\Sigma\,\beta_{1}\,$ as suggested in the fourth line of(85). From the relation (55), we have now

$\Gamma_{3}(W)=B\,\Sigma\,\Psi_{2}(W)+{1\over 12}\,B_{2}\,\Psi_{1}(W)-{1\over 6}\,B\,{\rm d}\Psi_{1}(W).\Gamma_{1}(W)$

$\quad\,=B\,\Sigma\,\beta_{2}\,W+{1\over 12}\,B_{2}\,\beta_{1}\,W-{1\over 6}\,B\,\beta_{1}.\alpha_{1}\,W$

\quad\,=\big{(}B\,\Sigma\,\beta_{2}+{1\over 12}\,B_{2}\,\beta_{1}-{1\over 6}\,B\,\beta_{1}.\alpha_{1}\big{)}\,W\,.

Then $\,\alpha_{3}=B\,\Sigma\,\beta_{2}+{1\over 12}\,B_{2}\,\beta_{1}-{1\over 6}\,B\,\beta_{1}.\alpha_{1}\,$ as proposed in the fifth line of (85). Observe now that due to (28), we have

$\partial^{2}\Psi_{1}.\,\Gamma_{1}(W)={\rm d}^{2}\Psi_{1}(W).(\Gamma_{1},\,\Gamma_{1})+{\rm d}\Psi_{1}(W).{\rm d}\Gamma_{1}(W).\,\Gamma_{1}$

$\quad\,={\rm d}\Psi_{1}(W).{\rm d}\Gamma_{1}(W).\,\Gamma_{1}$ because $\,\Psi_{1}\,$ is linear

$\quad\,=\beta_{1}\,\alpha_{1}^{2}\,W\,.$

Then due to (71), we have in this linear case,

$\Psi_{3}(W)={\rm d}\Phi(W).\Gamma_{3}(W)+\Sigma\,\,{\rm d}\Psi_{1}(W).\Gamma_{2}(W)+\Sigma\,\,{\rm d}\Psi_{2}(W).\Gamma_{1}(W)-D\,\Sigma\,\Psi_{2}(W)$

$\qquad\quad+{1\over 6}\,D\,{\rm d}\Psi_{1}(W).\Gamma_{1}(W)-{1\over 12}\,D_{2}\,\Psi_{1}(W)-{1\over 12}\,\partial^{2}\Psi_{1}(W).\Gamma_{1}(W)$

$\quad\,=E\,\alpha_{3}\,W+\Sigma\,\beta_{1}\,\alpha_{2}\,W+\Sigma\,\beta_{2}\,\alpha_{1}\,W-D\,\Sigma\,\beta_{2}\,W+{1\over 6}\,D\,\beta_{1}\,\alpha_{1}\,W$

$\qquad\quad-{1\over 12}\,D_{2}\,\beta_{1}\,W-{1\over 12}\,\beta_{1}\,\alpha_{1}^{2}\,W$

\quad\,=\big{(}E\,\alpha_{3}+\Sigma\,\beta_{1}\,\alpha_{2}+\Sigma\,\beta_{2}\,\alpha_{1}-D\,\Sigma\,\beta_{2}+{1\over 6}\,D\,\beta_{1}\,\alpha_{1}-{1\over 12}\,D_{2}\,\beta_{1}-{1\over 12}\,\beta_{1}\,\alpha_{1}^{2}\big{)}\,W

and $\,\beta_{3}=E\,\alpha_{3}+\Sigma\,\beta_{1}\,\alpha_{2}+\Sigma\,\beta_{2}\,\alpha_{1}-D\,\Sigma\,\beta_{2}+{1\over 6}\,D\,\beta_{1}\,\alpha_{1}-{1\over 12}\,D_{2}\,\beta_{1}-{1\over 12}\,\beta_{1}\,\alpha_{1}^{2}$. We have finally, due to (73),

$\Gamma_{4}(W)=B\,\Sigma\,\Psi_{3}(W)+{1\over 4}\,B_{2}\,\Psi_{2}(W)+{1\over 6}\,B\,D_{2}\,\Sigma\,\Psi_{1}-{1\over 6}\,A\,B\,\Psi_{2}(W)-{1\over 6}\,B\,{\rm d}\gamma_{1}(W).\Gamma_{2}(W)$

$\qquad\quad-{1\over 6}\,B\,{\rm d}\gamma_{2}(W).\Gamma_{1}(W)-{1\over 6}\,B\,\Sigma\,\partial^{2}\Psi_{1}(W).\Gamma_{1}(W)\,$

$\,\,\,=B\,\Sigma\,\beta_{3}\,W+{1\over 4}\,B_{2}\,\beta_{2}\,W+{1\over 6}\,B\,D_{2}\,\Sigma\,\beta_{1}\,W-{1\over 6}\,A\,B\,\beta_{2}\,W-{1\over 6}\,B\,E\,\alpha_{2}\,W$

$\qquad\quad-{1\over 6}\,B\,E\,\alpha_{1}\,W-{1\over 6}\,B\,\Sigma\,\beta_{1}\,\alpha_{1}^{2}\,W\,$

\,\,\,=\big{(}B\,\Sigma\,\beta_{3}+{1\over 4}\,B_{2}\,\beta_{2}+{1\over 6}\,B\,D_{2}\,\Sigma\,\beta_{1}-{1\over 6}\,A\,B\,\beta_{2}-{1\over 6}\,B\,E\,\alpha_{2}\,\alpha_{1}-{1\over 6}\,B\,E\,\alpha_{1}\,\alpha_{2}-{1\over 6}\,B\,\Sigma\,\beta_{1}\,\alpha_{1}^{2}\big{)}\,W

and the proposition is established. $\square$

**8)  Conclusion **

In this contribution, we have extended the Taylor expansion method of a multiple relaxation times lattice Boltzmann scheme up to fourth order accuracy. With this expansion, nonlinear partial differential equations for the conserved variables and various differential expressions for the nonconserved moments emerge naturally. This expansion has been validated with previous works on the subject.

Observe that the main result is implemented with this framework in the “pyLBM”’ software [31]. The next step is the introduction of source terms, in the spirit proposed in [26], and other geometries like triangles, in a way suggested in [24]. The adaptation of the Taylor expansion method to other variants of lattice Boltzmann schemes as recentered schemes [29, 20], regularized lattice Boltzmann [45, 49, 54] or entropic lattice Boltzmann schemes [42] are also into study. After a first step in [48], the link between the single scale Taylor expansion and the multiple scale Chapman-Enskog is also a project for the near future.

**Acknowledgments **

The author thanks Bruce Boghosian, Filipa Caetano, Loïc Gouarin, Benjamin Graille, Pierre Lallemand and Paulo Philippi for impressive disussions. The scientific exchanges with the referee greatly improved the matter of this contribution. All the formal computations have been implemented with SageMath [53].

**References **