# On the convergence of discontinuous Galerkin/Hermite spectral methods for the Vlasov-Poisson system

Marianne Bessemoulin-Chatard, Francis Filbet

► To cite this version:

Marianne Bessemoulin-Chatard, Francis Filbet. On the convergence of discontinuous Galerkin/Hermite spectral methods for the Vlasov-Poisson system. SIAM Journal on Numerical Analysis, 2023, 61 (4), pp.1664-1688. 10.1137/22M1518232. hal-03760711v2

HAL Id: hal-03760711

https://hal.science/hal-03760711v2

Submitted on 30 Jan 2023

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_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# ON THE CONVERGENCE OF DISCONTINUOUS GALERKIN/HERMITE SPECTRAL METHODS FOR THE VLASOV-POISSON SYSTEM

MARIANNE BESSEMOULIN-CHATARD AND FRANCIS FILBET

ABSTRACT. We prove the convergence of discontinuous Galerkin approximations for the Vlasov-Poisson system written as an hyperbolic system using Hermite polynomials in velocity. To obtain stability properties, we introduce a suitable weighted $L^2$ space, with a time dependent weight, and first prove global stability for the weighted $L^2$ norm and propagation of regularity. Then we prove error estimates between the numerical solution and the smooth solution to the Vlasov-Poisson system.

## CONTENTS

<table><tr><td>1. Introduction</td><td>1</td></tr><tr><td>2. Discontinuous Galerkin/Hermite spectral methods</td><td>3</td></tr><tr><td>2.1. Hermite spectral form</td><td>3</td></tr><tr><td>2.2. Spatial discretization</td><td>4</td></tr><tr><td>2.3. Discussion on the scaling function α and main results</td><td>6</td></tr><tr><td>3. Preliminary results</td><td>8</td></tr><tr><td>3.1. Propagation of the weighted Sobolev norms</td><td>8</td></tr><tr><td>3.2. Projection error of the Hermite decomposition</td><td>9</td></tr><tr><td>3.3. Projection error of the discontinuous Galerkin methods</td><td>12</td></tr><tr><td>3.4. Global projection error</td><td>12</td></tr><tr><td>4. Proof of Theorem 2.4</td><td>14</td></tr><tr><td>5. Numerical simulations</td><td>18</td></tr><tr><td>5.1. Test 1: order of convergence</td><td>18</td></tr><tr><td>5.2. Test 2: Bump-on-the-tail</td><td>19</td></tr><tr><td>6. Conclusion & Perspectives</td><td>21</td></tr><tr><td>Acknowledgement</td><td>21</td></tr><tr><td>References</td><td>21</td></tr></table>

## 1. INTRODUCTION

We consider a noncollisional plasma of charged particles (electrons and ions). For simplicity, we assume that the properties of the plasma are one dimensional and we take into account only the electrostatic forces, thus neglecting the electromagnetic effects. We denote by $f = f(t, x, v)$ the electron distribution function and by $E(t, x)$ the electrostatic field. The Vlasov-Poisson equations of the plasma in dimensionless variables can be rewritten as,

$$ (1.1) \qquad \left\{ \begin{array}{l} \displaystyle \frac{\partial f}{\partial t} + v \frac{\partial f}{\partial x} + E \frac{\partial f}{\partial v} = 0, \\[1.5ex] \displaystyle \frac{\partial E}{\partial x} = \rho - \rho_0, \\[1.5ex] \displaystyle f(t=0) = f_0, \end{array} \right. $$

<footnote n="2010">2010 Mathematics Subject Classification. Primary: 65N35, 82C40, Secondary: 65N08.
Key words and phrases. Convergence; Discontinuous Galerkin method; Hermite spectral method; Vlasov-Poisson.</footnote>
where the density $\rho$ is given by

$$ \rho(t, x) = \int_{\mathbb{R}} f(t, x, v) \, dv, \quad t \ge 0, x \in \mathbb{T}. $$

To ensure the well-posedness of the Poisson problem, we add the compatibility (or normalizing) condition

$$ (1.2) \qquad \int_{\mathbb{T}} \rho(t, x) \, dx = \int_{\mathbb{T}} \int_{\mathbb{R}} f(t, x, v) \, dv \, dx = \operatorname{mes}(\mathbb{T}) \rho_0, \quad \forall t \ge 0, $$

which is the condition for total charge neutrality. Let us notice that (1.2) express that the total charge of the system is preserved in time.

There is a wide variety of techniques to discretize the Vlasov-Poisson system. For instance, particle methods (PIC) consist in approximating the distribution function by a finite number of Dirac masses [6]. They allow to obtain satisfying results with a small number of discrete particles, hence these methods are very popular in the community of computational plasma physics, but a well-known drawback of this approach is the inherent numerical noise which only decreases in $1/\sqrt{N}$ when the number of discrete particles $N$ increases, preventing from getting an accurate description of the distribution function for some specific applications. To overcome this difficulty, Eulerian solvers have been applied. These methods discretize the Vlasov equation on a mesh of the phase space [14, 11, 35, 15]. Among them, we can mention finite volume methods [13] which are a simple and inexpensive option, but in general low order. Fourier-Fourier transform schemes [23] are based on a Fast Fourier Transform of the distribution function in phase space, but suffer from Gibbs phenomenon if other than periodic conditions are considered. Standard finite element methods [40, 41] have also been applied, but may present numerical oscillations when approximating the Vlasov equation. Later, semi-Lagrangian schemes have also been proposed [36], consisting in computing the distribution function at each grid point by following the characteristic curves backward. Despite these schemes can achieve high order allowing also for large time steps, they require high order interpolation to compute the origin of the characteristics, destroying the local character of the reconstruction. Many improvement have been proposed and studied to make this approach more efficient [4, 8, 39]. Finally, spectral Galerkin and spectral collocation methods for the asymmetric weighted Fourier-Hermite discretization have been proposed in [12, 25, 29]. In [7], the authors study a time implicit method allowing the exact conservation of charge, momentum and energy, and highlight that for some test cases, this scheme can be significantly more accurate than the PIC method.

In the present article, we focus on a class of Eulerian methods based on Hermite polynomials in the velocity variable, where the Vlasov-Poisson system (1.1) is written as an hyperbolic system. This idea of using Galerkin methods with a small finite set of orthogonal polynomials rather than discretizing the distribution function in velocity space goes back to the 60's [1, 22]. More recently, the merit to use rescaled orthogonal basis like the so-called scaled Hermite basis has been shown [12, 19, 34, 32, 37]. In [19], Holloway formalized two possible approaches. The first one, called symmetrically-weighted, is based on standard Hermite functions as the basis in velocity and as test functions in the Galerkin method. It appears that this symmetrically weighted method cannot simultaneously conserve the mass and the momentum. It makes up for this deficiency by correctly conserving the $L^2$ norm of the distribution function, ensuring the stability of the method. In the second approach, called asymmetrically-weighted, another set of test functions is used, leading to the simultaneous conservation of mass, momentum and total energy since the infinite hyperbolic system corresponds to the one satisfied by the moments of the distribution in the velocity space. However, this approach does not conserve the $L^2$ norm of the distribution function and is then not numerically stable. Recently in [5], we provide a stability analysis of the asymmetric Hermite method in a weighted $L^2$ space for the Vlasov-Poisson system. The main idea is to introduce a scaling function $t \mapsto \alpha(t)$ which is well adapted to the variation of the distribution function with respect to time. The aim of this work is to present a convergence analysis with error estimates based on the asymmetric weighted Hermite method with a discontinuous Galerkin method for the space discretization. It is worth to mention that the convergence of the symmetric weighted Fourier-Hermite method has been already studied in
[30] where the standard $L^2$ framework is well adapted. In [24], the authors study the conservation and $L^2$ stability properties of a generalized Hermite-Fourier semi-discretization, including as special cases the symmetric and asymmetric weighted approaches. Concerning discontinuous Galerkin methods, they are similar to finite elements methods but use discontinuous polynomials and are particularly well-adapted to handling complicated boundaries which may arise in many realistic applications. Due to their local construction, this type of methods provides good local conservation properties without sacrificing the order of accuracy. They were already used for the Vlasov-Poisson system in [18, 9]. Optimal error estimates and study of the conservation properties of a family of semi-discrete DG schemes for the Vlasov-Poisson system with periodic boundary conditions have been proved for the one [2] and multi-dimensional [3] cases. In all these works, the discontinuous Galerkin method is employed using a phase space mesh.

Here, we adopt this approach only in physical space, as in [16], with a Hermite approximation in the velocity variable. In [16], such schemes with discontinuous Galerkin spatial discretization are designed in such a way that the conservation of mass, momentum and total energy is rigorously provable. In the next section, we introduce the formulation of the Vlasov equation using the Hermite basis in velocity and a class of spatial discretizations based on discontinuous Galerkin approximations. Then we present our main result on error estimates (Theorem 2.4). In Section 3, we prove some preliminary results on approximation theory based on Spectral accuracy of Hermite spectral methods and remind some basic results on interpolation error for discontinuous Galerkin method. Then in Section 4, we prove an error estimate between the semi-discrete solution (time is continuous) and the exact smooth solution of the Vlasov-Poisson system. Finally in Section 5 we present numerical results in order to illustrate the order of convergence and the stability of our approach.

## 2. DISCONTINUOUS GALERKIN/HERMITE SPECTRAL METHODS

In this section, we present the Discontinuous Galerkin/Hermite spectral method. On the one hand, we focus on the velocity discretization by expanding the distribution function $f$ using Hermite polynomials. Then, we treat the space discretization using a discontinuous Galerkin method [16, 5].

### 2.1. Hermite spectral form.
For a given scaling positive function $t \mapsto \alpha(t)$ which will be determined later, we define the weight as

$$ (2.1) \qquad \omega(t, v) := \sqrt{2\pi} \exp \left( \frac{\alpha^2(t) |v|^2}{2} \right), $$

and the associated weighted $L^2$ space

$$ L^2(\omega(t) \, d\nu) := \left\{ g : \mathbb{R} \to \mathbb{R} : \int_{\mathbb{R}} |g(v)|^2 \omega(t, v) \, d\nu < +\infty \right\}, $$

with $\langle \cdot, \cdot \rangle_{L^2(\omega(t) \, d\nu)}$ the inner product and $\| \cdot \|_{L^2(\omega(t) \, d\nu)}$ the corresponding norm. As in [5], we choose the following basis of normalized scaled time dependent asymmetrically weighted Hermite functions:

$$ (2.2) \qquad \Psi_n(t, v) = \alpha(t) H_n (\alpha(t)v) \frac{e^{-(\alpha(t)v)^2/2}}{\sqrt{2\pi}}, $$

where $\alpha$ is a scaling function depending on time and $H_n$ are the Hermite polynomials defined by $H_{-1}(\xi) = 0$, $H_0(\xi) = 1$ and for $n \ge 1$, $H_n(\xi)$ has the following recursive relation

$$ \sqrt{n} H_n(\xi) = \xi H_{n-1}(\xi) - \sqrt{n-1} H_{n-2}(\xi), \quad \forall n \ge 1. $$

Let us also emphasize that $H'_n(\xi) = \sqrt{n} H_{n-1}(\xi)$ for all $n \ge 1$, and the set of functions $(\Psi_n)_n$ defined by (2.2) is an orthogonal system satisfying

$$ (2.3) \qquad \langle \Psi_n, \Psi_m \rangle_{L^2(\omega(t) \, d\nu)} = \alpha(t) \int_{\mathbb{R}} \Psi_n(v) H_m(\alpha(t)v) \, dv = \alpha(t) \delta_{n,m}, $$

where $\delta_{n,m}$ is the Kronecker delta function. Finally, for any integer $N \ge 1$ and $t \ge 0$, we introduce the space $V_N$ as the subspace of $L^2(\omega(t) \, d\nu)$ defined by

$$ (2.4) \qquad V_N := \operatorname{Span}\{\Psi_n(t), \quad 0 \le n \le N-1\}. $$
Then we look for an approximate solution $f_N$ of (1.1) as a finite sum which corresponds to a truncation of a series

$$ (2.5) \qquad f_N(t, x, v) = \sum_{n=0}^{N-1} C_n(t, x) \Psi_n(t, v), $$

where $N$ is the number of modes and $(C_n)_{0 \le n \le N-1}$ are computed using the orthogonality property (2.3), and taking $H_n(\alpha v)$ as test function in (1.1). Therefore, a system of evolution equations is obtained for the modes $(C_n)_{0 \le n < N}$ as in [5],

$$ (2.6) \qquad \left\{ \begin{array}{l} \partial_t C_n + \mathcal{T}_n[C] = S_n[C, E_N], \\ \mathcal{T}_n[C] = \frac{1}{\alpha} (\sqrt{n} \partial_x C_{n-1} + \sqrt{n+1} \partial_x C_{n+1}), \\ S_n[C, E_N] = \frac{\alpha'}{\alpha} (n C_n + \sqrt{(n-1)n} C_{n-2}) + E_N \alpha \sqrt{n} C_{n-1}, \end{array} \right. $$

with $C_n = 0$ when $n < 0$ and $n \ge N$, and the initial data $C_n(t = 0)$ is given by

$$ C_n(t = 0) = \frac{1}{\alpha(0)} \langle f_0, \Psi_n(0) \rangle_{L^2(\omega(0) \, d\nu)} . $$

Meanwhile, we observe that the density $\rho_N$ satisfies

$$ \rho_N = \int_{\mathbb{R}} f_N \, d\nu = C_0, $$

and then the Poisson equation becomes

$$ (2.7) \qquad \frac{\partial E_N}{\partial x} = C_0 - \rho_{0,N}, $$

with $\rho_{0,N}$ such that

$$ \int_{\mathbb{T}} (C_0 - \rho_{0,N}) \, d\nu = 0. $$

Observe that when we take $N = \infty$ in the expression (2.5), we get an infinite system (2.6)-(2.7) of equations for $(C_n)_{n \in \mathbb{N}}$ and $E_N$, which is formally equivalent to the Vlasov-Poisson system (1.1).

## 2.2. Spatial discretization.
As in [5], we consider a discontinuous Galerkin approximation for the Vlasov equation with Hermite spectral basis in velocity (2.6). Let us first introduce some notations and start with $\mathcal{J} = \{0, \dots, N_x - 1\}$ describing the set of subintervals and $\hat{\mathcal{J}} = \{0, \dots, N_x\}$ related to the number of edges, where $N_x \ge 1$ is an integer, then we consider the set $\{x_{j+1/2}\}_{j \in \hat{\mathcal{J}}}$, a partition of the torus $\mathbb{T}$, where each element is denoted as $I_j = [x_{j-1/2}, x_{j+1/2}]$ with its length $h_j$ for $j \in \hat{\mathcal{J}}$, and $h = \max_j h_j$. Finally, we introduce the parameter $\delta = (h, 1/N)$ related to the numerical discretization in space and velocity.

Given any $k \in \mathbb{N}$, we define a finite dimensional discrete piecewise polynomial space

$$ (2.8) \qquad X_h = \{u \in L^2(\mathbb{T}) : u|_{I_j} \in \mathcal{P}_k(I_j), \ j \in \mathcal{J}\}, $$

where the local space $\mathcal{P}_k(I)$ consists of polynomials of degree at most $k$ on the interval $I$. We further denote the jump $[u]_{j+1/2}$ and the average $\{u\}_{j+1/2}$ of $u$ at $x_{j+1/2}$ defined as

$$ [u]_{j+1/2} = u(x_{j+1/2}^+) - u(x_{j+1/2}^-) \quad \text{and} \quad \{u\}_{j+1/2} = \frac{1}{2} \left( u(x_{j+1/2}^+) + u(x_{j+1/2}^-) \right), \quad \forall j \in \hat{\mathcal{J}}, $$

where $u(x^{\pm}) = \lim_{\Delta x \to 0^{\pm}} u(x + \Delta x)$. We also denote

$$ u_{j+1/2} = u(x_{j+1/2}), \quad u_{j+1/2}^{\pm} = u(x_{j+1/2}^{\pm}), \quad \forall j \in \hat{\mathcal{J}}. $$
From these notations, we apply a semi-discrete discontinuous Galerkin method for (2.6) as follows. We look for an approximation $C_{\delta} = (C_{\delta,n})_{0 \le n \le N-1}$ with $C_{\delta,n}(t, \cdot) \in X_h$, such that for any $\varphi_n \in X_h$, we have

$$ (2.9) \quad \frac{d}{dt} \int_{I_j} C_{\delta,n} \varphi_n \, dx + \mathcal{A}_{n,j}(g_n(C_{\delta}), \varphi_n) = \int_{I_j} \mathcal{S}_n[C_{\delta}, E_{\delta}] \varphi_n \, dx, \quad \text{for } j \in \mathcal{J}, \quad 0 \le n \le N-1, $$

where $\mathcal{A}_{n,j}$ is defined by

$$ (2.10) \qquad \left\{ \begin{array}{l} \mathcal{A}_{n,j}(g_n, \varphi_n) = - \int_{I_j} g_n \varphi'_n \, dx + \hat{g}_{n,j+\frac{1}{2}} \varphi^-_{n,j+1/2} - \hat{g}_{n,j-1/2} \varphi^+_{n,j-1/2}, \\ g_n(C_{\delta}) = \frac{1}{\alpha} \left( \sqrt{n} C_{\delta,n-1} + \sqrt{n+1} C_{\delta,n+1} \right). \end{array} \right. $$

The numerical flux $\hat{g}_n$ in (2.10) is given by

$$ (2.11) \qquad \hat{g}_n = \frac{1}{2} [g_n^-(C_{\delta}) + g_n^+(C_{\delta}) - \frac{\nu_n}{\alpha} (C_{\delta,n}^+ - C_{\delta,n}^-)], $$

with the numerical viscosity coefficient $\nu_n$ such that $\nu_n \in [\underline{\nu}, \bar{\nu}]$ with $0 < \underline{\nu} \le \bar{\nu} < \infty$.

Therefore the approximate solution of (1.1) obtained using Hermite polynomials in velocity variable and discontinuous Galerkin discretization in space is then defined by

$$ (2.12) \qquad f_{\delta}(t, x, v) = \sum_{n=0}^{N-1} C_{\delta,n}(t, x) \Psi_{n}(t, v), $$

where $\delta = (1/N, h)$ is a small parameter, $(C_{\delta,n})_n$ satisfy (2.9) and $(\Psi_n)_n$ are the basis functions defined by (2.2).

We now deal with the approximation $E_{\delta}$ of the electric field. To this end, we consider the potential function $\Phi_{\delta}(t, x)$, such that

$$ (2.13) \qquad \left\{ \begin{array}{l} E_{\delta} = -\frac{\partial \Phi_{\delta}}{\partial x}, \\ \frac{\partial E_{\delta}}{\partial x} = C_{\delta,0} - \rho_{0,\delta}. \end{array} \right. $$

Hence we get the one dimensional Poisson equation

$$ -\frac{\partial^2 \Phi_{\delta}}{\partial x^2} = C_{\delta,0} - \rho_{0,\delta}, $$

with $\rho_{0,\delta}$ such that

$$ \int_{\mathbb{T}} (C_{\delta,0} - \rho_{0,\delta}) \, dx = 0. $$

We simply consider a conforming approximation of the electric potential, corresponding to a direct integration of this Poisson problem (2.13), which is straightforward in 1D.

In what follows, we study the scheme (2.9)-(2.13), where $\alpha$ is a time-dependent function defined in the next subsection by (2.15), and not a constant scaling parameter as usual. In [5], we provide a study of the conservation properties satisfied by this type of discontinuous Galerkin/Hermite spectral methods with a time-dependent scaling function $\alpha$. It appears that the conservation properties only rely on the choice of the spatial discretization and not on the definition of $\alpha$. Indeed, we proved in [5, Proposition 2.1] the conservation of mass, momentum and total energy for the Hermite velocity discretization (2.6). Then, concerning the spatial discretization, we established in [5, Theorem 3.4] the conservation of the discrete total energy for the scheme (2.9)-(2.11) with a centered numerical flux $\hat{g}_0$ (corresponding to $\nu_0 = 0$ in (2.11)) together with a discontinuous Galerkin approximation of the Poisson equation.
2.3. **Discussion on the scaling function α and main results.** Before to state an error estimate on the numerical solution to (2.9)-(2.13), let us introduce the suitable functional framework. We set $\mu_t$ the measure given as

$$ (2.14) \qquad \mathrm{d}\mu_t = \omega(t, v) \mathrm{d}x \mathrm{d}v $$

where the weight $\omega$ is provided in (2.1) and the following $L^2$ weighted space given by

$$ L^2(\mathrm{d}\mu_t) := \left\{ g : \mathbb{T} \times \mathbb{R} \to \mathbb{R} : \iint_{\mathbb{T} \times \mathbb{R}} |g(x, v)|^2 \mathrm{d}\mu_t < +\infty \right\}, $$

with $\langle \cdot, \cdot \rangle_{L^2(\mathrm{d}\mu_t)}$ the associated inner product, that is

$$ \langle f, g \rangle_{L^2(\mathrm{d}\mu_t)} = \iint_{\mathbb{T} \times \mathbb{R}} f(x, v) g(x, v) \mathrm{d}\mu_t, $$

and $\| \cdot \|_{L^2(\mathrm{d}\mu_t)}$ the corresponding norm.

Let us remark that for all functions $\alpha, \tilde{\alpha} : \mathbb{R}_+ \to \mathbb{R}_+$ such that $\tilde{\alpha}(t) \le \alpha(t)$ for all $t \ge 0$, the associated $\omega, \tilde{\omega}$ defined by (2.1) satisfy $\tilde{\omega}(t, v) \le \omega(t, v)$ for all $t \ge 0$, $v \in \mathbb{R}$. Then, the corresponding weighted $L^2$ norms verify

$$ \iint_{\mathbb{T} \times \mathbb{R}} |f(x, v)|^2 \tilde{\omega}(t, v) \mathrm{d}x \mathrm{d}v \le \iint_{\mathbb{T} \times \mathbb{R}} |f(x, v)|^2 \omega(t, v) \mathrm{d}x \mathrm{d}v. $$

In particular, taking $\tilde{\alpha} = 0$, the standard $L^2(\mathrm{d}x \mathrm{d}v)$ norm is controlled by the weighted $L^2$ norm:

$$ \|f\|_{L^2(\mathrm{d}x \mathrm{d}v)}^2 = \iint_{\mathbb{T} \times \mathbb{R}} |f(x, v)|^2 \mathrm{d}x \mathrm{d}v \le \iint_{\mathbb{T} \times \mathbb{R}} |f(x, v)|^2 \omega(t, v) \mathrm{d}x \mathrm{d}v = \|f\|_{L^2(\mathrm{d}\mu_t)}^2. $$

The first issue is to find the appropriate framework for the stability of approximations based on asymmetrically-weighted Hermite basis. Indeed, this choice fails to preserve the $L^2$ norm of the approximate solution, and therefore to ensure the long-time stability of the method. Consequently, we introduce a $L^2$ weighted space, with a time-dependent weight, allowing to prove the global stability of the solution in this space [5]. Actually, this idea has been already employed in [27, 28] to stabilize Hermite spectral methods for linear diffusion equations and nonlinear convection-diffusion equations in unbounded domains, yielding stability and spectral convergence of the considered methods. The main point now is to determine a function $\alpha$. We proved the following result in [5, Proposition 3.2].

**Proposition 2.1.** Let $(f_\delta, E_\delta)$ be the approximate solution defined by (2.9)-(2.13), with the scaling function $\alpha$ defined by

$$ (2.15) \qquad \alpha(t) := \alpha_0 \left( 1 + \gamma \alpha_0^2 \int_0^t \max(1, \|E_\delta(s)\|_{L^\infty}^2) \mathrm{d}s \right)^{-1/2}. $$

Assume that $\|f_\delta(0)\|_{L^2(\mathrm{d}\mu_0)} < +\infty$. Then, for any $t \ge 0$, we have

$$ (2.16) \qquad \frac{\mathrm{d}}{\mathrm{d}t} \|f_\delta(t)\|_{L^2(\mathrm{d}\mu_t)}^2 := \frac{\mathrm{d}}{\mathrm{d}t} \left( \alpha(t) \sum_{n=0}^{N-1} \int_{\mathbb{T}} |C_{\delta,n}|^2 \mathrm{d}x \right) \le - \sum_{n=0}^{N-1} \sum_{j \in \mathcal{J}} \nu_n [C_{\delta,n}]_{j-\frac{1}{2}}^2 + \frac{1}{2\gamma} \|f_\delta(t)\|_{L^2(\mathrm{d}\mu_t)}^2, $$

from which we deduce

$$ (2.17) \qquad \|f_\delta(t)\|_{L^2(\mathrm{d}\mu_t)} \le \|f_\delta(0)\|_{L^2(\mathrm{d}\mu_0)} e^{t/4\gamma}, $$

where $\gamma > 0$ is the fixed parameter, which can be chosen arbitrarily, appearing in the definition (2.15) of $\alpha$.

To study the convergence of the numerical method, we also need to establish a stability result for the exact solution $f$ in the weighted norm $L^2(\mathrm{d}\mu_t)$, where the weight depends on the approximate solution (see definition (2.15) of $\alpha$).
**Proposition 2.2.** Consider $(f, E)$ a smooth solution of the Vlasov-Poisson system (1.1). Assuming that the initial data $f_0$ belongs to $L^2(d\mu_0)$, then there exists $c_0 > 0$ such that the solution $f(t)$ satisfies for all $t \ge 0$:

$$ \|f(t)\|_{L^2(d\mu_t)} \le \|f_0\|_{L^2(d\mu_0)} e^{t/4\eta}, $$

where $d\mu_t$ is defined by (2.14), with $\alpha$ appearing in the weight $\omega$ given by (2.15).

The parameter $\eta > 0$ has to be chosen small enough, namely such that $\eta C/\gamma < 1$, where $C > 0$ is a constant depending only on the mass of $f_0$ such that $\|E(t)\|_{L^\infty}^2 \le C$ for all $t \ge 0$, and $\gamma > 0$ is the fixed parameter appearing in the definition (2.15) of $\alpha$.

*Proof.* Using the Vlasov equation (1.1) and the definition (2.1) of the weight $\omega$, we have

$$ \frac{1}{2} \frac{d}{dt} \|f(t)\|_{L^2(d\mu_t)}^2 = \frac{1}{2} \iint_{\mathbb{T} \times \mathbb{R}} f^2 (\alpha^2 E v + \alpha \alpha' |v|^2) \omega \, dx \, dv. $$

Applying now the Young inequality on the first term, we get for $\eta > 0$,

$$ \frac{1}{2} \frac{d}{dt} \|f(t)\|_{L^2(d\mu_t)}^2 \le \frac{1}{2} \left( \frac{\eta}{2} \|E\|_{L^\infty}^2 \alpha^4 + \alpha' \alpha \right) \iint_{\mathbb{T} \times \mathbb{R}} f^2 |v|^2 \omega \, dx \, dv + \frac{1}{4\eta} \|f(t)\|_{L^2(d\mu_t)}^2. $$

Then, using the definition (2.15) of $\alpha$, it is clear that

$$ \alpha' = -\frac{\gamma}{2} \max(1, \|E_\delta\|_{L^\infty}^2) \alpha^3, $$

leading to

$$ \frac{1}{2} \frac{d}{dt} \|f(t)\|_{L^2(d\mu_t)}^2 \le \frac{1}{2} \left( \frac{\eta}{2} \|E\|_{L^\infty}^2 - \frac{\gamma}{2} \max(1, \|E_\delta\|_{L^\infty}^2) \right) \alpha^4 \iint_{\mathbb{T} \times \mathbb{R}} f^2 |v|^2 \omega \, dx \, dv + \frac{1}{4\eta} \|f(t)\|_{L^2(d\mu_t)}^2. $$

In one space dimension, the $L^1$ estimate on $f$ can be used to bound the electric field, namely there exists a constant $C > 0$ depending only on the initial mass of $f$ such that for all $t \ge 0$, $\|E(t)\|_{L^\infty}^2 \le C$. Then, choosing $\eta > 0$ such that $\eta C/\gamma < 1$, the first term of the right-hand side is nonpositive, which concludes the proof. $\square$

Notice that the functional space $L^2(d\mu_t)$ depends on $\alpha$ given by (2.15), and then on the discretization parameter itself through the term $\|E_\delta(t)\|_{L^\infty}$ involved in this definition. Therefore, to establish a convergence result, it is mandatory to control $\alpha$ uniformly with respect to $\delta$. This control is achieved by bounding $\|E_\delta(t)\|_{L^\infty}$ uniformly with respect to $\delta$.

**Proposition 2.3.** We consider a solution $(f_\delta, E_\delta)$ of (2.9)-(2.13), with $\alpha$ defined by (2.15). We assume that $\|f_\delta(0)\|_{L^2(d\mu_0)} < +\infty$. Let $T > 0$ be a fixed final time. Then, there exists a constant $C_T > 0$ independent of the discretization parameter $\delta$ such that

$$ \|E_\delta(t)\|_{L^\infty} \le C_T, \quad \forall t \in [0, T]. $$

Moreover, the scaling function $\alpha$ satisfies

$$ (2.18) \qquad 0 < \underline{\alpha}_T \le \alpha(t) \le \alpha_0, \quad \forall t \in [0, T], $$

where the constant $\underline{\alpha}_T$ is independent of $\delta$ and given by

$$ \underline{\alpha}_T := \alpha_0 (1 + \gamma \alpha_0^2 (1 + C_T) T)^{-1/2}. $$

*Proof.* The proof of the uniform $L^\infty$ bound on $E_\delta$ is mainly based on the Sobolev and Poincaré-Wirtinger inequalities, together with the stability estimate (2.17), as detailed in [5, Theorem 3.6]. Thanks to this bound, it is straightforward to obtain the lower bound of $\alpha$ by using the definition (2.15). The upper bound $\alpha_0$ is also clear since $\alpha$ is a nonincreasing function. $\square$

We are now in position to state our main result.
**Theorem 2.4.** For any $t \in [0, T]$, consider the scaling function $\alpha$ defined by (2.15) and let $f(t, \cdot) \in H^m(\mathbf{d}\mu_t)$ be the solution of the Vlasov-Poisson system (1.1) where $m \ge k + 1$ and $f_\delta$ be the approximation defined by (2.9)-(2.12). Then there exists a constant $C > 0$, independent of $\delta = (h, 1/N)$, but depending on $T$, such that

$$ (2.19) \qquad \|f(t) - f_\delta(t)\|_{L^2(\mathbf{d}x \, \mathbf{d}v)} \le \|f(t) - f_\delta(t)\|_{L^2(\mathbf{d}\mu_t)} \le C \left[ \frac{1}{N^{(m-1)/2}} + h^{k+1/2} \right]. $$

Before to present the proof of this result, let us make some comments.

* On the one hand this result shows that the discretization in velocity using Hermite polynomials provides spectral accuracy. On the other hand, we get the classical order of convergence of the discontinuous Galerkin method for the space discretization.

* For the sake of clarity, we do not write explicitly how the error bound (2.19) depends on the scaling parameter $\alpha$, nevertheless in the proof we will follow carefully this dependence.

### 3. PRELIMINARY RESULTS

In the next section, we provide some results about the propagation of regularity of the solution to the Vlasov-Poisson system (1.1).

**3.1. Propagation of the weighted Sobolev norms.** The global existence of classical $\mathcal{C}^1$ solutions of the Vlasov-Poisson system has been obtained by Ukai-Okabe [38] in the case of space dimension two, and by Pfaffelmoser [31] and Lions-Perthame [26] independently in three dimensions. The key to obtain classical solutions of the Vlasov-Poisson system is to prove that the macroscopic (charge) density $\rho \in L^\infty([0, T] \times \mathbb{R}^d)$ for all $T > 0$. Following Ukai & Okabe who show the decay of the distribution function $f$ in the $v$ variable formulated in terms of a convenient weighted estimate, we get the propagation of $\mathcal{C}^m$ regularity for the solution $(f, E)$. More precisely, according for example to the article of Ukai & Okabe [38], the following result holds.

Let $\beta \in \mathcal{C}^m(\mathbb{R})$ be such that

$$ \beta \ge 0, \quad \beta' \le 0, \quad \text{and} \quad \beta(r) = O\left(\frac{1}{r^\lambda}\right), $$

with $\lambda > 1$.

**Theorem 3.1.** For any $m \ge 1$, assume that $f_0$ is nonnegative such that $f_0 \in \mathcal{C}^{m-1}(\mathbb{T} \times \mathbb{R})$ and for all $0 \le k \le m$,

$$ (|\partial_x^k f_0| + |\partial_v^k f_0|) (x, v) \le \beta(|v|). $$

Then there exists a unique classical solution to the Vlasov-Poisson system (1.1) satisfying for any $0 \le k \le m-1$

$$ (|\partial_x^k f| + |\partial_v^k f|) (t, x, v) = O(1/|v|^\lambda), \quad \text{as } |v| \to \infty, $$

uniformly in $(t, x) \in [0, T] \times \mathbb{T}$. Moreover the electric field $E$ satisfies for all $1 \le k \le m$

$$ (|E| + |\partial_x^k E|) (t, x) \le C_T. $$

The proof of this result is done in the two dimensional case in [38], but also applies in the simpler one dimensional case.

This latter result allows to obtain uniform estimates on the electric field and its space derivative in $\mathcal{C}^m([0, T] \times \mathbb{T})$, hence we can propagate the weighted Sobolev norm on the solution to the Vlasov-Poisson system (1.1)

**Corollary 3.2.** Under the assumption of Theorem 3.1 and assuming that

$$ \|f_0\|_{H^m(\mathbf{d}\mu_0)} < \infty, $$

we have that for any $t \in [0, T]$

$$ \|f(t)\|_{H^m(\mathbf{d}\mu_t)} \le \|f_0\|_{H^m(\mathbf{d}\mu_0)} \exp(Ct). $$
Proof. Theorem 3.1 ensures that the electric field is such that for all $t \in (0, T]$

$$ \|E(t)\|_{W^{m,\infty}} \leq C_T. $$

Then since the electric field is uniformly bounded, we propagate the $H^m(\mathrm{d}\mu_t)$ norms of $f(t)$ as we estimate the weighted $L^2(\mathrm{d}\mu_t)$ norm in Proposition 2.2. $\square$

3.2. **Projection error of the Hermite decomposition.** In this section we present some approximation properties of the chosen Hermite functions. Since the results presented here are very similar to those proposed in [17, Section 2], we only briefly outline the proofs.

For any integer $m \ge 1$, we define

$$ H^m(\omega(t) \mathrm{d}v) := \{g : \mathbb{R} \to \mathbb{R}; \partial_v^l g \in L^2(\omega(t) \mathrm{d}v), 0 \le l \le m\}, $$

with the following seminorm and norm:

$$ |g|_{H^m(\omega(t) \mathrm{d}v)} = \|\partial_v^m g\|_{L^2(\omega(t) \mathrm{d}v)}, \quad \|g\|_{H^m(\omega(t) \mathrm{d}v)} = \left( \sum_{l=0}^{m} |g|_{H^l(\omega(t) \mathrm{d}v)}^2 \right)^{1/2}. $$

Using the definition (2.2) of $\Psi_n$ and the properties of the Hermite polynomials $H_n$, one obtains that for all $n \ge 0$,

$$ (3.1) \qquad \left\{ \begin{array}{l} \partial_t \Psi_n = -\frac{\alpha'}{\alpha} \left(n \Psi_n + \sqrt{(n+1)(n+2)} \Psi_{n+2}\right), \\ \alpha v \Psi_n = \sqrt{n+1} \Psi_{n+1} + \sqrt{n} \Psi_{n-1}, \\ \partial_v \Psi_n = -\alpha \sqrt{n+1} \Psi_{n+1}. \end{array} \right. $$

Using the latter relation, the set $(\partial_v \Psi_n)_n$ is also an orthogonal system, namely

$$ (3.2) \qquad \langle \partial_v \Psi_n, \partial_v \Psi_m \rangle_{L^2(\omega(t) \mathrm{d}v)} = \alpha^3(t) (n+1) \delta_{n,m}. $$

Therefore, any $g \in L^2(\omega(t) \mathrm{d}v)$ can be expanded as

$$ (3.3) \qquad g(v) = \sum_{n \in \mathbb{N}} \hat{g}_n(t) \Psi_n(t, v), $$

and the Hermite coefficients are given by

$$ (3.4) \qquad \hat{g}_n(t) = \frac{1}{\alpha(t)} \langle g, \Psi_n(t) \rangle_{L^2(\omega(t) \mathrm{d}v)}. $$

Remark also that the following equalities hold for every $g \in L^2(\omega(t) \mathrm{d}v)$:

$$ (3.5) \qquad \|g\|_{L^2(\omega(t) \mathrm{d}v)}^2 = \alpha(t) \sum_{n \in \mathbb{N}} |\hat{g}_n(t)|^2 = \frac{1}{\alpha(t)} \sum_{n \in \mathbb{N}} \left| \langle g, \Psi_n(t) \rangle_{L^2(\omega(t) \mathrm{d}v)} \right|^2. $$

Finally, we also introduce $\mathcal{P}_{V_N}$ the orthogonal projection on $V_N$ such that we have

$$ \langle g - \mathcal{P}_{V_N} g, \varphi \rangle_{L^2(\omega(t) \mathrm{d}v)} = 0, \quad \forall \varphi \in V_N, $$

and

$$ \mathcal{P}_{V_N} g(t, v) = \sum_{n=0}^{N-1} \hat{g}_n(t) \Psi_n(t, v). $$

Following [17], we now establish some inverse inequalities and imbedding inequalities which are needed to analyze the spectral convergence property for the here considered Hermite method.

**Lemma 3.3.** For any $g \in V_N$,

$$ |g|_{H^1(\omega(t) \mathrm{d}v)} \le \alpha(t) \sqrt{N} \|g\|_{L^2(\omega(t) \mathrm{d}v)}. $$
Proof. Decomposing $g \in V_N$ as

$$g(t, v) = \sum_{n=0}^{N-1} \hat{g}_n(t) \Psi_n(t, v),$$

and using the orthogonality properties (2.3) and (3.2) yields

$$|g|_{H^1(\omega(t) \, d\nu)}^2 = \alpha^3(t) \sum_{n=0}^{N-1} (n+1) |\hat{g}_n(t)|^2 \le \alpha^3(t) N \sum_{n=0}^{N-1} |\hat{g}_n(t)|^2 = \alpha^2(t) N \|g\|_{L^2(\omega(t) \, d\nu)}^2.$$

□

Furthermore, we also show the following inequalities.

**Lemma 3.4.** For any $g \in H^1(\omega(t) \, d\nu)$, it holds

$$(3.6) \qquad \left\{ \begin{array}{l} \|g\|_{L^2(\omega(t) \, d\nu)} \le \frac{2}{\alpha(t)} |g|_{H^1(\omega(t) \, d\nu)}, \\ \|v g\|_{L^2(\omega(t) \, d\nu)} \le \frac{2}{\alpha^2(t)} |g|_{H^1(\omega(t) \, d\nu)}. \end{array} \right.$$

Proof. Using that $\partial_v\omega(t, v) = v\alpha^2(t)\omega(t, v)$, we have

$$\int_{\mathbb{R}} (v g(v))^2 \omega(t, v) \, d\nu = \frac{1}{\alpha^2(t)} \int_{\mathbb{R}} v g(v)^2 \partial_v \omega(t, v) \, d\nu.$$

By an integration by parts, one gets

$$\int_{\mathbb{R}} (v g(v))^2 \omega(t, v) \, d\nu = -\frac{1}{\alpha^2(t)} \int_{\mathbb{R}} g(v)^2 \omega(t, v) \, d\nu - \frac{2}{\alpha^2(t)} \int_{\mathbb{R}} v g(v) \partial_v g(v) \omega(t, v) \, d\nu.$$

Then, applying the Cauchy-Schwarz inequality to the second term of the right hand side, it yields

$$(3.7) \qquad \|v g\|_{L^2(\omega(t) \, d\nu)}^2 + \frac{1}{\alpha^2(t)} \|g\|_{L^2(\omega(t) \, d\nu)}^2 \le \frac{2}{\alpha^2(t)} \|v g\|_{L^2(\omega(t) \, d\nu)} \|g|_{H^1(\omega(t) \, d\nu)},$$

from which we deduce that

$$\|v g\|_{L^2(\omega(t) \, d\nu)} \le \frac{2}{\alpha^2(t)} |g|_{H^1(\omega(t) \, d\nu)}.$$

Now, using this later estimate in (3.7), we obtain that

$$\|g\|_{L^2(\omega(t) \, d\nu)}^2 \le \frac{4}{\alpha^2(t)} |g|_{H^1(\omega(t) \, d\nu)}^2,$$

which concludes the proof.

□

Still following [17], let us define the Fokker-Planck operator $\mathcal{F}$ as

$$(3.8) \qquad \mathcal{F}[g](v) = -\partial_v (\omega^{-1}(t) \partial_v (g\omega(t))),$$

where $\omega(t)$ is the weight defined in (2.1). It follows from Lemma 3.4 that $\mathcal{F}$ is a continuous mapping from $H^2(\omega(t) \, d\nu)$ to $L^2(\omega(t) \, d\nu)$, with $\|\mathcal{F}\|$ independent of $t$, as stated in the following lemma.

**Lemma 3.5.** For all $g \in H^2(\omega(t) \, d\nu)$, it holds

$$(3.9) \qquad \|\mathcal{F}[g]\|_{L^2(\omega(t) \, d\nu)} \le 7 |g|_{H^2(\omega(t) \, d\nu)}.$$

Proof. Using the definition (3.8), we obtain after some computations that for all $g \in H^2(\omega(t) \, d\nu)$,

$$\mathcal{F}[g](v) = -\partial_v^2 g - \alpha^2(t) v \partial_v g - \alpha^2(t) g.$$

Then, applying the triangle inequality together with the second inequality of Lemma 3.4 on the second term and twice the first inequality of Lemma 3.4 on the third term, we obtain the expected estimate.

□
Moreover, $\Psi_n$ is the n-th eigenfunction of the following singular Liouville problem:

$$ (3.10) \qquad -\mathcal{F}[g](v) + \lambda g(v) = 0, \quad v \in \mathbb{R}, $$

with corresponding eigenvalues $\lambda_n = \alpha^2(t)n$.

**Proposition 3.6.** Let $r \ge 0$. For any $g \in H^r(\omega(t) \, d\nu)$, it holds for all $N \ge 0$

$$ (3.11) \qquad \| g - \mathcal{P}_{V_N} g \|_{L^2(\omega(t) \, d\nu)} \le \frac{C}{(\alpha^2(t) N)^{r/2}} \| g \|_{H^r(\omega(t) \, d\nu)}, $$

with $C > 0$ independent of $N$ and $t$.

*Proof.* Throughout this proof, let $\mathcal{C}$ be a generic positive constant independent of $N$ and $t$, which may be different in different places. Using the orthogonal relation (2.3) and the definition of the orthogonal projection $\mathcal{P}_{V_N}$, we have

$$ \| g - \mathcal{P}_{V_N} g \|_{L^2(\omega(t) \, d\nu)}^2 = \alpha(t) \sum_{n \ge N} |\hat{g}_n(t)|^2. $$

Let us first treat the case where $r$ is an even integer. By the singular Liouville equation (3.10), we get

$$ \int_{\mathbb{R}} g \Psi_n \omega(t) \, d\nu = \frac{1}{\alpha^2(t) n} \int_{\mathbb{R}} g \mathcal{F}[\Psi_n] \omega(t) \, d\nu. $$

Then, using the definition of the Fokker-Planck operator $\mathcal{F}$ in (3.8) and performing two successive integrations by parts, it holds

$$ \begin{aligned} \int_{\mathbb{R}} g \Psi_n \omega(t) \, d\nu &= \frac{1}{\alpha^2(t) n} \int_{\mathbb{R}} \omega^{-1}(t) \partial_v (\Psi_n \omega(t)) \partial_v g \omega(t) \, d\nu, \\ &= -\frac{1}{\alpha^2(t) n} \int_{\mathbb{R}} \Psi_n \partial_v (\omega^{-1}(t) \partial_v (g \omega(t))) \omega(t) \, d\nu, \\ &= \frac{1}{\alpha^2(t) n} \int_{\mathbb{R}} \Psi_n \mathcal{F}[g] \omega(t) \, d\nu. \end{aligned} $$

Then by induction, we deduce that

$$ (3.12) \qquad \int_{\mathbb{R}} g \Psi_n \omega(t) \, d\nu = \frac{1}{(\alpha^2(t) n)^{r/2}} \int_{\mathbb{R}} \mathcal{F}^{r/2}[g] \Psi_n \omega(t) \, d\nu. $$

Consequently, using (3.4), it yields

$$ (3.13) \qquad |\hat{g}_n(t)| = \frac{1}{\alpha(t)} \frac{1}{(\alpha^2(t) n)^{r/2}} \left| \int_{\mathbb{R}} \mathcal{F}^{r/2}[g] \Psi_n \omega(t) \, d\nu \right|. $$

Furthermore, using (3.5) and Lemma 3.5, we get

$$ \begin{aligned} \| g - \mathcal{P}_{V_N} g \|_{L^2(\omega(t) \, d\nu)}^2 &= \alpha(t) \sum_{n \ge N} \frac{1}{\alpha^2(t) n} \frac{1}{(\alpha^2(t) n)^r} \left| \int_{\mathbb{R}} \mathcal{F}^{r/2}[g] \Psi_n \omega(t) \, d\nu \right|^2, \\ &\le \frac{1}{(\alpha^2(t) N)^r} \| \mathcal{F}^{r/2}[g] \|_{L^2(\omega(t) \, d\nu)}^2, \\ &\le \frac{C}{(\alpha^2(t) N)^r} \| g \|_{H^r(\omega(t) \, d\nu)}^2. \end{aligned} $$
Now, let $r$ be any odd integer. Applying (3.12) for $r-1$ (which is now even), using the Liouville equation (3.10) and an integration by parts, we have

$$ 
\begin{aligned}
\int_{\mathbb{R}} g \Psi_n \omega(t) \, dv &= \frac{1}{(\alpha^2(t) n)^{(r-1)/2}} \int_{\mathbb{R}} \mathcal{F}^{(r-1)/2}[g] \Psi_n \omega(t) \, dv \\
&= \frac{1}{(\alpha^2(t) n)^{(r+1)/2}} \int_{\mathbb{R}} \partial_v (\mathcal{F}^{(r-1)2}[g] \omega(t)) \partial_v (\Psi_n \omega(t)) \omega^{-1}(t) \, dv.
\end{aligned}
$$

On the one hand, by virtue of the two last relations in (3.1), we obtain

$$ \omega^{-1}(t) \partial_v (\Psi_n \omega(t)) = \alpha(t) \sqrt{n} \Psi_{n-1}. $$

On the other hand, we remark that

$$ \omega^{-1}(t) \partial_v (\mathcal{F}^{(r-1)/2}[g] \omega(t)) = \partial_v \mathcal{F}^{(r-1)/2}[g] + \alpha^2(t) v \mathcal{F}^{(r-1)/2}[g]. $$

Using these latter results, it yields that

$$ |\hat{g}_n| = \frac{1}{\alpha(t)} \frac{1}{(\alpha^2(t) n)^{r/2}} \left| \int_{\mathbb{R}} \left( \partial_v \mathcal{F}^{(r-1)/2}[g] + \alpha^2(t) v \mathcal{F}^{(r-1)/2}[g] \right) \Psi_{n-1} \omega(t) \, dv \right| $$

Then, proceeding as above and using (3.6), we get

$$ 
\begin{aligned}
\|g - \mathcal{P}_{V_N} g\|_{L^2(\omega(t) \, dv)}^2 &\le \frac{C}{(\alpha^2(t) N)^r} \left[ \|\partial_v \mathcal{F}^{(r-1)/2}[g]\|_{L^2(\omega(t) \, dv)}^2 + \alpha^4(t) \|v \mathcal{F}^{(r-1)/2}[g]\|_{L^2(\omega(t) \, dv)}^2 \right] \\
&\le \frac{C}{(\alpha^2(t) N)^r} \|\partial_v \mathcal{F}^{(r-1)/2}[g]\|_{L^2(\omega(t) \, dv)}^2 \\
&\le \frac{C}{(\alpha^2(t) N)^r} \|g\|_{H^r(\omega(t) \, dv)}^2,
\end{aligned}
$$

which concludes the proof. $\square$

**3.3. Projection error of the discontinuous Galerkin methods.** Now, let us recall a classical result concerning the spatial approximation (see for example [21, Lemma 2.1]).

**Proposition 3.7.** There exists a constant $\bar{C} > 0$, independent of $h$, such that for any $g \in H^{k+1}(\mathrm{d}x)$, the following inequality holds:

$$ (3.14) \qquad \|g - \mathcal{P}_{X_h} g\|_{L^2}^2 + h \|g - \mathcal{P}_{X_h} g\|_{L^2(\Gamma)}^2 \le \bar{C} h^{2(k+1)} \|g\|_{H^{k+1}}^2, $$

where $\| \cdot \|_{L^2(\Gamma)}$ is the norm over the mesh skeleton $\Gamma = (x_{j+1/2})_{j \in \hat{\mathcal{J}}}$ defined by

$$ (3.15) \qquad \|u\|_{L^2(\Gamma)}^2 = \sum_{j \in \hat{\mathcal{J}}} (|u_{j-1/2}^+|^2 + |u_{j-1/2}^-|^2). $$

**3.4. Global projection error.** In this section, we provide a global projection error on the combination of Hermite polynomial approximations and local discontinuous Galerkin interpolation in the suitable functional space. For any $t \ge 0$, we introduce $W_\delta$, the subspace of $L^2(\mathrm{d}\mu_t)$ defined by

$$ W_\delta := V_N \otimes X_h, $$

where $V_N$ is given by (2.4) and $X_h$ is defined in (2.8), and $\mathcal{P}_{W_\delta}$ the orthogonal projection on $W_\delta$. For $f \in L^2(\mathrm{d}\mu_t)$ which can be expanded as

$$ (3.16) \qquad f(x, v) = \sum_{n \in \mathbb{N}} \hat{C}_n(t, x) \Psi_n(t, v), $$

the projection $\mathcal{P}_{W_\delta}f$ is then given by

$$ (3.17) \qquad \mathcal{P}_{W_\delta}f = \sum_{n=0}^{N-1} \mathcal{P}_{X_h} \hat{C}_n(t, x) \Psi_n(t, v), $$
where $\mathcal{P}_{X_h}$ is the $L^2(\mathrm{d}x)$-orthogonal projection on $X_h$. We prove the following result.

**Theorem 3.8.** For any $t \ge 0$, we consider $g(t) \in H^r(\mathrm{d}\mu_t)$, with $r \ge k + 1$. Then we have

(i) the global projection error for all $t \in [0, T]$,

$$ \|g(t) - \mathcal{P}_{W_\delta} g(t)\|_{L^2(\mathrm{d}\mu_t)} \le C \left( \frac{1}{(\alpha^2(t)N)^{r/2}} + h^{k+1} \right) \|g(t)\|_{H^r(\mathrm{d}\mu_t)}, $$

(ii) the global projection error on the fluxes for all $t \in [0, T]$,

$$ \left( \int_{\mathbb{R}} \|(g - \mathcal{P}_{W_\delta} g)(t, \cdot, v)\|_{L^2(\Gamma)}^2 \omega(t) \mathrm{d}v \right)^{1/2} \le C \left( \frac{1}{(\alpha^2(t)N)^{r/2}} + h^{k+1/2} \right) \|g(t)\|_{H^r(\mathrm{d}\mu_t)}. $$

*Proof.* To prove (i), we first write $g(t) \in L^2(\mathrm{d}\mu_t)$ as

$$ (3.18) \qquad g(t, x, v) = \sum_{n \in \mathbb{N}} \hat{C}_n(t, x) \Psi_n(t, v), $$

where the modes $(\hat{C}_n)_{n \in \mathbb{N}}$ are computed from $g(t)$ using the orthogonality property (2.3). Hence, we get

$$ (3.19) \qquad g(t) - \mathcal{P}_{W_\delta} g(t) = \sum_{n \ge N} \hat{C}_n \Psi_n + \sum_{n=0}^{N-1} (\hat{C}_n - \mathcal{P}_{X_h} \hat{C}_n) \Psi_n. $$

On the one hand, observing that the first term of the right hand side is nothing else than $g(t) - \mathcal{P}_{V_N} g(t)$, it can be estimated thanks to Proposition 3.6

$$ \| \sum_{n \ge N} \hat{C}_n \Psi_n \|_{L^2(\mathrm{d}\mu_t)} \le \frac{C}{(\alpha^2(t)N)^{r/2}} \|g\|_{H^r(\mathrm{d}\mu_t)}. $$

On the other hand, the second term can be controlled by applying Proposition 3.7 with $g = \hat{C}_n$, for $0 \le n \le N - 1$, which yields

$$ \begin{aligned} \| \sum_{n=0}^{N-1} (\hat{C}_n - \mathcal{P}_{X_h} \hat{C}_n) \Psi_n \|_{L^2(\mathrm{d}\mu_t)}^2 &\le \bar{C} h^{2(k+1)} \alpha(t) \sum_{n=0}^{N-1} \| \hat{C}_n \|_{H^{k+1}(\mathrm{d}x)}^2, \\ &\le \bar{C} h^{2(k+1)} \| g \|_{H^{k+1}(\mathrm{d}\mu_t)}^2. \end{aligned} $$

Gathering the latter results and using that $r \ge k + 1$, we prove that there exists a constant $C > 0$, independent of the discretization parameter $\delta = (h, 1/N)$, such that for all $t \in [0, T]$

$$ \|g(t) - \mathcal{P}_{W_\delta} g(t)\|_{L^2(\mathrm{d}\mu_t)} \le C \left( \frac{1}{(\alpha^2(t)N)^{r/2}} + h^{k+1} \right) \|g(t)\|_{H^r(\mathrm{d}\mu_t)}. $$

The second estimate in (ii) is obtained using the same ideas, together with the definition (3.15) of the norm over the mesh skeleton and the classical trace theorem on Sobolev spaces. Indeed, from (3.19) and Proposition 3.6, we obtain that for all $x \in \mathbb{T}$,

$$ \begin{aligned} \int_{\mathbb{R}} |g(t, x, v) - \mathcal{P}_{W_\delta} g(t, x, v)|^2 \omega(t) \mathrm{d}v &\le \frac{C}{(\alpha^2(t)N)^r} \|g(t, x, \cdot)\|_{H^r(\omega(t)\mathrm{d}v)}^2 \\ &\quad + \alpha(t) \sum_{n=0}^{N-1} |\hat{C}_n - \mathcal{P}_{X_h} \hat{C}_n|^2(t, x). \end{aligned} $$

Furthermore, since $H^r(\mathbb{T}) \subset L^\infty(\mathbb{T})$ for $r \ge 1$ and thanks to Proposition 3.7 on the mesh skeleton $\Gamma$, it yields

$$ \int_{\mathbb{R}} \|(g - \mathcal{P}_{W_\delta} g)(t, \cdot, v)\|_{L^2(\Gamma)}^2 \omega(t) \mathrm{d}v \le C \left( \frac{1}{(\alpha^2(t)N)^r} + h^{2k+1} \right) \|g(t)\|_{H^r(\mathrm{d}\mu_t)}^2, $$

from which we deduce the second item.

□
4. PROOF OF THEOREM 2.4

From the previous stability analysis and global projection error, we are now ready to prove our main result on the convergence of the numerical solution $f_\delta$ given by (2.9)-(2.12) to the solution $f$ to the Vlasov-Poisson system (1.1). By the triangle inequality, we have

$$ (4.1) \qquad \|f(t) - f_\delta(t)\|_{L^2(\mathrm{d}\mu_t)} \le \|f(t) - P_{W_\delta}f(t)\|_{L^2(\mathrm{d}\mu_t)} + \|P_{W_\delta}f(t) - f_\delta(t)\|_{L^2(\mathrm{d}\mu_t)}, $$

where the first term on the right hand side is the projection error, which has already been estimated in Theorem 3.8, whereas the second one is the consistency error, which will be treated by considering its time derivative and using stability arguments together with interpolation properties.

We define the consistency error $\mathcal{B}$ for a given $C = (C_n)_{n \in \mathbb{N}}$ and a smooth test function $\varphi$ as

$$ (4.2) \qquad \mathcal{B}_n(C, E, \varphi) := \int_{\mathbb{T}} \partial_t C_n \varphi \, \mathrm{d}x + \sum_{j \in \mathcal{J}} \mathcal{A}_{n,j}(g_n(C), \varphi) - \int_{\mathbb{T}} \mathcal{S}_n[C, E] \varphi \, \mathrm{d}x, $$

where $(\mathcal{A}_{n,j}, g_n(C))$ are given in (2.10) and $\mathcal{S}_n[C, E]$ in (2.6).

On the one hand, from the modes $C_\delta = (C_{\delta,n})_{0 \le n < N}$ corresponding to $f_\delta$ satisfying (2.9), we construct $\bar{C}_\delta = (\bar{C}_{\delta,n})_{n \in \mathbb{N}}$ as

$$ \bar{C}_{\delta,n} = \begin{cases} C_{\delta,n}, & \text{if } 0 \le n \le N-1, \\ 0, & \text{else,} \end{cases} $$

which satisfies for any $n \in \mathbb{N}$

$$ (4.3) \qquad \mathcal{B}_n(\bar{C}_\delta, E_\delta, \varphi_n) = 0, \quad \forall \varphi_n \in X_h, $$

where $E_\delta$ is solution to (2.13).

On the other hand, by consistency of the numerical flux (2.11), the first modes $\hat{C} = (\hat{C}_n)_{n \in \mathbb{N}}$ corresponding to the exact continuous solution $f$ of (1.1) satisfies for all $n \in \mathbb{N}$

$$ (4.4) \qquad \mathcal{B}_n(\hat{C}, E, \varphi_n) = 0, \quad \forall \varphi_n \in X_h, $$

with $E$ given by the second equation (Poisson equation) of (1.1). Then we introduce, for any $t \ge 0$, $\eta_\delta(t) \in W_\delta$ as

$$ (4.5) \qquad \eta_\delta(t, x, v) := P_{W_\delta}f(t, x, v) - f_\delta(t, x, v) = \sum_{n \in \mathbb{N}} \eta_{\delta,n}(t, x) \Psi_n(t, v), $$

with

$$ \eta_{\delta,n} = P_{X_h} \hat{C}_n \mathbf{1}_{\{n<N\}} - \bar{C}_{\delta,n}, \quad n \in \mathbb{N}, $$

and $\xi_\delta(t) \in L^2(\mathrm{d}\mu_t)$ as

$$ (4.6) \qquad \xi_\delta(t, x, v) := f(t, x, v) - P_{W_\delta}f(t, x, v) = \sum_{n \in \mathbb{N}} \xi_{\delta,n}(t, x) \Psi_n(t, v), $$

with

$$ \xi_{\delta,n}(t, x) := \hat{C}_n - P_{X_h} \hat{C}_n \mathbf{1}_{\{n<N\}}, \quad n \in \mathbb{N}. $$

Since $\eta_{\delta,n} \in X_h$, by taking $\varphi_n = \eta_{\delta,n}$ in (4.3) and (4.4) and substrating the two equalities, we get

$$ (4.7) \qquad \mathcal{B}_n(\eta_\delta, E_\delta, \eta_{\delta,n}) = -\mathcal{B}_n(\xi_\delta, E_\delta, \eta_{\delta,n}) + \alpha \sqrt{n} \int_{\mathbb{T}} (E - E_\delta) \hat{C}_{n-1} \eta_{\delta,n} \, \mathrm{d}x. $$

Now, the aim is to estimate the consistency error defined for any $t \ge 0$ by

$$ \|P_{W_\delta}f(t) - f_\delta(t)\|_{L^2(\mathrm{d}\mu_t)}^2 = \sum_{n=0}^{N-1} \alpha(t) \int_{\mathbb{T}} |\eta_{\delta,n}|^2 \mathrm{d}x. $$

To do this, we compute the time derivative of $\|P_{W_\delta}f(t) - f_\delta(t)\|_{L^2(\mathrm{d}\mu_t)}^2$ given by

$$ \frac{1}{2} \frac{\mathrm{d}}{\mathrm{d}t} \left[ \sum_{n=0}^{N-1} \alpha(t) \int_{\mathbb{T}} |\eta_{\delta,n}(t)|^2 \mathrm{d}x \right] = \sum_{n=0}^{N-1} \int_{\mathbb{T}} \left[ \alpha(t) \eta_{\delta,n} \partial_t \eta_{\delta,n} + \frac{1}{2} \alpha'(t) |\eta_{\delta,n}|^2 \right] \mathrm{d}x, $$
hence using (4.7) together with the definition (4.2) of $\mathcal{B}_n$, we obtain

$$ (4.8) \qquad \frac{1}{2} \frac{d}{dt} \left[ \sum_{n=0}^{N-1} \alpha(t) \int_{\mathbb{T}} |\eta_{\delta,n}(t)|^2 \, dx \right] = \mathcal{I}_1 + \mathcal{I}_2 + \mathcal{I}_3, $$

where

$$ (4.9) \qquad \left\{ \begin{array}{l} \mathcal{I}_1 := - \sum_{n=0}^{N-1} \alpha(t) \left( \sum_{j \in \mathcal{J}} \mathcal{A}_{n,j}(g_n(\eta_\delta), \eta_{\delta,n}) - \int_{\mathbb{T}} \mathcal{S}_n[\eta_\delta, E_\delta] \eta_{\delta,n} \, dx - \frac{1}{2} \frac{\alpha'(t)}{\alpha(t)} \int_{\mathbb{T}} |\eta_{\delta,n}|^2 \, dx \right), \\ \mathcal{I}_2 := \sum_{n=0}^{N-1} \alpha^2(t) \sqrt{n} \int_{\mathbb{T}} (E - E_\delta) \hat{C}_{n-1} \eta_{\delta,n} \, dx, \\ \mathcal{I}_3 := - \sum_{n=0}^{N-1} \alpha(t) \mathcal{B}_n(\xi_\delta, E_\delta, \eta_{\delta,n}). \end{array} \right. $$

Let us now estimate each of these terms separately. Throughout the following computations, $\mathcal{C}$ will be a generic positive constant, depending on the $L^2(\mathrm{d}\mu_t)$ of the exact solution $f$ and its derivatives, but independent of $\delta = (h, 1/N)$, and which may be different in different places.

We proceed as in the stability analysis detailed in [5, Proposition 3.2] or in Proposition 2.1, replacing $C_\delta$ by $\eta_\delta$. We get

$$ (4.10) \qquad \mathcal{I}_1 \le -\frac{1}{2} \sum_{n=0}^{N-1} \sum_{j \in \mathcal{J}} \nu_n [\eta_{\delta,n}]_{j-1/2}^2 + \frac{1}{4\gamma} \|\eta_\delta(t)\|_{L^2(\mathrm{d}\mu_t)}^2. $$

For the second term $\mathcal{I}_2$, we apply the Cauchy-Schwarz inequality and write

$$ |\mathcal{I}_2| \le \|E - E_\delta\|_{L^\infty(\mathrm{d}x)} \left( \sum_{n=0}^{N-1} \alpha^3(t) n \int_{\mathbb{T}} |\hat{C}_{n-1}|^2 \mathrm{d}x \right)^{1/2} \left( \sum_{n=0}^{N-1} \alpha(t) \int_{\mathbb{T}} |\eta_{\delta,n}|^2 \mathrm{d}x \right)^{1/2}. $$

On the one hand, to estimate $\|E - E_\delta\|_{L^\infty(\mathrm{d}x)}$, we proceed as in the proof of [5, Proposition 2.3]. By the Sobolev and Poincaré-Wirtinger inequalities, there exists a constant $\mathcal{C} > 0$ such that

$$ \|E - E_\delta\|_{L^\infty}^2 \le \mathcal{C} \|\partial_x(E - E_\delta)\|_{L^2(\mathrm{d}x)}^2. $$

Substrating (2.13) and the second equation of (1.1), taking $\partial_x(E - E_\delta)$ as test function, and applying the Cauchy-Schwarz inequality, it yields

$$ \|\partial_x(E - E_\delta)\|_{L^2}^2 \le \|\partial_x(E - E_\delta)\|_{L^2} \|\hat{C}_0 - C_{\delta,0}\|_{L^2}, $$

and then

$$ \|E - E_\delta\|_{L^\infty}^2 \le \mathcal{C} \|\hat{C}_0 - C_{\delta,0}\|_{L^2}^2 \le \frac{\mathcal{C}}{\alpha(t)} \|f(t) - f_\delta(t)\|_{L^2(\mathrm{d}\mu_t)}^2. $$

On the other hand, we remark that using the decomposition (3.16) of $f$ and the third equality in (3.1), we have

$$ \partial_v f = - \sum_{n=0}^{+\infty} \alpha \sqrt{n} \hat{C}_{n-1} \Psi_n, $$

and then

$$ \sum_{n=0}^{N-1} \alpha^3(t) n \int_{\mathbb{T}} |\hat{C}_{n-1}|^2 \mathrm{d}x \le \sum_{n \in \mathbb{N}} \alpha(t) \int_{\mathbb{T}} |\alpha(t) \sqrt{n} \hat{C}_{n-1}|^2 \mathrm{d}x = \|\partial_v f\|_{L^2(\mathrm{d}\mu_t)}^2. $$

Gathering these results, we obtain an estimate on $\mathcal{I}_2$ as

$$ |\mathcal{I}_2| \le \frac{\mathcal{C}}{\sqrt{\alpha(t)}} \|f(t) - f_\delta(t)\|_{L^2(\mathrm{d}\mu_t)} \|\eta_\delta(t)\|_{L^2(\mathrm{d}\mu_t)}. $$
Since the total error is the sum of the projection and consistency errors, we get by the triangle inequality and the use of the first item of Theorem 3.8 that

$$|\mathcal{I}_2| \le \frac{C}{\sqrt{\alpha(t)}} \left[ \frac{1}{(\alpha^2 N)^{m/2}} + h^{k+1} + \|\eta_\delta\|_{L^2(\mathrm{d}\mu_t)} \right] \|\eta_\delta\|_{L^2(\mathrm{d}\mu_t)},$$

where the constant $C > 0$ depends on the weighted $H^m$ norm of $f(t)$. Finally applying the Young inequality, we conclude that

$$(4.11) \qquad |\mathcal{I}_2| \le \frac{C}{\sqrt{\alpha(t)}} \left( \frac{1}{(\alpha^2 N)^m} + h^{2(k+1)} + \|\eta_\delta\|_{L^2(\mathrm{d}\mu_t)}^2 \right).$$

Now we turn to the last term $\mathcal{I}_3$ and use the definition (4.2) of $B_n$ to obtain

$$\mathcal{I}_3 = \mathcal{I}_{31} + \mathcal{I}_{32} + \mathcal{I}_{33},$$

with

$$\left\{
\begin{aligned}
\mathcal{I}_{31} &:= - \sum_{n=0}^{N-1} \alpha(t) \int_{\mathbb{T}} \left( \partial_t \xi_{\delta,n} - \frac{\alpha'(t)}{\alpha(t)} \left( n \xi_{\delta,n} + \sqrt{(n-1)n} \xi_{\delta,n-2} \right) \right) \eta_{\delta,n} \mathrm{d}x, \\
\mathcal{I}_{32} &:= + \sum_{n=0}^{N-1} \alpha^2(t) \sqrt{n} \int_{\mathbb{T}} E_{\delta} \xi_{\delta,n-1} \eta_{\delta,n} \mathrm{d}x, \\
\mathcal{I}_{33} &:= - \sum_{n=0}^{N-1} \alpha(t) \sum_{j \in \mathcal{J}} \mathcal{A}_{n,j}(g_n(\xi_{\delta}), \eta_{\delta,n}).
\end{aligned}
\right.$$

We start with $\mathcal{I}_{31}$ and remark that

$$\sum_{n=0}^{N-1} \alpha(t) \int_{\mathbb{T}} \left[ \partial_t \xi_{\delta,n} - \frac{\alpha'(t)}{\alpha(t)} \left( n \xi_{\delta,n} + \sqrt{(n-1)n} \xi_{\delta,n-2} \right) \right]^2 \mathrm{d}x \le \|\partial_t f(t) - P_{W_\delta} \partial_t f(t)\|_{L^2(\mathrm{d}\mu_t)}^2.$$

Then, since $f(t) \in H^m(\mathrm{d}\mu_t)$ satisfies the Vlasov equation (1.1), and using the second estimate (3.6) of Lemma 3.4, we have

$$\begin{aligned}
\|\partial_t f\|_{H^{m-1}(\mathrm{d}\mu_t)} &\le \|v \partial_x f\|_{H^{m-1}(\mathrm{d}\mu_t)} + \|E\|_{L^\infty} \|\partial_v f\|_{H^{m-1}(\mathrm{d}\mu_t)} \\
&\le \frac{C}{\alpha^2(t)} \|f\|_{H^m(\mathrm{d}\mu_t)}.
\end{aligned}$$

Thus, applying the Cauchy-Schwarz inequality to $\mathcal{I}_{31}$ and using Theorem 3.8 to $\partial_t f$ with $r = m - 1$, we have

$$(4.12) \qquad |\mathcal{I}_{31}| \le \frac{C}{\alpha^2(t)} \left( \frac{1}{(\alpha^2(t)N)^{(m-1)/2}} + h^{k+1} \right) \|\eta_\delta\|_{L^2(\mathrm{d}\mu_t)}.$$

The estimate on $\mathcal{I}_{32}$ follows the same lines as the one for $\mathcal{I}_{31}$. Indeed, remarking again that

$$\sum_{n=0}^{N-1} \alpha(t) \int_{\mathbb{T}} |\alpha(t) \sqrt{n} \xi_{\delta,n-1}|^2 \mathrm{d}x \le \|\partial_v f - P_{W_\delta} \partial_v f\|_{L^2(\mathrm{d}\mu_t)}^2,$$

and applying the Cauchy-Schwarz inequality to $\mathcal{I}_{32}$ and Theorem 3.8 to $\partial_v f$ with $r = m - 1$, we get

$$|\mathcal{I}_{32}| \le C \|E_\delta\|_{L^\infty} \left( \frac{1}{(\alpha^2(t) N)^{(m-1)/2}} + h^{k+1} \right) \left( \alpha(t) \sum_{n=0}^{N-1} \int_{\mathbb{T}} |\eta_{\delta,n}|^2 \mathrm{d}x \right)^{1/2},$$

which can be written as

$$(4.13) \qquad |\mathcal{I}_{32}| \le \frac{C}{\sqrt{\alpha(t)}} \left( \frac{1}{(\alpha^2(t) N)^{(m-1)/2}} + h^{k+1} \right) \|\eta_\delta(t)\|_{L^2(\mathrm{d}\mu_t)}.$$
Finally, we turn to the estimation of $\mathcal{I}_{33}$ and split it as $\mathcal{I}_{33} = \mathcal{I}_{331} + \mathcal{I}_{332}$, with

$$ \begin{cases} \mathcal{I}_{331} := \sum_{n=0}^{N-1} \alpha(t) \sum_{j \in \mathcal{J}} \int_{I_j} g_n(\xi_\delta) \partial_x \eta_{\delta,n} \, d x, \\ \mathcal{I}_{332} := - \sum_{n=0}^{N-1} \alpha(t) \sum_{j \in \mathcal{J}} \left( (\hat{g}_n(\xi_\delta) \eta_{\delta,n}^-)_{j+1/2} - (\hat{g}_n(\xi_\delta) \eta_{\delta,n}^+)_{j-1/2} \right). \end{cases} $$

Using the definition (2.10) of $g_n(\xi_\delta)$, we have

$$ \mathcal{I}_{331} = \sum_{n=0}^{N-1} \sum_{j \in \mathcal{J}} \int_{I_j} \left( \sqrt{n} \xi_{\delta,n-1} + \sqrt{n+1} \xi_{\delta,n+1} \right) \partial_x \eta_{\delta,n} \, d x . $$

Since $\xi_{\delta,n} = \hat{C}_n - \mathcal{P}_{X_h} \hat{C}_n$ for $n = 0, \dots, N-1$ and $\partial_x \eta_{\delta,n} \in X_h$, we have by definition of the projection $\mathcal{P}_{X_h}$ that $\mathcal{I}_{331} = 0$. Now, it remains to estimate $\mathcal{I}_{332}$. Using the periodic boundary conditions, we may write it as

$$ \mathcal{I}_{332} = \sum_{n=0}^{N-1} \alpha(t) \sum_{j \in \mathcal{J}} (\hat{g}_n(\xi_\delta) [\eta_{\delta,n}])_{j-1/2} . $$

Then, applying the Young inequality, we obtain

$$ |\mathcal{I}_{332}| \le \frac{\alpha(t)}{2} \sum_{n=0}^{N-1} \sum_{j \in \mathcal{J}} \left( \frac{\alpha(t)}{\nu_n} |\hat{g}_n(\xi_\delta)_{j-1/2}|^2 + \frac{\nu_n}{\alpha(t)} [\eta_{\delta,n}]^2_{j-1/2} \right) . $$

The second term of the right hand side will be balanced with the dissipation term figuring in the estimate (4.10) of $\mathcal{I}_1$, whereas the first term is estimated as follows. Using the definition of the numerical flux $\hat{g}_n$ (2.11) and the artificial viscosity $\nu_n$, we have

$$ \frac{\alpha^2(t)}{2} \sum_{n=0}^{N-1} \sum_{j \in \mathcal{J}} \frac{1}{\nu_n} |\hat{g}(\xi_\delta)_{n,j-1/2}|^2 \le \sum_{n=0}^{N-1} \sum_{j \in \mathcal{J}} \frac{1}{\nu} \{\sqrt{n} \xi_{\delta,n-1} + \sqrt{n+1} \xi_{\delta,n+1}\}_{j-1/2}^2 \\ + \bar{\nu} \sum_{n=0}^{N-1} \sum_{j \in \mathcal{J}} [\xi_{\delta,n}]^2_{j-1/2} . $$

On the one hand, observing that

$$ v f - \mathcal{P}_{W_\delta}(v f) = \sum_{n \in \mathbb{N}} \frac{1}{\alpha} \left( \sqrt{n} \xi_{\delta,n-1} + \sqrt{n+1} \xi_{\delta,n+1} \right) \Psi_n, $$

we apply the second item of Theorem 3.8 to $g = v f$ with $r = m - 1$ and the second inequality (3.6) of Lemma 3.4,

$$ \sum_{n=0}^{N-1} \sum_{j \in \mathcal{J}} \frac{1}{\nu} \left\{ \sqrt{n} \xi_{\delta,n-1} + \sqrt{n+1} \xi_{\delta,n+1} \right\}_{j-1/2}^2 \le C \alpha^2 \int_{\mathbb{R}} \|v f - \mathcal{P}_{W_\delta}(v f)\|_{L^2(\Gamma)}^2 \omega(t) \, d v \\ \le C \left( \frac{1}{(\alpha^2(t)N)^{m-1}} + h^{2k+1} \right) \|f(t)\|_{H^m(\mathrm{d} \mu_t)}^2 . $$

The viscosity term is estimated in the same manner, hence we deduce that

$$ (4.14) \quad |\mathcal{I}_{33}| \le C \left( \frac{1}{(\alpha^2(t)N)^{m-1}} + h^{2k+1} \right) \|f(t)\|_{H^m(\mathrm{d} \mu_t)}^2 + \frac{1}{2} \sum_{n=0}^{N-1} \sum_{j \in \mathcal{J}} \nu_n [\eta_{n,h}]^2_{j-1/2} . $$
Gathering (4.12), (4.13) and (4.14), we obtain

$$ |\mathcal{I}_3| \le \frac{C}{\alpha^2(t)} \left( \frac{1}{(\alpha^2(t) N)^{(m-1)/2}} + h^{k+1} \right) \|\eta_\delta(t)\|_{L^2(d\mu_t)} \\ + C \left( \frac{1}{(\alpha^2(t) N)^{m-1}} + h^{2k+1} \right) + \frac{1}{2} \sum_{n=0}^{N-1} \sum_{j \in \mathcal{J}} \nu_n [\eta_{\delta,n}(t)]_{j-1/2}^2. $$

Finally applying the Young inequality to the first term, we conclude that

$$ (4.15) \quad |\mathcal{I}_3| \le \frac{C}{\alpha^4(t)} \left( \frac{1}{(\alpha^2(t) N)^{m-1}} + h^{2k+1} \right) + \|\eta_\delta(t)\|_{L^2(d\mu_t)}^2 + \frac{1}{2} \sum_{n=0}^{N-1} \sum_{j \in \mathcal{J}} \nu_n [\eta_{n,h}(t)]_{j-1/2}^2. $$

Collecting (4.10), (4.11) and (4.15) in (4.8), it yields

$$ \frac{1}{2} \frac{d}{dt} \| \mathcal{P}_{W_\delta} f(t) - f_\delta(t) \|_{L^2(d\mu_t)}^2 \le \frac{C}{\sqrt{\alpha(t)}} \| \mathcal{P}_{W_\delta} f(t) - f_\delta(t) \|_{L^2(d\mu_t)}^2 \\ + \frac{C}{\alpha^4(t)} \left( \frac{1}{(\alpha^2(t) N)^{m-1}} + h^{2k+1} \right). $$

We conclude using the Gronwall lemma that for any $t \ge 0$, we have

$$ \| \mathcal{P}_{W_\delta} f(t) - f_\delta(t) \|_{L^2(d\mu_t)}^2 \le \int_0^t \frac{C}{\alpha^4(s)} \left( \frac{1}{(\alpha^2(s) N)^{m-1}} + h^{2k+1} \right) ds \exp \left( \int_0^t \frac{C}{\sqrt{\alpha(\tau)}} d\tau \right). $$

Using the projection estimate of Theorem 3.8 together with the latter inequality, and since by Proposition 2.3 the function $\alpha$ is bounded from below and above on a finite interval of time, we have established (2.19).

## 5. NUMERICAL SIMULATIONS

It is worth to mention that the Hermite/discontinuous Galerkin method has already been validated in [16, 5] on classical numerical tests. Hence in this section, we perform complementary numerical simulations on the Vlasov-Poisson system (1.1) using the DG/Hermite Spectral method to illustrate our theoretical result and to investigate the impact of the choice of the free parameter $\gamma \ge 0$ which enters in the definition of the scaling function $\alpha$ in (2.15). We also refer to [37] for a discussion on the latter point.

5.1. **Test 1: order of convergence.** We take $N$ modes for the Hermite spectral bases and $N_x$ cells in space, and apply a third order Runge-Kutta scheme for the time discretization with a small time step $\Delta t = 0.001$ in order to neglect the time discretization error. The initial scaling parameter $\alpha(0)$ is chosen to be 1 and the Hou-Li filter with 2/3 dealiasing rule [20, 10] will be used.

We choose the following initial condition

$$ (5.1) \qquad f_0(x, v) = (1 + \delta \cos(kx)) \frac{1}{\sqrt{2\pi}} \exp\left(-\frac{v^2}{2}\right), $$

with $\delta = 0.01$, which corresponds to the Landau damping configuration. The background density is $\rho_0 = 1$, the length of the domain in the $x$-direction is $L = 4\pi$ (that is $k = \pi/6$) and the final time is $T = 0.1$. The free parameter $\gamma$ is chosen as $\gamma = 1$ and the errors are computed by comparing to a reference solution obtained using $N_x \times N = 512 \times 512$ with $P_2$ piecewise polynomial basis and Hermite polynomial in velocity. In Table 5.1, we show the weighted $L^2$ errors and orders for $P_k$ piecewise polynomials with $k = 1, 2$ respectively. Due to the fact that the time steps are smaller than the spatial mesh size, we can observe $(k+1)$-th order of convergence for $P_k$ polynomials respectively.
<table><caption>TABLE 5.1. Test 1: Numerical weighted $L^2$ errors and orders for Landau damping with initial distribution (5.1), $\delta = 0.01$ and $k = \pi/6$, $T = 0.5$.</caption><thead><tr><th rowspan="2">N<sub>x</sub> &times; N</th><th colspan="2">P<sub>1</sub></th><th colspan="2">P<sub>2</sub></th></tr><tr><th>L<sup>2</sup>(d &mu;<sub>t</sub>) error</th><th>Order</th><th>L<sup>2</sup>(d &mu;<sub>t</sub>) error</th><th>Order</th></tr></thead><tbody><tr><td>16 &times; 16</td><td>5.12E-4</td><td>-</td><td>1.44E-5</td><td>-</td></tr><tr><td>32 &times; 32</td><td>1.05E-4</td><td>2.28</td><td>1.68E-6</td><td>3.09</td></tr><tr><td>64 &times; 64</td><td>2.31E-5</td><td>2.18</td><td>2.05E-7</td><td>3.04</td></tr><tr><td>128 &times; 128</td><td>5.42E-6</td><td>2.09</td><td>2.48E-8</td><td>3.04</td></tr></tbody></table>

5.2. **Test 2: Bump-on-the-tail.** Now we investigate the impact of the parameter $\gamma$ appearing in the definition of scaling function $\alpha$ in (2.15). According to our analysis, when $\gamma$ decreases, the function $\alpha$ decays more slowly and the weighted norm $\|f_{\delta}(t)\|_{L^2(d\mu_t)}$ is bounded as

$$ (5.2) \qquad \|f_{\delta}(t)\|_{L^2(d\mu_t)} \le \|f_{\delta}(0)\|_{L^2(d\mu_0)} e^{t/4\gamma}. $$

For practical computation, we expect that the scaling function $\alpha$ follows the variation of the distribution function in the velocity space, hence it is crucial to have a good understanding of the impact of this free parameter. For these reasons we present a numerical example on the bump-on-the-tail [5, Section 5.2] where the distribution function strongly varies in $v$. We consider the initial distribution as

$$ (5.3) \qquad f(0, x, v) = f_b(v)(1 + \kappa \cos(k n x)), $$

where the bump-on-tail distribution is

$$ (5.4) \qquad f_b(v) = \frac{n_p}{\sqrt{\pi} v_p} e^{-v^2/v_p^2} + \frac{n_b}{\sqrt{\pi} v_b} e^{-(v-v_d)^2/v_b^2}. $$

We choose a strong perturbation with $\kappa = 0.04$, $n = 3$ and $k = 2\pi/L$ with $L = 62$ and the other parameters are set to be $n_p = 0.9$, $n_b = 0.1$, $v_d = 4.5$, $v_p = \sqrt{2}$, $v_b = \sqrt{2}/2$. The computational domain is $[0, L] \times [-8, 8]$. These settings have been used in [33] and [16, Section 4.3]. For this case, we take the initial scaling function to be $\alpha_0 = 5/7$. We take $N_x \times N = 64 \times 128$.

For the Vlasov-Poisson system, we know that the total energy and the $L^2$ norm of $f$ are exactly conserved, but this property is no more true at the discrete level. In Figure 5.1 (a)-(b), we present the time evolution of these quantities for different values of $\gamma$ corresponding to the definition (2.15). On the one hand, the amplitude of the variations of the total energy are of order $10^{-6}$ for the different values of $\gamma$ and are even smaller when $\gamma$ is small. On the other hand, the $L^2$ norm of $f$ oscillates around its initial value, but again the impact of the parameter $\gamma$ is negligible (see (b)). We also present the time evolution of the potential and kinetic energy in Figure 5.1 (c)-(d) for different values of $\gamma$. From these plots, we can observe that the impact of this free parameter is limited and does not affect the accuracy of the method.

Finally, in Figure 5.2, we present the time evolution of $\alpha$ and the corresponding weighted $L^2$ norm for different values of $\gamma$. Since the initial data $\alpha(0)$ is the same for our simulations, we know that when $\gamma \le \gamma'$, we have $\alpha_{\gamma'}(t) \le \alpha_{\gamma}(t)$ (since $\alpha_{\gamma'}$ decays faster than $\alpha_{\gamma}$) and then

$$ \|f(t)\|_{L^2(d\mu'_t)} \le \|f(t)\|_{L^2(d\mu_t)}, $$

where $\mu'_t$ represents the measure associated to $\alpha_{\gamma'}(t)$. We also notice that when $\gamma$ increases, the weighted $L^2$ norm $\|f(t)\|_{L^2(d\mu_t)}$ grows slowly, which is consistent with our estimate (5.2). Moreover, for a fixed $\gamma$, the time evolution of $\|f(t)\|_{L^2(d\mu_t)}$ first follows the growth of the potential energy which is almost exponentially fast, but when nonlinear effects dominate, it starts to oscillate and is stabilized. This last numerical result illustrates that the estimate (5.2) is certainly not optimal for large time...
<figure><img src="image_2.png" /><figcaption>FIGURE 5.1. **Test 2:** time evolution of (a) the variation of the total energy, (b) the standard $L^2$ norm of $f$, (c) the potential energy and (d) the kinetic energy for $\alpha$ given by (2.15), with several values of $\gamma = 0, \dots, 10^{-1}$.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>FIGURE 5.2. **Test 2:** time evolution of (a) the scaling function $\alpha$ and (b) the corresponding weighted $L^2$ norm of $f$, for $\alpha$ given by (2.15), with several values of $\gamma = 0, \dots, 10^{-1}$.</figcaption></figure>
## 6. CONCLUSION & PERSPECTIVES

In this article we investigate the convergence analysis of a spectral Hermite discretization of the Vlasov-Poisson system with a time-dependent scaling factor allowing to prove the stability and convergence of the numerical solution in an appropriate functional framework. The control of this scaling factor, and more precisely a positive lower bound, is crucial to ensure completely the stability of the method and our proof follows carefully this dependency. Our analysis is limited to the one dimensional case for which the control of the electric field is straightforward. In a future work we would like to adapt our approach to the multi-dimensional case following the ideas in [3] for the control of the $L^\infty$ norm of the electric field in two and three dimensions.

## ACKNOWLEDGEMENT

Both authors want to thank the two referees for their motivating comments to improve our analysis. Marianne Bessemoulin-Chatard is partially funded by the Centre Henri Lebesgue (ANR-11-LABX-0020-01) and ANR Project MoHyCon (ANR-17-CE40-0027-01). Francis Filbet is partially funded by the ANR Project Muffin (ANR-19-CE46-0004) and by the EUROfusion Consortium and has received funding from the Euratom research and training programme 2019-2020 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

## REFERENCES

[1] T. P. Armstrong. Numerical studies of the nonlinear Vlasov equation. *The Physics of Fluids*, 10(6):1269-1280, 1967.

[2] B. Ayuso, J. A. Carrillo, and C.-W. Shu. Discontinuous Galerkin methods for the one-dimensional Vlasov-Poisson system. *Kinet. Relat. Models*, 4(4):955-989, 2011.

[3] B. Ayuso, J. A. Carrillo, and C.-W. Shu. Discontinuous Galerkin methods for the multi-dimensional Vlasov-Poisson problem. *Mathematical Models and Methods in Applied Sciences*, 22(12):1250042, 2012.

[4] N. Besse. Convergence of a high-order semi-Lagrangian scheme with propagation of gradients for the one-dimensional Vlasov-Poisson system. *SIAM J. Numer. Anal.*, 46(2):639-670, 2008.

[5] M. Bessemoulin-Chatard and F. Filbet. On the stability of conservative discontinuous Galerkin/Hermite spectral methods for the Vlasov-Poisson system. *Journal of Computational Physics*, 451:110881, 2022.

[6] C. K. Birdsall and A. B. Langdon. *Plasma physics via computer simulation*. McGraw-Hill, New York, 1985.

[7] E. Camporeale, G. Delzanno, B. Bergen, and J. Moulton. On the velocity space discretization for the Vlasov-Poisson system: Comparison between implicit Hermite spectral and Particle-in-Cell methods. *Computer Physics Communications*, 198:47-58, 2016.

[8] F. Charles, B. Després, and M. Mehrenberger. Enhanced convergence estimates for semi-Lagrangian schemes: application to the Vlasov-Poisson equation. *SIAM J. Numer. Anal.*, 51(2):840-863, 2013.

[9] Y. Cheng, I. M. Gamba, and P. J. Morrison. Study of conservation and recurrence of Runge-Kutta discontinuous Galerkin schemes for Vlasov-Poisson systems. *Journal of Scientific Computing*, 56(2):319-349, 2013.

[10] Y. Di, Y. Fan, Z. Kou, R. Li, and Y. Wang. Filtered Hyperbolic Moment Method for the Vlasov Equation. *Journal of Scientific Computing*, 79(2):969-991, 2019.

[11] R. Duclous, B. Dubroca, F. Filbet, and V. Tikhonchuk. High order resolution of the Maxwell-Fokker-Planck-Landau model intended for ICF applications. *Journal of Computational Physics*, 228(14):5072-5100, 2009.

[12] F. Engelmann, M. Feix, E. Minardi, and J. Oxenius. Nonlinear effects from Vlasov's equation. *The Physics of Fluids*, 6(2):266-275, 1963.

[13] F. Filbet. Convergence of a finite volume scheme for the Vlasov-Poisson system. *SIAM J. Numer. Anal.*, 39(4):1146-1169, 2001.

[14] F. Filbet and E. Sonnendrücker. Comparison of Eulerian Vlasov solvers. *Comput. Phys. Comm.*, 150(3):247-266, 2003.

[15] F. Filbet and E. Sonnendrücker. Modeling and numerical simulation of space charge dominated beams in the paraxial approximation. *Mathematical Models and Methods in Applied Sciences*, 16(05):763-791, 2006.

[16] F. Filbet and T. Xiong. Conservative Discontinuous Galerkin/Hermite Spectral Method for the Vlasov-Poisson System. *Commun. Appl. Math. Comput.*, 2020.

[17] J. C. M. Fok, B. Guo, and T. Tang. Combined Hermite spectral-finite difference method for the Fokker-Planck equation. *Mathematics of computation*, 71(240):1497-1528, 2001.

[18] R. E. Heath, I. M. Gamba, P. J. Morrison, and C. Michler. A discontinuous Galerkin method for the Vlasov-Poisson system. *Journal of Computational Physics*, 231(4):1140-1174, 2012.

[19] J. P. Holloway. Spectral velocity discretizations for the Vlasov-Maxwell equations. *Transport theory and statistical physics*, 25(1):1-32, 1996.
[20] T. Y. Hou and R. Li. Computing nearly singular solutions using pseudo-spectral methods. *Journal of Computational Physics*, 226:379-397, 2007.

[21] C. Johnson and J. Pitkäranta. An analysis of the discontinuous Galerkin method for a scalar hyperbolic equation. *Mathematics of computation*, 46(173):1–26, 1986.

[22] G. Joyce, G. Knorr, and H. K. Meier. Numerical integration methods of the Vlasov equation. *Journal of Computational Physics*, 8(1):53–63, 1971.

[23] A. J. Klimas and W. M. Farrell. A splitting algorithm for Vlasov simulation with filamentation filtration. *J. Comput. Phys.*, 110(1):150–163, 1994.

[24] K. Kormann and A. Yurova. A generalized Fourier-Hermite method for the Vlasov-Poisson system. *BIT Numerical Mathematics*, pages 1–29, 2021.

[25] S. Le Bourdiec, F. De Vuyst, and L. Jacquet. Numerical solution of the Vlasov-Poisson system using generalized Hermite functions. *Computer physics communications*, 175(8):528–544, 2006.

[26] P.-L. Lions and B. Perthame. Propagation of moments and regularity for the 3-dimensional Vlasov-Poisson system. *Invent. Math.*, 105(2):415–430, 1991.

[27] H. Ma, W. Sun, and T. Tang. Hermite spectral methods with a time-dependent scaling for parabolic equations in unbounded domains. *SIAM journal on numerical analysis*, 43(1):58–75, 2005.

[28] H. Ma and T. Zhao. A stabilized Hermite spectral method for second-order differential equations in unbounded domains. *Numerical Methods for Partial Differential Equations: An International Journal*, 23(5):968–983, 2007.

[29] G. Manzini, G. L. Delzanno, J. Vencels, and S. Markidis. A Legendre-Fourier spectral method with exact conservation laws for the Vlasov-Poisson system. *Journal of Computational Physics*, 317:82–107, 2016.

[30] G. Manzini, D. Funaro, and G. L. Delzanno. Convergence of Spectral Discretizations of the Vlasov-Poisson System. *SIAM Journal on Numerical Analysis*, 55(5):2312–2335, 2017.

[31] K. Pfaffelmoser. Global classical solutions of the Vlasov-Poisson system in three dimensions for general initial data. *J. Differential Equations*, 95(2):281–303, 1992.

[32] J. W. Schumer and J. P. Holloway. Vlasov simulations using velocity-scaled Hermite representations. *Journal of Computational Physics*, 144(2):626–661, 1998.

[33] M. Shoucri and G. Knorr. Numerical integration of the Vlasov equation. *Journal of Computational Physics*, 14(1):84–92, 1974.

[34] J. W. Schumer and J. P. Holloway. Vlasov simulations using velocity-scaled Hermite representations. *Journal of Computational Physics*, 144(2):626–661, 1998.

[35] E. Sonnendrücker, F. Filbet, A. Friedman, E. Oudet, and J.-L. Vay. Vlasov simulations of beams with a moving grid. *Computer Physics Communications*, 164(1-3):390–395, 2004.

[36] E. Sonnendrücker, J. Roche, P. Bertrand, and A. Ghizzo. The semi-Lagrangian method for the numerical resolution of the Vlasov equation. *Journal of Computational Physics*, 149(2):201–220, 1999.

[37] T. Tang. The Hermite spectral method for Gaussian-type functions. *SIAM journal on scientific computing*, 14(3):594–606, 1993.

[38] S. Ukai and T. Okabe. On classical solutions in the large in time of two-dimensional Vlasov's equation. *Osaka Math. J.*, 15(2):245–261, 1978.

[39] C. Yang and M. Mehrenberger. Highly accurate monotonicity-preserving semi-Lagrangian scheme for Vlasov-Poisson simulations. *J. Comput. Phys.*, 446:Paper No. 110632, 33, 2021.

[40] S. Zaki, L. Gardner, and T. Boyd. A finite element code for the simulation of one-dimensional Vlasov plasmas. I. Theory. *Journal of Computational Physics*, 79(1):184–199, 1988.

[41] S. I. Zaki, T. Boyd, and L. Gardner. A finite element code for the simulation of one-dimensional Vlasov plasmas. II. Applications. *Journal of Computational Physics*, 79(1):200–208, 1988.

Email address: francis.filbet@math.univ-toulouse.fr

Email address: marianne.bessemoulin@univ-nantes.fr