Nonlinear fourth order Taylor expansion of lattice Boltzmann schemes

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.

Keywords

Partial differential equations asymptotic analysis

1. Introduction

The lattice Boltzmann schemes in their modern form have been developed with the contributions of d’Humières, Lallemand, Succi [13,37,43,50] 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 [4]. 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. [41] and to the prospective article of Succi [56]. 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 [8] expansion in the spirit of the book of Huang [38]. 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 [9] and by Qian and Zhou [51]. 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 [43], Philippi and Hegele [48], Shan et al. [53] and Otomo et al. [46] 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 [45] and Warming and Hyett [57]. 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,42]. It was extended to third order accuracy for thermal and fluid problems in [17], automated at an arbitrary order in the linear case in [3], 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. [1,55]) 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 [47]. 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. [58].

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 [13] 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 $R^{d}$ of dimension d, we consider a finite set of q discrete velocities $v_{j} \in V$ with components $v_{j}^{α}$ for $1 ⩽ α ⩽ 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}$ : $\begin{matrix} f_{j} = f_{j} (x, t), x \in R^{d}, t ⩾ 0, 0 ⩽ j < q . \end{matrix}$ The vector $f (x, t) \in R^{q}$ is constructed with the numbers $f_{j} (x, t)$ for $0 ⩽ j < q$ . We introduce a collision vector $R^{q} ∋ f ⟼ Q (f) \in R^{q}$ of components $Q_{j} (f)$ for $0 ⩽ j < q$ with given regular functions $Q_{j}$ . We introduce also a small parameter $ε > 0$ that can be interpreted as a Knudsen number in the context of gas dynamics (see e.g. [8]).

The Boltzmann model with discrete velocities takes the form $\begin{matrix} (1) & \frac{\partial f_{j}}{\partial t} + \sum_{1 ⩽ α ⩽ d} v_{j}^{α} \frac{\partial f_{j}}{\partial x_{α}} = \frac{1}{ε} Q_{j} (f), 0 ⩽ j < q . \end{matrix}$ It has been proposed by Carleman [7], Gross [32], Broadwell [6] and intensively developed by Gatignol [28] and her co-workers.

∙ Moments

Our framework concerns multi relaxation times: we introduce a constant invertible matrix M called “d’Humières matrix” [13] in this contribution. This matrix defines the vector of moments $m \in R^{q}$ by a simple product: $\begin{matrix} (2) & m_{k} \equiv \sum_{0 ⩽ j < q} M_{k j} f_{j} . \end{matrix}$ We introduce the number of conservation laws N ( $1 ⩽ N < q$ ) such that the N first moments of the collision kernel are equal to zero: $\begin{matrix} (3) & \sum_{0 ⩽ j < q} M_{k j} Q_{j} (f) = 0, \forall f \in R^{q}, 0 ⩽ k < N . \end{matrix}$ Then it is natural to divide the vector of moments into two families: $\begin{matrix} (4) & m \equiv (\begin{matrix} W \\ Y \end{matrix}) . \end{matrix}$ 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): $\begin{matrix} \frac{\partial W_{k}}{\partial t} + \sum_{1 ⩽ α ⩽ d} \sum_{0 ⩽ j < q} M_{k j} v_{j}^{α} \frac{\partial f_{j}}{\partial x_{α}} = 0, 0 ⩽ k < N . \end{matrix}$ Observe that the nonconserved moments or microscopic variables Y generate a linear space of dimension $q - N$ .

∙ Equilibrium states

We suppose that the equilibrim states $f^{eq}$ defined by the conditions $\begin{matrix} Q (f^{eq}) = 0 \end{matrix}$ are characterized with the help of a regular nonlinear vector field $Φ : R^{N} ⟶ R^{q - N}$ such that $\begin{matrix} (5) & f^{eq} = M^{- 1} (\begin{matrix} W \\ Φ (W) \end{matrix}) . \end{matrix}$ In other words, the vector field $W ⟼ Y^{eq} \equiv Φ (W)$ defines the set of equilibrium states.

We suppose moreover that the jacobian matrix $d Q (f^{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 $\begin{matrix} Z = diag (\frac{1}{τ_{1}}, \dots, \frac{1}{τ_{q - N}}) \end{matrix}$ of order $q - N$ with strictly positive coefficients $τ_{ℓ}$ such that $\begin{matrix} (6) & M d Q (f^{eq}) M^{- 1} = - (\begin{matrix} 0 & 0 \\ 0 & Z \end{matrix}) . \end{matrix}$ In other words, the jacobian operator $d Q (f^{eq})$ admits the matrix $M^{- 1}$ as a matrix of right eigenvectors and the associated eigenvalues are all non positive.

∙ Momentum-velocity operator

For a linear space of dimension d, we introduce the momentum-velocity operator matrix Λ defined by the relation $\begin{matrix} (7) & Λ = M diag (\sum_{1 ⩽ α ⩽ d} v^{α} \partial_{α}) M^{- 1} . \end{matrix}$ 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 Λ 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 $\begin{matrix} (8) & Λ \equiv (\begin{matrix} A & B \\ C & D \end{matrix}) . \end{matrix}$ Remember that in the following, the matrices A, B, C and D are matrices composed with first order space operators.

∙ Boltmann-BGK system with discrete velocities

We approximate now the previous discrete model with a BGK [4] 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: $\begin{matrix} Q (f) ≃ Q (f^{eq}) + d Q (f^{eq}) . (f - f^{eq}) . \end{matrix}$ Due to the condition $Q (f^{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 [4], a new system of partial differential equations: $\begin{matrix} (9) & \frac{\partial f_{j}}{\partial t} + \sum_{1 ⩽ α ⩽ d} v_{j}^{α} \frac{\partial f_{j}}{\partial x_{α}} = \frac{1}{ε} d Q_{j} (f^{eq}) . (f - f^{eq}), 0 ⩽ j < q . \end{matrix}$ Observe that $f^{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: $\begin{matrix} (10) & \frac{\partial W}{\partial t} + A W + B Y = 0, \frac{\partial Y}{\partial t} + C W + D Y = - \frac{1}{ε} Z (Y - Φ (W)) . \end{matrix}$

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

As ε 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: $\begin{matrix} (11) & \frac{\partial W}{\partial t} + Γ_{1} (W) + ε Γ_{2} (W) + ε^{2} Γ_{3} (W) + ε^{3} Γ_{4} (W) = O (ε^{4}) . \end{matrix}$ The nonlinear operators are related to the expansion of the microscopic moments: $\begin{matrix} (12) & Y = Φ (W) + ε Z^{- 1} Ψ_{1} (W) + ε^{2} Z^{- 1} Ψ_{2} (W) + ε^{3} Z^{- 1} Ψ_{2} (W) + O (ε^{4}) \end{matrix}$ and we have the interlaced relations $\begin{matrix} (13) & \{\begin{matrix} Γ_{1} (W) = A W + B Φ (W) \\ Ψ_{1} (W) = d Φ (W) . Γ_{1} (W) - (C W + D Φ (W)) \\ Γ_{2} (W) = B Z^{- 1} Ψ_{1} (W) \\ Ψ_{2} (W) = Z^{- 1} d Ψ_{1} (W) . Γ_{1} (W) + d Φ (W) . Γ_{2} (W) - D Z^{- 1} Ψ_{1} (W) \\ Γ_{3} (W) = B Z^{- 1} Ψ_{2} (W) \\ Ψ_{3} (W) = Z^{- 1} d Ψ_{1} (W) . Γ_{2} (W) + d Φ (W) . Γ_{3} (W) - D Z^{- 1} Ψ_{2} (W) \\ + Z^{- 1} d Ψ_{2} (W) . Γ_{1} (W) \\ Γ_{4} (W) = B Z^{- 1} Ψ_{3} (W) . \end{matrix} \end{matrix}$

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

Proof.
Proof of Proposition 1 First we inject the representation (12) inside the first equation of (10). It comes $\begin{matrix} \frac{\partial W}{\partial t} + A W + B (Φ (W) + ε Z^{- 1} Ψ_{1} (W) + ε^{2} Z^{- 1} Ψ_{2} (W) + ε^{3} Z^{- 1} Ψ_{2} (W)) = O (ε^{4}) \end{matrix}$ then after rearranging the terms, $\begin{matrix} \frac{\partial W}{\partial t} + A W + B Φ (W) + ε B Z^{- 1} Ψ_{1} (W) + ε^{2} B Z^{- 1} Ψ_{2} (W) + ε^{3} B Z^{- 1} Ψ_{3} (W) = O (ε^{4}) . \end{matrix}$ The confrontation of this partial differential equation with the Ansatz (11) shows that $\begin{matrix} \{\begin{matrix} Γ_{1} (W) = A W + B Φ (W) \\ Γ_{2} (W) = B Z^{- 1} Ψ_{1} (W) \\ Γ_{3} (W) = B Z^{- 1} Ψ_{2} (W) \\ Γ_{4} (W) = B Z^{- 1} Ψ_{3} (W) \end{matrix} \end{matrix}$ as proposed in (13).

Secondly, we write the second equation of (10) under the form $\begin{matrix} Y = Φ - ε Z^{- 1} (\partial_{t} Y + C W + D Y) . \end{matrix}$ We observe that $\begin{array}{l} - \partial_{t} Y = & d Y . (- \partial_{t} W) \\ = & d (Φ (W) + ε Z^{- 1} Ψ_{1} (W) + ε^{2} Z^{- 1} Ψ_{2} (W)) . (Γ_{1} + ε Γ_{2} (W) + ε^{2} Γ_{3} (W)) + O (ε^{3}) \\ = & d Φ (W) . Γ_{1} + ε (Z^{- 1} d Ψ_{1} (W) . Γ_{1} + d Φ (W) . Γ_{2}) \\ + ε^{2} (Z^{- 1} d Ψ_{2} (W) . Γ_{1} + Z^{- 1} d Ψ_{1} (W) . Γ_{2} + d Φ (W) . Γ_{3}) + O (ε^{3}) . \end{array}$ We report this expression inside the value of Y: $\begin{array}{l} Y = & Φ + ε Z^{- 1} (- \partial_{t} Y - C W - D Y) \\ = & Φ + ε Z^{- 1} [d Φ (W) . Γ_{1} + ε (Z^{- 1} d Ψ_{1} (W) . Γ_{1} + d Φ (W) . Γ_{2}) \\ + ε^{2} (Z^{- 1} d Ψ_{2} (W) . Γ_{1} + Z^{- 1} d Ψ_{1} (W) . Γ_{2} + d Φ (W) . Γ_{3}) - C W \\ - D (Φ (W) + ε Z^{- 1} Ψ_{1} (W) + ε^{2} Z^{- 1} Ψ_{2} (W))] + O (ε^{4}) \\ = & Φ + ε Z^{- 1} (d Φ (W) . Γ_{1} - C W - D Φ) + ε^{2} Z^{- 1} (Z^{- 1} d Ψ_{1} (W) . Γ_{1} + d Φ (W) . Γ_{2} - D Z^{- 1} Ψ_{1}) \\ + ε^{3} Z^{- 1} (Z^{- 1} d Ψ_{2} (W) . Γ_{1} + Z^{- 1} d Ψ_{1} (W) . Γ_{2} + d Φ (W) . Γ_{3} - D Z^{- 1} Ψ_{2}) + O (ε^{4}) . \end{array}$ Comparatively to the relation (12), we obtain $\begin{matrix} \{\begin{matrix} Ψ_{1} (W) = d Φ (W) . Γ_{1} - (C W + D Φ) \\ Ψ_{2} (W) = Z^{- 1} d Ψ_{1} (W) . Γ_{1} + d Φ (W) . Γ_{2} (W) - D Z^{- 1} Ψ_{1} (W) \\ Ψ_{3} (W) = Z^{- 1} d Ψ_{1} (W) . Γ_{2} + d Φ (W) . Γ_{3} - D Z^{- 1} Ψ_{2} + Z^{- 1} d Ψ_{2} (W) . Γ_{1} \end{matrix} \end{matrix}$ and the last relations of (13) are explicited. The proof is completed. □

3. Multiple relaxation times lattice Boltzmann schemes

We introduce now a regular cartesian lattice $L$ composed by vertices x separated by distances that are simple expressions of the space step $Δ x$ . A discrete time t is supposed to be an integer multiple of a time step $Δ 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 ⩽ j < q$ : $\begin{matrix} f_{j} = f_{j} (x, t), x \in L, t = n Δ t, n \in N, v_{j} \in V . \end{matrix}$

The discrete velocity set $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) ⟶ R (Δ t) f (t) \equiv f^{*}$ for the approximation of the system of ordinary differential equations $\begin{matrix} \frac{\partial f}{\partial t} = \frac{1}{ε} d Q (f^{eq}) . (f - f^{eq}) \end{matrix}$ and secondly a transport step $f^{*} ⟶ T (Δ t) f^{*} \equiv f (t + Δ t)$ for the free advection $\begin{matrix} \frac{\partial f_{j}}{\partial t} + \sum_{1 ⩽ α ⩽ d} v_{j}^{α} \frac{\partial f_{j}}{\partial x_{α}} = 0, 0 ⩽ j < q . \end{matrix}$

∙ Relaxation step

Due to the formulation (10), the relaxation scheme is simply described in the moment representation and we have $\begin{matrix} \frac{\partial W}{\partial t} = 0, \frac{\partial Y}{\partial t} = - \frac{1}{ε} Z (Y - Φ (W)) . \end{matrix}$ After relaxation, the vector of moments m is transformed into a new vector $m^{*}$ : $\begin{matrix} (14) & m^{*} \equiv (\begin{matrix} W^{*} \\ Y^{*} \end{matrix}) . \end{matrix}$ The conserved moments are not modified during this relaxation step: $\begin{matrix} (15) & W^{*} = W . \end{matrix}$ As remarked in [16], we use an explicit first order Euler scheme for the microscopic moments: $\begin{matrix} \frac{Y^{*} - Y}{Δ t} = - \frac{1}{ε} Z (Y - Φ (W)), \end{matrix}$ id est $\begin{matrix} (16) & Y^{*} = Y + S (Φ (W) - Y) \end{matrix}$ with $\begin{matrix} (17) & S = \frac{Δ t}{ε} Z . \end{matrix}$ Usually, as pointed by Lallemand and Luo [43], the relaxation matrix S is a diagonal matrix, in coherence with the hypothesis concerning the diagonalization of the operator $d Q (f^{eq})$ presented in the previous section: $\begin{matrix} S = diag (s_{1}, \dots, s_{q - N}) . \end{matrix}$ Moreover, the matrix S is a fixed invertible square matrix of dimension $(q - N) \times (q - N)$ . In other words, the time step $Δ t$ is if the order of the Knudsen number ε in our asymptotic study. Finally the particle distribution $f^{*}$ after relaxation is given according to $\begin{matrix} (18) & f^{*} = M^{- 1} m^{*} . \end{matrix}$

∙ Necessary conditions for stability

We deduce from the relation (16) $\begin{matrix} Y^{*} - Φ (W) = (I - S) (Y - Φ (W)) \end{matrix}$ and a natural stability condition takes the form $‖ I - S ‖ ⩽ 1$ , id est $\begin{matrix} 0 ⩽ s_{j} ⩽ 2, 0 ⩽ j < q - N . \end{matrix}$ 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.

∙ 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} Δ t$ of the lattice $L$ to the vertex $x \in L$ . In other words, the advection step occurs with a Courant–Friedrichs–Lewy number $c f l \equiv 1$ for all discrete velocities $v_{j}$ .

∙ Lattice Boltzmann scheme as a Lie–Trotter splitting scheme

The iteration of the lattice Boltzmann scheme is finally given by the relation $\begin{matrix} (19) & f (t + Δ t) = T (Δ t) R (Δ t) f (t) \end{matrix}$ and it is usefull to explicit it for each component: $\begin{matrix} (20) & f_{j} (x, t + Δ t) = f_{j}^{*} (x - v_{j} Δ t, t), v_{j} \in V, x \in L . \end{matrix}$ Observe that this framework is very general. The so-called “BGK version” of the lattice Boltzmann scheme [5,37,50] corresponds to the choice $M = I$ and $Z = \frac{1}{τ} I$ . Usual multi-relaxation lattice Boltzmann schemes [13,21,22,43] suppose in general some orthogonality property for the d’Humières matrix. This paradigm permits also the introduction of two particle distributions [2,21,34,40] or eventually more with the vectorial schemes developed by Graille [30] or a variant proposed in [18]. In those last cases, the numbering $j ⟼ 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 [54]: $\begin{matrix} f (t + Δ t) = R (\frac{Δ t}{2}) T (Δ t) R (\frac{Δ t}{2}) f (t) . \end{matrix}$ In this case, as put in evidence by Dellar [11], one can write $\tilde{f} (t) \equiv R {(\frac{Δ t}{2})}^{- 1} f (t)$ and the previous Strang scheme can be written as $\begin{matrix} \tilde{f} (t + Δ t) = T (Δ t) R (\frac{Δ t}{2}) R (\frac{Δ t}{2}) \tilde{f} (t) . \end{matrix}$ This scheme, up to the shift of one half relaxation step, is identical to the Lie–Trotter splitting (19), except that the product $R (\frac{Δ t}{2}) R (\frac{Δ t}{2})$ plays the role of $R (Δ t)$ . A first possibility is to divide the time step $Δ t$ by a factor 2 in the relations (16), (17): $\begin{matrix} \overline{Y} = Y + \frac{1}{2} S (Φ - Y), Y^{*} = \overline{Y} + \frac{1}{2} S (Φ - \overline{Y}) . \end{matrix}$ Then $\begin{matrix} Y^{*} = (1 - S + \frac{1}{4} S^{2}) Y + (S - \frac{1}{4} S^{2}) Φ (W) \end{matrix}$ and over-relaxation is no more possible because $s_{j} - \frac{1}{4} s_{j}^{2} ⩽ 1$ for all $s_{j} \in R$ . So such a numerical scheme is not used by the engineers at our knowledge.

A second possibility is to enforce the relation $R (\frac{Δ t}{2}) = \sqrt{R (Δ t)}$ , id est $\begin{matrix} \overline{Y} = \sqrt{I - S} Y + (I - \sqrt{I - S}) Φ (W) . \end{matrix}$ But this relation has a sense only if $s_{j} ⩽ 1$ and over-relaxation is again excluded.

An other idea is to exchange the roles of advection and relaxation in the Strang scheme: $\begin{matrix} f (t + Δ t) = T (\frac{Δ t}{2}) R (Δ t) T (\frac{Δ t}{2}) f (t) \end{matrix}$ 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 ): $\begin{matrix} (21) & m (x, t + Δ t) = exp (- Δ t Λ) m^{*} (x, t) . \end{matrix}$

Proof of Proposition 2.
We have simply a long sequence of identities: $\begin{array}{l} m_{k} (x, t + Δ t) \\ = \sum_{j} M_{k j} f_{j} (x, t + Δ t) due to (2) \\ = \sum_{j} M_{k j} f_{j}^{} (x - v_{j} Δ t, t) due to (20) \\ = \sum_{j ℓ} M_{k j} {(M^{- 1})}_{j ℓ} m_{ℓ}^{} (x - v_{j} Δ t, t) due to (18) \\ = \sum_{j ℓ} M_{k j} {(M^{- 1})}_{j ℓ} \sum_{n = 0}^{\infty} \frac{1}{n!} {(- Δ t \sum_{α} v_{j}^{α} \partial_{α})}^{n} m_{ℓ}^{} (x, t) Taylor expansion \\ = \sum_{ℓ} \sum_{n = 0}^{\infty} \frac{1}{n!} \sum_{j} M_{k j} {(- Δ t \sum_{α} v_{j}^{α} \partial_{α})}^{n} {(M^{- 1})}_{j ℓ} m_{ℓ}^{} (x, t) elementary algebra \\ = \sum_{ℓ} \sum_{n = 0}^{\infty} \frac{1}{n!} {(- Δ t)}^{n} {({(M (\sum_{α} v^{α} \partial_{α}) M^{- 1})}^{n})}_{k ℓ} m_{ℓ}^{} (x, t) because M is a constant matrix \\ = \sum_{ℓ} \sum_{n = 0}^{\infty} \frac{1}{n!} {(- Δ t)}^{n} {(Λ^{n})}_{k ℓ} m_{ℓ}^{} (x, t) due to (7) \\ = \sum_{ℓ} [\sum_{n = 0}^{\infty} \frac{1}{n!} {(- Δ t Λ)}_{k ℓ}^{n}] m_{ℓ}^{} (x, t) elementary algebra \\ = \sum_{ℓ} exp {(- Δ t Λ)}_{k ℓ} m_{ℓ}^{} (x, t) elementary algebra \\ = {(exp (- Δ t Λ) m^{*} (x, t))}_{k} elementary algebra \end{array}$ and the proof is completed. □

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 [5].

∙ Momentum-velocity operator for a D2Q9 lattice Boltzmann scheme

Fig. 1.
D2Q9 lattice Boltzmann scheme.

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

To fix the ideas, we use the moment matrix presented in Lallemand and Luo [43]: $\begin{matrix} (22) & M_{D 2 Q 9} = [\begin{matrix} 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\ 0 & λ & 0 & - λ & 0 & λ & - λ & - λ & λ \\ 0 & 0 & λ & 0 & - λ & λ & λ & - λ & - λ \\ - 4 λ^{2} & - λ^{2} & - λ^{2} & - λ^{2} & - λ^{2} & 2 λ^{2} & 2 λ^{2} & 2 λ^{2} & 2 λ^{2} \\ 0 & λ^{2} & - λ^{2} & λ^{2} & - λ^{2} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & λ^{2} & - λ^{2} & λ^{2} & - λ^{2} \\ 0 & - 2 λ^{3} & 0 & 2 λ^{3} & 0 & λ^{3} & - λ^{3} & - λ^{3} & λ^{3} \\ 0 & 0 & - 2 λ^{3} & 0 & 2 λ^{3} & λ^{3} & λ^{3} & - λ^{3} & - λ^{3} \\ 4 λ^{4} & - 2 λ^{4} & - 2 λ^{4} & - 2 λ^{4} & - 2 λ^{4} & λ^{4} & λ^{4} & λ^{4} & λ^{4} \end{matrix}] . \end{matrix}$ Observe that other choices are possible with the Hermite moments developed in [48,53]. The moments are associated with the lines of the matrix (22). They are denominated (in this order) by $ρ, J_{x}, J_{y}, ε, x_{x}, x_{y}, q_{x}, q_{y}, h$ . The lines of this invertible matrix are chosen orthogonal. For isothermal flows, the conserved moments $\begin{matrix} (23) & W = {(ρ, J_{x}, J_{y})}^{t} \end{matrix}$ correspond to the three first lines of the d’Humières matrix (22). The non-conserved moments complete the family: $\begin{matrix} (24) & Y = {(ε, x_{x}, x_{y}, q_{x}, q_{y}, h)}^{t} . \end{matrix}$ They are linked to the six last lines of the matrix defined in (22). With the velocities defined in Fig. 1 and the moment matrix (22), we construct without difficulty the operator matrix Λ for the thermal D2Q9 scheme:

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 $\begin{matrix} Λ^{2} \equiv (\begin{matrix} A_{2} & B_{2} \\ C_{2} & D_{2} \end{matrix}) = (\begin{matrix} A & B \\ C & D \end{matrix}) (\begin{matrix} A & B \\ C & D \end{matrix}) = (\begin{matrix} A^{2} + B C & A B + B D \\ C A + D C & C B + D^{2} \end{matrix}) \end{matrix}$ and to fix the ideas $\begin{matrix} Λ^{3} \equiv (\begin{matrix} A_{3} & B_{3} \\ C_{3} & D_{3} \end{matrix}) = (\begin{matrix} A_{2} & B_{2} \\ C_{2} & D_{2} \end{matrix}) (\begin{matrix} A & B \\ C & D \end{matrix}) = (\begin{matrix} A_{2} A + B_{2} C & A_{2} B + B_{2} D \\ C_{2} A + D_{2} C & C_{2} B + D_{2} D \end{matrix}) . \end{matrix}$
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.

∙ Hypotheses for a formal expansion

The formal Taylor expansion supposes some precise hypotheses. First, we adopt the acoustic scaling: the ratio $λ \equiv \frac{Δ x}{Δ t}$ is supposed to be constant in all this work. Moreover, the ratio $\frac{Δ t}{ε}$ 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 [36]. 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 Σ in this contribution. Roughtly speaking, it plays a role analogous to the matrix $Z^{- 1}$ in the relations (13). But it is defined by $\begin{matrix} (26) & Σ \equiv S^{- 1} - \frac{1}{2} I . \end{matrix}$ 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 Σ 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 $d ζ (W) . ξ$ for the action of the differential of some regular function $W ⟼ ζ (W)$ againts a test vector ξ. In particular, $\begin{matrix} (27) & γ_{j} (W) \equiv d Φ (W) . Γ_{j} (W), 1 ⩽ j ⩽ 3 . \end{matrix}$ For second order derivatives, we introduce $\begin{matrix} (28) & \partial^{2} Ψ . Γ_{1} \equiv d (d Ψ (W) . Γ_{1}) . Γ_{1} = d^{2} Ψ (W) . (Γ_{1}, Γ_{1}) + d Ψ (W) . d Γ_{1} (W) . Γ_{1} . \end{matrix}$ The asymptotic analysis occurs for a time step $Δ t$ (or a space step $Δ x$ ) tending to zero. Emerging partial differential equations for the conserved variables are denoted as $\begin{matrix} (29) & \partial_{t} W + Γ_{1} (W) + Δ t Γ_{2} (W) + Δ t^{2} Γ_{3} (W) + Δ t^{3} Γ_{4} (W) = O (Δ t^{4}) . \end{matrix}$ The vector $Γ_{j} (W)$ belongs in $R^{N}$ . It is obtained after j space derivations of the conserved moments W and of the equilibrium vector $Φ (W)$ . For this reason, we will speak in the following of e.g. $Γ_{1} (W)$ as “first order term”, even if it is a term of order zero relative to the infinitesimal $Δ t$ . For non-conserved moments (or microscopic variables), we suppose $\begin{matrix} (30) & \{\begin{matrix} Y = Φ (W) + S^{- 1} (Δ t Ψ_{1} (W) + Δ t^{2} Ψ_{2} (W) + Δ t^{3} Ψ_{3} (W)) + O (Δ t^{4}) \\ Y^{*} = Φ (W) + (Σ - \frac{1}{2} I) (Δ t Ψ_{1} (W) + Δ t^{2} Ψ_{2} (W) + Δ t^{3} Ψ_{3} (W)) + O (Δ t^{4}) \end{matrix} \end{matrix}$ where 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 $Φ_{j} (W)$ are analogous to $Γ_{j} (W)$ but not with the same dimension: $Γ_{j} (W) \in R^{N}$ and $Φ_{j} (W) \in R^{q - N}$ .

∙ 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 $Γ_{j} (W)$ and $Φ_{j} (W)$ up to fourth order accuracy: $\begin{matrix} (31) & \{\begin{matrix} Γ_{1} (W) = A W + B Φ (W) \\ Ψ_{1} (W) = d Φ (W) . Γ_{1} (W) - (C W + D Φ (W)) \\ Γ_{2} (W) = B Σ Ψ_{1} (W) \\ Ψ_{2} (W) = Σ d Ψ_{1} (W) . Γ_{1} (W) + d Φ (W) . Γ_{2} (W) - D Σ Ψ_{1} (W) \\ Γ_{3} (W) = B Σ Ψ_{2} (W) + \frac{1}{12} B_{2} Ψ_{1} (W) - \frac{1}{6} B d Ψ_{1} (W) . Γ_{1} (W) \\ Ψ_{3} (W) = Σ d Ψ_{1} (W) . Γ_{2} (W) + d Φ (W) . Γ_{3} (W) - D Σ Ψ_{2} (W) + Σ d Ψ_{2} (W) . Γ_{1} (W) \\ + \frac{1}{6} D d Ψ_{1} (W) . Γ_{1} (W) - \frac{1}{12} D_{2} Ψ_{1} (W) - \frac{1}{12} \partial^{2} Ψ_{1} (W) . Γ_{1} (W) \\ Γ_{4} (W) = B Σ Ψ_{3} (W) + \frac{1}{4} B_{2} Ψ_{2} (W) + \frac{1}{6} B D_{2} Σ Ψ_{1} (W) - \frac{1}{6} A B Ψ_{2} (W) \\ - \frac{1}{6} B (d γ_{1} (W) . Γ_{2} (W) + d γ_{2} (W) . Γ_{1} (W)) - \frac{1}{6} B Σ \partial^{2} Ψ_{1} (W) . Γ_{1} (W) . \end{matrix} \end{matrix}$ 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 $Γ_{2} = B Σ Ψ_{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 Σ 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 $Δ t$ tends to zero, the microscopic moments are close to their equilibrium value: $\begin{matrix} (32) & Y = Φ (W) + O (Δ t), Y^{*} = Φ (W) + O (Δ t) . \end{matrix}$

Proof of Proposition 3.
We expand one iteration of the scheme (21) at order zero: $\begin{matrix} (33) & m + O (Δ t) = m^{} + O (Δ t^{3}) . \end{matrix}$ Then we deduce from (33): $\begin{matrix} Y - Y^{} = O (Δ t) . \end{matrix}$ Then due to the iteration (16) and the fact that the matrix S is supposed fixed, we have $\begin{matrix} (34) & Y = Φ (W) + O (Δ t) . \end{matrix}$ 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. □

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 $Δ t$ tends to zero, the macroscopic moments W satisfy asymptotically the following system of N first order equations: $\begin{matrix} (35) & \partial_{t} W + Γ_{1} (W) = O (Δ t) . \end{matrix}$ Moreover, the vector $Γ_{1} (W)$ for the first order dynamics satisfy $\begin{matrix} (36) & Γ_{1} (W) = A W + B Φ (W) . \end{matrix}$

Proof of Proposition 4.
We expand one iteration of the scheme (21) at first order: $\begin{matrix} (37) & m + Δ t \partial_{t} m + O (Δ t^{2}) = m^{} - Δ t Λ m^{} + O (Δ t^{2}) . \end{matrix}$ We can replace the vector m by its two components W and Y. We have for the first line $\begin{matrix} W + Δ t \partial_{t} W + O (Δ t^{2}) = W - Δ t (A W + B Y^{*}) + O (Δ t^{2}) . \end{matrix}$ We simplify by the constant W term, divide by $Δ t$ and take into account the expansion (32): $\begin{matrix} \partial_{t} W + O (Δ t) = - (A W + B Φ (W)) + O (Δ t) . \end{matrix}$ This relation is exactly the expansion (35) with the first order dynamics $Γ_{1} (W)$ evaluated according to the relation (36). □
∙ Example of the D2Q9 fluid scheme

The results presented here are essentially a reformulation of the classic work of Lallemand and Luo [43]. With the moments introduced in (22), (23) and (24), we denote by $Φ_{ε}$ , $Φ_{x x}$ and $Φ_{x y}$ 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, $\begin{matrix} Γ_{1} = (\begin{matrix} \partial_{x} J_{x} + \partial_{y} J_{y} \\ \frac{2}{3} λ^{2} \partial_{x} ρ + \frac{1}{6} \partial_{x} Φ_{ε} + \frac{1}{2} \partial_{x} Φ_{x x} + \partial_{y} Φ_{x y} \\ \frac{2}{3} λ^{2} \partial_{y} ρ + \frac{1}{6} \partial_{y} Φ_{ε} - \frac{1}{2} \partial_{y} Φ_{x x} + \partial_{x} Φ_{x y} \end{matrix}) . \end{matrix}$ These relations can be easily fitted with the Euler equations of gas dynamics in two space dimensions $\begin{matrix} \{\begin{matrix} \partial_{t} ρ + \partial_{x} J_{x} + \partial_{y} J_{y} = 0 \\ \partial_{t} J_{x} + \partial_{x} (\frac{1}{ρ} J_{x}^{2} + p) + \partial_{y} (\frac{1}{ρ} J_{x} J_{y}) = 0 \\ \partial_{t} J_{y} + \partial_{x} (\frac{1}{ρ} J_{x} J_{y}) + \partial_{y} (\frac{1}{ρ} J_{y}^{2} + p) = 0 . \end{matrix} \end{matrix}$

We identify the two expressions with $J_{x} \equiv ρ u$ and $J_{y} \equiv ρ v$ : $\begin{matrix} \frac{2}{3} λ^{2} ρ + \frac{1}{6} Φ_{ε} + \frac{1}{2} Φ_{x x} = ρ u^{2} + p, Φ_{x y} = ρ u v \frac{2}{3} λ^{2} ρ + \frac{1}{6} Φ_{ε} - \frac{1}{2} Φ_{x x} = ρ v^{2} + p, \end{matrix}$ and we deduce the expression for the three first noncenserved moments: $\begin{matrix} Φ_{ε} = 6 p - 4 λ^{2} ρ + 3 ρ (u^{2} + v^{2}), Φ_{x x} = ρ (u^{2} - v^{2}), Φ_{x y} = ρ u v . \end{matrix}$ 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 $Δ t$ (or a space step $Δ x$ ) tending to zero. The expansion of the microscopic variables takes the form $\begin{matrix} (38) & Y = Φ (W) + Δ t S^{- 1} Ψ_{1} (W) + O (Δ t^{2}) \end{matrix}$ and the vector $Ψ_{1} (W)$ satisfies the relation $\begin{matrix} (39) & Ψ_{1} (W) = d Φ (W) . Γ_{1} (W) - (C W + D Φ (W)) . \end{matrix}$ A set of second order partial differential equations emerges $\begin{matrix} (40) & \partial_{t} W + Γ_{1} (W) + Δ t Γ_{2} (W) = O (Δ t^{2}) . \end{matrix}$ The vector $Γ_{1} (W)$ has been precised in ( 36 ) the vector $Γ_{2} (W)$ is obtained after two derivations of the conserved moments W and the equilibrium vector $Φ (W)$ : $\begin{matrix} (41) & Γ_{2} (W) = B Σ Ψ_{1} (W) . \end{matrix}$ The expression of the second order term puts in evidence the importance of taking over-relaxation in the applications where Σ should also be as small as possible.

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 $\begin{matrix} Y + Δ t \partial_{t} Y + O (Δ t^{2}) = Y^{} - Δ t (C W + D Y^{}) + O (Δ t^{2}) . \end{matrix}$ Then $\begin{matrix} Y - Y^{} = - Δ t (\partial_{t} Y + (C W + D Φ (W))) + O (Δ t^{2}) . \end{matrix}$ The explicitation of $\partial_{t} Y$ is given by the chain rule from the expansion (34): $\begin{array}{l} \partial_{t} Y & = \partial_{t} (Φ (W) + O (Δ t)) \\ = d Φ (W) . \partial_{t} W + O (Δ t) \\ = d Φ (W) . (- Γ_{1}) + O (Δ t) . \end{array}$ Then $\begin{matrix} S (Y - Φ (W)) = Δ t (d Φ (W) . Γ_{1} - (C W + D Φ)) + O (Δ t^{2}) . \end{matrix}$ The relation (38) is established, and the expression of $Ψ_{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 $\begin{array}{l} Y^{} & = Y + S (Φ (W) - Y) \\ = Y - (Δ t Ψ_{1} + O (Δ t^{2})) due to (38) \\ = Φ (W) + (Σ + \frac{1}{2} I) (Δ t Ψ_{1} (W) + O (Δ t^{2})) - (Δ t Ψ_{1} + O (Δ t^{2})) due to (26) and (38) \\ = Φ (W) + (Σ - \frac{1}{2} I) Δ t Ψ_{1} (W) + O (Δ t^{2}) \end{array}$ and we have the relation $\begin{matrix} (42) & Y^{} = Φ (W) + (Σ - \frac{1}{2} I) Δ t Ψ_{1} (W) + O (Δ t^{2}) . \end{matrix}$ The two relations (38) and (42) are nothing else that the first order terms of the general expansion (30).

∙ Second order partial differential equations

We expand now one iteration of the scheme (21) at second order accuracy: $\begin{matrix} (43) & m + Δ t \partial_{t} m + \frac{1}{2} Δ t^{2} \partial_{t}^{2} m + O (Δ t^{3}) = m^{} - Δ t Λ m^{} + \frac{1}{2} Δ t^{2} Λ^{2} m^{} + O (Δ t^{3}) . \end{matrix}$ We replace the vector m by its two components W and Y. We use the decomposition (8) of the momentum-velocity operator matrix Λ. We have for the first component $\begin{matrix} (44) & \{\begin{matrix} W + Δ t \partial_{t} W + \frac{1}{2} Δ t^{2} \partial_{t}^{2} W + O (Δ t^{3}) \\ = W - Δ t (A W + B Y^{}) + \frac{1}{2} Δ t^{2} (A_{2} W + B_{2} Y^{}) + O (Δ t^{3}) . \end{matrix} \end{matrix}$ We explicit some terms present in the relation (44). We have $\begin{matrix} \partial_{t} W = - Γ_{1} (W) + O (Δ t) = - (A W + B Φ (W)) + O (Δ t) \end{matrix}$ Then, using the formal rule stated at the beginning of Section 4, we keep the order of accuracy after derivation: $\begin{array}{l} \partial_{t}^{2} W & = - \partial_{t} (Γ_{1} + O (Δ t)) \\ = - d Γ_{1} . \partial_{t} W + O (Δ t) \\ = d Γ_{1} . Γ_{1} + O (Δ t) \\ = A Γ_{1} + B d Φ . Γ_{1} + O (Δ t) . \end{array}$

We deduce from (44): $\begin{matrix} \partial_{t} W = - \frac{1}{2} Δ t \partial_{t}^{2} W - (A W + B Y^{}) + \frac{1}{2} Δ t (A_{2} W + B_{2} Y^{*}) + O (Δ t^{2}) . \end{matrix}$ We order the various terms by powers of $Δ t$ : $\begin{array}{l} \partial_{t} W = & - A W - B (Φ + (Σ - \frac{1}{2} I) Δ t Ψ_{1}) - \frac{1}{2} Δ t (A Γ_{1} + B d Φ . Γ_{1}) \\ + \frac{1}{2} Δ t ((A^{2} + B C) W + (A B + B D) Φ) + O (Δ t^{2}) \\ = & - A W - B Φ + Δ t [- B Σ Ψ_{1} + \frac{1}{2} B (d Φ . Γ_{1} - C W - D Φ) - \frac{1}{2} A (A W + B Φ) \\ - \frac{1}{2} B d Φ . Γ_{1} + \frac{1}{2} (A^{2} + B C) W + \frac{1}{2} (A B + B D) Φ] + O (Δ t^{2}) \\ = & - Γ_{1} (W) - Δ t B Σ Ψ_{1} (W) + O (Δ t^{2}) \end{array}$ and the relation (40) is proven with $Γ_{2} (W) = B Σ Ψ_{1}$ as proposed in (41). □
∙ Link with the original Taylor expansion method

In [15 ,16], the moments at equilibrium are denoted $m_{ℓ}^{eq}$ : $\begin{matrix} m_{ℓ}^{eq} (W) = \{\begin{matrix} W_{ℓ} & if 0 ⩽ ℓ < N \\ {(Φ (W))}_{ℓ - N} & if ℓ ⩾ N, \end{matrix} \end{matrix}$ the coefficients ${\tilde{Λ}}_{i ℓ}^{β}$ are given according to $\begin{matrix} (45) & {\tilde{Λ}}_{i ℓ}^{β} = \sum_{j} M_{i j} v_{j}^{β} {(M^{- 1})}_{j ℓ}, 1 ⩽ β ⩽ d, 0 ⩽ i, ℓ < q, \end{matrix}$ are closely related to the operator matrix Λ introduced in (7). The coefficients of the Hénon’s matrix (26) have a shifted numbering $\begin{matrix} (46) & {\tilde{σ}}_{ℓ} = σ_{ℓ - N} \equiv \frac{1}{s_{ℓ - N}} - \frac{1}{2}, ℓ ⩾ N \end{matrix}$ and the defect of conservation ${\tilde{θ}}_{k}$ is defined according to $\begin{matrix} (47) & {\tilde{θ}}_{k} (W) \equiv \partial_{t} m_{k}^{eq} + \sum_{0 ⩽ ℓ < q} \sum_{1 ⩽ β ⩽ d} {\tilde{Λ}}_{k ℓ}^{β} \partial_{β} m_{ℓ}^{eq} (W), k ⩾ N . \end{matrix}$ 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 ]: $\begin{matrix} (48) & \partial_{t} W_{i} + {\tilde{Λ}}_{i ℓ}^{β} \partial_{β} m_{ℓ}^{eq} (W) = Δ t \sum_{k ⩾ N} {\tilde{Λ}}_{i k}^{β} {\tilde{σ}}_{k} \partial_{β} {\tilde{θ}}_{k} (W) + O (Δ t^{2}) . \end{matrix}$

Proof of Proposition 6.
Remark first that $\begin{matrix} Λ_{i k} = \sum_{j = 0}^{q - 1} M_{i j} (\sum_{1 ⩽ β ⩽ d} v_{j}^{β} \partial_{β}) {(M^{- 1})}_{j k} = \sum_{1 ⩽ β ⩽ d} (\sum_{j = 0}^{q - 1} M_{i j} v_{j}^{β} {(M^{- 1})}_{j k}) \partial_{β} = \sum_{1 ⩽ β ⩽ d} {\tilde{Λ}}_{i k}^{β} \partial_{β} . \end{matrix}$ Then on one hand, the component number i of the first order term $(A W + B Φ (W))$ can be expanded in the following way: $\begin{array}{l} {(A W + B Φ (W))}_{i} & = \sum_{k < N} {\tilde{Λ}}_{i k}^{β} \partial_{β} W_{k} + \sum_{k ⩾ N} {\tilde{Λ}}_{i k}^{β} \partial_{β} Φ_{k} \\ = \sum_{k < N} {\tilde{Λ}}_{i k}^{β} \partial_{β} m_{k}^{eq} + \sum_{k ⩾ N} {\tilde{Λ}}_{i k}^{β} \partial_{β} m_{k}^{eq} \\ = \sum_{0 ⩽ k < q} {\tilde{Λ}}_{i k}^{β} \partial_{β} m_{k}^{eq} \end{array}$ and the first order terms of (40) and (48) are identical. Secondly, due to (47), we have for $k ⩾ N$ $\begin{array}{l} {\tilde{θ}}_{k} & = \partial_{t} m_{k}^{eq} + \sum_{ℓ β} {\tilde{Λ}}_{k ℓ}^{β} \partial_{β} m_{ℓ}^{eq} \\ = \partial_{t} Φ_{k - N} + {(Λ . m^{eq})}_{k} \\ = {(d Φ . \partial_{t} W)}_{k - N} + {(Λ . m^{eq})}_{k} \\ = {(d Φ . (Γ_{1} + O (Δ t)))}_{k - N} - {(C W + D Φ)}_{k - N} \end{array}$ and due to (39), $\begin{matrix} {\tilde{θ}}_{k} = - {(Ψ_{1})}_{k - N} + O (Δ t), k ⩾ N . \end{matrix}$ We deduce that for second order terms with $0 ⩽ i < N$ : $\begin{array}{l} {(B \tilde{σ} Ψ_{1})}_{i} & = \sum_{0 ⩽ ℓ < q - N} B_{i ℓ} {\tilde{σ}}_{ℓ + N} {(Ψ_{1} (W))}_{ℓ} \\ = \sum_{k ⩾ N} \sum_{0 ⩽ β ⩽ d} {\tilde{Λ}}_{i k}^{β} \partial_{β} (σ_{k} (- {\tilde{θ}}_{k} + O (Δ t))) with k = ℓ + N \\ = - \sum_{k ⩾ N} \sum_{0 ⩽ β ⩽ d} {\tilde{Λ}}_{i k}^{β} {\tilde{σ}}_{k} (\partial_{β} {\tilde{θ}}_{k}) + O (Δ t) \end{array}$ because the coefficients ${\tilde{σ}}_{k}$ are supposed constant. Then the relation (48) is established. □
∙ Fluid D2Q9 diffusive tensor

If we compute the second order term $Γ_{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 μ and a given bulk viscosity ζ, we write the viscous terms of the lattice Boltzmann expansion as $\begin{matrix} (49) & - Δ t Γ_{2} = (\begin{matrix} 0 \\ \partial_{j} τ_{x j} \equiv \partial_{x} (2 μ \partial_{x} u + (ζ - μ) (\partial_{x} u + \partial_{y} v)) + \partial_{y} (μ (\partial_{x} v + \partial_{y} u)) \\ \partial_{j} τ_{y j} \equiv \partial_{x} (μ (\partial_{x} v + \partial_{y} u)) + \partial_{y} ((ζ - μ) (\partial_{x} u + \partial_{y} v) + 2 μ \partial_{y} v) \end{matrix}) + (\begin{matrix} 0 \\ Ψ_{x} \\ Ψ_{y} \end{matrix}), \end{matrix}$ with the notation $Ψ_{x}$ , $Ψ_{y}$ for the error for the dissipation of momentum.

The problem to enforce $Ψ_{x} = Ψ_{y} = 0$ is ill-posed, as recognized by Dellar [10,12], Házi and Kávrán [35] and Prasianakis at al. [49] 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 $\begin{matrix} p (ρ) = \frac{λ^{2}}{3} ρ, \end{matrix}$ id est a sound velocity equal to $\frac{1}{\sqrt{3}}$ in the velocity unit of the lattice Boltzmann scheme. Nevertheless, with the notation $Φ_{q x}$ and $Φ_{q y}$ for the equilibrium fonctions of the moments $q_{x}$ and $q_{y}$ and the coefficients of the Hénon matrix Σ fixed as $\begin{matrix} (50) & Σ = diag (σ_{e}, σ_{x}, σ_{x}, σ_{q}, σ_{q}, σ_{h}), \end{matrix}$ the choice $\begin{matrix} \{\begin{matrix} Φ_{q x} = - ρ λ^{2} u + 3 ρ (u^{2} + v^{2}) u \\ Φ_{q y} = - ρ λ^{2} v + 3 ρ (u^{2} + v^{2}) v, \end{matrix} \end{matrix}$ conducts to the relation (49) with associated viscosities fixed in a very classical way: $\begin{matrix} (51) & μ = \frac{λ}{3} ρ σ_{x} Δ x, ζ = \frac{λ}{3} ρ σ_{e} Δ x . \end{matrix}$ The equilibrium value $Φ_{h}$ of the last fourth order moment h (see the relation (24)) has no incidence on the value of the diffusive term $Γ_{2}$ . The error is then reduced to third order terms relative to the velocity and we have precisely $\begin{matrix} \{\begin{matrix} Ψ_{x} = σ_{x} Δ t [\partial_{x} (u^{3} \partial_{x} ρ - v^{3} \partial_{y} ρ + 3 ρ (u^{2} \partial_{x} u - v^{2} \partial_{y} v)) \\ + \partial_{y} (- v^{3} \partial_{x} ρ - u^{3} \partial_{y} ρ - 3 ρ (u^{2} \partial_{y} u + v^{2} \partial_{x} v))] \\ Ψ_{y} = σ_{x} Δ t [\partial_{x} (- v^{3} \partial_{x} ρ - u^{3} \partial_{y} ρ - 3 ρ (u^{2} \partial_{y} u + v^{2} \partial_{x} v) \\ + \partial_{y} (- u^{3} \partial_{x} ρ + v^{3} \partial_{y} ρ + 3 ρ (- u^{2} \partial_{x} u + v^{2} \partial_{y} v))] . \end{matrix} \end{matrix}$
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 $Δ t$ has an infinitesimal value, the expansion of the microscopic variables takes the form $\begin{matrix} (52) & Y = Φ (W) + Δ t S^{- 1} Ψ_{1} (W) + Δ t^{2} S^{- 1} Ψ_{2} (W) + O (Δ t^{3}) . \end{matrix}$ The vector $Ψ_{1} (W)$ is evaluated at the relation ( 39 ) and $Ψ_{2} (W)$ satisfies the relation $\begin{matrix} (53) & Ψ_{2} (W) = Σ d Ψ_{1} (W) . Γ_{1} (W) + d Φ (W) . Γ_{2} (W) - D Σ Ψ_{1} (W) . \end{matrix}$ It is possible to precise a system of third order partial differential equations $\begin{matrix} (54) & \partial_{t} W + Γ_{1} (W) + Δ t Γ_{2} (W) + Δ t^{2} Γ_{3} (W) = O (Δ t^{3}) . \end{matrix}$ The vectors $Γ_{1} (W)$ and $Γ_{2} (W)$ have been precised at the relations ( 36 ) and ( 41 ) respectively. We have moreover $\begin{matrix} (55) & Γ_{3} (W) = B Σ Ψ_{2} (W) + \frac{1}{12} B_{2} Ψ_{1} (W) - \frac{1}{6} B d Ψ_{1} (W) . Γ_{1} (W) . \end{matrix}$

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 $\begin{array}{l} Y + Δ t \partial_{t} Y + \frac{1}{2} Δ t^{2} \partial_{t}^{2} Y + O (Δ t^{3}) = Y^{} - Δ t (C W + D Y^{}) + \frac{1}{2} Δ t^{2} (C_{2} W + D_{2} Y^{}) + O (Δ t^{3}) \end{array}$ and, due to the relaxation process $Y - Y^{} \equiv S (Y - Φ (W))$ , we have $\begin{matrix} (56) & S (Y - Φ (W)) = - Δ t (C W + D Y^{} + \partial_{t} Y) + \frac{1}{2} Δ t^{2} (C_{2} W + D_{2} Y^{} - \partial_{t}^{2} Y) + O (Δ t^{3}) . \end{matrix}$ We can precise the two first time derivatives of the non-conserved moments. The relation (38) can also be written as $\begin{matrix} Y = Φ (W) + (Σ + \frac{1}{2} I) Δ t Ψ_{1} (W) + O (Δ t^{2}) . \end{matrix}$ We differentiate this previous relation relative to time: $\begin{array}{l} \partial_{t} Y & = d Φ . \partial_{t} W + (Σ + \frac{1}{2} I) Δ t d Ψ_{1} . \partial_{t} W + O (Δ t^{2}) \\ = - d Φ . (Γ_{1} (W) + Δ t Γ_{2} (W) + O (Δ t^{2})) + Δ t (Σ + \frac{1}{2} I) d Ψ_{1} . (- Γ_{1} + O (Δ t)) + O (Δ t^{2}) \end{array}$ and $\begin{matrix} (57) & \partial_{t} Y = - d Φ (W) . Γ_{1} (W) - Δ t [d Φ (W) . Γ_{2} (W) + (Σ + \frac{1}{2} I) d Ψ_{1} (W) . Γ_{1} (W)] + O (Δ t^{2}) . \end{matrix}$ 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 $\begin{matrix} (58) & \partial_{t}^{2} Y = d γ_{1} . Γ_{1} + O (Δ t) . \end{matrix}$ We precise now the expansions of $C W + D Y^{}$ and $C_{2} W + D_{2} Y^{}$ . We have $\begin{array}{l} C W + D Y^{} = C W + D [Φ (W) + (Σ - \frac{1}{2} I) Δ t Ψ_{1}] + O (Δ t^{2}) \end{array}$ and $\begin{matrix} (59) & C W + D Y^{} = γ_{1} (W) - Ψ_{1} (W) + Δ t D (Σ - \frac{1}{2} I) Ψ_{1} (W) + O (Δ t^{2}) . \end{matrix}$ With a new set of operator matrices, $\begin{array}{l} C_{2} W + D_{2} Y^{} & = C_{2} W + D_{2} (Φ (W) + O (Δ t)) \\ = (C A + D C) W + (C B + D^{2}) Φ (W) + O (Δ t) \\ = C Γ_{1} + D (d Φ . Γ_{1} - Ψ_{1}) + O (Δ t) \end{array}$ and because $d Φ . Γ_{1} = d γ_{1} . Γ_{1} - (C Γ_{1} + D d Φ . Γ_{1})$ , $\begin{matrix} (60) & C_{2} W + D_{2} Y^{} = d (γ_{1} - Ψ_{1}) . Γ_{1} (W) - D Ψ_{1} (W) + O (Δ t) . \end{matrix}$ Then due to (56), (57), (58), (59) and (60), we have $\begin{array}{l} S (Y - Φ (W)) \\ = - Δ t (C W + D Y^{} + \partial_{t} Y) + \frac{1}{2} Δ t^{2} (C_{2} W + D_{2} Y^{} - \partial_{t}^{2} Y) + O (Δ t^{3}) \\ = - Δ t (γ_{1} - Ψ_{1} + Δ t D (Σ - \frac{1}{2} I) Ψ_{1}) - Δ t (- d Φ . Γ_{1} - Δ t [d Φ . Γ_{2} + (Σ + \frac{1}{2} I) d Ψ_{1} . Γ_{1}]) \\ + \frac{1}{2} Δ t^{2} (d (γ_{1} - Ψ_{1}) . Γ_{1} - D Ψ_{1}) - \frac{1}{2} Δ t^{2} d γ_{1} . Γ_{1} + O (Δ t^{3}) \\ = Δ t Ψ_{1} + Δ t^{2} (- D Σ Ψ_{1} + \frac{1}{2} D Ψ_{1} + d Φ . Γ_{2} + Σ d Ψ_{1} . Γ_{1} - \frac{1}{2} D Ψ_{1}) + O (Δ t^{3}) \\ = Δ t Ψ_{1} + Δ t^{2} (- D Σ Ψ_{1} + d Φ . Γ_{2} + Σ d Ψ_{1} . Γ_{1}) + O (Δ t^{3}) \end{array}$ and the relations (52) and (53) are established.

∙ Third order partial differential equations

We consider now the Taylor expansion of the relation (21) at third order: $\begin{matrix} (61) & \{\begin{matrix} m + Δ t \partial_{t} m + \frac{1}{2} Δ t^{2} \partial_{t}^{2} m + \frac{1}{6} Δ t^{3} \partial_{t}^{3} m + O (Δ t^{4}) \\ = m^{} - Δ t Λ m^{} + \frac{1}{2} Δ t^{2} Λ^{2} m^{} - \frac{1}{6} Δ t^{3} Λ^{3} m^{} + O (Δ t^{4}) . \end{matrix} \end{matrix}$ We consider the first conserved components of the moments: $\begin{matrix} \{\begin{matrix} W + Δ t \partial_{t} W + \frac{1}{2} Δ t^{2} \partial_{t}^{2} W + \frac{1}{6} Δ t^{3} \partial_{t}^{3} W + O (Δ t^{4}) \\ = W - Δ t (A W + B Y^{}) + \frac{1}{2} Δ t^{2} (A_{2} W + B_{2} Y^{}) - \frac{1}{6} Δ t^{3} (A_{3} W + B_{3} Y^{}) + O (Δ t^{4}) . \end{matrix} \end{matrix}$ We simplify the constant term W and divide by $Δ t$ . We deduce $\begin{matrix} (62) & \{\begin{matrix} \partial_{t} W = - A W - B Y^{} + \frac{1}{2} Δ t (A_{2} W + B_{2} Y^{} - \partial_{t}^{2} W) \\ - \frac{1}{6} Δ t^{2} (A_{3} W + B_{3} Y^{} + \partial_{t}^{3} W) + O (Δ t^{3}) . \end{matrix} \end{matrix}$ We explicit the partial derivatives $\partial_{t}^{2} W$ and $\partial_{t}^{3} W$ at orders one and zero respectively. From the relation $\begin{matrix} \partial_{t} W + Γ_{1} + Δ t Γ_{2} = O (Δ t^{2}) \end{matrix}$ we have $\begin{array}{l} \partial_{t}^{2} W & = \partial_{t} (- Γ_{1} - Δ t Γ_{2} + O (Δ t^{2})) \\ = d (- Γ_{1} - Δ t Γ_{2} + O (Δ t^{2})) . (\partial_{t} W) \\ = d (- Γ_{1} - Δ t Γ_{2} + O (Δ t^{2})) . (- Γ_{1} - Δ t Γ_{2} + O (Δ t^{2}) \end{array}$ and $\begin{matrix} (63) & \partial_{t}^{2} W = d Γ_{1} . Γ_{1} + Δ t (d Γ_{1} . Γ_{2} + d Γ_{2} . Γ_{1}) + O (Δ t^{2}) . \end{matrix}$ Then $\begin{array}{l} \partial_{t}^{3} W & = \partial_{t} (\partial_{t}^{2} W) \\ = \partial_{t} (d Γ_{1} . Γ_{1}) + O (Δ t) \\ = d (d Γ_{1} . Γ_{1}) . \partial_{t} W + O (Δ t) \\ = - d (d Γ_{1} . Γ_{1}) . Γ_{1} + O (Δ t) \\ (64) & \partial_{t}^{3} W & = - \partial^{2} Γ_{1} . Γ_{1} + O (Δ t) . \end{array}$ In order to compute the coefficient $Γ_{3}$ , we precise now the expansions of the quantities $A_{j} W + B_{j} Y^{}$ for $j = 1$ , 2 and 3. We have $\begin{array}{l} A W + B Y^{} \\ = A W + B (Φ (W) + (Σ - \frac{1}{2} I) (Δ t Ψ_{1} + Δ t^{2} Ψ_{2}) + O (Δ t^{3}) \\ = A W + B Φ + Δ t B (Σ - \frac{1}{2} I) Ψ_{1} + Δ t^{2} B (Σ - \frac{1}{2} I) Ψ_{2} + O (Δ t^{3}) \\ A_{2} W + B_{2} Y^{} \\ = A_{2} W + B_{2} (Φ (W) + (Σ - \frac{1}{2} I) Δ t Ψ_{1}) + O (Δ t^{2}) \\ = (A^{2} + B C) W + (A B + B D) Φ + Δ t B_{2} (Σ - \frac{1}{2} I) Ψ_{1} + O (Δ t^{2}) \\ = A Γ_{1} + B (γ_{1} - Ψ_{1}) + Δ t B_{2} (Σ - \frac{1}{2} I) Ψ_{1} + O (Δ t^{2}) because C W + D Ψ_{1} = γ_{1} - Ψ_{1} \\ = d Γ_{1} . Γ_{1} - B Ψ_{1} + Δ t B_{2} (Σ - \frac{1}{2} I) Ψ_{1} + O (Δ t^{2}) because d Γ_{1} . Γ_{1} = A Γ_{1} + B γ_{1}, \\ A_{3} W + B_{3} Y^{} \\ = A_{3} W + B_{3} Φ (W) + O (Δ t) \\ = (A_{2} A + B_{2} C) W + (A_{2} B + B_{2} D) Φ + O (Δ t) \\ = A_{2} Γ_{1} + B_{2} (γ_{1} - Ψ_{1}) + O (Δ t) \\ = (A^{2} + B C) Γ_{1} + (A B + B D) γ_{1} - B_{2} Ψ_{1} + O (Δ t) \\ = A (A Γ_{1} + B γ_{1}) + B (C Γ_{1} + D γ_{1}) - B_{2} Ψ_{1} + O (Δ t) \\ = A d Γ_{1} . Γ_{1} + B (C Γ_{1} + D γ_{1}) - B_{2} Ψ_{1} + O (Δ t) because d Γ_{1} . Γ_{1} = A Γ_{1} + B γ_{1} \\ = A d Γ_{1} . Γ_{1} + B d (γ_{1} - Ψ_{1}) . Γ_{1} - B_{2} Ψ_{1} + O (Δ t) because C Γ_{1} + D γ_{1} = d (γ_{1} - Ψ_{1}) . Γ_{1} \\ = \partial^{2} Γ_{1} . Γ_{1} - B d Ψ_{1} . Γ_{1} - B_{2} Ψ_{1} + O (Δ t) \\ because \partial^{2} Γ_{1} . Γ_{1} \equiv d (d Γ_{1} . Γ_{1}) . Γ_{1} = A d Γ_{1} . Γ_{1} + B d γ_{1} . Γ_{1} . \end{array}$

We deduce now from (62), (63), (64) and the previous expressions: $\begin{array}{l} \partial_{t} W = & - A W - B Y^{} + \frac{1}{2} Δ t (A_{2} W + B_{2} Y^{} - \partial_{t}^{2} W) - \frac{1}{6} Δ t^{2} (A_{3} W + B_{3} Y^{*} + \partial_{t}^{3} W) + O (Δ t^{3}) \\ = & - (A W + B Φ + Δ t B (Σ - \frac{1}{2} I) Ψ_{1} + Δ t^{2} B (Σ - \frac{1}{2} I) Ψ_{2}) \\ + \frac{1}{2} Δ t (d Γ_{1} . Γ_{1} - B Ψ_{1} + Δ t B_{2} (Σ - \frac{1}{2} I) Ψ_{1} - d Γ_{1} . Γ_{1} - Δ t (d Γ_{1} . Γ_{2} + d Γ_{2} . Γ_{1})) \\ - \frac{1}{6} Δ t^{2} (\partial^{2} Γ_{1} . Γ_{1} - B d Ψ_{1} . Γ_{1} - B_{2} Ψ_{1} - \partial^{2} Γ_{1} . Γ_{1}) + O (Δ t^{3}) \\ = & - Γ_{1} - Δ t B Σ Ψ_{1} + Δ t^{2} (- B Σ Ψ_{2} + \frac{1}{2} B Ψ_{2} + \frac{1}{2} B_{2} Σ Ψ_{1} - \frac{1}{4} B_{2} Ψ_{1} - \frac{1}{2} (d Γ_{1} . Γ_{2} + d Γ_{2} . Γ_{1}) \\ + \frac{1}{6} B d Ψ_{1} . Γ_{1} + \frac{1}{6} B_{2} Ψ_{1}) + O (Δ t^{3}) because \partial^{2} Γ_{1} . Γ_{1} = A d Γ_{1} . Γ_{1} + B d γ_{1} . Γ_{1} \\ = & - Γ_{1} - Δ t B Σ Ψ_{1} + Δ t^{2} (- B Σ Ψ_{2} + \frac{1}{2} B (Σ d Ψ_{1} . Γ_{1} + γ_{2} - D Σ Ψ_{1}) + \frac{1}{2} (A B + B D) Σ Ψ_{1} \\ - \frac{1}{12} B_{2} Ψ_{1} - \frac{1}{2} (A Γ_{2} + B γ_{2}) - \frac{1}{2} B Σ d Ψ_{1} . Γ_{1} + \frac{1}{6} B d Ψ_{1} . Γ_{1}) + O (Δ t^{3}) \\ = & - Γ_{1} - Δ t Γ_{2} - Δ t^{2} (B Σ Ψ_{2} + \frac{1}{12} B_{2} Ψ_{1} - \frac{1}{6} B d Ψ_{1} . Γ_{1}) + O (Δ t^{3}) \end{array}$ and the relations (54) and (55) are established. □

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 $\begin{matrix} M_{0 j} = 1, M_{α j} = v_{j}^{α}, 0 ⩽ j < q, 1 ⩽ α ⩽ d . \end{matrix}$ The first moments are denoted by $ρ \equiv \sum_{j} f_{j}$ and $J_{α} \equiv \sum_{j} v_{j}^{α} f_{j}$ . For the scalar advection diffusion ( $N = 1$ ), the third order equivalent partial differential equation ( 54 ) ( 55 ) can be written $\begin{matrix} (65) & \partial_{t} ρ + \partial_{α} J_{α}^{eq} - Δ t {\tilde{σ}}_{α} \partial_{α} {\tilde{θ}}_{α} + Δ t^{2} [Λ_{β ℓ}^{γ} ({\tilde{σ}}_{β} {\tilde{σ}}_{ℓ} - \frac{1}{12}) \partial_{β} \partial_{γ} {\tilde{θ}}_{ℓ} + ({\tilde{σ}}_{β}^{2} - \frac{1}{6}) \partial_{t} (\partial_{β} {\tilde{θ}}_{β})] = O (Δ t^{3}) \end{matrix}$ with $J_{α}^{eq} \equiv Φ_{α} (ρ)$ and the defects of conservation ${\tilde{θ}}_{k} (W)$ defined in ( 47 ). For the fluid system ( $N = d + 1$ ), we have $\begin{matrix} (66) & \{\begin{matrix} \partial_{t} ρ + \partial_{α} J_{α} - \frac{1}{12} Δ t^{2} Λ_{β ℓ}^{γ} \partial_{β} \partial_{γ} {\tilde{θ}}_{ℓ} = O (Δ t^{3}) \\ \partial_{t} J_{α} + Λ_{α k}^{β} \partial_{β} m_{k}^{eq} - Δ t Λ_{α k}^{β} {\tilde{σ}}_{k} \partial_{β} {\tilde{θ}}_{k} \\ + Δ t^{2} [Λ_{α k}^{β} Λ_{k ℓ}^{γ} ({\tilde{σ}}_{k} {\tilde{σ}}_{ℓ} - \frac{1}{12}) \partial_{β} \partial_{γ} {\tilde{θ}}_{ℓ} + Λ_{α k}^{β} ({\tilde{σ}}_{k}^{2} - \frac{1}{6}) \partial_{β} \partial_{t} {\tilde{θ}}_{k}] = O (Δ t^{3}) . \end{matrix} \end{matrix}$ The relations ( 65 ) and ( 66 ) are axactly the ones proposed with numbers (35), (40) and (41) in the reference [ 17 ].

Proof of Proposition 8.
We first consider the conservation defect ${\tilde{θ}}_{k} (W) \equiv \partial_{t} m_{k}^{eq} + \sum_{ℓ β} {\tilde{Λ}}_{k ℓ}^{β} \partial_{β} m_{ℓ}^{eq}$ introduced in (47). We define a vector conservation defect θ of dimension $q - N$ just obtained by a shift of component numbering: $\begin{matrix} θ_{k} = {\tilde{θ}}_{k + N} (W), 0 ⩽ k < q - N . \end{matrix}$ We have for $0 ⩽ k < q - N$ : $\begin{array}{l} θ_{k} & = {(d Φ . \partial_{t} W)}_{k} + {(\tilde{Λ} . m^{eq})}_{k + N} \\ = {(d Φ . \partial_{t} W)}_{k} + {(C W + D Φ (W))}_{k} \\ = {(d Φ . (- Γ_{1} - Δ t Γ_{2}))}_{k} + {(C W + D Φ (W))}_{k} + O (Δ t^{2}) due to the expansion (29) \\ = - {(Ψ_{1} + Δ t d Φ . Γ_{2})}_{k} + O (Δ t^{2}) due to the relation (31) \end{array}$ and $\begin{matrix} θ = - Ψ_{1} - Δ t d Φ . Γ_{2} + O (Δ t^{2}) . \end{matrix}$ In particular, $\partial_{t} θ = - \partial_{t} Ψ_{1} + O (Δ t) = - d Ψ_{1} . (- Γ_{1}) + O (Δ t)$ and $\begin{matrix} \partial_{t} θ = d Ψ_{1} . Γ_{1} + O (Δ t) . \end{matrix}$ Due to (41) and (55), we have the following calculus: $\begin{array}{l} Γ_{2} + Δ t Γ_{3} = & B Σ Ψ_{1} + Δ t (B Σ Ψ_{2} - \frac{1}{12} B_{2} Ψ_{1} + \frac{1}{6} B d Ψ_{1} . Γ_{1}) \\ = & B Σ (- θ + Δ t d Φ . Γ_{2} + O (Δ t^{2})) + Δ t B Σ (D Σ Ψ_{1} - d Φ . Γ_{2} - Σ d Ψ_{1} . Γ_{1}) \\ - \frac{1}{12} Δ t B_{2} Ψ_{1} + \frac{1}{6} Δ t B d Ψ_{1} . Γ_{1} \\ = & - B Σ θ + Δ t [(B Σ D Σ - \frac{1}{12} B_{2}) Ψ_{1} + B (\frac{1}{6} - Σ^{2}) d Ψ_{1} . Γ_{1}] + O (Δ t^{2}) \\ = & - B Σ θ + Δ t [- (B Σ D Σ - \frac{1}{12} B_{2}) θ + B (Σ^{2} - \frac{1}{6}) \partial_{t} θ] + O (Δ t^{2}) . \end{array}$ The third-order partial equivalent equations (54) can be written under the form $\begin{matrix} \partial_{t} W + \tilde{Λ} m^{eq} + Δ t B Σ θ + Δ t^{2} [(B Σ D Σ - \frac{1}{12} B_{2}) θ - B (Σ^{2} - \frac{1}{6}) \partial_{t} θ] = O (Δ t^{3}) . \end{matrix}$

∙ We explicit now the following cartesian components: $\begin{array}{l} {(B Σ D Σ θ)}_{i} = - {\tilde{Λ}}_{i k}^{β} {\tilde{σ}}_{k} (- {\tilde{Λ}}_{k ℓ}^{γ}) {\tilde{σ}}_{ℓ} \partial_{β} \partial_{γ} {\tilde{θ}}_{ℓ} = {\tilde{Λ}}_{i k}^{β} {\tilde{σ}}_{k} {\tilde{Λ}}_{k ℓ}^{γ} {\tilde{σ}}_{ℓ} \partial_{β} \partial_{γ} {\tilde{θ}}_{ℓ}, \\ {(B_{2} θ)}_{i} = ({\tilde{Λ}}_{i k}^{β} \partial_{β}) ({\tilde{Λ}}_{k ℓ}^{γ} \partial_{γ}) {\tilde{θ}}_{ℓ} = {\tilde{Λ}}_{i k}^{β} {\tilde{Λ}}_{k ℓ}^{γ} \partial_{β} \partial_{γ} {\tilde{θ}}_{ℓ}, \\ {((B Σ D Σ - \frac{1}{12} B_{2}) \tilde{θ})}_{i} = {\tilde{Λ}}_{i k}^{β} {\tilde{Λ}}_{k ℓ}^{γ} ({\tilde{σ}}_{k} {\tilde{σ}}_{ℓ} - \frac{1}{12}) \partial_{β} \partial_{γ} {\tilde{θ}}_{ℓ} \end{array}$ and $\begin{array}{l} - {(B (Σ^{2} - \frac{1}{6}) \partial_{t} \tilde{θ})}_{i} = {\tilde{Λ}}_{i k}^{β} ({\tilde{σ}}_{k}^{2} - \frac{1}{6}) \partial_{β} \partial_{t} {\tilde{θ}}_{k} . \end{array}$

We deduce a new form of the third-order partial equivalent equations for $0 ⩽ i < N$ : $\begin{matrix} (67) & \{\begin{matrix} \partial_{t} W_{i} + {\tilde{Λ}}_{i k}^{α} \partial_{α} m_{k}^{eq} - Δ t {\tilde{Λ}}_{i k}^{β} {\tilde{σ}}_{k} \partial_{β} {\tilde{θ}}_{k} \\ + Δ t^{2} [{\tilde{Λ}}_{i k}^{β} {\tilde{Λ}}_{k ℓ}^{γ} ({\tilde{σ}}_{k} {\tilde{σ}}_{ℓ} - \frac{1}{12}) \partial_{β} \partial_{γ} {\tilde{θ}}_{ℓ} + {\tilde{Λ}}_{i k}^{β} ({\tilde{σ}}_{k}^{2} - \frac{1}{6}) \partial_{β} \partial_{t} {\tilde{θ}}_{k}] = O (Δ t^{3}) . \end{matrix} \end{matrix}$

∙ We focus on the mass conservation ( $i = 0$ ). We have $\begin{array}{l} {\tilde{Λ}}_{0 ℓ}^{α} = M_{0 j} v_{j}^{α} {(M^{- 1})}_{j ℓ} = M_{α j} {(M^{- 1})}_{j ℓ} = δ_{α ℓ} . \end{array}$ Then $\begin{array}{l} {\tilde{Λ}}_{0 k}^{β} {\tilde{σ}}_{k} \partial_{β} {\tilde{θ}}_{k} = δ_{β k} {\tilde{σ}}_{k} \partial_{β} {\tilde{θ}}_{k} = {\tilde{σ}}_{β} \partial_{β} {\tilde{θ}}_{β} . \end{array}$

In the thermal case, we have only one conservation law and ${\tilde{σ}}_{β} \neq 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, ${\tilde{σ}}_{β} = 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 $\begin{array}{l} {\tilde{Λ}}_{0 k}^{β} {\tilde{Λ}}_{k ℓ}^{γ} ({\tilde{σ}}_{k} {\tilde{σ}}_{ℓ} - \frac{1}{12}) \partial_{β} \partial_{γ} θ_{ℓ} = δ_{β k} {\tilde{Λ}}_{k ℓ}^{γ} ({\tilde{σ}}_{k} {\tilde{σ}}_{ℓ} - \frac{1}{12}) \partial_{β} \partial_{γ} θ_{ℓ} = {\tilde{Λ}}_{β ℓ}^{γ} ({\tilde{σ}}_{β} {\tilde{σ}}_{ℓ} - \frac{1}{12}) \partial_{β} \partial_{γ} θ_{ℓ} . \end{array}$

In the termal case, ${\tilde{σ}}_{β} \neq 0$ and this term is exactly the first of the two third order terms of the relation (65). In the fluid case, ${\tilde{σ}}_{β} = 0$ and this term has a more compact expression. Finally, $\begin{matrix} (68) & {\tilde{Λ}}_{0 k}^{β} {\tilde{Λ}}_{k ℓ}^{γ} ({\tilde{σ}}_{k} {\tilde{σ}}_{ℓ} - \frac{1}{12}) \partial_{β} \partial_{γ} θ_{ℓ} = \{\begin{matrix} {\tilde{Λ}}_{β ℓ}^{γ} ({\tilde{σ}}_{β} {\tilde{σ}}_{ℓ} - \frac{1}{12}) \partial_{β} \partial_{γ} θ_{ℓ} & thermics \\ - \frac{1}{12} {\tilde{Λ}}_{β ℓ}^{γ} \partial_{β} \partial_{γ} θ_{ℓ} & fluid . \end{matrix} \end{matrix}$

For the second term relative to third order, we have $\begin{array}{l} {\tilde{Λ}}_{0 k}^{β} ({\tilde{σ}}_{k}^{2} - \frac{1}{6}) \partial_{β} \partial_{t} θ_{k} = ({\tilde{σ}}_{β}^{2} - \frac{1}{6}) \partial_{t} (\partial_{β} θ_{β}) \end{array}$

In the thermal case, ${\tilde{σ}}_{β} \neq 0$ and ${\tilde{Λ}}_{0 k}^{β} ({\tilde{σ}}_{k}^{2} - \frac{1}{6}) \partial_{β} \partial_{t} θ_{k} = ({\tilde{σ}}_{β}^{2} - \frac{1}{6}) \partial_{t} (\partial_{β} θ_{β})$ . In the tluid case, the momentum is conserved and $θ_{β} = O (Δ t)$ and ${\tilde{Λ}}_{0 k}^{β} ({\tilde{σ}}_{k}^{2} - \frac{1}{6}) \partial_{β} \partial_{t} θ_{k} = O (Δ t)$ . In a synthetic way, $\begin{matrix} (69) & {\tilde{Λ}}_{0 k}^{β} ({\tilde{σ}}_{k}^{2} - \frac{1}{6}) \partial_{β} \partial_{t} θ_{k} = \{\begin{matrix} ({\tilde{σ}}_{β}^{2} - \frac{1}{6}) \partial_{t} (\partial_{β} θ_{β}) & thermics \\ O (Δ t) & fluid . \end{matrix} \end{matrix}$ We deduce from (68) and (69) the expansion (65) in the thermal case and the first relation of (66) in the fluid case.

∙ For the momentum equation of the fluid case, $1 ⩽ i = α ⩽ d$ and the relation (67) can be rewritten as $\begin{matrix} \{\begin{matrix} \partial_{t} J_{α} + {\tilde{Λ}}_{α k}^{α} \partial_{α} m_{k}^{eq} - Δ t {\tilde{Λ}}_{α k}^{β} {\tilde{σ}}_{k} \partial_{β} {\tilde{θ}}_{k} \\ + Δ t^{2} [{\tilde{Λ}}_{α k}^{β} {\tilde{Λ}}_{k ℓ}^{γ} ({\tilde{σ}}_{k} {\tilde{σ}}_{ℓ} - \frac{1}{12}) \partial_{β} \partial_{γ} {\tilde{θ}}_{ℓ} + {\tilde{Λ}}_{α k}^{β} ({\tilde{σ}}_{k}^{2} - \frac{1}{6}) \partial_{β} \partial_{t} {\tilde{θ}}_{k}] = O (Δ t^{3}) . \end{matrix} \end{matrix}$ This relation is exactly the second relation of (66). □
Proposition 9 (Nonlinear expansion at fourth order).

When the time step $Δ t$ has an infinitesimal value, the expansion of the microscopic variables takes the form $\begin{matrix} (70) & Y = Φ (W) + Δ t S^{- 1} Ψ_{1} (W) + Δ t^{2} S^{- 1} Ψ_{2} (W) + Δ t^{3} S^{- 1} Ψ_{3} (W) + O (Δ t^{4}) . \end{matrix}$ The vectors $Ψ_{1} (W)$ and $Ψ_{2} (W)$ have been evaluated at the relations ( 39 ) and ( 53 ). We have $\begin{matrix} (71) & \{\begin{matrix} Ψ_{3} (W) = d Φ (W) . Γ_{3} (W) + Σ d Ψ_{1} (W) . Γ_{2} (W) + Σ d Ψ_{2} (W) . Γ_{1} (W) - D Σ Ψ_{2} (W) \\ + \frac{1}{6} D d Ψ_{1} (W) . Γ_{1} (W) - \frac{1}{12} D_{2} Ψ_{1} (W) - \frac{1}{12} \partial^{2} Ψ_{1} (W) . Γ_{1} (W) . \end{matrix} \end{matrix}$ A system of fourth order partial differential equations is asymptotically satisfied by the lattice Boltzmann scheme: $\begin{matrix} (72) & \partial_{t} W + Γ_{1} (W) + Δ t Γ_{2} (W) + Δ t^{2} Γ_{3} (W) + Δ t^{3} Γ_{4} (W) = O (Δ t^{4}) . \end{matrix}$ The vectors $Γ_{1} (W)$ , $Γ_{2} (W)$ and $Γ_{3} (W)$ have been precised at the relations ( 36 ), ( 41 ) and ( 55 ) respectively. We have now $\begin{matrix} (73) & \{\begin{matrix} Γ_{4} (W) = B Σ Ψ_{3} (W) + \frac{1}{4} B_{2} Ψ_{2} (W) + \frac{1}{6} B D_{2} Σ Ψ_{1} - \frac{1}{6} A B Ψ_{2} (W) \\ - \frac{1}{6} B d γ_{1} (W) . Γ_{2} (W) - \frac{1}{6} B d γ_{2} (W) . Γ_{1} (W) - \frac{1}{6} B Σ \partial^{2} Ψ_{1} (W) . Γ_{1} (W) . \end{matrix} \end{matrix}$

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 $\begin{array}{l} \{\begin{matrix} Y + Δ t \partial_{t} Y + \frac{1}{2} Δ t^{2} \partial_{t}^{2} Y + \frac{1}{6} Δ t^{3} \partial_{t}^{3} Y + O (Δ t^{4}) \\ = Y^{} - Δ t (C W + D Y^{}) + \frac{1}{2} Δ t^{2} (C_{2} W + D_{2} Y^{}) - \frac{1}{6} Δ t^{3} (C_{3} W + D_{3} Y^{}) + O (Δ t^{4}) . \end{matrix} \end{array}$

The relaxation process can be written $Y - Y^{} \equiv S (Y - Φ (W))$ . We deduce $\begin{matrix} (74) & \{\begin{matrix} S (Y - Φ (W)) = - Δ t (C W + D Y^{} + \partial_{t} Y) + \frac{1}{2} Δ t^{2} (C_{2} W + D_{2} Y^{} - \partial_{t}^{2} Y) \\ - \frac{1}{6} Δ t^{3} (C_{3} W + D_{3} Y^{} + \partial_{t}^{3} Y) + O (Δ t^{4}) . \end{matrix} \end{matrix}$ 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: $\begin{matrix} Y = Φ (W) + (Σ + \frac{1}{2} I) Δ t Ψ_{1} (W) + (Σ + \frac{1}{2} I) Δ t^{2} Ψ_{2} (W) + O (Δ t^{3}) . \end{matrix}$ We differentiate this previous relation relative to time: $\begin{array}{l} \partial_{t} Y = & d Φ . \partial_{t} W + (Σ + \frac{1}{2} I) Δ t d Ψ_{1} . \partial_{t} W + (Σ + \frac{1}{2} I) Δ t^{2} d Ψ_{2} . \partial_{t} W + O (Δ t^{3}) \\ = & - d Φ . (Γ_{1} (W) + Δ t Γ_{2} (W) + Δ t^{2} Γ_{3} (W) + O (Δ t^{3})) \\ - Δ t (Σ + \frac{1}{2} I) d Ψ_{1} . (Γ_{1} + Δ t Γ_{2} + O (Δ t^{2})) \\ + Δ t^{2} (Σ + \frac{1}{2} I) d Ψ_{2} . (- Γ_{1} + O (Δ t)) + O (Δ t^{3}) \end{array}$ and $\begin{matrix} (75) & \{\begin{matrix} \partial_{t} Y = - d Φ (W) . Γ_{1} (W) - Δ t [d Φ (W) . Γ_{2} (W) + (Σ + \frac{1}{2} I) d Ψ_{1} (W) . Γ_{1} (W)] \\ - Δ t^{2} [d Φ (W) . Γ_{3} (W) + (Σ + \frac{1}{2} I) d Ψ_{1} (W) . Γ_{2} (W) \\ + (Σ + \frac{1}{2} I) d Ψ_{2} (W) . Γ_{1} (W)] + O (Δ t^{3}) . \end{matrix} \end{matrix}$ 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 $\begin{array}{l} \partial_{t}^{2} Y = - d γ_{1} . (- Γ_{1} - Δ t Γ_{2}) - Δ t [- d γ_{2} . Γ_{1} - (Σ + \frac{1}{2} I) \partial^{2} Φ (W) . Γ_{1} (W)] + O (Δ t^{2}) \end{array}$ and $\begin{matrix} (76) & \partial_{t}^{2} Y = d γ_{1} . Γ_{1} + Δ t [d γ_{1} . Γ_{2} + d γ_{2} . Γ_{1} + (Σ + \frac{1}{2} I) \partial^{2} Ψ_{1} (W) . Γ_{1} (W)] + O (Δ t^{2}) . \end{matrix}$ Finally, $\begin{matrix} (77) & \partial_{t}^{3} Y = - \partial^{2} γ_{1} (W) . Γ_{1} (W) + O (Δ t) . \end{matrix}$ 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 $\begin{array}{l} C W + D Y^{} = C W + D [Φ (W) + (Σ - \frac{1}{2} I) Δ t Ψ_{1} + (Σ - \frac{1}{2} I) Δ t^{2} Ψ_{2}] + O (Δ t^{2}) \end{array}$ and $\begin{matrix} (78) & \{\begin{matrix} C W + D Y^{} = γ_{1} (W) - Ψ_{1} (W) + Δ t D (Σ - \frac{1}{2} I) Ψ_{1} (W) \\ + Δ t^{2} D (Σ - \frac{1}{2} I) Ψ_{2} (W) + O (Δ t^{3}) . \end{matrix} \end{matrix}$ For the second order term, $\begin{array}{l} C_{2} W + D_{2} Y^{} = C_{2} W + D_{2} [Φ (W) + Δ t (Σ - \frac{1}{2} I) Ψ_{1}] + O (Δ t^{2}) \end{array}$ and due to (60), $\begin{matrix} (79) & C_{2} W + D_{2} Y^{} = d (γ_{1} - Ψ_{1}) . Γ_{1} (W) - D Ψ_{1} (W) + Δ t D_{2} (Σ - \frac{1}{2} I) Ψ_{1} + O (Δ t^{2}) . \end{matrix}$ We have finally for the third order terms $\begin{array}{l} C_{3} W + D_{3} Y^{} & = C_{3} W + D_{3} Φ (W) + O (Δ t) \\ = (C_{2} A + D_{2} C) W + (C_{2} B + D_{2} D) Φ (W) + O (Δ t) \\ = C_{2} Γ_{1} + D_{2} (γ_{1} - Ψ_{1}) + O (Δ t) \\ = (C A + D C) Γ_{1} + (C B + D^{2}) γ_{1} - D_{2} Ψ_{1} + O (Δ t) \\ = C d γ_{1} . Γ_{1} + D d (γ_{1} - Ψ_{1}) . Γ_{1} - D_{2} Ψ_{1} + O (Δ t) \\ = C d γ_{1} . Γ_{1} + D d γ_{1} . Γ_{1} - D d Ψ_{1} . Γ_{1} - D_{2} Ψ_{1} + O (Δ t) \end{array}$ and $\begin{matrix} (80) & C_{3} W + D_{3} Y^{} = \partial^{2} (γ_{1} - Ψ_{1}) . Γ_{1} - D d Ψ_{1} . Γ_{1} - D_{2} Ψ_{1} + O (Δ t) \end{matrix}$ Then due to (74), (75), (76), (77), (78), (79) and (80), we have $\begin{array}{l} S (Y - Φ (W)) \\ = - Δ t (C W + D Y^{} + \partial_{t} Y) + \frac{1}{2} Δ t^{2} (C_{2} W + D_{2} Y^{} - \partial_{t}^{2} Y) \\ - \frac{1}{6} Δ t^{3} (C_{3} W + D_{3} Y^{} + \partial_{t}^{3} Y) + O (Δ t^{4}) \\ = - Δ t (γ_{1} - Ψ_{1} + Δ t D (Σ - \frac{1}{2} I) Ψ_{1} + Δ t^{2} D (Σ - \frac{1}{2} I) Ψ_{2}) \\ - Δ t [- d Φ . Γ_{1} - Δ t (d Φ . Γ_{2} + (Σ + \frac{1}{2} I) d Ψ_{1} . Γ_{1}) \\ - Δ t^{2} (d Φ . Γ_{3} + (Σ + \frac{1}{2} I) d Ψ_{1} . Γ_{2} + (Σ + \frac{1}{2} I) d Ψ_{2} . Γ_{1})] \\ + \frac{1}{2} Δ t^{2} (d (γ_{1} - Ψ_{1}) . Γ_{1} - D Ψ_{1} + Δ t D_{2} (Σ - \frac{1}{2} I) Ψ_{1}) \\ - \frac{1}{2} Δ t^{2} (d γ_{1} . Γ_{1} + Δ t [d γ_{1} . Γ_{2} + d γ_{2} . Γ_{1} + (Σ + \frac{1}{2} I) \partial^{2} Ψ_{1} . Γ_{1}]) \\ - \frac{1}{6} Δ t^{3} (\partial^{2} (γ_{1} - Ψ_{1}) . Γ_{1} - D d Ψ_{1} . Γ_{1} - D_{2} Ψ_{1}) - \frac{1}{6} Δ t^{3} (- \partial^{2} γ_{1} . Γ_{1}) + O (Δ t^{4}) \\ = Δ t Ψ_{1} + Δ t^{2} Ψ_{2} + Δ t^{3} [- D (Σ - \frac{1}{2} I) Ψ_{2} + d Φ . Γ_{3} + (Σ + \frac{1}{2} I) d Ψ_{1} . Γ_{2} \\ + (Σ + \frac{1}{2} I) d Ψ_{2} . Γ_{1} + \frac{1}{2} D_{2} (Σ - \frac{1}{2} I) Ψ_{1} - \frac{1}{2} d γ_{1} . Γ_{2} - \frac{1}{2} d γ_{2} . Γ_{1} - \frac{1}{2} (Σ + \frac{1}{2} I) \partial^{2} Ψ_{1} . Γ_{1} \\ + \frac{1}{6} \partial^{2} Ψ_{1} . Γ_{1} + \frac{1}{6} D d Ψ_{1} . Γ_{1} + \frac{1}{6} D_{2} Ψ_{1}] + O (Δ t^{4}) \\ = Δ t Ψ_{1} + Δ t^{2} Ψ_{2} + Δ t^{3} [- D Σ Ψ_{2} + \frac{1}{2} D (Σ d Ψ_{1} . Γ_{1} + γ_{2} - D Σ Ψ_{1}) + d Φ . Γ_{3} + Σ d Ψ_{1} . Γ_{2} \\ + \frac{1}{2} d Ψ_{1} . Γ_{2} + Σ d Ψ_{2} . Γ_{1} + \frac{1}{2} (Σ \partial^{2} Ψ_{1} . Γ_{1} + d γ_{2} . Γ_{1} - D Σ d Ψ_{1} . Γ_{1}) + \frac{1}{2} (C B + D^{2}) Σ Ψ_{1} \\ - \frac{1}{4} D_{2} Ψ_{1} - \frac{1}{2} d γ_{1} . Γ_{2} - \frac{1}{2} d γ_{2} . Γ_{1} - \frac{1}{2} Σ \partial^{2} Ψ_{1} . Γ_{1} - \frac{1}{12} \partial^{2} Ψ_{1} . Γ_{1} + \frac{1}{6} D d Ψ_{1} . Γ_{1} + \frac{1}{6} D_{2} Ψ_{1}] \\ + O (Δ t^{4}) \\ = Δ t Ψ_{1} + Δ t^{2} Ψ_{2} + Δ t^{3} [d Φ . Γ_{3} - D Σ Ψ_{2} + \frac{1}{2} D γ_{2} + Σ d Ψ_{1} . Γ_{2} + \frac{1}{2} d Ψ_{1} . Γ_{2} + Σ d Ψ_{2} . Γ_{1} \\ + \frac{1}{2} C B Σ Ψ_{1} - \frac{1}{12} D_{2} Ψ_{1} - \frac{1}{2} d γ_{1} . Γ_{2} - \frac{1}{12} \partial^{2} Ψ_{1} . Γ_{1} + \frac{1}{6} D d Ψ_{1} . Γ_{1}] + O (Δ t^{4}) \\ = Δ t Ψ_{1} + Δ t^{2} Ψ_{2} \\ + Δ t^{3} [d Φ . Γ_{3} - D Σ Ψ_{2} + \frac{1}{2} D γ_{2} + Σ d Ψ_{1} . Γ_{2} + \frac{1}{2} (- C Γ_{2} - D γ_{2} + d γ_{1} . Γ_{2}) \\ + Σ d Ψ_{2} . Γ_{1} + \frac{1}{2} C B Σ Ψ_{1} - \frac{1}{12} D_{2} Ψ_{1} - \frac{1}{2} d γ_{1} . Γ_{2} - \frac{1}{12} \partial^{2} Ψ_{1} . Γ_{1} + \frac{1}{6} D d Ψ_{1} . Γ_{1}] + O (Δ t^{4}) \\ = Δ t Ψ_{1} + Δ t^{2} Ψ_{2} \\ + Δ t^{3} [d Φ . Γ_{3} - D Σ Ψ_{2} + Σ d Ψ_{1} . Γ_{2} + Σ d Ψ_{2} . Γ_{1} - \frac{1}{12} D_{2} Ψ_{1} - \frac{1}{12} \partial^{2} Ψ_{1} . Γ_{1} + \frac{1}{6} D d Ψ_{1} . Γ_{1}] \\ + O (Δ t^{4}) . \end{array}$ The relations (70) and (71) are established.

∙ Fourth order partial differential equations

We consider now the Taylor expansion of the relation (21) at fourth order: $\begin{array}{l} m + Δ t \partial_{t} m + \frac{1}{2} Δ t^{2} \partial_{t}^{2} m + \frac{1}{6} Δ t^{3} \partial_{t}^{3} m + \frac{1}{24} Δ t^{4} \partial_{t}^{4} m + O (Δ t^{5}) \\ = m^{} - Δ t Λ m^{} + \frac{1}{2} Δ t^{2} Λ^{2} m^{} - \frac{1}{6} Δ t^{3} Λ^{3} m^{} + \frac{1}{24} Δ t^{4} Λ^{4} m^{} + O (Δ t^{5}) . \end{array}$ For the first components (conserved moments): $\begin{array}{l} W + Δ t \partial_{t} W + \frac{1}{2} Δ t^{2} \partial_{t}^{2} W + \frac{1}{6} Δ t^{3} \partial_{t}^{3} W + \frac{1}{24} Δ t^{4} \partial_{t}^{4} W + O (Δ t^{5}) \\ = W - Δ t (A W + B Y^{}) + \frac{1}{2} Δ t^{2} (A_{2} W + B_{2} Y^{}) - \frac{1}{6} Δ t^{3} (A_{3} W + B_{3} Y^{}) \\ + \frac{1}{24} Δ t^{4} (A_{4} W + B_{4} Y^{}) + O (Δ t^{5}) . \end{array}$ We simplify the constant term W and divide by $Δ t$ . We deduce $\begin{matrix} (81) & \{\begin{matrix} \partial_{t} W = - (A W + B Y^{}) + \frac{1}{2} Δ t (A_{2} W + B_{2} Y^{} - \partial_{t}^{2} W) \\ - \frac{1}{6} Δ t^{2} (A_{3} W + B_{3} Y^{} + \partial_{t}^{3} W) + \frac{1}{24} Δ t^{3} (A_{4} W + B_{4} Y^{} - \partial_{t}^{4} W) + O (Δ t^{4}) . \end{matrix} \end{matrix}$ 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 $\begin{matrix} \partial_{t} W + Γ_{1} + Δ t Γ_{2} + Δ t^{2} Γ_{3} = O (Δ t^{3}) \end{matrix}$ where $Γ_{1}$ , $Γ_{2}$ and $Γ_{3}$ have been evaluated previously. We have $\begin{array}{l} \partial_{t}^{2} W & = \partial_{t} (- Γ_{1} - Δ t Γ_{2} - Δ t^{2} Γ_{3} + O (Δ t^{3})) \\ = d (- Γ_{1} - Δ t Γ_{2} - Δ t^{2} Γ_{3} + O (Δ t^{3})) . (\partial_{t} W) \\ = d (- Γ_{1} - Δ t Γ_{2} - Δ t^{2} Γ_{3} + O (Δ t^{2})) . (- Γ_{1} - Δ t Γ_{2} - Δ t^{2} Γ_{3}) + O (Δ t^{3}) \end{array}$ and $\begin{matrix} (82) & \partial_{t}^{2} W = d Γ_{1} . Γ_{1} + Δ t (d Γ_{1} . Γ_{2} + d Γ_{2} . Γ_{1}) + Δ t^{2} (d Γ_{1} . Γ_{3} + d Γ_{2} . Γ_{2} + d Γ_{3} . Γ_{1}) + O (Δ t^{2}) . \end{matrix}$ Then $\begin{array}{l} \partial_{t}^{3} W = & \partial_{t} (\partial_{t}^{2} W) \\ = & \partial_{t} (A Γ_{1} + B d Φ . Γ_{1} + Δ t \partial_{t} (d Γ_{2} . Γ_{1} + d Γ_{1} . Γ_{2})) + O (Δ t^{2}) \\ = & (A d Γ_{1} + B d (d Φ . Γ_{1})) . (- Γ_{1} - Δ t Γ_{2} + O (Δ t^{2})) \\ + Δ t \partial_{t} (B Σ d Ψ_{1} . Γ_{1} + A Γ_{2} + B d Φ . Γ_{2}) + O (Δ t^{2}) \\ = & \partial_{t} (d Γ_{1} . Γ_{1} + Δ t (d Γ_{1} . Γ_{2} + d Γ_{2} . Γ_{1})) + O (Δ t^{2}) \\ = & d (d Γ_{1} . Γ_{1} + Δ t (d Γ_{1} . Γ_{2} + d Γ_{2} . Γ_{1})) . \partial_{t} W + O (Δ t^{2}) \\ = & d (d Γ_{1} . Γ_{1} + Δ t (d Γ_{1} . Γ_{2} + d Γ_{2} . Γ_{1})) . (- Γ_{1} - Δ t Γ_{2}) + O (Δ t^{2}) \\ = & - d (d Γ_{1} . Γ_{1}) . Γ_{1} - Δ t [d (A Γ_{2} + B γ_{2} + B Σ d Ψ_{1} . Γ_{1}) . Γ_{1} + d (A Γ_{1} + B γ_{1}) . Γ_{2}] + O (Δ t^{2}) \\ = & - \partial^{2} Γ_{1} . Γ_{1} - Δ t [A d Γ_{2} . Γ_{1} + B d γ_{2} . Γ_{1} + B Σ \partial^{2} Ψ_{1} . Γ_{1} + A d Γ_{1} . Γ_{2} + B d γ_{1} . Γ_{2}] + O (Δ t^{2}) \end{array}$ and we have at first order for this third order derivative $\begin{matrix} (83) & \{\begin{matrix} \partial_{t}^{3} W = - \partial^{2} Γ_{1} . Γ_{1} - Δ t [A (d Γ_{2} . Γ_{1} + d Γ_{1} . Γ_{2}) \\ + B (d γ_{2} . Γ_{1} + d γ_{1} . Γ_{2} + Σ \partial^{2} Ψ_{1} . Γ_{1})] + O (Δ t^{2}) . \end{matrix} \end{matrix}$ Last but not least, we have $\begin{array}{l} \partial_{t}^{4} W & = \partial_{t} (- \partial^{2} Γ_{1} . Γ_{1} + O (Δ t)) \\ = - \partial_{t} (A d Γ_{1} . Γ_{1} + B d γ_{1} . Γ_{1}) + O (Δ t)) \\ = d (A d Γ_{1} . Γ_{1} + B d γ_{1} . Γ_{1}) . (Γ_{1} + O (Δ t)) \end{array}$ and finally $\begin{matrix} (84) & \partial_{t}^{4} W = A \partial^{2} Γ_{1} . Γ_{1} + B \partial^{2} γ_{1} . Γ_{1} + O (Δ t) . \end{matrix}$ For the determination of the coefficient $Γ_{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: $\begin{array}{l} A W + B Y^{} \\ = A W + B [Φ (W) + (Σ - \frac{1}{2} I) (Δ t Ψ_{1} + Δ t^{2} Ψ_{2} + Δ t^{3} Ψ_{3})] + O (Δ t^{4}) \\ = A W + B Φ + Δ t B (Σ - \frac{1}{2} I) Ψ_{1} + Δ t^{2} B (Σ - \frac{1}{2} I) Ψ_{2} + Δ t^{3} B (Σ - \frac{1}{2} I) Ψ_{3} + O (Δ t^{4}) \\ A_{2} W + B_{2} Y^{} \\ = A_{2} W + B_{2} [Φ (W) + (Σ - \frac{1}{2} I) (Δ t Ψ_{1} + Δ t^{2} Ψ_{2})] + O (Δ t^{3}) \\ = (A^{2} + B C) W + (A B + B D) Φ + Δ t B_{2} (Σ - \frac{1}{2} I) Ψ_{1} + Δ t^{2} B_{2} (Σ - \frac{1}{2} I) Ψ_{2} + O (Δ t^{3}) \\ = d Γ_{1} . Γ_{1} - B Ψ_{1} + Δ t B_{2} (Σ - \frac{1}{2} I) Ψ_{1} + Δ t^{2} B_{2} (Σ - \frac{1}{2} I) Ψ_{2} + O (Δ t^{3}) \\ A_{3} W + B_{3} Y^{} \\ = A_{3} W + B_{3} [Φ (W) + Δ t (Σ - \frac{1}{2} I) Ψ_{1}] + O (Δ t^{2}) \\ = (A_{2} A + B_{2} C) W + (A_{2} B + B_{2} D) Φ + Δ t B_{3} (Σ - \frac{1}{2} I) Ψ_{1} + O (Δ t^{2}) \\ = \partial^{2} Γ_{1} . Γ_{1} - B d Ψ_{1} . Γ_{1} - B_{2} Ψ_{1} + Δ t B_{3} (Σ - \frac{1}{2} I) Ψ_{1} + O (Δ t^{2}) \\ A_{4} W + B_{4} Y^{} \\ = A_{4} W + B_{4} Φ (W) + O (Δ t) \\ = (A_{3} A + B_{3} C) W + (A_{3} B + B_{3} D) Φ + O (Δ t) \\ = A_{3} Γ_{1} + B_{3} (γ_{1} - Ψ_{1}) + O (Δ t) \\ = (A_{2} A + B_{2} C) Γ_{1} + (A_{2} B + B_{2} D) γ_{1} - B_{3} Ψ_{1} + O (Δ t) \\ = A_{2} d Γ_{1} . Γ_{1} + B_{2} d (γ_{1} - Ψ_{1}) . Γ_{1} - B_{3} Ψ_{1} + O (Δ t) \\ = (A^{2} + B C) d Γ_{1} . Γ_{1} + (A B + B D) d γ_{1} . Γ_{1} - B_{2} d Ψ_{1} . Γ_{1} - B_{3} Ψ_{1} + O (Δ t) \\ = A \partial^{2} Γ_{1} . Γ_{1} + B \partial^{2} (γ_{1} - Ψ_{1}) . Γ_{1} - B_{2} d Ψ_{1} . Γ_{1} - B_{3} Ψ_{1} + O (Δ t) \\ = A \partial^{2} Γ_{1} . Γ_{1} + B \partial^{2} γ_{1} . Γ_{1} - B \partial^{2} Ψ_{1} . Γ_{1} - B_{2} d Ψ_{1} . Γ_{1} - B_{3} Ψ_{1} + O (Δ t) . \end{array}$ We deduce now from (81), (82), (83), (84) and the previous expressions: $\begin{array}{l} \partial_{t} W = & - (A W + B Y^{}) + \frac{1}{2} Δ t (A_{2} W + B_{2} Y^{} - \partial_{t}^{2} W) - \frac{1}{6} Δ t^{2} (A_{3} W + B_{3} Y^{} + \partial_{t}^{3} W) \\ + \frac{1}{24} Δ t^{3} (A_{4} W + B_{4} Y^{} - \partial_{t}^{4} W) + O (Δ t^{4}) \\ = & - (A W + B Φ + Δ t B (Σ - \frac{1}{2} I) Ψ_{1} + Δ t^{2} B (Σ - \frac{1}{2} I) Ψ_{2} + Δ t^{3} B (Σ - \frac{1}{2} I) Ψ_{3}) \\ + \frac{1}{2} Δ t [d Γ_{1} . Γ_{1} - B Ψ_{1} + Δ t B_{2} (Σ - \frac{1}{2} I) Ψ_{1} + Δ t^{2} B_{2} (Σ - \frac{1}{2} I) Ψ_{2} \\ - (d Γ_{1} . Γ_{1} + Δ t (d Γ_{1} . Γ_{2} + d Γ_{2} . Γ_{1}) + Δ t^{2} (d Γ_{1} . Γ_{3} + d Γ_{2} . Γ_{2} + d Γ_{3} . Γ_{1}))] \\ - \frac{1}{6} Δ t^{2} [\partial^{2} Γ_{1} . Γ_{1} - B d Ψ_{1} . Γ_{1} - B_{2} Ψ_{1} + Δ t B_{3} (Σ - \frac{1}{2} I) Ψ_{1} \\ - \partial^{2} Γ_{1} . Γ_{1} - Δ t (A (d Γ_{2} . Γ_{1} + d Γ_{1} . Γ_{2}) + B (d γ_{2} . Γ_{1} + d γ_{1} . Γ_{2} + Σ \partial^{2} Ψ_{1} . Γ_{1}))] \\ + \frac{1}{24} Δ t^{3} [A \partial^{2} Γ_{1} . Γ_{1} + B \partial^{2} γ_{1} . Γ_{1} - B \partial^{2} Ψ_{1} . Γ_{1} - B_{2} d Ψ_{1} . Γ_{1} - B_{3} Ψ_{1} \\ - (A \partial^{2} Γ_{1} . Γ_{1} + B \partial^{2} γ_{1} . Γ_{1})] + O (Δ t^{4}) \\ = & - Γ_{1} - Δ t Γ_{2} - Δ t^{2} Γ_{3} + Δ t^{3} [- B (Σ - \frac{1}{2} I) Ψ_{3} + \frac{1}{2} B_{2} (Σ - \frac{1}{2} I) Ψ_{2} \\ - \frac{1}{2} (d Γ_{1} . Γ_{3} + d Γ_{2} . Γ_{2} + d Γ_{3} . Γ_{1}) + \frac{1}{6} (- B_{3} (Σ - \frac{1}{2} I) Ψ_{1} + A (d Γ_{2} . Γ_{1} + d Γ_{1} . Γ_{2}) \\ + B (d γ_{2} . Γ_{1} + d γ_{1} . Γ_{2} + Σ \partial^{2} Ψ_{1} . Γ_{1})) + \frac{1}{24} (A \partial^{2} Γ_{1} . Γ_{1} + B \partial^{2} γ_{1} . Γ_{1} \\ - B \partial^{2} Ψ_{1} . Γ_{1} - B_{2} d Ψ_{1} . Γ_{1} - B_{3} Ψ_{1} - A \partial^{2} Γ_{1} . Γ_{1} - B \partial^{2} γ_{1} . Γ_{1})] + O (Δ t^{4}) \\ = & - Γ_{1} - Δ t Γ_{2} - Δ t^{2} Γ_{3} - Δ t^{3} [B Σ Ψ_{3} - \frac{1}{2} B (γ_{3} + Σ d Ψ_{1} . Γ_{2} + Σ d Ψ_{2} . Γ_{1} - D Σ Ψ_{2} \\ + \frac{1}{6} D d Ψ_{1} . Γ_{1} - \frac{1}{12} D_{2} Ψ_{1} \\ - \frac{1}{12} \partial^{2} Ψ_{1} . Γ_{1}) - \frac{1}{2} (A B + B D) Σ Ψ_{2} + \frac{1}{4} B_{2} Ψ_{2} \\ + \frac{1}{2} (A Γ_{3} + B γ_{3} + B Σ d Ψ_{1} . Γ_{1} + d Γ_{3} . Γ_{1}) + \frac{1}{6} B_{3} Σ Ψ_{1} - \frac{1}{12} B_{3} Ψ_{1} - \frac{1}{6} A B Σ d Ψ_{1} . Γ_{1} \\ - \frac{1}{6} A (A Γ_{2} + B γ_{2}) - \frac{1}{6} B (d γ_{2} . Γ_{1} + d γ_{1} . Γ_{2} + Σ \partial^{2} Ψ_{1} . Γ_{1}) + \frac{1}{24} B \partial^{2} Ψ_{1} . Γ_{1} \\ + \frac{1}{24} B_{2} d Ψ_{1} . Γ_{1} + \frac{1}{24} B_{3} Ψ_{1}] + O (Δ t^{4}) \\ = & - Γ_{1} - Δ t Γ_{2} - Δ t^{2} Γ_{3} - Δ t^{3} [B Σ Ψ_{3} - \frac{1}{2} B Σ d Ψ_{2} . Γ_{1} - \frac{1}{12} B D d Ψ_{1} . Γ_{1} + \frac{1}{24} B D_{2} Ψ_{1} \\ + \frac{1}{24} B \partial^{2} Ψ_{1} . Γ_{1} - \frac{1}{2} A B Σ Ψ_{2} + \frac{1}{4} B_{2} Ψ_{2} + \frac{1}{2} A (B Σ Ψ_{2} - \frac{1}{6} B d Ψ_{1} . Γ_{1} + \frac{1}{12} B_{2} Ψ_{1}) \\ + \frac{1}{2} (B Σ d Ψ_{2} . Γ_{1} - \frac{1}{6} B \partial^{2} Ψ_{1} . Γ_{1} + \frac{1}{12} B_{2} d Ψ_{1} . Γ_{1}) + \frac{1}{6} B_{3} Σ Ψ_{1} - \frac{1}{24} B_{3} Ψ_{1} \\ - \frac{1}{6} A B Σ d Ψ_{1} . Γ_{1} \\ - \frac{1}{6} A^{2} Γ_{2} - \frac{1}{6} A B γ_{2} - \frac{1}{6} B (d γ_{2} . Γ_{1} + d γ_{1} . Γ_{2} + Σ \partial^{2} Ψ_{1} . Γ_{1}) + \frac{1}{24} B \partial^{2} Ψ_{1} . Γ_{1} \\ + \frac{1}{24} B_{2} d Ψ_{1} . Γ_{1}] + O (Δ t^{4}) \\ = & - Γ_{1} - Δ t Γ_{2} - Δ t^{2} Γ_{3} - Δ t^{3} [B Σ Ψ_{3} - \frac{1}{12} B_{2} d Ψ_{1} . Γ_{1} \\ + \frac{1}{24} B D_{2} Ψ_{1} + \frac{1}{4} B_{2} Ψ_{2} + \frac{1}{24} A B_{2} Ψ_{1} \\ + \frac{1}{6} B_{3} Σ Ψ_{1} - \frac{1}{24} B_{3} Ψ_{1} - \frac{1}{6} A B Σ d Ψ_{1} . Γ_{1} - \frac{1}{6} A^{2} B Σ Ψ_{1} - \frac{1}{6} A B γ_{2} \\ - \frac{1}{6} B (d γ_{2} . Γ_{1} + d γ_{1} . Γ_{2} + Σ \partial^{2} Ψ_{1} . Γ_{1}) + \frac{1}{12} B_{2} d Ψ_{1} . Γ_{1}] + O (Δ t^{4}) \\ = & - Γ_{1} - Δ t Γ_{2} - Δ t^{2} Γ_{3} - Δ t^{3} [B Σ Ψ_{3} + \frac{1}{4} B_{2} Ψ_{2} + \frac{1}{6} B_{3} Σ Ψ_{1} - \frac{1}{6} A B Σ d Ψ_{1} . Γ_{1} \\ - \frac{1}{6} A^{2} B Σ Ψ_{1} - \frac{1}{6} A B (Ψ_{2} - Σ d Ψ_{1} . Γ_{1} + D Σ Ψ_{1}) - \frac{1}{6} B (d γ_{2} . Γ_{1} + d γ_{1} . Γ_{2} + Σ \partial^{2} Ψ_{1} . Γ_{1})] \\ + O (Δ t^{4}) because A B_{2} + B D_{2} = B_{3} \\ = & - Γ_{1} - Δ t Γ_{2} - Δ t^{2} Γ_{3} - Δ t^{3} [B Σ Ψ_{3} + \frac{1}{4} B_{2} Ψ_{2} + \frac{1}{6} B D_{2} Σ Ψ_{1} - \frac{1}{6} A B Ψ_{2} \\ - \frac{1}{6} B (d γ_{2} . Γ_{1} + d γ_{1} . Γ_{2} + Σ \partial^{2} Ψ_{1} . Γ_{1})] + O (Δ t^{4}) \end{array}$ because $B_{3} - A (A B + B D) = B D_{2}$ . The relations (72) and (73) are established and the proposition is proven. □

7. Revisiting the “Berlin algorithm” in the linear case

In [3], 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: $\begin{matrix} Φ (W) = E W \end{matrix}$ with a given fixed $(q - N) \times N$ rectangle matrix E. Then the nonequilibium moments can be expanded as $\begin{matrix} Y = [E + S^{- 1} (Δ t β_{1} + Δ t^{2} β_{2} + Δ t^{3} β_{3})] W + O (Δ t^{4}) . \end{matrix}$ The $(q - N) \times N$ matrices $β_{j}$ are linear operators of order j. The equivalent partial equivalent equations are linear and can be written as $\begin{matrix} \partial_{t} W + [α_{1} + Δ t α_{2} + Δ t^{2} α_{3} + Δ t^{3} α_{4}] W + O (Δ t^{4}), \end{matrix}$ and $α_{j}$ is a $N \times N$ differential matrix of order j. In other terms, $Φ_{j} (W) = β_{j} W$ and $Γ_{j} (W) = α_{j} W$ are linear functions of the conserved moments. The coefficent matrices $α_{j}$ and $β_{j}$ are computed according to a linearized version of ( 31 ): $\begin{matrix} (85) & \{\begin{matrix} α_{1} = A + B E \\ β_{1} = E α_{1} - C - D E \\ α_{2} = B Σ β_{1} \\ β_{2} = Σ β_{1} α_{1} + E α_{2} - D Σ β_{1} \\ α_{3} = B Σ β_{2} + \frac{1}{12} B_{2} β_{1} - \frac{1}{6} B β_{1} α_{1} \\ β_{3} = Σ β_{1} α_{2} + E α_{3} - D Σ β_{2} + Σ β_{2} α_{1} + \frac{1}{6} D β_{1} α_{1} - \frac{1}{12} D_{2} β_{1} - \frac{1}{12} β_{1} α_{1}^{2} \\ α_{4} = B Σ β_{3} + \frac{1}{4} B_{2} β_{2} + \frac{1}{6} B D_{2} Σ β_{1} - \frac{1}{6} A B β_{2} - \frac{1}{6} B E α_{1} α_{2} - \frac{1}{6} B E α_{2} α_{1} - \frac{1}{6} B Σ β_{1} α_{1}^{2} . \end{matrix} \end{matrix}$

Proof of Proposition 10.
First observe that $d Φ . ξ = E ξ$ for any test vector ξ. Due to the relation (36), we have $Γ_{1} = A W + B Φ (W) = A W + B E W = (A + B E) W$ and $α_{1} = A + B E$ . Due to (39), we have moreover $\begin{array}{l} Ψ_{1} & = d Φ (W) . Γ_{1} (W) - (C W + D Φ (W)) \\ = E α_{1} W - C W - D E W \\ = (E α_{1} - C - D E) W \end{array}$ and the second relation of (85) is established. Due to (41), $\begin{array}{l} Γ_{2} (W) & = B Σ Ψ_{1} (W) \\ = B Σ β_{1} W . \end{array}$ Then $α_{2} = B Σ β_{1}$ and the third line of (85) is proven. We have now from (53), $\begin{array}{l} Ψ_{2} (W) & = Σ d Ψ_{1} (W) . Γ_{1} (W) + d Φ (W) . Γ_{2} (W) - D Σ Ψ_{1} (W) \\ = Σ β_{1} α_{1} W + E α_{2} W - D Σ β_{1} W \\ = (Σ β_{1} α_{1} + E α_{2} - D Σ β_{1}) W \end{array}$ and $β_{2} = Σ β_{1} α_{1} + E α_{2} - D Σ β_{1}$ as suggested in the fourth line of(85). From the relation (55), we have now $\begin{array}{l} Γ_{3} (W) & = B Σ Ψ_{2} (W) + \frac{1}{12} B_{2} Ψ_{1} (W) - \frac{1}{6} B d Ψ_{1} (W) . Γ_{1} (W) \\ = B Σ β_{2} W + \frac{1}{12} B_{2} β_{1} W - \frac{1}{6} B β_{1} . α_{1} W \\ = (B Σ β_{2} + \frac{1}{12} B_{2} β_{1} - \frac{1}{6} B β_{1} . α_{1}) W . \end{array}$ Then $α_{3} = B Σ β_{2} + \frac{1}{12} B_{2} β_{1} - \frac{1}{6} B β_{1} . α_{1}$ as proposed in the fifth line of (85). Observe now that due to (28), we have $\begin{array}{l} \partial^{2} Ψ_{1} . Γ_{1} (W) & = d^{2} Ψ_{1} (W) . (Γ_{1}, Γ_{1}) + d Ψ_{1} (W) . d Γ_{1} (W) . Γ_{1} \\ = d Ψ_{1} (W) . d Γ_{1} (W) . Γ_{1} because Ψ_{1} is linear \\ = β_{1} α_{1}^{2} W . \end{array}$ Then due to (71), we have in this linear case, $\begin{array}{l} Ψ_{3} (W) = & d Φ (W) . Γ_{3} (W) + Σ d Ψ_{1} (W) . Γ_{2} (W) + Σ d Ψ_{2} (W) . Γ_{1} (W) - D Σ Ψ_{2} (W) \\ + \frac{1}{6} D d Ψ_{1} (W) . Γ_{1} (W) - \frac{1}{12} D_{2} Ψ_{1} (W) - \frac{1}{12} \partial^{2} Ψ_{1} (W) . Γ_{1} (W) \\ = & E α_{3} W + Σ β_{1} α_{2} W + Σ β_{2} α_{1} W - D Σ β_{2} W + \frac{1}{6} D β_{1} α_{1} W \\ - \frac{1}{12} D_{2} β_{1} W - \frac{1}{12} β_{1} α_{1}^{2} W \\ = & (E α_{3} + Σ β_{1} α_{2} + Σ β_{2} α_{1} - D Σ β_{2} + \frac{1}{6} D β_{1} α_{1} - \frac{1}{12} D_{2} β_{1} - \frac{1}{12} β_{1} α_{1}^{2}) W \end{array}$ and $β_{3} = E α_{3} + Σ β_{1} α_{2} + Σ β_{2} α_{1} - D Σ β_{2} + \frac{1}{6} D β_{1} α_{1} - \frac{1}{12} D_{2} β_{1} - \frac{1}{12} β_{1} α_{1}^{2}$ . We have finally, due to (73), $\begin{array}{l} Γ_{4} (W) = & B Σ Ψ_{3} (W) + \frac{1}{4} B_{2} Ψ_{2} (W) + \frac{1}{6} B D_{2} Σ Ψ_{1} - \frac{1}{6} A B Ψ_{2} (W) - \frac{1}{6} B d γ_{1} (W) . Γ_{2} (W) \\ - \frac{1}{6} B d γ_{2} (W) . Γ_{1} (W) - \frac{1}{6} B Σ \partial^{2} Ψ_{1} (W) . Γ_{1} (W) \\ = & B Σ β_{3} W + \frac{1}{4} B_{2} β_{2} W + \frac{1}{6} B D_{2} Σ β_{1} W - \frac{1}{6} A B β_{2} W - \frac{1}{6} B E α_{2} W \\ - \frac{1}{6} B E α_{1} W - \frac{1}{6} B Σ β_{1} α_{1}^{2} W \\ = & (B Σ β_{3} + \frac{1}{4} B_{2} β_{2} + \frac{1}{6} B D_{2} Σ β_{1} - \frac{1}{6} A B β_{2} - \frac{1}{6} B E α_{2} α_{1} - \frac{1}{6} B E α_{1} α_{2} - \frac{1}{6} B Σ β_{1} α_{1}^{2}) W \end{array}$ and the proposition is established. □

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 [19,29], regularized lattice Boltzmann [44,48,53] or entropic lattice Boltzmann schemes [39] are also into study. After a first step in [47], the link between the single scale Taylor expansion and the multiple scale Chapman–Enskog is also a project for the near future.

Footnotes

Acknowledgements

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 [].

References

Agarwal,

K.Y.

Yun and

Balakrishnan, Beyond Navier Stokes – Burnett equations for flow simulations in continuum-transition regime, Physics of Fluids 13 (2001), 3061–3085. doi:10.1063/1.1397256.

Alexander,

Chen and

Sterling, Lattice Boltzmann thermohydrodynamic, Physical Review E 47 (1993), R2249–R2252.

Augier,

Dubois,

Graille and

Lallemand, On rotational invariance of lattice Boltzmann schemes, Computers and Mathematics with Applications 67 (2014), 239–255. doi:10.1016/j.camwa.2013.06.009.

Bhatnagar,

Gross and

Krook, A model for collision processes in gases I. Small amplitude processes in charged and neutral one-component systems, Physical Review 94 (1954), 511–525. doi:10.1103/PhysRev.94.511.

Boghosian and

Coveney, Inverse Chapman–Enskog derivation of the thermohydrodynamic lattice-BGK model for the ideal gas, International Journal of Modern Physics C 9 (1998), 1231–1245. doi:10.1142/S0129183198001114.

Broadwell, Study of rarefied shear flow by the discrete velocity method, Journal of Fluid Mechanics 19 (1964), 401–414. doi:10.1017/S0022112064000817.

Carleman, Problèmes mathématiques dans la théorie cinétique des gaz, Publications Scientifiques de l’Institut Mittag-Leffler, Almqvist & Wiksells Boktryckeri AB, Uppsala University, 1957.

Chapman and

T.G.

Cooling, The Mathematical Theory of Non-uniform Gases, Cambridge University Press, 1939.

Chen and

G.D.

Doolen, Lattice Boltzmann method for fluid flows, Annual Review of Fluid Mechanics 30 (1998), 329–364. doi:10.1146/annurev.fluid.30.1.329.

10.

Dellar, Non-hydrodynamic modes and a priori construction of shallow water lattice Boltzmann equations, Physical Review E 65 (2002), 036309. doi:10.1103/PhysRevE.65.036309.

11.

Dellar, An interpretation and derivation of the lattice Boltzmann method using Strang splitting, Computers & Mathematics with Applications 65 (2013), 129–141. doi:10.1016/j.camwa.2011.08.047.

12.

Dellar, Lattice Boltzmann algorithms without cubic defects in Galilean invariance on standard lattices, Journal of Computational Physics 259 (2014), 270–283. doi:10.1016/j.jcp.2013.11.021.

13.

d’Humières, Generalized lattice–Boltzmann equations, in: Rarefied Gas Dynamics: Theory and Simulations, AIAA Progress in Astronautics and Astronautics, Vol. 159, 1992, pp. 450–458.

14.

Drui,

Franck,

Helluy and

Navoret, An analysis of over-relaxation in a kinetic approximation of systems of conservation laws, Comptes Rendus Mécanique 347 (2019), 259–269. doi:10.1016/j.crme.2018.12.001.

15.

Dubois, Une introduction au schéma de Boltzmann sur réseau, ESAIM Proceedings and Surveys 18 (2007), 181–215. doi:10.1051/proc:071815.

16.

Dubois, Equivalent partial differential equations of a Boltzmann scheme, Computers and mathematics with applications 55 (2008), 1441–1449. doi:10.1016/j.camwa.2007.08.003.

17.

Dubois, Third order equivalent equation of lattice Boltzmann scheme, Discrete and Continuous Dynamical Systems-Series A 23(1/2) (2009), 221–248.

18.

Dubois, Simulation of strong nonlinear waves with vectorial lattice Boltzmann schemes, International Journal of Modern Physics C 25(12) (2014), 1441014. doi:10.1142/S0129183114410149.

19.

Dubois,

Février and

Graille, Lattice Boltzmann schemes with relative velocities, Communications in Computational Physics 17(4) (2015), 1088–1112. doi:10.4208/cicp.2014.m394.

20.

Dubois,

Février and

Graille, Third order equivalent equation for the relative velocity lattice Boltzmann schemes with one conservation law, arxiv.org/abs/1502.02143, 7 February 2015.

21.

Dubois,

Graille and

Lallemand, Recovering the full Navier Stokes equations with lattice Boltzmann schemes, in: 30th International Symposium on Rarefied Gas Dynamics, Victoria, BC, Canada, Jul 2016, Vol. 1786, American Institute of Physics Proceedings, 2016, pp. 40003–40003.

22.

Dubois and

Lallemand, Towards higher order lattice Boltzmann schemes, Journal of Statistical Mechanics: Theory and Experiment P (2009), 06006.

23.

Dubois and

Lallemand, Quartic parameters for acoustic applications of lattice Boltzmann scheme, Computers and Mathematics with Applications 61 (2011), 3404–3416. doi:10.1016/j.camwa.2011.01.011.

24.

Dubois and

Lallemand, On triangular lattice Boltzmann schemes for scalar problems, Communications in Computational Physics 13 (2013), 649–670. doi:10.4208/cicp.381011.270112s.

25.

Dubois and

Lallemand, Comparison of simulations of convective flows, Communications in Computational Physics 17 (2015), 1169–1184. doi:10.4208/cicp.2014.m400.

26.

Dubois,

Lallemand and

M.M.

Tekitek, Taylor expansion method for linear lattice Boltzmann schemes with an external drift. Application to boundary conditions, in: High Order Nonlinear Numerical Schemes for Evolutionary PDEs,

Abgrall et al., eds, Springer Lecture Notes in Computational Science and Engineering, Vol. 99, 2014, pp. 89–107. doi:10.1007/978-3-319-05455-1_6.

27.

Fletcher, Semi-Lagrangian Advection Methods and Their Applications in Geoscience, Elsevier Science Publishing Co, United States, 2019.

28.

Gatignol, Théorie cinétique des gaz à répartition discrète de vitesses, Springer Lecture Notes in Physics 36 (1975).

29.

M.C.

Geier,

Greiner and

J.G.

Korvink, Cascaded digital lattice Boltzmann automata for high Reynolds number flow, Physical Review E 73 (2006), 066705. doi:10.1103/PhysRevE.73.066705.

30.

Graille, Approximation of mono-dimensional hyperbolic systems: A lattice Boltzmann scheme as a relaxation method, Journal of Computational Physics 266 (2014), 74–88. doi:10.1016/j.jcp.2014.02.017.

31.

Graille and

Gouarin, pyLBM, see www.imo.universite-paris-saclay.fr/~graille/pylbm.php and pypi.org/project/pylbm, 2017.

32.

Gross, Recent investigations of the Boltzmann equation, in: Proceedings of the First International Symposium Held in Nice,

F.M.

Devienne, ed., Rarefied Gas Dynamics, Vol. 3, Pergamon Press, London, 1960.

33.

Guo and

Shu, Lattice Boltzmann Method and Its Applications in Engineering, World Scientific, 2013.

34.

Guo,

Zheng and

Shi, Thermal lattice Boltzmann equation for low Mach number flows: Decoupling model, Physical Review E 75 (2007), 036704. doi:10.1103/PhysRevE.75.036704.

35.

Házi and

Kávrán, On the cubic velocity deviations in lattice Boltzmann methods, Journal of Physics A: Mathematical and General 39(12) (2006), 3127–3136. doi:10.1088/0305-4470/39/12/019.

36.

Hénon, Viscosity of a lattice gas, Complex Systems 1 (1987), 763–789.

37.

F.J.

Higuera,

Succi and

Benzi, Lattice gas dynamics with enhanced collisions, Europhysics Letters 9(4) (1991), 345–349. doi:10.1209/0295-5075/9/4/008.

38.

Huang, Statistical Mechanics, John Wiley & Sons Inc, New York, 1963.

39.

I.V.

Karlin,

A.N.

Gorban,

Succi and

Boffi, Maximum entropy principle for lattice kinetic equations, Physical Review Letters 81 (1998), 6–9. doi:10.1103/PhysRevLett.81.6.

40.

I.V.

Karlin,

Sichau and

S.S.

Chikatamarla, Consistent two-population lattice Boltzmann model for thermal flows, Physical Review E 88 (2013), 063310. doi:10.1103/PhysRevE.88.063310.

41.

Krüger,

Kusumaatmaja,

Kuzmin,

Shardt,

Silva and

E.M.

Viggen, The Lattice Boltzmann Method: Principles and Practice, Springer Verlag, 2017.

42.

Lallemand and

Dubois, Comparison of simulations of convective flows, Communications in Computational Physics 17 (2015), 1169–1184. doi:10.4208/cicp.2014.m400.

43.

Lallemand and

L.-S.

Luo, Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, Galilean invariance, and stability, Physical Review E 61 (2000), 6546–6562. doi:10.1103/PhysRevE.61.6546.

44.

Latt and

Chopard, Lattice Boltzmann method with regularized pre-collision distribution functions, Mathematics and Computers in Simulation 72 (2006), 165–168. doi:10.1016/j.matcom.2006.05.017.

45.

Lerat and

Peyret, Noncentered schemes and shock propagation problems, Computers and Fluids 2 (1974), 35–52. doi:10.1016/0045-7930(74)90004-8.

46.

Otomo,

B.M.

Boghosian and

Dubois, Two complementary lattice-Boltzmann-based analyses for nonlinear systems, Physica A: Statistical Mechanics and its Applications 486 (2017), 1000–1011. doi:10.1016/j.physa.2017.06.010.

47.

Otomo,

B.M.

Boghosian and

Dubois, Efficient lattice Boltzmann models for the Kuramoto–Sivashinsky equation, Computers and Fluids 172 (2018), 683–688. doi:10.1016/j.compfluid.2018.01.036.

48.

Philippi and

Hegele, From the continuous to the lattice Boltzmann equation: The discretization problem and thermal models, Physical Review E 73(5) (2006), 056702. doi:10.1103/PhysRevE.73.056702.

49.

N.I.

Prasianakis,

Karlin,

Mantzaras and

Boulouchos, Lattice Boltzmann method with restored Galilean invariance, Physical Review E 79 (2009), 066702. doi:10.1103/PhysRevE.79.066702.

50.

Y.-H.

Qian,

d’Humières and

Lallemand, Lattice BGK models for Navier–Stokes equation, Europhysics Letters 17(6) (1992), 479–484. doi:10.1209/0295-5075/17/6/001.

51.

Y.H.

Qian and

Zhou, Higher-order dynamics in lattice-based models using the Chapman–Enskog method, Physical Review E 61(2) (2000), 2103–2106. doi:10.1103/PhysRevE.61.2103.

52.

SageMath, the Sage Mathematics Software System (Version 7.5.1), The Sage Developers, http://www.sagemath.org, 2017.

53.

Shan,

X.F.

Yuan and

Chen, Kinetic theory representation of hydrodynamics: A way beyond the Navier–Stokes equation, Journal of Fluid Mechanics 550 (2006), 413–441. doi:10.1017/S0022112005008153.

54.

Strang, On the construction and comparison of difference schemes, SIAM Journal on Numerical Analysis 5 (1968), 506–517. doi:10.1137/0705041.

55.

Struchtrup, Failures of the Burnett and super-Burnett equations in steady state processes, Continuum Mechanics and Thermodynamics 17 (2005), 43–50. doi:10.1007/s00161-004-0186-0.

56.

Succi, Lattice Boltzmann 2038, Europhysics Letters 109 (2015), 50001. doi:10.1209/0295-5075/109/50001.

57.

R.F.

Warming and

B.J.

Hyett, The modified equation approach to the stability and accuracy analysis of finite difference methods, Journal of Computational Physics 14 (1974), 159–179. doi:10.1016/0021-9991(74)90011-4.

58.

W.-A.

Yong,

Zhao and

L.-S.

Luo, Theory of the lattice Boltzmann method: Derivation of macroscopic equations via the Maxwell iteration, Physical Review E 93 (2016), 033310. doi:10.1103/PhysRevE.93.033310.