Alain Blaustein

Francis Filbet
email: francis.filbet@math.univ-toulouse.fr

A STRUCTURE AND ASYMPTOTIC PRESERVING SCHEME FOR THE VLASOV-POISSON-FOKKER-PLANCK MODEL

Keywords: 2010 Mathematics Subject Classification. Primary: 82C40, Secondary: 65N08, 65N35 Hermite spectral method, Vlasov-Poisson-Fokker-Planck, Hypocoercive estimates

We propose a numerical method for the Vlasov-Poisson-Fokker-Planck model written as an hyperbolic system thanks to a spectral decomposition in the basis of Hermite functions with respect to the velocity variable and a structure preserving finite volume scheme for the space variable. On the one hand, we show that this scheme naturally preserves both stationary solutions and linearized free-energy estimate. On the other hand, we adapt previous arguments based on hypocoercivity methods to get quantitative estimates ensuring the exponential relaxation to equilibrium of the discrete solution for the linearized Vlasov-Poisson-Fokker-Planck system, uniformly with respect to both scaling and discretization parameters. Finally, we perform substantial numerical simulations for the nonlinear system to illustrate the efficiency of this approach for a large variety of collisional regimes (plasma echos for weakly collisional regimes and trend to equilibrium for collisional plasmas) and to highlight its robustness (unconditional stability, asymptotic preserving properties).

Introduction

The Vlasov-Poisson-Fokker-Planck system provides a kinetic description of a gas constituted of charged particles, let us say electrons and heavy positive ions, interacting through a mean electrostatic field:

(1.1)

       ∂f ∂t + v • ∇ x f + q m E • ∇ v f = 1 τ div v (vf + T 0 ∇ v f ) , E = -∇ x Φ ; -ε 0 ∆ x Φ = q (n -n i ) ; n = R 3 f dv .
In (1.1), f (t, x, v) is the distribution of electrons over the phase space (x, v) ∈ T 3 × R 3 at time t ≥ 0. Field interactions are taken into account thanks to a coupling between kinetic and Poisson equations (first and second line of (1.1) respectively). The coupling displays several constants: vacuum permittivity ε 0 , elementary charge q, mass m of an electron as well as macroscopic densities n(t, x) of electrons and n i (x) of ions. Thermodynamic effects are taken into account thanks to a Fokker-Planck operator, on the right-hand side of the kinetic equation, where appears the relaxation time τ > 0 of electrons due to their collisions with the ionic surrounding bath described by a spatially homogeneous temperature T 0 > 0.

To investigate the physical behavior of the solution to system (1.1), two important quantities will be considered [START_REF] Filbet | Fokker-Planck multi-species equations in the adiabatic asymptotics[END_REF]. The first one is ε the square root of the ratio between the mass of electrons and positive heavy ions, given by

ε := m e m i 1,
whereas the second one is τ (ε) > 0, the ratio between the elapsed time between two collisions of electrons and the observable time. We focus on a situation where the one dimensional Vlasov-Poisson-Fokker-Planck system can be reformulated using these parameters as

(1.2)        ε ∂ t f + v ∂ x f -∂ x φ ∂ v f = ε τ (ε) ∂ v (v f + T 0 ∂ v f ) , -∂ 2 x φ = ρ -ρ i , ρ(t, x) = R f (t, x, v) dv .
This system is completed with the following condition, which ensures uniqueness of the electrical potential φ

(1.3) T φ(t, x) dx = 0 .
Since we focus on the situation where the electron to ion mass ratio is very small, it allows to describe ions as a steady thermal bath. More precisely, ions are supposed to be fixed with a prescribed temperature T 0 > 0 and a density ρ i > 0, which is an integrable function over T. Furthermore, the following quasi-neutrality assumption is satisfied for all time t ≥ 0 T ρ(t, x) dx = T ρ i (x) dx as soon as this condition is initially verified. As already mentioned, the other scaling parameter τ (ε) > 0 stands for the ratio between the time which separates two collisions of an electron with the ionic background and the timescale of observation. In this work, we suppose τ (ε) = τ 0 ε.

Therefore, as ε goes to zero, we expect that the couple (f, φ), solution to (1.2)-(1.3), converges to (f ∞ , φ ∞ ) given by

       v ∂ x f ∞ -∂ x φ ∞ ∂ v f ∞ = 1 τ 0 ∂ v (v f ∞ + T 0 ∂ v f ∞ ) , -∂ 2 x φ ∞ = ρ ∞ -ρ i , ρ ∞ (x) = R f ∞ (x, v) dv .
Actually, the equilibrium state f ∞ is uniquely determined as follows

f ∞ (x, v) = ρ ∞ (x) M(v) ,
where M denotes the Maxwellian with temperature T 0

(1.4) M(v) = 1 √ 2π T 0 exp - |v| 2 2 T 0 ,
and where ρ ∞ solves (1.5)

     ρ ∞ = c ∞ exp - φ ∞ T 0 , -∂ 2 x φ ∞ = ρ ∞ -ρ i , completed with the following condition on the constant c ∞ > 0 c ∞ T exp - φ ∞ T 0 dx = T ρ i dx .
Therefore, the electrons' temperature relaxes to the background temperature whereas their spatial distribution converges to a Maxwell-Boltzmann equilibrium.

Let us point out that there is an attractive version of the Vlasov-Poisson-Fokker-Planck system (1.2)- (1.3) which is also widely used in stellar physics. For that case the repulsive electrostatic force is replaced by the attractive gravitational force, responsible for a change of sign in the Poisson equation. In this paper, we only consider the repulsive case. The numerical resolution of the system (1.2)-(1.3) shares the same difficulties as most of kinetic equations: high dimensionality, presence of various scales, etc... Several numerical methods have been developed for (1.2)-(1.3), we mention for instance [START_REF] Havlak | On deterministic particle methods for solving Vlasov-Poisson-Fokker-Planck systems[END_REF][START_REF] Wollman | Numerical approximation of the Vlasov-Poisson-Fokker-Planck system in one dimension[END_REF][START_REF] Crestetto | Kinetic/fluid micro-macro numerical schemes for Vlasov-Poisson-BGK equation using particles[END_REF][START_REF] Degond | Asymptotic-preserving methods and multiscale models for plasma physics[END_REF][START_REF] Carrillo | Variational Asymptotic Preserving Scheme for the Vlasov-Poisson-Fokker-Planck System[END_REF][START_REF] Jang | High order asymptotic preserving DG-IMEX schemes for discrete-velocity kinetic equations in a diffusive scaling[END_REF]. These numerical schemes are either deterministic or stochastic, with an effort to capture some physical phenomena associated to weakly collisional plasmas such as Landau damping or two-stream instability, occurring for short time range, before being canceled by collisions. More recently, dynamical low-rank algorithms have been proposed [START_REF] Einkemmer | An asymptotic-preserving dynamical low-rank method for the multi-scale multidimensional linear transport equation[END_REF][START_REF] Coughlin | Efficient dynamical low-rank approximation for the Vlasov-Ampère-Fokker-Planck system[END_REF], they decouple the dimensions of the phase space allowing to reduce the computational cost.

Here, we want to design a numerical scheme preserving the large-time behavior of solutions to the Vlasov-Poisson-Fokker-Planck equation or related nonlinear kinetic models. More precisely, we want to describe accurately the transient regime in which transport effects dominate and where the distribution function exhibits filamentation, but also large time scales in which collisions take over, forcing the relaxation of the distribution to thermal equilibrium, given here by a Maxwell-Boltzmann distribution. Let us remind that when we neglect both the coupling with the Poisson equation and the free transport term, a fully discrete finite difference scheme for the homogeneous Fokker-Planck equation has been proposed in the pioneering work of Chang and Cooper [11]. This scheme has nice properties since it preserves the stationary solution and the entropy decay of the numerical solution. Later, finite volume schemes preserving the exponential trend to equilibrium have been studied for nonlinear convection-diffusion equations (see for example [START_REF] Bessemoulin-Chatard | A finite volume scheme for nonlinear degenerate parabolic equations[END_REF][START_REF] Burger | A mixed finite element method for nonlinear diffusion equations[END_REF][START_REF] Chainais-Hillairet | Asymptotic behaviour of a finite-volume scheme for the transient drift-diffusion model[END_REF] and more recently [START_REF] Georgoulis | Hypocoercivity-compatible finite element methods for the long-time computation of kolmogorov's equation[END_REF][START_REF] Dujardin | Coercivity, hypocoercivity, exponential time decay and simulations for discrete Fokker-Planck equations[END_REF][START_REF] Bessemoulin-Chatard | Hypocoercivity and diffusion limit of a finite volume scheme for linear kinetic equations[END_REF] in the frame of hypocoercive methods). Let us also mention that in [START_REF] Pareschi | Residual equilibrium schemes for time dependent partial differential equations[END_REF], authors investigate the question of describing correctly the equilibrium state of nonlinear diffusion and kinetic models for high order schemes.

More recently, in our previous work [START_REF] Blaustein | On a discrete framework of hypocoercivity for kinetic equations[END_REF], we have proposed and carefully studied a numerical scheme based on Hermite polynomials in the velocity space for the linear Vlasov-Fokker-Planck equation when the electric field E ∞ = -∂ x φ ∞ is prescribed. This approach preserves the equilibrium and thanks to a discrete hypocoercive method, it is shown that numeric solutions relax exponentially fast to thermodynamic equilibrium. When collisions are neglected, Hermite decomposition has also been successfully applied to approximate the Vlasov-Poisson system on a finite time interval [START_REF] Filbet | Conservative Discontinuous Galerkin/Hermite Spectral Method for the Vlasov-Poisson System[END_REF][START_REF] Bessemoulin-Chatard | On the stability of conservative discontinuous Galerkin/Hermite spectral methods for the Vlasov-Poisson system[END_REF][START_REF] Bessemoulin-Chatard | On the convergence of discontinuous Galerkin/Hermite spectral methods for the Vlasov-Poisson system[END_REF]. Therefore, our goal is to provide here an efficient approximation of the nonlinear Vlasov-Poisson-Fokker-Planck system (1.2)-(1.3) able to describe accurately both weakly and strongly collisional regimes and to preserve the long time behavior of the solution. To illustrate the feature of our numerical scheme, we will analyze the asymptotic behavior of the solution given by the discrete linearized system and then propose various numerical experiments in various regimes.

In order to design a well-balanced approximation for (1.2)-(1.3), we consider an equivalent reformulation where we define a potential ψ = φ -φ ∞ and replace the electrical potential φ in the equation, it yields that ρ i is now replaced by the quantities at equilibrium as follows,

(1.6)        ε ∂ t f + v ∂ x f -∂ x φ ∞ ∂ v f -∂ x ψ ∂ v f = 1 τ 0 ∂ v (v f + T 0 ∂ v f ) , -∂ 2 x ψ = ρ -ρ ∞ , ρ(t, x) = R f (t, x, v) dv ,
coupled with the condition on ψ,

(1.7) T ψ(t, x) dx = 0 . The equilibrium to (1.6)-(1.7) is now characterized by (f ∞ , ψ ∞ ) where ψ ∞ ≡ 0.
The key-estimate to prove the trend to equilibrium of solutions to (1.6)-(1.7) is given by

d dt H(f, f ∞ ) = - 1 ε τ 0 I(f, f ∞ ),
where H(f, f ∞ ) denotes the free energy

H(f, f ∞ ) := T×R ln f f ∞ f dx dv + 1 2 T 0 ∂ x ψ 2 L 2 (T) ,
and

I(f, f ∞ ) is the entropy dissipation I(f, f ∞ ) := 4 T 0 T×R ∂ v f f ∞ 2 f ∞ dx dv .
The free energy estimate is said to be "hypocoercive" since we have

I(f (t), f ∞ ) = 0 ⇐⇒ f (t, x, v) = ρ(t, x) M(v) , ∀(x, v) ∈ T × R ,
meaning that the entropy dissipation I controls at most the distance between f and its associated local equilibrium ρ M, giving no straightforward information on the long time behavior of ρ.

In the linear case, corresponding to ψ ≡ 0 in (1.6), L. Desvillettes and C. Villani have proposed a general method, based on the latter relative entropy estimate and logarithmic Sobolev inequalities, to overcome this degeneracy in the position variable [START_REF] Desvillettes | On the trend to global equilibrium in spatially inhomogeneous entropy-dissipating systems: the linear Fokker-Planck equation[END_REF]. This approach has been widely explored in the last decade [START_REF] Villani | Hypocoercivity[END_REF][START_REF] Dolbeault | Hypocoercivity for linear kinetic equations conserving mass[END_REF]. However, when the Vlasov-Fokker-Planck equation is coupled with the Poisson equation for the electrical potential, a different functional framework, based on weighted L 2 spaces, is applied motivated by the following estimate when f is near equilibrium [START_REF] Hérau | On global existence and trend to the equilibrium for the Vlasov-Poisson-Fokker-Planck system with exterior confining potential[END_REF][START_REF] Herda | Large-time behavior of solutions to Vlasov-Poisson-Fokker-Planck equations: from evanescent collisions to diffusive limit[END_REF]. Indeed, plugging the following formal expansion into the free energy

f ln f f ∞ ∼ f →f∞ f -f ∞ + |f -f ∞ | 2
2 f ∞ and using that mass is conserved for solutions to (1.6)-(1.7), we define a new functional, named the linearized free energy, as

(1.8) E(t) = f (t) -f ∞ 2 L 2 (f -1 ∞ ) + 1 T 0 ∂ x ψ(t) 2 L 2 (T) .
Unfortunately, this functional is not dissipated for the solution to the nonlinear system (1.6)-(1.7), but only for its linearized version given by (1.9)

       ε ∂ t f + v ∂ x f -∂ x φ ∞ ∂ v f -∂ x ψ ∂ v f ∞ = 1 τ 0 ∂ v (v f + T 0 ∂ v f ) , -∂ 2 x ψ = ρ -ρ ∞ , ρ(t, x) = R f (t, x, v) dv ,
coupled with the condition on ψ,

(1.10) T ψ(t, x) dx = 0 .
This yields for the solution (f, ψ) to (1.9)-(1.10) (see Proposition 2.1 for a complete proof)

(1.11) 1 2 d dt E(t) = - T 0 ε τ 0 T×R ∂ v f (t) f ∞ 2 f ∞ dx dv .
The purpose of this paper is to design a numerical scheme for the nonlinear Vlasov-Poisson-Fokker-Planck system (1.6)-(1.7) for which such estimate occurs on its linearized version (1.9)-(1.10). To this aim, we propose a simple time splitting scheme, where the first stage consists in solving the linearized system (1.9)-(1.10), whereas the second stage solves the remaining quadratic part of (1.6)-(1.7), that is

ε ∂ t f -∂ x ψ ∂ v (f -f ∞ ) = 0 , for which ψ is unchanged.
This approach has several advantages from the computational and stability point of view. Indeed, both steps will be fully implicit in time allowing to use a large time step, uniformly with respect to the parameter ε. Moreover, the linearized equation is autonomous, hence it requires to solve the same linear system at each time step, which can be done using a direct solver. Furthermore, solving the time dependent implicit second step is in fact negligible in terms of computational costs since the associated system is trivially invertible due to the Hermite discretization. Finally, the numerical approximation of the linearized system allows to capture a consistent asymptotic profile when ε → 0 and it also preserves the free energy estimate as in [START_REF] Blaustein | On a discrete framework of hypocoercivity for kinetic equations[END_REF], which treats the case of a Vlasov-Fokker-Planck equation without a coupling with Poisson.

In Section 2 we propose a numerical discretization of the full model (1.6)-(1.7) based on Hermite's decomposition in the velocity space and finite volume scheme for the space discretization. Then, in Section 3, we prove quantitative properties on its linearized version (1.9)-(1.10). More precisely, we first prove a discrete version of the free energy estimate (1.11) and then, using discrete hypocoercive methods, we prove the exponential trend to equilibrium with rate 1/ε, this uniformly with respect to discretization parameters. This result constitutes a theoretical proof of the asymptotic-preserving properties of the method at the linearized level. Finally, in Section 4, we carry out numerical experiments which illustrate the robustness of our scheme in a wide variety of situations ranging from near-to-collisionless regime 1 τ 0 to the stiff limit ε ∼ 0 and including inhomogeneous ionic background. In particular, we highlight the asymptotic preserving properties of the method. Furthermore, we observe formation of nonlinear echoes and study their suppression in weakly collisional settings as well as simultaneous vortex/filamentation formation for inhomogeneous ionic background. These phenomena have drawn intense mathematical interest in the kinetic community over the past decade [START_REF] Bedrossian | Suppression of Plasma Echoes and Landau Damping in Sobolev Spaces by Weak Collisions in a Vlasov-Fokker-Planck Equation[END_REF][START_REF] Grenier | Plasma Echoes Near Stable Penrose Data[END_REF][START_REF] Chaturvedi | The Vlasov-Poisson-Landau system in the weakly collisional regime[END_REF]. With this work, we aim at taking part in these efforts by proposing numerical methods capable to capture these phenomena efficiently.

Numerical scheme

In this section, we will introduce our numerical method. We follow the lines of our previous work for the linear Vlasov-Fokker-Planck equation [START_REF] Blaustein | On a discrete framework of hypocoercivity for kinetic equations[END_REF] and then propose a discretization of the Poisson equation allowing to preserve the energy estimate for the linearized problem (1.9)-(1.10). Finally, the nonlinear part is discretized using an implicit scheme.

2.1. Hermite's decomposition for the velocity variable. Let us first focus on the discretization of the velocity variable. It consists in performing a spectral decomposition of the distribution f into its Hermite modes (Ψ k ) k∈N defined as

Ψ k (v) = H k v √ T 0 M(v) ,
and which constitute an orthonormal system for the inverse Gaussian weight since it holds

R Ψ k (v) Ψ l (v) M -1 (v)dv = δ k,l ,
where M is the Maxwellian corresponding to the stationary state of the Fokker-Planck operator (1.4) and δ k,l the Kronecker symbol (δ k,l = 1 when k = l and δ k,l = 0 otherwise). In the latter definition, (H k ) k∈N stands for the family of Hermite polynomials defined recursively as follows H -1 = 0, H 0 = 1 and

ξ H k (ξ) = √ k H k-1 (ξ) + √ k + 1 H k+1 (ξ) , ∀ k ≥ 0 .
Let us also point out that Hermite's polynomials verify the following relation

H k (ξ) = √ k H k-1 (ξ) , ∀ k ≥ 0 .
The Hermite system arises naturally in our context since it offers a simple discrete reformulation of the L 2 f -1 ∞ -norm which appears in the key estimate (1.11), indeed it holds

f (t) 2 L 2 (f -1 ∞ ) = k∈N C k (t) 2 L 2 (ρ -1 ∞ ) , where C = (C k ) k∈N stand for the Hermite components of f f (t, x, v) = k∈N C k (t, x) Ψ k (v) .
As one can see in the latter relation, each term of the sequence C = (C k ) k∈N naturally belongs to the weighted space L 2 ρ -1 ∞ . From the numerical point of view, working in weighted spaces induces difficulties when it comes to integro/differential manipulation such as integration by part. This is the reason why rather than discretizing coefficients C = (C k ) k∈N , we consider their re-normalized versions D = (D k ) k∈N defined as

(2.1) f (t, x, v) = √ ρ ∞ (x) k∈N D k (t, x) Ψ k (v) .
According to latter considerations, renormalized Hermite coefficients D = (D k ) k∈N verify

f (t) 2 L 2 (f -1 ∞ ) = k∈N D k (t) 2 L 2 (T) .
To sum up, renormalized Hermite coefficients play a fundamental role in our analysis for two reasons: they offer a discrete reformulation of the key quantity E(t) given by (1.8) and they belong to the unweighted L 2 -Lebesgue space over T. Moreover, there is another benefit coming out of this choice: thanks to the properties of Hermite polynomials, one can see that Hermite functions diagonalize the Fokker-Planck operator since it holds

∂ v [ v Ψ k + T 0 ∂ v Ψ k ] = -k Ψ k .
Therefore, following [START_REF] Blaustein | On a discrete framework of hypocoercivity for kinetic equations[END_REF], we substitute the decomposition (2.1) in the first line of (1.6) : using the identities

E ∞ = -∂ x φ ∞ and ρ ∞ E ∞ = T 0 ∂ x ρ ∞ , we get that D = (D k ) k∈N satisfies the following system (2.2)      ε ∂ t D k + √ k A D k-1 - √ k + 1 A D k+1 + k T 0 ∂ x ψ D k-1 = - k τ 0 D k , ∀ k ∈ N , D k (t = 0) = D 0,ε k , ∀ k ∈ N ,
where we set D -1 = 0. Due to our unweighted L 2 -framework, the latter formulation enjoys a nice duality structure since A and A are adjoint operators in L 2 (T). More precisely, A and A are given by

       A u = + T 0 ∂ x u - E ∞ 2 √ T 0 u , A u = -T 0 ∂ x u - E ∞ 2 √ T 0 u .
Notice that both operators A and A may also be rewritten as follows

           A u = + T 0 ρ ∞ ∂ x u √ ρ ∞ , A u = - T 0 ρ ∞ ∂ x √ ρ ∞ u .
To conclude, we denote by D ∞ the Hermite decomposition of the equilibrium f ∞ . It is determined for all k ∈ N by

(2.3) D ∞,k = √ ρ ∞ , if k = 0 , 0, else .

Poisson equation formulated in the Hermite framework.

To compute the electrical potential ψ, we will reformulate the Poisson equation in such a way that the free energy estimate (1.11) for the linearized system (1.9)-(1.10) is satisfied. The main idea to preserve such estimate is to construct a scheme for which it is possible to perform at the discrete level analogous calculations as the ones needed at the continuous level to derive the free energy estimate (1.11). To this aim, we introduce a modified potential ω given by

ω = √ ρ ∞ T 0 ψ .
Using the definition of the operator A, the electric field

E = -∂ x ψ is now given by √ ρ ∞ E = - √ ρ ∞ ∂ x ψ = -T 0 A ω
and the modified potential ω solves the modified Poisson equation

(2.4) A ρ -1 ∞ A ω = D 0 - √ ρ ∞ .
This new formulation will allow us to easily construct a numerical scheme for the Poisson equation preserving the key energy estimate (1.11), where in this framework, the linearized system (1.9)-(1.10) rewrites

(2.5)

             ∂ t D k + 1 ε √ k A D k-1 - √ k + 1 A D k+1 + A ω δ k,1 = - k ε τ 0 D k , ∀ k ∈ N , A ρ -1 ∞ A ω = D 0 - √ ρ ∞ , D k (t = 0) = D 0,ε k , ∀ k ∈ N , completed with the condition (2.6) T ω √ ρ ∞ dx = 0 ,
and where the linearized free energy E reads (2.7)

E(t) = 1 2 D(t) -D ∞ 2 L 2 + A ω √ ρ ∞ 2 L 2 (T)
, where • L 2 stands for the overall L 2 -norm with no weight

D 2 L 2 = k∈N D k 2 L 2 (T) .
From this reformulated equation, we first prove the free energy estimate on the linearized system (1.9)-(1.10).

Proposition 2.1. Consider a formal solution (D k ) k∈N to (2.5) such that for all t > 0,

k≥0 k D k (t) 2 L 2 < ∞.
The following free energy estimate holds for all t ≥ 0

(2.8) d dt E(t) + 1 ε τ 0 k∈N k D k (t) 2 L 2 (T) = 0 .
Proof. We multiply equation (2.5) by D k -D ∞,k , sum over all k ∈ N, integrate in x ∈ T and then rearrange all the terms, this yields

1 2 d dt D -D ∞ 2 L 2 + 1 ε T A ω D 1 dx = - 1 ε τ 0 k∈N * k D k 2 L 2 (T) .
We rewrite the integral term in the latter estimate using the duality structure and the equation on

D 0 T A ω D 1 dx = T ω A D 1 dx = ε T ω ∂ t (D 0 - √ ρ ∞ ) dx.
Using the reformulated Poisson equation in (2.5) and the duality structure again, we deduce

ε T ω ∂ t (D 0 - √ ρ ∞ ) dx = ε T ω ∂ t A ρ -1 ∞ A ω dx = ε 2 d dt A ω √ ρ ∞ 2 L 2 (T)
.

It finally yields the free energy estimate

d dt E(t) = - 1 ε τ 0 k∈N k D k (t) 2 L 2 (T) .
Let us now draw two conclusions from these considerations. On the one hand, we observe that the estimate given in Proposition 2.1 is not sufficient to prove convergence of the solution to the linearized system (2.5) towards the stationary state (2.3) because of the lack of coercivity. To bypass this difficulty, we will define a modified relative energy H as (2.9)

H(t) = E(t) + β 0 A D 1 (t), u(t) L 2 (T) ,
where β 0 > 0 is a small free parameter and u is solution to

     A A u = D 0 - √ ρ ∞ , T u √ ρ ∞ dx = 0 .
To get the convergence to the solution to the linearized system (2.5) to the stationary state, the strategy consists in proving that H and E are equivalent and that there exists a constant κ > 0 such that

d dt H(t) ≤ - κ ε min τ 0 , τ -1 0 H(t).
We do not detail these results in the continuous case since they constitute the object of Section 3 in the discrete framework.

On the other hand, from the reformulated Poisson equation for the linearized system (1.6)-(1.7), we are now able to write the Hermite method for the full Vlasov-Poisson-Fokker-Planck system, which reads

               ε ∂ t D k + √ k A D k-1 - √ k + 1 A D k+1 + k ρ ∞ A ω D k-1 = - k τ 0 D k , ∀ k ∈ N , A ρ -1 ∞ A ω = D 0 - √ ρ ∞ , D k (t = 0) = D 0,ε k , ∀ k ∈ N , completed with the condition T ω √ ρ ∞ dx = 0 .
This latter formulation is equivalent to the initial system (2.2). The Hermite spectral method now consists in considering a finite number of modes N H . Therefore we solve for k = 0, . . . , N H , (2.10)

               ε ∂ t D k + √ k A D k-1 - √ k + 1 A D k+1 + k ρ ∞ A ω D k-1 = - k τ 0 D k , A ρ -1 ∞ A ω = D 0 - √ ρ ∞ , D k (t = 0) = D 0,ε k , where we set D -1 = D N H +1 = 0.
It is worth to mention that Proposition 2.1 also holds true for (2.10). Indeed, this truncation does not modify the above computations on the free energy estimate.

In the following subsection, we design a well-balanced finite volume discretization of the nonlinear system (2.10) such that the associated approximation of the linearized system (2.5) satisfies the estimate given in Proposition 2.1 (Section 3).

2.3.

Finite volume discretization for the space variable. We now turn to the phase space and time discretization of the system (2.10) and propose a well-balanced time splitting scheme. To discretize the phase space domain, we fix a number of Hermite modes N H ∈ N * . Then, we consider an interval (a, b) of R and for N x ∈ N , introduce the set J = {1, . . . , N x } and a family of control volumes (K j ) j∈J such that K j = x j-1/2 , x j+1/2 with x j the middle of the interval K j and

a = x 1/2 < x 1 < x 3/2 < ... < x j-1/2 < x j < x j+1/2 < ... < x Nx < x Nx+1/2 = b . Let us set ∆x j = x j+1/2 -x j-1/2 , for j ∈ J , ∆x j+1/2 = x j+1 -x j , for 1 ≤ j ≤ N x -1 .
We also introduce the parameter h such that h = max j∈J ∆x j .

To discretize the time variable, we fix a time step ∆t and we set t n = n∆t with n ∈ N. Our time discretization of R + is then given by the increasing sequence of (t n ) n∈N .

(2.11)

                             ε D n+1/2 1 -D n 1 ∆t + A h D n+1/2 0 - √ 2 A h D n+1/2 2 + A h ω n+1/2 h = - 1 τ 0 D n+1/2 1 , ε D n+1/2 k -D n k ∆t + √ k A h D n+1/2 k-1 - √ k + 1 A h D n+1/2 k+1 = - k τ 0 D n+1/2 k , for k = 1 , (A h ρ -1 ∞ A h ) ω n+1/2 h = D n+1/2 0 - √ ρ ∞ , j∈J ∆x j ω n+1/2 j √ ρ -1 ∞,j = 0 , for k ∈ {0, . . . , N H } and D n+1/2 k = 0 when k > N H and k = -1. Moreover, operator A h (resp. A h
) is an approximation of the operator A (resp. A ) given by (2.12)

A h = (A j ) j∈J and A h = (A j ) j∈J
and where for D = (D j ) j∈J it holds (2.13)

       A j D = T 0 D j+1 -D j-1 2∆x j - E ∞,j 2 √ T 0 D j , j ∈ J , A j D = -T 0 D j+1 -D j-1 2∆x j - E ∞,j 2 √ T 0 D j , j ∈ J ,
whereas the discrete electric field E ∞,j is given by

(2.14) E ∞,j = - φ ∞,j+1 -φ ∞,j-1 2∆x j = 2 T 0 √ ρ ∞,j √ ρ ∞,j+1 - √ ρ ∞,j-1 2 ∆x j ,
where ρ ∞,j is an approximation of the stationary density ρ ∞ on the cell K j . This latter formula is consistent with the definition of √ ρ ∞ = exp (-φ ∞ /(2T 0 )) and the fact that

1 2 T 0 ∂ x φ ∞ = - 1 √ ρ ∞ ∂ x √ ρ ∞ ,
furthermore, our discretization of the field E ∞ allows to preserve the equilibrium since it ensures

(2.15) A j √ ρ ∞ = T 0 √ ρ ∞,j+1 - √ ρ ∞,j-1 2∆x j - E ∞,j 2 √ T 0 √ ρ ∞,j = 0 , ∀j ∈ J .
This first step requires the numerical resolution of a linear system which does not depend on the time index n. Hence a direct solver based on LU factorization is applied to get the solution D n+1/2 , ω n+1/2 h .

On the other hand, we solve the nonlinear part using again a fully implicit Euler scheme (2.16) do not change during this second step, the latter system is trivially invertible and hence does not require any linear solver.

                 D n+1 0 = D n+1/2 0 , ω n+1 h = ω n+1/2 h , ε D n+1 k -D n+1/2 k ∆t + k ρ ∞ A h ω n+1 h D n+1 k-1 -D ∞,k-1 = 0 , if k ≥ 1 ,
In the next section we analyze the linearized step (2.11) for which we prove exponential relaxation towards equilibrium in the long time regime. Then, in Section 4, we will perform numerical simulations on the full nonlinear scheme (2.11)-(2.16).

Trend to equilibrium of the discrete linearized system

In this section, we only consider the numerical scheme applied to the linearized system (2.5) corresponding to the first step in the splitting method (2.11)-(2.16), that is, 2)

(3.1)                            ε D n+1 1 -D n 1 ∆t + A h D n+1 0 - √ 2 A h D n+1 2 + A h ω n+1 h = - 1 τ 0 D n+1 1 , ε D n+1 k -D n k ∆t + √ k A h D n+1 k-1 - √ k + 1 A h D n+1 k+1 = - k τ 0 D n+1 k , for k = 1 , (A h ρ -1 ∞ A h ) ω n+1 h = D n+1 0 - √ ρ ∞ , j∈J ∆x j ω n+1 j √ ρ -1 ∞,j =
E n = 1 2 D n -D ∞ 2 l 2 + A h ω n h √ ρ ∞ 2 l 2 (T)
, where

D 2 l 2 = N H k=0 D k 2 l 2 (T) , and D k 2 l 2 (T) = j∈J |D k,j | 2 ∆x j ,
and where D ∞ is defined as

D ∞,k = √ ρ ∞ , if k = 0 , 0, else ,
and recall [7, Lemma 3.3] that A h and A h are adjoint operators in l 2 (T).

Then, we prove that the solution to the fully discrete system (3.1) converges exponentially fast to its discrete equilibrium, which is consistent with the continuous system. Theorem 3.1. Consider the solution (D n ) n∈N to (3.1). Then the following discrete energy estimate holds for all n ≥ 0

E n ≤ 3 1 + κ ε min τ 0 , τ -1 0 ∆t -n E 0 ,
where κ > 0 depends only on ρ i , T 0 and |b -a|.

To show this result, we couple a discrete version of the free energy estimate in Proposition 2.1 with hypocoercive estimates for the discrete version of the modified relative entropy functional defined in (2.9).

Before to give the proof, let us comment this result. On the one hand, we emphasize that the rate of convergence is uniform with respect to discretization parameters. On the other hand, the convergence rate is proportional to 1/ε, regardless of discretization parameters. This last property ensures that the scheme for the linearized model is asymptotic preserving in the long time regime. It also ensures that our method applied to the linearized model is unconditionally stable, regardless of both scaling and discretization parameters. To conclude, we emphasize that the explicit dependence with respect to τ 0 is coherent with the results for the continuous model [START_REF] Herda | Large-time behavior of solutions to Vlasov-Poisson-Fokker-Planck equations: from evanescent collisions to diffusive limit[END_REF].

3.1.

A priori estimates. In this section, we prove a discrete free energy estimate on (E n ) n∈N analogous to the one in Proposition 2.1. Proposition 3.2. Consider the solution (D n ) n∈N to (3.1). The following discrete energy estimate holds for all n ≥ 0

E n+1 -E n ∆t + ∆t R n h = - 1 ε τ 0 N H k=1 k D n+1 k 2 l 2 (T) ,
where R n h is the following positive remainder due to numeric dissipation

R n h = 1 2   D n+1 -D n ∆t 2 l 2 + A h ω n+1 h -ω n h ∆t √ ρ ∞ 2 l 2 (T)   .
Proof. To compute the variations of D n -D ∞ 2 l 2 between time step n and n + 1, we take the l 2 (T) scalar product between the first line in (3.1) and

D n+1 1 -D ∞,1 , the second and D n+1 k -D ∞,k and sum over all k ∈ {0, . . . , N H }, this yields 1 2 ∆t D n+1 -D ∞ 2 l 2 -D n -D ∞ 2 l 2 + D n+1 -D n 2 l 2 = I 1 + I 2 ,
where I 1 and I 2 are given by

(3.3)          I 1 := - 1 ε τ 0 N H k=0 k D n+1 k 2 l 2 (T) - 1 ε N H k=0 √ k A h D n+1 k-1 - √ k + 1 A h D n+1 k+1 , D n+1 k -D ∞,k l 2 (T)
,

I 2 := - 1 ε A h ω n+1 h , D n+1 1 l 2 (T)
, where I 2 stands for the contribution of the electric field and I 1 gathers all the other terms.

First, we rewrite I 1 using that A h and A h are adjoint in l 2 (T) and that A h D ∞,0 = 0 according to (2.15)

I 1 = - 1 ε τ 0 N H k=0 k D n+1 k 2 l 2 (T) - 1 ε N H k=0 √ k A h D n+1 k-1 , D n+1 k l 2 (T) - √ k + 1 D n+1 k+1 , A h D n+1 k l 2 (T) .
Hence, splitting and re-indexing the second sum in the latter relation, we see that it is in fact zero. Therefore, I 1 rewrites as follows

I 1 = - 1 ε τ 0 N H k=0 k D n+1 k 2 l 2 (T) .
Furthermore, considering the case k = 0 in the second equation of system (3.1), we deduce that I 2 rewrites as follows

I 2 = -ω n+1 h , D n+1 0 -D n 0 ∆t l 2 (T)
.

Then, we replace D n+1 0 -D n 0 in the latter relation using the third line in (3.1), it yields

I 2 = -ω n+1 h , (A h ρ -1 ∞ A h ) ω n+1 h -ω n h ∆t l 2 (T)
.

Using that A h and A h are adjoint in l 2 (T), we deduce the following relation

I 2 = - 1 2 ∆t   A h ω n+1 h √ ρ ∞ 2 l 2 (T) - A h ω n h √ ρ ∞ 2 l 2 (T) + A h ω n+1 h -ω n h √ ρ ∞ 2 l 2 (T)   ,
which concludes the proof.

3.2. Proof of Theorem 3.1. We are now ready to proceed to the proof of Theorem 3.1. To do so, we develop a discrete hypocoercive technique which consists in introducing the analog H n of the continuous modified entropy functional defined in (2.9). It reads as follows

(3.4) H n = E n + β 0 A h D n 1 , u n h l 2 (T)
, where β 0 is a positive constant and where u n h is solution the solution to

(3.5)        A h A h u n h = D n 0 - √ ρ ∞ , j∈J ∆x j u n j √ ρ ∞,j = 0 .
Let us first recall useful estimates on u n h [7, Lemma 3.5] Lemma 3.3. For each n ∈ N, the solution u n h to (3.5) satisfies

(3.6)          A 2 h u n h l 2 (T) + A h u n h l 2 (T) ≤ C D n 0 - √ ρ ∞ l 2 (T) , A h u n+1 h -u n h ∆t l 2 (T) ≤ min 1 ε D n+1 1 l 2 (T) , C D n+1 0 -D n 0 ∆t l 2 (T)
, for some constant C > 0 only depending on ρ i and T 0 .

Among other, the latter Lemma ensures that the modified relative entropy functional H n is in fact equivalent to the energy E n for β 0 small enough. Lemma 3.4. For all ρ i and T 0 > 0, there exists a positive constant β 0 such that for all β 0 ∈ (0, β 0 ), one has

(3.7) 1 4 E n ≤ H n ≤ 3 4 E n , ∀ n ∈ N,
where E n and H n are given by (3.2) and (3.4).

Proof. Since A h and A h are adjoint in l 2 (T), applying Cauchy-Schwarz inequality and the first line in (3.6), we obtain

A h D n 1 , u n h l 2 (T) ≤ C D n -D ∞ 2
l 2 , which allows to bound the additional term in the definition of H n and therefore conclude the proof.

Building on the latter lemmas, we now prove that H n verifies a dissipation relation Proposition 3.5. Consider the solution (D n ) n∈N to (3.1). The modified relative entropy functional defined by (3.4) verifies for all n ≥ 0

H n+1 -H n ∆t ≤ - κ ε min τ 0 , τ -1 0 H n+1 ,
for some positive constant κ depending only on ρ i and T 0 and |b -a|.

Proof. We first focus on the additional term A h D n 1 , u n h l 2 (T) in the definition of H n . We write its discrete time derivative as follows

1 ∆t A h D n+1 1 , u n+1 h l 2 (T) -A h D n 1 , u n h l 2 (T) = A h D n+1 1 -D n 1 ∆t , u n+1 h l 2 (T) + A h D n 1 , u n+1 h -u n h ∆t l 2 (T)
and replace the discrete time derivative of D n 1 in the latter right-hand side thanks to the first line in (3.1)

(3.8) 1 ∆t A h D n+1 1 , u n+1 h l 2 (T) -A h D n 1 , u n h l 2 (T) = J 1 + J 2 + J 3 , where J 1 , J 2 and J 3 are given by                      J 1 := - 1 ε A h A h D n+1 0 - √ ρ ∞ , u n+1 h l 2 (T) , J 2 := - 1 ε A h A h ω n+1 h , u n+1 h l 2 (T) , J 3 := + 1 ε √ 2 (A h ) 2 D n+1 2 - 1 τ 0 A h D n+1 1 , u n+1 h l 2 (T) + A h D n 1 , u n+1 h -u n h ∆t l 2 (T)
.

First, J 1 is the desired dissipation term: since A h and A h are adjoint in l 2 (T) and according to (3.5), it holds

J 1 = - 1 ε D n+1 0 - √ ρ ∞ , A h A h u n+1 h l 2 (T) = - 1 ε D n+1 0 - √ ρ ∞ 2 l 2 (T)
. Then, J 2 takes into account the contribution of the electric field. As it turns out, it is also a dissipation term: thanks to (3.1) (third line) and (3.5), we have

A h A h u n+1 h = (A h ρ -1 ∞ A h ) ω n+1 h
, which yields

J 2 = - 1 ε ω n+1 h , (A h ρ -1 ∞ A h ) ω n+1 h l 2 (T) = - 1 ε A h ω n+1 h √ ρ ∞ 2 l 2 (T)
.

Finally, J 3 gathers all terms without good sign. Since A h and A h are adjoint in l 2 (T), it rewrites

J 3 = √ 2 ε D n+1 2 , (A h ) 2 u n+1 h l 2 (T) - 1 τ 0 ε D n+1 1 , A h u n+1 h l 2 (T) + D n+1 1 , A h u n+1 h -u n h ∆t l 2 (T) + D n 1 -D n+1 1 , A h u n+1 h -u n h ∆t l 2 (T)
.

We estimate the first two terms in the latter right-hand side using Young inequality and the first line in (3.6). The last two terms are estimated applying Cauchy-Schwarz inequality and the second line in (3.6). Hence, we get

J 3 ≤ 1 ε C τ 0 η D n+1 0 -D ∞,0 2
l 2 (T) + C η ε D n+1 2 2 l 2 (T) + 1 τ 0 D n+1 1 2 l 2 (T) + C ∆t D n+1 -D n 2 l 2 .
for any positive η, for some positive constant C depending only on ρ i and T 0 and with τ 0 = min(1, τ 0 ). Taking the sum between (3.8) multiplied by β 0 and the estimate in Proposition 3.2 and replacing J 1 , J 2 and J 3 with the latter estimates, it yields

H n+1 -H n ∆t + ∆t (1 -C β 0 ) R n h ≤ - 1 ε τ 0   1 - C η τ 0 β 0 τ 0 N H k=1 D n+1 k 2 l 2 (T) + β 0 τ 0 1 - C τ 0 η D n+1 0 -D ∞,0 2
l 2 (T) + β 0 τ 0 A h ω n+1 h √ ρ ∞ 2 l 2 (T)   .
Choosing η = τ 0 /(2C) and

β 0 τ 0 = 2τ 2 0 /(τ 2 0 + 4C 2 ) it holds 1 - C η τ 0 β 0 τ 0 N H k=1 D n+1 k 2 l 2 (T) + β 0 τ 0 1 - C τ 0 η D n+1 0 -D ∞,0 2 l 2 (T) = τ 2 0 4 C 2 + τ 2 0 D n+1 -D ∞ 2 l 2 ,
hence the latter inequality becomes

H n+1 -H n ∆t + ∆t 1 -C 2τ 2 0 τ 0 (τ 2 0 + 4C 2 ) R n h ≤ - 1 ετ 0 2τ 2 0 4 C 2 + τ 2 0 E n+1 .
According to Lemma 3.4 and since β 0 τ 0 = 2τ 2 0 /(τ 2 0 + 4C 2 ), we may replace E n with H n for C ≥ 0 great enough in the latter estimate. Hence, after simple computations, we deduce the result

H n+1 -H n ∆t ≤ - κ ε min τ 0 , τ -1 0 H n+1 ,
for some κ > 0 depending only on ρ i and T 0 and |b -a|.

We now conclude the proof of Theorem 3.1. First, from Proposition 3.5, it is straightforward to obtain

H n ≤ 1 + κ ε min τ 0 , τ -1 0 ∆t -n H 0 .
Then, we apply Lemma 3.4 on each side of the latter inequality and obtain the result.

Numerical simulations

For numerical experiments, we apply a slight modification of the scheme (2.11)-(2.16) since a Strang splitting scheme with a second order implicit Runge-Kutta scheme is used to get second order accuracy in time.

In our simulations, we fix the the temperature T 0 to 1 and numerical parameters as follows: N x = 128, ∆t = 0.1 and we adapt the number of Hermite modes depending on the collisional regime. In Section 4.1 we highlight the asymptotic preserving properties of our scheme and investigate the behavior of solutions when ε 1 and τ 0 is fixed. Then, from Section 4.2 to Section 4.4, we will fix ε = 1 in (2.11)-(2.16) and investigate the robustness of the scheme in different collisional regimes ranging from weakly collisional regime when τ 0 1 to strongly collisional plasmas when τ 0 1.

4.1.

Asymptotic-preserving properties. In this first test, we illustrate the robustness and the asymptotic preserving property of the scheme in the limit ε → 0. To do so, we keep ∆t = 0.1 fixed and perform numerical simulations with ε = 10 -k , for k ranging from 0 to 6. Since ε only appears in front of the time derivative, it can be interpreted as a time scaling parameter. We emphasize that the extreme case ε = 10 -6 corresponds to taking time step of the order ∆t/ε = 10 5 : we expect that the solution approaches f ∞ in one time step! In this study, the numerical parameters are fixed and we only consider a few number of Hermite modes taking N H = 80. We choose the following spatially inhomogeneous equilibrium

ρ ∞ (x) = c ∞ exp (-φ ∞ (x)) , x ∈ (-L, L) ,
where the potential φ ∞ is given by

φ ∞ (x) = 0.2 sin(k x) , x ∈ (-L, L),
with k = π/L and L = 6 whereas the constant c ∞ ensures that

1 2L L -L ρ ∞ (x) dx = 1.
Thus, we take the initial distribution function as a perturbation of this steady state, that is,

f (t = 0, x, v) = 1 √ 2π (ρ ∞ (x) + δ cos(k x)) exp - v 2 2 , (x, v) ∈ (-L, L) × R,
where δ = 0.01 and consider the case τ 0 = 10 5 , which corresponds to a weakly collisional regime. On Figure 4.1, we represent the time evolution of the potential energy and the quantity

(4.1) L 2 (t) := f (t) -f ∞ L 2 (f -1 ∞ ) .
L 2 f -1 ∞ (in log scale).
As predicted by our analysis of the linearized model, we observe exponential relaxation towards equilibrium at a rate which is proportional to 1/ε. The scheme is uniformly stable with respect to ε and the solution converges to the discrete equilibrium f ∞ when ε → 0. The left chart in Figure 4.2 represents the time evolution of the potential energy in log-scale for t/ε ∈ [0, 100]. We plot the approximations corresponding to ε = 10 -k , for k ranging from 0 to 3 and, in each case, re-scale the time variable as (s ← t/ε) to compare solutions. We can see that the case ε = 10 -1 (which corresponds to taking a time step ∆t/ε = 1) fits very well with the approximation obtained when ε = 1. Both the first phase for t/ε ∈ [0, 20] corresponding to fast oscillations and steep descent and the second phase for t/ε ∈ [20, 100] corresponding to slower oscillations without damping are transcribed correctly in this case. When ε = 10 -2 , we do not expect the solution to be precise during the first phase t/ε ∈ [0, 20] since the corresponding time step ∆t/ε = 10, is greater than the time period of the oscillations. However, we see that in the second phase t/ε ∈ [20, 100], the approximation catches up and even seems to capture the oscillatory behavior of the solution, despite the fact that time step and oscillation period have the same order of magnitude [START_REF] Chainais-Hillairet | Asymptotic behaviour of a finite-volume scheme for the transient drift-diffusion model[END_REF]. To conclude, the case ε = 10 -3 corresponds to taking a time step ∆t/ε = 100. Therefore, the approximation at time t/ε = 100 is surprisingly accurate considering that it was calculated in only a single iteration !

The right chart in Figure 4.2 represents the time evolution of the potential energy in log-scale for t/ε ∈ [0, 10 5 ]. We plot the approximations corresponding to ε = 10 -k , for k ranging from 3 to 6 and once again, in each case, re-scale the time variable as (s ← t/ε) to compare solutions. We do not observe oscillations anymore since collisions effects take over transport phenomena at this time scale. Therefore, we compare the exponential decay rate of our approximations in order to validate at the nonlinear level our theoretical result, which holds for the linearized scheme. We observe that in all cases k = 3, . . . , 6, approximations present similar decay rate, even when k = 6, which is surprising considering the fact that in this case the approximation at time t/ε = 10 5 was obtained in one iteration only. All these results are very satisfying! Indeed, for the differents regimes, corresponding to the values of ε, our numerical schemes is able to describe correctly the different phases: an oscillatory behavior when ε ≥ 10 -3 and an exponential decay to equilibrium when ε ≤ 10 -3 .

To conclude this section, it is worth to mention that in this test only, we tuned the expert options of the Super LU library [START_REF] Li | An overview of SuperLU: Algorithms, implementation, and user interface[END_REF] used as a direct solver in our code. More precisely, for small values of ε, the system associated to (2.11) may be ill conditioned, hence we disabled the equilibration option and tuned the threshold used for a diagonal entry to be an acceptable pivot in the factorization. 4.2. Perturbation of non uniform density. For this second numerical test, we consider the same initial condition as in the preceding Section 4.1. However, we now fix ε to 1 and perform simulations with variable τ 0 . To enforce numerical convergence, we have chosen a large number of Hermite modes N H = 400 when the plasma is weakly collisional, that is when τ 0 10 2 , since filamentation may occur in phase space whereas N H = 50 is enough when collisions dominate.

On the one hand, we take τ 0 = 10 4 corresponding to the weakly collisional regime and compare two solutions, one is obtained using (3.1) corresponding to the linearized Vlasov-Poisson-Fokker-Planck system (1.9) and the second one is given by (2.11)-(2.16) corresponding to the nonlinear Vlasov-Poisson-Fokker-Planck system (1.2). Our results show that both solutions have the same behavior, which means that, for such a small perturbation, the linear regime governs the dynamics. To illustrate this observation, we report in Figure 4. In Figure 4.4, we show several snapshots of the difference between the distribution function f and its equilibrium f ∞ for t ∈ [START_REF] Bessemoulin-Chatard | On the stability of conservative discontinuous Galerkin/Hermite spectral methods for the Vlasov-Poisson system[END_REF]70]. As expected, thin filaments propagate in phase space for large velocities but surprisingly we also observe that a vortex is generated in the region where |v| ≤ 1. For large time, this vortex remains and continues to rotate around the point (x C , v C ) = (-3, 0). For such a regime, where collisions are almost negligible, the amplitude of the perturbation does not vanish even when t 70 and transport phenomena dominate. On the other hand, we study the influence of the collision frequency τ 0 and perform several numerical simulations for the nonlinear Vlasov-Poisson-Fokker-Planck system (1.2) with the same initial data for τ 0 = 10 k , with k = 0, . . . , 4 (see Figure 4.5). For a weakly collision regime, that is k ≥ 3, we again observe oscillations of the potential energy and an exponential decay. However, when collisions dominate, fast oscillations only occur for short time, then the potential energy decays rapidly to zero without any oscillations. This trend to equilibrium can be also viewed on the distribution function as shown on Figure 4.5-(b), where we present the time evolution of the quantity L 2 defined in (4.1).

Finally in Figure 4.6, we again present snapshots of f -f ∞ at different time when τ 0 = 10 2 . In this situation, collisions dissipate thin filaments generated by the transport term and the amplitude of the perturbation vanishes. As a consequence, we do not detect anymore the vortex structure on the difference f -f ∞ and the solution converges exponentially fast to the equilibrium as it has been shown in Theorem 3.1 for the linearized system. 4.3. Plasma echo. We now investigate a much more intricate problem where the non-linearity plays the main role. Following the work [START_REF] Gould | Plasma wave echo[END_REF][START_REF] Malmberg | Plasma wave echo experiment[END_REF] or more recently [START_REF] Bedrossian | Suppression of Plasma Echoes and Landau Damping in Sobolev Spaces by Weak Collisions in a Vlasov-Fokker-Planck Equation[END_REF][START_REF] Grenier | Plasma Echoes Near Stable Penrose Data[END_REF], we will consider a perturbation of an homogeneous Maxwellian distribution f ∞ (x, v) := M(v) where

M(v) = 1 √ 2π exp - v 2 2 , v ∈ R ,
on the space interval [-L, L], with L = 6. Of course, this homogeneous Maxwellian is stable: for a small perturbation in high order Sobolev norms, we expect to observe the well known Landau damping on the macroscopic quantities. However, our aim is to investigate a transient regime where a plasma echo occurs for a well chosen perturbation. This echo appears when two waves with distinct frequencies interact. For a large time period, the effect is very small but at certain particular times, the interaction becomes strong: this is known in plasma physics as the plasma echo, and can be thought of as a kind of resonance [START_REF] Gould | Plasma wave echo[END_REF][START_REF] Malmberg | Plasma wave echo experiment[END_REF]. In all this section, we fix ε to 1. To build our initial condition, we proceed in two steps [START_REF] Filbet | Comparison of Eulerian Vlasov solvers[END_REF]. We first solve numerically the Vlasov-Poisson system with almost no collisions (τ 0 = 10 6 ) on the time interval [-30, 0] with an initial data at time t 0 = -30 given by fin (x, v)

= (1 + δ cos(k 1 x)) M(v) ,
where δ = 0.01 and k 1 = π/L. This choice induces an exponential decay of the potential energy by Landau damping at the rate associated to the perturbed mode k 1 , hence it gives a distribution function which is almost space homogeneous with small and fast oscillations in velocity. Then, at t = 0, we consider this solution, denoted by f (0), and choose it as initial data with a perturbation of the mode k 2 := 2k 1 . More precisely, we take

f 0 (x, v) = f (0, x, v) + δ cos(k 2 x) M(v)
. This initial data is represented in Figure 4.9. Then, starting from f 0 , we solve the Vlasov-Poisson-Fokker-Planck system on the time interval [0, 120].

On the one hand, we take τ 0 = 10 6 , which corresponds to a weakly collisional regime generating an oscillatory solution in velocity. For this reason, we choose a large number of Hermite modes N H = 8000 in such a way that our numerical results are not anymore sensitive to the numerical parameters. Let us emphasize that we compare our numerical results with those obtained using a semi-Lagrangian method [START_REF] Filbet | Comparison of Eulerian Vlasov solvers[END_REF][START_REF] Sonnendrücker | Vlasov simulations of beams with a moving grid[END_REF][START_REF] Filbet | Convergence of a finite volume scheme for the Vlasov-Poisson system[END_REF], which also requires such a large number of points in order to reach convergence. The phenomena is so intricate that we want to be sure that numerical parameters do not produce any artefact... Now let us comment our numerical results. In this weakly collisional regime, we expect that this initial data will induce a first Landau damping phenomena due to the perturbation of the second mode k 2 . However, after a time much longer than the inverse Landau damping rate, a new wave, called "echo" will appear as a modulation of the density at the wave number k echo = k 2 -k 1 = k 1 . This echo is due to the interaction between modes and is essentially a phenomenon of beating between two waves. First, we will see that the nonlinearity is crucial here. Indeed, in Figure 4.7, we show the time evolution of the potential energy and the first modes of the electric field obtained by applying the scheme (3.1) corresponding to the linearized Vlasov-Poisson-Fokker-Planck system and the scheme (2.11)-(2.16) corresponding to the nonlinear system. The numerical solution corresponding to the linearized system exhibits a simple Landau damping, when t ≥ 5, with a decay rate corresponding to the predicted value γ L = 0.355, whereas the numerical solution corresponding to the nonlinear system differs completely. It exhibits a first fast decay as for the linearized system, but when t ≥ 15, it produces an echo on the potential energy and the first mode (see right-hand side of Figure 4.7). The echo reaches its maximal amplitude at t = 30 for which we report the snapshots of the f -f ∞ in Figure 4.9. The first damping of the perturbed mode k 2 for short time t ≤ 5 and the subsequent echo are accurately reproduced. From [START_REF] Filbet | Comparison of Eulerian Vlasov solvers[END_REF], the echo wave number is indeed expected to be k echo = k 1 the first echo time is predicted at

t echo = t 0 + k 2 k 2 -k 1 (0 -t 0 ) = 30 ,
which corresponds very well with the numerical value we obtain here. From time t = 0 to t 20, the second wave corresponding to the mode k 2 has no effect on the first mode k 1 of the electric field, but at time t = 30, it is strongly perturbed by the echo effect. Actually, our numerical results illustrate the complexity of these phenomena.

On the one hand, we notice that echoes are repeated through time. On the potential energy (see the top right chart of Figure 4.7), a Landau damping is observed when 30 ≤ t ≤ 80, then we again discern a new echo around time t = 90 followed by a new damping. We also remark that this second damping (t ≥ 90) unfolds with a smaller decay rate than the first one (30 ≤ t ≤ 80). This repetition of echoes can be also perceived on the modes of the electric field (see the bottom right chart of Figure 4.7).

On the other hand, we report the time evolution of the first modes of the electric field (see the bottom right chart of Figure 4.7) and notice that other modes are also subjected to echoes but at times which differ from the "macroscopic" echo time t echo = 30. These echoes are not visible on the potential energy since their amplitude is smaller than the one of the potential energy by several order of magnitude. However, we point out that the third mode is subjected to a dramatic echo whose maximal amplitude, reached at time t = 15, is greater than the initial amplitude by a factor almost 10 5 . A careful inspection allows to distinguish the effect of this echo on the overall amplitude of the potential energy (compare bottom and top right charts of Figure 4.7). It is worth to mention that all modes corresponding to the linearized system are subjected to the classical Landau damping without any echo. This is a nice example where nonlinear effects, even small, induce intricate oscillatory behaviors. For this weakly collisional regime (τ 0 = 10 6 ), we also report (on Figure 4.8) the time evolution of the quantity

L 2 (t) := f -f ∞ L 2 (f -1 ∞ )
. First, we point out that unlike for potential energies, we do not observe any difference between the behavior of the linearized (3.1) and the nonlinear (2.11)-(2.16) solutions at the level of distribution functions on these charts. Second, we notice that on this time interval, collisions are negligible and we clearly see that the distribution function f does not yet converge to f ∞ . Figure 4.8 also shows that unlike in the case of strong Landau damping, variations of the spatial distribution occurs at an amplitude which is way smaller than the error between kinetic distributions. This can be also viewed in Figure 4.9, where we report the snapshots of the difference between the distribution function f solution to the nonlinear system (2.11)-(2.16) and the equilibrium f ∞ . We first observe the projection of the initial data which exhibits oscillations in velocity and a smooth perturbation in x with a small amplitude of order 10 -3 . At time t = 30 when the echo occurs, we clearly see that the first mode k 1 = 1 dominates, then the solution continues to oscillate due to the transport operator in a periodic domain in space.

-f∞ L 2 (f -1 ∞ ) , f -ρ M L 2 (f -1 ∞ ) and ρ -ρ∞ L 2 (ρ - 1
A natural question in physics is "how to cancel plasma echo?" for which a natural answer is that collisions may play a role, as shown in several recent articles [START_REF] Bedrossian | Suppression of Plasma Echoes and Landau Damping in Sobolev Spaces by Weak Collisions in a Vlasov-Fokker-Planck Equation[END_REF][START_REF] Chaturvedi | The Vlasov-Poisson-Landau system in the weakly collisional regime[END_REF]. To illustrate this phenomenon, we perform new numerical simulations passing from weakly to strongly collisional regime and again compare the two solutions corresponding to the the linearized system (3.1) and the nonlinear one (2.11)-(2.16). The results are now reported in Figure 4.10. Roughly speaking, when τ 0 > 10 2 , the nonlinear system exhibits a plasma echo whereas when collisions dominate, the electric field is rapidly damped and the solution converges to its equilibrium f ∞ exponentially fast as predicted by our analysis. It is worth to mention again that at the level of the distribution function, the L 2 time evolution of linearized (3.1) and nonlinear (2.11)-(2.16) solutions are globally the same for various regimes independently of the plasma echo. 4.4. Two-stream. In this last experiment, we fix ε to 1 and consider the equilibrium

φ ∞ (x) = 0.1 (1 -cos(k x)) , x ∈ (-L, L), with k = π/L with L = 6. The equilibrium is therefore given by f ∞ (x, v) = c 0 exp -φ ∞ (x) - |v| 2 2 ,
where c 0 is renormalizing constant. Furthermore, we consider the initial distribution function as

f (t = 0, x, v) = 1 6 √ 2π (1 + δ cos(k x)) (1 + 5v 2 ) exp - v 2 2 , (x, v) ∈ (-L, L) × R,
where δ = 10 -2 . These conditions can be viewed as a perturbation of data leading to the well-known two-stream instability when φ ∞ ≡ 0. For this case, we fix the number of Hermite modes at N H = 800 and consider the solution f to the nonlinear scheme (2.11)- (2.16). The purpose of this experiment is to compare our results with the classical two-stream instability which is usually performed with an homogeneous background distribution and to study the influence of collisions on our results. Our first comment is that unlike in the classical two-stream instability [START_REF] Bessemoulin-Chatard | On the stability of conservative discontinuous Galerkin/Hermite spectral methods for the Vlasov-Poisson system[END_REF], it is not clear that the electric field develops an exponential growth in this case. This may be observed on the left plot of where a non-constant stationary force field was applied in strongly collisional settings. However we deal here with a self induced force field whereas [START_REF] Blaustein | On a discrete framework of hypocoercivity for kinetic equations[END_REF] focuses on the linear case. These oscillations seem robust enough to persist in the present situation.

Conclusion and perspectives

In this work, we proposed a numerical scheme for the Vlasov-Poisson-Fokker-Planck model. On the one hand, we proved that our method is asymptotic preserving in the long time regime for the linearized model. To do so, we derived the exponential relaxation of the numerical solution towards its equilibrium with rates independent of scaling and discretization parameters. On the other hand, we tested the robustness of the method in various numerical experiments. These experiments show the accuracy of our method in both weakly collisional regime where small scales of the system are captured, allowing to reproduce filamentation, vortex formation as well as fine nonlinear phenomena such as plasma echoes but also in strong collisional regime, where we witness exponential trend to equilibrium, as predicted by our analysis of the linearized model. Many interesting perspectives arise from this work. On the theoretical view point, an important continuation consists in extending our theoretical results, which apply for a linear coupling with the Poisson equation, to the nonlinear scheme (2.11)-(2.16) by proving its asymptotic preserving properties and exponential trend towards equilibrium of discrete solutions. This might be doable in a perturbative setting by controlling the nonlinear contribution using discrete Sobolev inequalities. Carrying such proof in higher dimensions d = 2, 3 would be a great challenge and would probably require new theoretical tools. Indeed, equivalent studies on the continuous model in the literature rely on propagation of regularity [START_REF] Hérau | On global existence and trend to the equilibrium for the Vlasov-Poisson-Fokker-Planck system with exterior confining potential[END_REF][START_REF] Hwang | On the Vlasov-Poisson-Fokker-Planck equation near Maxwellian[END_REF][START_REF] Liu | Hypocoercivity based sensitivity analysis and spectral convergence of the stochastic galerkin approximation to collisional kinetic equations with multiple scales and random inputs[END_REF][START_REF] Herda | Large-time behavior of solutions to Vlasov-Poisson-Fokker-Planck equations: from evanescent collisions to diffusive limit[END_REF]. In our case it would require to propagate regularity at the discrete level. The groundwork towards such result has been laid in [START_REF] Blaustein | On a discrete framework of hypocoercivity for kinetic equations[END_REF], where we propagated discrete H 1 norms in the linear setting.

Another important continuation of this work is to incorporate nonlinear collisions to the model. Let us first observe that in [START_REF] Filbet | Fokker-Planck multi-species equations in the adiabatic asymptotics[END_REF], the Hermite spectral method is applied to a nonlinear Fokker-Planck operator conserving mass, momentum and energy. However, extending our analysis of the longtime regime at the discrete level to this case may require modifications and further investigations have to be done. L 2 -hypocoercivity methods have been applied in the case of nonlinear BGK and linearized Boltzmann operators at the continuous level [START_REF] Hérau | Hypocoercivity and exponential time decay for the linear inhomogeneous relaxation Boltzmann equation[END_REF][START_REF] Achleitner | On multi-dimensional hypocoercive BGK models[END_REF], however such analysis at the discrete level is not available in the literature in the framework of Hermite decomposition.

Regarding simulations, the study of echoes also raises interesting perspectives. In [START_REF] Grenier | Plasma Echoes Near Stable Penrose Data[END_REF] were constructed theoretical solutions to the Vlasov-Poisson equation which display infinite cascades of echoes and for which Landau damping is therefore not verified. Furthermore, sharp joint conditions on the collision frequency and the size of the initial data were obtained in order to ensure suppression of these echoes in [START_REF] Bedrossian | Suppression of Plasma Echoes and Landau Damping in Sobolev Spaces by Weak Collisions in a Vlasov-Fokker-Planck Equation[END_REF][START_REF] Chaturvedi | The Vlasov-Poisson-Landau system in the weakly collisional regime[END_REF]. Constructing such numerical solutions and illustrating the threshold obtained in these analysis would be of great interest. Another possible continuation would consist in finding a non-homogeneous background configuration where damping phenomena occur as in the homogeneous case, and then construct an experiment where nonlinear effect play the main role, even for small perturbation as for plasma echoes.

  for k ∈ {0, . . . , N H } and D n+1 k = 0 when k > N H . Observe that since ω

  0 , for k ∈ {0, . . . , N H } and D n+1 k = 0 when k > N H and k = -1. Moreover, we define the discrete free energy of the solution D n = (D n k ) k∈N to (3.1) as follows (3.

Figure 4 . 1 .

 41 Figure 4.1. Asymptotic-preserving properties: time development of (a) the potential energy (b) f -f∞

Figure 4 . 2 .

 42 Figure 4.2. Asymptotic-preserving properties: re-scaled time development (s ← ε t) of the potential energy (in log scale).

  3 the time evolution of the potential energy E p (t) := T |∂ x ψ(t, x)| 2 dx obtained using (3.1) in (a) and (2.11)-(2.16) in (b). Both solutions first produce fast damped oscillations up to time t ≤ 20 and then oscillate with a lower frequency while converging exponentially fast to zero with the same convergence rate γ L 0.004.

Figure 4 . 3 .

 43 Figure 4.3. Perturbation of non uniform density for τ 0 = 10 4 (weakly collisional regime): time development of the potential energy in log scale (for (a) the linearized Vlasov-Poisson-Fokker-Planck system and (b) the nonlinear Vlasov-Poisson-Fokker-Planck system.

Figure 4 . 4 .

 44 Figure 4.4. Perturbation of non uniform density (weakly collisional regime, τ 0 = 10 4 ): snapshots of the difference between the solution f and the equilibrium f∞ at time t = 4, 8, 16, 30, 40 and 70.

Figure 4 . 5 . 1 ∞

 451 Figure 4.5. Perturbation of non uniform density: time development of (a) the potential energy (b) f -f∞ L 2 f -1 ∞

Figure 4 . 6 .

 46 Figure 4.6. Perturbation of non uniform density (moderate collisional regime, τ 0 = 10 2 ): snapshots of the difference between the solution f and the equilibrium f∞ at time t = 4, 8, 16, 30, 40 and 70.

Figure 4 . 7 .

 47 Figure 4.7. Plasma echo for τ 0 = 10 6 (weakly collisional regime): time development of the potential energy (top) and square of the k-th mode of the electric field for k = k 1 , ..., 4k 1 in log scale (bottom) for (a) the linearized Vlasov-Poisson-Fokker-Planck system and (b) the nonlinear Vlasov-Poisson-Fokker-Planck system.

Figure 4 . 8 .

 48 Figure 4.8. Plasma echo for τ 0 = 10 6 (weakly collisional regime): time development of f-f∞ L 2 (f -1 ∞ ) , f -ρ M L 2 (f -1 ∞ ) and ρ -ρ∞ L 2 (ρ -1 ∞ ) in logscale for (a) the linearized Vlasov-Poisson-Fokker-Planck system and (b) the nonlinear Vlasov-Poisson-Fokker-Planck system.

  Figure 4.8. Plasma echo for τ 0 = 10 6 (weakly collisional regime): time development of f-f∞ L 2 (f -1 ∞ ) , f -ρ M L 2 (f -1 ∞ ) and ρ -ρ∞ L 2 (ρ -1 ∞ ) in logscale for (a) the linearized Vlasov-Poisson-Fokker-Planck system and (b) the nonlinear Vlasov-Poisson-Fokker-Planck system.

  Figure 4.11 considering the curves of the quantity E -E ∞ L 2 (T) in weak and intermediate collisional regimes τ 0 = 10 k ,

Figure 4 . 9 .

 49 Figure 4.9. Plasma echo for τ 0 = 10 6 (weakly collisional regime): snapshots of the difference between the solution f and the equilibrium f∞ at time t = 0, 20, 30, 40 and 50.

Figure 4 . 10 .

 410 Figure 4.10. Plasma echo for τ 0 = 1, . . . , 10 6 (various regimes): time development of the potential energy (top) and f -f∞ L 2 (f -1 ∞ ) (bottom) in log scale for (a) the linearized Vlasov-Poisson-Fokker-Planck system and (b) the nonlinear Vlasov-Poisson-Fokker-Planck system.

Figure 4 . 11 .

 411 Figure 4.11. Two-stream : time development of (a) E -E∞ L 2 (T) (b) and f -f∞ L 2 (f -1 ∞ ) for various τ 0 = 10 2 , . . . , 10 6 (in log scale).

  (a) τ 0 = 1 (b) τ 0 = 10 (c) τ 0 = 100 (d) τ 0 = 1000

Figure 4 . 12 .

 412 Figure 4.12. Two-stream : time development of f-f∞ L 2 (f -1 ∞ ) , f -ρ M L 2 (f -1 ∞ ) and ρ -ρ∞ L 2 (ρ -1 ∞ )for various τ 0 (in log scale).

  (a) τ 0 = 10 4 (b) τ 0 = 10 3 (c) τ 0 = 10 2

Figure 4 . 13 .

 413 Figure 4.13. Two-stream : snapshots of the distribution function f at time t = 8, 16, 30 and 60 for various τ 0 .

Acknowledgement

Both authors are partially funded by the ANR Project Muffin (ANR-19-CE46-0004). This work has been initiated during the semester "Frontiers in kinetic theory: connecting microscopic to macroscopic scales " at the Isaac Newton Institute for Mathematical Sciences, Cambridge, UK.