On the Use of Hermite Functions for the Vlasov–Poisson System

We apply a second-order semi-Lagrangian spectral method for the Vlasov–Poisson system, by implementing Hermite functions in the approximation of the distribution function with respect to the velocity variable. Numerical tests are performed on a standard benchmark problem, namely the two-stream instability test case. The approach uses two independent sets of Hermite functions, based on Gaussian weights symmetrically placed with respect to the zero velocity level. An experimental analysis is conducted to detect a reasonable location of the two weights in order to improve the approximation properties.

velocity domain. The solution at each space-velocity node is traced back along the characteristic curve originating backward from that node. In [8] a high-order Taylor expansion of the characteristic curves is used to trace back the solution in time, which is then approximated by spectral interpolation. Such a method guarantees the conservation of the main physical quantities (charge, mass, and momentum).
The first attempt in using Hermite polynomials to solve the Vlasov equation dates back to the work [10], where the Hermite basis is used in the velocity variable to describe a plasma in a physical state near the thermodynamic equilibrium. Within this approach, exact discrete conservation laws can be constructed [7,13,14,20,21]. The weight function of the Hermite basis can be generalized by introducing a parameter α in such a way that it becomes exp(−α 2 v 2 ). A proper choice of this parameter can significantly improve the convergence [2,3,19]. This fact was also confirmed in earlier works on plasmas physics based on Hermite spectral methods (see [11,17] and more recently [4]).
The paper is organized as follows. In Sect. 2, we present the continuous model, i.e., the 1D-1V Vlasov equation. In Sect. 3, we introduce the spectral approximation in the phase space. In Sect. 4, we present the semi-Lagrangian schemes based on an approximation of the characteristic curves coupled with a second-order backward differentiation formula (BDF). In Sect. 5, we numerically assess the performance of the method for a standard test case, and we show how the solution's behavior can be affected by the choice of a certain parameter β, acting on the location of Hermite weight function.

The Continuous Model
We deal with the 1D-1V Vlasov equation defined in the domain = x × R, with x ⊆ R. The unknown f = f (t, x, v) denotes the probability of finding negative charged particles at the location x with velocity v. This is solution of the problem where ρ denotes the electron charge density. System (1)- (2) in the unknowns f and E is a simplification of the Vlasov-Poisson equations in two or three dimensional space domains. Uniqueness of the solution is ensured by imposing that where | x | is the size of x . We assume periodic boundary conditions in the variable x and a suitable exponential decay at infinity for the variable v. After integration and by using the boundary constraints, we obtain the conservation of mass When f and E are smooth enough, for a sufficiently small δ > 0, the local system of characteristics associated with (1) is given by the curves (X(τ ), V (τ )) solving with the condition that (X(t), V (t)) = (x, v) when τ = t. With this setting we have in mind that for τ > 0 we proceed backward. Under suitable regularity assumptions, there exists a unique solution of the Vlasov-Poisson problem (1)-(2) which is formally obtained by propagating the initial condition along the characteristic curves described by (5), i.e. we have where we recall thatf is the initial datum. By using the first-order approximation the Vlasov equation is satisfied up to an error decaying as |τ − t|, for τ tending to t.

Phase-Space Discretization
We briefly recall the construction of the approximation method proposed in [8].
where δ ij is the usual Kronecker symbol. We recall that Hermite functions are obtained from Hermite polynomials after multiplication by the weight ω(v) = e −v 2 . We also define the discrete spaces Any function f N,M that belongs to Y N,M can be represented as where the coefficients of such an expansion are given by In the following, the matrices d As a special case, we set d To evaluate a function f N,M ∈ Y N,M at the new points (x nm ,ṽ nm ) through the coefficients c ij , we use a Taylor expansion in time. By omitting the terms in t of order higher than one, we get By substituting (11) in (9), we obtain the approximation which is the main building block for more advanced schemes.

Discretization of the Vlasov Equation
Given the time instants t k = k t = k T /K for any integer k = 0, 1, . . . , K, we consider the approximation of the unknowns f and E of problem (1)- (2), given by where the function f Hence, at any time step k, we express f where c Suppose that E (k) N is given at step k. According to [8], we write where the discrete Fourier coefficientsâ The distribution function f is expected to remain constant along the characteristics. The most straightforward discretization method is obtained by advancing the coefficients according to the approximation This states that the value of f (k+1) N,M , at the grid points and time step (k + 1) t, is assumed to correspond to the previous value at time k t, recovered by going backwards along the characteristics. To computeṽ nm , we should use E N (x n ). However, the distance between these two quantities is of the order of t, so that the replacement has no practical effects on the accuracy of firstorder methods. Between each step k and the successive one, we need to update the electric field. This can be done by using the Gaussian quadrature formula in (14), so obtaining where w j , for j = 1, . . . , M − 1, are the quadrature weights. Afterwards, in order to compute the new point-values E (k+1) N (x n ) of the electric field, it is necessary to integrate ρ (k) N . By using approximation (12) in (17), we end up with the first-order explicit scheme of Euler type: where The parameter t must satisfy a suitable CFL condition, which is obtained by requiring that the point (x nm ,ṽ nm ) falls inside the box ]x n−1 , A straightforward way to increase the time accuracy is to use a multistep discretization scheme as the second-order accurate two-step BDF scheme. We have where (x nm ,ṽ nm ) is the point obtained from (x n , v m ) going back of one step t along the characteristic lines. Similarly, the point (x nm ,ṽ nm ) is obtained by going two steps back along the characteristic lines, i.e., by using 2 t instead of t when computingx nm andṽ nm . Despite the fact that a BDF scheme is commonly presented as an implicit technique, in our context (f constant along the characteristics) it assumes the form of an explicit method. In terms of the coefficients, we end up with the scheme From theoretical considerations and the experiments in [8], it turns out that the above method is actually second-order accurate in t. Higher order schemes can be obtained with similar principles. All the above schemes guarantee mass conservation (see (4) for the continuous case), which is a crucial physical property. For practical purposes, it is advisable to make the change of variable f (t, A generalization consists in introducing a real parameter α and assuming that the weight function is ω(v) = exp(−α 2 v 2 ). The approximation scheme can be easily adjusted by modifying nodes and weights of the Gaussian formula, through a multiplication by suitable constants. The difficulty in the implementation is practically the same, but, as observed in [9], the results are quite sensitive to the variation of α.

Numerical Experiments
The numerical scheme here proposed is validated in the standard two-stream instability benchmark test. We consider the Vlasov-Poisson problem (1)-(2) where we set 5]. The initial solution is given bȳ where G R (v) = e −α 2 (v−β) 2 and G R (v) = e −α 2 (v+β) 2 are two Gaussians centered symmetrically at the points v = ±β. The parameters for (24) are: a = 1/ √ 8, In all the experiments that follow, we integrate up to time T = 30 using the second-order BDF scheme with a suitably small time step, in order to guarantee stability and a good accuracy. In this way we can concentrate our attention to the spectral approximation in the variable x and v. A study of the convergence rate in time of the proposed numerical scheme can be found in [8]. First of all, in Fig. 1 we show the results at time T = 30 of the solution recovered by the Fourier-Fourier method, by choosing N = 2 5 , M = 2 6 and time step equal to t = 0.00125. This will be the referring figure for the successive comparisons. Besides we show the corresponding time evolution of |â (k) 1 |, the first Fourier mode of the electric field E (k) N in (16). The behavior of this last quantity is predicted by theoretical  (16) considerations, and the slope of the "segment" starting at T = 15 agrees with the expectancy [1, Chapter 5]. As done in [9], we perform a series of experiments using less degrees of freedom than those actually necessary to resolve accurately the equation. In practice, we set N = M = 2 4 . In this way, we could for instance detect what happens by varying the parameters α and β. Of course, if we increase the number of degrees of freedom, the numerical solution improves and cannot be distinguished from the referring one shown in Fig. 1. The purpose in [9] was to check what happens by varying the parameter α in the Hermite weight exp (−α 2 v 2 ). The conclusions are that the approximate solution is very sensitive to the choice of α and that there are values of α that perform better than others. In general these values are those belonging to a neighbourhood of α = 1. Moreover, in [9], we note that keeping α constantly equal to the value that better fits the initial datum (i.e. α =ᾱ = 2 for (24)) may create instability as time increases. For such motivations, since at the moment a practical algorithm able to vary α in a dynamical way during the computations is not available, in the numerical experiments that follow we fix α = 1, while play with β.
Due to the particular initial condition, we adopt a two-species decomposition of the Vlasov equation, where the distribution function is given by the sum of two electron distribution functions, i.e., f = f R + f L . These distribution functions refer to the two initial electron distributions, so that f R = p R G R and f L = p L G L , where p L and p R are given polynomials. We consider the two systems of electrons described by the distribution functions f L and f R at the initial time as distinct plasma species that maintain their diversity throughout the whole numerical simulation. Therefore, we can split the Vlasov equation into two equations that are still of Vlasov type and are solvable independently, although they are coupled through the same electric field, which depends on the total charge density. This amounts to approximate two independent equations of the same type of that given The two unknowns are then coupled through the density function as in (2). The plots of Fig. 2 (16). The distribution functions presented in the left column of Fig. 2 are visibly and significantly different depending on β, while the first Fourier mode of the electric field shown in the right column seems to be less affected. These differences practically confirm that the choice of the Hermite weight functions ω(v) = exp(−α 2 (v ± β) 2 ) is a crucial aspect of the method (see also [11,12,17,22]). This conclusion is heuristic. Unfortunately, there is no space enough for a deeper quantitative analysis in these pages. The question deserves however further investigation. Moreover, it would be advisable to develop appropriate algorithms allowing for the automatic adjustment of both parameters α and β during the time advancing procedure, in order to optimize the performance.
Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence and indicate if changes were made.
The images or other third party material in this chapter are included in the chapter's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.