Lattice Boltzmann schemes are a class of computational methods used to simulate systems of conservation laws under the form of Partial Differential Equations (PDEs). Their basic way of working is the following: instead of takingPDEs and directly discretize them, a lattice Boltzmann scheme enlarges the size of the problem from to and treats it in a kinetic-like fashion. This means that the new variables undergo, at each time step, a local collision phase where different particle distribution functions interact, followed by a lattice-constrained stream phase where no interaction is possible. The advantage of such idiosyncratic approach compared to more traditional numerical methods (e.g. Finite Difference, Finite Volume, Finite Elements, etc.) is that the local nature of the collision phase allows for massive parallelization of the method and the lattice-constrained stream can be computationally implemented as a pointer shift. Although this way of proceeding is highly beneficial from a computational perspective, it yields a deficient structure to construct a clear and rigorous notion of consistency with respect to the target equations, as well as a rigorous theory of stability. Indeed, only formal procedures, either based on the Chapman-Enskog expansion  or on the equivalent equations by Dubois [16, 18]
are currently available to study the consistency of lattice Boltzmann schemes. As far as stability is concerned, most of the studies rely on the linear stability analysis of the eigenvalues of the system, see[4, 37].
In order to bridge the gap between the lattice Boltzmann methods and the traditional approaches known to numerical analysts, the aim of the present contribution is to show that any lattice Boltzmann scheme can be rewritten as a corresponding multi-step Finite Difference scheme on the conserved variables, regardless of the linearity of the equilibria. This is made possible by developing an appropriate formalism based on commutative algebra and therefore yields a proper notion of consistency with respect to the target equations, which is that of Finite Difference schemes (see any standard textbook such as ). Furthermore, we confirm that the customary von Neumann analysis used for lattice Boltzmann schemes is equivalent to performing the same analysis on the corresponding Finite Difference scheme and is consequently particularly relevant. The price to pay for passing from an explicit scheme with variables and utilizing information only at the previous time-step to a method with variables is to increase the number of previous time-steps the new solution depends on, yielding a multi-step Finite Difference scheme.
In the past, few authors have noticed that for some particular lattice Boltzmann schemes, one has a corresponding (sometimes called “equivalent”) Finite Difference formulation on the conserved variables. Despite this, no general theory has been formulated. For instance: Suga  derives by direct computations a three-stages Finite Difference scheme from a uni-dimensional three-velocities scheme,111It is customary to call a scheme in a -dimensional space using discrete velocities. limiting the computations to a linear framework with one relaxation parameter (SRT). Dellacherie  derives a two-stages Finite Difference scheme for the lattice Boltzmann scheme. Again, this is limited to one spatial dimension and to a linear framework. A higher level of generality has been reached by the works of Ginzburg and collaborators, see  for a recap. They succeeded, using a link formalism, in writing a class of Lattice Boltzmann schemes as Finite Difference schemes 
. With their highly constrained link structure to be enforced, the resulting Finite Difference scheme with three stages is valid regardless of the spatial dimension and the choice of discrete velocities. The limitations are that the choice of moments is heavily constrained and only the case of one conserved moment is handled. Moreover, the evolution equation of the moving particles can depend on the distribution of the still particles onlyvia the conserved moment the equilibria depend upon and the schemes must be two-relaxation time (TRT) models with “magic parameter” equal to one-fourth for any link. The difficulty in establishing a general result comes from the coupling between spatial operators and time shifts. We must mention that during the drafting of the present contribution, an interesting work by Fuc̆ik and Straka  has been published covering the very same subject and essentially coming to the same conclusion as our paper. Their focus is different than ours since they adopt a purely algorithmic approach rather than a precise algebraic characterization of lattice Boltzmann schemes. We actually provide more insight into the bound on the number of time steps of the corresponding Finite Difference scheme and our formalism, based on polynomials, aims at providing a direct link with the classical tools for the stability analysis and allows to establish a link with the Taylor expansions from , as introduced in . In , the authors rely on a decomposition of the scheme using an hollow matrix222Matrix with zero entries on the diagonal. yielding an equivalent form of the scheme with the diagonal non-equilibrium part, after a finite number of steps of their algorithm. However, to the best of our understanding, the origin of such algorithm is not fully clear. In their work, the spatial shifts of data introduced by the stream phase are taken into account using a rather cumbersome system of indices, whereas we rely on an straightforward algebraic characterization of the stream phase.
Our paper is structured as follows: in Section 2
, we introduce – in guise of friendly introduction – the link of our problem with Ordinary Differential Equations (ODEs). The right formalism to make lattice Boltzmann schemes looking very close to a system of ODEs is provided inSection 3 and allows to prove the main results of the work showcased in Section 4. We devote Section 5 to discuss examples, possible simplifications of the problem and particular cases deserving particular attention. In Section 6, we prove the equivalence of the von Neumann analysis for lattice Boltzmann and Finite Difference schemes. In Section 7, we show how the well-known tools for Finite Difference schemes can be used to prove convergence theorems for lattice Boltzmann schemes. We corroborate our claim via numerical simulations. We eventually conclude in Section 8.
2 The example of Ordinary Differential Equations
Since our way of reducing any lattice Boltzmann scheme to a multi-step Finite Difference scheme has been originally inspired by an analogy with systems of ODEs, let us introduce this way of reasoning with the following example. Consider the system of ODEs of size with matrix given by
Transforming a single equation of higher order into a system of first order equations like Equation 1 by considering the companion matrix is a current practice, which unsurprisingly makes the problem more handy from the computational standpoint. Though being the analogous of what we aim at doing of lattice Boltzmann schemes, the other way around, passing from a system of first order to a single equation of higher order, seems to be seldom considered. We proceed like in . By iterating, we have that for .333We shall consistently use the notation for and . Let be real coefficients, then write . Taking as the coefficients of the characteristic polynomial444In the whole work, the indeterminate of any polynomial shall be denoted by . of , by virtue of the Cayley-Hamilton theorem, we deduce the corresponding equation on the first variable (playing the role of the conserved moment), given by
This provides a systematic way of performing the transformation without having to rely on hand computations and substitutions. To give an example, consider
Hence, the corresponding ODE on the first variable is given by .
3 Algebraic form of lattice Boltzmann schemes
Now that the reader is familiar – through a simple example – with the main idea and the final aim of the present contribution, we introduce the general framework of lattice Boltzmann schemes and the right formalism to treat them almost as systems of ODEs.
3.1 Spatial and temporal discretization
We set the problem in any spatial dimension considering the whole space , because we are not interested in studying boundary conditions. The space is discretized by a -dimensional lattice of constant step in all direction. The time is uniformly discretized with step . The discrete instants of time shall be indexed by the integer indices so that the corresponding time is . We finally introduce the so-called lattice velocity defined by . Observe that the developing theory is totally discrete and thus fully independent from the scaling between and .
3.2 Discrete velocities and shift operators
The first choice to be made when devising a lattice Boltzmann scheme concerns the discrete velocities with , which are multiples of the lattice velocity, namely for any with . Therefore, particles are stuck to move – at each time step – on the lattice . We denote the distribution density of the particles moving with velocity by for every . The shift operators associated with the discrete velocities are an important element of the following analysis.
Definition 1 (Shift operator).
Let , then the associated shift operator on the lattice , denoted , is defined in the following way. Take be any function defined on the lattice,555The function could take values in any ring, see . then the action of is
We also introduce .
The shift yields information sought in the upwind direction with respect to the considered velocity. Let us introduce the natural binary operation between shifts.
Definition 2 (Product).
Let the “product” be the binary operation defined as , for any .
Henceforth, the product is understood whenever no ambiguity is possible. This operation provides an algebraic structure to the shifts, directly inherited from that of .
forms an Abelian group.
Moreover, there is only “one movement” for each Cartesian direction which “generates” the shifts. More precisely
where is the customary notation for the generating set of a group. We can add one more binary operation, which is non-internal to . This yields the cornerstone of this work, namely the set of Finite Difference operators, finite combinations of weighted shifts operators via a sum. It is defined as follows, see Chapter 3 of .
Definition 3 (Finite Difference operators).
The set of Finite Difference operators on the lattice is defined as
the group ring (or group algebra) of over . The sum the product666Which interestingly corresponds to the discrete convolution product. of two elements are defined by
Furthermore, the product of with elements of is given by
With the two binary operations, behaves closely to , or as stated by the following result, see .
Proposition 2 (Ring of Finite Difference operators).
is a commutative ring.777It also an (Hopf) algebra over and can also be viewed as a free module where the scalars belong to and the basis are the elements of the group .
Observe that is not a field: not every element of has multiplicative inverse, take for example the centered approximation of the derivative along : and see for instance the concept of indefinite sum in the calculus of Finite Differences [33, 32]. The elements having inverse are called “units” and divide all the other elements. It can be easily seen that the units are the product of a non-zero real number and a shift in . Indeed for any and . The inverse of a unit shall also be denoted by a bar.
One can see as the ring of Laurent polynomials of variables over the field , where the indeterminates are , and . For example, for , the identification holds. This automatically implies that is more than a commutative ring, namely a unique factorization domain.
The reals can be identified with the subring .
3.3 Lattice Boltzmann algorithm: collide and stream
Any lattice Boltzmann scheme consists in an algorithm made up of two phases: a local collision phase performed on each site of the lattice and a stream phase, where particles are exchanged between different sites of the lattice. Let us introduce each of them.
3.3.1 Collision phase
We adopt the point of view of the multiple-relaxation-times (MRT) schemes, where it is customary to consider the collision written as a diagonal relaxation in the moments basis, see . For this reason, we introduce a change of basis called moment matrix . The entries of can depend on and/or on but cannot be a function of the space and time variables. Gathering the distributions into , the moments are recovered by . We also introduce
which is the identity matrix of size;
the matrix is the relaxation matrix which is a singular with , where is the number of conserved moments:
where the first entries are zero888This is not always the case in literature but shall be used consistently in this paper. We put them at the beginning for the sake of presentation. and correspond to the conserved moments, the following are such that for , see .
We employ the notation for , where are possibly non-linear functions of the conserved moments. Since these equilibria are then multiplied by , the first components do not need to be defined.
The collision phase reads, denoting by any post-collision state
In the collision phase Equation 4, the entries of can depend on or , but not on space and time. The equilibria are allowed to follow the same dependencies plus those on space and time and can also depend on some “external variable” like in the case of vectorial schemes .
3.3.2 Stream phase
The stream phase is diagonal in the space of the distributions. It can be written as
where for the first time, the matrices have entries in a commutative ring, see  and , instead than in the field . The set of square matrices of size with entries belonging to forms a ring under the usual operations between matrices. Even if is commutative from Proposition 2, is not commutative for , as for real matrices and matrices of first-order differential operators .
3.3.3 Monolithic scheme
The stream phase Equation 5 can be rewritten in a non-diagonal form in the space of moments as done by [18, 20] by introducing the matrix and merged with the collision phase Equation 4 to obtain the scheme
where and . In the sequel, we shall not indicate the spatial variable for the sake of readability.
We observe that the operators are the eigenvalues of the matrix . However, they are not the eigenvalues of the matrix . Indeed, it is general false that the eigenvalues of belong to the space . It is interesting to interpret the lattice Boltzmann scheme under the form Equation 6 as discrete-time linear control system with matrices on a commutative ring . The moments are the state of the system evolving via the matrix , whereas the equilibria are the control via being a feedback observing only a part of the state, namely the conserved moments.
We introduce our example of choice, which shall be used through the whole paper.
4 Main result of the paper
With a new way of writing any lattice Boltzmann scheme using Definition 3 and thanks to Proposition 2, which provides the ideal setting to generalize the Cayley-Hamilton theorem, we can proceed like in Section 2 to prove the main result of the paper: any lattice Boltzmann can be viewed as a multi-step Finite Difference scheme on the conserved variables.
4.1 Characteristic polynomial and Cayley-Hamilton theorem
Polynomials with coefficients in and matrices with entries in play a central role in what we are going to develop.
Definition 4 (Characteristic polynomial).
Let be a commutative ring and for some . The characteristic polynomial of , denoted , is given by , where is the determinant and is the identity matrix.
The naive computation of the characteristic polynomial using its definition via the determinant could be computationally expensive, especially when dealing with symbolic computations like in our case. For this reason, we employed the Faddeev-Leverrier algorithm  which is of polynomial complexity, generally lower than that of the pivot method.
The process is detailed in Algorithm 1 and only uses matrix-matrix multiplications and the computation of the trace, denoted by .
A central result used in this work is the Cayley-Hamilton theorem for matrices over a commutative ring, see  for the proof, generalizing the same result holding for matrices on a field utilized in Section 2.
Theorem 3 (Cayley-Hamilton).
Let be a commutative ring and for some . Then is a monic polynomial in the ring in the indeterminate , under the form with . Then999Sometimes, we shall indulge to the notation . .
This result states that any square matrix with entries in a commutative ring verifies its characteristic equation.
4.2 Corresponding Finite Difference schemes
4.2.1 One conserved moment
We first analyze the case of one conserved moment, namely , to keep the presentation as simple as possible. We shall eventually deal with once the principles are established.
Proposition 4 (Corresponding Finite Difference scheme for ).
Let , then the lattice Boltzmann scheme Equation 6 corresponds to a multi-step explicit Finite Difference scheme on the conserved moment under the form
where are the coefficients of , the characteristic polynomial of .
This result means that the conserved moment satisfies an explicit multi-step Finite Difference scheme with at most steps, thus involving discrete time instants, see Figure 1. The maximal size of spatial influence at each past time step can be deduced by looking at Algorithm 1, derived from the Newton’s identities.
It is interesting to observe that also the non-conserved moments satisfy a Finite Difference numerical scheme, see the following proof. However, these schemes would depend on the conserved moment via the equilibria and are therefore not independent from the rest of the system.
Let . Then for any , applying Equation 6 recursively we have
We perform a temporal shift in order to fix the first term on the right hand side regardless of the value of . Introduce , therefore
This holds true, in particular, for any . We can then consider the coefficients of the characteristic polynomial of and write
The last sum can start from . Performing a change of indices in the last double sum yields the result.
One could think of allowing and/or to depend on the space and time variables. This would imply to consider weights made up of functions instead of the real numbers in Definition 3. However, would no longer be commutative, because the multiplication by a function does not commute with shifts (not shift-invariant according to ). For example, take and a function , then
for every and for any function . The right-hand sides are not equal in general, except if is constant.
4.2.2 Several conserved moments and vectorial schemes
Consider now to deal with multiple conservation laws, namely . We select a conserved moment and we consider the other conserved moments as “slave” variables as the equilibria have been until so far, for , because they imply variables that we eventually want to keep. In particular, we utilize different polynomials for different conserved moments to obtain the Finite Difference schemes. To formalize this concept, for any square matrix , consider for any , corresponding to the matrix where only the rows and columns of indices are conserved and the remaining ones are set to zero. We can also consider the matrix obtained by keeping only the rows and the columns indexed in . A useful corollary of Theorem 3 and of the Laplace formula for the determinant is the following.
Let and , then one has that . Moreover, the polynomial annihilates .
This means that the characteristic polynomial of is directly linked to that of the smaller matrix , which is thus faster to compute, and that the latter is an annihilator for the first matrix.
For any conserved moment indexed by we introduce the matrix and . Notice that we have the decomposition . Indeed, we “save” the conserved moments other than the by placing them into , which shall not participate in the computation of the characteristic polynomial. With this notations, we have generated a family of problems from Equation 6 under the form
It is useful to stress that the term in Equation 9 does not involve any conserved moment other than the . Conversely, does not involve any function except the conserved moments other than the . Then, the corresponding Finite Difference schemes come under the form stated by the following Proposition.
Proposition 6 (Corresponding Finite Difference scheme for ).
This Proposition states that for each conserved moment, the corresponding Finite Difference scheme has at most steps, thus involves discrete times. This result encompasses and generalizes Proposition 4. The proof is the same than that of Proposition 4 by taking advantage of Corollary 5. We show in another contribution  that the result of Proposition 6 is the right one to bridge between the consistency analysis of Finite Difference schemes and the Taylor expansions on the lattice Boltzmann schemes for proposed by .
Example 4 ( for two conservation laws).
Consider the scheme  with , and , and
thus having . This scheme can be used to simulate the system of conservation laws and under the acoustic scaling . Using Proposition 6 we have
One could remark that the linear part is different from one scheme to the other, since we have used different polynomials for each conserved moment.
4.3 Initialization schemes
In the corresponding Finite Difference schemes in Proposition 4 and Proposition 6, the only remaining freedom is to devise the initialization schemes for the multi-step schemes at regime, analogously to Equation 2 in Section 2. This is the counter-part of the freedom of choice on the initial data for the original lattice Boltzmann scheme, which are not necessarily taken at equilibrium, see . By applying Equation 6 to the initial data as many times as needed, one progressively obtains the initialization schemes, as function of the initial datum. It is worthwhile observing that the choice of initial datum does not play any role in the previous procedure and does not influence the stability analysis of Section 6. It only comes into play during the consistency analysis of the numerical method, which is not investigated in this paper, in particular, as far as time boundary layers are concerned, see [40, 34].
5 Examples, simplifications and particular cases
Now that the main results of the paper, namely Proposition 4 and Proposition 6, have been stated and proved, we can analyze and comment some particular cases which deserve a closer look. More examples are available in the Appendices.
Example 5 (ODEs).
To illustrate some basic peculiarities that easily transpose to lattice Boltzmann schemes, we introduce the following matrices extending the discussion of Section 2.
For , we have , corresponding to . However, contrarily to in Section 2, the characteristic polynomial does not correspond to the minimal polynomial . Thus in this case, we could use the latter to obtain Equation 2 having . This phenomenon is studied in Section 5.1. It indicates that we can achieve a more compact corresponding ODE by using the annihilating polynomial of smallest degree on every variable. This does not change the core of the strategy.
For , we obtain , corresponding to . However, by inspecting , one notices that the first two equations do not depend on the last variable . For this reason, we could have considered the matrix obtained from by removing the last row and column. In this case , corresponding to the equation . This kind of situation for lattice Boltzmann schemes is investigated in Section 5.2. It is interesting to observe once more that divides . This shows that an initial inspection of the matrix can yield a reduction of the size of the problem that can be achieved by a simple trimming operation, which eliminates some variable from the problem but treats the remaining ones as usual.
Finally, consider . In this case the characteristic polynomial and the minimal polynomial coincide corresponding to the equation . However, if we take the polynomial such that divides