A generalized Fourier–Hermite method for the Vlasov–Poisson system

A generalized Fourier–Hermite semi-discretization for the Vlasov–Poisson equation is introduced. The formulation of the method includes as special cases the symmetrically-weighted and asymmetrically-weighted Fourier–Hermite methods from the literature. The numerical scheme is formulated as a weighted Galerkin method with two separate scaling parameters for the Hermite polynomial and the exponential part of the new basis functions. Exact formulas for the error in mass, momentum, and energy conservation of the method depending on the parameters are devised and L2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2$$\end{document} stability is discussed. The numerical experiments show that an optimal choice of the additional parameter in the generalized method can yield improved accuracy compared to the existing methods, but also reveal the distinct stability properties of the symmetrically-weighted method.


Introduction
The Vlasov-Poisson equations models the evolution of charged particles in their selfconsistent electric field. In this work, we consider the spectral discretization of the Vlasov-Poisson system. Compared to the widely used Particle-In-Cell (PIC) methods, which employ macroparticles that move throughout the computational mesh according to Newton's equations, spectral methods are not prone to statistical noise. Spectral methods naturally introduce increased computational costs per degree of freedom in exchange for improved accuracy. That is why, even though the first advancements in this area were made as early as 1963 [9], these methods became popular only recently due to the rapid increase of the computational power available. For test cases that are periodic, the Fourier basis is the natural choice for the spatial discretization. For the velocity discretization, one of the most common choices are Hermite-type basis functions. In this work, we introduce Hermite-type functions in a generalized form and build a mathematical framework for the solution in this basis. Alternatively, Fourier [18], Chebyshev [27], or Legendre [22] bases have also been used in velocity.

Previous Fourier-Hermite discretizations of the Vlasov equation
Hermite-type basis functions have been proposed already in early 1d1v Vlasov-Poisson simulations [2,3,15,17]. Later Holloway formalized two possible velocity discretizations of the Vlasov equation based on Hermite polynomials [16]. His first approach was based on standard Hermite functions as the basis in velocity. Then, the standard Galerkin method was used with Hermite functions as the test functions. This method is called symmetrically-weighted (SW) Hermite method. However, it turned out that mass and momentum cannot be conserved simultaneously for this method. To overcome that pitfall, the so-called asymmetricallyweighted (AW) Hermite basis was introduced In this case, in order to preserve orthogonality, another set of functions was used as test functions, For this method, it turned out to be possible to conserve mass, momentum, and energy exactly. For both methods scaling of the argument of the basis functions was considered. Certain choices of the scaling parameter proved to provide significant improvements in the quality of the results. It is consistent with the result of Boyd [5] that scaling of the series of Hermite functions is beneficial for accuracy. In addition, a shift in velocity is sometimes added depending on the initial value.
In the follow up work [26], Schumer & Holloway carried out a thorough numerical study of both methods which indicated that even though the AW method preserves mass, momentum, and energy, the SW method is more robust and better suited for long-time simulations. Despite that, most of the further developments were focusing on the AW method with a stabilization that is achieved through filtering of high modes (cf. e.g. [16,Sec. 3.3.3]). In [20,21] both spectral Galerkin and spectral collocation methods for the AW Fourier-Hermite discretization have been considered. Around the same time, the Hermite-based solution of the linearized Vlasov-Maxwell model has been investigated in [7]. A multi-dimensional spectral Vlasov-Maxwell solver based on the AW discretization has been proposed by Delzanno [8]. Camporeale, Delzanno, Bergen, & Moulton [6] demonstrated that for certain test cases the AW Fourier-Hermite method can be significantly more accurate than the PIC method. The spectral solver has been further enhanced by adding an adaptive strategy for regulating the number of basis functions [29], which is the basis for the SpectralPlasmaSolver code [30]. Moreover, the AW method was also considered for gyrokinetics [25].
On the other hand, some analysis was carried out for the SW case: Gibelli & Shizgal studied the convergence of the expansion of the distribution functions via Hermite functions in [12]. Convergence theory for the SW Fourier-Hermite method was provided in [23]. This study is, however, limited to a finite velocity interval.

Our contribution
We introduce a generalized Hermite basis based on results from the doctoral thesis [31], that contains the symmetric and the asymmetric basis as special cases. The contribution of this work is twofold: Firstly, a generalization of the method is introduced by the general formulation that includes two independent scaling parameters for the Gaussian and the polynomial part, which allows to include intermediate cases in the description. Secondly, we systematically study the conservation properties as well as the initialization procedure of the generalized Fourier-Hermite method with the aim of providing a mathematical foundation of the methods.

Organization of the paper
In the next section, we define our generalized Hermite basis. Section 3 introduces the Vlasov-Poisson model and describes its phase-space discretization in the generalized Fourier-Hermite basis. In Sect. 4, we find analytic expressions for the error in mass, momentum, and energy for the general method and discuss L 2 norm conservation. Finally, a numerical study of the method is provided in Sect. 5.

Generalized Hermite basis
As a basis in velocity space, we introduce the generalized Hermite basis of the form where ε > 0 is a suitable scaling parameter and γ > 0 defines the discretization method. We refer to the PhD thesis [31] for a detailed analysis of this class of functions. Note that we have slightly changed the definition such that our parameter γ corresponds to γ ε in [31]. Other than the standard Hermite functions, these functions are not orthonormal in the standard L 2 (R) inner product. Instead, the basis is orthonormal in the following weighted space, and the inner product We later refer to ·, · L 2 ω (R d ) as ·, · ω . Let us note that the symmetrically and asymmetrically weighted Hermite bases, known from the literature, are included in our generalized setup. The symmetric basis is obtained-up to a constant scaling factor, that does not influence the corresponding spectral Galerkin scheme-by the choice The asymmetrically-weighted basis is obtained for From the orthogonality and completeness of the Hermite functions in L 2 , the orthogonality and completeness of the generalized Hermite functions (2.1) in the weighted L 2 space can be deduced as stated in the following Lemma. From the recursion relation of the Hermite polynomials (see [11, §6.4, p. 186]), we deduce the following recursion relations for the generalized Hermite functions: Let us now derive a couple of other properties that involve integrals of the basis function over R. These formulas will be useful for the computation of the observables. Analogous formulas for the SW and AW cases were considered in [16,26]. We now derive them for the generalized setup.

Lemma 2.2 Denote
Then, the following relations hold: 1.
Proof First, we note that even (odd) Hermite polynomials are of even (odd) functions and so are the generalized Hermite functions. Therefore, even moments of the basis integrate to zero for odd indices, and odd moments for even indices.
ε . For 2 , > 0, we use the following result from [1,Expr. 22.13.17], With the coordinate transformationṽ = εv, we can use this result to find With the analogous expression for I 2 +2 , the fraction (2.5) can be found. 2. The formula for J 1 follows from [14, 3.
Applying the same formula to J −1 , we find This gives an expression for I in terms of J −1 which can be inserted to the first part of (2.6) to obtain its second part. 3. The formula forĪ 0 follows from [14, 3.

Remark 2.1
We note that the AW case has the very special property that the mth order moment of H 1,ε vanishes for all > m. This property is particularly useful when filtering of higher modes is applied (cf. [7,Sec. 3]), since then the low order moments of the solution are not affected. In the same way, a velocity shift can be added to the Hermite polynomials without destroying conservation (cf. [6]).

Phase-space discretization with the generalized Hermite-Fourier method
In this section, we first introduce the Vlasov-Poisson model, followed by a discussion of the phase-space discretization and the projection of the initial value.

The Vlasov-Poisson model
The d-dimensional Vlasov-Poisson system for the distribution function f : Ω ×R d × R + → R for electrons in neutralizing background on the domain Ω ⊂ R d reads as, where the electric field E(x, t) is computed via the Poisson equation and the density is given by Note that we use dimensionless units. In this paper, we will restrict ourselves to the one dimensional case on a periodic domain [0, L x ] in x and the whole space R in velocity. To close the system, some initial distribution f 0 (x, v) need to be specified. The Vlasov-Poisson system exhibits a lot of structure and many quantities are conserved over time (see, for example, [28, § 3.2.2]). In particular, the following observables are conserved over time: We now introduce a spectral discretization and discuss the conservation of these observables on the semi-discrete level.

Discretization in velocity
We now proceed to the discretization of the Vlasov equation in velocity. We look for an approximation f N v of the distribution function f in the approximation space In order to treat the coefficients with the indices outside of the range 0, . . . , N v − 1, we need a closure scheme. The most obvious choice is to assign them to zero This is the most common closure scheme but the alternatives can e.g. be found in [9]. Inserting expression (3.2) into the Vlasov equation (3.1) and exploiting the relations (2.3) and (2.4), yields We now discretize the Vlasov equation by a Galerkin method in the weighted space also as test functions and using expression (3.3) together with the orthonormality of the basis, we arrive to the following system of partial differential equations for the coefficients c (x, t), = 0, . . . ,

Discretization in space
We now discretize our system in space via the Fourier basis. For all = 0 . . . N v − 1, we take the following ansatz, For the convenience of notation, we consider only an odd number of Fourier modes. For an even number of modes, when k = −N x , . . . , N x +1, the coefficients corresponding to the indices 0 and N x + 1 have to be tackled separately to preserve the symmetricity, but all our results equally apply. The electric field E can be represented in the Fourier basis as The product E(x, t)c (x, t) is given in terms of the corresponding Fourier coefficients as k=−N x and the vector E(t) is additionally padded with zeros for all other indices. Here * denotes the convolution With a standard Galerkin method in Fourier space, we get the following system of differential equations for the coefficients: We will later refer to the discretization (3.7) as the generalized Fourier-Hermite discretization.

Computation of the electric field
To complete the formulation of the method, it is now left to compute the coefficients {E m (t)} N x m=−N x of the electric field E in the Fourier basis. From the representation of the distribution function in (3.5), we get A spectral Galerkin scheme for the Poisson equation gives the following Fourier representation of the electric field and the constant mode should vanish for the electric field E 0 (t) = 0.

Representation of the initial distribution in the generalized Fourier-Hermite basis
For the complete setup of the method, we need to discuss the representation of the initial distribution function in the generalized Fourier-Hermite basis. This is obtained by a (weighted) Galerkin projection to the discrete spaces. In many cases, the initial distribution is separable in space and velocity variable, and the spatial part is given by a trigonometric function and the velocity part by a (sum of) shifted Gaussians. Representing a trigonometric function in Fourier space is trivial. Let us therefore consider the weighted Galerkin projection velocity space for a function of the form, This task is particularly simple for a single Gaussian centered at zero, which is exactly represented by the 0th basis function if the width of the basis is chosen as ε = 1 For the general case, the generalized Fourier coefficients can be computed as a linear combination of the coefficients for the individual Gaussians according to the following proposition.

Proposition 3.1 For a Gaussian g n
, the generalized Fourier coefficients of can be computed according to the following formula Proof With the definition of the function H γ,ε and some simplification, we get We observe that In order to evaluate the integral, we use the formulas [14, §7.374, Expr. 8 and 6] for into (3.9), together with the definition ofγ n , yields the assertion.
Note that we only consider the case of a single set of basis functions as opposed to an ansatz with multiple sets of basis functions for each Gaussian proposed e.g. in [6]. For this reason, we do not consider a shift in velocity space.
The optimal choice of the parameter with respect to the initial value is not necessarily the best choice of the parameter after time has evolved (cf. the findings in [10] in the context of a semi-Lagrangian Hermite discretization). The accuracy of the method can therefore be improved by readapting the parameter to time evolution of the velocity moments of the solution.

Temporal discretization
The main focus of our paper is on the spatial discretization. For the temporal discretization, we choose the following form of a mid-point rule to compute the numerical approximation c k,m at time t m from the approximation at time t m−1 = t m − Δt: This form of the temporal discretization was also chosen in [7] for the AW method. The advantage of this scheme is that the conservation of energy translates from the semi-discrete to the fully discrete system as we will see in Sect. 4.6.

Computation of the observables and conservation properties
In this section, we provide formulas for mass, momentum, energy, and L 2 norm for the numerical scheme and analyze how they evolve over time.

Telescopic sums
As a preparation, we first compute some telescopic sums in our generalized setup that will be useful in the sequel. Note that Holloway [16] already used the idea of telescoping sums to show conservation for his two special cases.

Lemma 4.1 (Telescopic sums) Consider a sequence {a }
Let S J and SĪ be the same sums as above but with the terms J orĪ instead of I , respectively. Then, the following expressions hold for these sums Then, telescoping the sum above and noting that we arrive to the formula (4.1).
The formula for S J can be proved in similar fashion, but this time we express J through I −1 with the help of (2.6).
We note that the first term is a telescopic sum and can be computed as a I in this case. On the other hand, for odd N v , we can borrow the missing term −I N v −1 a N v −1 from the result of the telescopic sum. Putting it all together we arrive at the expression (4.2).
For the last sum SĪ let us first use the representation (2.7) ofĪ through J −1 . This yields two terms similar to the ones obtained for S I and S J that can then be shown to vanish with the same arguments.

Mass
For a solution represented in the generalized Fourier-Hermite basis, the semi-discrete mass is defined as Analytically, the mass is conserved and for the semi-discrete version we have the following result.
Proof Using the representation (4.4) together with expression (3.7), we get: where we have used the notation from (3.6) for the convolution. Using expression (4.1) for the sum above with {β 0 (t)} in place of the sequence {a }, we get Keeping in mind that the closure scheme implies c 0 N v (t) = 0 for all t ∈ R + , and hence β 0 N v (t) = 0, for all t, we arrive to the expression (4.5). Remark 4.1 Note that an analogous expression also holds on the semi-discrete level in phase space, i.e. after the velocity discretization while the spatial variable is left continuous.

Momentum
For our generalized Fourier-Hermite basis representation, the semi-discrete momentum is given by J c 0 (t). (4.6) and evolves according to the following theorem.

Theorem 4.2 (Fourier-Hermite momentum evolution)
The momentum is conserved for even numbers of basis functions. Moreover, Proof Using the representation (4.6) together with the expression (3.7), we get Using (4.2) for the sum above with {β 0 (t)} in place of the sequence {a }, we get Since the closure scheme implies c 0 N v (t) = 0, and hence β 0 Let us now consider the sum N v −1 =0 β 0 (t)I . Using formula (3.8) for the Fourier coefficients E k (t) and E 0 = 0, we get

Energy
Similarly to the other observables, the kinetic and the potential energies of the approximated solution can be computed as

Theorem 4.3 (Fourier-Hermite energy evolution) For odd N v , the energy is conserved,
i.e.

For even N v , the time derivative of the energy can be computed as
Proof For the derivative of the electric energy, we need to consider the derivative of the Fourier coefficients of the electric field. Using expression (3.8) for the coefficients {E k } and equation (3.7) for the coefficients c k (t), we get, for k = 0, Using (4.1) for the second sum, we get, for k = 0, where we used the fact that according to the closure scheme c k N v (t) = 0 for all t > 0, k = −N x , . . . , N x . As for the first sum, using the expression (2.6) of I through J −1 and J +1 , we get where we used that and we added the term J 0 c k 0 (t) to the sum for both odd and even N v since J 0 = 0. The term J N v −1 c k N v −1 (t) we added for odd N v only, since J N v −1 = 0 for even N v . Introducing the notation ΔM k N v (t) and putting everything together yields for even N v Let us now consider the kinetic energy. Using (3.7), the representation can be expressed as Using the expression (4.3) for the sum above with {β 0 (t)} in place of the sequence {a }, we get where we used again that β 0 N v (t) = 0 for all t to eliminate the first term of the telescopic sum for odd N v .
Summing up the terms for kinetic and potential energy, we see that they cancel out for odd N v and find expression (4.8) for even N v .

L 2 norm
Finally, let us consider the L 2 norm. Compared to the moments of the distribution function, that we have considered so far, the computation of the L 2 norm is more involved, since we have a quadratic dependence on the distribution function. Let us therefore first derive an expression for the semi-discrete L 2 norm. (4.10) with Γ and 2 F 1 being the gamma and hypergeometric functions, respectively.
Proof By the bilinearity of the inner product, we get where we have used the orthogonality of the Fourier basis in the second step. Now, we only need to compute the products H For all other cases, we transform the integral to If 1 + 2 is odd, the integrand is odd and, hence, these terms vanish. When 1 + 2 is even, the expression (4.10) can be derived from [4,Eq. (1.5)].

Remark 4.2
If the first two arguments of the hypergeometric function are negative integers, the following form can be derived from [14, 9.14] where (·) n denotes the Pochhammer symbol. This gives the following expression for the L 2 product of the generalized Hermite functions For small values of ξ 1, we observed the following asymptotic behavior (4.11) The limit ξ → 0 corresponds to γ → √ 2, i.e., the limit of the SW case.
The conservation of the L 2 norm is crucial for the numerical stability of the method. Let us first consider the special case of γ = √ 2, which corresponds to the symmetrically weighted case. Since then the weight is constant, a function f is in Proof Let us denote the differential operator of the Vlasov equation by L, i.e., we write the Vlasov equations as ∂ t f + L f = 0. By an integration by parts argument, one can show that Following [13, § 2, Eq. (2.5)], we write the Galerkin approximation as where P N v ,N x is the Galerkin projection operator that is selfadjoint in L 2 ((0, L x )) × L 2 ω SW (R), and hence in L 2 ((0, L x )) × L 2 (R). Therefore, we find The proof relies on the two major properties f , L f L 2 ((0,L x ))×L 2 (R) = 0 and . If we take the weighted Galerkin projection with a non-constant weight instead, which is the case for all setups except for the SW one, we cannot transfer from the L 2 norm of f N v ,N x to the weighted one by means of multiplication with a constant. Therefore, the Galerkin projection operator is no longer self-adjoint. Indeed, the numerical results show that the L 2 norm is not conserved for all other choices of γ . Alternatively, one could prove the L 2 norm conservation for the SW case starting from formula (4.9) using the differential equation (3.7) to replace the derivative. Then the first part vanishes by a telescoping argument similar to the proofs of the other conservation properties. The second part is not a telescopic sum but the two parts sum up to a term multiplied by γ 2 − 2. Hence, also this part vanishes for the SW case. However, for all other cases, the second term remains. Moreover, we will have contributions from the off-diagonal terms in (4.9). Using (4.11) one can show that the terms with 1 = 2 and | 1 − 2 | = 2 have a linear dependence on the deviation of γ from γ SW in the limit γ → γ SW , when neglecting the dependence of the coefficients on γ . Indeed, the numerical results reported in Sect. 5.4 verify the linear decay close to the SW case.

Conservation properties of the fully discrete system
Conservation of mass and momentum for N v odd or even, respectively, readily translates to a time-discrete approximation. Energy conservation requires more structure but is satisfied for the scheme proposed in Sect. 3.6 as the following theorem shows: (3.11) conserves energy, i.e.,

Theorem 4.4 For odd values of N v the Fourier-Hermite discretization with a mid point temporal discretization given by
Proof With the same arguments as in the proof of Theorem 4.3, we find that and E m+1 With the latter expression for E m+1 k , we find This proves the assertion.

Comparative summary
For the generalized Fourier-Hermite discretization of the Vlasov-Poisson equation with a zero closure of the truncated expansion, all even (odd) moments of the distribution function are conserved over time for an odd (even) number of basis functions. In terms of the conservation of moments, the asymmetrically-weighted setup is special, since mass, momentum, and energy are conserved due to the fact that I = J =Ī = 0 for all indices > 2. In particular, the mass, momentum, and energy are defined by the 0th, 1st, or 0th and 2nd AW Hermite function, respectively. For this reason, it is only possible for the asymmetrically weighted method to damp higher modes or add a velocity shift without affecting the conservation properties. On the other hand, the L 2 norm is only conserved for the symmetrically-weighted case, where the discretization corresponds to a non-weighted Galerkin scheme. This is of special importance, since L 2 norm conservation is linked to numerical stability. For all other choices of the method, a numerical stabilization is required by e.g. adding (hyper)viscosity or filtering of high modes. Since the effect of such a stabilization is unphysical, it is desirable to keep its amount as small as possible. From the dependence of the L 2 norm deviation on the parameter γ , it is clear that the norm error decreases to zero as γ = √ 2 is approached. This fact is also confirmed in our numerical experiments. For this reason, the further γ from the symmetric case, the more stabilization needs to be added. The choice of an intermediate value of γ between γ AW = 1 and γ SW = √ 2 offers the possibility to optimize the method parameter in such a way that a good compromise between the amount of stabilization that is necessary is balanced against the amount of violation of the conservation properties.

Numerical results
The generalized Fourier-Hermite method has been implemented in MATLAB. The parameter ε of the generalized Hermite basis, corresponding to the width of the Gaussians, is chosen based on the initial conditions of the corresponding test cases. The parameter γ , however, is left free and its influence is studied. For propagation in time with the conservative mid-point scheme as described in 3.6 we use a time step of Δt = 0.01 in all our simulations.

Test cases
First, we introduce three standard test cases that we will use for our numerical study.

Landau damping
The initial distribution for this case takes the following form For our numerical experiments, we take α = 0.05 and k = 0.5. For the verification of the results, we compare the damping rate to the rate of −0.1533 predicted by linear theory.
Since we only have one Gaussian, we fix ε = 1 √ 2 for the Hermite basis in order to match the width of the Gaussian. Then, the initial value is exactly represented with coefficients For the Landau damping, we choose a resolution of N x = 32-that is 65 modes-in space and N v = 32 or N v = 33.

Two-stream instability
The following initial distribution yields the two-stream instability test case with α = 0.001, k = 0.2 and v 0 = 3. In this case, the linear theory predicts a growth rate of 0.569. We set ε = 1 √ 2 that matches the width of the Gaussians and use Proposition 4.1 to approximate the initial value. We choose a resolution of N x = 32 in space and N v = 64 or N v = 65 in velocity.

Bump-on-tail instability
Finally, we consider the initial distribution , and v b = 4.5. For the initialization we use the result of Lemma 3.1 with the corresponding parameters. We set the width of the Gaussian to ε = 0.83, which proved to give good results in this case. The resolution is chosen as for the two-stream instability.

Verification of the method for an intermediate
First, we set γ = √ 2 − 0.1, which corresponds to an intermediate case between γ AW = 1 and γ SW = √ 2, and compare the results to linear theory. In Fig. 1a, the evaluation of the electric energy over time is shown for simulations of the bump-ontail instability with N v = 64 and N v = 65, respectively. During the linear phase both simulations match the results predicted by linear theory (cf. [24]). Moreover, we compare the deviation in mass, momentum, and energy with the values obtained by the analytic predictions derived in (4.5), (4.7), and (4.8), which we have propagated with the same temporal approximation as the original problem. Figures 1b-d show that the observables are indeed preserved up to machine precision where predicted and otherwise match the analytical formulas (up to machine precision). We have repeated the same experiments for the other two test cases as well with the same results (not shown in the paper).
For the Landau damping test case, we have shown in Fig. 2a the evolution of the electric energy for the SW, the AW, and an intermediate method. We observe that the exact damping rate is recovered the longest for the intermediate value of γ , while the well-known numerical recurrence phenomenon occurs the latest for the SW case. This shows that the choice of γ influences the quality of the solution and intermediate values can be advantageous. Note that the results can be improved by adding an articifical collision term as shown in [6]. Figure 2b shows the evolution of the L 2 norm for the Landau damping test case for three values of γ . As predicted, the L 2 norm is only conserved in the SW case, however, no numerical stability issues were encountered for this test case.

L 2 stability
For the other two test cases, the choice of the parameter γ clearly affected the numerical stability of the algorithm. In order to illustrate this, we run simulations for γ ∈ [ε √ 2 − 0.4, ε √ 2 + 0.25] with a step of 0.02 and report in Figs. 3a, b the time t max ∈ (0, 100] that could be reached before the error in the L 2 norm exceeded 10 (with analytic values between 2 and 3). One can see in the figures that the closer we

Conservation of the observables
In this section, we study the conservation properties of the method depending on γ . We simulate until time 30 and report the maximum deviation in mass, momentum, energy, and L 2 norm in this time interval as a function of γ . Here, we vary γ ∈ [−0.3, 0.05] in the two-stream instability test case and γ ∈ [−0.15, 0.05] in the bump-on-tail test case with a step of 0.01 in both cases. In light of the results of the previous section, we consider the simulation to be numerically stable in this range. Note that the range does not include the AW case for the bump-on-tail test. Figures 4a and 5a show the results for the two test cases with an even number of basis functions. As predicted, the odd moment vanishes. The even moments are conserved for the AW case and, remarkably, also for values close to this case, but show a more or less monotonic increase towards the SW case and beyond. For an odd number of basis functions, the roles of odd and even moments are swapped as shown in Figs. 4b and 5b. Only for the two-stream instability, where the momentum is zero due to the symmetry in the initial value, the momentum is numerically conserved as well.
The behavior of the L 2 norm shows a moderate increase for values of γ away from the SW case. However, close to the SW case it shows a sharp decay, since the SW case is the only choice with L 2 norm conservation. We therefore show a zoom of the L 2 norm conservation in Fig. 6 in a double logarithmic plot. Close to the SW case,      Finally, let us observe that the conservation of mass and energy is still dominated by roundoff erros for γ around 1.2 while the L 2 norm is already several orders of magnitude improved for this range of values of γ for the two-stream instability test in Fig. 4a.

Conclusions and outlook
In this paper, we have studied a generalized Fourier-Hermite discretization of the Vlasov-Poisson equation. We have derived exact formulas for the error in mass,  The parameter ε, corresponding to the width of the Gaussian in the basis, is usually adapted to the initial distribution. The parameter γ , however, allows to vary the scaling of the argument of the Hermite polynomials in the basis, which in practice means that it controls the ratio between the scaling of the exponent and the Hermite polynomial in the basis. The AW and SW methods are just cases of specific values of γ . The general theoretical framework derived in this paper also highlights the special properties of the SW and AW cases. The SW setup that implies working in a standard L 2 (R) space is the only method offering exact L 2 norm conservation. Moreover, the new generalized approach devised in this paper allowed us to demonstrate that this is indeed a distinct property of the SW method and even a small deviation yields substantial loss in L 2 norm conservation. On the other hand, only the AW method allows for simultaneous conservation of mass, momentum, and energy. However, as we have seen in our numerical study, deviation from the AW method does not immediately yield considerable losses in terms of conservation. Moreover, the analytical formulas for the loss in conservation can potentially be used to adapt the size of the generalized Hermite basis during simulation and for an optimization of the method parameter. Monitoring the loss in L 2 norm conservation, on the other hand, may be used to automatically adjust the hyperviscosity in stabilizing methods other than the SW method.
Further directions of future research include an extension of the method to multiple dimension based on the anisotropic extension of the basis that was introduced in [19] in the context of radial basis function stabilization. The anisotropy is supposed to be beneficial in a setup with a guide field causing an anisotropy in the velocities parallel and perpendicular to the field.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, 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 article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article'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. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.