Finite Difference formulation of any lattice Boltzmann scheme

Abstract

Lattice Boltzmann schemes rely on the enlargement of the size of the target problem in order to solve PDEs in a highly parallelizable and efficient kinetic-like fashion, split into a collision and a stream phase. This structure, despite the well-known advantages from a computational standpoint, is not suitable to construct a rigorous notion of consistency with respect to the target equations and to provide a precise notion of stability. In order to alleviate these shortages and introduce a rigorous framework, we demonstrate that any lattice Boltzmann scheme can be rewritten as a corresponding multi-step Finite Difference scheme on the conserved variables. This is achieved by devising a suitable formalism based on operators, commutative algebra and polynomials. Therefore, the notion of consistency of the corresponding Finite Difference scheme allows to invoke the Lax-Richtmyer theorem in the case of linear lattice Boltzmann schemes. Moreover, we show that the frequently-used von Neumann-like stability analysis for lattice Boltzmann schemes entirely corresponds to the von Neumann stability analysis of their Finite Difference counterpart. More generally, the usual tools for the analysis of Finite Difference schemes are now readily available to study lattice Boltzmann schemes. Their relevance is verified by means of numerical illustrations.

Appendices

Appendix A: A different reduction strategy: mathematical framework

According to the discussion presented in Sect. 5.3 through an example and again in the case $N= 1$, we define the polynomial annihilating all the first row of the matrix $\varvec{A}$, except the very first element.

Definition 8

We call $\nu _{\varvec{A}} \in {\mathcal {D}}^{d}_{\Delta x}[X]$ “minimal polynomial annihilating most of the first row” (MPAMFR) of $\varvec{A}$ the monic polynomial of minimal degree under the form

$$\begin{aligned} \nu _{\varvec{A}} = X^{\text {deg}(\nu _{\varvec{A}})} + \psi _{\text {deg}(\nu _{\varvec{A}})-1} X^{\text {deg}(\nu _{\varvec{A}})-1} + \dots + \psi _{1} X+ \psi _0, \end{aligned}$$

with $(\psi _{k})_{k= 0}^{k= \text {deg}(\nu _{\varvec{A}})} \subset {\mathcal {D}}^{d}_{\Delta x}$ such that for every $j\in \llbracket 2, q \rrbracket$

$$\begin{aligned} (\varvec{A}^{\text {deg}(\nu _{\varvec{A}})})_{1j} + \psi _{\text {deg}(\nu _{\varvec{A}})-1} (\varvec{A}^{\text {deg}(\nu _{\varvec{A}})-1})_{1j} + \dots + \psi _1 (\varvec{A})_{1j} = 0. \end{aligned}$$

By seeing the coefficients of this unknown polynomial as the unknowns of a linear system, the problem of finding $\nu _{\varvec{A}}$ can be rewritten in terms of matrices.¹³ Let $K \in \llbracket 1, \text {deg}(\mu _{\varvec{A}}) \rrbracket$ and construct the matrix of variable size

$$\begin{aligned} \varvec{V}_{K}&= \begin{pmatrix} (\varvec{A})_{12} &{} \cdots &{} (\varvec{A}^{K-1})_{12} \\ \vdots &{} &{} \vdots \\ (\varvec{A})_{1, Q+1} &{} \cdots &{} (\varvec{A}^{K-1})_{1, Q+1} \\ \end{pmatrix}\in {\mathcal {M}}_{Q\times (K - 1)}({\mathcal {D}}^{d}_{\Delta x}), \\ \varvec{r}_K&= \begin{pmatrix} -(\varvec{A}^K)_{12} \\ \vdots \\ -(\varvec{A}^K)_{1, Q+1} \end{pmatrix} \in ({\mathcal {D}}^{d}_{\Delta x})^{Q}. \nonumber \end{aligned}$$

(A1)

Therefore, we want to find the smallest $K \in \llbracket 1, \text {deg}(\mu _{\varvec{A}}) \rrbracket$ such that $\varvec{V}_{K} (\psi _1, \dots , \psi _{K-1})^{\intercal } = \varvec{r}_K$ On the other hand, it should be observed that the zero order coefficient $\psi _0$ remains free. This underdetermination comes from the fact that we do not request that $\nu _{\varvec{A}}$ annihilates the whole first row.

Proposition 14

Let $N= 1$, then the lattice Boltzmann scheme (6) corresponds to a multi-step explicit Finite Difference scheme on the conserved moment $m_1$ under the form

$$\begin{aligned} {m}_1^{{n} + 1} =&- \sum _{k= 1}^{\text {deg}(\nu _{\varvec{A}}) - 1} \psi _{k} {m}_1^{{n} +1 - \text {deg}(\nu _{\varvec{A}}) + k} + \left( \sum _{k= 1}^{\text {deg}(\nu _{\varvec{A}})} \psi _{k} (\varvec{A}^{k})_{11} \right) m_1^{n+1- \text {deg}(\nu _{\varvec{A}})} \nonumber \\&+\left( \sum _{k= 0}^{\text {deg}(\nu _{\varvec{A}}) - 1} \left( \sum _{\ell = 0}^{k} \psi _{\text {deg}(\nu _{\varvec{A}}) + \ell - k} \varvec{A}^{\ell } \right) \varvec{B}\varvec{m}^{\text {eq}} \mid ^{n- k} \right) _1, \end{aligned}$$

(A2)

where $(\psi _{k})_{k= 1}^{k= \text {deg}(\nu _{\varvec{A}})} \subset {\mathcal {D}}^{d}_{\Delta x}$ are the coefficients of $\nu _{\varvec{A}} = \sum _{k= 0}^{k= \text {deg}(\nu _{\varvec{A}})} \psi _{k} X^{k}$.

Proof

By the choice of polynomial, we have that

$$\begin{aligned} \left( \sum _{k= 0}^{\text {deg}(\nu _{\varvec{A}})} \psi _{k} \varvec{A}^{k} \right) _{1\cdot } = \left( \psi _0 + \sum _{k= 1}^{\text {deg}(\nu _{\varvec{A}})} \psi _{k} (\varvec{A}^k)_{11}, 0, \dots , 0 \right) . \end{aligned}$$

Restarting from the proof of Proposition 4, we have

$$\begin{aligned}&\sum _{k= 0}^{\text {deg}(\nu _{\varvec{A}})} \psi _{k} {m}_1^{{\tilde{n}} + k}= {m}_1^{{\tilde{n}} + \text {deg}(\nu _{\varvec{A}})} + \sum _{k= 1}^{\text {deg}(\nu _{\varvec{A}}) - 1} \psi _{k} {m}_1^{{\tilde{n}} + k} + \psi _0 {m}_1^{{\tilde{n}}}, \\&\quad = \left( \left( \sum _{k= 0}^{\text {deg}(\nu _{\varvec{A}})} \psi _{k} \varvec{A}^{k} \right) \varvec{m}^{{\tilde{n}}} \right) _1 + \sum _{k= 1}^{\text {deg}(\nu _{\varvec{A}})} \psi _{k} \left( \sum _{\ell = 0}^{k-1} \varvec{A}^{\ell } \varvec{B}\varvec{m}^{\text {eq}} \mid ^{{\tilde{n}} + k- 1 - \ell } \right) _1, \\&= \psi _0 m^{{\tilde{n}}} + \left( \sum _{k= 1}^{\text {deg}(\nu _{\varvec{A}})} \psi _{k} (\varvec{A}^{k})_{11} \right) m_1^{{\tilde{n}}} + \sum _{k= 1}^{\text {deg}(\nu _{\varvec{A}})} \psi _{k} \left( \sum _{\ell = 0}^{k-1} \varvec{A}^{\ell } \varvec{B}\varvec{m}^{\text {eq}} \mid ^{{\tilde{n}} + k- 1 - \ell } \right) _1, \end{aligned}$$

therefore

$$\begin{aligned} {m}_1^{{\tilde{n}} + \text {deg}(\nu _{\varvec{A}})}&= - \sum _{k= 1}^{\text {deg}(\nu _{\varvec{A}}) - 1} \psi _{k} {m}_1^{{\tilde{n}} + k} + \left( \sum _{k= 1}^{\text {deg}(\nu _{\varvec{A}})} \psi _{k} (\varvec{A}^{k})_{11} \right) m_1^{{\tilde{n}}} \nonumber \\&\quad \,\, + \sum _{k= 1}^{\text {deg}(\nu _{\varvec{A}})} \psi _{k} \left( \sum _{\ell = 0}^{k-1} \varvec{A}^{\ell } \varvec{B}\varvec{m}^{\text {eq}} \mid ^{{\tilde{n}} + k- 1 - \ell } \right) _1. \end{aligned}$$

(A3)

Performing the usual change of indices yields the result. $\square$

Looking at Eq. (A2), we see that we do not need the value of $\psi _0$ to reduce the scheme, neither to reduce $\varvec{A}$ nor to deal with the equilibria through $\varvec{B}$. Changing time indices and putting everything on the left hand side

$$\begin{aligned} \sum _{k= 1}^{\text {deg}(\nu _{\varvec{A}})} \psi _{k} {m}_1^{{\tilde{n}} + k}&- \left( \sum _{k= 1}^{\text {deg}(\nu _{\varvec{A}})} \psi _{k} (\varvec{A}^k)_{11} \right) m^{{\tilde{n}}} = \sum _{k= 0}^{\text {deg}(\nu _{\varvec{A}})} {\tilde{\psi }}_{k} {m}_1^{{\tilde{n}} + k}, \\&= \sum _{k= 1}^{\text {deg}(\nu _{\varvec{A}})} \psi _{k} \left( \sum _{\ell = 0}^{k-1} \varvec{A}^{\ell } \varvec{B}\varvec{m}^{\text {eq}} \mid ^{{\tilde{n}} + k- 1 - \ell } \right) _1, \\&= \sum _{k= 1}^{\text {deg}(\nu _{\varvec{A}})} {\tilde{\psi }}_{k} \left( \sum _{\ell = 0}^{k-1} \varvec{A}^{\ell } \varvec{B}\varvec{m}^{\text {eq}} \mid ^{{\tilde{n}} + k- 1 - \ell } \right) _1, \end{aligned}$$

where we have defined

$$\begin{aligned} {\tilde{\psi }}_{k} = {\left\{ \begin{array}{ll} \psi _{k}, \qquad &{}k\in \llbracket 1, \text {deg}(\nu _{\varvec{A}}) \rrbracket , \\ -\sum _{\ell = 1}^{\ell = \text {deg}(\nu _{\varvec{A}})} \psi _{\ell } (\varvec{A}^{\ell })_{11}, \qquad &{}k= 0. \end{array}\right. } \end{aligned}$$

This generates a polynomial, which is indeed $\nu _{\varvec{A}}$ but with a precise choice of $\psi _0$. We will soon give a precise characterization of this particular polynomial.

Definition 9

We call ${\tilde{\nu }}_{\varvec{A}} \in {\mathcal {D}}^{d}_{\Delta x}[X]$ “minimal polynomial annihilating the first row” (MPAFR) of $\varvec{A}$ the monic polynomial of minimal degree under the form

$$\begin{aligned} {\tilde{\nu }}_{\varvec{A}} = X^{\text {deg}({\tilde{\nu }}_{\varvec{A}})} + {\tilde{\psi }}_{\text {deg}({\tilde{\nu }}_{\varvec{A}})-1} X^{\text {deg}({\tilde{\nu }}_{\varvec{A}})-1} + \dots + {\tilde{\psi }}_{1} X^{1} + {\tilde{\psi }}_0, \end{aligned}$$

such that for every $j\in \llbracket 1, q \rrbracket$

$$\begin{aligned} (\varvec{A}^{\text {deg}({\tilde{\nu }}_{\varvec{A}})})_{1j} + {\tilde{\psi }}_{\text {deg}({\tilde{\nu }}_{\varvec{A}})-1} (\varvec{A}^{\text {deg}({\tilde{\nu }}_{\varvec{A}})-1})_{1j} + \dots + {\tilde{\psi }}_1 (\varvec{A})_{1j} + {\tilde{\psi }}_0 = 0. \end{aligned}$$

(A4)

Compared to Definition 8, we are just asking the property to hold also for the very first element of the first row, namely for $j= 1$. This polynomial is $\nu _{\varvec{A}}$ for a particular choice of $\psi _0$. It has been deduced from the reduction of the lattice Boltzmann scheme.

Lemma 15

The polynomial of degree $\text {deg}(\nu _{\varvec{A}})$ given by

$$\begin{aligned} {\tilde{\nu }}_{\varvec{A}} = X^{\text {deg}(\nu _{\varvec{A}})} + \psi _{\text {deg}(\nu _{\varvec{A}})-1} X^{\text {deg}(\nu _{\varvec{A}})-1} + \dots + \psi _1 X-\sum _{l = 1}^{\text {deg}(\nu _{\varvec{A}})} \psi _l (\varvec{A}^l)_{11}, \end{aligned}$$

where $(\psi _{k})_{k= 1}^{k= \text {deg}(\nu _{\varvec{A}})} \subset {\mathcal {D}}^{d}_{\Delta x}$ are the coefficients of a MPAMFR $\nu _{\varvec{A}}$ of $\varvec{A}$ being $\nu _{\varvec{A}} = \sum _{k= 0}^{k= \text {deg}(\nu _{\varvec{A}})} \psi _{k} X^{k}$, is the MPAFR ${\tilde{\nu }}_{\varvec{A}}$ of $\varvec{A}$.

Proof

We are only left to check Eq. (A4) for $j= 1$. $\square$

So in order to reduce the lattice Boltzmann scheme to a Finite Difference scheme using the new strategy, considering a MPAMFR or the MPAFR is exactly the same thing. Moreover, the MPAFR (but not the more general MPAMFR) can be linked to the minimal/characteristic polynomial.¹⁴

Lemma 16

Let $\mu _{\varvec{A}} \in {\mathcal {D}}^{d}_{\Delta x}[X]$ be the minimal polynomial of $\varvec{A}$, then ${\tilde{\nu }}_{\varvec{A}}$ exists and divides the minimal polynomial $\mu _{\varvec{A}}$. Moreover $\text {deg}({\tilde{\nu }}_{\varvec{A}}) = \text {deg}(\nu _{\varvec{A}}) \le \text {deg}(\mu _{\varvec{A}})$.

Proof

The proof goes like the standard one of Lemma 7. Consider $\mu _{\varvec{A}} = X^{\text {deg}(\mu _{\varvec{A}})} + \omega _{\text {deg}(\mu _{\varvec{A}})-1}X^{\text {deg}(\mu _{\varvec{A}})-1} + \dots + \omega _1 X+ \omega _0$. Consider the Euclidian division between $\mu _{\varvec{A}}$ and ${\tilde{\nu }}_{\varvec{A}}$: there exist $Q, R \in {\mathcal {D}}^{d}_{\Delta x}[X]$ such that

$$\begin{aligned} \mu _{\varvec{A}} = {\tilde{\nu }}_{\varvec{A}}Q + R, \end{aligned}$$

with either $0< \text {deg}(R) < \text {deg}({\tilde{\nu }}_{\varvec{A}})$ or $\text {deg}(R) =0$ (constant reminder polynomial). Let us indeed write

$$\begin{aligned} Q&= q_{\text {deg}(\mu _{\varvec{A}}) - \text {deg}({\tilde{\nu }}_{\varvec{A}})} X^{\text {deg}(\mu _{\varvec{A}}) - \text {deg}({\tilde{\nu }}_{\varvec{A}})} + \dots + q_1 X+ q_0, \\ R&= r_{\text {deg}(R)} X^{\text {deg}(R)} + \dots + r_1 X+ r_0, \\ \end{aligned}$$

Suppose that $R \not \equiv 0$, then we have for every $j\in \llbracket 1, q \rrbracket$

$$\begin{aligned} \overbrace{(\varvec{A}^{\text {deg}(\mu _{\varvec{A}})})_{1j} + \omega _{\text {deg}(\mu _{\varvec{A}})-1}(\varvec{A}^{\text {deg}(\mu _{\varvec{A}})-1})_{1j} + \dots + \omega _1 (\varvec{A})_{1j} + \omega _0 \delta _{1j}}^{=0} \\ = r_{\text {deg}(R)} (\varvec{A}^{\text {deg}(R})_{1j} + \dots + r_1 (\varvec{A})_{1j} + r_0\delta _{1j} + \\ \underbrace{\left( (\varvec{A}^{\text {deg}({\tilde{\nu }}_{\varvec{A}})})_{1j} + \psi _{\text {deg}({\tilde{\nu }}_{\varvec{A}})-1} (\varvec{A}^{\text {deg}({\tilde{\nu }}_{\varvec{A}})-1})_{1j} + \dots + \psi _1 (\varvec{A})_{1j} + \psi _0 \delta _{1j}\right) }_{= 0} \\ \times \left( q_{\text {deg}(\mu _{\varvec{A}}) - \text {deg}(\nu _{\varvec{A}})} (\varvec{A}^{\text {deg}(\mu _{\varvec{A}}) - \text {deg}(\nu _{\varvec{A}})})_{1j} + \dots + q_1 (\varvec{A})_{1j}+ q_0 \delta _{1j}\right) , \end{aligned}$$

thus

$$\begin{aligned} r_{\text {deg}(R)} (\varvec{A}^{\text {deg}(R})_{1j} + \dots + r_1 (\varvec{A})_{1j} + r_0 \delta _{1j} = 0, \qquad j\in \llbracket 1, q \rrbracket , \end{aligned}$$

with $0<\text {deg}(R) < \text {deg}({\tilde{\nu }}_{\varvec{A}})$, which contradicts the minimality of ${\tilde{\nu }}_{\varvec{A}}$. Thus necessarily $\text {deg}(R) = 0$ so the polynomial is constant, but to have the previous property, the constant must be zero, thus $R \equiv 0$. $\square$

Appendix B: Additional examples

In this Section, we gather more examples concerning the application of our theory to lattice Boltzmann schemes which can be found in the literature.

1.1 B. 1 $\text {D}_{1}\text {Q}_{2}$ with one conservation law

Consider the scheme by [9, 22] taking $d= 1$, $q= 2$ and $N= 1$ with $c_1 = 1$ and $c_2 = -1$ and

$$\begin{aligned} \varvec{M}= \left( \begin{matrix} 1 &{} 1 \\ \lambda &{} - \lambda \end{matrix} \right) , \qquad \text {and} \quad \varvec{S}= \text {diag}(0, s_2), \quad \text {with} \quad s_2\,\in ]0, 2]. \end{aligned}$$

(B5)

We recall that $m_2^{\text {eq}}: {\mathbb {R}}\rightarrow {\mathbb {R}}$ is a given function of $m_1$. The scheme can be used to simulate a non-linear scalar conservation law such as Eq. (7) using an acoustic scaling and a non-linear diffusion equation with a parabolic scaling, namely when $\lambda \propto 1/\Delta x$. However, the scheme is not rich enough to simulate more complex equations. As already pointed out in the introduction, the Finite Difference equivalent of this scheme has already been studied by [9] in the case where the equilibria are linear functions. It can be easily seen, even by hand since dealing with a $2\times 2$ matrix, that

The minimal polynomial coincides with the characteristic polynomial. This can be seen, as usual, by trying to consider $\alpha _0$ and $\alpha _1$ such that

$$\begin{aligned} \alpha _0 \varvec{I} + \alpha _1 \varvec{A}= \begin{pmatrix} \alpha _0 + \frac{({\mathsf {x}}+ \overline{{\mathsf {x}}})}{2}\alpha _1 &{} \frac{(1-s_2)({\mathsf {x}}- \overline{{\mathsf {x}}})}{2\lambda }\alpha _1 \\ \frac{\lambda ({\mathsf {x}}- \overline{{\mathsf {x}}})}{2}\alpha _1 &{} \alpha _0 + \frac{(1-s_2)({\mathsf {x}}+ \overline{{\mathsf {x}}})}{2}\alpha _1 \end{pmatrix} = \begin{pmatrix} 0 &{} 0 \\ 0 &{} 0 \end{pmatrix}. \end{aligned}$$

The only way of annihilating the first entry is to take $\alpha _0 = \alpha _1 = 0$, which is trivial. Thus the minimal polynomial is of degree two and then coincides with the characteristic polynomial. The equivalent Finite Difference scheme is

$$\begin{aligned} m_1^{n+ 1} = \frac{(2-s_2)}{2} ({\mathsf {x}}+ \overline{{\mathsf {x}}})m_1^{n} - (1-s_2)m_1^{n- 1} + \frac{s_2}{2\lambda } ({\mathsf {x}}- \overline{{\mathsf {x}}})m_2^{\text {eq}} \mid ^{n}. \end{aligned}$$

This scheme is a linear combination of the Lax-Friedrichs scheme (for $s_2 = 1$) with weight $2 - s_2$ and the leap-frog scheme (for $s_2 = 2$) with weight $-(1-s_2)$. It is a $\theta$-scheme between them with $\theta = 2 - s_2$ as long as $s_2 \in [1, 2]$.

1.2 B. 2 $\text {D}_{1}\text {Q}_{3}$ MRT for one conservation law

Consider the $\text {D}_{1}\text {Q}_{3}$ MRT scheme by [12], which reads with our notations $d= 1$, $q= 3$ and $N= 1$ with $c_1 = 0$, $c_2 = 1$ and $c_3 = -1$ and

$$\begin{aligned} \varvec{M}= \begin{pmatrix} 1 &{} 1 &{} 1 \\ 0 &{} \lambda &{} -\lambda \\ 0 &{} \lambda ^2 &{} \lambda ^2 \end{pmatrix}, \qquad \text {and} \quad \varvec{S}= \text {diag}(0, s_2, s_3), \quad \text {with} \quad s_2, s_3 \in ]0, 2], \end{aligned}$$

and $m_2^{\text {eq}}, m_3^{\text {eq}}: {\mathbb {R}}\rightarrow {\mathbb {R}}$ being given functions of $m_1$. The characteristic polynomial, corresponding to the minimal polynomial is

Then the corresponding Finite Difference scheme is

$$\begin{aligned} m_1^{n+ 1}&= m_1^{n} - \frac{(s_2 + s_3 - 2)}{2}({\mathsf {x}}+ \overline{{\mathsf {x}}}) m_1^{n} - (1-s_2) (1- s_3 ) m_1^{n- 1} \\&\quad + \frac{(s_2+ s_3-2)}{2}({\mathsf {x}}+ \overline{{\mathsf {x}}}) m_1^{n- 1} + (1-s_2)(1-s_3) m_1^{n- 2} \\&\quad + \frac{s_2}{2\lambda }({\mathsf {x}}- \overline{{\mathsf {x}}})m_2^{\text {eq}} \mid ^{n} - \frac{s_2(1-s_3)}{2\lambda } ({\mathsf {x}}- \overline{{\mathsf {x}}})m_2^{\text {eq}} \mid ^{n- 1} \\&\quad + \frac{s_3}{2\lambda ^2} ({\mathsf {x}}- 2 + \overline{{\mathsf {x}}})m_3^{\text {eq}} \mid ^{n} + \frac{s_3(1-s_2)}{2\lambda ^2} ({\mathsf {x}}- 2 + \overline{{\mathsf {x}}})m_3^{\text {eq}} \mid ^{n- 1}, \end{aligned}$$

This coincides with the one found by [12]. The SRT version comes by considering $s_2 = s_3$, thus having

$$\begin{aligned} m_1^{n+ 1}&= m_1^{n} + (1-s_2)({\mathsf {x}}+ \overline{{\mathsf {x}}}) m_1^{n} - (1-s_2)^2 m_1^{n- 1} -(1-s_2) ({\mathsf {x}}+ \overline{{\mathsf {x}}}) m_1^{n- 1} \\&\quad + (1-s_2)^2m_1^{n- 2} + \frac{s_2}{2\lambda }({\mathsf {x}}- \overline{{\mathsf {x}}})m_2^{\text {eq}} \mid ^{n} - \frac{s_2(1-s_2)}{2\lambda } ({\mathsf {x}}- \overline{{\mathsf {x}}})m_2^{\text {eq}} \mid ^{n- 1} \\&\quad + \frac{s_2}{2\lambda ^2} ({\mathsf {x}}- 2 + \overline{{\mathsf {x}}})m_3^{\text {eq}} \mid ^{n} + \frac{s_2(1-s_2)}{2\lambda ^2} ({\mathsf {x}}- 2 + \overline{{\mathsf {x}}})m_3^{\text {eq}} \mid ^{n- 1}, \end{aligned}$$

coinciding with the one found by [12] and by [8].

1.3 B. 3 $\text {D}_{2}\text {Q}_{4}$ for one conservation law

Consider $d= 2$, $q= 4$ and $N= 1$ with $\varvec{c}_1 = (1, 0)^{\intercal }$, $\varvec{c}_2 = (0, 1)^{\intercal }$, $\varvec{c}_3 = (-1, 0)^{\intercal }$ and $\varvec{c}_4 = (0, -1)^{\intercal }$ and

$$\begin{aligned} \varvec{M}= \left( \begin{matrix} 1 &{} 1 &{} 1 &{} 1\\ \lambda &{} 0 &{} -\lambda &{} 0 \\ 0 &{} \lambda &{} 0 &{} -\lambda \\ \lambda ^2 &{} -\lambda ^2 &{} \lambda ^2 &{} -\lambda ^2 \end{matrix} \right) , \qquad {\text {and}} \quad \varvec{S}= \text {diag}(0, s_2, s_2, 1), \quad \text {with} \quad s_2 \in ]0, 2]. \end{aligned}$$

(B6)

We have decided to set $s_3 = s_2$ not to favor one particular direction between x and y and $s_4 = 1$ to keep the scheme simpler. Moreover, we take $m_4^{\text {eq}} \equiv 0$ for simplicity and $m_2^{\text {eq}}, m_3^{\text {eq}} : {\mathbb {R}}\rightarrow {\mathbb {R}}$ are given functions of $m_1$. This can be used, for example, coupled with other schemes of the same nature (building what we call a “vectorial scheme” [39]) to easily simulate systems of non-linear conservation laws for $d= 2$, see [40], under acoustic scaling. After some computation, the characteristic polynomial of $\varvec{A}$ reads

where we have introduced the short-hands ${\mathsf {A}}_{\text {a}}{:}{=}({\mathsf {x}}+ \overline{{\mathsf {x}}} + {\mathsf {y}}+ \overline{{\mathsf {y}}})/4 \in {\mathcal {D}}^{d}_{\Delta x}$ and ${\mathsf {A}}_{\text {d}}{:}{=}({\mathsf {x}}{\mathsf {y}}+ {\mathsf {x}}\overline{{\mathsf {y}}} + \overline{{\mathsf {x}}} {\mathsf {y}}+ \overline{{\mathsf {x}}} \overline{{\mathsf {y}}})/4 \in {\mathcal {D}}^{d}_{\Delta x}$, yielding respectively the average between neighbors along the axis and along the diagonals. One can check as usual that it coincides with the minimal polynomial. The equivalent Finite Difference scheme is

$$\begin{aligned} m_1^{n+1}&= (3-2s_2){\mathsf {A}}_{\text {a}}m_1^{n} - (1-s_2) m_1^{n-1} - (1-s_2)(2-s_2){\mathsf {A}}_{\text {d}}m_1^{n-1} \\&+(1-s_2)^2 {\mathsf {A}}_{\text {a}}m_1^{n-2} +\frac{s_2}{2\lambda }({\mathsf {x}}- \overline{{\mathsf {x}}})m_2^{\text {eq}} \mid ^{n} + \frac{s_2}{2\lambda }({\mathsf {y}}- \overline{{\mathsf {y}}})m_3^{\text {eq}} \mid ^{n} \\&-\frac{s_2(1-s_2)}{4\lambda }\left( {\mathsf {y}}{({\mathsf {x}}- \overline{{\mathsf {x}}})} + \overline{{\mathsf {y}}} {({\mathsf {x}}- \overline{{\mathsf {x}}})}\right) m_2^{\text {eq}} \mid ^{n-1} \\&-\frac{s_2(1-s_2)}{4\lambda }\left( {\mathsf {x}}{({\mathsf {y}}- \overline{{\mathsf {y}}})} + \overline{{\mathsf {x}}} {({\mathsf {y}}- \overline{{\mathsf {y}}})}\right) m_3^{\text {eq}} \mid ^{n-1}, \end{aligned}$$

It is interesting to observe how the scheme computes the first-order derivatives in space of $m_2^{\text {eq}}$ and $m_3^{\text {eq}}$ at time $t^{n- 1}$.

1.4 B. 4 $\text {D}_{2}\text {Q}_{5}$ for one conservation law

In this example, we consider the scheme taken from [41]. Let $d= 2$, $q= 5$ and $N= 1$ with $\varvec{c}_1 = (0, 0)^{\intercal }$, $\varvec{c}_2 = (1, 0)^{\intercal }$, $\varvec{c}_3 = (0, 1)^{\intercal }$, $\varvec{c}_4 = (-1, 0)^{\intercal }$ and $\varvec{c}_5 = (0, -1)^{\intercal }$ and

$$\begin{aligned} \varvec{M}= \begin{pmatrix} 1 &{} 1 &{} 1 &{} 1 &{} 1 \\ 0 &{} \lambda &{} 0 &{} -\lambda &{} 0 \\ 0 &{} 0 &{} \lambda &{} 0 &{} -\lambda \\ -4\lambda ^2 &{} \lambda ^2 &{} \lambda ^2 &{} \lambda ^2 &{} \lambda ^2 \\ 0 &{} \lambda ^2 &{} -\lambda ^2 &{} \lambda ^2 &{} -\lambda ^2 \end{pmatrix}, \qquad \varvec{S}= \text {diag}(0, s_2, s_3, s_4, s_5), \end{aligned}$$

with $s_2, s_3, s_4, s_5 \in ]0, 2]$. One does not favor the x direction compared to the y direction, see [41], thus we set $s_2 = s_3$. Moreover, the assumptions on the moments at equilibrium are $m_2^{\text {eq}} \equiv m_3^{\text {eq}} \equiv m_5^{\text {eq}} \equiv 0$. The only one which does not vanish is $m_4^{\text {eq}} : {\mathbb {R}}\rightarrow {\mathbb {R}}$. The scheme can be used to solve thermic problems. Then, the characteristic polynomial of $\varvec{A}$ is given by

where the coefficients are given by

$$\begin{aligned} \gamma _4&= \frac{s_4}{5} (4 + {\mathsf {A}}_{\text {a}}) + {s_5}{\mathsf {A}}_{\text {a}}+ 2{s_2} {\mathsf {A}}_{\text {a}}- 1 - 4{\mathsf {A}}_{\text {a}}.\\ \gamma _3&= \frac{s_4 s_5}{5}(4 {\mathsf {A}}_{\text {a}}+ {\mathsf {A}}_{\text {d}}) + \frac{s_2 s_4 }{5} ( 1 + 8 {\mathsf {A}}_{\text {a}}+ {\mathsf {A}}_{\text {d}}) - \frac{s_4}{5} (1 + 17 {\mathsf {A}}_{\text {a}}+ 2 {\mathsf {A}}_{\text {d}}) \\&\quad + { s_2 s_5} ( 1 + {\mathsf {A}}_{\text {d}}) - s_5( 1 + {\mathsf {A}}_{\text {a}}+ 2{\mathsf {A}}_{\text {d}}) + s_2^2 {\mathsf {A}}_{\text {d}}- 2s_2 ( 1 + {\mathsf {A}}_{\text {a}}+ 2 {\mathsf {A}}_{\text {d}}) \\&\quad + 2 + 4{\mathsf {A}}_{\text {a}}+ 4 {\mathsf {A}}_{\text {d}}.\\ \gamma _2&= \frac{2s_2 s_4 s_5 }{5} ( 2 + {\mathsf {A}}_{\text {a}}+ 2 {\mathsf {A}}_{\text {d}}) - \frac{s_4 s_5}{5} (4 + 2 {\mathsf {A}}_{\text {a}}+ 9 {\mathsf {A}}_{\text {d}}) + \frac{ s_2^2 s_4}{5} ({\mathsf {A}}_{\text {a}}+ 4{\mathsf {A}}_{\text {d}}) \\&\quad -\frac{s_2 s_4 }{5} (9 + 4{\mathsf {A}}_{\text {a}}+ 17{\mathsf {A}}_{\text {d}}) +\frac{3 s_4}{5} ( 3 + {\mathsf {A}}_{\text {a}}+ 6 {\mathsf {A}}_{\text {d}}) +s_2^2 s_5 {\mathsf {A}}_{\text {a}}\\&\quad - s_2 s_5 ( 1 + 4{\mathsf {A}}_{\text {a}}+ {\mathsf {A}}_{\text {d}}) + {s_5}( 1 + 3{\mathsf {A}}_{\text {a}}+ 2{\mathsf {A}}_{\text {d}}) - s_2^2 ( 2 {\mathsf {A}}_{\text {a}}+ {\mathsf {A}}_{\text {d}}) \\&\quad + 2{s_2} ( 1 + 3 {\mathsf {A}}_{\text {a}}+ 2{\mathsf {A}}_{\text {d}}) - 2 - 4{\mathsf {A}}_{\text {a}}- 4{\mathsf {A}}_{\text {d}}. \\ \gamma _1&= -\frac{(1-s_2)}{20} \bigl ( 4 s_2 s_4 s_5 (1 + 4{\mathsf {A}}_{\text {a}}) -4 s_4 s_5 ( 1 + 14 {\mathsf {A}}_{\text {a}}) - 4s_2 s_4 (1 + 9{\mathsf {A}}_{\text {a}}) \\&\quad + 4s_4 ( 1 + 19 {\mathsf {A}}_{\text {a}}) - 20 s_2 s_5 (1 + {\mathsf {A}}_{\text {a}}) + 20 s_5 (1 + 3{\mathsf {A}}_{\text {a}}) + 20 s_2 ( 1 + 2{\mathsf {A}}_{\text {a}}) \\&\quad - 20 ( 1 + 4{\mathsf {A}}_{\text {a}}) \bigr ), \end{aligned}$$

and finally

$$\begin{aligned} \gamma _0 = - (1-s_2 )^2 (1-s_4) (1-s_5). \end{aligned}$$

Therefore, the corresponding Finite Difference scheme on the conserved variable $m_1$ reads

$$\begin{aligned} m_1^{n+ 1} = -\sum _{k= 0}^{4} \gamma _{4 - k} m_1^{n- k} + \sum _{k= 0}^3 \upsilon _{4 - k} m_4^{\text {eq}} \mid ^{n- k}, \end{aligned}$$

where the operators applied on the equilibrium $m_4^{\text {eq}}$ are given by

$$\begin{aligned} \upsilon _4&= \frac{s_4}{5 \lambda ^2} (- 1 + {\mathsf {A}}_{\text {a}}).\\ \upsilon _3&= \frac{s_4 s_5}{5 \lambda ^2} ( -{\mathsf {A}}_{\text {a}}+ {\mathsf {A}}_{\text {d}}) + \frac{s_2 s_4 }{5\lambda ^2} (1 - 2{\mathsf {A}}_{\text {a}}+ {\mathsf {A}}_{\text {d}}) - \frac{s_4}{ 5 \lambda ^2} (1 - 3{\mathsf {A}}_{\text {a}}+ 2{\mathsf {A}}_{\text {d}}).\\ \upsilon _2&= -\frac{s_2 s_4 s_5 }{5 \lambda ^2} (1 - 2{\mathsf {A}}_{\text {a}}+ {\mathsf {A}}_{\text {d}}) + \frac{s_4 s_5}{5\lambda ^2} (1 - 2{\mathsf {A}}_{\text {a}}+ {\mathsf {A}}_{\text {d}}) - \frac{ s_2^2 s_4}{5 \lambda ^2} ( -{\mathsf {A}}_{\text {a}}+ {\mathsf {A}}_{\text {d}}) \\&\quad +\frac{ s_2 s_4 }{5 \lambda ^2}(1 - 4{\mathsf {A}}_{\text {a}}+ 3{\mathsf {A}}_{\text {d}}) -\frac{s_4}{5\lambda ^2} (1 - 3{\mathsf {A}}_{\text {a}}+ 2{\mathsf {A}}_{\text {d}}), \end{aligned}$$

and finally

$$\begin{aligned} \upsilon _1&= -\frac{ s_2^2 s_4 s_5}{5\lambda ^2} (-1 + {\mathsf {A}}_{\text {a}}) + \frac{2s_2 s_4 s_5 }{5 \lambda ^2} ( -1 + {\mathsf {A}}_{\text {a}})- \frac{s_4 s_5}{5 \lambda ^2} (-1 + {\mathsf {A}}_{\text {a}}) \\&\quad + \frac{s_2^2 s_4 }{5\lambda ^2} (-1 + {\mathsf {A}}_{\text {a}}) - \frac{2s_2 s_4 }{5\lambda ^2} ( -1 + {\mathsf {A}}_{\text {a}}) + \frac{s_4}{5\lambda ^2} ( -1 + {\mathsf {A}}_{\text {a}}). \end{aligned}$$

This example shows that for large lattice Boltzmann schemes, the corresponding Finite Difference scheme can be rather complicated and not simple to analyze, especially for large $Q$.