# 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

(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

(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

(3a) | |||

(3b) | |||

(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),

(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

(5a) | |||

(5b) | |||

(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,

(6a) | |||

(6b) |

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

(7a) | |||

(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),

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

(8) |

with

(9) |

Here,

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

(10) |

becomes, in Fourier space,

(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

(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

(13) |

Our expansion of over the Hermite functions reads

(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 ,

(15) | ||||

(16) |

while

(17) |

We can write (15) in tensorial form as

(19) |

where,

(20) |

while, when ,

(21) |

and we have introduced

(22) | ||||

(23) | ||||

(24) |

Equations (20–24) 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)

(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

(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

(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

(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

(29a) | |||

They satisfy the orthogonality condition | |||

(29b) |

and therefore

(30) |

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

(31a) | ||||

(31b) |

Using the identity

(33) |

we can show

(34) |

## Appendix B Representation of translation operator

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

(35) |

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

(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 ,

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.

The truncated Fourier-Hermite expansions read

(37a) | |||

(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 ,

(38) |

while

(39) |

and we have introduced