A Hermite basis

# Stability of nonlinear Vlasov-Poisson equilibria through spectral deformation and Fourier-Hermite expansion

## Abstract

We study the stability of spatially periodic, nonlinear Vlasov-Poisson equilibria as an eigenproblem in a Fourier-Hermite basis (in the space and velocity variables, respectively) of finite dimension, . When the advection term in Vlasov equation is dominant, the convergence with of the eigenvalues is rather slow, limiting the applicability of the method. We use the method of spectral deformation introduced in [J. D. Crawford and P. D. Hislop, Ann. Phys. 189, 265 (1989)] to selectively damp the continuum of neutral modes associated with the advection term, thus accelerating convergence. We validate and benchmark the performance of our method by reproducing the kinetic dispersion relation results for linear (spatially homogeneous) equilibria. Finally, we study the stability of a periodic Bernstein-Greene-Kruskal mode with multiple phase space vortices, compare our results with numerical simulations of the Vlasov-Poisson system and show that the initial unstable equilibrium may evolve to different asymptotic states depending on the way it was perturbed.

###### pacs:
52.25.Dg, 52.35.Fp, 52.35.Sb

## I Introduction

The stability of stationary nonlinear electrostatic waves, the so-called Bernstein-Greene-Kruskal (BGK) modes introduced in Ref. Bernstein et al. (1957), is a very old and basic problem that is still of interest Paškauskas and De Ninno (2009); Khain and Friedland (2010); Manfredi and Bertrand (2000). It appears, moreover, to be closely related to the saturation of stimulated Raman scattering (SRS), an issue that motivated the present work. Indeed, for physical conditions typical of those met in inertial confinement fusion, the electron plasma wave (EPW) that grows unstable due to SRS sees its amplitude change over space and time scales much larger than the Debye length and the plasma period, respectively Bénisti et al. (2010). In a sense, this EPW is therefore close to a BGK mode. Moreover, recent one-dimensional (-D) Vlasov simulations of SRS Brunner and Valeo (2004) indicate that Raman reflectivity stopped increasing monotonically with time due to the growth of sidebands, resulting from a purely electrostatic instability similar to that introduced in Ref. Kruer et al. (1969). More recently, the so-called vortex fusion instability Ghizzo et al. (1988), which we will further detail in this paper, was invoked in Ref. Ghizzo et al. (2009) to explain why Raman reflectivity stopped growing monotonically.

In this paper, we strictly restrict to BGK equilibria, and describe a systematic and very efficient method to address their stability properties. Compared to a purely numerical approach consisting in integrating the Vlasov-Poisson system, our method not only very precisely predicts the growth rate of the instability but also the functional form of the few fastest growing modes. This allows one to illuminate the physics behind the instability, to discern different pathways for subsequent evolution depending on the mode triggered and to devise viable control strategies Sipp et al. (2010). Although many approaches Goldman (1970); Santini (1970); Lewis and Symon (1979); Ghizzo et al. (1988); Paškauskas and De Ninno (2009) have been developed to study the stability of nonlinear electrostatic waves, to the best of our knowledge none has provided such a precise description of the unstable modes by using a very general formalism and with very moderate computational cost.

In order to determine the functional dependence of the unstable modes, an eigenproblem formulation is required. It is derived by linearizing the governing equations around any equilibrium distribution function . This leads to the following general formulation

 ∂f1∂t=Af1, (1)

where is a linear operator which depends on , while is an infinitesimal perturbation. The eigenvalues of  determine the stability of the equilibrium characterized by . This equilibrium is unstable if some of the eigenvalues of have a strictly positive real part, the largest of which being the growth rate of the instability.

For spatially homogeneous equilibria, the eigenproblem (1) has been treated by Van Kampen Van Kampen (1955) and Case Case (1959). For spatially inhomogeneous BGK equilibria, characterized by a non-vanishing electric field, we propose here an approximate resolution of the eigenproblem (1) by making use of the Galerkin spectral method Boyd (2001). Hence, we expand the total distribution function over a finite set of global smooth functions which fulfill the boundary conditions of our problem, and which are moreover chosen to be orthogonal. Then, the operator  is approximated by a finite dimensional matrix.

In many areas of physics, most notably fluid mechanics, spectral methods have been particularly successful in solving eigenproblems of the form (1), with exponentially fast convergence of the result as a function of the number of orthogonal functions retained in the expansion for a sufficiently smooth (see Ref. Boyd (2001), for example). This hallmark of spectral methods makes them attractive for the study of the stability of Vlasov-Poisson and Vlasov-Maxwell equilibria, especially in dimensions higher than one. However, a convincing application of spectral methods to the problem of stability of BGK waves is still, to the best of our knowledge, lacking. Problems arise owing to the very nature of the operator . In particular the linear, advection term in Vlasov equation [see Eq. (3a) below] is responsible for the transfer of energy to very fine velocity scales, a phenomenon known as velocity space filamentation (see, for example, Ref. Manfredi (1997)). When the contribution of the advection term is significant, eigenfunctions of  are expected to involve high order velocity modes (fine velocity scales) leading to slow, power law convergence of the expansions on the orthogonal functions rather than the usual exponential convergence. An extreme example is the continuum of singular, neutral van Kampen modes whose approximation by means of smooth functions is clearly out of reach. Although we will not be interested in the continuum of neutral modes, it still has to be accounted for in the numerical approximation of  and can interfere with the determination of the unstable modes of interest, particularly the weakly unstable ones.

To ensure fast convergence of the eigenvalue calculation, it is clear that the role of the neutral modes has to be undermined. To this end, we use the method of spectral deformation, originally developed for quantum mechanical problems and introduced in the study of the Vlasov-Poisson system by Crawford and Hislop Crawford and Hislop (1989); Hislop and Crawford (1989). The method introduces an operator

 B(θ)=U(θ)AU−1(θ), (2)

with non-unitary for complex . It can be proved Hislop and Crawford (1989) that the eigenvalues of  with nonzero real part remain unchanged under suitably chosen transformations , while the continuum of neutral modes is damped. Our central observation is that if the corresponding damping rate is chosen so that the continuum spectrum is well separated from the discrete eigenvalues of interest, then exponential convergence with the truncation order can be recovered.

We choose to expand the distribution function in Fourier series (in space) and Hermite functions (in velocity). This decomposition was first introduced by Grant and Feix Grant and Feix (1967) in the study of stability of spatially homogeneous Vlasov-Poisson equilibria and recently generalized to the full Vlasov-Maxwell system and inhomogeneous equilibria by Camporeale et al. Camporeale et al. (2006). Even more recently Paškauskas and De Ninno Paškauskas and De Ninno (2009) revisited the method from a nonlinear dynamics perspective, studying homogeneous equilibria of the Hamiltonian mean field model. We show through numerical examples that, with the introduction of spectral deformation, Fourier-Hermite expansions converge fast enough to be useful for the determination of unstable modes, even for strongly inhomogeneous equilibria.

The organization of this paper is as follows. We introduce the Vlasov-Poisson system and its linear approximation around an inhomogeneous equilibrium in Sec. II. For simplicity, we restrict to one space and one velocity dimension. Spectral deformation is introduced and applied to the Vlasov-Poisson system in Sec. III. In Sec. IV we derive the representation of  in the Fourier-Hermite basis. We illustrate the treatment of Landau damped modes by our method in Sec. V.1.1. In Sec. V.1.2 we compare our results for a spatially homogeneous, bump-on-tail, distribution function against those obtained by numerically solving Landau’s dispersion relation. We also use this test problem to illustrate some of the convergence issues that can arise as the wavelength of perturbations decreases and the advection term becomes more significant, and their resolution through spectral deformation. In Sec. V.2 we study a nonlinear example, namely a BGK mode with multiple phase-space vortices. Our results for the growth rate and the unstable modes agree with numerical simulations of the (nonlinear) Vlasov-Poisson system. At a qualitative level our calculations show that the collective modes trigger the vortex fusion observed in numerical simulations. We discuss our findings and the potential for optimization, as well as our future studies based on this work, in Sec. VI. In the appendices we provide some technical details on the properties of the Hermite basis used here (Appendix A), the representation of in Fourier-Hermite basis (Appendix B), the truncation of the representation of  (Appendix C) and the calculation of the Fourier-Hermite coefficients (Appendix D).

## Ii Linearization of the Vlasov-Poisson system

We consider the one-dimensional motion of electrons in an unmagnetized plasma, with immobile ions forming a neutralizing background. We restrict to situations that can be modeled using periodic boundary conditions in a spatial domain . The motion is described in terms of the Vlasov-Poisson system

 ∂f∂t+v∂f∂x+E∂f∂v=0, (3a) ∂E∂x=(∫+∞−∞fdv−1), (3b) ∫L0dxE=0, (3c)

where velocity is normalized to the thermal one , space to the Debye length , time to the inverse of the electron plasma frequency and electric field to , where is the electron charge.

Let be an equilibrium solution of (3a),

 v∂f0∂x+E0∂f0∂v=0. (4)

The Vlasov equation is Galilean invariant and therefore any traveling wave solution can be reduced to an equilibrium solution without loss of generality.

Substituting in (3), where and are infinitesimal perturbations, and accounting only for first order terms in , we get

 ∂f1∂t=−v∂f1∂x−(E0∂f1∂v+E1∂f0∂v), (5a) ∂E1∂x=∫+∞−∞f1dv, (5b) ∫L0dxE1=0, (5c)

which describe the evolution of in the linear neighborhood of .

Owing to periodic boundary conditions in the space variable, the distribution function and electric field may be expressed in terms of Fourier series,

 f0(x,v)=+∞∑r=−∞fr0(v)Φr(x), (6a) E0(x)=+∞∑r=−∞Er0Φr(x), (6b)

where , , and similar expansions hold for and . Plugging system (6) into system (5) and using the standard Fourier basis orthogonality relations, we get

 ∂∂tfk1(v,t)=−ikk0vfk1(v,t)+ik0∞∑∑′r=−∞1r∂∂vfk−r0(v)∫+∞−∞dvfr1(v,t)+ik0∞∑∑′r=−∞1r∂∂vfk−r1(v,t)∫+∞−∞dvfr0(v), (7a) Ek1(t)={0,if k=0,−ikk0∫+∞−∞dvfk1(v,t),if k≠0, (7b)

where the prime in summations indicates that we omit the term, as we have incorporated (7b) into (7a) to eliminate the electric field. The restriction follows from condition (5c) on the electric field. Equation (7a) is of the form (1),

 ∂f1∂t=Af1,

where  is a linear integro-differential operator that depends on . For the rest of this paper we will use the distribution function alone to refer to solutions, keeping in mind that the electric field is determined self-consistently through Poisson’s equation.

Owing to the Hamiltonian structure of the Vlasov-Poisson system Morrison (1980), eigenvalues of the real operator  come into quartets [see Fig. 1(a)]. Moreover, characteristically has a continuum spectrum on the imaginary axis.

If for some , a perturbation in the ’th eigenspace grows in modulus as , while its phase oscillates at frequency . Then, the norm of a generic perturbation having non-zero components along all eigendirections would grow asymptotically in time as . In this case, we will say that the equilibrium is unstable.

If no eigenvalue with strictly positive real part exists, for all , the eigenvalues form a continuoum which coincides with the imaginary axis and the corresponding eigenmodes are undamped (neutrally stable). These modes are singular (described by generalized functions or distributions) and do not represent physically observable modes of the system. However, their presence is connected to the collisionless damping of generic electric field perturbations. Solving the initial value problem for small amplitude electrostatic waves in a Maxwellian plasma, Landau Landau (1946) famously showed that the electric field amplitude vanishes at an exponential rate. As shown by Van Kampen Van Kampen (1955) and generalized to more general spatially homogeneous equilibria by Case Case (1959), Landau damping can be understood as a destructive interference effect (known as phase mixing) of the neutral modes.

In this paper, we will adopt the convention of sorting our eigenvalues by decreasing real part, so that for unstable equilibria corresponds to the largest growth rate and the eigenfunction is the fastest growing mode. Moreover, we will not use as is traditionally the case in plasma physics when writing solutions of (5). Instead, we conform to the convention usually employed in the study of linear ordinary differential equations, which is more natural to a spectral discretization of the eigenproblem for . As a result, the continuous part of the spectrum, , lies on the imaginary axis, as illustrated in Fig. 1(a).

## Iii Spectral deformation

As already stated in the introduction, in order to alleviate the difficulty due to the continuum of neutral modes in the spectrum of , we do not address the stability of a BGK equilibrium by solving the eigenproblem for , but for the transformed operator

 B(θ)=U(θ)AU−1(θ), (8)

with

 U(θ)hk(v)=hk(v+θk). (9)

Here,

 θk={sgn(k)θ,k≠0,0,k=0,

and is the th Fourier mode of a given function . As shown by Hislop and Crawford in Refs. Crawford and Hislop (1989); Hislop and Crawford (1989), the main merit of this transformation is that for the continuous spectrum becomes damped, while the eigenvalues with remain unchanged [see Fig. 1(b) and Ref. Crawford and Hislop (1989) for a homogenous equilibrium]. When eigenvalues with of multiplicity higher than one exist, we can talk of discrete eigenvalues “embedded” in the continuous spectrum (see Ref. Crawford and Hislop (1989)). These embedded discrete eigenvalues also remain unchanged under spectral deformation. Hence, the stability issue of a BGK equilibrium may be equivalently addressed by calculating the eigenvalues of  or of , except that these are much more easily and accurately estimated for , since the damped continuous spectrum is much easier to account for numerically.

One may however wonder how Landau damping, which results from the phase mixing of the neutral modes, could be recovered using spectral deformation. As shown in Ref. Crawford and Hislop (1989), for a homogeneous equilibrium, when the neutral modes of become damped by . When this is larger in absolute value than the damping rate obtained through Landau’s dielectric tensor formalism, a new pair of discrete eigenvalues with [and equal to the frequency predicted by Landau] appears in the spectrum of  [see Fig. 1(b)]. The complex conjugate eigenvalue corresponds to wavelength . Hence, Landau damping is indeed recovered but, unlike for the original Vlasov-Poisson system, it now appears as being due to the damping of an eigenmode (corresponding to the largest negative eigenvalue) for the dissipative dynamics represented by the operator . This will be discussed in more detail in Sec. V.1.1.

Let us now, as for the original Vlasov-Poisson system, write the eigenproblem for  in Fourier space. Using, (7a) and , one easily finds that the equation

 ∂g1∂t=B(θ)g1 (10)

becomes, in Fourier space,

 ∂∂tgk1(v,t)=−ikk0(v+θk)gk1(v,t)+ik0∞∑∑′r=−∞1r∂∂vgk−r0(v+θk−θk−r)∫+∞−∞dvgr1(v,t)+ik0∞∑∑′r=−∞1rU(θk−θk−r)∂∂vgk−r1(v,t)∫+∞−∞dvgr0(v). (11)

In (11) the introduction of dissipation through spectral deformation with becomes apparent through the presence of in the advection term. The invariance of the discrete spectrum with is not obvious from (11) but has been established in Ref. Hislop and Crawford (1989). The unstable eigenmodes of the initial Vlasov-Poisson problem can be recovered through the inverse transformation, , since and is time independent.

## Iv Hermite expansion

In order to approximately solve the eigenproblem for , we now write as a sum of orthonormal Hermite functions. Such expansions have been first introduced in numerical studies of the Vlasov-Poisson system by Grant and Feix Grant and Feix (1967) and Armstrong Armstrong (1967) and present various advantages. Hermite functions decay as at large and therefore allow one to treat boundary conditions correctly without truncating the infinite interval. Moreover, convenient three term relations (see Appendix A) of the Hermite functions will result in an explicit, sparse representation of the operator . Hermite functions are related to velocity moments, the few first of which are directly linked to physical observables Grant and Feix (1967); Camporeale et al. (2006). As pointed out by Paškauskas and De Ninno Paškauskas and De Ninno (2009), this allows one, at least in principle, to naturally separate thermal and filamentation scale effects. As we will see in our numerical examples of Sec. V, filamentation scale often strikes back, rendering Hermite expansions problematic, a shortcoming pointed out a long time ago Grant and Feix (1967) and overcome here by the introduction of spectral deformation.

Here we consider the so-called asymmetrically weighted Hermite basis Schumer and Holloway (1998). Denoting by the basis functions and by the weight functions we have

 Ψn(v)=Cne−v2Hn(v),Ψn(v)=CnHn(v), (12)

where are Hermite polynomials and . More details on the properties of Hermite polynomials used here can be found in Appendix A. We note the important orthonormality relation

 ∫+∞−∞Ψm(v)Ψn(v)dv=δmn. (13)

Our expansion of over the Hermite functions reads

 gk(v,t)=+∞∑s=0gks(t)Ψs(u), (14)

where , with an arbitrary velocity scale factor whose importance will be discussed at the end of this section. Plugging Eq. (14) into (11), multiplying by , integrating over and using (31), (34) and the orthonormality of the Hermite basis we get, when ,

 ddtgkj1 =−ikk0[vs(√j2gk,j−11+√j+12gk,j+11)+θkgkj1] (15) −iπ1/4k0∞∑∑′r=−∞√2jr¯¯¯gk−r,j−10gr01−iπ1/4k0∞∑∑′r=−∞gr00rj−1∑n=0Kj,n+1(θk−θk−r)gk−r,n1, (16)

while

 dgk01dt=−ikk0[vs√2gk11+θkgk01]. (17)

In Eq. (15) we have introduced the notations

 Kjn(y)=⎧⎪ ⎪⎨⎪ ⎪⎩(−2)j−n√2nCjCn(j)n(yvs)j−n,if jn, (18)

and .

We can write (15) in tensorial form as

 ddtgkj1=∑l,mBkjlmglm1, (19)

where,

 Bk0lm=−ikk0δkl[vs√2δ1m+θkδ0m], (20)

while, when ,

 Bkjlm≡Dkjlm+Fkjlm+Gkjlm, (21)

and we have introduced

 Dkjlm =−ikk0δkl[vs(√j2δj,m+1+√j+12δj,m−1)+θkδjm], (22) Fkjlm ={−iπ1/4k0√2jlδ0m¯¯¯gk−l,j−10,if l≠0,0,if l=0\,, (23) Gkjlm =⎧⎨⎩−iπ1/4k0gk−l,00k−lKj,m+1(θk−θl),if l≠k and % m≤j−1,0,otherwise. (24)

Equations (2024) provide the representation of the linear operator in the Fourier-Hermite basis. In practice, we find it more convenient and efficient to compute the Fourier-Hermite coefficients of rather than that of . Then, and we can compute as described in Appendix B.

In computations, has to be truncated to finite order by setting for some cutoff values and , see Appendix C. Physically, such a truncation holds provided that and can be well described by functions that do not oscillate too rapidly with space and velocity. The truncated matrix, shown in Fig. 2(a), is sparse and the computation of its first few eigenvalues with the largest real parts can be efficiently handled by iterative schemes, such as Arnoldi iteration Trefethen and Bau (1997). For , contains even fewer non-zero elements, see Appendix C and Fig. 2(b).

The velocity factor reflects the freedom to rescale the infinite domain in which Hermite functions are defined. This freedom has been exploited both in numerical simulations Schumer and Holloway (1998) and eigenproblem calculations Camporeale et al. (2006); Paškauskas and De Ninno (2009). Paškauskas and De Ninno use to minimize the number of modes needed for a good approximation of the unperturbed distribution function . Camporeale et al., on the other hand, optimize (and also allow for a shift of origin in ) iteratively so as to reduce the quadratic error in the eigenvectors , given an initial approximation computed with non-optimal . We have found that for the problems of interest to us, i.e. strongly inhomogeneous equilibria, varying can be used to ensure fast convergence of the expansion of . However, in cases of power-law convergence of the eigenvalue computation, we have found that tuning is of limited help as it would not yield exponentialy fast convergence. Therefore, in our numerical examples we will fix to a value that provides fast convergence of the expansion of and resort to spectral deformation to ensure exponential convergence rate of the eigenvalue calculation.

## V Numerical Results

### v.1 Comparison with dielelectric tensor formalism

For spatially homogeneous unperturbed distribution functions, the equilibrium electric field vanishes. The problem decouples in Fourier space and the growth (or damping) rate and frequency vary continuously with and can be computed by Landau’s dielelectric tensor formalism Landau (1946). We will test the validity of our expansions by comparing our results against the predictions of Landau’s theory for two examples. In the first example, is a Maxwellian, which lets us discuss the qualitative difference in addressing Landau damping when diagonalizing versus . In the second example, we extensively compare the rate of convergence and the accuracy obtained with and without spectral deformation for an unstable “bump-on-tail” distribution.

#### Landau Damping

Since we approximate the linear operator  by a finite dimensional matrix, the continuous, neutral spectrum is represented by a discrete set of eigenvalues, the separation of which decreases as we increase . As shown by Grant and Feix Grant and Feix (1967), Landau damping is then only approximately recovered when solving the initial value problem for the electric field. Indeed, if we consider a Maxwellian and an initial small amplitude electric field perturbation of the form , which is not an eigenmode of the problem, then the solution of the initial value problem in the linear approximation yields for the first Fourier component, , of the electric field Grant and Feix (1967)

 E11(t)=E11(0)Nv∑j=0(a−1)j0a0jeiωjt, (25)

where the ’s are the eigenvalues of [given by (38)], the matrix of its column eigenvectors, and the matrix of row vectors of the dual basis. Hence, Landau damping of the EPW occurs as the consequence of the destructive interference effect of the neutral eigenmodes, but only over a finite time, because the sum (25) is finite. This is illustrated in Figure 3(a) plotting the time evolution of the EPW amplitude as predicted by (25), for . The electrostatic wave is indeed initially damped and the damping rate and frequency of oscillation can be determined from the slope and local maxima distance in Fig. 3(a), , . As shown in Fig. 3(b) and (c), agreement with the results from Landau’s analysis is excellent. However, after the electric field amplitude grows again to finite magnitude, a fact known as recurrence. The recurrence time increases as we decrease the scale factor or we increase (see Ref. Schumer and Holloway (1998)) indicating that recurrence is ultimately related to the approximation of the continuous spectrum by a finite set of eigenvalues (see also Ref. Knorr and Shoucri (1974)).

With spectral deformation, , the physics of the eigenproblem becomes dissipative: the Landau damped mode appears in the spectrum as a true eigenmode with , accompanied by and recurrence is absent. At the same time some of the collisionless physics is sacrificed, namely phenomena that depend on the reversible structure of Vlasov equation and the creation of fine velocity scales in the distribution function, such as plasma echos (See Ref. Chen (1995), for example). We do not plot the results obtained with spectral deformation in Fig. 3(b) and (c) since, within the resolution of this plot, they perfectly overlap with the results.

#### Bump-on-tail instability

Let us now compare the results of our method against those obtained from the numerical resolution of Landau’s dielectric tensor formalism for a bump-on-tail distribution of the form

 f(v)=np√2πe−v2/2+nb√2πe−(v−vb)2/2, (26)

with and , see Fig. 4(a).

The probability distribution function (26) can be effectively approximated with Hermite functions, as indicated by the exponential decay of the magnitude of Hermite coefficients in Fig. 4(a), which fall bellow roundoff for (with ). Note that both axes are logarithmic in the inset of Fig. 4(a) and exponential convergence corresponds to an envelope of the with an ever-increasing negative slope. For , the residual in is of the order of . In Fig. 4(b) we compare our results for to numerical values obtained through Landau’s analysis. Without spectral deformation, and and , agreement is good for small wavelengths, but deteriorates as is increased. In particular, the mere addition of one Hermite term in the series changes the results significantly, and for odd Grant and Feix (1967) the unstable range is narrower than predicted by Landau’s analysis. Beyond the cutoff wavelength, the damping rate would have to be determined as in Sec. V.1.1, but we will not get into this, as we are primarily interested in growing modes.

To study how the accuracy of our calculations is affected by we fix and vary the number of Hermite polynomials up to , see Fig. 4(c). As a measure of convergence we plot in Fig. 5(a) , i.e. the relative change in with the addition of a new term in the Hermite series, while in Fig. 5(b) we compare with the exact value of derived from Landau’s analysis. Without spectral deformation we observe slow, power-law convergence. Beyond the expansion of only adds numerical noise to the eigenvalue computation; however the relative change in eigenvalues is of the order of , with pronounced oscillations with odd and even order truncation in , see also Fig. 4(c).

This even-odd order oscillation behavior has been observed by other authors in similar Grant and Feix (1967) or more general settings Camporeale et al. (2006). While some (see Ref. Paškauskas and De Ninno (2009)) overcome such difficulties by averaging the eigenvalues computed over different values of , such an approach does not warrant convergence and cannot justify the choice of an expansion over Hermite polynomials, as opposed to an estimation with a low order method, for instance with finite differences.

Equation (7) shows that the relative importance of the advection term, responsible for the poor convergence of the method, increases with . For relatively large ’s, spectral deformation is called for, and with and , agreement with Landau’s analysis is recovered for all ’s in Fig. 4(b). The unstable range in is now accuratelly retrieved, and the stability threshold is crossed smoothly, since damped modes are represented as true eigemodes of  rather than through the interference of neutral modes.

A closer look at the convergence rate for different values of in Fig. 5(a) reveals that for small values of the convergence still obeys a power-law, yet steeper than the one for . For large enough values of , convergence becomes exponential and there appears to be no practical advantage in further increase of , since the convergence rate remains practically the same. From the results plotted in Fig. 4(c), we could say that the eigenvalue computation follows an overdamped oscillation pattern, with very fast relaxation towards the exact value. As the expansion coefficients fall below roundoff at , further precision gain with an increase in ceases, see Fig. 5(a) and (b).

It is interesting at this point to examine how rapidly the expansion coefficients of the computed eigenfunction fall off. With either or , coefficients fall off as a power-low, rather than exponentially as for , see Fig. 5(c). This should have been anticipated: owing to the action of the advection term the eigenfunctions span all velocity scales, up to the filamention scale. We therefore need non-vanishing contributions from higher order Hermite functions to capture such thin scale effects in the eigenfunction. The payback of using spectral deformation is a (relatively) accurate computation of high spectral coefficients. On the contrary, with higher spectral coefficients are subject to strong even-odd order oscillations, see Fig. 5(c).

### v.2 Stability of BGK waves

Bernstein-Greene-Kruskal (BGK) modes Bernstein et al. (1957) are stationary, nonlinear electrostatic waves accounting for the presence of both trapped and untrapped particles. They have attracted a lot of interest because of their resemblance with the saturated state of nonlinear processes in plasmas. We will study the stability of a particular BGK equilibrium introduced by Ghizzo et al. Ghizzo et al. (1988) to model processes that involve vortex fusion. Although no claim can be made regarding the generality of the shape of prescribed in Ref. Ghizzo et al. (1988), it serves well to demonstrate how our stability calculation can be used to predict vortex fusion.

The chosen equilibrium distribution function is

 f0(H)=μ√2π2−2ξ3−2ξ(1+H1−ξ)e−H, (27)

where is the total energy, a parameter that controls inhomogeneity ( corresponds to the homogeneous case), and a parameter that controls the depth of the distribution function’s “depression” or “hole” at (see Ref. Ghizzo et al. (1988) for details). Using Poisson equation, one then finds that the potential solves

 Φ′′(x)=−μ3−2ξ+2Φ(x)3−2ξe−Φ(x)+1. (28)

In order to allow for subharmonic perturbations, Ghizzo et al. study in Ref. Ghizzo et al. (1988) the stability of the equilibrium (27) for a physical system whose periodicity is times the period, , of the BGK mode. This system may be viewed as an “-cell replica” of the basic cell. Using the marginal stability analysis developped in Santini (1970), Ghizzo et al. show that when the equilibrium is unstable. This theoretical prediction is then tested against long-time numerical integrations of the Vlasov-Poisson system, with an initial state consisting of the equilibrium (27) perturbed by a small amplitude monochromatic wave whose wavelength is . For , the numerical results suggest that the equilibrium is stable. For , instability is clearly demonstrated in Ref. Ghizzo et al. (1988), as the perturbed equilibrium evolves towards a different final state through hole merging.

Here, we compute the unstable modes and their growth rates for the BGK equilibrium (27) using our Galerkin projection method. Unless otherwise noted, for all Fourier-Hermite calculations that follow, we use points/cell, , in the range , and or . Since growth rates are not computed in Ref. Ghizzo et al. (1988), we compare our results with the growth rates obtained numerically from the resolution of the Vlasov-Possion system with the 1D Eulerian Vlasov code VADOR Filbet et al. (2001).

Following Ref. Ghizzo et al. (1988), we solve Eq. (28) numerically for . The spatial period of the solution is fixed to by specifying the initial condition . Phase space portraits of the distribution function are characterized by one vortex or “hole” for each spatial period, as plotted in Fig. 6. Fourier-Hermite expansion coefficients of fall off exponentially, see Fig. 7(a).

For , the Fourier-Hermite method converges rapidly to the growth rate when we use spectral deformation with , with the relative error falling bellow for , see Fig. 8. The computation on the other hand suffers from large amplitude oscillations between even and odd order in and the convergence rate is slow.

The unstable eigenfunction is plotted in Fig. 6(b). In Fig. 6(c) we present the predicted linear evolution in the neighborhood of the BGK mode, showing that the instability causes the holes to approach each other. This triggers the eventual merging of the two holes, illustrated in Fig. 6(d), under the nonlinear dynamics.

The distribution function corresponding to the unstable mode determined by using the Fourier-Hermite expansion with and for and is shown in Fig. 9(a). We can see that the even/odd oscillations in Fig. 8(a) translate into rather large deviations in the shape of the eigenfunctions.

Our results for and the unstable mode, denoted as , are compared against those obtained through Vlasov simulations. We initialize the code VADOR with a perturbed distribution function of the form , with , and follow the evolution of the system in the neighborhood of until assumes a constant shape, plotted in Fig. 9(b), and only changes in norm. We estimate a growth rate from the rate of change of the distance . Note that the resolution used in VADOR, , , is much higher than what we needed in the Galerkin method, while the precision is lower, due to the error in graphically estimating the growth rate of the perturbations. We compare the distribution functions corresponding to the predicted unstable mode with , to the one determined by numerical integration of a small sinusoidal perturbation with code VADOR in Fig. 9(b). We observe that the two profiles agree rather well, except for some small-scale structure not captured by our method. However, the resolution is too different in the two methods for a direct comparison at this scale to be meaningfull.

In Fig. 7(b), we observe that the high-order Hermite function components of the unstable eigenmode do not fall off rapidly, indicating that the eigenfunctions we approximate involve fine velocity scales. The utility of the combination of spectral deformation and Galerkin projection is that it allows the accurate representation of the thermal scale effects of the instability, independently of filamentation scale effects. By contrast, not using spectral deformation () results in a much coarser approximation of the eigenmodes, with large fluctuations between odd and even , see Fig. 9(a). However, both and calculations provide very accurate approximation of the electric field, Fig. 9(c), as the differences in the distribution function are smoothed out when one considers its lower moments.

The -cell system has a slightly smaller positive eigenvalue , with an eigenmode that leads to two-hole fusion, as shown in Fig. 10. The numerical simulations in Ref. Ghizzo et al. (1988) show that this is indeed the case, with the third hole subsequently merged with the other two, leading to an asymptotic one-hole state. We will not present here any detailed convergence study or comparisons with numerical simulations, as the results are qualitatively similar to the -cell case.

For the -cell state a new possibility arises, since we can think of it as two -cell systems stacked together. Thus, the -cell, , unstable mode, denoted as , is still present, with eigenvalue . There is yet another unstable mode, , with smaller eigenvalue , which only exists for perturbations of wavelength . The fastest growing mode, is expected to prevail in the case of a broadband initial fluctuation spectrum. However it is possible to excite each mode independently, by imposing perturbations of appropriate wavelength. Excitation of triggers binary fusion of neighboring holes, leading to a final state with two holes shown in Fig. 11(a)-(c) which, owing to the periodicity of the system, is just the same final state as that obtained when . By contrast, excitation of leads to a more asymmetric hole-fusion senario, with subsequent nonlinear evolution leading to a one-hole state, see Fig. 11(d)-(e).

## Vi Conclusions

We have shown that the combination of spectral deformation and Fourier-Hermite expansion can be used to compute the stability of nonlinear Vlasov-Poisson waves in an efficient and systematic manner. The computation of unstable eigenfunctions for the BGK mode of Sec. V.2 illuminated the role of linear instability in triggering the subsequent vortex fusion. The ability to detect sub-dominant unstable modes permits a direct assesment of the role of perturbations of different wavelength in the evolution of the system. Our method will be used in a future paper to test the relevance of considering BGK equilibria in order to address Raman saturation.

Spectral deformation was introduced to handle filamentation scale effects through damping of the advective term in Vlasov equation. This was the key to achieving exponentially fast convergence of the eigenvalue computation and justify the choice of spectral methods. In contrast to estimates of the growth rate based on numerical integration with a Vlasov-Poisson solver, the combination of spectral deformation and Galerkin projection scales well with the dimension of phase space as the exponentially convergent series can be truncated very early, while still providing reliable results. We therefore find our method to be more suitable for extention to stability calculations in dimensions higher than -D, than an approach based on direct integration with a Vlasov code. Moreover, there is in principle no obstacle to generalizing the present method to the Vlasov-Maxwell system and to different geometries, as the choice of basis and transformed operator  can be adjusted to the problem at hand. Although our emphasis is on manipulating the structure of the linear operator itself, rather than on a careful choice of basis, we expect that the method can be further optimized, if needed, through the techniques introduced in Ref. Camporeale et al. (2006).

## Appendix A Hermite basis

The Hermite polynomials we use here follow the standarization for large , see Ref. Abramowitz and Stegun (1964). Their explicit definition is

 Hs(u)=(−1)seu2dsduse−u2. (29a) They satisfy the orthogonality condition ∫+∞−∞Hm(u)Hn(u)e−u2du=δmn√π2nn!, (29b)

and therefore

 ∫+∞−∞Ψm(u)Ψn(u)du=δmn. (30)

We note the following relations Abramowitz and Stegun (1964):

 H′n(u) =2nHn−1(u), (31a) Hn+1(u) =2uHn(u)−2nHn−1(u). (31b)

Using (31) we can show that the orthonormalized Fourier-Hermite basis (12) satisfies

 Ψ′n(u)=−√2(n+1)Ψn+1(u), (32a) uΨn(u)=√n+12Ψn+1(u)+√n2Ψn−1(u). (32b)

Using the identity

 Hj(u+c)=j∑k=0(jk)Hk(u)(2c)j−k, (33)

we can show

 ∫+∞−∞duΨj(u)Ψn(u+c)=⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩0,jn. (34)

## Appendix B Representation of translation operator

The velocity translation operator is represented in Fourier-Hermite basis by the matrix elements

 Ukjmn(θ)=∫L/2−L/2dx∫+∞−∞duΦk(x)Ψj(u)U(θm)Φm(x)Ψn(u). (35)

A direct calculation gives, with the help of (34),

 Ukjmn(θ)=δkm∫+∞−∞duΨj(u)Ψn(u+θm/vs)=δkm⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩0,jn. (36)

However, in our computations we only need the action of on a vector representing a function . We can then exploit the sparsity of the representation of the generator of velocity translations ,

 (∂v)kjlm=−v−1s√2(m+1)δklδj,m+1,

and employ Krylov subspace approximations Trefethen and Bau (1997) to , to compute in a stable and efficient manner.

## Appendix C Truncation

In practice the infinite ladder of equations (15) has to be truncated by setting for some cutoff values and . We present the truncated form of (21) to emphasize the treatment of boundary terms in Fourier space and to introduce the notations required to discuss the computation of the expansion coefficients in Appendix D. For clarity, we only present the case here; extension to the spectrally deformed case is straightforward.

 f(xj,vm,t)=Nx/2∑r=−Nx/2+1Nv∑s=0frsΦr(xj)Ψs(um), (37a) E(xj,t)=Nx/2∑r=−Nx/2+1ErΦr(xj), (37b)

and , and . The calculation of coefficients is presented in Appendix D.

Note that when taking odd derivatives of (37a) with respect to reality of the result is not ensured, since the term is not included in the sum. Thus, we have to set the (presumably small) term in the discrete Fourier transform of odd derivatives equal to zero, see Ref. (Trefethen, 2000, Chapter 3). Then, proceeding as in Sec. IV with we get, when ,

 Akjlm=Kkjlm+Lkjlm+Mkjlm, (38)

while

 Ak0lm≡{−ikk0δklvs√2δ1m,k≠Nx/2,0,k=Nx/2, (39)

and we have introduced

 Kkjlm =⎧⎪⎨⎪⎩−ikk0vsδkl(√j2δj,m+1+√j+12δj,m−1),k≠Nx/2,0,k