# Construction and performance of kinetic schemes for linear systems of conservation laws

Emmanuel Audusse, Sébastien Boyaval, Virgile Dubos, Minh-Hoang Le

<figure><img src="image_1.png" /></figure>

## To cite this version:

Emmanuel Audusse, Sébastien Boyaval, Virgile Dubos, Minh-Hoang Le. Construction and performance of kinetic schemes for linear systems of conservation laws. 2025. hal-05405231

HAL Id: hal-05405231

https://enpc.hal.science/hal-05405231v1

Preprint submitted on 8 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# CONSTRUCTION AND PERFORMANCE OF KINETIC SCHEMES FOR LINEAR SYSTEMS OF CONSERVATION LAWS

EMMANUEL AUDUSSE, SÉBASTIEN BOYAVAL, VIRGILE DUBOS, AND MINH LE

RÉSUMÉ. We describe a methodology to build vectorial kinetic schemes, targeting the numerical solution of linear symmetric-hyperbolic systems of conservation laws – a minimal application case for those schemes. Precisely, we fully detail the construction of kinetic schemes that satisfy a discrete equivalent to a convex extension (an additional non-trivial conservation law) of the target system – the (linear) acoustic and elastodynamics systems, specifically –. Then, we evaluate numerically the convergence of various possible kinetic schemes toward smooth solutions, in comparison with standard finite-difference and finite-volume discretizations on Cartesian meshes. Our numerical results confirm the interest of ensuring a discrete equivalent to a convex extension, and show the influence of remaining parameter variations in terms of error magnitude, both for “first-order” and “second-order” kinetic schemes : the parameter choice with largest CFL number (equiv., smallest spurious diffusion in the equivalent equation analysis) has the smallest discretization error.

## 1. INTRODUCTION

Systems of Conservation Laws (SCL) with smooth fluxes $F_a(q) \in \mathbb{R}^L$, $1 \le a \le d$,

$$ (1) \qquad \partial_t q^l + \sum_{a=1}^{d} \partial_a F_a^l(q) = 0 \quad 1 \le l \le L $$

that satisfy an additional non-trivial conservation law

$$ (2) \qquad \partial_t \eta(q) + \sum_{a=1}^{d} \partial_a G_a(q) = 0 $$

where $\eta(q)$ is a scalar $C^2(\mathbb{R}^L)$ functional strictly convex in $q$ and $G_a$ are smooth “entropy” fluxes, are widely spread in the physical sciences. They have a symmetric-hyperbolic formulation so solutions $q \in C^0([0, T], [H^s(\mathbb{R}^d)]^L)$ are unequivocally defined on small times $T$ given initial conditions $q(t=0) = q^0 \in [H^s(\mathbb{R}^d)]^L$, $s > 1+\frac{d}{2}[12, 5]$. The particular case of linear symmetric-hyperbolic systems $F_a^l(q) = \sum_{m=1}^L A_a^{lm} q_m$, $A_a^{lm} = A_a^{ml}$, where $\eta(q) = \frac{1}{2}|q|^2$ and $G_a(q) = \frac{1}{2} \sum_{m,l=1}^L A_a^{lm} q_m q_l$, is a minimal benchmark to evaluate the accuracy of numerical approximations to $q(t \equiv x_0, x_1, \dots, x_d)$. Noticeably, solutions are then well-defined whatever $T > 0$ and whatever $s \in \mathbb{R}$. The same remarkable properties extend to Friedrichs symmetrizable linear systems i.e. linear systems (1) for which a symmetric, positive-definite

<footnote>Thanks to ANR for its grant 15-CE01-0013 to the SEDIFLO project. Thanks to Lucas Brélivet for its help as an ENPC student in numerical simulations.</footnote>

<footnote>AMS classification : 76A10; 35L45; 74D10.</footnote>
matrix $S_0$ exists so that $S_0A_a = A_a^T S_0$ ($S_0^\frac{1}{2}q$ is solution to a symmetric-hyperbolic system then, with $\eta(q) = q^T S_0 q/2$). Friedrichs symmetrizable linear systems are physically relevant : here, we concentrate on the linear flux case having in mind the evaluation of discretization methods for linear *elastodynamics*.

Many numerical methods have been proposed to compute solutions $q(t, x_1, \dots, x_d)$ to (1). Various Finite-Difference (FD) [19] and Finite-Volume (FV) methods [23, 18] have proved able to approximate well solutions. Yet, as demonstrated e.g. in [21], standard techniques leave room to reach accuracy at a reasonable computational price, especially for multidimensional solutions to nonlinear systems. See also [9] for a recent study of Discontinuous Galerkin (DG) methods toward solutions to linear (symmetric-hyperbolic) systems, or [22] for the long-time behaviour of FV schemes (the particular DG-0 case).

In this work, we study the construction and numerical performance of *vector(ial) kinetic schemes* as defined in Section 2, for linear symmetric-hyperbolic systems of conservation laws, with a view to founding the development of the kinetic approach to further cases. After recalling some elements of analysis about vectorial kinetic schemes in Section 2, we develop in Section 3 a *construction method* for the paradigmatic (linear) acoustics and elastodynamics systems. Finally, the construction method is first based on ensuring a discrete equivalent to a convex extension of the SCL at, as well as preserving some (rotational) symmetries of the SCL. Next, it leaves one with various choices for parameters. In Section 4, we numerically evaluate the performance of specific schemes (with specific parameter values) in comparison with standard FD and FV schemes.

Our numerical results first confirm the interest of ensuring a discrete equivalent to a convex extension, both for “first-order” and “second-order” kinetic schemes : ensuring a discrete equivalent to a convex extension through adequate parameter choices appears as a sufficient and necessary condition of stability in our numerical tests. Moreover, both first and second order schemes do have first and second order convergence rates. Last, the remaining parameter choices (within the numerical stability region) show a strong influence on the magnitude of the discretization error. The parameter choice with largest CFL number (equiv., the one with smallest spurious diffusion in the equivalent equation analysis) has the smallest discretization error, smaller than standard second order FD and FV schemes.

Our numerical results are an incentive to further studies about the convergence and the error of kinetic schemes for (symmetric-hyperbolic) SCL

## 2. VECTORIAL KINETIC SCHEMES

Kinetic schemes use $Q$ kinetic variables $f_\zeta$, $\zeta \in \{0, 1, \dots, Q-1\}$, which are solutions to (discretized) Boltzmann-BGK equations with speeds $c_\zeta \in \mathbb{R}^d$ and a relaxation time $\varepsilon \ll T$. Here we consider *vectorial* (or vector) kinetic schemes, with vectorial unknown $f_\zeta \in \mathbb{R}^L$ solutions to Boltzmann-BGK equations

$$ (3) \qquad \partial_t f_{\zeta,l} + \sum_{a=1}^{d} c_{\zeta}^{a} \partial_a f_{\zeta,l} = \frac{\bar{f}_{\zeta,l}(Mf) - f_{\zeta,l}}{\varepsilon} $$

for $l \in \{1, \dots, L\}$, using “Maxwellians” $\bar{f}_{\zeta,l}$ functions of moments of $f_{\zeta,l}$ like

$$ (4) \qquad M_0 f_l := \sum_{\zeta} f_{\zeta,l}. $$
Here, we consider $Q \ge d$ such that the speeds $c_\zeta \in \mathbb{R}^d$ define lattices of $\mathbb{R}^d$ parametrized by a uniform time-step $k = T/N > 0$ ($N \in \mathbb{N}^*$): for all $\zeta = 0, \dots, Q-1$ it holds $a_{ij} + kc_\zeta = a_{i'j'}$ for some $i', j' \in \mathbb{Z}$. Such (vectorial) kinetic schemes are sometimes referred to as (vector) “Lattice Boltzmann Methods” (LBMs).

In this work, we consider square lattices equiv. Cartesian grids with $c_0 = 0$

$$ (5) \qquad \Delta = \inf_{ij \neq i'j'} |a_{ij} - a_{i'j'}| = k \min_{\zeta>0} |c_\zeta|, $$

$a_{ij}$ has coordinates $(i\Delta, j\Delta)$. Denoting $f_{\zeta,l}^n(a_{ij})$ the kinetic variables for computing $q_{ij}^n \approx q(kn, a_{ij})$, $n = 1, \dots, N$, we are interested in kinetic schemes as e.g. [10, 2], which read on 2D grids of points $a_{ij}$ with coordinates $(\Delta i, \Delta j)$

$$ (6) \qquad f_{\zeta,l}^{n+1}(a_{ij} + kc_{\zeta}) = (1-\omega)f_{\zeta,l}^{n}(a_{ij}) + \omega \bar{f}_{\zeta,l}(M_0 f^{n}(a_{ij})) $$

where for all $\zeta$ and $ij$, $a_{i'j'} = a_{ij} + kc_\zeta$ is uniquely defined on the lattice, while $\omega \in (1, 2)$ ranges for “first-order” schemes when $\omega = 1$ to “second-order” schemes when $\omega = 2$. The latter schemes coincide with a two-step time-integration of (3):

$$ (7) \qquad f_{\zeta,l}^{n+\frac{1}{2}}(a_{ij}) = (1-\omega)f_{\zeta,l}^{n}(a_{ij}) + \omega \bar{f}_{\zeta,l}(M_0 f^{n}(a_{ij})) $$

$$ (8) \qquad f_{\zeta,l}^{n+1}(a_{ij} + kc_{\zeta}) = f_{\zeta,l}^{n+\frac{1}{2}}(a_{ij}) $$

resulting from a Projection–Transport operator–splitting on $t \in [nk, (n+1)k)$.

The approximation of (1) by (6) is not yet fully understood, even in the case of (smooth) solutions to linear SCL, to our knowledge. Various stability conditions have been proposed in the literature [7, 20], as well as elements of error analysis (i.e. heuristics for a future rigorous proof) [14, 26, 4, 3], which presently guide the construction of a kinetic scheme. But we are not aware of sufficient conditions, about finitely-many kinetic speeds or “celerities” $c_\zeta$ and associated vector Maxwellians $\bar{f}_\zeta$, that guarantee a priori the convergence of (6) to (1). That is why here, focusing on the most simple SCL, the linear acoustics/elastodynamics systems in particular, we explore some sensible choices of speeds and Maxwellians numerically (see Section 4). To that aim, let us first recall what are sensible choices presently.

First, the (moments $M_0 f^n$ of the) solutions to (6) are known to solve equivalent equations consistent with (1) provided the Maxwellians $\bar{f}$ satisfy

$$ (9) \qquad M_0 \bar{f}(q) = q $$

$$ (10) \qquad M_a \bar{f}(q) = F_a(q) $$

for all $q \in \mathbb{R}^L$, on denoting $M_0$ the operator defined in (4), and $M_a : \mathbb{R}^{Q \times L} \to \mathbb{R}^L$

$$ (11) \qquad M_a f_l := \sum_{\zeta=0}^{Q-1} c_\zeta^a f_{\zeta,l} $$

see e.g. [3] for a recent statement of that result. Thus, (9) and (10) are standardly required in the construction of kinetic schemes for (1). Choosing speeds and Maxwellians that satisfy (9) and (10) is sensible, it provides one with $(d+1)L$ equations for $QL$ unknowns $\bar{f}$. But (9) and (10) are not enough in general to define uniquely the Maxwellians $\bar{f}$ (e.g. when using standard lattices with $Q > d+1$ celerities and linear Maxwellians). Moreover, letting alone the fact that (9) and (10) are not enough to guarantee convergence a priori, (9) and (10) do also not obviously provide one with a (simple, pointwise) error analysis in the case of smooth solutions – unlike standard FD and FV schemes as Yee’s scheme or the upwind scheme (see
Sections 4, A.1 and A.2). Now, an error analysis is one important step in proving convergence, independent of stability, which still needs studying [4]. So, with a view to complementing (9) and (10) here, let us now provide some elements for a detailed error analysis in the case of smooth solutions, inspired by the heuristics [7, Proposition 2.1] leading to an equivalent equations when $\omega = 1$.

**Proposition 1.** Given $\{c_\zeta, \zeta = 0 \dots Q - 1\}$ consider Maxwellians $\bar{f}_\zeta : \mathbb{R}^L \to \mathbb{R}^L$ linear, that satisfy (9)-(10) and

$$ (*) \quad \sum_{\zeta=0}^{Q} \max_{\|q\|_{\infty}=1} |\bar{f}_{\zeta,l}(q)| \le 1 \quad \forall l = 1 \dots L. $$

Then, for any smooth solution

$$ (12) \quad q \in C_t^0([0, T], [H^2(\mathbb{R}^d)]^L) \cap C_t^1([0, T], [H^1(\mathbb{R}^d)]^L) \cap C_t^2([0, T], [H^0(\mathbb{R}^d)]^L) $$

to linear conservation laws (1) with initial condition $q(t=0) = q^0 \in [H^2(\mathbb{R}^d)]^L$, there exists $C_T > 0$ such that it holds $\forall k = T/N, N \in \mathbb{N}^*$:

$$ (13) \quad \|q(T) - M_0 f^N\|_0 := \sqrt{\int_{\mathbb{R}^d} d\mathbf{a} \sum_{l=1}^L |q_l(T, \mathbf{a}) - M_0 f_l^N(\mathbf{a})|^2} \le kC_T $$

when, for each lattice $\{\mathbf{a}_{ij}\}$, $f_{\zeta,l}^n(\mathbf{a}_{ij})$ is computed by (6) with $\omega = 1$, and

$$ (14) \quad f^0 = \bar{f}(q^0). $$

Denoting $F_a^l(q) = \sum_{m=1}^L A_a^{lm} q_m$ the flux, the constant $C_T$ can be chosen as follows :

$$ C_T = 2T \left( \max_{\zeta} |c_{\zeta}| + \max_a \|A_a\|_{\infty} \right)^2 \max_{a,b} \sup_{s \in (0,T)} \|\partial_{ab}^2 q\|_0(s). $$

**Démonstration.** First, for $\omega = 1$ the “transport-projection” scheme (6) rewrites

$$ (15) \quad f_{\zeta,l}^{n+1}(\mathbf{a}) = \bar{f}_{\zeta,l}(M_0\{f_{\zeta}^n(\mathbf{a} - k c_{\zeta})\}) $$

whatever $k > 0$. Then, since $M_0 f^0 = q^0 \in H^2(\mathbb{R}^d)$, one can establish by recurrence for $n = 1 \dots N$ that $M_0 f^n = \sum_\zeta \bar{f}_\zeta(M_0\{f_\zeta^{n-1}(\cdot - k c_\zeta)\}) \in H^2(\mathbb{R}^d)$, and $\max_{a,b} \|\partial_{ab}^2 M_0 f^n\|_0 \le \max_{a,b} \|\partial_{ab}^2 q^0\|_0$ by assumptions on $\bar{f}_\zeta$.

Moreover, $f_{\zeta,l}^{n+1}(\mathbf{a}')$ in (8) is equivalently $f_{\zeta,l}(t^{n+1}, \mathbf{a}')$ the solution to

$$ (16) \quad (\partial_t + c_\zeta \cdot \nabla) f_{\zeta,l} = 0 \quad \forall t \in (t^n, t^{n+1}), \mathbf{a} \in \mathbb{R}^d \quad ; f_{\zeta,l}(t^n, \mathbf{a}) = f_{\zeta,l}^{n+\frac{1}{2}}(\mathbf{a}) $$

which possesses the regularity (12) because $f_{\zeta,l}^{n+\frac{1}{2}} = \bar{f}_{\zeta,l}(M_0 f^n) \in H^2(\mathbb{R}^d)$. Then, using (16) to Taylor expand $f_{\zeta,l}^{n+1}(\mathbf{a}) = f_{\zeta,l}(t^n + k, \mathbf{a})$ implies whatever $\mathbf{a} \in \mathbb{R}^d$

$$ (17) \quad f_{\zeta,l}^{n+1}(\mathbf{a}) = f_{\zeta,l}^{n+\frac{1}{2}}(\mathbf{a}) - k(c_{\zeta} \cdot \nabla) f_{\zeta,l}^{n+\frac{1}{2}}(\mathbf{a}) + \int_{t^n}^{t^{n+1}} (t^{n+1} - s) \partial_{tt}^2 f_{\zeta,l}(s, \mathbf{a}) ds $$
thus, using (9)–(10) and recalling the projection step $f_{\zeta,l}^{n+\frac{1}{2}}(\mathbf{a}) = \bar{f}_{\zeta,l}(M_0 f^n(\mathbf{a}))$,

$$ (18) \quad M_0 f_l^{n+1}(\mathbf{a}) = M_0 f_l^n(\mathbf{a}) - k \sum_{a=1}^d \partial_a F_a^l (M_0 f^n(\mathbf{a})) \\ \qquad + \int_{t^n}^{t^{n+1}} (t^{n+1} - s) \partial_{tt}^2 M_0 f_l(s, \mathbf{a}) ds. $$

Eq. (18) can be compared with a Taylor expansion of $q(nk + k, \mathbf{a})$ using (1) i.e.

$$ (19) \quad q_l((n+1)k, \mathbf{a}) = q_l(nk, \mathbf{a}) - k \sum_{a=1}^{d} \partial_a F_a^l (q(nk, \mathbf{a})) \\ \qquad + \int_{t^n}^{t^{n+1}} (t^{n+1} - s) \partial_{tt}^2 q_l(s, \mathbf{a}) ds $$

where $\partial_{tt}^2 q_l = \sum_{a,b=1}^d \sum_{m,m'=1}^L A_a^{lm} A_b^{mm'} \partial_{ab}^2 q_{m'}$ here in the linear flux case.
Subtracting (19) from (18), moreover using (16) we obtain

$$ (20) \quad q_l((n+1)k, \mathbf{a}) - M_0 f_l^{n+1}(\mathbf{a}) \\ = q_l(nk, \mathbf{a}) - M_0 f_l^n(\mathbf{a}) - k \sum_{a=1}^d \sum_{m=1}^L \partial_m F_a^l (\partial_a q_m(nk, \mathbf{a}) - \partial_a M_0 f_m^n(\mathbf{a})) + R_l^n(\mathbf{a}) $$

with a remainder $R_l^n(\mathbf{a})$ that satisfies by Minkowski's inequality

$$ (21) \quad \|R^n\|_0 \le k^2 \max_{a,b} \left( \max_{\zeta} |c_{\zeta}|^2 \|\partial_{ab}^2 M_0 f^n\|_0 + \max_a \|A_a\|_{\infty}^2 \sup_{s \in (t^n, t^{n+1})} \|\partial_{ab}^2 q\|_0(s) \right). $$

Using also Minkowski's inequality with (20) yields

$$ (22) \quad \|q((n+1)k, \cdot) - M_0 f^{n+1}(\cdot)\|_0 \le \|q(nk, \cdot) - M_0 f^n(\cdot)\|_0 \\ \qquad + k \max_{a,l,m} |\partial_m F_a^l| \sum_{a=1}^d \|\partial_a q(nk, \cdot) - \partial_a M_0 f^n(\cdot)\|_0 + \|R^n\|_0. $$

On the other hand, the following Taylor expansions

$$ (23) \quad M_0 f_l^{n+1}(\mathbf{a}) = M_0 f_l^n(\mathbf{a}) + \int_{t^n}^{t^{n+1}} \partial_t M_0 f_l(s, \mathbf{a}) ds $$

$$ (24) \quad q_l((n+1)k, \mathbf{a}) = q_l(nk, \mathbf{a}) + \int_{t^n}^{t^{n+1}} \partial_t q_l(s, \mathbf{a}) ds $$

yield, using differentiation by $\partial_a$, subtraction and Minkowski's inequality

$$ (25) \quad \|\partial_a q(nk, \cdot) - \partial_a M_0 f^n(\cdot)\|_0 \le \|\partial_a q((n-1)k, \cdot) - \partial_a M_0 f^{n-1}(\cdot)\|_0 \\ \qquad + k \max_b \left( \max_{\zeta} |c_{\zeta}| \|\partial_{ab}^2 M_0 f^{n-1}\|_0 + \max_{\alpha} \|A_{\alpha}\|_{\infty} \sup_{s \in (t^{n-1}, t^n)} \|\partial_{ab}^2 q\|_0(s) \right). $$

Using (25) in (22) and iterating for $n = N - 1 \dots 0$ yields

$$ (26) \quad \|q(T, \cdot) - M_0 f^N(\cdot)\|_0 \le \|q(0, \cdot) - M_0 f^0(\cdot)\|_0 \\ \qquad + (kN \equiv T) \max_a \|A_a\|_{\infty} \sum_{a=1}^d \|\partial_a q(0, \cdot) - \partial_a M_0 f^0(\cdot)\|_0 + C_T k $$
that is (13), since (14) and (9) imply $M_0f^0 = q^0$.

<figure><img src="image_3.png" /></figure>

The proof above shows the consistency of kinetic schemes under the assumption (*) to ensure the “luxurious” $H^s$ stability ($s = 2$) of the schemes (namely $\max_{a,b} \|\partial_{ab}^2 M_0 f^n\|_0 \le \max_{a,b} \|\partial_{ab}^2 q^0\|_0$). That condition (*) is too strong : combined with (9) it leads to $\bar{f}_{\zeta,l}(q) = \sum_m \Omega_{\zeta,m}^l q_m$ with $\sum_{\zeta} \sum_m |\Omega_{\zeta,m}^l| \in (0,1)$ for all $l$ which strongly limits the fluxes captured by (10). That is why we do not claim here to have proved any (interesting) convergence result about kinetic scheme and we do not proceed to a (cumbersome) extension of Prop. 1 to the case $\omega = 2$. However, the analysis above can guide one further in the construction of kinetic schemes. Indeed note that the approximation of (24) by (23) is in fact one order higher than (25) (one can Taylor expand one order higher (24) and (23), next get further cancellation using (10)). Thus, a higher-order approximation than the first-order rate (13) can be reached when $\omega = 1$: it suffices that the $\frac{d(d+1)}{2}$ conditions

$$ (27) \quad M_{ab}\bar{f}_l(q) := \sum_{\zeta=0}^{Q-1} c_{\zeta}^a c_{\zeta}^b \bar{f}_{\zeta,l}(q) = \sum_{m,m'=1}^{L} A_a^{lm} A_b^{mm'} q_{m'} \quad \forall a, b = 1 \dots d $$

be satisfied so as to cancel the first remaining term in the difference of the Taylor expansions (19) and (18). Then one may want to require the $\frac{d(d+1)}{2}$ conditions in the construction of convergent kinetic schemes. In fact, satisfying the $\frac{d(d+1)}{2}$ conditions is impossible in general : see e.g. the next Section 3 where a practical application case is investigated. From the analysis of the main error term, a sensible choice to construct a kinetic scheme is then : minimize

$$ (28) \quad B_{ab,l}(q) := M_{ab}\bar{f}_l(q) - \sum_{m,m'=1}^{L} A_a^{lm} A_b^{mm'} q_{m'} $$

in some norm on $\mathbb{R}^{\frac{d(d+1)}{2}}$. Such a requirement for the construction of kinetic schemes seems new to us in the literature, although the matrix $B(q) \in \mathbb{R}^{d \times L, d \times L}$ has already appeared (for instance, in [1, Prop. 1], the positivity of $B(q)$ is shown a sufficient condition for the diffusivity of the equation formally satisfied by the approximation $M_0f$ of $q$). In Section 4, we next investigate numerically the interest of such a requirement.

On the other hand, to construct kinetic schemes, one would also like to ensure sufficient stability conditions for the scheme (6) to converge. Now, other stability conditions less strong than (*) can be found in the literature e.g. [6, 7, 8, 15, 16] which have already proved useful numerically in the case of nonlinear fluxes (even though rigorous convergence has not yet been proved to our knowledge), and which turn out compatible with (9) and (10). Precisely, denoting $\eta^*(\varphi) := \sup_q (\varphi^T q - \eta(q))$ the Legendre transform of $\eta$, a $C^2(\mathbb{R}^l)$ strictly convex functional of the entropy variable $\varphi^l(q) := \partial_{q^l} \eta(q)$ while $\partial_{q^l} \eta$ is a bijection onto $\mathbb{R}^L$, a “H-theorem” can further constrain $\bar{f}$ [6, 7] :

**Proposition 2.** Assume there exist $C^1$ convex functionals $h_\zeta^*$ of $\varphi$ such that

$$ (29) \quad \sum_{\zeta=0}^{Q-1} h_{\zeta}^{*}(\varphi) = \eta^{*}(\varphi) \quad \sum_{\zeta=0}^{Q-1} c_{\zeta}^{a} h_{\zeta}^{*}(\varphi) = G_{a}^{*}(\varphi) $$
hold, using $G_a^*(\varphi) := \varphi^l F_a^l(q(\varphi)) - G_a(q(\varphi))$. Then (9)-(10) hold on requiring

$$ (30) \qquad \bar{f}_{\zeta,l}(q) = \frac{\partial h_{\zeta}^{*}}{\partial \varphi^{l}}(\varphi(q)), $$

as well as, on denoting $h_{\zeta}(f_{\zeta}) = \sup_{\varphi} (\varphi^l f_{\zeta,l} - h_{\zeta}^*(\varphi))$ for solutions to (3):

$$ (31) \qquad \partial_t \left( \sum_{\zeta=0}^{Q-1} h_{\zeta}(f) \right) + \sum_{a=1}^{d} \partial_a \left( \sum_{\zeta=0}^{Q-1} c_{\zeta}^{a} h_{\zeta}(f) \right) \le 0. $$

The proof of Prop. 2 is straightforward by direct computations, follow e.g. [16] (in the present smooth case). The key point in (31) is that whatever $\zeta$, for all $l$ it holds $\partial_{f_{\zeta,l}} h_{\zeta}(f) = \varphi^l$ for all $\zeta$, with $\varphi$ such that $f_{\zeta,l} = \partial_{\varphi,l} h_{\zeta}^*(\varphi)$, in particular $\partial_{f_{\zeta,l}} h_{\zeta}(\bar{f}(q)) = \varphi^l(q) = \partial_{q^l} \eta(q)$ whatever $q$.

The importance of requiring (a discrete equivalent to) Prop. 2 for numerical stability does not seem well established yet. But it can quite easily enforced here, see Section 3, and we shall therefore consider it here for the practical construction of kinetic schemes. Precisely, for linear symmetric-hyperbolic systems of conservation laws (1) such that $q \in \mathbb{R}^L$ are actually the variables after symmetrization (i.e. $F_a^l(q) = A_a^{lm} q_m$ with $A_a^{lm} = A_a^{ml}$ for all $a = 1 \dots d$, $1 \le m, l \le L$) in particular, we require linear Maxwellians $\bar{f}_{\zeta,l}(q) = \sum_{m=1}^L \Omega_{\zeta,l}^m q_m$ that satisfy (9) and (10) with symmetric positive matrices $\Omega_{\zeta} \in \mathbb{R}^{L \times L}$ for all $\zeta = 0 \dots Q-1$, further denoting $[\Omega_{\zeta}]_{lm} = \Omega_{\zeta,l}^m$ their entries. On noting (30), this is equivalent to require $h_{\zeta}(f_{\zeta}) = \frac{1}{2} \sum_{l,m=1}^L f_{\zeta,l} [\Omega_{\zeta}^{-1}]_{lm} f_{\zeta,m}$ where $\Omega_{\zeta}^{-1}$ denotes the matrix inverse of $\Omega_{\zeta}$ (standard inverse if $\Omega_{\zeta}$ is strictly positive equiv. positive definite, or generalized well-defined by e.g. SVD in the semi-definite case). Then indeed, the following discrete version of Prop. 2 holds:

**Proposition 3.** Assume $h_{\zeta}(f_{\zeta}) = \frac{1}{2} f_{\zeta,l} [\Omega_{\zeta}^{-1}]_{lm} f_{\zeta,m}$ then it holds for the kinetic scheme (6) if $\omega = 1$, or with equality if $\omega = 2$:

$$ \sum_{\zeta=0}^{Q-1} h_{\zeta}(f_{\zeta}^{n+1}) + \int_{t^n}^{t^{n+1}} \sum_{a=1}^{d} \partial_a \left( \sum_{\zeta=0}^{Q-1} c_{\zeta}^{a} h_{\zeta}(f_{\zeta}) \right) ds \le \sum_{\zeta=0}^{Q-1} h_{\zeta}(f_{\zeta}^{n}) $$

where $f_{\zeta}$ is the solution to the transport step of $f^{n+\frac{1}{2}}$ to $f^{n+1}$.

**Démonstration.** The transport step of $f^{n+\frac{1}{2}}$ to $f^{n+1}$ satisfies

$$ \partial_t h_{\zeta}(f_{\zeta}) + \sum_{a=1}^{d} \partial_a (c_{\zeta}^{a} h_{\zeta}(f_{\zeta})) = 0 $$

thus

$$ h_{\zeta}(f_{\zeta}^{n+1}) + \int_{t^n}^{t^{n+1}} \sum_{a=1}^{d} \partial_a (c_{\zeta}^{a} h_{\zeta}(f_{\zeta})) ds = h_{\zeta}(f_{\zeta}^{n+\frac{1}{2}}) $$

while, using the convexity of $h_{\zeta}$ and $f_{\zeta,l}^{n+\frac{1}{2}} = \sum_{m'=1}^L [\Omega_{\zeta}]_{lm'} M_0 f_{m'}^n$, it holds:

$$ \begin{aligned} & (32) \quad h_{\zeta}(f_{\zeta}^{n+\frac{1}{2}}) - h_{\zeta}(f_{\zeta}^{n}) \le \partial_{f_{\zeta}} h_{\zeta}(f_{\zeta}^{n+\frac{1}{2}}) \cdot (f_{\zeta}^{n+\frac{1}{2}} - f_{\zeta}^{n}) \\ & = \sum_{m=1}^{L} f_{\zeta,l}^{n+\frac{1}{2}} [\Omega_{\zeta}^{-1}]_{lm} \left( \sum_{m'=1}^{L} [\Omega_{\zeta}]_{mm'} M_0 f_{m'}^n - f_{\zeta,m}^n \right) = \sum_{m=1}^{L} M_0 f_m^n \left( \sum_{m'=1}^{L} [\Omega_{\zeta}]_{mm'} M_0 f_{m'}^n - f_{\zeta,m}^n \right) \end{aligned} $$
where the last term in (32) satisfies

$$ (33) \qquad \sum_{\zeta} \left( \sum_{m'=1}^{L} [\Omega_{\zeta}]_{mm'} M_0 f_{m'}^n - f_{\zeta,m}^n \right) = 0 $$

which concludes the case $\omega = 1$. On the other hand, using $x^2 - y^2 = (x+y)(x-y)$ and $f_{\zeta,l}^{n+\frac{1}{2}} + f_{\zeta,l}^n = 2 \sum_{m'=1}^L [\Omega_\zeta]_{lm'} M_0 f_{m'}^n$, it holds when $\omega = 2$:

$$ (34) \qquad \begin{aligned} h_{\zeta}(f_{\zeta}^{n+\frac{1}{2}}) - h_{\zeta}(f_{\zeta}^n) &= \frac{1}{2} \Omega_{\zeta}^{-1} (f_{\zeta,l}^{n+\frac{1}{2}} + f_{\zeta,l}^n) \cdot (f_{\zeta}^{n+\frac{1}{2}} - f_{\zeta}^n) \\ &= \sum_{m=1}^{L} M_0 f_m^n \left( \sum_{m'=1}^{L} [\Omega_{\zeta}]_{mm'} M_0 f_{m'}^n - f_{\zeta,m}^n \right) \end{aligned} $$

and one again finishes the proof on noting that the last term in (34) satisfies (33). □

<figure><img src="image_4.png" /></figure>

A similar entropy stability as in Prop. 3 is proved in [20, Chapter 6] for the D1Q2 kinetic scheme applied to the scalar transport equation with $\omega \in (1, 2)$.

### 3. CONSTRUCTION OF KINETIC SCHEMES FOR 2D ELASTODYNAMICS

The 2-dimensional linear (Hookean) elastodynamics system reads

$$ (35) \qquad \begin{cases} \partial_t \mathbb{F}_a^m - \partial_a u^m = 0 \\ \partial_t u^m - c^2 \partial_a \mathbb{F}_a^m = 0 \end{cases}. $$

It governs the position $\sum_{m=1}^2 x_m(t, \sum_{a=1}^2 a_a E^a) e^m$ of a 2-dimensional body, in the Euclidean physical space equipped with a Cartesian coordinate system $\{e^1, e^2\}$, using Lagrangian coordinates in a Cartesian system $\{E^1, E^2\}$ equipping a reference configuration of the body, assuming $\mathbb{F}_a^m = \partial_a x^m$, $u^m = \partial_t x^m$.

It is a paradigmatic *Friedrichs symmetrizable linear system of conservation laws*:

$$ q = (\sqrt{c}\mathbb{F}_1^1, \sqrt{c}\mathbb{F}_2^1, u^1/\sqrt{c}, \sqrt{c}\mathbb{F}_1^2, \sqrt{c}\mathbb{F}_2^2, u^2/\sqrt{c}) $$

is solution to a symmetric-hyperbolic system of conservation laws. Noticeably, the linear homogeneous system (35), symmetrized in $q$, can be interpreted as a vectorial version of the (symmetrized) 2D linear acoustics system for $r = (v, w, p) \in \mathbb{R}^3$

$$ (36) \qquad \partial_t r + A_1 \partial_1 r + A_2 \partial_2 r = 0 $$

that governs velocities $\sum_{a=1}^2 \sqrt{c} F_a^m E^a \sim vE^1 + wE^2$ and pressures $u^m/\sqrt{c} \sim p$ in two fluids $m = 1, 2$ with (same) sound celerity $c > 0$ parametrizing

$$ A_1 = \begin{pmatrix} 0 & 0 & -c \\ 0 & 0 & 0 \\ -c & 0 & 0 \end{pmatrix}, \quad A_2 = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & -c \\ 0 & -c & 0 \end{pmatrix}. $$

Here, we develop our construction methodology of a kinetic scheme (based on elements from the previous sections) specifically for the symmetric system on 2D Cartesian grids with grid-size $\Delta > 0$. We consider either D2Q5 or D2Q9 sets of speeds $c_\zeta$ (see details below). We require *linear* Maxwellians $\tilde{f}_{\zeta,l}(q) = \Omega_{\zeta,l}^m q_m$ satisfying (9), (10) and (29), which is natural insofar as $\eta^*, G_1^*, G_2^*$ are quadratic functionals of $\varphi$. Moreover we require that the kinetic scheme preserves the spatial
isotropy of the 2D elastodynamics (i.e. ($F_1^2, F_2^2, u^2$) satisfies exactly the same dynamics as ($F_1^1, F_2^1, u^1$) when initial conditions are the same). So we look only for decompositions (29) of the form

$$ (37) \qquad \partial_{\varphi^l} h_{1,\zeta}^*(\varphi) = \begin{pmatrix} \Omega_\zeta & 0 \\ 0 & \Omega_\zeta \end{pmatrix} q^T $$

using symmetric-positive matrices $\Omega_\zeta \in \mathbb{R}^{3 \times 3}$ which obviously ensure :

**Proposition 4.** The numerical approximation $\sum_\zeta f_{\zeta,l}, l \in \{1, 2, 3\}$, to $(cF_1^1, cF_2^1, u^1/c)$ satisfies the same dynamics as the numerical approximation $\sum_\zeta f_{\zeta,l}, l \in \{4, 5, 6\}$, to $(cF_1^2, cF_2^2, u^2/c)$ when initial conditions are the same, and (37) holds.

We consider classical sets of kinetic speeds : either the standard D2Q5 case :

$$ c_0 = (0,0); \quad c_1 = \lambda(1,0); \quad c_2 = \lambda(0,1); \quad c_3 = \lambda(-1,0); \quad c_4 = \lambda(0,-1) $$

with $Q = 5$ speeds, or the standard D2Q9 case with four additional speeds :

$$ c_5 = \lambda(1,1); \quad c_6 = \lambda(-1,1); \quad c_7 = \lambda(-1,-1); \quad c_8 = \lambda(1,-1) $$

which both use grids with step size $\Delta = k\lambda$. Both cases possess symmetries : given one speed $c_1$ in the material space $\mathbb{R}^2$ equipped with the Cartesian basis $\{E^1, E^2\}$,

$$ c_\zeta = R_{\pi/2}^{\zeta-1} c_1 $$

holds for $\zeta \in \{2, 3, 4\}$, and

$$ c_\zeta = R_{\pi/2}^{\zeta-5} c_5 $$

for $\zeta \in \{6, 7, 8\}$ when $Q = 9$ and $c_5 = c_1 + c_2$, denoting

$$ (38) \qquad R_{\pi/2} = \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix} $$

the (matrix representation of) $\pi/2$-rotation of $\mathbb{R}^2$. Now, for the numerical simulation on a Cartesian grid of homogeneous physics characterized by a uniform constant $c^2$, we naturally also require the preservation of the ambient-space symmetries when possible at discrete level, e.g. the $\pi/2$-rotation.

Precisely, still denoting $R_{\pi/2}$ the $3 \times 3$ matrix with $2 \times 2$ upper-left block (38) complemented by identity (and $R_{-\pi/2}$ its inverse), we require that $f_\zeta$ should satisfy exactly the same dynamics as $R_{\pi/2}^{\zeta-1} f_\zeta$ when initial conditions are also rotated :

$$ (39) \qquad \Omega_\zeta = R_{\pi/2}^{\zeta-1} \Omega_1 R_{-\pi/2}^{\zeta-1} $$

for $\zeta \in \{2, 3, 4\}$ and

$$ (40) \qquad \Omega_\zeta = R_{\pi/2}^{\zeta-5} \Omega_5 R_{-\pi/2}^{\zeta-5} $$

for $\zeta \in \{6, 7, 8\}$ when $Q = 9$.

**Proposition 5.** The numerical approximations ($\sum_\zeta f_{\zeta,l}$) to $(cF_1^1, cF_2^1, u^1/c)$ (with $l \in \{1, 2, 3\}$) and to $(cF_1^2, cF_2^2, u^2/c)$ (with $l \in \{4, 5, 6\}$) are exactly the same as $R_{-\pi/2}(\sum_\zeta f_{\zeta,l})$ after rotating the initial conditions similarly.
### 3.1. D2Q5 scheme.
Applying the construction principles above to approximate linear elastodynamics with D2Q5 speeds, it remains to specify

$$ \Omega_1 = \begin{pmatrix} \alpha & x & y \\ x & \beta & z \\ y & z & \gamma \end{pmatrix} $$

insofar as $\Omega_{2,3,4}$ can next be computed from $\Omega_1$ using (39) : (41)

$$ \Omega_2 = \begin{pmatrix} \beta & -x & -z \\ -x & \alpha & y \\ -z & y & \gamma \end{pmatrix}, \quad \Omega_3 = \begin{pmatrix} \alpha & x & -y \\ x & \beta & -z \\ -y & -z & \gamma \end{pmatrix}, \quad \Omega_4 = \begin{pmatrix} \beta & -x & z \\ -x & \alpha & -y \\ z & -y & \gamma \end{pmatrix}. $$

The consistency relations (29), equivalent to (9–10) with linear Maxwellians, impose

$$ \lambda(\Omega_1 - \Omega_3) = \begin{pmatrix} 0 & 0 & -c \\ 0 & 0 & 0 \\ -c & 0 & 0 \end{pmatrix}, \quad \lambda(\Omega_2 - \Omega_4) = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & -c \\ 0 & -c & 0 \end{pmatrix} $$

i.e.

$$ y = -\frac{c}{2\lambda}, \quad z = 0, $$

and

$$ (42) \qquad \Omega_0 = \begin{pmatrix} 1-2(\alpha+\beta) & 0 & 0 \\ 0 & 1-2(\alpha+\beta) & 0 \\ 0 & 0 & 1-4\gamma \end{pmatrix}. $$

So there remain infinitely many choices for a scheme like (6).

To ensure numerical stability, we moreover require Prop. 3, that is : all the matrices $\Omega_\zeta$ should be positive. Here, in view of (41) we need only check that $\Omega_1$ and $\Omega_0$ are positive, i.e. all their principal minors are positive :

$$ (43) \quad \alpha\beta \ge x^2, \quad \alpha\gamma \ge \frac{c^2}{4\lambda^2}, \quad \alpha, \beta, \gamma \ge 0, \quad \alpha + \beta \le 1/2, \quad \gamma \le 1/4. $$

It is remarkable that here, for a 3 × 3 SCL in two dimensions, the “(kinetic) entropy” conditions (43) imply the positivity of the symmetric matrix $B(q)$ defined in (28), which specifically reads in $\mathbb{R}^6$ here as

$$ \left[ \frac{M_{11}\bar{f} - A_1 A_1}{-A_2 A_1} \quad \frac{-A_1 A_2}{M_{22}f - A_2 A_2} \right]. $$

Indeed, our conditions (43) imply

$$ \alpha\gamma \ge \frac{c^2}{4\lambda^2}, \quad \frac{1}{4\gamma} \ge 1 \Rightarrow \alpha \ge \frac{c^2}{\lambda^2} $$

and

$$ \alpha\gamma \ge \frac{c^2}{4\lambda^2}, \quad \frac{1}{2\alpha} \ge 1 \Rightarrow \gamma \ge \frac{c^2}{2\lambda^2} $$

i.e. $B(q) \ge 0$. Now, although established for (discrete velocity) kinetic approximations of scalar equations [1, 20], that implication has not yet been proved generically, for kinetic approximations of any (symmetric-hyperbolic, multidimensional) SCL. Recall that one consequence of the positivity of $B(q)$ is that the Chapman-Enskog expansion of the (discrete velocity) kinetic scheme is actually a diffusive approximation of the target SCL, which is a desired stability requirement, see [1].
One noticeable choice then corresponds to optimizing the accuracy of the scheme on aiming at (27) as “second-order” consistency conditions (equivalent to minimum numerical diffusion in the standard analysis by equivalent equations), where

$$M_{11}\bar{f} \equiv \lambda^2(\Omega_1 + \Omega_3) = 2\lambda^2 \begin{pmatrix} \alpha & x & 0 \\ x & \beta & 0 \\ 0 & 0 & \gamma \end{pmatrix} \qquad M_{12}\bar{f} = 0 = M_{21}\bar{f}$$

$$M_{22}\bar{f} \equiv \lambda^2(\Omega_2 + \Omega_4) = 2\lambda^2 \begin{pmatrix} \beta & -x & 0 \\ -x & \alpha & 0 \\ 0 & 0 & \gamma \end{pmatrix} \qquad A_1 A_1 = \begin{pmatrix} c^2 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & c^2 \end{pmatrix}$$

$$A_2 A_2 = \begin{pmatrix} 0 & 0 & 0 \\ 0 & c^2 & 0 \\ 0 & 0 & c^2 \end{pmatrix} \qquad A_1 A_2 = \begin{pmatrix} 0 & c^2 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} \qquad A_2 A_1 = \begin{pmatrix} 0 & 0 & 0 \\ c^2 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}$$

The “second-order” conditions (27) cannot be fully fullfilled, because $M_{12}\bar{f} = 0 = M_{21}\bar{f}$. But one can minimize numerical diffusion using $\beta = 0 = x$ and e.g.

$$\underset{\{0 \le \alpha \le 1/2\}}{\operatorname{argmin}} \quad \underset{\{0 \le \gamma \le 1/4\}}{\lambda^2 \ge \frac{c^2}{4\alpha\gamma}} \quad J(\alpha, \gamma, \lambda) := (|2\lambda^2\alpha - c^2|^2 + |2\lambda^2\gamma - c^2|^2)$$

that is,

$$(44) \qquad \gamma = \frac{1}{4} \qquad \alpha = \frac{1}{2} \qquad \lambda = \sqrt{2c}$$

on noting $J$ is smooth, coercive, constrained on bounded compact domain, with extrema reached on boundaries $\alpha = \frac{1}{2}$, $\gamma = \frac{1}{4}$ using “Lagrange multiplyers” $\mu, \nu$

$$\gamma\lambda^2 = \frac{c^2}{2} + \frac{c^2\mu}{\lambda^2} \qquad \alpha\lambda^2 = \frac{c^2}{2} + \frac{c^2\nu}{\lambda^2}$$

compatible with $\lambda \ge \sqrt{2c}$. In the sequel, we shall compare (44) with a more standard CFL $\frac{1}{2}$, i.e. $\lambda = 2c$ (so $\alpha\gamma \ge 1/16$), and three extremal values for $\alpha, \gamma$:

$$\gamma = \frac{1}{4} \qquad \alpha = \frac{1}{2},$$

$$\text{or} \quad \gamma = \frac{1}{4} \qquad \alpha = \frac{1}{4},$$

$$\text{and also} \quad \gamma = \frac{1}{8} \qquad \alpha = \frac{1}{2}.$$

Note that when $\gamma = \frac{1}{4}$, $\alpha = \frac{1}{2}$, the D2Q5 scheme is in fact a D2Q4 since $\Omega_0 = 0$.

3.2. **D2Q9 scheme.** Applying the construction principles above with D2Q9 speeds, we have only to specify

$$\Omega_1 = \begin{pmatrix} \alpha & x & y \\ x & \beta & z \\ y & z & \gamma \end{pmatrix} \qquad \Omega_5 = \begin{pmatrix} \tilde{\alpha} & \tilde{x} & \tilde{y} \\ \tilde{x} & \tilde{\beta} & \tilde{z} \\ \tilde{y} & \tilde{z} & \tilde{\gamma} \end{pmatrix}$$
but not less since it does not seem obvious to link $\Omega_1$ and $\Omega_5$ through symmetries of the Cartesian grid. The consistency relations (29) (or (9) and (10)) impose

$$ (45) \quad \lambda (\Omega_1 - \Omega_3) + \sqrt{2}\lambda (\Omega_5 - \Omega_6 - \Omega_7 + \Omega_8) = \begin{pmatrix} 0 & 0 & -c \\ 0 & 0 & 0 \\ -c & 0 & 0 \end{pmatrix} $$ 
 $$ \lambda (\Omega_2 - \Omega_4) + \sqrt{2}\lambda (\Omega_5 + \Omega_6 - \Omega_7 - \Omega_8) = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & -c \\ 0 & -c & 0 \end{pmatrix} $$

i.e.

$$ (46) \quad 2\lambda(y + \sqrt{2}\tilde{y} - \sqrt{2}\tilde{z}) = -c $$

$$ (47) \quad 2\lambda(z + \sqrt{2}\tilde{z} + \sqrt{2}\tilde{y}) = 0 $$

$$ (48) \quad 2\lambda(y + \sqrt{2}\tilde{z} - \sqrt{2}\tilde{y}) = -c $$

$$ (49) \quad 2\lambda(z + \sqrt{2}\tilde{y} + \sqrt{2}\tilde{z}) = 0 $$

so

$$ \tilde{y} = \tilde{z} = 0 = z \quad y = -\frac{c}{2\lambda} $$

and

$$ (50) \quad \Omega_0 = \begin{pmatrix} 1 - 2(\alpha + \beta + \tilde{\alpha} + \tilde{\beta}) & 0 & 0 \\ 0 & 1 - 2(\alpha + \beta + \tilde{\alpha} + \tilde{\beta}) & 0 \\ 0 & 0 & 1 - 4(\gamma + \tilde{\gamma}) \end{pmatrix} $$

Furthermore, the positivity of $\Omega_0$, $\Omega_1$ and $\Omega_5$ requires :

$$ (51) \quad \alpha\beta \ge x^2, \quad \alpha\gamma \ge \frac{c^2}{4\lambda^2}, \quad \tilde{\alpha}\tilde{\beta} \ge \tilde{x}^2, \quad \tilde{\alpha}\tilde{\gamma} \ge 0, \\ \qquad \tilde{\alpha}, \tilde{\beta}, \tilde{\gamma} \ge 0, \alpha, \beta, \gamma \ge 0, \alpha + \beta + \tilde{\alpha} + \tilde{\beta} \le 1/2, \gamma + \tilde{\gamma} \le 1/4. $$

Finally, on computing

$$ M_{11}\bar{f} \equiv \lambda^2(\Omega_1+\Omega_3)+2\lambda^2(\Omega_5+\Omega_6+\Omega_7+\Omega_8) = 2\lambda^2 \begin{pmatrix} \alpha+2(\tilde{\alpha}+\tilde{\beta}) & x & 0 \\ x & \beta+2(\tilde{\alpha}+\tilde{\beta}) & 0 \\ 0 & 0 & \gamma+2\tilde{\gamma} \end{pmatrix} $$

$$ M_{12}\bar{f} = M_{21}\bar{f} \equiv 2\lambda^2(\Omega_5 - \Omega_6 + \Omega_7 - \Omega_8) = 2\lambda^2 \begin{pmatrix} 2(\tilde{\alpha} - \tilde{\beta}) & 4\tilde{x} & 0 \\ 4\tilde{x} & 2(\tilde{\beta} - \tilde{\alpha}) & 0 \\ 0 & 0 & 0 \end{pmatrix} $$

$$ M_{22}\bar{f} \equiv \lambda^2(\Omega_2+\Omega_4)+2\lambda^2(\Omega_5+\Omega_6+\Omega_7+\Omega_8) = 2\lambda^2 \begin{pmatrix} \beta+2(\tilde{\alpha}+\tilde{\beta}) & -x & 0 \\ -x & \alpha+2(\tilde{\alpha}+\tilde{\beta}) & 0 \\ 0 & 0 & \gamma+2\tilde{\gamma} \end{pmatrix} $$

minimizing numerical diffusion leads, using (51), to

$$ \beta = x = 0 = \tilde{\alpha} = \tilde{\beta} = \tilde{\gamma} = \tilde{x}, \quad \lambda^2 = 2c^2, \quad \gamma = \frac{1}{4}, \quad \alpha = \frac{1}{2} $$

that is the D2Q5 scheme already constructed in the previous section! Nevertheless, to see the effect of a few additional celerities (i.e. kinetic directions), we could also
do the following choice :

$$
(52) \qquad \beta = 0 \quad \gamma = \frac{1}{4} \quad \tilde{\alpha} = \tilde{\beta} = \tilde{x} = \frac{c^2}{16\lambda^2} \quad \lambda = \frac{3c}{2} \quad \text{and} \quad \alpha = \frac{4}{9}.
$$

4. NUMERICAL EVALUATION

Here we evaluate numerically our vectorial kinetic schemes built in Section 3,
for the approximation in $L^2(\mathcal{D})$, $\mathcal{D} := \{\|\mathbf{a}\| \le R\}$, at $t = T$ of solutions to (35)
computed on $\mathcal{D}_T \supset \{\|\mathbf{a}\| \le R + cT\}, t \in (0, T)$. Note that the conditions used at
the boundary $\partial \mathcal{D}_T$ are not influential insofar as the simulation domain $\mathcal{D}_T$ contains
the dependence domain of $\mathcal{D}$ at $t = T$. Simulations can be thus compared on $\mathcal{D}$
with smooth analytical solutions to (35) on $\mathbb{R}^d$, which can be defined e.g. from

$$
(53) \qquad \left( \partial_{tt}^2 - c^2 \sum_{\alpha=1}^{d} \right) u^m = 0 \quad \text{for } m = 1, 2
$$

complemented with $u^m|_{t=0}$ and $\partial_t u^m|_{t=0} = c^2 \sum_{\alpha=1}^d \partial_\alpha F_\alpha^m|_{t=0}$, also computable as

$$
(54) \quad u^m(t, \mathbf{a}) = \partial_t \frac{1}{2\pi} \int_{D_{t,a}} \frac{u^m|_{t=0}(\mathbf{b})}{\sqrt{c^2 t^2 - |\mathbf{a} - \mathbf{b}|^2}} d\mathbf{b} - \frac{1}{2\pi} \int_{D_{t,a}} \frac{\partial_t u^m|_{t=0}(\mathbf{b})}{\sqrt{c^2 t^2 - |\mathbf{a} - \mathbf{b}|^2}} d\mathbf{b}
$$

where $D_{t,a} := \{\mathbf{b}, c^2t^2 \ge |\mathbf{a} - \mathbf{b}|^2\}$ is the domain of influence on the solution at $(t, \mathbf{a})$. Focusing here on axisymmetric solutions (that depend only on $r > 0$ the distance to the origin), let us introduce the (zero-order) Hankel transform of $u^m$

$$
\widetilde{u}^m(k, t) = \int_0^\infty r J_0(kr) u^m(r, t) dr
$$

where $J_0$ is zero-order Bessel function, also often denoted $H_0(u^m)(k, t)$. It satisfies

$$
\int_{0}^{\infty} \partial_{r} (r \partial_{r} u^{m}(r, t)) J_{0}(kr) dr = -k^{2} \tilde{u}^{m}(k, t)
$$

so one can also compute solutions to (53) by

$$
(55) \quad u^m(t, r) = \int_0^\infty \left( \widetilde{u^m}_{|t=0}(k) \cos(ckt) + \partial_t \widetilde{u^m}_{|t=0}(k) \frac{\sin(ckt)}{c} \right) J_0(kr) kdk.
$$

For instance, recalling $r^2 = \|\mathbf{a}\|^2$, choosing $\partial_t u^m|_{t=0} = 0 = \sum_{\alpha=1}^d \partial_\alpha F_\alpha^m|_{t=0}$ and
$u^m|_{t=0}(r) = \kappa e^{-\mu r^2}$ yields [13, Chap.7] $\widetilde{u^m}|_{t=0}(k) = \frac{\kappa}{2\mu} e^{-\frac{k^2}{4\mu}}$ so

$$
(56) \quad u^m(t, r) = \frac{\kappa}{2\mu} \int_0^{+\infty} e^{-\frac{k^2}{4\mu}} \cos(ckt) J_0(kr) kdk = \mathcal{H}_0 \left( \frac{\kappa}{2\mu} e^{-\frac{k^2}{4\mu}} \cos(ckt) \right) (t, r).
$$

The latter analytic solution can be numerically computed with high precision given
$r$, see e.g. the solution plot in Fig. 1 using the open-source gsl<fnref n="1" />.

For benchmarking purposes, we also compute numerical approximations using
other discretization methods on the same uniform Cartesian grid of $\mathcal{D}_T$, with square
cells $\Delta \times \Delta$ centered at $\mathbf{a}_{ij}$, $ij \in \mathcal{I} := \{-I, -I+1, \dots, I\}$ ($I \in \mathbb{N}$), ensuring a finite
propagation speed of the information through *time-explicit* schemes. We impose as
initial condition $u_0(r) = \kappa e^{-\mu r^2}$ for $\kappa = 1, \mu = 2$, with $r = \sqrt{x_i^2 + y_j^2}$ at cell centers
$(x_i, y_j) \in \mathcal{D}_T := (-4, 4)^2$ in the FV case (i.e. simply grid points in the FD case).

<footnote n="1">1. The discrete Hankel transform, with $M = 2^{12}$ zeros of $J_0$ and $X = 4$ precisely in Fig. 1,
using the notations of https://www.gnu.org/software/gsl/doc/html/dht.html</footnote>
<figure><img src="image_5.png" /><figcaption>FIGURE 1. Solution (56) for $\kappa = 1, \mu = 2$ computed at $t \in \{0; 0.5; 1; 2\}$ as a function of $r$ using gsl.</figcaption></figure>

A mesh convergence study is plotted on Fig. 2 when $R = 2$, $c = \frac{1}{\sqrt{2}}$ at $T = 1$. Our proposed D2Q5 vectorial kinetic schemes are compared one another, when $\omega = 1$ and when $\omega = 2$, with first and second order FV schemes detailed in Appendix A.2.

First, our numerical experiments show the convergence of the vectorial kinetic schemes, at rates $O(\Delta)$ and $O(\Delta^2)$ respectively when $\omega = 1$ and when $\omega = 2$, whenever the parameters $\lambda, \alpha, \gamma$ satisfy (43). It is remarkable that when only the less stringent conditions $\alpha \ge \frac{c^2}{\lambda^2}, \gamma \ge \frac{c^2}{2\lambda^2}$ (so $B(q) \ge 0$) are satisfied strictly (i.e. $\lambda < c\sqrt{2}$ and (43) are not all satisfied) then we observe numerical instabilities. For instance, our D2Q5 schemes all blow up with $\Delta = 2^{-9} \equiv .001953125$ when $\lambda = c(\sqrt{2} - .1)$ whatever $\omega \in (1, 2)$, and for $\omega \ge 1.9$ when $\lambda = c(\sqrt{2} - .01)$. Such numerical results are not surprising in view of previous results [1, 20]. A full rigorous theory regarding the necessary/sufficient conditions for convergence of the kinetic schemes toward SCL however remains needed.

Second, we observe that the kinetic schemes with largest $CFL = 1/\sqrt{2}$, i.e. using $\lambda = \sqrt{2c}, \gamma = 1/4, \alpha = 1/2$ are the most accurate at both first and second order ($\omega = 1, \omega = 2$). That parameter choice was motivated by minimizing the numerical diffusion when $\omega = 1$ (recall the proof of Prop. 1). The D2Q9 schemes using a slightly more diffusive CFL coefficient $\lambda = 3c/2$ from (52) are second best, next come the cases $\lambda = 2c$ corresponding to a 1/2 CFL condition. Moreover, when simulated with $\omega = 1$ or 2 in the stability region of Prop. 3, different (admissible) values of $\alpha, \gamma$ have remarkable impact on the accuracy of numerical solution. This also calls for further studies (of the discretization error), theoretically and numerically (note that $\beta \ne 0, x \ne 0$ in the Maxwellians has not been tested here.)

Interestingly, we observe that the convergences of the first order FV scheme and that of the first order kinetic scheme with $\lambda = 2c, \gamma = 1/4 = \alpha$ are superimposed,
which is not surprising insofar as, recalling $\lambda = \frac{\Delta}{k}$, the scheme

$$q_{i,j}^{n+\frac{1}{2}} = \Omega_0 q_{i,j}^{n-\frac{1}{2}} + \Omega_1 q_{i-1,j}^{n-\frac{1}{2}} + \Omega_2 q_{i,j-1}^{n-\frac{1}{2}} + \Omega_3 q_{i+1,j}^{n-\frac{1}{2}} + \Omega_4 q_{i,j+1}^{n-\frac{1}{2}}$$

that results from LBM for the first-order moments when $\lambda = 2c$ i.e.

$$\Omega_1 = \begin{pmatrix} \frac{1}{4} & 0 & -\frac{1}{4} \\ 0 & 0 & 0 \\ -\frac{1}{4} & 0 & \frac{1}{4} \end{pmatrix} \quad \Omega_0 = \begin{pmatrix} \frac{1}{4} & 0 & 0 \\ 0 & \frac{1}{4} & 0 \\ 0 & 0 & 0 \end{pmatrix}$$

then exactly coincides with the first-order FV scheme at CFL $\frac{1}{2}$, see (62–63–64) in Appendix A.2. Similarly, the second-order numerical results for the Yee scheme look very much like those for the kinetic scheme with $\lambda = 2c, \gamma = 1/8, \alpha = 1/2$.

<figure><img src="image_6.png" /><figcaption>FIGURE 2. Errors in $L_2(\{\|a\| \le 2\})$ of approximations of (56) at $t = 1$ when $\kappa = 1, \mu = 2, c = \frac{1}{\sqrt{2}}$, by (i) various D2Q5, D2Q9 kinetic schemes, with various values of parameters $\alpha, \gamma, \lambda, \omega$, and (ii) the FD, FV schemes with CFL = 1/2 – on $D_T = (-4, 4)^2$ –.</figcaption></figure>

## ACKNOWLEDGMENTS

We would like to acknowledge the help of ANR with grant ANR-15-CE01-0013 for the SEDIFLO project. We also thank Lucas Brélivet for its help in a preliminary version of our numerical implementation with basilisk<fnref n="2" />.

<footnote n="2">2. http://basilisk.fr/</footnote>
ANNEXE A. FINITE DIFFERENCE AND FINITE VOLUME SCHEMES

A.1. Finite-Difference scheme. Each “acoustic” sub-system $(cF_1^m, cF_2^m, u^m) \equiv (v, w, p)$ for $m = 1, 2$ can be numerically simulated using a standard (leap-frog, centered) Finite-Difference scheme, often named after Yee [25]. Denoting staggered nodal approximations

$$p_{ij}^n \approx p(t^n, i\Delta, j\Delta) \quad w_{i\Delta+\frac{1}{2},j}^{n+\frac{1}{2}} \approx w(t^n+\frac{k}{2}, i\Delta, j\Delta+\frac{\Delta}{2}) \quad v_{i+\frac{1}{2},j}^{n+\frac{1}{2}} \approx w(t^n+\frac{k}{2}, i\Delta+\frac{\Delta}{2}, j\Delta)$$

the following Yee scheme writes

$$(57) \qquad v_{i+\frac{1}{2},j}^{n+\frac{1}{2}} = v_{i+\frac{1}{2},j}^{n-\frac{1}{2}} + \frac{ck}{\Delta}(p_{i+1,j}^n - p_{ij}^n)$$

$$(58) \qquad w_{i\Delta+\frac{1}{2},j}^{n+\frac{1}{2}} = w_{i\Delta+\frac{1}{2},j}^{n-\frac{1}{2}} + \frac{ck}{\Delta}(p_{i\Delta+1,j}^n - p_{ij}^n)$$

$$(59) \qquad p_{ij}^{n+1} = p_{ij}^{n} + \frac{ck}{\Delta}(v_{i+\frac{1}{2},j}^{n+\frac{1}{2}} - v_{i-\frac{1}{2},j}^{n+\frac{1}{2}} + w_{i\Delta+\frac{1}{2},j}^{n+\frac{1}{2}} - w_{i\Delta-\frac{1}{2},j}^{n+\frac{1}{2}})$$

yields approximations such as the space-continuous reconstructions

$$\sum_{ij \in \mathcal{I}} p_{ij}^n \psi_{ij} \approx p(t^n) \quad \sum_{ij \in \mathcal{I}} w_{i\Delta+\frac{1}{2},j}^n \psi_{ij+\frac{1}{2}} \approx w(t^n) \quad \sum_{ij \in \mathcal{I}} v_{i+\frac{1}{2},j}^n \psi_{i+\frac{1}{2},j} \approx v(t^n)$$

using the tensor-product $\mathbb{Q}_{1,2}$ Finite-Element basis functions $\psi_{ij}, \psi_{i+\frac{1}{2},j}, \psi_{i+\frac{1}{2},j}$. The latter converges in $L^2(\mathcal{D})$ to smooth $(v, w, p) \in H^2(\mathbb{R}^2)$ with an error rate $O(\Delta^2)$ as $\Delta = ck \to 0$ [24]. A proof for $p$ consists in bounding above the RHS of

$$(60) \quad \|p(t^n) - \sum_{ij \in \mathcal{I}} p_{ij}^n \psi_{ij}\|_{L^2(\mathcal{D})} \le \|p(t^n) - \sum_{ij \in \mathcal{I}} p(t^n, i\Delta, j\Delta)\psi_{ij}\|_{L^2(\mathcal{D})} \\ \qquad + \sup_{ij \in \mathcal{I}} \|\psi_{ij}\|_{L^2(\mathcal{D})} \sqrt{\sum_{ij \in \mathcal{I}} |p(t^n, i\Delta, j\Delta) - p_{ij}^n|^2}$$

by the sum of an upper-bound for each term : the first term in RHS of (60) can be bounded above by $C_1\Delta^2\|p(t^n)\|_{H^2(\mathcal{D})}$ with a constant $C_1$ function of $p$ but independent of $\Delta$ using a standard interpolation result and the second term in RHS of (60) can be bounded above by $C_2\Delta^2 \sup_{t \in [0, t^n]} \|p(t)\|_{H^2(\mathcal{D})}$ using a standard a priori analysis for the pointwise errors like $p(t^n, i\Delta, j\Delta) - p_{ij}^n$ which are solution to (57–58–59) with additional error terms due to Taylor approximation (the bound holds thanks to the discrete $\ell^2$ stability of the scheme and the consistency of the initial conditions). Smoothness of the reference solution is key to the error estimation.

A.2. Finite-Volume methods. Explicit-in-time Finite-Volume (FV) methods (using time discretization like the Euler forward method) are also standardly used to discretize (35). Unlike the FD method above, they provide one with co-located cell-averaged approximations $q_{ij}^n$ for $q(t^n, x, y) = (F_1^m, F_2^m, u^m)$ over the control volumes $C_{ij} = (i\Delta \pm \frac{\Delta}{2}, j\Delta \pm \frac{\Delta}{2})$ defined by

$$q_{ij}^n \approx \frac{1}{\Delta^2} \int_{C_{ij}} q(t^n, x, y) dxdy.$$

The generic FV scheme writes

$$(61) \qquad q_{ij}^{n+1} = q_{ij}^{n} + \frac{k}{\Delta}(\mathcal{F}_{i,j-\frac{1}{2}} + \mathcal{F}_{i-\frac{1}{2},j} - \mathcal{F}_{i,j+\frac{1}{2}} - \mathcal{F}_{i+\frac{1}{2},j})$$
which requires formulas for the numerical fluxes $\mathcal{F}_{i\pm\frac{1}{2},j}$ and $\mathcal{F}_{i,j\pm\frac{1}{2}}$ approximating the fluxes $F_a(q)$ of the conservation laws (1) at cell faces

$$ \mathcal{F}_{i\pm\frac{1}{2},j} \approx \frac{1}{k\Delta} \int_{t^n}^{t^{n+1}} \int_{y_{j-\frac{1}{2}}}^{y_{j+\frac{1}{2}}} F_1(q(x_{i\pm\frac{1}{2}}, y, t)) dydt, $$

$$ \mathcal{F}_{i,j\pm\frac{1}{2}} \approx \frac{1}{k\Delta} \int_{t^n}^{t^{n+1}} \int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} F_2(q(x, y_{j\pm\frac{1}{2}}, t)) dxdt. $$

A standard Two-Points (TP) flux is the upwind flux which yields for each “acoustic” sub-system with variable $(F_1^m, F_2^m, u^m) \equiv (v, w, p), m \in \{1, 2\}$

$$ (62) \qquad v_{ij}^{n+1} = v_{ij}^n + \frac{kc}{2\Delta}(p_{i+1,j}^n - p_{i-1,j}^n - 2v_{ij}^n + v_{i+1,j}^n + v_{i-1,j}^n) $$

$$ (63) \qquad w_{ij}^{n+1} = w_{ij}^n + \frac{kc}{2\Delta}(p_{ij+1}^n - p_{ij-1}^n - 2w_{ij}^n + w_{ij+1}^n + w_{ij-1}^n) $$

$$ (64) \qquad \begin{aligned} p_{ij}^{n+1} = p_{ij}^n + \frac{kc}{2\Delta}(v_{i+1,j}^n - v_{i-1,j}^n + w_{ij+1}^n - w_{ij-1}^n - 4p_{ij}^n + p_{ij+1}^n + p_{ij-1}^n + p_{i+1,j}^n + p_{i-1,j}^n) \end{aligned} $$

where, in comparison with the “staggered” scheme (57–59), new terms appear that are not obviously consistent but can be interpreted as artificial diffusion, small as $O(\Delta^2)$ noticeably. Now, as is well-known, the latter “non-consistent” terms do not prevent convergence : in the case of smooth solutions, the $O(\Delta^2)$ scaling of those terms even preserves for (62–64) the convergence rate observed for (57–59). A proof is as follows : Consistency errors $V_{ij}^n = v_{ij}^n - v(kn, \Delta i, \Delta j)$, $W_{ij}^n = w_{ij}^n - w(kn, \Delta i, \Delta j)$, $P_{ij}^n = p_{ij}^n - p(kn, \Delta i, \Delta j)$ satisfy

$$ (65) \qquad V_{ij}^{n+1} = V_{ij}^n + \frac{kc}{2\Delta}(P_{i+1,j}^n - P_{i-1,j}^n - 2V_{ij}^n + V_{i+1,j}^n + V_{i-1,j}^n) + R V_{ij}^n $$

$$ (66) \qquad W_{ij}^{n+1} = W_{ij}^n + \frac{kc}{2\Delta}(P_{ij+1}^n - P_{ij-1}^n - 2W_{ij}^n + W_{ij+1}^n + W_{ij-1}^n) + R W_{ij}^n $$

$$ (67) \qquad \begin{aligned} P_{ij}^{n+1} = P_{ij}^n + \frac{kc}{2\Delta}(V_{i+1,j}^n - V_{i-1,j}^n + W_{ij+1}^n - W_{ij-1}^n - 4P_{ij}^n + P_{ij+1}^n + P_{ij-1}^n + P_{i+1,j}^n + P_{i-1,j}^n) + R P_{ij}^n \end{aligned} $$

using Taylor formulas, with remainder terms $O(k^2) + O(\Delta^2)$

$$ (68) \qquad \begin{aligned} RV_{ij}^n = & - \int_{t^n}^{t^{n+1}} (t^{n+1} - t) \partial_{tt}^2 v(t, i\Delta, j\Delta) dt \\ & + \frac{kc}{2\Delta} \int_{i\Delta}^{(i+1)\Delta} ((i+1)\Delta - x) \partial_{xx}^2 p(kn, x, j\Delta) dx - \frac{kc}{2\Delta} \int_{i\Delta}^{(i-1)\Delta} ((i-1)\Delta - x) \partial_{xx}^2 p(kn, x, j\Delta) dx \\ & + \frac{kc}{2\Delta} \int_{i\Delta}^{(i+1)\Delta} ((i+1)\Delta - x) \partial_{xx}^2 v(kn, x, j\Delta) dx + \frac{kc}{2\Delta} \int_{i\Delta}^{(i-1)\Delta} ((i-1)\Delta - x) \partial_{xx}^2 v(kn, x, j\Delta) dx \end{aligned} $$
$$
(69) \quad RW_{ij}^n = - \int_{t^n}^{t^{n+1}} (t^{n+1} - t) \partial_{tt}^2 w(t, i\Delta, j\Delta) dt \\
\qquad + \frac{kc}{2\Delta} \int_{j\Delta}^{(j+1)\Delta} ((j+1)\Delta - y) \partial_{yy}^2 p(kn, i\Delta, y) dy - \frac{kc}{2\Delta} \int_{j\Delta}^{(j-1)\Delta} ((j-1)\Delta - y) \partial_{yy}^2 p(kn, i\Delta, y) dy \\
\qquad + \frac{kc}{2\Delta} \int_{j\Delta}^{(j+1)\Delta} ((j+1)\Delta - y) \partial_{yy}^2 v(kn, i\Delta, y) dy + \frac{kc}{2\Delta} \int_{j\Delta}^{(j-1)\Delta} ((j-1)\Delta - y) \partial_{yy}^2 v(kn, i\Delta, y) dy
$$

$$
(70) \quad RP_{ij}^n = - \int_{t^n}^{t^{n+1}} (t^{n+1} - t) \partial_{tt}^2 p(t, i\Delta, j\Delta) dt \\
\qquad + \frac{kc}{2\Delta} \int_{(i\Delta}^{(i+1)\Delta} ((i+1)\Delta - x) \partial_{xx}^2 v(kn, x, j\Delta) dx - \frac{kc}{2\Delta} \int_{(i\Delta}^{(i-1)\Delta} ((i-1)\Delta - x) \partial_{xx}^2 v(kn, x, j\Delta) dx \\
\qquad + \frac{kc}{2\Delta} \int_{(i\Delta}^{(i+1)\Delta} ((i+1)\Delta - x) \partial_{xx}^2 p(kn, x, j\Delta) dx + \frac{kc}{2\Delta} \int_{(i\Delta}^{(i-1)\Delta} ((i-1)\Delta - x) \partial_{xx}^2 p(kn, x, j\Delta) dx \\
\qquad + \frac{kc}{2\Delta} \int_{(j\Delta}^{(j+1)\Delta} ((j+1)\Delta - y) \partial_{yy}^2 w(kn, i\Delta, y) dy - \frac{kc}{2\Delta} \int_{(j\Delta}^{(j-1)\Delta} ((j-1)\Delta - y) \partial_{yy}^2 w(kn, i\Delta, y) dy \\
\qquad + \frac{kc}{2\Delta} \int_{(j\Delta}^{(j+1)\Delta} ((j+1)\Delta - y) \partial_{yy}^2 p(kn, i\Delta, y) dy + \frac{kc}{2\Delta} \int_{(j\Delta}^{(j-1)\Delta} ((j-1)\Delta - y) \partial_{yy}^2 p(kn, i\Delta, y) dy
$$

as long as $(v, w, p)$ are $C^2$ smooth solutions. Then, on a square domain $[0, 1]^2$ with periodic boundary conditions say (for the sake of simplicity : see our numerical Section 4 for an illustration when that situation applies), one can study convergence using the Fourier modes $\hat{Q}_{gh}^n$ of the error vector $Q_{ij}^n = \sum_{g,h \in \mathbb{N}^2} \hat{Q}_{gh}^n e^{2i\pi(ig+hj)} \in \mathbb{R}^3$ with components $V_{ij}^n = v_{ij}^n - v(kn, \Delta i, \Delta j)$, $W_{ij}^n = w_{ij}^n - w(kn, \Delta i, \Delta j)$, $P_{ij}^n = p_{ij}^n - p(kn, \Delta i, \Delta j)$ following Von Neumann. The Fourier modes satisfy

$$
\hat{Q}_{gh}^{n+1} = C(g, h) \hat{Q}_{gh}^{n} + \hat{R}_{gh}^{n}
$$

with $C(g, h)$ the amplification matrix of the linear scheme (65–66–68), and $C(g, h)$
is power-bounded as long as $k \in (\epsilon_0, \frac{2}{c})\Delta$ see e.g. [17, 11] : then the upwind scheme
is $\ell^2$ stable and convergence to smooth solutions holds at “first order” i.e.

$$
\sqrt{\sum_{ij} |Q_{ij}^N|^2 = \sum_{gh} |\hat{Q}_{ij}^N|^2 \leq C_N \Delta}
$$

with $C_N$ independent of $\Delta$, insofar as $R_{gh}^n = O(k^2 + \Delta^2)$ and $\sum_{ij} |Q_{ij}^0|^2 = O(\Delta^2)$.
To improve the convergence rate, one classically considers

$$
(71) \quad q_{ij}^{(1)} = q_{ij}^{n} + \frac{k}{\Delta} (\mathcal{F}_{i,j-\frac{1}{2}} + \mathcal{F}_{i-\frac{1}{2},j} - \mathcal{F}_{i,j+\frac{1}{2}} - \mathcal{F}_{i+\frac{1}{2},j}),
$$

$$
(72) \quad q_{ij}^{n+1} = \frac{q_{ij}^{n}}{2} + \frac{q_{ij}^{(1)}}{2} + \frac{k}{2\Delta} \left( \mathcal{F}_{i,j-\frac{1}{2}}^{(1)} + \mathcal{F}_{i-\frac{1}{2},j}^{(1)} - \mathcal{F}_{i,j+\frac{1}{2}}^{(1)} - \mathcal{F}_{i+\frac{1}{2},j}^{(1)} \right)
$$

based on a second-order time-integrator like RK2 method instead of (61), with
numerical fluxes $\mathcal{F}_{i,j \pm \frac{1}{2}}^{(1)}$, $\mathcal{F}_{i \pm \frac{1}{2},j}^{(1)}$, consistently improved, i.e. using at cell faces a pie-
cewise bilinear reconstruction of the fields $q$. Precisely, the numerical fluxes are
computed from the value of reconstructed variables $\tilde{q}$ evaluated at cell faces rather
than the cell value $q$. A linear reconstruction $\tilde{q}$ on a cell $C_{ij}$ has the form

$$ (73) \qquad \tilde{q}(x, y) = c_0 + c_1\varphi_1(x, y) + c_2\varphi_2(x, y) $$

with $c_0, c_1, c_2$ coefficients to be determined, using basis functions with zero mean :

$$ \varphi_1 = \frac{x - x_i}{\Delta}, \quad \varphi_2 = \frac{y - y_j}{\Delta}. $$

The coefficients $c_k$ are typically defined such that $\tilde{q}$ conserves exactly cell-averaged value $q_0 = q_{ij}^n$ on the central cell $C_0 = C_{ij}$, so $c_0 = q_0$, while it matches some averaged state for a given set $\{C_l\}_l$ of neighbour cells in a least-square sense :

$$ \tilde{q} = \arg \min \sum_{C_l} \left( q_l - \frac{1}{\Delta^2} \int_{C_l} \tilde{q} dxdy \right)^2 . $$

Consequently, the coefficients $c_1$ and $c_2$ are least-square solution of the (overdetermined) linear system

$$ \sum_{k=1}^{2} c_k \bar{\varphi}_{kl} = q_l - q_0 $$

with $\bar{\varphi}_{kl} = \frac{1}{\Delta^2} \int_{C_l} \varphi_k dxdy$, i.e. the averaged value of the basis function $\varphi_k$ over $C_l$. Introducing the matrix $\mathbf{A}$ whose coefficients are $A_{kl} = \bar{\varphi}_{kl}$ and denoting $\mathbf{B} = (\mathbf{A}^T \mathbf{A})^{-1} \mathbf{A}^T$ (easily pre-computed for the case of Cartesian mesh) one obtains

$$ (74) \qquad c_k = \sum_l B_{kl}(q_l - q_0). $$

In the D2Q5 case with neighbour cells $\{C_l\}_l = \{C_1, C_2, C_3, C_4\}$ we obtain

$$ \mathbf{A} = \begin{pmatrix} 1 & 0 \\ 0 & 1 \\ -1 & 0 \\ 0 & -1 \end{pmatrix}, \quad \mathbf{B} = \begin{pmatrix} \frac{1}{2} & 0 & -\frac{1}{2} & 0 \\ 0 & \frac{1}{2} & 0 & -\frac{1}{2} \end{pmatrix}, \\ c_1 = \frac{q_1 - q_3}{2}, \quad c_2 = \frac{q_2 - q_4}{2}. $$

The resulting linear reconstruction $\tilde{q}$ is none other than the classical centered discretization of gradient in each direction.

## RÉFÉRENCES

[1] Denise Aregba-Driollet and Roberto Natalini. Discrete kinetic schemes for multidimensional systems of conservation laws. *SIAM J. Numer. Anal.*, 37(6):1973–2004, May 2000.

[2] Hubert Baty, Florence Drui, Philippe Helluy, Emmanuel Franck, Christian Klingenberg, and Lukas Thanhäuser. A robust and efficient solver based on kinetic schemes for magnetohydrodynamics (MHD) equations. *Appl. Math. Comput.*, 440:Paper No. 127667, 18, 2023.

[3] Thomas Bellotti. Truncation errors and modified equations for the lattice Boltzmann method via the corresponding finite difference schemes. *ESAIM Math. Model. Numer. Anal.*, 57(3):1225–1255, 2023.

[4] Thomas Bellotti, Benjamin Graille, and Marc Massot. Finite difference formulation of any lattice Boltzmann scheme. *Numer. Math.*, 152(1):1–40, 2022.

[5] S. Benzoni-Gavage and D. Serre. *Multidimensional hyperbolic partial differential equations*. Oxford Mathematical Monographs. The Clarendon Press, Oxford University Press, Oxford, 2007. First order systems and applications.

[6] François Bouchut. Construction of BGK models with a family of kinetic entropies for a given system of conservation laws. *Journal of Statistical Physics*, 95:113–170, 1999. 10.1023/A:1004525427365.
[7] François Bouchut. Entropy satisfying flux vector splittings and kinetic BGK models. *Numerische Mathematik*, 94 :623–672, 2003. 10.1007/s00211-002-0426-9.

[8] François Bouchut. *Nonlinear stability of finite volume methods for hyperbolic conservation laws and well-balanced schemes for sources*. Frontiers in Mathematics. Birkhäuser Verlag, Basel, 2004.

[9] Daniele Corallo, Willy Dörfler, and Christian Wieners. Space-time discontinuous Galerkin methods for weak solutions of hyperbolic linear symmetric Friedrichs systems. *J. Sci. Comput.*, 94(1):Paper No. 27, 30, 2023.

[10] David Coulette, Clémentine Courtès, Emmanuel Franck, and Laurent Navoret. Vectorial kinetic relaxation model with central velocity. Application to implicit relaxations schemes. *Commun. Comput. Phys.*, 27(4):976–1013, 2020.

[11] Jean-François Coulombel. On the strong stability of finite difference schemes for hyperbolic systems in two space dimensions. *Calcolo*, 51(1):97–108, 2014.

[12] Constantine M. Dafermos. *Hyperbolic conservation laws in continuum physics*, volume 325 of *Grundlehren der Mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences]*. Springer-Verlag, Berlin, Berlin, 2000.

[13] Lokenath Debnath and Dambaru Bhatta. *Integral transforms and their applications*. CRC Press, Boca Raton, FL, third edition, 2015.

[14] François Dubois. Equivalent partial differential equations of a lattice Boltzmann scheme. *Comput. Math. Appl.*, 55(7):1441–1449, 2008.

[15] François Dubois. Stable lattice Boltzmann schemes with a dual entropy approach for monodimensional nonlinear waves. *Comput. Math. Appl.*, 65(2):142–159, 2013.

[16] François Dubois. Simulation of strong nonlinear waves with vectorial lattice boltzmann schemes. *International journal of modern physics. C, Computational physics, physical computation*, 25(12):1441014, 2014.

[17] Edwige Godlewski and Pierre-Arnaud Raviart. *Numerical approximation of hyperbolic systems of conservation laws*, volume 118 of *Applied Mathematical Sciences*. Springer-Verlag, New York, 1996.

[18] Edwige Godlewski and Pierre-Arnaud Raviart. *Numerical approximation of hyperbolic systems of conservation laws*, volume 118 of *Applied Mathematical Sciences*. Springer-Verlag, New York, [2021] ©2021. Second edition [of 1410987].

[19] Bertil Gustafsson. *High order difference methods for time dependent PDE*, volume 38 of *Springer Series in Computational Mathematics*. Springer-Verlag, Berlin, 2008.

[20] Romane Hélie. *Schéma de relaxation pour la simulation de plasmas dans les tokamaks*. Theses, Université de Strasbourg, March 2023.

[21] J. Hesthaven. *Numerical methods for conservation laws*. SIAM, 2018.

[22] Jonathan Jung and Vincent Perrier. Long time behavior of finite volume discretization of symmetrizable linear hyperbolic systems. *IMA J. Numer. Anal.*, 43(1):326–356, 2023.

[23] Randall J. LeVeque. *Finite volume methods for hyperbolic problems*. Cambridge Texts in Applied Mathematics. Cambridge University Press, Cambridge, 2002.

[24] Peter Monk and Endre Süli. A convergence analysis of Yee's scheme on nonuniform grids. *SIAM J. Numer. Anal.*, 31(2):393–412, 1994.

[25] Kane S. Yee. Numerical solution of initial boundary value problems involving maxwell's equations in isotropic media. *IEEE Transactions on Antennas and Propagation*, 14:302–307, 1966.

[26] Wen-An Yong, Weifeng Zhao, and Li-Shi Luo. Theory of the lattice Boltzmann method : derivation of macroscopic equations via the Maxwell iteration. *Phys. Rev. E*, 93(3):033310, 11, 2016.

UNIVERSITÉ PARIS 13, VILLETANEUSE, FRANCE

LHSV, ENPC, INSTITUT POLYTECHNIQUE DE PARIS, EDF R&D, CHATOU (SEBASTIEN.BOYAVAL@ENPC.FR)
& INRIA, FRANCE

LHSV, ENPC, INSTITUT POLYTECHNIQUE DE PARIS, EDF R&D, CHATOU, FRANCE