# Discontinuous-galerkin methods for a kinetic model of self-organized dynamics

Francis Filbet, Chi-Wang Shu

► To cite this version:

Francis Filbet, Chi-Wang Shu. Discontinuous-galerkin methods for a kinetic model of self-organized dynamics. Mathematical Models and Methods in Applied Sciences, In press, 28 (6), pp.1171-1197. 10.1142/S0218202518500318. hal-01526278

HAL Id: hal-01526278

https://hal.science/hal-01526278

Submitted on 22 May 2017

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.

Public Domain
# DISCONTINUOUS-GALERKIN METHODS FOR A KINETIC MODEL OF SELF-ORGANIZED DYNAMICS

FRANCIS FILBET AND CHI-WANG SHU

**ABSTRACT.** This paper deals with the numerical resolution of kinetic models for systems of self-propelled particles subject to alignment interaction and attraction-repulsion. We focus on the kinetic model considered in [18, 17] where alignment is taken into account in addition of an attraction-repulsion interaction potential. We apply a discontinuous Galerkin method for the free transport and non-local drift velocity together with a spectral method for the velocity variable. Then, we analyse consistency and stability of the semi-discrete scheme. We propose several numerical experiments which provide a solid validation of the method and its underlying concepts.

**Key words:** Self-propelled particles, alignment dynamics, kinetic model, discontinuous Galerkin method

**AMS Subject classification:** 35L60, 35K55, 35Q80, 82C05, 82C22, 82C70, 92D50.

## CONTENTS

<table><tr><td>1. Introduction</td><td>1</td></tr><tr><td>1.1. Agent-based model of self-alignment with attraction-repulsion</td><td>3</td></tr><tr><td>1.2. Kinetic model of self-alignment with attraction-repulsion</td><td>3</td></tr><tr><td>2. Numerical Methods</td><td>4</td></tr><tr><td>2.1. Notations</td><td>4</td></tr><tr><td>2.2. The semi-discrete discontinuous Galerkin method</td><td>5</td></tr><tr><td>2.3. Temporal discretizations</td><td>7</td></tr><tr><td>3. Conservation and stability</td><td>7</td></tr><tr><td>4. Proof of Theorem 2.2</td><td>10</td></tr><tr><td>4.1. Basic results</td><td>10</td></tr><tr><td>4.2. Error estimates ||f - f<sub>h</sub>||<sub>L<sup>2</sup></sub></td><td>15</td></tr><tr><td>5. Numerical simulations</td><td>16</td></tr><tr><td>5.1. Accuracy test</td><td>16</td></tr><tr><td>5.2. Taylor-Green vortex problem</td><td>17</td></tr><tr><td>5.3. Formation of bands problem</td><td>17</td></tr><tr><td>6. Conclusion and perspective</td><td>19</td></tr><tr><td>Acknowledgements</td><td>19</td></tr><tr><td>References</td><td>20</td></tr></table>

## 1. INTRODUCTION

Theoretical and mathematical biology communities have paid a great deal of attention to explain large scale structures in animal groups. Coherent structures appearing from seemingly direct interactions between individuals have been reported in many different species like fishes, birds, and insects [39, 44, 42, 5, 15, 6] and many others, see also the reviews [10, 8, 40]. There has been an intense literature about the modeling of interactions between individuals among
animal societies such as fish schools, bird flocks, herds of mammalians, etc. We refer for instance to [1, 2, 15, 37] but an exhaustive bibliography is out of reach. Among these models, the Vicsek model [52] has received particular attention due to its simplicity and the universality of its qualitative features. This model is an individual based model or agent-based model which consists of a time-discretized set of Ordinary Differential Equations for the particle positions and velocities. A time-continuous version of this model and its kinetic formulation are available in [18]. A rigorous derivation of this kinetic model from the time-continuous Vicsek model can be found in [4] and in [17] when adding an attraction-repulsion force.

On the other hand, hydrodynamic models are attractive over particle ones due to their computational efficiency. For this reason, many such models have been proposed in the literature [9, 12, 23, 42, 50, 51]. However, most of them are phenomenological. For instance in [18], the authors propose one of the first rigorous derivations of a hydrodynamic version of the Vicsek model (see also [41, 45, 46] for phenomenological derivations). It has been expanded in [19] to account for a model of fish behavior where particles interact through curvature control, and in [20] to include diffusive corrections. Other variants have also been investigated [16, 32, 33]. For instance, [32] studies the influence of a vision angle and of the dependency of the alignment frequency upon the local density, whereas in [16, 33], the authors study a modification of the model which results in phase transitions from disordered to ordered equilibria as the density increases and reaches a threshold, in a way similar to polymer models [24, 43].

In this paper, we will focus on the numerical approximation of a kinetic model for self-propelled particles. The self-propulsion speed is supposed to be constant and identical for all the particles. Therefore, the velocity variable reduces to its orientation in the $(d-1)$-dimensional sphere $S^{d-1}$. The particle interactions consist in two parts:

* an alignment rule which tends to relax the particle velocity to the local average orientation;

* an attraction-repulsion rule which makes the particles move closer or farther away from each other.

This model is inspired both by the Vicsek model [52] and the Couzin model [2, 15] describing interactions at the microscopic level. This approach has led to different types of models for swarming: microscopic models and macroscopic models involving macroscopic quantities (e.g. mass, flux). Here we study an intermediate approach, called the mesoscopic scale, where we investigate the time evolution of a distribution function of particles $f(t, \mathbf{x}, \mathbf{v})$, depending on time $t \ge 0$, position $\mathbf{x} \in \Omega \subset \mathbb{R}^d$ and velocity $\mathbf{v} \in S^{d-1}$. This distribution function is solution to a kinetic equation describing the motion of particles and their change of directions. Several works have already studied kinetic models for swarming [7, 38, 9], but few have done a numerical investigation. For macroscopic and microscopic models, we refer to [18, 21, 22] for discrete particle approximation of microscopic models and finite volume or finite difference approximations [22] where complex structures may be observed as in [27, 22]. Concerning kinetic model, we refer to [34] for the first work on this topic, where the authors propose a spectral discretization of the operator describing the change of direction with the flavor of what has been done for the Boltzmann equation [47, 29]. This spectral discretization is coupled with a finite volume approximation for the transport [25, 26, 28, 31]. For the study of bacteria's motion, we also mention [30] where similar structures have been observed (band's formation). Here we present a local discontinuous Galerkin method for computing the approximate solution to the kinetic model and to get high order approximations in time, space and velocity. Indeed, the preservation of high order accuracy allows to investigate complex structures in space as it has already been observed for macroscopic models [21, 22].
The paper is organized as follows: we first present precisely the kinetic model and give the main assumptions on the regularity of the unique solution to prove convergence and error estimates on the approximation to the exact smooth solution. Then, in Section 2 we develop a numerical scheme (local discontinuous Galerkin method) for the kinetic model. In Sections 3 and 4, we perform a stability and convergence analysis of the proposed numerical methods. Numerical investigations are presented in Section 5 where the order accuracy is verified and we observe the formation of complex structures.

1.1. Agent-based model of self-alignment with attraction-repulsion. The starting point of this study is an Individual-Based Model of particles interacting through self-alignment [52] and attraction-repulsion [2, 15]. Specifically, we consider N particles $x_i \in \mathbb{R}^d$, with $d = 2$ or 3, moving at a constant speed $v_i \in \mathbb{S}^{d-1}$. Each particle adjusts its velocity to align with its neighbors and to get closer or further away. Therefore, the evolution of each particle is modeled by the following dynamics: for any $1 \le i \le N$

$$ (1.1) \qquad \left\{ \begin{array}{l} \displaystyle \frac{d\mathbf{x}_i}{dt} = \mathbf{v}_i, \\[1em] \displaystyle d\mathbf{v}_i = \mathbf{P}_{\mathbf{v}_i^{\perp}} (\bar{\mathbf{v}}_i dt + \sqrt{2\nu} \, d\mathbf{B}_i^i), \end{array} \right. $$

where $\mathbf{B}_i^i$ is a Brownian motion and $d$ represents the noise intensity whereas $\mathbf{P}_{\mathbf{v}_i^{\perp}}$ is the projection matrix onto the normal plane to $\mathbf{v}_i$:

$$ \mathbf{P}_{\mathbf{v}^{\perp}} = \mathrm{Id} - \mathbf{v} \otimes \mathbf{v}, $$

which ensures that $\mathbf{v}_i$ stays of norm 1.

Both the alignment and attraction-repulsion rules are taken into account in the macroscopic velocity $\bar{\mathbf{v}}_i \in \mathbb{S}^{d-1}$:

$$ \bar{\mathbf{v}}_i = \frac{1}{|\mathbf{J}_i + \mathbf{R}_i|} (\mathbf{J}_i + \mathbf{R}_i), $$

where $\mathbf{J}_i$ counts for the alignment and $\mathbf{R}_i$ for the attraction-repulsion:

$$ (1.2) \qquad \mathbf{J}_i = \sum_{j=1}^{N} k(|\mathbf{x}_j - \mathbf{x}_i|) \mathbf{v}_j, \qquad \mathbf{R}_i = -\sum_{j=1}^{N} \nabla_{\mathbf{x}_i} \phi(|\mathbf{x}_j - \mathbf{x}_i|), $$

where the kernel $k$ is a positive function, $\phi'$ can be both negative (repulsion) and positive (attraction) and for simplicity, we will assume that both $k$ and $\phi$ are compactly supported in $[0, \infty)$.

1.2. Kinetic model of self-alignment with attraction-repulsion. When the number of particles becomes large, that is $N \to \infty$, one can formally derive a Vlasov type equation. It describes the evolution of a system of particles under the effects of external and self-consistent fields. The unknown $f(t, \mathbf{x}, \mathbf{v})$, depending on the time $t$, the position $\mathbf{x}$, and the velocity $\mathbf{v}$, represents the distribution of particles in phase space for each species with $(\mathbf{x}, \mathbf{v}) \in \Omega \times \mathbb{S}^{d-1}$, $d = 1, \dots, 3$, where $\Omega \subset \mathbb{R}^d$. Its behaviour is given by the Vlasov equation [4, 18, 49, 17],

$$ (1.3) \qquad \frac{\partial f}{\partial t} + \mathbf{v} \cdot \nabla_{\mathbf{x}} f = -\mathrm{div}_{\mathbf{v}} [\mathbf{P}_{\mathbf{v}^{\perp}} \mathbf{v}_f f - \nu \nabla_{\mathbf{v}} f], $$
where $\nu > 0$ and

$$ (1.4) \qquad \left\{ \begin{array}{l} \mathbf{v}_f = \frac{1}{|\mathbf{J}_f + \mathbf{R}_f|} (\mathbf{J}_f + \mathbf{R}_f), \\[1.5em] \mathbf{J}_f = \int_{\Omega \times \mathbb{S}^{d-1}} k(|\mathbf{x} - \mathbf{x}'|) \mathbf{v}' f(t, \mathbf{x}', \mathbf{v}') d\mathbf{x}' d\mathbf{v}', \\[1.5em] \mathbf{R}_f = -\nabla_{\mathbf{x}} \int_{\Omega \times \mathbb{S}^{d-1}} \phi(|\mathbf{x} - \mathbf{x}'|) f(t, \mathbf{x}', \mathbf{v}') d\mathbf{x}' d\mathbf{v}'. \end{array} \right. $$

In general, the function $\phi$ is such that $\phi(r) \to 0$ when $r \to \infty$, but here we will assume that both $k$ and $\phi$ are nonnegative functions which satisfy

$$ (1.5) \qquad k, \phi \in C_c^p([0, \infty)), \quad \text{with } p \ge 2 $$

and for periodic boundary conditions in space, we have

$$ \left\{ \begin{array}{l} \mathbf{J}_f(t, \mathbf{x}) = \int_{\text{supp}(k)} k(|\mathbf{y}|) \rho \mathbf{u}(t, \mathbf{x} + \mathbf{y}) \, d\mathbf{y}, \\ \\ \mathbf{R}_f(t, \mathbf{x}) = -\nabla_{\mathbf{x}} \int_{\text{supp}(\phi)} \phi(|\mathbf{y}|) \rho(t, \mathbf{x} + \mathbf{y}) \, d\mathbf{y}. \end{array} \right. $$

Furthermore we assume that the system (1.3)-(1.4) has a smooth solution such that

$$ f \in H^{k+2}([0, T] \times \Omega \times \mathbb{S}^{d-1}), $$

with $\mathbf{J}_f$ and $\mathbf{R}_f$ such that for any $T > 0$, there exists a constant $\xi_T > 0$ such that for all $(t, \mathbf{x}) \in [0, T] \times \Omega$

$$ (1.6) \qquad |\mathbf{J}_f(t, \mathbf{x}) + \mathbf{R}_f(t, \mathbf{x})| \ge \xi_T. $$

Using the distribution $f$ and a rescaling of (1.3), it is possible to identify the asymptotic behavior of the model in different regimes as in [17] and to recover some classical hydrodynamic model for the self-organized dynamics.

The model studied in this paper is a generalization of the model of [18] with the addition of an attraction-repulsion interaction potential [17].

## 2. NUMERICAL METHODS

In this section, we will introduce the discontinuous Galerkin algorithm for the system (1.3)-(1.4). Discontinuous Galerkin methods are particularly suited for transport type equations with several attractive properties, such as their easiness for adaptivity and parallel computation, and their nice stability properties. We refer to the survey paper [14] and the references therein for an introduction to discontinuous Galerkin methods. For discontinuous Galerkin methods solving kinetic type equations we refer to [11, 3]. We consider an open set $\Omega \subset \mathbb{R}^d$ where all boundary conditions are periodic, and $f(t, \mathbf{x}, \mathbf{v})$ is assumed to be in the unit sphere $\mathbb{S}^{d-1}$.

### 2.1. Notations.
Let $\mathcal{T}_h^x = \{K_x\}$ and $\mathcal{T}_h^v = \{K_v\}$ be partitions of $\Omega$ and $\mathbb{S}^{d-1}$, respectively, with $K_x$ and $K_v$ being Cartesian elements; then

$$ \mathcal{T}_h = \{K \subset \Omega \times \mathbb{S}^{d-1}; \ K = K_x \times K_v \in \mathcal{T}_h^x \times \mathcal{T}_h^v \} $$

defines a partition of $\Omega \times \mathbb{S}^{d-1}$.
Let $\mathcal{E}$ be the set of the edges of $\mathcal{T}_h$ will be $\mathcal{E} = \mathcal{E}_x \cup \mathcal{E}_v$ as

$$ 
\left\{ 
\begin{array}{l}
\mathcal{E}_x = \{ \sigma_x \times K_v : \sigma_x \in \partial K_x, K_v \in \mathcal{T}_h^v \}, \\
\mathcal{E}_v = \{ K_x \times \sigma_v : K_x \in \mathcal{T}_h^x, \sigma_v \in \partial K_v \}.
\end{array}
\right.
$$

Next we define the discrete spaces

$$ (2.1) \qquad \mathcal{G}_h^k = \{ g \in L^2(\Omega \times \mathbb{S}^{d-1}) : g|_K \in P^k(K), K \in \mathcal{T}_h \}, $$

and

$$ \mathcal{U}_h^k = \{ \mathbf{U} \in [L^2(\Omega \times \mathbb{S}^{d-1})]^{d-1} : \mathbf{U}|_K \in [P^k(K)]^{d-1}, K \in \mathcal{T}_h \}, $$

where $P^k(K)$ denotes the set of polynomials of total degree at most $k$ on $K$, and $k$ is a nonnegative integer.

Note the space $\mathcal{G}_h^k$, which we use to approximate $f$, is called P-type, and it can be replaced by the tensor product of P-type spaces in $\mathbf{x}$ and $\mathbf{v}$,

$$ \{ g \in L^2(\Omega \times \mathbb{S}^{d-1}) : g|_K \in P^k(K_{\mathbf{x}}) \times P^k(K_{\mathbf{v}}), K = K_{\mathbf{x}} \times K_{\mathbf{v}} \in \mathcal{T}_h \}, $$

or by the tensor product space in each variable, which is called Q-type

$$ \{ g \in L^2(\Omega \times \mathbb{S}^{d-1}) : g|_K \in Q^k(K_{\mathbf{x}}) \times Q^k(K_{\mathbf{v}}), K = K_{\mathbf{x}} \times K_{\mathbf{v}} \in \mathcal{T}_h \}. $$

Here $Q^k(K)$ denotes the set of polynomials of degree at most $k$ in each variable on $K$. The numerical methods formulated in this paper, as well as the conservation, stability, and error estimates, hold when any of the spaces above is used to approximate $f$.

**Remark 2.1.** In our simulations of Section 5, we use the P-type of (2.1) as it is the smallest and therefore renders the most cost efficient algorithm.

For piecewise functions defined with respect to $\mathcal{T}_h^x$ or $\mathcal{T}_h^v$, we further introduce the jumps and averages as follows. For $\alpha \in \{\mathbf{x}, \mathbf{v}\}$ and for any edge $\sigma = \{K_\alpha^+ \cap K_\alpha^-\} \in \mathcal{E}_\alpha$, with $\mathbf{n}_\alpha^\pm$ as the outward unit normal to $\partial K_\alpha^\pm$, $g^\pm = g|_{K_\alpha^\pm}$, the jumps across $\sigma$ and the averages are defined as

$$ (2.2) \qquad [g]_\alpha = g^+ - g^-, \quad \{g\}_\alpha = \frac{1}{2}(g^+ + g^-), \quad \alpha \in \{\mathbf{x}, \mathbf{v}\}. $$

2.2. The semi-discrete discontinuous Galerkin method. The numerical methods proposed in this section are formulated for the system (1.3)-(1.4). Given $k, r \ge 0$, the semi-discrete discontinuous Galerkin methods for the system (1.3)-(1.4) are defined by the following procedure: for any $K = K_{\mathbf{x}} \times K_{\mathbf{v}} \in \mathcal{T}_h$, we look for $(f_h, q_h) \in \mathcal{G}_h^k \times \mathcal{U}_h^k$, $v_{fh} \in \mathcal{U}_h^r$, such that for all $g \in \mathcal{G}_h^k$,

$$ (2.3) \qquad \int_K \frac{\partial f_h}{\partial t} g \,dx\,dv - \int_K f_h \mathbf{v} \cdot \nabla_{\mathbf{x}} g \,dx\,dv - \int_K (\mathbf{P}_{\mathbf{v}\perp} \mathbf{v}_{fh} - \nu \mathbf{q}_h) \cdot \nabla_{\mathbf{v}} g \,dx\,dv \\ + \int_{\sigma_{\mathbf{x}}} \widehat{f_h \mathbf{v}} \cdot \mathbf{n}_{\mathbf{x}} g^{-} \,ds_{\mathbf{x}}\,d\mathbf{v} + \int_{\sigma_{\mathbf{v}}} \widehat{(f_h \mathbf{P}_{\mathbf{v}\perp} \mathbf{v}_{fh} - \nu \mathbf{q}_h)} \cdot \mathbf{n}_{\mathbf{v}} g^{-} \,ds_{\mathbf{v}}\,d\mathbf{x} = 0, $$

where $\mathbf{n}_{\mathbf{x}}$ and $\mathbf{n}_{\mathbf{v}}$ are outward unit normals of $\partial K_{\mathbf{x}}$ and $\partial K_{\mathbf{v}}$, respectively, whereas $\mathbf{q}_h$ is given by

$$ (2.4) \qquad \int_K \mathbf{q}_h \cdot \mathbf{u} \,dx\,dv + \int_K f_h \operatorname{div}_{\mathbf{v}} \mathbf{u} \,dx\,dv - \int_{\sigma_{\mathbf{v}}} \widehat{f_h \mathbf{n}_{\mathbf{v}}} \cdot \mathbf{u}^{-} \,dx\,ds_{\mathbf{v}} = 0, \forall \mathbf{u} \in \mathcal{U}_h^k. $$

Furthermore, the velocity $\mathbf{v}_{fh} \in L^\infty(\Omega)$ with $\|\mathbf{v}_{fh}\| = 1$, and

$$ (2.5) \qquad \mathbf{v}_{fh}(t, \mathbf{x}) = \frac{1}{\|\mathbf{J}_h(t, \mathbf{x}) + \mathbf{R}_h(t, \mathbf{x})\|} (\mathbf{J}_h(t, \mathbf{x}) + \mathbf{R}_h(t, \mathbf{x})), $$
with $\mathbf{J}_h$ and $\mathbf{R}_h$ computed by

$$ (2.6) \qquad \left\{ \begin{array}{l} \mathbf{J}_h(t, \mathbf{x}) = \int_{\Omega} k(|\mathbf{x} - \mathbf{x}'|) \rho_h \mathbf{u}_h(t, \mathbf{x}') \, d\mathbf{x}', \\ \\ \mathbf{R}_h = \int_{\Omega} \nabla_{\mathbf{x}} \phi(|\mathbf{x} - \mathbf{x}'|) \rho_h(t, \mathbf{x}') \, d\mathbf{x}', \end{array} \right. $$

where $\rho_h$ and $\mathbf{u}_h$ are defined from the distribution function $f_h$, by

$$ \rho_h = \int_{\mathbb{S}^{d-1}} f_h d\mathbf{v}, \quad \rho_h \mathbf{u}_h = \int_{\mathbb{S}^{d-1}} \mathbf{v} f_h d\mathbf{v}. $$

All hat functions are numerical fluxes that are determined by upwinding for convection and local DG alternating for diffusion, i.e. for the convective terms in (2.3)

$$ (2.7) \qquad \left\{ \begin{array}{l} \widehat{f_h \mathbf{v}} \cdot \mathbf{n_x} = \mathbf{v} \cdot \mathbf{n_x} \{f_h\}_x - \frac{|\mathbf{v} \cdot \mathbf{n_x}|}{2} [f_h]_x, \\ \\ f_h \widehat{\mathbf{P}_{v\perp}} \mathbf{v}_{f_h} \cdot \mathbf{n_v} = \mathbf{P}_{v\perp} \mathbf{v}_{f_h} \cdot \mathbf{n_v} \{f_h\}_v - \frac{|\mathbf{P}_{v\perp} \mathbf{v}_{f_h}|}{2} [f_h]_v, \end{array} \right. $$

and for the diffusive terms in (2.4), we apply

$$ (2.8) \qquad \widehat{\mathbf{q}_h} \cdot \mathbf{n_v} = \{\mathbf{q}_h\} \cdot \mathbf{n_v} + \frac{C_{11}}{2} [f_h]_v, \qquad \widehat{f_h} \mathbf{n_v} = \{f_h\} \mathbf{n_v} + \frac{C_{22}}{2} [\mathbf{q}_h]_v $$

where $C_{11}, C_{22} > 0$.

This completes the definition of our Discontinuous-Galerkin method, but to facilitate its study, we recast its formulation. We sum (2.3) and (2.4) over all elements and define $a_h(\cdot)$ and $b_h(\cdot)$ such that for $(f_h, \mathbf{q}_h) \in \mathcal{G}_h^k \times \mathcal{U}_h^k$ and $g \in \mathcal{G}_h^k$

$$ (2.9) \qquad a_h(f_h, \mathbf{q}_h, \mathbf{v}_{f_h}, g) := \int_{\Omega \times \mathbb{S}^{d-1}} \left( \frac{\partial f_h}{\partial t} g - f_h \mathbf{v} \cdot \nabla_{\mathbf{x}} g \right) d\mathbf{x} d\mathbf{v} \\ - \int_{\Omega \times \mathbb{S}^{d-1}} (\mathbf{P}_{v\perp} \mathbf{v}_{f_h} f_h - \nu \mathbf{q}_h) \cdot \nabla_{\mathbf{v}} g \, d\mathbf{x} d\mathbf{v} \\ - \sum_{\sigma_x \in \mathcal{E}_x} \int_{\sigma_x} \widehat{f_h \mathbf{v}} \cdot \mathbf{n}_x [g]_x \, ds_x \, d\mathbf{v} \\ - \sum_{\sigma_v \in \mathcal{E}_v} \int_{\sigma_v} \widehat{(f_h \mathbf{P}_{v\perp} \mathbf{v}_{f_h} - \nu \widehat{\mathbf{q}_h})} \cdot \mathbf{n}_v [g]_v \, ds_v \, d\mathbf{x}, $$

and for $\mathbf{u} \in \mathcal{U}_h^k$

$$ (2.10) \qquad b_h(f_h, \mathbf{q}_h, u) := \int_{\Omega \times \mathbb{S}^{d-1}} (\mathbf{q}_h \cdot \mathbf{u} + f_h \operatorname{div}_v \mathbf{u}) \, d\mathbf{x} d\mathbf{v} \\ \qquad \qquad + \sum_{\sigma_v \in \mathcal{E}_v} \int_{\sigma_v} \widehat{f_h} \mathbf{n}_v \cdot [\mathbf{u}]_v \, d\mathbf{x} \, ds_v, $$

where we notice that $a_h$ (resp. $b_h$) is linear with respect to $(f_h, \mathbf{q}_h)$ and $g$ (resp. $\mathbf{u}$).

We prove the following boundedness and error estimate results.

**Theorem 2.2.** Assume that the solution to (1.3)-(1.4) is such that $f \in H^{k+2}([0, T] \times \Omega \times \mathbb{S}^{d-1})$ with the two hypothesis (1.5) and (1.6). We also assume that the initial data $f_h(0)$ is uniformly bounded in $L^2(\Omega \times \mathbb{S}^{d-1})$ and for $k \ge 0$, we consider the numerical solution $(f_h, \mathbf{q}_h) \in \mathcal{G}_h^k \times \mathcal{U}_h^k$
given by (2.3)-(2.8) supplemented with periodic boundary conditions. Then, for any $h_0 > 0$, there exists $C_T > 0$, depending on $f$, $T$ and $h_0$, such that for $h < h_0$

$$ \|f_h(t)\|_{L^2}^2 + \int_0^t \|\mathbf{q}_h(s)\|_{L^2}^2 ds \le C_T, \quad t \in [0, T] $$

and

$$ \|f(t) - f_h(t)\|_{L^2}^2 + \left( \int_0^t \|\mathbf{q}(s) - \mathbf{q}_h(s)\|_{L^2}^2 ds \right)^{1/2} \le C_T h^{k+1/2}, \quad t \in [0, T]. $$

**2.3. Temporal discretizations.** We use total variation diminishing (TVD) high-order Runge-Kutta methods to solve the method of lines ordinary differential equation resulting from the semi-discrete discontinuous Galerkin scheme,

$$ \frac{df_h}{dt} = \mathcal{R}(f_h). $$

Such time stepping methods are convex combinations of the Euler forward time discretization. The commonly used third-order TVD Runge-Kutta method is given by

$$ (2.11) \qquad \begin{cases} f_h^{(1)} = f_h^n + \Delta t \mathcal{R}(f_h^n), \\ f_h^{(2)} = \frac{1}{4} \left( 3 f_h^n + f_h^{(1)} + \Delta t \mathcal{R}(f_h^{(1)}) \right), \end{cases} $$

and

$$ (2.12) \qquad f_h^{n+1} = \frac{1}{3} \left( f_h^n + 2 f_h^{(2)} + 2 \Delta t \mathcal{R}(f_h^{(2)}) \right), $$

where $f_h^n$ represents a numerical approximation of the solution at discrete time $t_n$.

A detailed description of the TVD Runge-Kutta method can be found in [48]; see also [35] and [36] for strong-stability-preserving methods.

### 3. CONSERVATION AND STABILITY

In this section, we will establish conservation and stability properties of the semi-discrete discontinuous Galerkin methods. In particular, we prove that for periodic boundary condition, the total density (mass) is always conserved. We also show that $f_h$ is $L^2$ stable, which facilitates the error analysis of Section 4.

**Lemma 3.1 (Mass conservation).** Consider the numerical solution $(f_h, \mathbf{q}_h) \in \mathcal{G}_h^k \times \mathcal{U}_h^k$ for $k \ge 0$ given by (2.3)-(2.8) supplemented with periodic boundary conditions. Then it satisfies

$$ (3.1) \qquad \frac{d}{dt} \int_{\Omega \times \mathbb{S}^{d-1}} f_h d\mathbf{x} d\mathbf{v} = 0. $$

Equivalently, for $\rho_h(\mathbf{x}, t)$, for any $t > 0$, the following holds:

$$ (3.2) \qquad \int_{\Omega} \rho_h(t, \mathbf{x}) d\mathbf{x} = \int_{\Omega} \rho_h(0, \mathbf{x}) d\mathbf{x}. $$

*Proof.* Let $g(\mathbf{x}, \mathbf{v}) = 1$ and noticing that $g \in \mathcal{G}_h^k$, for any $k \ge 0$, is continuous $\nabla_{\mathbf{x}}g = 0$ and $\nabla_{\mathbf{v}}g = 0$. Taking this $g$ as the test function in (2.3), one has

$$ \frac{d}{dt} \int_K f_h d\mathbf{x} d\mathbf{v} + \int_{\sigma_{\mathbf{x}}} \widehat{f_h \mathbf{v}} \cdot \mathbf{n}_x ds_{\mathbf{x}} d\mathbf{v} + \int_{\sigma_{\mathbf{v}}} \left( f_h \widehat{\mathbf{P}_{\mathbf{v} \perp} \mathbf{v}}_{f_h} - \nu \widehat{\mathbf{q}_h} \right) \cdot \mathbf{n}_{\mathbf{v}} ds_{\mathbf{v}} d\mathbf{x} = 0. $$
Then summing over $K \in \mathcal{T}_h$ and thanks to the periodic boundary conditions, we get the conservation of total mass for any $t \ge 0$,

$$ \frac{d}{dt} \int_{\Omega \times \mathbb{S}^{d-1}} f_h(t) d\mathbf{x} d\mathbf{v} = 0. $$

Finally from the definition of $\rho_h$ and integrating in time from 0 to $t$, it gives (3.2). $\square$

Finally, we can obtain the $L^2$-stability result for $f_h$. This result will be used in the error analysis of Section 4.

**Lemma 3.2** ($L^2$-stability of $f_h$). Assume that the initial data $f_h(0)$ is uniformly bounded in $L^2(\Omega \times \mathbb{S}^{d-1})$ and for $k \ge 0$, consider the numerical solution $(f_h, \mathbf{q}_h) \in \mathcal{G}_h^k \times \mathcal{U}_h^k$ given by (2.3)-(2.8) supplemented with periodic boundary conditions. Then $(f_h, \mathbf{q}_h)$ satisfies for any $t \ge 0$

$$ \begin{aligned} & \frac{1}{2} \frac{d}{dt} \int_{\Omega \times \mathbb{S}^{d-1}} |f_h|^2 d\mathbf{x} d\mathbf{v} \\ & \quad + \nu \int_{\Omega \times \mathbb{S}^{d-1}} |\mathbf{q}_h|^2 d\mathbf{x} d\mathbf{v} \\ & \quad + \frac{1}{2} \sum_{\sigma_{\mathbf{v}} \in \mathcal{E}_{\mathbf{v}}} \int_{\sigma_{\mathbf{v}}} (|\mathbf{P}_{\mathbf{v} \perp} \mathbf{v}_{f_h} \cdot \mathbf{n}_{\mathbf{v}}| + \nu C_{11}) [f_h]_{\mathbf{v}}^2 + \nu C_{22} [\mathbf{q}_h]_{\mathbf{v}}^2 ds_{\mathbf{v}} d\mathbf{x} \\ & \quad + \frac{1}{2} \sum_{\sigma_{\mathbf{x}} \in \mathcal{E}_{\mathbf{x}}} \int_{\sigma_{\mathbf{x}}} |\mathbf{v} \cdot \mathbf{n}_{\mathbf{x}}| [f_h]_{\mathbf{x}}^2 ds_{\mathbf{x}} d\mathbf{v} \\ & = -\frac{1}{2} \int_{\Omega \times \mathbb{S}^{d-1}} f_h^2 \operatorname{div}_{\mathbf{v}} (\mathbf{P}_{\mathbf{v} \perp} \mathbf{v}_{f_h}) d\mathbf{x} d\mathbf{v}. \end{aligned} $$

*Proof.* Observing that

$$ \begin{aligned} \int_{K_{\mathbf{v}}} f_h \mathbf{P}_{\mathbf{v} \perp} \mathbf{v}_{f_h} \cdot \nabla_{\mathbf{v}} f_h \, d\mathbf{v} &= \frac{1}{2} \int_{K_{\mathbf{v}}} \mathbf{P}_{\mathbf{v} \perp} \mathbf{v}_{f_h} \cdot \nabla_{\mathbf{v}} f_h^2 \, d\mathbf{v}, \\ &= \frac{1}{2} \int_{K_{\mathbf{v}}} \operatorname{div}_{\mathbf{v}} (\mathbf{P}_{\mathbf{v} \perp} \mathbf{v}_{f_h} f_h^2) - f_h^2 \operatorname{div}_{\mathbf{v}} (\mathbf{P}_{\mathbf{v} \perp} \mathbf{v}_{f_h}) \, d\mathbf{v}, \\ &= \frac{1}{2} \int_{\partial K_{\mathbf{v}}} \mathbf{P}_{\mathbf{v} \perp} \mathbf{v}_{f_h} \cdot \mathbf{n}_{\mathbf{v}} |f_h|^{-2} ds_{\mathbf{v}} - \frac{1}{2} \int_{K_{\mathbf{v}}} f_h^2 \operatorname{div}_{\mathbf{v}} (\mathbf{P}_{\mathbf{v} \perp} \mathbf{v}_{f_h}) \, d\mathbf{v}, \end{aligned} $$

and summing over all control volume $K \in \mathcal{T}^h$ and recasting the edges, it yields

$$ \begin{aligned} \int_{\Omega \times \mathbb{S}^{d-1}} f_h \mathbf{P}_{\mathbf{v} \perp} \mathbf{v}_{f_h} \cdot \nabla_{\mathbf{v}} f_h \, d\mathbf{v} d\mathbf{x} &= -\frac{1}{2} \int_{\Omega \times \mathbb{S}^{d-1}} f_h^2 \operatorname{div}_{\mathbf{v}} (\mathbf{P}_{\mathbf{v} \perp} \mathbf{v}_{f_h}) \, d\mathbf{v} d\mathbf{x} \\ &\quad - \sum_{\sigma_{\mathbf{v}} \in \mathcal{E}_{\mathbf{v}}} \int_{\sigma_{\mathbf{v}}} \mathbf{P}_{\mathbf{v} \perp} \mathbf{v}_{f_h} \cdot \mathbf{n}_{\mathbf{v}} \{f_h\}_{\mathbf{v}} [f_h]_{\mathbf{v}} ds_{\mathbf{v}} d\mathbf{x}. \end{aligned} $$

Moreover, we have after recasting

$$ \int_{\Omega \times \mathbb{S}^{d-1}} f_h \mathbf{v} \cdot \nabla_{\mathbf{x}} f_h \, d\mathbf{v} d\mathbf{x} = - \sum_{\sigma_{\mathbf{x}} \in \mathcal{E}_{\mathbf{x}}} \int_{\sigma_{\mathbf{x}}} \mathbf{v} \cdot \mathbf{n}_{\mathbf{x}} \{f_h\}_{\mathbf{x}} [f_h]_{\mathbf{x}} ds_{\mathbf{v}} d\mathbf{x}. $$

Then we take $g = f_h$ in (2.9), it gives after an integration by part on each control volume

$$ (3.3) \quad \frac{1}{2} \frac{d}{dt} \int_{\Omega \times \mathbb{S}^{d-1}} |f_h|^2 d\mathbf{x} d\mathbf{v} + \frac{1}{2} \int_{\Omega \times \mathbb{S}^{d-1}} f_h^2 \operatorname{div}_{\mathbf{v}} (\mathbf{P}_{\mathbf{v} \perp} \mathbf{v}_{f_h}) d\mathbf{x} d\mathbf{v} + I_1 + I_2 + I_3 = 0, $$
with

$$
\left\{
\begin{aligned}
I_1 &= - \sum_{\sigma_x \in \mathcal{E}_x} \int_{\sigma_x} (\widehat{f_h \mathbf{v}} \cdot \mathbf{n}_x - \mathbf{v} \cdot \mathbf{n}_x \{f_h\}_x) [f_h]_x \, ds_x \, d\mathbf{v}, \\
I_2 &= - \sum_{\sigma_v \in \mathcal{E}_v} \int_{\sigma_v} (f_h \widehat{\mathbf{P}_{v\perp} \mathbf{v}_{f_h}} \cdot \mathbf{n}_v - \mathbf{P}_{v\perp} \mathbf{v}_{f_h} \cdot \mathbf{n}_v \{f_h\}_v) [f_h]_v \, ds_v \, d\mathbf{x}, \\
I_3 &= \nu \left( \int_{\Omega \times \mathbb{S}^{d-1}} \mathbf{q}_h \cdot \nabla_{\mathbf{v}} f_h \, d\mathbf{v} \, d\mathbf{x} + \sum_{\sigma_v \in \mathcal{E}_v} \int_{\sigma_v} \widehat{\mathbf{q}_h} \cdot \mathbf{n}_v [f_h]_v \, ds_v \, d\mathbf{x} \right).
\end{aligned}
\right.
$$

Let us prove that each term $I_k$, for $1 \le k \le 3$, is nonnegative. On the one hand using the definition of the upwinding flux (2.7), we simply have

$$
(3.4) \qquad I_1 = \frac{1}{2} \sum_{\sigma_x \in \mathcal{E}_x} \int_{\sigma_x} |\mathbf{v} \cdot \mathbf{n}_x| [f_h]^2_x ds_x d\mathbf{v} \ge 0
$$

and

$$
(3.5) \qquad I_2 = \frac{1}{2} \sum_{\sigma_v \in \mathcal{E}_v} \int_{\sigma_x} |\mathbf{P}_{v\perp} \mathbf{v}_{f_h} \cdot \mathbf{n}_v| [f_h]^2_v ds_v d\mathbf{x} \ge 0.
$$

On the other hand, to deal with the last term $I_3$, we choose $\mathbf{u} = \mathbf{q}_h$ in (2.10), hence we get

$$
\nu \left( \int_{\Omega \times \mathbb{S}^{d-1}} (|\mathbf{q}_h|^2 + f_h \operatorname{div}_v \mathbf{q}_h) \, d\mathbf{x} \, d\mathbf{v} + \sum_{\sigma_v \in \mathcal{E}_v} \int_{\sigma_v} \hat{f}_h \mathbf{n}_v \cdot [\mathbf{q}_h]_v \, d\mathbf{x} \, ds_v \right) = 0.
$$

Then performing an integration by part in velocity of the second term and using the definition
of $I_3$, we have

$$
\begin{align*}
I_3 ={}& \nu \int_{\Omega \times \mathbb{S}^{d-1}} |\mathbf{q}_h|^2 \, d\mathbf{x} \, d\mathbf{v} \\
& + \nu \sum_{\sigma_v \in \mathcal{E}_v} \int_{\sigma_v} (f_h^- \mathbf{q}_h^- - f_h^+ \mathbf{q}_h^+) \cdot \mathbf{n}_v + \widehat{\mathbf{q}}_h \cdot \mathbf{n}_v [f_h]_v + \hat{f}_h \mathbf{n}_v \cdot [\mathbf{q}_h]_v \, d\mathbf{x} \, ds_v.
\end{align*}
$$

Therefore from the definition of the “alternating fluxes” (2.8), we finally get that

$$
(3.6) \quad I_3 = \nu \left( \int_{\Omega \times \mathbb{S}^{d-1}} |\mathbf{q}_h|^2 \, d\mathbf{x} \, d\mathbf{v} + \frac{1}{2} \sum_{\sigma_v \in \mathcal{E}_v} \int_{\sigma_v} C_{11} [f_h]^2_{\mathbf{v}} + C_{22} [\mathbf{q}_h]^2_{\mathbf{v}} \right) \ge 0.
$$

Gathering (3.3) together with (3.4)-(3.6), we obtain the result

$$
\begin{align*}
& \frac{1}{2} \frac{d}{dt} \int_{\Omega \times \mathbb{S}^{d-1}} |f_h|^2 \, d\mathbf{x} \, d\mathbf{v} + \nu \int_{\Omega \times \mathbb{S}^{d-1}} |\mathbf{q}_h|^2 \, d\mathbf{x} \, d\mathbf{v} \\
& + \frac{1}{2} \sum_{\sigma_v \in \mathcal{E}_v} \int_{\sigma_v} (|\mathbf{P}_{v\perp} \mathbf{v}_{f_h} \cdot \mathbf{n}_v| + \nu C_{11}) [f_h]^2_{v} + \nu C_{22} [\mathbf{q}_h]^2_{v} \, ds_v \, d\mathbf{x} \\
& + \frac{1}{2} \sum_{\sigma_x \in \mathcal{E}_x} \int_{\sigma_x} |\mathbf{v} \cdot \mathbf{n}_x| [f_h]^2_x \, ds_x \, d\mathbf{v} = -\frac{1}{2} \int_{\Omega \times \mathbb{S}^{d-1}} f_h^2 \operatorname{div}_v (\mathbf{P}_{v\perp} \mathbf{v}_{f_h}) \, d\mathbf{x} \, d\mathbf{v}.
\end{align*}
$$

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

From this Lemma, we get $L^2$ boundedness estimates on $(f_h, \mathbf{q}_h)$ and on the macroscopic
quantities.
**Proposition 3.3.** Under the assumptions of Lemma 3.2, consider the numerical solution $(f_h, \mathbf{q}_h) \in \mathcal{G}_h^k \times \mathcal{U}_h^k$ given by (2.3)-(2.8) supplemented with periodic boundary conditions. Then for any $t \ge 0$

$$ (3.7) \qquad \|f_h(t)\|_{L^2}^2 + 2\nu e^t \int_0^t \|\mathbf{q}_h(s)\|_{L^2}^2 ds \le \|f_h(0)\|_{L^2}^2 e^t. $$

Furthermore $(\rho_h, \rho_h \mathbf{u}_h)$ computed from the distribution function $f_h$ satisfies for any $t \ge 0$

$$ (3.8) \qquad \left\{ \begin{array}{l} \|\rho_h(t)\|_{L^2(\Omega)} \le \text{vol}(\mathbb{S}^{d-1})^{1/2} \|f_h(0)\|_{L^2} e^t, \\ \|\rho_h \mathbf{u}_h(t)\|_{L^2(\Omega)} \le \text{vol}(\mathbb{S}^{d-1})^{1/2} \|f_h(0)\|_{L^2} e^t \end{array} \right. $$

and $\|\mathbf{u}_h\|_{L^\infty(\Omega)} \le 1$.

*Proof.* Starting from Lemma 3.2 and using the fact that $\|\mathbf{v}_{f_h}\| \le 1$, we have

$$ \frac{d}{dt} \int_{\Omega \times \mathbb{S}^{d-1}} |f_h|^2 d\mathbf{x}d\mathbf{v} + 2\nu \int_{\Omega \times \mathbb{S}^{d-1}} |\mathbf{q}_h|^2 d\mathbf{x}d\mathbf{v} \le \int_{\Omega \times \mathbb{S}^{d-1}} f_h^2 d\mathbf{x}d\mathbf{v}, $$

or

$$ \frac{d}{dt} \left( e^{-t} \int_{\Omega \times \mathbb{S}^{d-1}} |f_h|^2 d\mathbf{x}d\mathbf{v} \right) + 2\nu \int_{\Omega \times \mathbb{S}^{d-1}} |\mathbf{q}_h|^2 d\mathbf{x}d\mathbf{v} \le 0. $$

Hence after integration, it yields to the $L^2$ estimate for any $t \ge 0$,

$$ \int_{\Omega \times \mathbb{S}^{d-1}} |f_h(t)|^2 d\mathbf{x}d\mathbf{v} + 2\nu e^t \int_0^t \int_{\Omega \times \mathbb{S}^{d-1}} |\mathbf{q}_h(s)|^2 d\mathbf{x}d\mathbf{v} ds \le e^t \int_{\Omega \times \mathbb{S}^{d-1}} |f_h(0)|^2 d\mathbf{x}d\mathbf{v}. $$

The estimates on $\rho_h(t)$ and $\rho_h \mathbf{u}_h(t)$ now easily follow since $\mathbf{v} \in \mathbb{S}^{d-1}$ and by application of the Cauchy-Schwarz inequality

$$ \left\{ \begin{array}{l} \|\rho_h(t)\|_{L^2(\Omega)} \le \text{vol}(\mathbb{S}^{d-1})^{1/2} \|f_h(t)\|_{L^2}, \\ \|\rho_h \mathbf{u}_h(t)\|_{L^2(\Omega)} \le \text{vol}(\mathbb{S}^{d-1})^{1/2} \|f_h(t)\|_{L^2}. \end{array} \right. $$

Hence, we conclude the proof of (3.8) from (3.7). $\square$

## 4. PROOF OF THEOREM 2.2

For any nonnegative integer $m$, $H^m(\Omega)$ denotes the $L^2$-Sobolev space of order $m$ with the standard Sobolev norm and for $m = 0$, we use $H^0(\Omega) = L^2(\Omega)$.

For any nonnegative integer $k$, let $\Pi_h$ be the $L^2$ projection onto $\mathcal{G}_h^k$, then we have the following classical result [13].

### 4.1. Basic results.

**Lemma 4.1** (Approximation properties). There exists a constant $C > 0$, such that for any $g \in H^{m+1}(\Omega)$, the following hold:

$$ (4.1) \qquad \|g - \Pi^m g\|_{L^2(K)} + h_K^{1/2} \|g - \Pi^m g\|_{L^2(\partial K)} \le C h_K^{m+1} \|g\|_{H^{m+1}(K)}, \quad \forall K \in \mathcal{T}_h, $$

where the constant $C$ is independent of the mesh sizes $h_K$ but depends on $m$ and the shape regularity parameters $\sigma_\mathbf{x}$ and $\sigma_\mathbf{v}$ of the mesh.

Moreover, we also remind the classical inverse inequality [13]
**Lemma 4.2** (Inverse inequality). There exists a constant $C > 0$, such that for any $g \in P^m(K)$ or $P^m(K_x) \times P^m(K_v)$ with $K = (K_x \times K_v) \in \mathcal{T}_h$, the following holds:

$$ \|\nabla_x g\|_{L^2(K)} \le C h_{K_x}^{-1} \|g\|_{L^2(K)}, \quad \|\nabla_v g\|_{L^2(K)} \le C h_{K_v}^{-1} \|g\|_{L^2(K)}, $$

where the constant $C$ is independent of the mesh sizes $h_{K_x}$, $h_{K_v}$, but depends on $m$ and the shape regularity parameters $\sigma_x$ and $\sigma_v$ of the mesh.

Now let us start the error estimate analysis and consider $(f, \mathbf{q} = \nabla_v f)$ the exact solution to the kinetic equation (1.3) and $(f_h, \mathbf{q}_h)$ the approximated solution given by (2.3)-(2.8).

We introduce the consistency error function $\delta_h$ and the projected error $\xi_h \in \mathcal{G}_h^k \times \mathcal{U}_h^k$ such that

$$ (4.2) \qquad \left\{ \begin{array}{l} \xi_h = (\xi_{h,1}, \xi_{h,2}) = (f - \Pi_h f, \mathbf{q} - \Pi_h \mathbf{q}), \\ \\ \delta_h = (\delta_{h,1}, \delta_{h,2}) = (\Pi_h f - f_h, \Pi_h \mathbf{q} - \mathbf{q}_h), \end{array} \right. $$

where $\Pi_h$ represents the $L^2$ projection onto $\mathcal{G}_h^k$ and $\mathcal{U}_h^k$. Hence, we define the total error $\varepsilon_h = \delta_h + \xi_h$.

We are now ready to prove the following Lemma

**Lemma 4.3** (Estimate of $\delta_h$). Consider the numerical solution $(f_h, \mathbf{q}_h) \in \mathcal{G}_h^k \times \mathcal{U}_h^k$ for $k \ge 0$ given by (2.3)-(2.8) supplemented with periodic boundary conditions. Then for any $h_0 > 0$, there exists a constant $C > 0$ depending on $f$ and $h_0$, such that for $h \le h_0$,

$$ (4.3) \qquad \frac{1}{2} \frac{d}{dt} \|\delta_{h,1}\|_{L^2}^2 + \nu \|\delta_{h,2}\|_{L^2}^2 \le C \left[ \|\delta_{h,1}\|_{L^2}^2 + h^{2k+1} + \|\mathbf{v}_f - \mathbf{v}_{f_h}\|_{L^\infty} \|\delta_{h,1}\|_{L^2} \right]. $$

*Proof.* On the one hand the numerical approximation $(f_h, \mathbf{q}_h, \mathbf{v}_{f_h})$ given by (2.9)-(2.10) satisfies

$$ (4.4) \qquad \left\{ \begin{array}{l} a_h(f_h, \mathbf{q}_h, \mathbf{v}_{f_h}, g) = 0, \quad \forall g \in \mathcal{G}_h^k, \\ b_h(f_h, \mathbf{q}_h, \mathbf{u}) = 0, \quad \forall \mathbf{u} \in \mathcal{U}_h^k. \end{array} \right. $$

On the other hand since the numerical fluxes of (2.3)-(2.8) are consistent, the exact solution $(f, \mathbf{q}, \mathbf{v}_f)$ satisfies

$$ (4.5) \qquad \left\{ \begin{array}{l} a_h(f, \mathbf{q}, \mathbf{v}_f, g) = 0, \quad \forall g \in \mathcal{G}_h^k, \\ b_h(f, \mathbf{q}, \mathbf{u}) = 0, \quad \forall \mathbf{u} \in \mathcal{U}_h^k. \end{array} \right. $$

Then we notice that $\delta_h \in \mathcal{G}_h^k \times \mathcal{U}_h^k$; by taking $g = \delta_{h,1}$ and $\mathbf{u} = \delta_{h,2}$ in (4.4) and (4.5) and subtracting the two equalities, one has

$$ (4.6) \qquad \left\{ \begin{array}{l} a_h(\delta_h, \mathbf{v}_{f_h}, \delta_{h,1}) = -a_h(\xi_h, \mathbf{v}_{f_h}, \delta_{h,1}) + R_0, \\ b_h(\delta_h, \delta_{h,2}) + b_h(\xi_h, \delta_{h,2}) = 0, \end{array} \right. $$

where $R_0$ contains the nonlinear terms

$$ R_0 := \int_{\Omega \times \mathbb{S}^{d-1}} f \mathbf{P}_{\mathbf{v}^{\perp}} (\mathbf{v}_f - \mathbf{v}_{f_h}) \cdot \nabla_{\mathbf{v}} \delta_{h,1} \,dx\,d\mathbf{v} \\ + \sum_{\sigma_{\mathbf{v}} \in \mathcal{E}_{\mathbf{v}}} \int_{\sigma_{\mathbf{v}}} f \mathbf{P}_{\mathbf{v}^{\perp}} \widehat{(\mathbf{v}_f - \mathbf{v}_{f_h})} \cdot \mathbf{n}_{\mathbf{v}} [\delta_{h,1}]_{\mathbf{v}} \,ds_{\mathbf{v}}\,d\mathbf{x}. $$
Following the same lines as in the proof of Lemma 3.2 and using the definition of $b_h$, we get

$$
\begin{align*}
a_h(\delta_h, \mathbf{v}_{f_h}, \delta_{h,1}) &= \frac{1}{2} \frac{d}{dt} \int_{\Omega \times \mathbb{S}^{d-1}} |\delta_{h,1}|^2 d\mathbf{x}d\mathbf{v} + \nu \int_{\Omega \times \mathbb{S}^{d-1}} |\delta_{h,2}|^2 d\mathbf{x}d\mathbf{v} \\
&\quad + \frac{1}{2} \sum_{\sigma_{\mathbf{v}} \in \mathcal{E}_{\mathbf{v}}} \int_{\sigma_{\mathbf{v}}} (|\mathbf{P}_{\mathbf{v}^{\perp}} \mathbf{v}_{f_h} \cdot \mathbf{n}_{\mathbf{v}}| + \nu C_{11}) [\delta_{h,1}]_{\mathbf{v}}^2 + \nu C_{22} [\delta_{h,2}]_{\mathbf{v}}^2 ds_{\mathbf{v}} d\mathbf{x} \\
&\quad + \frac{1}{2} \sum_{\sigma_{\mathbf{x}} \in \mathcal{E}_{\mathbf{x}}} \int_{\sigma_{\mathbf{x}}} |\mathbf{v} \cdot \mathbf{n}_{\mathbf{x}}| [\delta_{h,1}]_{\mathbf{x}}^2 ds_{\mathbf{x}} d\mathbf{v} \\
&\quad + \frac{1}{2} \int_{\Omega \times \mathbb{S}^{d-1}} |\delta_{h,1}|^2 \mathrm{div}_{\mathbf{v}} (\mathbf{P}_{\mathbf{v}^{\perp}} \mathbf{v}_{f_h}) d\mathbf{x}d\mathbf{v} + \nu b_h(\xi_h, \delta_{h,2}).
\end{align*}
$$

Therefore, using (4.6), it yields

$$
\begin{align*}
(4.7) \quad & \frac{1}{2} \frac{d}{dt} \int_{\Omega \times \mathbb{S}^{d-1}} |\delta_{h,1}|^2 d\mathbf{x} d\mathbf{v} + \nu \int_{\Omega \times \mathbb{S}^{d-1}} |\delta_{h,2}|^2 d\mathbf{x} d\mathbf{v} \\
& + \frac{1}{2} \sum_{\sigma_{\mathbf{v}} \in \mathcal{E}_{\mathbf{v}}} \int_{\sigma_{\mathbf{v}}} (|\mathbf{P}_{\mathbf{v}^{\perp}} \mathbf{v}_{f_h} \cdot \mathbf{n}_{\mathbf{v}}| + \nu C_{11}) [\delta_{h,1}]_{\mathbf{v}}^2 + \nu C_{22} [\delta_{h,2}]_{\mathbf{v}}^2 ds_{\mathbf{v}} d\mathbf{x} \\
& + \frac{1}{2} \sum_{\sigma_{\mathbf{x}} \in \mathcal{E}_{\mathbf{x}}} \int_{\sigma_{\mathbf{x}}} |\mathbf{v} \cdot \mathbf{n}_{\mathbf{x}}| [\delta_{h,1}]_{\mathbf{x}}^2 ds_{\mathbf{x}} d\mathbf{v} \\
& \leq \frac{1}{2} \int_{\Omega \times \mathbb{S}^{d-1}} |\delta_{h,1}|^2 |\mathrm{div}_{\mathbf{v}} (\mathbf{P}_{\mathbf{v}^{\perp}} \mathbf{v}_{f_h})| d\mathbf{x} d\mathbf{v} + \nu |b_h(\xi_h, \delta_{h,2})| + |R_0| + |a_h(\xi_h, \mathbf{v}_{f_h}, \delta_{h,1})|,
\end{align*}
$$

where we need to evaluate the right hand side.

On the one hand, we observe that $|\mathrm{div}_{\mathbf{v}}(\mathbf{P}_{\mathbf{v}^{\perp}}\mathbf{v}_{f_h})| \le C$, so that the first term is straightforward.

Then, we evaluate the second term $|b_h|$ in (4.7). Using that $\xi_{h,2} = \nabla_{\mathbf{v}}f - \Pi_h \nabla_{\mathbf{v}}f$ the integral on the control volume vanishes and applying the Young inequality and Lemma 4.1, we have

$$
\begin{align*}
(4.8) \quad |b_h(\xi_h, \delta_{h,2})| & \le \sum_{\sigma_{\mathbf{v}} \in \mathcal{E}_{\mathbf{v}}} \int_{\sigma_{\mathbf{v}}} |\widehat{\xi}_{h,1} \mathbf{n}_{\mathbf{v}} \cdot [\delta_{h,2}]_{\mathbf{v}}| d\mathbf{x} ds_{\mathbf{v}}, \\
& \le \frac{C_{22}}{2} \sum_{\sigma_{\mathbf{v}} \in \mathcal{E}_{\mathbf{v}}} \int_{\sigma_{\mathbf{v}}} [\delta_{h,2}]_{\mathbf{v}}^2 d\mathbf{x} ds_{\mathbf{v}} + \frac{C}{2C_{22}} \|f\|_{H^{k+1}}^2 h^{2k+1}.
\end{align*}
$$

Now we treat the third term $|R_0|$ in (4.7) and use the fact that $f$ and $\mathbf{P}_{\mathbf{v}^{\perp}}$ are continuous and $\mathbf{v}_f - \mathbf{v}_{f_h}$ does not depend on $\mathbf{v}$, hence after an integration by part and by consistency of the flux, all the integrals over $\sigma_{\mathbf{v}} \in \mathcal{E}_{\mathbf{v}}$ vanish and there exists a constant $C > 0$, only depending on $\|\nabla_{\mathbf{v}}f\|_{W^{1,\infty}}$, such that

$$
(4.9) \qquad |R_0| \le C \|\mathbf{v}_f - \mathbf{v}_{f_h}\|_{L^\infty} \|\delta_{h,2}\|_{L^2}.
$$

Finally the last term in (4.7) is $a_h(\xi_h, \mathbf{v}_{f_h}, \delta_{h,1})$, we split it in two parts

$$
(4.10) \qquad a_h(\xi_h, \mathbf{v}_{f_h}, \delta_{h,1}) = R_1 + R_2,
$$
where $R_1$ and $R_2$ are given by

$$ egin{cases} R_1 = \int_{\Omega \times \mathbb{S}^{d-1}} \left( \frac{\partial \xi_{h,1}}{\partial t} \delta_{h,1} - \xi_{h,1} \mathbf{v} \cdot \nabla_{\mathbf{x}} \delta_{h,1} - (\mathbf{P}_{\mathbf{v}^{\perp}} \mathbf{v}_{f_h} \xi_{h,1} - \nu \xi_{h,2}) \cdot \nabla_{\mathbf{v}} \delta_{h,1} \right) d\mathbf{x} d\mathbf{v}, \\ R_2 = - \sum_{\sigma_{\mathbf{x}} \in \mathcal{E}_{\mathbf{x}}} \int_{\sigma_{\mathbf{x}}} \widehat{\xi_{h,1} \mathbf{v}} \cdot \mathbf{n}_{\mathbf{x}} [\delta_{h,1}]_{\mathbf{x}} ds_{\mathbf{x}} d\mathbf{v} \\ - \sum_{\sigma_{\mathbf{v}} \in \mathcal{E}_{\mathbf{v}}} \int_{\sigma_{\mathbf{v}}} \left( \xi_{h,1} \widehat{\mathbf{P}_{\mathbf{v}^{\perp}} \mathbf{v}_{f_h}} - \nu \widehat{\xi_{h,2}} \right) \cdot \mathbf{n}_{\mathbf{v}} [\delta_{h,1}]_{\mathbf{v}} ds_{\mathbf{v}} d\mathbf{x}. \end{cases} $$

Let us first evaluate the term $R_1$ and decompose it as $R_1 = R_{11} + R_{12} + R_{13}$, with

$$ |R_{11}| := \left| \int_{\Omega \times \mathbb{S}^{d-1}} \frac{\partial \xi_{h,1}}{\partial t} \delta_{h,1} d\mathbf{x} d\mathbf{v} \right| \le C h^{k+1} \|\partial_t f\|_{H^{k+1}} \|\delta_{h,1}\|_{L^2} $$

and $R_{12}$ is

$$ R_{12} := \int_{\Omega \times \mathbb{S}^{d-1}} \xi_{h,1}(\mathbf{v} - \mathbf{v}_0) \cdot \nabla_{\mathbf{x}} \delta_{h,1} d\mathbf{x} d\mathbf{v} + \int_{\Omega \times \mathbb{S}^{d-1}} \xi_{h,1} \mathbf{v}_0 \cdot \nabla_{\mathbf{x}} \delta_{h,1} d\mathbf{x} d\mathbf{v}, $$

where $\mathbf{v}_0$ be the $L^2$ projection of the function $\mathbf{v}$ onto the piecewise constant space with respect to $\mathcal{T}_h$. Hence, by definition of the $L^2$ projection $\xi_{h,1} = f - \Pi_h f$ the last term vanishes and we have from Lemma 4.1 and since $\|\mathbf{v} - \mathbf{v}_0\|_{L^\infty} \le C h$,

$$ |R_{12}| \le C \|\mathbf{v} - \mathbf{v}_0\|_{L^\infty} h^k \|f\|_{H^{k+1}} \|\delta_{h,1}\|_{L^2} \le C h^{k+1} \|f\|_{H^{k+1}} \|\delta_{h,1}\|_{L^2}. $$

Finally, we evaluate $R_{13}$ defined as

$$ R_{13} := - \int_{\Omega \times \mathbb{S}^{d-1}} \mathbf{P}_{\mathbf{v}^{\perp}} (\mathbf{v}_{f_h} - \mathbf{v}_f) \xi_{h,1} \cdot \nabla_{\mathbf{v}} \delta_{h,1} d\mathbf{x} d\mathbf{v} \\ - \int_{\Omega \times \mathbb{S}^{d-1}} \mathbf{P}_{\mathbf{v}^{\perp}} \mathbf{v}_f \xi_{h,1} \cdot \nabla_{\mathbf{v}} \delta_{h,1} d\mathbf{x} d\mathbf{v} \\ + \nu \int_{\Omega \times \mathbb{S}^{d-1}} \xi_{h,2} \cdot \nabla_{\mathbf{v}} \delta_{h,1} d\mathbf{x} d\mathbf{v}. $$

Using the definition of the $L^2$ projection $\xi_{h,2} = \mathbf{q} - \Pi_h \mathbf{q}$, we first observe that the last term vanishes

$$ \int_{\Omega \times \mathbb{S}^{d-1}} \xi_{h,2} \cdot \nabla_{\mathbf{v}} \delta_{h,1} d\mathbf{x} d\mathbf{v} = 0, $$

then we proceed as on the estimate of $R_{12}$ by introducing the $L^2$ projection of the function $\mathbf{P}_{\mathbf{v}^{\perp}} \mathbf{v}_f$ onto the piecewise constant space with respect to $\mathcal{T}^h$, whereas we apply the Cauchy-Schwarz inequality to treat the first term. It yields that there exists a constant $C > 0$,

$$ |R_{13}| \le C \left( \|\mathbf{v}_f - \mathbf{v}_{f_h}\|_{L^\infty} \|\xi_{h,1}\|_{L^2} \|\nabla_{\mathbf{v}} \delta_{h,1}\|_{L^2} + h^{k+1} \|f\|_{H^{k+1}} \|\delta_{h,1}\|_{L^2} \right). $$

Again we apply the two Lemmas 4.2 and 4.1, which gives that

$$ |R_{13}| \le C h^k \left( \|\mathbf{v}_f - \mathbf{v}_{f_h}\|_{L^\infty} + h \right) \|f\|_{H^{k+1}} \|\delta_{h,1}\|_{L^2}. $$

Gathering these results on $R_{11}$, $R_{12}$ and $R_{13}$, we get the following estimate on $R_1$,

$$ (4.11) \qquad |R_1| \le C h^k \left( \|\mathbf{v}_f - \mathbf{v}_{f_h}\|_{L^\infty} + h \right) \|f\|_{H^{k+1}} \|\delta_{h,1}\|_{L^2}. $$
Now we want to estimate the term $R_2$ containing the fluxes such that $R_2 = R_{21} + R_{22} + R_{23}$, with

$$ egin{aligned} |R_{21}| &:= \left| \sum_{\sigma_x \in \mathcal{E}_x} \int_{\sigma_x} \widehat{\xi_{h,1}} \mathbf{v} \cdot \mathbf{n}_x [\delta_{h,1}]_x ds_x d\mathbf{v} \right| \\ &\le \frac{1}{2} \sum_{\sigma_x \in \mathcal{E}_x} \int_{\sigma_x} |\mathbf{v} \cdot \mathbf{n}_x| [\delta_{h,1}]_x^2 ds_x d\mathbf{v} + \frac{1}{2} \sum_{\sigma_x \in \mathcal{E}_x} \int_{\sigma_x} |\xi_{h,1}|^2 ds_x d\mathbf{v} \\ &\le \frac{1}{2} \sum_{\sigma_x \in \mathcal{E}_x} \int_{\sigma_x} |\mathbf{v} \cdot \mathbf{n}_x| [\delta_{h,1}]_x^2 ds_x d\mathbf{v} + C \|f\|_{H^{k+1}}^2 h^{2k+1}, \end{aligned} $$

whereas from Lemma 4.1

$$ egin{aligned} |R_{22}| &:= \left| \sum_{\sigma_v \in \mathcal{E}_v} \int_{\sigma_v} \xi_{h,1} \widehat{\mathbf{P}_{v\perp}} \mathbf{v}_{f_h} \cdot \mathbf{n}_v [\delta_{h,1}]_v ds_v d\mathbf{x} \right|, \\ &\le \frac{1}{2} \sum_{\sigma_v \in \mathcal{E}_v} \int_{\sigma_v} |\mathbf{P}_{v\perp} \mathbf{v}_{f_h} \cdot \mathbf{n}_v| [\delta_{h,1}]_v^2 ds_v d\mathbf{x} + \frac{1}{2} \sum_{\sigma_v \in \mathcal{E}_v} \int_{\sigma_v} |\mathbf{P}_{v\perp} \mathbf{v}_{f_h} \cdot \mathbf{n}_v| |\xi_{h,1}|^2 ds_v d\mathbf{x}, \\ &\le \frac{1}{2} \sum_{\sigma_v \in \mathcal{E}_v} \int_{\sigma_v} |\mathbf{P}_{v\perp} \mathbf{v}_{f_h} \cdot \mathbf{n}_v| [\delta_{h,1}]_v^2 ds_v d\mathbf{x} + C \|f\|_{H^{k+1}}^2 h^{2k+1}. \end{aligned} $$

Finally, from the Young inequality and Lemma 4.1, we get that for any $\eta > 0$,

$$ egin{aligned} |R_{23}| &:= \nu \left| \sum_{\sigma_v \in \mathcal{E}_v} \int_{\sigma_v} \widehat{\xi_{h,2}} \cdot \mathbf{n}_v [\delta_{h,1}]_v ds_v d\mathbf{x} \right|, \\ &\le \frac{\nu \eta}{2} \sum_{\sigma_v \in \mathcal{E}_v} \int_{\sigma_v} [\delta_{h,1}]_v^2 ds_v d\mathbf{x} + \frac{C \nu}{2\eta} \|f\|_{H^{k+1}}^2 h^{2k+1}. \end{aligned} $$

Gathering these results on $R_{21}$, $R_{22}$ and $R_{23}$, we get the following estimate on $R_2$ with $\eta = C_{11}$,

$$ egin{aligned} |R_2| &\le \frac{1}{2} \sum_{\sigma_x \in \mathcal{E}_x} \int_{\sigma_x} |\mathbf{v} \cdot \mathbf{n}_x| [\delta_{h,1}]_x^2 ds_x d\mathbf{v} \\ &\quad + \frac{1}{2} \sum_{\sigma_v \in \mathcal{E}_v} \int_{\sigma_v} |\mathbf{P}_{v\perp} \mathbf{v}_{f_h} \cdot \mathbf{n}_v| [\delta_{h,1}]_v^2 ds_v d\mathbf{x} \\ &\quad + \frac{C_{11} \nu}{2} \sum_{\sigma_v \in \mathcal{E}_v} \int_{\sigma_v} [\delta_{h,1}]_v^2 ds_v d\mathbf{x} + C \|f\|_{H^{k+1}}^2 h^{2k+1}. \end{aligned} $$

To conclude the proof, we consider again (4.7) and use the estimates obtained in (4.8)-(4.12), it yields

$$ egin{aligned} &\frac{1}{2} \frac{d}{dt} \int_{\Omega \times \mathbb{S}^{d-1}} |\delta_{h,1}|^2 d\mathbf{x} d\mathbf{v} + \nu \int_{\Omega \times \mathbb{S}^{d-1}} |\delta_{h,2}|^2 d\mathbf{x} d\mathbf{v} \\ &\le C \left( \|\delta_{h,1}\|_{L^2}^2 + h^{2k+1} + \left( (1+h^k)\|\mathbf{v}_f - \mathbf{v}_{f_h}\|_{L^\infty} + h^{k+1} \right) \|\delta_{h,1}\|_{L^2} \right). \end{aligned} $$

Finally for any $h_0 > 0$, it yields to the result (4.3) for $h \le h_0$. $\square$

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

Now to complete the proof of convergence it remains to estimate the error on the velocity field $\mathbf{v}_{f_h}$.
**Lemma 4.4** (Estimate of $v_{f_h}$). Consider the numerical solution $(f_h, q_h) \in \mathcal{G}_h^k \times \mathcal{U}_h^k$ for $k \ge 0$ given by (2.3)-(2.8) supplemented with periodic boundary conditions. Then there exists a constant $C > 0$ such that

$$ (4.13) \qquad \|v_f - v_{f_h}\|_{L^\infty} \le C \|\varepsilon_{h,1}\|_{L^2}. $$

*Proof.* Thanks to (3.8) in Proposition 3.3, there exists a constant $C > 0$, only depending on the dimension $d$, $\Omega$ and the support of $k$ and $\phi$ such that

$$ \|\mathbf{J}_f - \mathbf{J}_h\|_{L^\infty} + \|\mathbf{R}_f - \mathbf{R}_h\|_{L^\infty} \le C (\|\phi'\|_{L^2} + \|k\|_{L^2}) \|\varepsilon_1\|_{L^2(\Omega \times \mathbb{S}^{d-1})}. $$

Then we evaluate the error $v_f - v_{f_h}$ by

$$ |v_f - v_{f_h}| = \left| \frac{(\mathbf{J}_f + \mathbf{R}_f) - (\mathbf{J}_h + \mathbf{R}_h)}{|\mathbf{J}_f + \mathbf{R}_f|} + (\mathbf{J}_h + \mathbf{R}_h) \frac{|\mathbf{J}_h + \mathbf{R}_h| - |\mathbf{J}_f + \mathbf{R}_f|}{|\mathbf{J}_h + \mathbf{R}_h| |\mathbf{J}_f + \mathbf{R}_f|} \right|. $$

Hence from assumption (1.6) on $\mathbf{J}_f$ and $\mathbf{R}_f$, there exists a new constant $C_f > 0$ depending on the exact solution $f$, $k$, $\phi$ and $\Omega$ such that

$$ \|v_f - v_{f_h}\|_{L^\infty} \le \frac{2}{\xi_T} (\|\mathbf{J}_f - \mathbf{J}_h\|_{L^\infty} + \|\mathbf{R}_f - \mathbf{R}_h\|_{L^\infty}) \le C_f \|\varepsilon_{h,1}\|_{L^2}. $$

□

**4.2. Error estimates** $\|f - f_h\|_{L^2}$. To prove Theorem 2.2, we first obtain the $L^2$ estimate on $(f_h, q_h)$, which is a direct consequence of the stability estimate proven in Lemma 3.2: there exists a constant $C_T > 0$ such that for any $t \in [0, T]$

$$ \|f_h(t)\|_{L^2} + \|q_h(t)\|_{L^2} \le C_T. $$

The error estimate follows by applying Lemma 4.3 with the estimate on $\|v_f - v_{f_h}\|_{L^\infty}$ given in Lemma 4.4, it yields

$$ \frac{1}{2} \frac{d}{dt} \|\delta_{h,1}\|_{L^2}^2 + \nu \|\delta_{h,2}\|_{L^2}^2 \le C \left[ \|\delta_{h,1}\|_{L^2}^2 + h^{2k+1} + \|\varepsilon_{h,1}\|_{L^2} \|\delta_{h,1}\|_{L^2} \right]. $$

Then we remind that $\varepsilon_h = \delta_h + \xi_h$, where $\xi_{h,1} = f - \Pi_h f$ satisfies from Lemma 4.1

$$ (4.14) \qquad \|\xi_{h,1}\|_{L^2} = \|f - \Pi^m f\|_{L^2(K)} \le C h_K^{k+1} \|f\|_{H^{k+1}(K)}, \quad \forall K \in \mathcal{T}_h, $$

hence for any $h_0 > 0$, there exists another constant $C > 0$, depending on $f$ and $h_0$, such that $h \le h_0$ and

$$ \frac{1}{2} \frac{d}{dt} \|\delta_{h,1}\|_{L^2}^2 + \nu \|\delta_{h,2}\|_{L^2}^2 \le C \left[ \|\delta_{h,1}\|_{L^2}^2 + h^{2k+1} \right]. $$

Applying the Gronwall's Lemma, we get that there exists a constant $C_T > 0$, depending on $f$, $T$ and $h_0$, such that $h \le h_0$ and for all $t \in [0, T]$

$$ (4.15) \qquad \|\delta_{h,1}(t)\|_{L^2} + \left( \int_0^t \|\delta_{h,2}\|_{L^2}^2 ds \right)^{1/2} \le C_T h^{k+1/2}. $$

Finally gathering (4.14) and (4.15) and using the same kind of estimate as (4.14) for $\nabla_v f - \Pi_h \nabla_v f$, we get that for $h \le h_0$ and for all $t \in [0, T]$,

$$ \|f(t) - f_h(t)\|_{L^2} + \left( \int_0^t \|\nabla_v f(s) - q_h(s)\|_{L^2}^2 ds \right)^{1/2} \le C_T h^{k+1/2}. $$
5. NUMERICAL SIMULATIONS

We now present several numerical experiments and simply choose $C_{11} = C_{22} = 1$. We first propose an accuracy test to verify the order of accuracy of the method and then give two examples on creation of vortices and band formation.

5.1. **Accuracy test.** We first consider the model (1.3), where the velocity $\mathbf{v}_f$ is fixed and given by $\mathbf{v}_f = \mathbf{x}t$. The initial datum is

$$f_0(\mathbf{x}, \mathbf{v}) = \frac{1}{2\pi \nu} \exp \left( -\frac{\|\mathbf{x}\|^2}{2\nu} \right), \quad \mathbf{x} \in \Omega,$$

where the computational domain is chosen as $\Omega = [-1, 1]^2$. Hence the exact solution is given by

$$f(t, \mathbf{x}, \mathbf{v}) = \frac{1}{2\pi \nu} \exp \left( -\frac{\|\mathbf{x} - \mathbf{v}t\|^2}{2\nu} \right).$$

In the numerical simulations, uniform meshes are used, with $N$ cells in each direction. In addition, the third order TVD RungeKutta method is applied in time, with the CFL number for the upwind and alternating flux in $P_1$ and $P_2$ cases. In Tables 1 and 2, we present the error $\varepsilon_N^1$ (resp. $\varepsilon_N^\infty$) on the exact solution for $L^1$ (resp. $L^\infty$) norm for $k = 1$ and $2$ with

$$\varepsilon_N^1(t) = \int_{\Omega \times \mathbb{S}^1} |f(t) - f_h(t)| d\mathbf{x} d\mathbf{v}, \quad \varepsilon_N^\infty(t) = \sup_{(\mathbf{x}, \mathbf{v}) \in \Omega \times \mathbb{S}^1} |f(t) - f_h(t)|.$$

We observe that the schemes with the upwind and alternating fluxes achieve optimal $(k+1)$-th order accuracy in approximating the solution compared to $(k+1/2)$-th order of accuracy established in the previous section.

<table><caption>TABLE 1. **Accuracy test.** Error norm $\varepsilon_N^1$ and $\varepsilon_N^\infty$ for $k = 1$ where $N$ represents the number of points in each direction.</caption><thead><tr><th></th><th>N</th><th>L<sup>1</sup> error</th><th>order</th><th>L<sup>∞</sup> error</th><th>order</th></tr></thead><tbody><tr><td rowspan="5">k = 1</td><td>16</td><td>3.09753e-00</td><td>–</td><td>1.74318e-00</td><td>–</td></tr><tr><td>24</td><td>1.57566e-00</td><td>1.7</td><td>9.76815e-01</td><td>1.45</td></tr><tr><td>32</td><td>9.21703e-01</td><td>1.7</td><td>5.97741e-01</td><td>1.54</td></tr><tr><td>48</td><td>4.03124e-01</td><td>1.95</td><td>2.73951e-01</td><td>1.83</td></tr><tr><td>64</td><td>2.18455e-01</td><td>2.00</td><td>1.51392e-01</td><td>1.98</td></tr></tbody></table>

<table><caption>TABLE 2. **Accuracy test.** Error norm $\varepsilon_N^1$ and $\varepsilon_N^\infty$ for $k = 2$ where $N$ represents the number of points in each direction.</caption><thead><tr><th></th><th>N</th><th>L<sup>1</sup> error</th><th>order</th><th>L<sup>∞</sup> error</th><th>order</th></tr></thead><tbody><tr><td rowspan="5">k = 2</td><td>16</td><td>8.61814e-01</td><td>–</td><td>5.67255e-01</td><td>–</td></tr><tr><td>24</td><td>2.71516e-01</td><td>2.85</td><td>1.95118e-01</td><td>2.63</td></tr><tr><td>32</td><td>1.13208e-01</td><td>2.93</td><td>8.55585e-02</td><td>2.73</td></tr><tr><td>48</td><td>3.18561e-02</td><td>3.09</td><td>2.53507e-02</td><td>2.94</td></tr><tr><td>64</td><td>1.41038e-02</td><td>3.00</td><td>1.05810e-02</td><td>3.01</td></tr></tbody></table>
5.2. **Taylor-Green vortex problem.** We now consider the model (1.3)-(1.4) with periodic boundary conditions in $\Omega = (0, 10)^2$, where the velocity $\mathbf{v}_f$ is given as $\mathbf{R}_f \equiv 0$ and

$$ \mathbf{J}_f(t, \mathbf{x}) = \int_{\Omega \times \mathbb{S}^{d-1}} k(|\mathbf{x} - \mathbf{x}'|) \mathbf{v}' f(t, \mathbf{x}', \mathbf{v}') d\mathbf{x}' d\mathbf{v}', $$

with $k(r) = \exp(-r^2/(2\sigma^2))$ and $\sigma = 0.1$. Here we neglect the repulsion force $\mathbf{R}_f$ and only take into account the alignment of particles with averaged velocity $\mathbf{v}_f$.

We compare the numerical solutions provided by the local discontinuous Galerkin method with the one obtained with the particle method in [21, 22]. The initial data are

$$ f_0(\mathbf{x}, \mathbf{v}) = \rho_0 (2 + v_x \Omega_x(\mathbf{x}) + v_y \Omega_y(\mathbf{x})), $$

where $\mathbf{v} = (\cos \theta, \sin \theta)$, and

$$ \left\{ \begin{aligned} \Omega_x &= +\frac{1}{3} \left[ \sin\left(\frac{\pi x}{5}\right) \cos\left(\frac{\pi y}{5}\right) + \sin\left(\frac{3\pi x}{10}\right) \cos\left(\frac{3\pi y}{10}\right) + \sin\left(\frac{\pi x}{2}\right) \cos\left(\frac{\pi y}{2}\right) \right] \\ \Omega_y &= -\frac{1}{3} \left[ \cos\left(\frac{\pi x}{5}\right) \sin\left(\frac{\pi y}{5}\right) + \cos\left(\frac{3\pi x}{10}\right) \sin\left(\frac{3\pi y}{10}\right) + \cos\left(\frac{\pi x}{2}\right) \sin\left(\frac{\pi y}{2}\right) \right] \end{aligned} \right. $$

with $\mathbf{x} = (x, y) \in (0, 10)^2$.

This model is supplemented by periodic boundary conditions in both directions. The numerical parameters for the kinetic model (1.3)-(1.4) are: $\Delta x = \Delta y = 0.2$, $\Delta t = 0.01$. In Figure 1, we report the density $\rho$ and the flux direction $\mathbf{U}$ at different time $t \in (0, 30)$ given by

$$ \rho(t, \mathbf{x}) = \int_{\mathbb{S}^1} f(t, \mathbf{x}, \mathbf{v}) \, d\mathbf{v}, \quad \rho \mathbf{U}(t, \mathbf{x}) = \int_{\mathbb{S}^1} \mathbf{v} f(t, \mathbf{x}, \mathbf{v}) \, d\mathbf{v}. $$

We find a very good agreement with the results in [21, 22] for agent based models (1.1) and macroscopic models in spite of the quite complex structure of the solution (see Figures 4 and 5 in [21] for short time $t = 5$). In our simulation, we also present simulations for large time and observe the time evolution of vortices.

Finally we also propose in Figure 2 the time evolution of the local averaged velocity $\mathbf{v}_f$ given in (1.4) and the persistence of several vortices for large time $t \simeq 30$.

5.3. **Formation of bands problem.** We still consider the kinetic model (1.3)-(1.4) but with a different scaling for $\varepsilon > 0$,

$$ \frac{\partial f}{\partial t} + \mathbf{v} \cdot \nabla_{\mathbf{x}} f = -\frac{1}{\varepsilon} \mathrm{div}_{\mathbf{v}} [\mathbf{P}_{\mathbf{v}^{\perp}} \mathbf{v}_f f - \nu \nabla_{\mathbf{v}} f], $$

where $\mathbf{v}_f$ is defined as previously. We set periodic boundary conditions in $\Omega = (-1/2, 1/2) \times (0, 1)$, and the initial data is given by

$$ f_0(\mathbf{x}, \theta) = \left(1 + \frac{1}{2} \cos(\theta)\right) \left(1 + \frac{3}{5} \sin(2\pi x) + \frac{3}{10} \cos(2\pi y)\right), \quad \mathbf{x} = (x, y) \in \Omega, \theta \in (0, 2\pi). $$

We choose $\varepsilon = 0.25$ and $\nu = 0.005$ and we investigate the long time behavior of the numerical solution. On the one hand, we report the time evolution of the density $\rho$ in Figure 3 and observe after time $t$ larger than 15, the formation of a band which propagates with an horizontal velocity of speed $\simeq 1$. Such a behaviour has been already observed for numerical simulations of stochastic models with only local alignment interactions [37] as (1.1). These moving structures appear for large enough systems after some transient. Then, they extend transversally with
<figure><img src="image_3.png" /><figcaption>FIGURE 1. **Taylor-Green vortex problem.** Numerical solution obtained with a (a) density $\rho$, (b) mean velocity **u** at time $t = 5$, $t = 15$ and $t = 30$.</figcaption></figure>

respect to the mean direction of motion. The advantage of kinetic models as (1.3)-(1.4) is that bands can be described quantitatively through local quantities, such as the local density $\rho$, but also the mena velocity **u** and the local averaged mean velocity $\mathbf{v}_f$.

We finally propose in Figure 4, a snapshot of the mean velocity **u** and the local averaged velocity $\mathbf{v}_f$ given in (1.4) at the final time $t = 30$.
<figure><img src="image_4.png" /></figure>

$t = 5$

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

$t = 15$

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

$t = 30$

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

$t = 50$

FIGURE 2. Taylor-Green vortex problem. Numerical solution $v_f$ at time $t = 5, t = 15, t = 20$ and $t = 30$.

## 6. CONCLUSION AND PERSPECTIVE

In this paper we proposed a discontinuous Galerkin discretization technique for a kinetic model of self-alignment introduced in [18, 17]. The main feature of this approach is to guarantee the accuracy and stability for the $L^2$ norm. Furthermore, we performed a complete analysis to get error estimates for smooth solutions. The scheme has been tested using an exact solution where the order of accuracy has been verified. The proposed method has been applied to study the long time dynamics of this system and can be further investigated to improve the model.

## ACKNOWLEDGEMENTS

Francis Filbet acknowledges the Division of Applied Mathematics, Brown University for the invitation in January/February 2016, where the present work has been initiated. The research of Chi-Wang Shu is partially supported by DOE grant DE-FG02-08ER25863 and NSF grant DMS-1418750.
<figure><img src="image_8.png" /><figcaption>FIGURE 3. Formation of bands problem. Numerical solution of the density $\rho$ at time $t = 2, t = 9, t = 13, t = 17.25, t = 22.75$ and $t = 30$.</figcaption></figure>

REFERENCES

[1] Aldana, M., Huepe, C.: Phase transitions in self-driven many-particle systems and related non-equilibrium models: a network approach. J. Stat. Phys., **112** (2003), pp. 135–153.

[2] Aoki, I.: A simulation study on the schooling mechanism in fish. Bulletin of the Japan Society of Scientific Fisheries, **48** (1982), pp. 1081–1088.

[3] Ayuso de Dios, B., Carrillo, J.A., Shu, C.-W.: *Discontinuous Galerkin methods for the multi-dimensional Vlasov-Poisson problem*. Mathematical Models and Methods in Applied Sciences (*M*³AS), **22** (2012), 1250042.
<figure><img src="image_9.png" /><figcaption>FIGURE 4. Formation of bands problem. Numerical solution (a) **u** and (b) **v** at time $t = 30$.</figcaption></figure>

[4] Bolley, F., Cañizo, J.A., Carrillo, J.A.: Mean-field limit for the stochastic Vicsek model. arXiv preprint 1102.1325.

[5] Bonabeau, E., Dorigo, M., Theraulaz, G.: Intelligence: From Natural to Artificial Systems. Oxford University Press, New York, 1999.

[6] Camazine, S., Deneubourg, J.-L., Franks, N.R., Sneyd, J., Theraulaz, G., Bonabeau, E.: Self-Organization in Biological Systems. Princeton University Press, 2003.

[7] Cañizo, J.A., Carrillo, J.A., Rosado, J.: A well-posedness theory in measures for some kinetic models of collective motion, Math. Mod. Meth. Appl. Sci., **21** (2011), pp. 515-539.

[8] Cañizo, J.A., Carrillo, J.A., Rosado, J.: Collective Behavior of Animals: Swarming and Complex Patterns. Arbor, **186** (2010), pp. 1035-1049.

[9] Carrillo, J.A., D'Orsogna, M.R., Panferov, V.: Double milling in self-propelled swarms from kinetic theory. Kinetic and Related Models, **2** (2009), pp. 363-378.

[10] Carrillo, J.A., Fornasier, M., Toscani, G., Vecil, F.: Particle, Kinetic, and Hydrodynamic Models of Swarming. In Naldi, G., Pareschi, L., Toscani, G. (eds.), Mathematical Modeling of Collective Behavior in Socio-Economic and Life Sciences, Series: Modelling and Simulation in Science and Technology, Birkhauser, 2010, pp. 297-336.

[11] Cheng, Y., Gamba, I.M., Majorana, A., Shu, C.-W.: A discontinuous Galerkin solver for Boltzmann-Poisson systems in nano devices. Computer Methods in Applied Mechanics and Engineering, **198** (2009), pp. 3130-3150.

[12] Chuang, Y.L., D'Orsogna, M.R., Marthaler, D., Bertozzi, A.L., Chayes, L.: State transitions and the continuum limit for a 2D interacting, self-propelled particle system. Physica D, **232** (2007), pp. 33-47.

[13] Ciarlet, P.-G. The Finite Element Methods for Elliptic Problems. North-Holland, Amsterdam, 1975.

[14] Cockburn, B., Shu, C.-W.: Runge-Kutta discontinuous Galerkin methods for convection-dominated problems. Journal of Scientific Computing, **16** (2001), pp. 173-261.

[15] Couzin, I.D., Krause, J., Franks, N.R., Levin, S.A.: Effective leadership and decision making in animal groups on the move. Nature, **433** (2005), pp. 513-516.

[16] Degond, P., Frouvelle, A., Liu, J.-G.: Macroscopic limits and phase transition in a system of self-propelled particles. J. Nonlinear Sci., **23** (2013), pp. 427-456.

[17] Degond, P., Liu, J.-G., Motsch S., Panferov, V.: Hydrodynamic models of self-organized dynamics: derivation and existence theory. Methods Appl. Anal. **20** (2013), pp. 89-114.

[18] Degond, P., Motsch, S.: Continuum limit of self-driven particles with orientation interaction. Math. Models Methods Appl. Sci., **18** (2008), pp. 1193-1215.

[19] Degond, P., Motsch, S.: A macroscopic model for a system of swarming agents using curvature control. J. Stat. Phys., (2011), available online (DOI 10.1007/s10955-011-0201-3).

[20] Degond, P., Yang, T.: Diffusion in a continuum model of self-propelled particles with alignment interaction, Math. Models Methods Appl. Sci., 20, Suppl. (2010), pp. 1459-1490.
[21] Degond, P., Dimarco, G., Mac, Th. B.N., Wang, N.: Macroscopic models of collective motion with repulsion. Commun. Math. Sci. **13** (2015), 1615–1638.

[22] Dimarco, G. and Motsch, S.: *Self-alignment driven by jump processes: Macroscopic limit and numerical investigation*, Math. Models Methods Appl. Sci. **26** (2016), no. 7, 13851410

[23] D'Orsogna, M.R., Chuang, Y.L., Bertozzi, A.L., Chayes, L.: *Self-propelled particles with soft-core interactions: patterns, stability, and collapse*. Phys. Rev. Lett., **96** (2006), 104302

[24] Doi M., Edwards, S.F.: The Theory of Polymer Dynamics. Clarendon Press, 1999.

[25] F. Filbet, Convergence of a finite volume scheme for the one dimensional Vlasov-Poisson system, *SIAM J. Numer. Analysis*, **39**, pp. 1146–1169 (2001).

[26] Filbet, F.; Russo, G. *High order numerical methods for the space non-homogeneous Boltzmann equation*, J. Comput. Phys. **186** (2003), no. 2, 457480

[27] Filbet, F. and Shu, C.-W. *Approximation of hyperbolic models for chemosensitive movement*, SIAM J. Sci. Comput. **27** (2005), no. 3, 850872.

[28] Filbet, F. and Mouhot, C. *Analysis of spectral methods for the homogeneous Boltzmann equation*, Trans. Amer. Math. Soc. **363** (2011), no. 4, 19471980.

[29] Filbet, F.; Hu, J. and Jin, S. *A numerical scheme for the quantum Boltzmann equation with stiff collision terms*, ESAIM Math. Model. Numer. Anal. **46** (2012), no. 2, 443463.

[30] Filbet, F. and Yang, C., *An inverse Lax-Wendroff method for boundary conditions applied to Boltzmann type models*, J. Comput. Phys. **245** (2013), 4361.

[31] Filbet, F. and Yang, C. *Numerical simulations of kinetic models for chemotaxis*, SIAM J. Sci. Comput. **36** (2014), no. 3, B348B366.

[32] Frouvelle, A.: *A continuous model for alignment of self-propelled particles with anisotropy and density-dependent parameters*. preprint arXiv 0912.0594.

[33] Frouvelle, A., Liu, J.-G.: *Dynamics, in a kinetic model of oriented particles with phase transition*. preprint arXiv 1101.2380.

[34] Gamba, I. M., Haack, J. R. and Motsch, S. *Spectral method for a kinetic swarming model*, J. Comput. Phys. **297** (2015), 3246.

[35] Gottlieb, S., Shu, C.-W.: *Total variation diminishing Runge-Kutta schemes*. Math. Comput. **67**, (1998), pp. 73–85.

[36] Gottlieb, S., Shu, C.-W., Tadmor, E.: *Strong stability preserving high order time discretization methods*. SIAM Review, **43**, (2001), pp. 89–112.

[37] Grégoire, G., Chaté, H.: *Onset of collective and cohesive motion*. Phy. Rev. Lett., **92** (2004), 025702.

[38] Ha, S.-Y., Tadmor, E.: *From particle to kinetic and hydrodynamic descriptions of flocking*. Kinetic and Related Models, **1** (2008), pp. 415–435.

[39] Huth, A. and Wissel, C.: *The simulation of the movement of fish schools*. Journal of Theoretical Biology, **152** (1992), pp. 365–385.

[40] Kolokolnikov, T., Carrillo, J.A., Bertozzi, A.L., Fetecau, R., Lewis, M.: *Emergent behaviour in multi-particle systems with non-local interactions*. Phys. D **260** (2013), pp. 1–4.

[41] Kulinskii, V.L., Ratushnaya, V.I., Zvelindovsky, A.V., Bedeaux, D.: *Hydrodynamic model for a system of self-propelling particles with conservative kinematic constraints*. Europhys. Lett., **71** (2005), pp. 207–213.

[42] Mogilner, A., Edelstein-Keshet, L., Bent, L., Spiros, A.: *Mutual interactions, potentials, and individual distance in a social aggregation*. J. Math. Biol., **47** (2003), pp. 353–389.

[43] Onsager, L.: *The effects of shape on the interaction of colloidal particles*. Annals of the New York Academy of Sciences, **51** (1949), pp. 627–659.

[44] Parrish, J., Edelstein-Keshet, L.: *Complexity, pattern, and evolutionary trade-offs in animal aggregation*. Science, **294** (1999), pp. 99–101.

[45] Ratushnaya, V.I., Bedeaux, D., Kulinskii, V.L., Zvelindovsky, A.V.: *Collective behaviour of self propelling particles with kinematic constraints; the relations between the discrete and the continuous description*. Physica A, **381** (2007), pp. 39–46.

[46] Ratushnaya, V.I., Kulinskii, V.L., Zvelindovsky, A.V., Bedeaux, D.: *Hydrodynamic model for the system of self propelling particles with conservative kinematic constraints; two dimensional stationary solutions*. Physica A, **366** (2006), pp. 107–114.

[47] Rey, T. and Tan, C. *An exact rescaling velocity method for some kinetic flocking models*, SIAM J. Numer. Anal. **54** (2016), no. 2, 641664.

[48] Shu, C.-W., Osher, S.: *Efficient implementation of essentially non-oscillatory shock-capturing schemes*. J. Comput. Phys. **77**, (1988), pp. 439–471.
[49] Sznitman, A.S.: Topics in propagation of chaos, École d'été de probabilités de Saint-Flour XIX-1989. Lecture Notes in Math, 1464:165–251, 1989.

[50] Topaz, C.M., Bertozzi, A.L.: Swarming patterns in a two-dimensional kinematic model for biological groups. SIAM J. Appl. Math, 65 (2004), pp. 152–174.

[51] Topaz, C.M., Bertozzi, A.L., Lewis, M.A: A nonlocal continuum model for biological aggregation. Bull. Math. Biol., 68 (2006), pp. 1601–1623.

[52] Vicsek, T., Czirók, A., Ben-Jacob, E., Cohen, I., Shochet, O.: Novel type of phase transition in a system of self-driven particles. Phys. Rev. Lett., 75 (1995), pp. 1226–1229.

FRANCIS FILBET

UNIVERSITÉ DE TOULOUSE III & IUF

INSTITUT DE MATHÉMATIQUES DE TOULOUSE,

118, ROUTE DE NARBONNE

F-31062 TOULOUSE CEDEX, FRANCE

E-MAIL: francis.filbet@math.univ-toulouse.fr

CHI-WANG SHU

DIVISION OF APPLIED MATHEMATICS

BROWN UNIVERSITY

PROVIDENCE, RI 02912, USA

E-MAIL: shu@dam.brown.edu