Long-time behavior of second order linearized Vlasov-Poisson equations near a homogeneous equilibrium

The asymptotic behavior of the solutions of the second order linearized Vlasov-Poisson system around homogeneous equilibria is derived. It provides a fine description of some nonlinear and multidimensional phenomena such as the existence of Best frequencies. Numerical results for the 1D x 1D and 2D x 2D Vlasov-Poisson system illustrate the effectiveness of this approach.

1. Introduction. We consider potentials φ = φ(t, x) : R×T d → R and distribution Here periodic boundary conditions being used, T d is a d dimensional torus: there exist L 1 , . . . , L d > 0 such that T d = (R/L 1 Z) × · · · × (R/L d Z). Furthermore, doing an assumption of neutrality, we only consider solutions of (VP) such that In this paper, we aim at exhibiting nonlinear and multidimensional phenomena of solutions of (VP), pursuing a first preliminary work [4] on this subject. Beyond their physical interest, these phenomena can be relevant to evaluate the performances and the qualitative properties of numerical methods.
Since the very first developments of numerical methods for solving VP (we refer to [15], for a review; the literature is particularly huge in 1D × 1D and we can mention [14], as one of the earliest works in 2D × 2D), the numerical solutions are compared to the solutions of the Vlasov-Poisson system linearized around a homogeneous equilibria f eq ≡ f eq (v) 1 . It consists in looking for solutions of (VP) of the type f = f eq + εg φ = 0 + εψ. Neglecting second order terms, g is formally a solution of This equation being linear and homogeneous, it is natural to try to solve it realizing a Fourier transform we respect to the variable x. Thus, we get where k ∈ T d = (2π/L 1 )Z × · · · × (2π/L d )Z and the Fourier transform with respect to the space variable is defined for u ∈ L 1 (T d ) and k ∈ T d by It is relevant to notice on (VPLF) that there is no energy exchange between space modes at the linear level. In other words, if g is a solution of (VPL) such that g 0 ≡ g 0 (v)e ik·x then it is of the form g(t, x, v) = g(t, v)e ik·x . As a consequence, the linear analysis is not well suited to describe multidimensional phenomena that could be confronted with numerical simulations. Since (VPLF) is linear and autonomous, it is natural to solve it with the Laplace transform. This transform is defined for functions u : R * + → C such that there exists λ ∈ R satisfying ue −λt ∈ L ∞ (R * + ) and z ∈ C such that z > λ by Thus, it can be proven that solutions of (VPLF) are given by (1) and for z large enough where N k and D k are holomorphic functions defined when z is large enough by Thus to get a solution g of (VPL) by (1) we just have to solve the equation (2) (called dispersion relation) determining an inverse Laplace transform.
Up to some strong assumptions on f eq and g 0 (k) (precised later), it can be proven that N k and D k are entire functions and that, for all λ ∈ R, the number of zeros of D k with an imaginary part larger than λ is finite (see Remark 2). Thus, using the formula and realizing precise estimates of remainder terms, we can prove that (2) has an analytic solution ψ whose analytic expansion is given, for all λ ∈ R, by where P ω,k is the polynomial such that M k (z) = z→ω L[P ω,k (t)e −iωt ](z) + O(1) is the expansion of M k (z) in ω.
Such an analysis was first realized by Landau [10], in 1946. It has been done rigorously and generalized in 1986 by Degond [7]. It gave a partial explanation to the phenomenon of Landau damping. This latter corresponds to the dynamic of (VP) when for all k ∈ T d , D k (z) does not vanish if z ≥ 0. In this case, the electric potential goes exponentially fast to zero as t goes to +∞. In 2011, Mouhot and Villani proved the existence of this phenomenon for the nonlinear Vlasov-Poisson equation (VP) in [11].
As we have just seen, due to the absence of energy exchange between the spaces modes at the linear level, the linearization is not relevant to explain really multidimensional phenomena. Furthermore, of course, it can not explain nonlinear phenomena. This motivates thus the study of the dynamic of the second order term in the expansion of f as powers of ε. More precisely, we look for a solution of (VP) under the form f = f eq + εg + ε 2 h + o(ε 2 ), where h(t = 0) ≡ 0. Neglecting the third order terms, it can be proven formally We recognize the linearized Vlasov-Poisson equations, with an initial condition equal to zero but with a source term. In that case, we refer to Denavit [8], for one of the first works on the subject, in 1965. Different second order oscillations appear and have been studied by physicists (see for example [13] and references therein; there are many references especially from the 1960s and 1970s). Our aim is to make here a rigorous mathematical study of the asymptotical behavior of the solutions of these equations, which has, to the best of our knowledge, not already been performed. This expansion in power of ε is quite natural because if we assume that f , the solution of (VP), is a smooth function of ε then g and h are just its first and second Taylor coefficients around ε = 0. We admit that there is a priori no reason that this second order approximation provides a more accurate approximation of f for long times than the usual first order approximation. Nevertheless, numerical results suggest that this second order approximation is relevant for long times both if the equilibrium is stable or not (see [4], [13] and Section 5). Furthermore, to the best of our knowledge, the asymptotic expansion of the solutions of (VPL2) we realize in this paper (see Theorem 1.4) is much more precise than what is known for the asymptotic behavior of the solution of (VP) (see [5], [11]).
Here, as for the linear case, we solve (VPL2) using a Laplace transform for the time variable and a Fourier transform for the space variable. More precisely, some calculations prove that (VPL2) is equivalent to and when z is large enough where D k is given by (3) and N k is a meromorphic function on C explicitly known. As previously, there is just to invert a Laplace transform to solve (VPL2). As for the linear case, a precise study of M k and its poles gives a solution µ of (6) and an asymptotic expansion of the form where Q ω,k is the polynomial such that The poles of M k are of two kinds: they can be zeros of D k (generating the same frequencies as at the first order) or poles of N k . The study of the poles is technical because N k is defined from the solution of (VPL). However, the asymptotic expansion of ψ (see (5)) enables a decomposition of N k in more elementary terms whose poles can be determined.
In order to give an intuition of these poles, we consider a term that is very representative 2 of this decomposition: is the Fourier transform of f eq . The later being defined for u ∈ L 1 (R d ) and ξ ∈ R d by , it can be proven that its poles are given by the asymptotic expansion of F (rep) k (t) with the formula (4). As it is suggested by the formula (8), the behavior of this later is quite different if the set of the points (τ, s) such that τ k 1 + sk 2 = 0 is a line segment (resonant case) or a point (non-resonant case).
In the non-resonant case, there exists a constant c > 0 such that So, assuming that f eq is regular enough so that F [f eq ](ξ) decreases faster than any exponential as |ξ| goes to +∞ (for example like a Gaussian), we can prove that the integral in (8) converges faster than any exponential as t goes to +∞. As a consequence, we get a constant a ∈ C such that In the resonant case, there exists γ ∈ (0, 1) such that Realizing a natural change of coordinates in (8), we get Thus, assuming that f eq is regular enough so that F [f eq ](ξ) decrease faster than any exponential as |ξ| goes to +∞, we have and we can prove that the third term decreases faster than any exponential. Thus, this decomposition provides the following asymptotic expansion where a, b ∈ C and ω b = (1 − γ)ω 1 = (|k|/|k 1 |)ω 1 is the Best frequency (according to [13]). As suggested by this sketch of proof, we can prove that M k have three kinds of poles. More precisely, if ω is a pole of M k it satisfies one of the following conditions We recall that these poles drive the asymptotic behavior of µ(k) through formula (7). The frequencies (I) and (II) have already been identified in our preliminary work on this subject [4], but not the frequency (III). We emphasize that all the three type of frequencies are listed in [13], which makes our analysis coherent with the physics litterature.
To conclude this presentation, we are going to state a precise theorem giving the asymptotic behavior of the solutions of (VPL2). To this end, we need to introduce some notations. Definition 1.1. Let E (R d ) be the subspace of the Schwartz space S (R d ), of functions f , whose Fourier transform, F f , extends to an entire function on C d and such that where | · | denotes the canonical Hermitian norm of C d .
This assumption is probably not optimal but it is crucial in our proof in order to invert easily some Laplace transforms (see , D k and N k are entire functions and for all λ ∈ R, the number of zeros of D k with an imaginary part larger than λ is finite (proof will be given in Corollary 3 and Proposition 7). Appendix 6.3 provides an algorithm to computate the zeros of D k .
Definition 1.2. If k ∈ T d , n k,ω denotes the multiplicity of ω as zero of D k , i.e.
Most of the result of this paper will require that g 0 is supported on a finite number of spatial modes whose set is denoted K ⊂ T d \ {0}. More precisely, they require the following assumption This assumption seems clearly not optimal but it is general enough to exhibit the relevant phenomena and it corresponds to the usual initial data used for numerical simulations. Furthermore, it simplifies most of the proof avoiding several problems of convergences.
We can now state the main result of this paper: the following theorem proves the existence of smooth solutions of (VPL) and (VPL2) and describes their asymptotic behavior.

Remark 3.
We have e −iωt = e (ω)t−i (ω)t , so that all the terms of the sum except maybe the remainder R are of the form e ik·x P (t)eλ t+iβt , with P (t) ∈ C[t]. If one of this term satisfiesλ < λ, it can be put in the remainder term R.

Remark 4.
Taking λ decreasing to −∞ makes the sum larger, but it always remain finite, for a fixed λ, since K, K + K are finite together with the zeros (see Remark 2). We warn the reader that, a priori, the expansion does not converge as λ goes to −∞.
Remark 5. It may exist some degenerate cases for which the four types of functions introduced in the second part of Theorem 1.4 are non distinct. In such a case, the numbers σ k1,k2 ω1,ω2 and ν k1,k2 ω1,ω2 do not vanish and we have where 1 P denotes the characteristic function of the property P .

Remark 6.
In the case where d k2 k1 = 0, which we will call resonant case, where the Best frequency, that is the term B appears, p can a priori be ≥ 1. For the J and I terms, the multiplicity can be equal to one, corresponding to m = = 0.
Remark 7. It is quite direct to extend the classical linear analysis of the Vlasov Poisson with 1 specie to the case of multi-species charged particles (see § 3.1.2 in [4]). Similarly, we expect that it would be possible to extend our second order analysis to the multi-species case. Actually, such an extension is discussed in our preliminary work ( § 5.1 in [4]) to explain more precisely the numerical results associated with a multi-species test case introduced in [2].
Remark 8. It may be interesting to try to extend this second order analysis to non-homogeneous equilibria. However such an analysis seems much more involved. Indeed, in this case, to carry out an analysis of the linearized equation similar to the analysis of the homogeneous case, it is usual to use action-angle variables (see for example [3], [9]). However, this change of coordinates generates some singularities and boundary effects leading to an algebraic decay of the electric field.
The fifth section of this paper is devoted to some numerical experiments. They principally aim at highlighting the Best's waves because most of the other phenomena associated with second order terms have been studied numerically in the proceedings [4]. Unlike the linear case, it seems that there is no elementary way to determine a priori the coefficients associated with the asymptotic expansion of µ. Indeed, they depend non trivially on the solution of (VPL) (and not only on its asymptotic expansion). Consequently, we use here least squares procedures, which permit to have a simple and quick way to find these coefficients.
There are some difficulty arising of these computations because, as we compare the solution of the second order expansion to the solution of (VP), this gives a constraint on ε and the final time that should be small enough. As we have seen, the final time should also not be too small in order to be in the asymptotic regime, and this is also true for ε (which is here put to the square, as we consider second order expansion) due to the limits imposed by machine precision.
We admit that for the numerical checking of codes, second order terms have not gained much popularity, maybe as the linear terms generally already give the main phenomena. We emphasize that we are here able to identify the contributions of the different frequencies, and thus do an effective comparison with, as already told, multidimensional and nonlinear features.
Some remarks about the notations. In order to keep proofs as readable as possible, we do some classical abuses of notation for integral transforms. For example, the Fourier transform on R d is always associated with the variable v, it means that if Outline of the work . In Section 2, we derive some integral equations (called dispersion relations) satisfied by solutions ψ, µ of (VPL) and (VPL2). Then we prove that it is enough to solve these dispersion relations to get solutions for (VPL) and (VPL2). The next two sections are devoted to the resolution of these dispersion relations and to the asymptotic expansions of their solutions: Section 3 is for the first order expansion and Section 4 is for the second order expansion. Finally in Section 5, we give some numerical results.
2. Derivation of the dispersion relations.
2.1. Dispersion relations for first and second order. In the following propositions, we give the dispersion relations, that are obtained through Fourier and Laplace transforms. Note that we have an expression for both the electric potentials ψ resp. µ and the distribution function g resp. h of the first resp. second order dispersion relations.
Proposition 1. Assume f eq ∈ E (R d ) and g 0 satisfies Assumption 1.3. Assume there exists a C ∞ function on R * + × T d , denoted ψ, and there exists λ 0 > 0 such that for z > λ 0 . If we define g by then g is a C ∞ function on R * + × T d × R d , continuous on R + × T d × R d and (g, ψ) is a solution of (VPL).
Proposition 2. Assume f eq ∈ E (R d ) and g 0 satisfies Assumption 1.3. Assume there exists a solution of (VPL) as in Proposition (1). Assume there exists a C ∞ function on R * + × T d , denoted µ, and there exists for z > λ 1 .
is a solution of (VPL2).

A general linearized Vlasov-Poisson equation.
In order to prove Propositions 1 and 2, as (VPL) and (VPL2) share the same structure, we focus on a general linearized Vlasov Poisson equation In the following proposition, we derive a general dispersion relation satisfied by u.
We first do not consider the coupling with the Poisson equation.
Furthermore, for all z ∈ C with (z) > λ 0 we have Proof. First, applying a space Fourier transform to (13), we get for Consequently, applying Duhamel formula, we get for We deduce of this last estimation, that for any fixed v ∈ R d and for any λ > λ 0 , e −λt g(t, k, v) is continuous and bounded on R + . Consequently, we can apply a Laplace transform on (16) and get for all z ∈ C such that z > λ 0 and v ∈ R d , Since z > 0, this relation can be divided by Finally we conclude this proof integrating with respect to v and applying Fubini Theorem (with the control (14)) to get for all z ∈ C with z > λ 0 If we want to get a closed equation on u, we have to use Poisson equation Formally, applying a space Fourier transform and a Laplace transform we would get Consequently, applying (15), we should get the following dispersion relation where D k is defined by (3).
is a solution of (VPLG).
Proof. By construction of g through Duhamel formula (20), g is obviously a con- Furthermore, we may verify by a straightforward calculation that g is solution of the Vlasov equation (13). Consequently, we just have to prove that g, u is solution of Poisson equation (18).
However u and g satisfy assumptions of Proposition 3, so we can apply it. Consequently, we know that if λ > λ 0 then e −λt g(t, k, v)dv is continuous and bounded and that its Laplace transform satisfies (15). But since L[ u(t, k)] is a solution of the dispersion relation (10), we deduce that for all z ∈ C such that z > λ 0 we have But it is well known that Laplace transform is injective on continuous functions with an exponential order (i.e. bounded by an exponential function), see Theorem 1.7.3 in [1]. Consequently, we have for all t > 0 Since space Fourier transform is also injective on regular functions, we have proven that u, g is a solution of Poisson equation (18).

Proof of Propositions 1 and 2.
We now apply Proposition 4 for the proof of Propositions 1 and 2.
But, we have proven in Lemma 6 that D k (z) = 0 if z is large enough. Consequently, L ψ(t, k) (z) = 0 if z is large enough. So we deduce by a classical criterion about Laplace transform (see Theorem 1.7.3 in [1]) that ψ(t, k) = 0.

JOACKIM BERNIER AND MICHEL MEHRENBERGER
We observe on (11) that g is clearly a C ∞ function on R * + × T d × R d . Finally, we just need to apply Proposition 4 to prove that (g, ψ) is a solution of (VPL).
Proof of Proposition 2. Let S be defined by Since, space Fourier transform of ψ is supported by K (see proof of Proposition 1), its space Fourier transform is supported by K + K. Furthermore, since λ 1 > 2λ 0 , we can construct, by a straightforward estimation, a continuous function In particular, this estimation proves that the right member of (12) is well defined if z ≥ λ 1 . Now, as in Proposition 1, we can first prove that space Fourier transform of µ is supported by (K + K) ∪ {0}, then observe that h is a C ∞ function and conclude that (h, µ) is a solution of (VPL) by Proposition 4.

Resolution and expansion of the linearized equation.
3.1. Introduction and statement of the result. In Proposition 1, we have proved that it is enough to solve dispersion relation (10) to get a solution (g, ψ) to linearized Vlasov-Poisson equation (VPL). So the aim of this section is to solve this dispersion relation introducing most of the theoretical tools useful in the resolution of the second order relation (12). In particular, many of them deal with analytic function defined on open sectors, denoted Σ α , with α ∈ (0, π), and defined by The result we are going to establish in this section is the following.
Proposition 5. Assume f eq ∈ E (R d ) and g 0 satisfies Assumption 1.3. For all λ ∈ R, for all k ∈ K, for all zero point ω of D k there exists a polynomial, denoted P k,ω , whose degree is strictly smaller than the multiplicity of ω, α ∈ (0, π 2 ) and there exists r k,λ an analytic and bounded function on Σ α such that the following expansion defines a solution of the dispersion relation (10) This proposition will be proven at the end of this section. First, we introduce some notations and many useful theoretical tools.
3.2. Definition of N k and theoretical tools. The right member of the dispersion relation (10) is very important in our study. We denote it N k (z). More precisely, it is an analytic function defined, when z > 0 by In the first part of this proof we study the regularity and the behavior of D k and N k . However, we need to introduce some classical results on Laplace transform. First, consider the following Theorem that is very useful to invert Laplace transforms and to control it. Theorem 3.1. (Analytic representation) Let α ∈ (0, π 2 ), λ 0 ∈ R and q : i(λ 0 , ∞) → C. The following assertions are equivalent: (i) There exists a holomorphic function f : Σ α → C such that (ii) The function q has a holomorphic extension q : Proof. See Theorem 2.6.1 in [1] page 87.
There is a direct corollary of the proof of Theorem 3.1 that is useful in our study.
Corollary 1. Assume that conclusion of Theorem 3.1 holds. Then for all 0 < γ < β < α, we have Then, we observe that D k and N k are defined through a integral operator whose kernel is 1 v·k−z . The following lemma links this operator to more classical ones.
dv is a continuous and bounded function, so its Laplace transform is well defined if z > 0. Now, consider the following function Since f ∈ L 1 (R d ), it is a regular function and we have But, since z > 0 we observe that F (t) goes to 0 when t goes to +∞. Consequently, we get 3.3. Estimations for D k and N k . With Lemma 3.2, we can write D k and N k as Laplace transforms. So, in the following proposition, we can prove their analyticity using the analytic representation theorem (Theorem 3.1). In particular, we prove and extend Remark 2.

JOACKIM BERNIER AND MICHEL MEHRENBERGER
• there exists α ∈ (0, π 2 ) such that for all 0 < γ < α and for all λ 0 ∈ R there exists C > 0 satisfying Consequently, D k is well defined as an analytic function on {z ∈ C | z > 0}. Furthermore, we can apply Lemma 3.2 to get for z > 0 Then we define e k = k |k| and we get by the change Now, using Theorem 3.1, we are going to prove this Laplace transform defines an entire function and we are going to control it. Since f eq ∈ E (R d ), it extends to an analytic function and there exists α ∈ (0, π 2 ) such that for all β ∈ (0, α) and for all λ ∈ R, there exist C > 0 such that Consequently, we get Finally, we have proven that for all λ 0 ∈ R, there exists a constant M > 0 (independent of e k ) such that Applying Theorem 3.1 and its corollary, we have proven that is an entire function and that for all γ ∈ (0, α) and all λ 0 ∈ R, there exists M > 0 (associated to β = α+γ 2 ) such that Finally, we deduce directly the result from formula (21): , D k is an entire function and for all λ ∈ R, {ω ∈ C | D k (ω) = 0 and ω ≥ λ} is a finite set.
Proof. In Proposition 6, we have proved that D k is an entire function and it can be directly deduced from its Corollary 2 that {z ∈ C | D k (z) = 0 and ω ≥ λ} is bounded. Consequently, since zero points of D k are isolated, it is a finite set.
It is very natural to adapt this result to N k . More precisely, we deduce the following proposition.
Since v → g 0 (k, v) ∈ S (R d ), N k defines naturally an analytic function for z > 0. Furthermore, applying Lemma 3.2 we know that for z > 0 is an entire function and there exists α ∈ (0, π 2 ) such that for all β ∈ (0, α) and all λ 0 ∈ R, t → |e λ0t F [ g 0 (k, v)] (kt)| is bounded on Σ β . Consequently, we deduce of Theorem 3.1, that for all λ 0 ∈ R, N k has a holomorphic extension on iλ 0 + iΣ α+ π 2 such that ∀β ∈ (0, α), sup Finally, since N k is analytic on iλ 0 +iΣ α+ π 2 for all λ 0 ∈ R, it is an entire function. 3.4. A theoretical tool for the control of L −1 [N k /D k ]. Now, we introduce a general criterion to invert Laplace transform and get an asymptotic expansion. Lemma 3.3. Let R ∈ H(C) be an entire function and N be a meromorphic function defined on C. If there exists α ∈ (0, π 2 ) such that ∃C > 0, sup then for any λ ∈ R, there exist β ∈ (0, α) and a function r ∈ H(Σ β ) analytic and bounded on Σ β such that if z is large enough then where Z is the set of poles of N (z) 1−R(z) , is the polynomial whose coefficients are defined by the expansion of Remark 10. In the application, for the proof of Proposition 5, N will be entire (thus meromorphic), but for the second order case, in the next section, we will really need that N is meromorphic.
Proof of Lemma 3.3. Many geometrical objects are going to be introduced in this proof. The reader can refer to Figure 1 to an illustration of these constructions.
First observe that to prove the lemma, we can assume that λ is negative enough. In particular we assume that λ < −2C.
Proof of Proposition 5. In Proposition 6 and 7 we have proven that we can apply Lemma 3.3 with D k = 1 − R and N = N k , taking λ 0 = 0. But the result of this lemma is exactly the expansion of Proposition 5.
Remark 11. As we use only λ 0 = 0 in Proposition 6 and 7 for the proof of Proposition 5, we may wonder of we could use a weaker assumption on f eq for getting the estimate on D k for example. Indeed, that estimate derives from Theorem 3.1 for λ 0 = 0 and so the weaker assumption could be the hypothesis of (i) in Theorem 3.1 for λ 0 = 0. However, we also need to have that D k is entire, and there we have used Theorem 3.1 for all λ 0 ∈ R.

4.
Resolution and expansion of the second order equation.

4.1.
Introduction and statement of the result. In Proposition 2, we have proven that it is enough to solve dispersion relation (12) to get a solution (h, µ) to second order linearized Vlasov-Poisson equation (VPL2). So this section is devoted to the resolution of this second order dispersion relation, following the strategy established for the first order dispersion relation in the previous section, by proving the following proposition, which permits to complete the proof of our main result, Theorem 1.4.
Proposition 8. Assume f eq ∈ E (R d ) and g 0 satisfies Assumption 1.3. Consider the solution (g, ψ) of (VPL) given by Proposition 5 and Proposition 1. Then there exists a solution µ of the dispersion relation (12) whose expansion is detailed in Theorem 1.4.

4.2.
Definition of N 1 k and N 2 k (the right hand side). We will first look for the right hand side of the second order dispersion relation, that was N k in the first order case.
Consequently, we can determine more precisely the right member of (12). However, to be rigorous we need to prove that our integrals are convergent. Indeed, as in Proposition 2, there exists λ 0 ∈ R such that for any λ > λ 0 , we can construct a continuous and integrable function d . Consequently, if z > 2λ 0 , we can apply Lemma 3.2 to prove that where we have used the change of variable τ + t ← s and the notation Furthermore, using definition of g (see (11)), we can precise F g(τ, k 2 , k(t − τ )). Indeed, we start from Consequently, we get So we have two numerators to study for this dispersion relation. On the one hand, we have with On the other hand, we have with So with these notations,the dispersion relation (8) may be written as, for all z such that z > 2λ 0 , 4.3. Estimates for N 1 k and N 2 k . We are going to apply the same strategy as for the resolution of the first order dispersion relation. It will be solve using Lemma 3.3. The denominator has been studied in Proposition 6. The following lemma describes the regularity and the behavior of the numerators N 1 k and N 2 k . Lemma 4.1. The function N 1 k , N 2 k have a meromorphic continuation and there exist λ ∈ R and β ∈ (0, π 2 ) such that sup • If there exists γ ∈ (0, 1) such k 2 = −γk 1 then the poles of N 1 k are the points ω ∈ C such that ω = ω 1 |k| |k1| where D k1 (ω 1 ) = 0 and its multiplicity is smaller than or equal to n k1,ω1 + 1. N 2 k has two kinds of poles. On the one hand, there are the points ω ∈ C such that ω = ω 1 + ω 2 where D k1 (ω 1 ) = D k2 (ω 2 ) = 0. On the other hand, there are the points ω ∈ C such that ω = ω 1 |k| |k1| where D k1 (ω 1 ) = 0. The multiplicity of a pole belonging to the two families is smaller than or equal to n k1,ω1 + n k2,ω2 + 1. Else the multiplicity of a pole of the first kind is smaller than or equal to n k1,ω1 + n k2,ω2 − 1 and the multiplicity of a pole of the second kind is smaller than or equal to n k1,ω1 + 1.
• Else N 1 k is an entire function and the poles of N 2 k are the points ω ∈ C such that ω = ω 1 + ω 2 where D k1 (ω 1 ) = D k2 (ω 2 ) = 0 and its multiplicity is smaller than or equal to n k1,ω1 + n k2,ω2 − 1. Now, we focus on proving Lemma 4.1. However, using analytic representation Theorem 3.1, it is directly deduced of the two following lemmas (Lemma 4.2 and Lemma 4.3) involving properties of F 1 k and F 2 k . Lemma 4.2. For all λ ∈ R there exist β ∈ (0, π 2 ) and an analytic and bounded function on Σ β denoted r such that • if k 2 = −γk 1 , γ ∈ (0, 1), then for all t > 0 with R ω1 a polynomial of degree smaller than or equal to n k1,ω1 , • else, for all t > 0 F 1 k (t) = e λt r(t).
Remark 12. In Lemma 4.1, we need that the inequality is true for a givenλ, in order to apply Lemma 3.3. However, applying Theorem 3.1, we deduce from Lemmae 4.2 and 4.3 that the inequality is true for allλ ∈ R. On the other hand, we have needed that Lemma 4.2 and 4.3 are true for all λ ∈ R, in order to prove that N 1 k and N 2 k are meromorphic.
We are going to prove these lemmas distinguishing the non resonant case from the resonant case (when there exists γ ∈ (0, 1) such that k 2 = −γk 1 ). In order to get proofs as clear as possible we do not prove that the remainder term can be extended on complex cones and we only control them on R * + . Indeed, there are no real issues to extend them and the arguments to control them on Σ α or R * + are the same. Furthermore, the notations induced for the complex extensions are quite heavy and so do not help to understand the ideas. However, in the first proof, to give an example, we prove the analytic extension and we really estimate it.
Indeed, in the resonant case there exists γ ∈ (0, 1) such that k 2 = −γk 1 = γ(k 2 −k), so that (1 − γ)k 2 + γk = 0. We have proven in Proposition 5 that there exists α ∈ (0, π 2 ) such that ψ(t, k 1 ) extends to an analytic function on Σ α and that there exists λ 0 ∈ R and M > 0 such that Furthermore, since v → g 0 (k 2 , v) ∈ E (R d ), its Fourier transform extends to an entire function on C d and we can assume (choosing α small enough) that for all λ 2 ∈ R there exists C λ2 > 0 such that Now observe that by a change of variable, F 1 k (t) can be written as Consequently, F 1 k (t) naturally extends to an analytic function on Σ α . Now, we have to control F 1 k (z)e −λz on Σ α for any λ ∈ R. Indeed, we have, for z ∈ Σ α , as we can assume λ 2 ≤ 0, So this quantity is bounded uniformly with respect to z ∈ Σ α if λ 2 < λ−|λ0| δ .

4.5.
Proof of Lemma 4.2 in the resonant case. As explained before, from now, we do not pay attention to the analytic extension anymore. First, we use the resonance to give a more adapted expression of F 1 making the change of variable s ← (1 − γ)t − τ . We want to expand ψ, so we introduce the dependency of F 1 k with respect to t → ψ(t, k 1 ) by denoting F 1 k [ ψ(t, k 1 )](t). Consequently, using the expansion of ψ of Proposition 5, for any λ 1 ∈ R, we get where r k1,λ1 is a bounded function on R * + . First, we are going to control the remainder term F 1 k [e λ1t r k1,λ1 (t)](t). Using the same control of the Fourier transform as previously (see (29)), we have, as we can assume λ 1 , λ 2 ≤ 0, So this quantity is bounded uniformly with respect to t ∈ R * + if (1 − γ)λ 1 < λ and λ 2 |k 1 | < λ 1 . Now, we are going to study one leading term of the type F 1 k [t n e −iω1t ](t). So, we are doing a new expansion.
where b 0 , . . . , b n+1 are real numbers. Here we recognise the leading terms of (25) since, by construction, 1 − γ = |k| |k1| . Then, observe that since the right integral is convergent (see (29)), there exists A ∈ C such that for any t ∈ R * + , we have The complex number A is the leading term of this integral whereas the other ones are remainder terms. So we just have to control them. Indeed, we have as we can assume λ ≤ 0 and since t ≤ s γ . Consequently, it is bounded uniformly with respect to t if |k 1 |γλ 2 < λ − ω 1 . The estimation of the third integral can be realized with the same ideas.
As we have a term in t n+1 , we see that R ω1 is of degree ≤ n k1,ω1 , since P k1,ω1 is of degree ≤ n k1,ω1 − 1.

4.6.
Proof of Lemma 4.3 in the non resonant case. First, operating the change of variable τ = t − τ , s = t − s, we can write F 2 k as In order to get notations general enough but compact, we denote We define, for continuous functions φ 1 , φ 2 with an exponential order, a bilinear operator q by With these notations, we have Consequently, using the expansions of ψ(t, k 1 ) and ψ(t, k 2 ) established in Proposition 5, we get 3 for λ 1 , λ 2 ∈ R, where r k1,λ1 and r k2,λ2 are respectively bounded by constants C λ1 and C λ2 . Furthermore, we can also assume that there exists λ 0 ∈ R and M > 0 such that Finally, since we are treating the non-resonant case, we may assume that there exists δ > 0 such that So first, we are going to control the remainder terms of (31). For example, we consider q[e λ1t r k1,λ1 , ψ(t, So, this quantity is bounded uniformly with respect to t > 0 if λ 1 < λ − λ 0 and λ 3 < λ1+λ0 δ < λ δ . Similarly, we could prove that if λ 2 is chosen negative enough then we could control q[ ψ(t, k 1 ), e λ2t r k2,λ2 ](t)e −λt uniformly with respect to t, and also q[e λ1t r k1,λ1 , e λ2t r k2,λ2 ]. Now, we consider a generic leading terms of (31) of the type q[t n1 e −iω1t , t n2 e −iω2t ]. So first, we can expand it where b ∈ R 0,n1 × 0,n2 are some real coefficients. 3 Realizing a decomposition of the form We observe that this last integral converge when t goes to +∞. Indeed, we have s 0 τ j1+1 s j2 e iω1τ +iω2s u(kτ + k 2 (s − τ ))dτ ≤ C λ3 s j1+j2+2 e (|ω1|+|ω2|)s e λ3δs Consequently, there exists a complex constant A ∈ C such that 0≤τ ≤s≤t τ j1+1 s j2 e iω1τ +iω2s u(kτ + k 2 (s − τ ))dsdτ This complex number A generates the term of frequency ω 1 + ω 2 in (28). So we just need to prove that the other term is a remainder term controlling it. Indeed, we have ). Concerning the degree, we see that it is ≤ n k1,ω1 − 1 + n k2,ω2 − 1, since n 1 ≤ n k1,ω1 − 1 and n 2 ≤ n k2,ω2 − 1, which corresponds to what is expected. 4.7. Proof of Lemma 4.3 in the resonant case. We consider now the last case, which is the most complex. We keep the notations of the previous subsection but we need a new expression of q adapted to the resonance: The term τ + γs is quite heavy for our estimations, so we introduce a last notation Consequently, we can expand q[φ 1 , φ 2 ] as follow We also introduce 4 a new expansion of F 2 k more adapted to the resonance Now, we are going to study each one of the terms of this expansion. Last term. First, we control the last remainder term, q[e λ1t r k1,λ1 , ψ(t, k 2 )]. Indeed, if t > 0 we have So this last quantity is bounded uniformly with respect to t if λ 1 and λ 3 are chosen negative enough. More precisely, we need (1 − γ)λ 1 < λ and λ 3 |k 1 | < λ 1 . Second term. Now, we study the behavior of the second kind of term in the expansion of F 2 k . Expanding P k1,ω1 (t−τ −γs), we can write q[P k1,ω1 e −iω1t , e λ2t r k2,λ2 ](t) as a linear combination of term of the type t j q m l+1 [e −iω1t , e λ2t r k2,λ2 ](t) and t j q m+1 l [e −iω1t , e λ2t r k2,λ2 ](t), where B λ2,p is well defined if λ 2 is negative enough (i.e. λ 2 < −γ|λ 0 |). Consequently, we get ( where C p m = m p is a binomial coefficient. Here there are three kinds of terms. The first one is one of expected leading term. The two others are remainder terms. So we have to control them. So this last quantity is finite if λ 2 is negative enough. Then we control the last kind of term. If t > 0 then So this last quantity is bounded uniformly with respect to t if λ 2 < λ − | ω 1 |γ and λ 3 |k 1 | < ω 1 ) − |λ2| γ . Of course, we could control the other remainder term (with R + ) in a similar way.
Concerning the degree, it is smaller or equal than the degree of P k1,ω1 , that is ≤ n k1,ω1 − 1, as j + m ≤ deg P k1,ω1 . This is for the moment one degree less than what is expected in the Lemma 4.3.
Remark 13. Note the term B λ2,p in (33) is not explicit, as it relies on a remainder term of the first order dispersion relation. It is worth mentioning that this term contributes to the second order expansion, and not as a remainder term.
First term. Finally we study the first kind of terms in the expansion of F 2 k . These terms are of the type q[P k1,ω1 e −iω1t , P k2,ω2 e −iω2t ]. By a straightforward calculation, as in the previous case, it can be extended as a linear combination of terms of the type t j q m l+1 [e −iω1t , t n e −iω2t ] and t j q m+1 l [e −iω1t , t n e −iω2t ] with j + l + m = deg P k1,ω1 and n ≤ deg P k2,ω2 .
In order to pursue the proof for this first kind of terms, in the following elementary lemma, we introduce a useful algebraic decomposition. It is proven in Appendix 6.2. If ω = 0 then deg Q m,n,ω = m and deg R m,n,ω = n. If ω = 0 then Q m,n,ω = 0 and deg R m,n,ω = m + n + 1.
Furthermore, using the previous constructions, we introduce Finally we just have to prove that this last integral is a remainder term. Indeed, we have So this last quantity is finite if λ 3 is negative enough.
Concerning the degree, we consider first the case γω 1 + ω 2 = 0. As B is of degree ≤ n and R m,n,γω1+ω2 is of degree ≤ n. So we get, as j ≤ deg P k1,ω1 and n ≤ deg P k2,ω2 , that Q k1,k2 ω1,ω2 is of degree ≤ n k1,ω1 −1+n k2,ω2 −1, which is the expected value. Now, as we can have a q m+1 l term, leading to Q m+1,n,γω1+ω2 which is of degree ≤ m + 1 and as m can be chosen ≤ deg P k1,ω1 , R ω1 k1,k2 is of degree ≤ n k1,ω1 , which is now the expected value. We consider finally the case γω 1 + ω 2 = 0, so that the terms e −i(ω1+ω2)t and e −iω1 |k| |k 1 | t are the same. The terms of highest degree is then R m+1,n,γω1+ω2 which is here of degree ≤ m+n+2, that is ≤ n k1,ω1 −1+n k2,ω2 −1+2.
All the values found are thus those that are expected.

Proof of Proposition 8.
Proof of Proposition 8. In Proposition 6 and 4.1 we have proven that we can apply Lemma 3.3 with 1 − R(z) = D k (z + i λ) and N (z) = N 1 k (z + i λ) + N 2 k (z + i λ), taking λ 0 = λ/|k| in Proposition 6: we get from Proposition 6 and Lemma 4.1 But the result of this lemma is that for all λ ∈ R, we have with a function r ∈ H(Σβ) analytic and bounded on Σβ, for z large enough, with someβ satisfying 0 <β < γ < β and P ω is the polynomial such that Thus, we have we get (24), which is (12). We finally have the expansion of Theorem 1.4. Concerning the multiplicity, if one pole is common to N 1 k + N 2 k and D −1 k we have to sum up the multiplicity, leading to add n k,ω1+ω2 − 1 to the range for and n k, |k| |k 1 | ω1 − 1 to the range for p. The other concerns about the multiplicity follow from Lemmae 4.2 and 4.3, and the condition k · k 1 = 0 directly follows from the factor k · k 1 in front of (22) and (23). Note also that R * + ⊂ Σ β , so that r is bounded on R * + as stated in Theorem 1.4.

5.
Numerical results. Simulations have already been performed for multi-species and multi-dimensional simulations in [4], highlighting the relevance of second order expansion. We focus here more specifically on exhibiting a case where the Best frequency, that corresponds to the terms B in Theorem 1.4, appears.
5.1. First example. We consider the one dimensional case (d = 1 and L 1 = 2π) and solve numerically (VP) with a Semi-Lagrangian scheme and an adapted 6-th order splitting [6]. 1D periodic centered Lagrange interpolation of degree 17 is used in both x and v directions and the periodic Poisson solver is solved with fast Fourier transform.
Initial condition is and σ 2 = 2 1/4 , σ 3 = √ π/2 and ε = 0.001. We take v ∈ [−v max , v max ], with v max = 10. Numerical parameters are: the number of uniform cells in x (resp. v) that are N x (resp. N v ) and the time step ∆t ∈ R * + , leading to a grid which will be referred as N x × N v × ∆t grid. The first Fourier mode E 1,num (t) of the electric field E := −∇Φ is computed from the simulation at each time step t = t n = n∆t, using a discrete Fourier transform.
We first compute the zeros of D k = D −k (see Remark 15), for |k| = 1, 2, 3 with greatest imaginary part that are ω 1,± ±2.511728081 − 0.4796966410i, The second frequency of the mode 1 is ω (2) 1,± ±3.498058625 − 2.374303389i. Such zeros can be computing with a symbolic calculus software. An example using Maple is provided in the Appendix. Here the modes that are initialized are k 1 , k 2 ∈ {±2, ±3}. The main term is for k = k 1 + k 2 = ±1, with k 1 = ∓2 and k 2 = ±3, as ω ±1 has the greatest imaginary part among the ω k1+k2 , with k 1 , k 2 ∈ {±2, ±3}. For having k 2 = −γk 1 , with γ ∈ (0, 1), we have to take k 2 = ±2 and k 1 = ∓3, so that the Best frequencies ω b,± of greatest imaginary part are defined by In order to see such term, we have to remove the main part coming from ω ±1 . The procedure is detailed as follows. From Theorem 1.4, we look here for with ω 1 = ω 1,+ or ω 1 = ω 1,− , as it leads to the same value, and similarly for ω b . We estimate z by using a least square procedure: we first define and then define z by minimizing this quantity, that is, χ 2 (z) = min y∈C χ 2 (y), which is explicitely given by as solution of with A a matrix given by its 2 columns and b a vector, all the three vectors being indexed by j that goes through all the values such that t min ≤ t j ≤ t max . Once z is found, we estimate z 1 and z 2 using again a least square procedure on the remainder: defining this timẽ z 1 and z 2 are obtained by minimizing this quantity, that is, Again the solution is explicitely given, the matrix A being here On Figure 2, we represent the time evolution of the real part of the first Fourier mode | ( E 1,num )(t)| in absolute value, together with | ( E 1,num )(t) − ze −iω1t |, that is the quantity where we have removed the main part (it is a term J in Theorem 1.4); the latter is compared to | (z 1 + t j z 2 )e −iω b tj | that corresponds to the Best term. The parameters t min , t max ,t min andt max are chosen properly so that, in the corresponding interval, the approximation is valid. Note that a too low value is not good, as the expansion is only asymptotic and we consider only one term which is the main term asymptotically. A too high value is also not good, as we have to face with the round off or numerical error and the nonlinear behavior (note that we do not solve here the second linearized equation but the full nonlinear equation). We observe a well agreement, which is even better, by refining the grid, so that we can claim that we have exhibited the Best frequency in the numerical results, which is fully coherent with the theoretical results.

5.2.
Another case where the Best frequency is almost dominant on a spatial mode. Now we consider again d = 1 (dimension 1), but we change the spatial length of the domain L 1 = 20π, and take and σ 2 = 2 1/4 , σ 3 = √ π/2 and ε = 0.001. Now the modes that are initialized are k 1 , k 2 ∈ {±1, ±0.1}. We now need to know (we already have the value of ω 1,± from the previous subsection) ω 0.1,± ±1.592755970 + 3.218848582 · 10 −52 i, ω 0.2,± ±1.621955006 − 2.569883158 · 10 −12 i, ω 0.9,± ±2.382548194 − 0.3594880484i, The second frequency of the mode 0.9 is ω (2) 0.9,± ±3.181466437−2.102684847i. The possible values of k = k 1 + k 2 are in the set {±0.2, ±0.9, ±1.1, ±2}. The first order expansion already gives a term that is not damped (the imaginary part is almost equal to zero). We also have terms on the second order expansion that are not damped (for k = ±0.2). Nevertheless, if one consider the mode k = ±0.9, one can look at | ( E 0.9,num )(t)|. From Theorem 1.4, we look thus here for an approximation of ε −2 ( E 0.9,num )(t) in the form with z = (z 1 , z 2 , z 3 , z 4 , z 5 ) ∈ C 5 , using again ω = ω ,+ or ω = ω ,− , for ∈ R, as it leads to the same result. In order to estimate z, we compute that is attained for y = z, by using the least square method as previously. Note that we add here the weight e λt , with λ = 0.48 and then we look for all the coefficients in one step. The choice of the value of λ is coherent with the fact that from Theorem 1.4, the function e λt ε −2 ( E 0.9,num )(t) − E(t, z) should be bounded. Numerical results are shown on Figure 3. We use t min = 0 and t max = 30 for the coarse grid and have increased t max to 35 for the fine grid (for the fine grid, we could even increase this value, which was not possible for the coarse grid: the results were worse, as the solution is not precise enough for the coarse grid on late times, as shown on Figure 3). For the fine grid, we could also not really increase further than around t max = 50, as we are limited, with nonlinear effects, convergence and/or machine precision; we have also preferred not to go until t max = 50, as it leads to a worser matching, since the least square procedure tends to match for values around 50, where the matching is less good. We could also change the initial time, but it has not so much impact, as it was the case for the previous subsection, since we have added here a weight function in the least square procedure. We emphasize that we can again exhibit the Best frequency and also the two other types of frequencies, which are all in the same range, for this example. In order to get this results, we note that we had to adapt he strategy concerning the least square method that was presented for the first example; this is due to the fact the several modes are in a similar range, and it was not easy to use the first procedure (used for the first example) to catch the different frequencies.
The values of z are given here for coarse and fine mesh:  Figure 3. Time evolution of if m 1 = 0 and n 1 = 0. If m 1 or m 2 = 0, we get m 1 = m 2 = 0, and similarly for n 1 and n 2 . In order to have a "real" 2D case, we can suppose that m 1 = 0 and n 1 = 0. We have −m 2 m 1 = −n 2 n 1 = 1 − γ = p q , p, q ∈ N, p < q, p ∧ q = 1.

5.3.2.
A 2D test case with Best frequency. We choose here L 1 = L 2 = L, m 1 = n 1 = 3, m 2 = n 2 = −2, so that and We will write ω , instead of ω ,+ or ω ,− , when we can either use ω ,+ or ω ,− . We will need for this subsection and the next one, the following values (note that the values are here not the same as in the one dimensional case, since the dispersion relation is not the same, as we have considered here a normalized Maxwellian): Also, the second frequency of the mode √ 2/10 is ω (2) √ 2/10,± ±0.5196579915 − 0.2520173386i.

3
(the last one is the Best frequency). In that case, we expect a similar behavior as the test case of the first subsection.
We take v 1 , v 2 ∈ [−6, 6], 32 cells in x 1 and x 2 directions, 64 cells in v 1 and v 2 directions; time step is fixed to ∆t = 0.1, leading to a 32 × 32 × 64 × 64 × 0.1 grid. The diagnostics are here obtained form the charge density ρ(t, x 1 , x 2 ) = R 2 f dv (computed from trapezoidal rule): we define ρ 1, 2,num (t) the Discrete Fourier Transform of the charge density at time t = t n = n∆t. Results are given on Figure 4. The least square procedure is here applied to minimize: min y∈C 3 tmin≤tj ≤tmax e λtj E(t j , y) − ε −2 ρ 1,1,num (t j ) 2 , with E(t, y) = y 1 e −iω √ 2/10 t + (y 2 t + y 3 )e −i 1 and is attained for y j = z j , j = 1, 2, 3, where the z j are given in Figure 4. We clearly see on Figure 4, that the Best frequency is needed: with the combination of the main frequency ω √ 2 2π L the simulated mode ρ 1,1,num is accurately asymptotically described. In that case, we see both frequencies are useful; the main frequency is not enough as we can see it on Figure 4. Indeed both modes (main and Best) are shown (they are shifted towards bottom of the Figure in order to see them better), and we see that the combination of the modes is needed to describe the simulated mode.

2π
L ,± (these frequencies were defined in the previous subsection), as we have this time and the initial data being changed to g 0 (x, v) = (cos(0.3x 1 + 0.2x 2 ) + cos(0.2x 1 + 0.1x 2 )) f eq (v), and we have still L = 20π. For the least square procedure, we consider the minimization problem min y∈C 3 tmin≤tj ≤tmax e λtj E(t j , y) − ε −2 ρ 1,1,num (t j ) 2 , with E(t, y) = y 1 e −iω √ 2/10 t + y 2 e −i(ω √ 5/10,+ +ω √ 13/10,− )t + y 3 e −i(ω √ 5/10,+ +ω √ 13/10,+ )t , attained for y j = z j , j = 1, 2, 3, where the z j are given in Figure 5. We remark here that we have some unexpected frequency at the beggining which might be interpreted as a Best frequency (the simu-first approx curve), but such one is damped and we get the right asymptotic behavior, which shows that we cannot get a Best frequency in the asymptotic limit, which is fully consistant with Theorem 1.4.
6. Appendix. Many details about these spaces can be found in [12], in particular these spaces are stable by multiplication by a polynomial or a trigonometric polynomial, derivation  The parameters λ = 0.09 and [t min , t max ] = [0, 60] are used for the least square procedure to fit (simu) by (approx) and leads to z 1 0.036159 + 0.042602i, z 2 −0.0031761 − 0.00089598i and z 3 0.010351 − 0.046355i. and the natural action of the affine group of R d . Furthermore, we obviously have F S β α (R d ) = S α β (R d ).
To get some other example, we remark that E (R d ) is clearly stable by multiplication by a trigonometric polynomial, derivation and the natural action of the affine group of R d . Furthermore, it enjoys the following tensor product property.   Proof of Lemma 4.4. If ω = 0, we get the result by expanding the polynomial (t − s) n . So we suppose now that ω = 0. Since we recognize a convolution product, we apply a Laplace transform. So we get We can apply a partial fraction decomposition to get some complex coefficients (a j ) j=0,...,n and (b j ) j=0,...,m such that n!m!(−i) n+m+2 (z + ω) m+1 z n+1 = n j=0 a j z j+1 + m j=0 b j (z + ω) j+1 .