Landau Damping

Let us begin our study of the Vlasov equation by examining what appears, at first sight, to be a fairly simple and straightforward problem: namely, the propagation of small amplitude plasma waves through a uniform plasma possessing no equilibrium magnetic field. For the sake of simplicity, we shall only consider electron motion, assuming that the ions form an immobile, neutralizing background. The ions are also assumed to be singly charged. We shall search for electrostatic plasma waves of the type discussed in Section 5.7. Such waves are longitudinal in nature (i.e., ${\bf E}\propto {\bf k}$), and possess a perturbed electric field, but no perturbed magnetic field.

Our starting point is the Vlasov equation for an unmagnetized, collisionless plasma:

$$\frac{\partial f_e}{\partial t} + {\bf v}\cdot\nabla f_e - \frac{e}{m_e}\, {\bf E}\cdot\nabla_v f_e = 0,$$ (7.1)

where $f_e({\bf r},{\bf v}, t)$ is the ensemble-averaged electron distribution function. The electric field satisfies

$${\bf E} = -\nabla{\mit\Phi}.$$ (7.2)

where

$$\nabla^2{\mit\Phi} = -\frac{e}{\epsilon_0}\left(n-\int \!f_e\,d^3{\bf v}\right).$$ (7.3)

Here, $n$ is the number density of ions (which is the same as the equilibrium number density of electrons).

Because we are dealing with small amplitude waves, it is appropriate to linearize the Vlasov equation. Suppose that the electron distribution function is written

$$f_e({\bf r}, {\bf v}, t) = f_0({\bf v}) + f_1({\bf r}, {\bf v}, t).$$ (7.4)

Here, $f_0$ represents the equilibrium electron distribution, whereas $f_1$ represents the small perturbation due to the wave. Of course, $\int f_0\,d^3{\bf v} =n$, otherwise the equilibrium state would not be quasi-neutral. The electric field is assumed to be zero in the unperturbed state, so that ${\bf E}({\bf r}, t)$ can be regarded as a small quantity. Thus, linearization of Equations (7.1) and (7.3) yields

$$\frac{\partial f_1}{\partial t} + {\bf v}\cdot\nabla f_1 - \frac{e}{m_e}\,{\bf E}\cdot\nabla_v f_0 =0,$$ (7.5)

and

$$\nabla^2{\mit\Phi} = \frac{e}{\epsilon_0}\int\! f_1\,d^3{\bf v},$$ (7.6)

respectively.

Let us now follow the standard procedure for analyzing small amplitude waves, by assuming that all perturbed quantities vary with ${\bf r}$ and $t$ like $\exp[\,{\rm i}\,({\bf k}\cdot{\bf r}-\omega\,t)]$. Equations (7.5) and (7.6) reduce to

$$-{\rm i}\,(\omega - {\bf k}\cdot{\bf v})\, f_1 + {\rm i}\,\frac{e}{m_e}\,{\mit\Phi} \,{\bf k}\cdot\nabla_v f_0 = 0,$$ (7.7)

and

$$-k^2\,{\mit\Phi} = \frac{e}{\epsilon_0} \int\! f_1\,d^3{\bf v},$$ (7.8)

respectively. Solving the first of these equations for $f_1$, and substituting into the integral in the second, we conclude that if ${\mit\Phi}$ is non-zero then we must have

$$1 + \frac{e^2}{\epsilon_0\,m_e\,k^2} \int\frac{ {\bf k}\cdot\nabla_v f_0} {\omega - {\bf k}\cdot{\bf v}}\,d^3{\bf v} = 0.$$ (7.9)

We can interpret Equation (7.9) as the dispersion relation for electrostatic plasma waves, relating the wavevector, ${\bf k}$, to the frequency, $\omega$. However, in doing so, we run up against a serious problem, because the integral has a singularity in velocity space, where $\omega={\bf k}\cdot{\bf v}$, and is, therefore, not properly defined.

The way to resolve this problem was first explained by Landau in a very influential paper that was the foundation of much subsequent work on plasma oscillations and instabilities (Landau 1946). Landau showed that, instead of simply assuming that $f_1$ varies in time as $\exp(-{\rm i}\,\omega\,t)$, the problem must be regarded as an “initial value problem” in which $f_1$ is specified at $t=0$, and calculated at later times. We may still Fourier analyze with respect to ${\bf r}$, so we write

$$f_1({\bf r}, {\bf v}, t) = f_1({\bf v},t)\,\exp\left(\,{\rm i}\,{\bf k}\cdot {\bf r}\right).$$ (7.10)

It is helpful to define $u$ as the velocity component along ${\bf k}$ (i.e., $u= {\bf k}\cdot{\bf v}/k$), and to also define $F_0(u)$ and $F_1(u,t)$ as the integrals of $f_0({\bf v})$ and $f_1({\bf v},t)$, respectively, over the velocity components perpendicular to ${\bf k}$. Thus, Equations (7.5) and (7.6) yield

$$\frac{\partial F_1}{\partial t} + {\rm i}\,k\,u\,F_1 - \frac{e}{m_e}\,E\, \frac{\partial F_0}{\partial u} = 0,$$ (7.11)

and

$${\rm i}\,k\,E =-\frac{e}{\epsilon_0} \int_{-\infty}^{\infty} F_1(u)\,du,$$ (7.12)

respectively, where ${\bf E}= E\,{\bf k}/k$.

In order to solve Equations (7.11) and (7.12) as an initial value problem, we introduce the Laplace transform of $F_1$ with respect to $t$ (Riley 1974):

$$\bar{F}_1(u,p) = \int_0^\infty F_1(u,t)\,{\rm e}^{-p\,t}\,dt.$$ (7.13)

If the rate of increase of $F_1$ with increasing $t$ is no faster than exponential then the integral on the right-hand side of the previous equation converges, and defines $\bar{F}_1$ as an analytic function of $p$, provided that the real part of $p$ is sufficiently large.

Noting that the Laplace transform of $\partial F_1/\partial t$ is $p\,\bar{F}_1-F_1(u,t=0)$ (as is easily shown by integration by parts), we can Laplace transform Equations (7.11) and (7.12) to obtain

$$p\,\bar{F}_1 + {\rm i}\,k\,u\,\bar{F}_1 = \frac{e}{m_e}\,\bar{E}\,\frac{\partial F_0}{\partial u} + F_1(u,t=0),$$ (7.14)

and

$${\rm i}\,k\,\bar{E} =-\frac{e}{\epsilon_0}\int_{-\infty}^{\infty}\!\bar{F}_1(u)\, du,$$ (7.15)

respectively. The previous two equations can be combined to give

$${\rm i}\,k\,\bar{E} = -\frac{e}{\epsilon_0}\int_{-\infty}^{\infty... ...\partial u} {p + {\rm i}\,k\,u} + \frac{F_1(u,t=0)}{p+{\rm i}\,k\,u}\right] du,$$ (7.16)

yielding

$$\bar{E}(p) = -\frac{(e/\epsilon_0)}{{\rm i}\,k\,\epsilon(k,p)} \int_{-\infty}^\infty \frac{F_1(u,t=0)}{p+{\rm i}\,k\,u}\,du,$$ (7.17)

where

$$\epsilon(k,p) = 1 + \frac{e^2}{\epsilon_0\,m_e\,k} \int_{-\infty}^{\infty} \frac{\partial F_0/\partial u}{{\rm i}\,p- k\,u}\,du.$$ (7.18)

The function $\epsilon(k,p)$ is known as the plasma dielectric function. Of course, if $p$ is replaced by $-{\rm i}\,\omega$ then the dielectric function becomes equivalent to the left-hand side of Equation (7.9). However, because $p$ possesses a positive real part, the integral on the right-hand side of the previous equation is well defined.

Figure 7.1: The Bromwich contour.

According to Equations (7.14) and (7.17), the Laplace transform of the distribution function is written

$$\bar{F}_1 = \frac{e}{m_e}\,\bar{E}\, \frac{\partial F_0/\partial u} {p+ {\rm i}\,k\,u} + \frac{F_1(u,t=0)}{p+ {\rm i}\,k\,u},$$ (7.19)

or

$$\bar{F}_1(u,p) = - \frac{e^2}{\epsilon_0\,m_e\,{\rm i}\,k} \frac{... ...ac{F_1(u',t=0)} {p+ {\rm i}\,k\,u'}\,du' + \frac{F_1(u,t=0)}{p+ {\rm i}\,k\,u}.$$ (7.20)

Having found the Laplace transforms of the electric field and the perturbed distribution function, we must now invert them to obtain $E$ and $F_1$ as functions of time. The inverse Laplace transform of the distribution function is given by (Riley 1974)

$$F_1(u,t) = \frac{1}{2\pi\,{\rm i}}\int_C \bar{F}_1(u,p)\,{\rm e}^{p\,t}\,dp,$$ (7.21)

where $C$—the so-called Bromwich contour—is a contour running parallel to the imaginary axis, and lying to the right of all singularities (otherwise known as poles) of $\bar{F}_1$ in the complex-$p$ plane. (See Figure 7.1.) There is an analogous expression for the parallel electric field, $E(t)$.

Rather than trying to obtain a general expression for $F_1(u,t)$, from Equations (7.20) and (7.21), we shall concentrate on the behavior of the perturbed distribution function at large times. Looking at Figure 7.1, we note that if $\bar{F}_1(u,p)$ has only a finite number of simple poles in the region ${\rm Re}(p)>-\sigma$ (where $\sigma$ is real and positive) then we may deform the contour as shown in Figure 7.2, with a loop around each of the singularities. A pole at $p_0$ gives a contribution that varies in time as ${\rm e}^{p_0\,t}$, whereas the vertical part of the contour gives a contribution that varies as ${\rm e}^{-\sigma\,t}$. For sufficiently large times, the latter contribution is negligible, and the behavior is dominated by contributions from the poles furthest to the right.

Figure 7.2: The distorted Bromwich contour.

Equations (7.17), (7.18), and (7.20) all involve integrals of the form

$$\int_{-\infty}^\infty \frac{G(u)}{u-{\rm i}\,p/k}\,du.$$ (7.22)

Such integrals become singular as $p$ approaches the imaginary axis. In order to distort the contour $C$, in the manner shown in Figure 7.2, we need to continue these integrals smoothly across the imaginary $p$-axis. As a consequence of the way in which the Laplace transform was originally defined—that is, for ${\rm Re}(p)$ sufficiently large—the appropriate way to do this is to take the values of these integrals when $p$ lies in the right-hand half-plane, and to then find the analytic continuation into the left-hand half-plane (Flanigan 2010).

If $G(u)$ is sufficiently well-behaved that it can be continued off the real axis as an analytic function of a complex variable $u$ then the continuation of (7.22) as the singularity crosses the real axis in the complex $u$-plane, from the upper to the lower half-plane, is obtained by letting the singularity take the contour with it, as shown in Figure 7.3 (Cairns 1985).

Figure 7.3: The Bromwich contour for Landau damping.

Note that the ability to deform the Bromwich contour into that of Figure 7.3, and so to find a dominant contribution to $E(t)$ and $F_1(u,t)$ from a few poles, depends on $F_0(u)$ and $F_1(u,t=0)$ having smooth enough velocity dependences that the integrals appearing in Equations (7.17), (7.18), and (7.20) can be analytically continued sufficiently far into the lower half of the complex $u$-plane (Cairns 1985).

If we consider the electric field given by the inversion of Equation (7.17) then we see that its behavior at large times is dominated by the zero of $\epsilon(k,p)$ that lies furthest to the right in the complex $p$-plane. According to Equations (7.20) and (7.21), $F_1(u,t)$ has a similar contribution, as well as a contribution that varies in time as ${\rm e}^{-{\rm i}\,k\,u\,t}$. Thus, for sufficiently long times after the initial excitation of the wave, the electric field depends only on the positions of the roots of $\epsilon(k,p)=0$ in the complex $p$-plane. The distribution function, on the other hand, has corresponding components from these roots, as well as a component that varies in time as ${\rm e}^{-{\rm i}\,k\,u\,t}$. At large times, the latter component of the distribution function is a rapidly oscillating function of velocity, and its contribution to the charge density, obtained by integrating over $u$, is negligible.

As we have already noted, the function $\epsilon(k,p)$ is equivalent to the left-hand side of Equation (7.9), provided that $p$ is replaced by $-{\rm i}\,\omega$. Thus, the dispersion relation, (7.9), obtained via Fourier transformation of the Vlasov equation, gives the correct behavior at large times, as long as the singular integral is treated correctly. Adapting the procedure that we discovered using the complex variable $p$, we see that the integral is defined as it is written for ${\rm Im}(\omega)>0$, and analytically continued, by deforming the contour of integration in the $u$-plane (as shown in Figure 7.3), into the region ${\rm Im}(\omega)<0$. The simplest way to remember how to do the analytic continuation is to observe that the integral is continued from the part of the $\omega$-plane corresponding to growing perturbations to that corresponding to damped perturbations. Once we know this rule, we can obtain kinetic dispersion relations in a fairly direct manner, via Fourier transformation of the Vlasov equation, and there is no need to attempt the more complicated Laplace transform solution.

In Chapter 5, where we investigated the cold-plasma dispersion relation, we found that for any given $k$ there were a finite number of values of $\omega$, say $\omega_1$, $\omega_2$, $\cdots$, and a general solution was a linear superposition of functions varying in time as ${\rm e}^{-{\rm i}\,\omega_1\,t}$, ${\rm e}^{-{\rm i}\,\omega_2\,t}$, et cetera. The set of values of $\omega$ corresponding to a given value of $k$ is called the spectrum of the wave. It is clear that the cold-plasma equations yield a discrete wave spectrum. On the other hand, in the kinetic problem, we obtain contributions to the distribution function that vary in time as ${\rm e}^{-{\rm i}\,k\,u\,t}$, with $u$ taking any real value. In other words, the kinetic equation yields a continuous wave spectrum. All of the mathematical difficulties of the kinetic problem arise from the existence of this continuous spectrum (Cairns 1985). At short times, the behavior is very complicated, and depends on the details of the initial perturbation. It is only asymptotically that a mode varying in time as ${\rm e}^{-{\rm i}\,\omega\,t}$ is obtained, with $\omega$ determined by a dispersion relation that is solely a function of the unperturbed state. As we have seen, the emergence of such a mode depends on the initial velocity disturbance being sufficiently smooth.

Suppose, for the sake of simplicity, that the background plasma state is a Maxwellian distribution. Working in terms of $\omega$, rather than $p$, the kinetic dispersion relation for electrostatic waves takes the form

$$\epsilon(k,\omega) = 1 + \frac{e^2}{\epsilon_0\,m_e\,k} \int_{-\infty}^{\infty} \frac{\partial F_0/\partial u}{\omega- k\,u}\,du=0,$$ (7.23)

where

$$F_0(u) = \frac{n}{(2\pi\,T_e/m_e)^{1/2}}\,\exp\left(-\frac{m_e\,u^2}{2\,T_e}\right).$$ (7.24)

Suppose that, to a first approximation, $\omega$ is real. Letting $\omega$ tend to the real axis from the domain ${\rm Im}(\omega)>0$, we obtain

$$\int_{-\infty}^{\infty} \frac{\partial F_0/\partial u}{\omega- k\... ...ac{{\rm i}\,\pi}{k} \left(\frac{\partial F_0}{\partial u} \right)_{u=\omega/k},$$ (7.25)

where $P$ denotes the Cauchy principal part of the integral (Flanigan 2010). The origin of the two terms on the right-hand side of the previous equation is illustrated in Figure 7.4. The first term—the principal part—is obtained by removing an interval of length $2\,\epsilon$, symmetrical about the pole, $u=\omega/k$, from the range of integration, and then letting $\epsilon\rightarrow 0$. The second term comes from the small semi-circle linking the two halves of the principal part integral. Note that the semi-circle deviates below the real $u$-axis, rather than above, because the integral is calculated by letting the pole approach the axis from the upper half-plane in $u$-space.

Incidentally, because Equation (7.25) holds for any well-behaved distribution function, it follows that

$$\frac{1}{\omega-k\,u} = P\,\frac{1}{\omega -k\,u} -{\rm i}\,\pi\,\delta(\omega-k\,u).$$ (7.26)

This famous expression is known as the Plemelj formula (Plemelj 1908).

Figure 7.4: Integration path about a pole.

Suppose that $k$ is sufficiently small that $\omega\gg k\,u$ over the range of $u$ where $\partial F_0/\partial u$ is non-negligible. It follows that we can expand the denominator of the principal part integral in a Taylor series:

$$\frac{1}{\omega-k\,u} \simeq \frac{1}{\omega}\left(1+ \frac{k\,u}{\omega} + \frac{k^2\,u^2}{\omega^2} + \frac{k^3\, u^3}{\omega^3} + \cdots\right).$$ (7.27)

Integrating the result term by term, and remembering that $\partial F_0/\partial u$ is an odd function, Equation (7.23) reduces to

$$1-\frac{{\mit\Pi}_e^{2}}{\omega^2} - 3\,k^2\,\frac{T_e\,{\mit\Pi}... ...\,\pi}{k^2} \left(\frac{\partial F_0}{\partial u}\right)_{u=\omega/k} \simeq 0,$$ (7.28)

where ${\mit\Pi}_e=(n\,e^2/\epsilon_0\,m_e)^{1/2}$ is the electron plasma frequency. Equating the real part of the previous expression to zero yields

$$\omega^2 \simeq {\mit\Pi}_e^{2}\left[1+ 3\,(k\,\lambda_D)^2\right],$$ (7.29)

where $\lambda_D =(T_e/m_e\,{\mit\Pi}_e^{2})^{1/2}$ is the Debye length, and it is assumed that $k\,\lambda_D\ll 1$. We can regard the imaginary part of $\omega$ as a small perturbation, and write $\omega=\omega_0+\delta\omega$, where $\omega_0$ is the root of Equation (7.29). It follows that

$$2\,\omega_0\,\delta\omega \simeq \omega_0^{2}\, \frac{e^2}{\epsil... ...c{{\rm i}\,\pi}{k^2} \left(\frac{\partial F_0}{\partial u}\right)_{u=\omega/k},$$ (7.30)

and so

$$\delta\omega \simeq \frac{{\rm i}\,\pi}{2} \frac{e^2\,\omega_0}{\epsilon_0\,m_e\,k^2} \left(\frac{\partial F_0}{\partial u}\right)_{u=\omega/k},$$ (7.31)

giving

$$\delta\omega \simeq - {\rm i}\sqrt{\frac{\pi}{8}} \frac{{\mit\Pi}_e} {(k\,\lambda_D)^3} \exp\left[-\frac{1}{2\,(k\,\lambda_D)^2}\right].$$ (7.32)

If we compare the previous results with those for a cold plasma, where the dispersion relation for an electrostatic plasma wave was found to be simply $\omega^2={\mit\Pi}_e^{2}$ (see Section 5.7), we see, first, that $\omega$ now depends on $k$, according to Equation (7.29), so that, in a warm plasma, the electrostatic plasma wave is a propagating mode, with a non-zero group-velocity. Such a mode is known as a Langmuir wave. Second, we now have an imaginary part to $\omega$, given by Equation (7.32), corresponding, because it is negative, to the damping of the wave in time. This damping is generally known as Landau damping. If $k\,\lambda_D\ll 1$ (i.e., if the wavelength is much larger than the Debye length) then the imaginary part of $\omega$ is small compared to the real part, and the wave is only lightly damped. However, as the wavelength becomes comparable to the Debye length, the imaginary part of $\omega$ becomes comparable to the real part, and the damping becomes strong. Admittedly, the approximate solution given previously is not very accurate in the short wavelength case, but it is nevertheless sufficient to indicate the existence of very strong damping.

There are no dissipative effects explicitly included in the collisionless Vlasov equation. Thus, it can easily be verified that if the particle velocities are reversed at any time then the solution up to that point is simply reversed in time. At first sight, this reversible behavior does not seem to be consistent with the fact that an initial perturbation dies out. However, we should note that it is only the electric field that decays in time. The distribution function contains an undamped term varying in time as ${\rm e}^{-{\rm i}\,k\,u\,t}$. Furthermore, the decay of the electric field depends on there being a sufficiently smooth initial perturbation in velocity space. The presence of the ${\rm e}^{-{\rm i}\,k\,u\,t}$ term means that, as time advances, the velocity space dependence of the perturbation becomes more and more convoluted. It follows that if we reverse the velocities after some time then we are not starting with a smooth distribution. Under these circumstances, there is no contradiction in the fact that, under time reversal, the electric field grows initially, until the smooth initial state is recreated, and subsequently decays away (Cairns 1985).

Landau damping was first observed experimentally in the 1960s (Malmberg and Wharton 1964; Malmberg and Wharton 1966; Derfler and Simonen 1966).