Lecture #2

Mathematical Methods

# Summary of Lecture #1

Plasmas are composed of ensembles of charged (and neutral particles)
which may, or may not, be in charge and thermal equilibrium. →
Maxwellian distribution.

As a result of the free motion of charged particles, Debye shielding
will occur → Debye length.

The fields produced by density perturbations cause a restoring force,
which leads to plasma oscillations →plasma frequency.

Overview of Lecture #2

Mathematical fundamentals

Distribution function

Vlasov/Boltzmann equation

Methods of solution

Moments

Recall the definition of a distribution function from lecture #1:

> f(x,v) = probability density that particles exist between
> [x]{.underline}, [v]{.underline}, and [x]{.underline} +
> d[x]{.underline}, [v]{.underline} + d[v]{.underline}

This concept is actually based on a more formal expansion of the
probabilities of finding clusters of particles, i.e. correlations
between particles, which we shall mostly ignore in these notes.

The most probable distribution is the Maxwell / Boltzmann, given by
(lecture #1)

$f(v) = \frac{\eta_{o}}{(2\pi\frac{\text{kT}}{m})^{\frac{3}{2}}}e^{- \frac{m\left( v_{x}^{2} + v_{y}^{2} + v_{z}^{2} \right)}{2\text{kT}}}$

though frequently the actual plasma distribution is not in equilibrium
and may, for instance, exhibit multiple temperatures.

$f(v_{\text{//}},v_{}) = \frac{n_{0}}{(\frac{2\pi k}{m})^{\frac{3}{2}}T_{}T_{\text{//}^{\frac{1}{2}}}}e^{- \frac{\text{mv}_{}^{2}}{2\text{kT}_{}} - \frac{\text{mv}_{\text{//}}^{2}}{2{kΤ}_{\text{//}}}}$

Free Energy of a Plasma

It is conceptually useful to consider the concept of [free
energy]{.underline}, which is the energy of a system available to effect
changes in the system. For instance, the Maxwell-Boltzmann distribution
is the one that occurs at a minimum of the free energy. The free energy
of a plasma is defined by

F = W -- TS

energy entropy

temperature

# F is defined by 

> $$F = - \frac{\text{kT}}{V}\text{ln}\int_{}^{}{\text{exp}( - \frac{1}{\text{kT}}}\sum_{i,k}^{}W_{\text{ik}})\text{dx}_{1}\text{.}\text{.}\text{.}\text{dx}_{N}$$

W~ik~ is the interaction energy, which for the Coulomb interaction is

$W_{\text{ik}} = \frac{q_{i}q_{k}}{\left| {\underline{x}}_{i} - {\underline{x}}_{k} \right|}$

Note: for an interaction free ensemble F=F~o~

Important note: Neglecting correlations is not the same as neglecting
interactions entirely. The potential at a point x is given by
[the]{.underline} [collective]{.underline} potential due to all other
particles in the system.

$\varphi_{p} = \sum_{j}^{}\frac{q_{j}}{\left| {\underline{x}}_{p} - {\underline{x}}_{j} \right|}$

this can be approximated by (using the smoothness criterion)

$\varphi_{p} = \int_{}^{}\frac{d{\underline{x}}_{j}\text{qn}\left( {\underline{x}}_{j} \right)}{\left| {\underline{x}}_{p} - {\underline{x}}_{j} \right|}$

In order to handle the plasma medium as a continuum, it is necessary to
have many particles within a Debye sphere.

A brief (and superficial) summary of Plasma Statistical Mechanics.

We have already introduced the concept of a [distribution
function.]{.underline}

Formally, this is the probability that a system occupies a given state
in 6-D (3-p/3-q) phase space.

Actually, it is the [joint]{.underline} probability that the system
falls within an infinitesimal volume in phase space, namely that
particle 1 is the volume dx~1~ dv~1~, [and]{.underline} particle 2 is in
dx~2~ dv~2~, and so on. We will ignore higher-order terms and just
consider the [reduced probability]{.underline} that is given by

$f_{k}\left( x_{},v_{},\text{.}\text{.}\text{.},x_{k},v_{k} \right) = \int_{}^{}{\text{dx}_{k + 1}\text{dv}_{k + 1},\text{.}\text{.}\text{.}\text{dx}_{n}}\text{dv}_{n}f_{\text{jo}\text{int}}$

This is what we have been calling our distribution function.

We will actually only care about average values, such as

$\overline{Q} = \frac{1}{T}\int_{- T2}^{T2}{Q(q,p,t)\text{dt}}$ time
average

$\left\langle Q \right\rangle = \frac{1}{v}\int_{}^{}{Q(q,p,t)\text{fdq}_{}}\text{.}\text{.}\text{.}\text{dp}_{N}$
ensemble average

Now the number of particles in phase space must be conserved (apart from
specific sources and sinks).

Thus in a differential volume, we consider first motion along the q~1~
axis.

> ![](media/image1.jpeg){width="3.051388888888889in"
> height="1.9104166666666667in"}

The number of systems passing q~0~ per unit time is

$f\left( q_{0},p,t \right){\overset{˙}{q}}_{1}{\Delta q}_{2}\text{.}\text{.}\text{.}{\Delta q}_{N}{\Delta p}_{2}\text{.}\text{.}\text{.}{\Delta p}_{N}$

and at the $q_{0} + {\Delta q}_{1}$ surface, the number leaving the
volume per unit time is

> $$\left( f + \frac{\partial f}{\partial q_{1}} \right)\left( \overset{\text{.}}{{\overset{˙}{q}}_{1} + \frac{\partial{\overset{˙}{q}}_{1}}{\partial q_{1}}{\Delta q}_{1}} \right){\Delta q}_{2}\text{.}\text{.}\text{.}{\Delta q}_{N}{\Delta p}_{1}\text{.}\text{.}\text{.}{\Delta p}_{n}$$

which means, to lowest order in ${\Delta q}_{1}$, the change in the
number of states in $\Delta V$ is

$- \left( f\frac{\partial{\overset{˙}{q}}_{1}}{\partial q_{1}} + {\overset{˙}{q}}_{1}\frac{\partial f}{\partial q_{1}} \right){\Delta q}_{1}\text{.}\text{.}\text{.}{\Delta q}_{N}{\Delta p}_{1}\text{.}\text{.}\text{.}{\Delta p}_{N}$

But the same result is true for every p~i~ and q~i~ coordinate, thus

$\frac{\partial f}{\partial t} = - \sum_{i}^{}\left( f\left\lbrack \frac{\partial{\overset{˙}{q}}_{i}}{\partial q_{i}} + \frac{\partial{\overset{˙}{p}}_{i}}{\partial p_{i}} \right\rbrack + \left\lbrack \frac{\partial f}{\partial q_{i}}{\overset{˙}{q}}_{i} + \frac{\partial f}{\partial p_{i}}{\overset{˙}{p}}_{i} \right\rbrack \right)$

But recall that for any system, Hamilton's equations are given by

${\overset{˙}{q}}_{i} = \frac{\partial H}{\partial p_{i}}$
;${\overset{˙}{p}}_{i} = - \frac{\partial H}{\partial q_{i}}$

where H is the [Hamiltonian]{.underline}.

Thus, we have

$\frac{\partial{\overset{˙}{q}}_{i}}{\partial q_{i}} = - \frac{\partial{\overset{˙}{p}}_{i}}{\partial p_{i}}$

and the first term on the RHS of the $\frac{\partial f}{\partial t}$
equation = 0.

$\frac{\partial f}{\partial t} = - \sum_{i}^{}\left( \frac{\partial f}{\partial q_{i}}\frac{\partial H}{\partial p_{i}} - \frac{\partial f}{\partial p_{i}}\frac{\partial H}{\partial q_{i}} \right)$

But this equation implies that

$\frac{\partial f}{\partial t} + \sum_{i}^{}\left( \frac{\partial f}{\partial q_{i}}\frac{\partial H}{\partial p_{i}} - \frac{\partial f}{\partial p_{i}}\frac{\partial H}{\partial q_{i}} \right) \equiv \frac{\text{df}}{\text{dt}} = 0$

This is know as Liouville's Theorem and may be stated

-   Systems of an ensemble are conserved [along a trajectory of the
    system.]{.underline}

What does this mean?

-   individual particles, and hence the system, must obey Newton's laws

-   as the system evolves, phase space is conserved

This is a very general result and describes the behavior, in all its
complexity, of virtually all many-body systems. But it is also why
'plasma-like phenomena' can be found throughout physics.

We note that all the features peculiar to a plasma are imbedded in the
force term.

The general force on a charged particle is given by the [Lorentz force
law]{.underline}

$F = q\left( \underline{E} + \underline{v} \times \underline{B} \right)$

(we will neglect gravitational forces here), which has the equivalent
(non-relativistic) Hamiltonian

> $$H = \frac{1}{2m}\left( \underline{P} - q\underline{A} \right)^{2} - q\varphi$$

where p is the momentum, A is the vector potential, and φ is the scalar
potential. (note that q here means charge, not position coordinate)

So finally, we may write the general (non-relativistic) form of the
[Vlasov Equation]{.underline} as

$\frac{\partial f_{\alpha}}{\partial t} + \underline{v} \cdot \frac{\partial f}{\partial\underline{x}} + \frac{q_{\alpha}}{m_{\alpha}}\left( \underline{E} + \underline{v} \times \underline{B} \right) \cdot \frac{\partial f_{\alpha}}{\partial\underline{v}} = 0$

α= species

In general, this is a very difficult equation to solve, as it is
manifestly nonlinear. Moreover, it does not include any effects of
particle interaction other than through the long-range fields E and B.

The simplest corrections to this picture are two-body collisions, and
this can be written schematically as

$\frac{\partial f_{\alpha}}{\partial t}|_{\text{coll}} = \sum_{\beta}^{}{\int_{}^{}\left( f_{\alpha^{\prime}}f_{\beta^{\prime}} - f_{\alpha}f_{\beta} \right)}v_{\text{rel}}\text{bdbd}{\theta d}^{3}v$

The primes correspond to distributions following a collision, the
unprimed quantities to distributions before a collision. B and θ are
geometric parameters of the collision itself, which we will examine more
in later detail.

This slightly more general result is known as the [Boltzmann
Equation.]{.underline}

$\frac{\partial f_{\alpha}}{\partial t} + \underline{v} \cdot \frac{\partial f_{\alpha}}{\partial\underline{x}} + \frac{q_{\alpha}}{m_{\alpha}}\left( \underline{Ε} + \underline{v} \times \underline{B} \right) \cdot \frac{\partial f_{\alpha}}{\partial\underline{v}} = \frac{\partial f_{\alpha}}{\partial t}|_{\text{coll}}$

![](media/image2.png)

These equations contain virtually all the complexities of plasma
physics. To begin to understand their properties, the usual procedure is
to linearize them, and then find solutions in some simple geometric
configuration. They are still hard to solve, and can involve
infinite-order differential equations. (what does that mean?)

We will find some simple linear solutions, and will find nonlinear
solutions for some particular cases. This is where the real interest of
applying plasma physics concepts to beams is: in the detailed
understanding of the Vlasov/Boltzmann equation and in methods of its
solution. There is a wealth of new understanding to be gained.

We note that, in general, the Vlasov equation is a statement of the
constancy of the [total derivative]{.underline} of the distribution
function [along a particle trajectory]{.underline}. Thus if we construct
a distribution function out of [constants of the motion]{.underline}, we
have found stationary, or equilibrium, solutions to the equation.

Example:

Equilibrium solutions when E=B=0

In this case the energy ½ mv^2^ and the momentum mv are constants of the
motion. Thus [any]{.underline} function of velocity is an equilibrium
solution.

> $f_{\text{eq}} = f\left( v_{x},v_{y},v_{z} \right)$

Equilibrium solutions when E=0, B=constant

In this case, the longitudinal momentum mv~z~ and the transverse energy
$\frac{1}{2}m\left( v_{x}^{2} + v_{y}^{2} \right)$ are constants of the
motion. Thus equilibrium solutions are of the form

$f_{\text{eq}} = f\left( v_{},v_{z} \right)$

Equilibrium solutions when B=0 $Ε = Ε(x)\overset{}{x}$

In this case, the longitudinal momentum and the transverse energy given
by $\frac{1}{2}\text{mv}_{x^{2}} + q\varphi$ are constants of the
motion.

$f_{\text{eq}} = f\left( v_{x}^{2} + \frac{2q\varphi}{m},v_{y},v_{z} \right)$

Example: Electrostatic Waves

One of the immediate consequences of the Vlasov equation is the
existence of [electrostatic waves]{.underline} in a plasma. These are
waves which can exist which have no magnetic field and where energy is
exchanged between the particles and the electric field, producing wave
motion.

The general concept of linearization can be demonstrated by the
following:

Let the distribution be given by an equilibrium part, and by a
fluctuating part that is assumed to be small, i.e.

$f = f_{0} + f_{}(x,t)$ $\left| f_{} \right|\text{<<}f_{0}$

Then consider only small sinusoidal perturbations (equivalent to taking
a Fourier transform of the Vlasov equation).

$- {i\omega f}_{}(\omega,k) + \text{ikvf}_{}(\omega,k) = - \frac{q}{m}Ε(\omega,k)\frac{\partial f_{0}}{\partial v}$

In the product term of E and $\frac{\partial f}{\partial v}$, we only
use the equilibrium part of $f$ so that the term is of first order.

This can be solved for $f_{1}$ in terms of E as

$f_{} = - \frac{iqm}{\omega - \text{kv}}Ε\frac{\partial f_{0}}{\partial v}$

As before in these notes, we use Poisson's equation to relate the field
to the perturbation of the particle density which produced it

> $$\text{ik}Ε = \frac{e}{Ε_{o}}\left( n_{i} - n_{e} \right) = \frac{e}{Ε_{o}}\int_{}^{}\left( f_{i} - f_{e} \right)\text{dv}$$

Using this equation, we can eliminate E from the above equations,
leading to the [dispersion relation]{.underline} of electrostatic waves.

$1 + \frac{e^{2}}{m_{e}Ε_{0}k}\int_{- \infty}^{\infty}\left( \frac{m_{e}}{m_{i}}\frac{\frac{\partial f_{0i}}{\partial v}}{\omega - \text{kv}} - \frac{\frac{\partial f_{0e}}{\partial v}}{\omega - \text{kv}} \right)\text{dv} = 0$

This dispersion relation has many interesting properties, like Landau
damping, which we shall examine more later. For now, we note the
similarity to the relation for longitudinal oscillations in a coasting
beam.

$1 - Z(\omega)\frac{i\left( {e\omega}_{0^{2}} \right)}{2\pi}N\int_{\infty}^{\infty}\frac{\frac{\partial f_{0}}{\partial Ε}}{\omega - {m\omega}_{0} + k_{0}Ε} = 0$

These look the same because they come from the same Vlasov equation. The
only difference is the force law that relates the density perturbation
to the induced field. In the case of the beam it is the wakefield, i.e.
the wall currents and their retarded effect on the passing beam. In the
case of the plasma, (for this unmagnetized case), it is the space charge
force of neighboring charges.

Note that the impedance in the beam case takes on the role of a mass,
and can be negative, or complex. Thus a beam can exhibit plasma
oscillations, but the associated mass is complex.

The resonant denominator of the dispersion integral evokes a particular
physical picture. It suggests a resonance occurs when ω = kv namely when
particles are just traveling with the velocity of a wave which has the
phase velocity. This is, in fact the resonant energy exchange which
leads to Landau damping.

Moments of the Vlasov Equation

A connection to the classical theory of fluids can be drawn by
considering [moments]{.underline} of the Vlasov equation. If we assume
the equilibrium function is a known function, say a Maxwell-Boltzmann
distribution, then we may define the following average over the
distribution function:

$n = \int_{}^{}f_{0}d^{3}v$ density

$\left\langle m\underline{v} \right\rangle = \int_{}^{}{m\underline{v}f_{0}}d^{3}v$
momentum

$\left\langle 12m{\underline{v}}^{2} \right\rangle = 12\int_{}^{}{m\underline{v}}^{2}f_{0}d^{3}v$
energy

$\left\langle 12m{\underline{v}}^{2}\underline{v} \right\rangle = 12\int_{}^{}{m\underline{v}}^{2}\underline{v}f_{0}d^{3}v$
energy flux etc.

These are the elements of classical fluid theory and represent an
approximation to the solution of the Vlasov equation. Information about
the internal structure of the distribution function is almost completely
lost by the averaging process.

To derive the moment equations, we multiply the Vlasov equation, given
by

$\frac{\partial f_{\alpha}}{\partial t} + v_{i}\frac{\partial f_{\alpha}}{\partial x_{i}} + a_{i}\frac{\partial f_{\alpha}}{\partial v_{i}} = 0$

by a function $Q\left( \underline{x},\underline{v},t \right)$ and
integrate over velocity. The first term is given by

$\int_{}^{}\frac{\partial f_{\alpha}}{\partial t}\text{Qd}^{3}v = \frac{\partial}{\partial t}n_{\alpha}\left\langle Q \right\rangle_{\alpha} - n_{\alpha}\left\langle \frac{\partial Q}{\partial t} \right\rangle_{\alpha}$

where α is a species tag and

$\left\langle Q \right\rangle_{\alpha} \equiv \frac{1}{n_{\alpha}}\int_{}^{}f_{\alpha}\text{Qd}^{3}v$

In the second term we get

$$\int_{}^{}v_{i}\frac{\partial f_{\alpha}}{\partial x_{i}}\text{Qd}^{3}v = \frac{\partial}{\partial x_{i}}n_{\alpha}\left\langle \text{Qv}_{i} \right\rangle_{\alpha} - n_{\alpha}\left\langle v_{i}\frac{\partial Q}{\partial x_{i}} \right\rangle_{\alpha}$$

And in the third term we get

$\int_{}^{}a:\frac{\partial f_{\alpha}}{\partial v_{i}}\text{Qd}^{3}v = \int_{}^{}\frac{\partial}{\partial v_{i}}\left( a_{i}f_{\alpha}Q \right)d^{3}v - \int_{}^{}a_{i}f_{\alpha}\frac{\partial Q}{\partial v_{i}}d^{3}v$

but

$\int_{}^{}\frac{\partial}{\partial v_{i}}\left( a_{i}f_{\alpha}Q \right) = 0$

$\int_{}^{}a_{i}\frac{\partial f_{\alpha}}{\partial v_{i}}\text{Qd}^{3}v = - n_{\alpha}\left\langle a_{i}\frac{\partial Q}{\partial v_{i}} \right\rangle_{\alpha}$

Thus, the general [transport equation]{.underline} is given by

$\begin{array}{r}
\frac{\partial}{\partial t}n_{\alpha}\left\langle Q \right\rangle_{\alpha} - n_{\alpha}\left\langle \frac{\partial Q}{\partial t} \right\rangle_{\alpha} + \frac{\partial}{\partial x_{i}}n_{\alpha}\left\langle \text{Qv}_{i} \right\rangle_{\alpha} \\
 - n_{\alpha}\left\langle v_{i}\frac{\partial Q}{\partial x_{i}} \right\rangle_{\alpha} - n_{\alpha}\left\langle a_{i}\frac{\partial Q}{\partial v_{i}} \right\rangle_{\alpha} = 0
\end{array}$

And if we consider collisions, (Boltzmann equation), the RHS becomes

$\text{RHS} \equiv {\delta Q}_{\alpha} = \int_{}^{}{Q\frac{\partial f_{\alpha}}{\partial t}}|_{\text{coll}}d^{3}v$

Setting Q=m~α~ give an equation for [mass conservation]{.underline}

$\left\langle Q \right\rangle_{\alpha} = m_{\alpha}$

$\begin{array}{r}
\frac{\partial Q}{\partial t} = \frac{\partial Q}{\partial x_{i}} = \frac{\partial Q}{\partial v_{i}} = 0 \\
 \Rightarrow \left\langle \text{Qv}_{i} \right\rangle_{\alpha} = m_{\alpha}u_{\alpha i}
\end{array}$

where

$\begin{array}{r}
{\underline{u}}_{\alpha} = \frac{1}{n_{\alpha}}\int_{}^{}f_{\alpha}\underline{v}d^{3}v \\
 \Rightarrow \frac{\partial n}{\partial t} + \nabla \cdot \left( n_{\alpha}{\underline{u}}_{\alpha} \right) = 0
\end{array}$

This is also known as the [continuity equation]{.underline}. At the next
higher velocity moment, we find an equation for momentum (force) balance

$\begin{array}{r}
Q = m\left( v_{k} - U_{\alpha}k \right) \\
 \Rightarrow m_{\alpha}n_{\alpha}\left( \frac{\partial}{\partial t} + U_{\alpha i}\frac{\partial}{\partial x_{i}} \right)U_{\alpha k} = m_{\alpha}n_{\alpha}\left\langle \text{ak} \right\rangle_{\alpha} - \frac{\partial p_{\alpha\text{ik}}}{\partial x_{i}}
\end{array}$

convective derivative pressure term

> $\rho_{\alpha\text{ik}} =$ pressure tensor
> $= n_{\alpha}\left\langle v_{j}\left( v_{k} - U_{\alpha k} \right) \right\rangle$

if the pressure is a scalar, then
$p_{\alpha} = n_{\alpha}\text{kT}_{\alpha}$

For the case of electrostatic waves considered previously in the context
of the Vlasov equation, a fluid model for electrostatic waves using
velocity moments can also be found. The corresponding dispersion
relation is given by

$\omega^{2} = \omega_{p^{2}}\left( 1 + 2v_{\text{th}}\frac{k}{\omega} + 3v_{\text{th}}^{2}\frac{k^{2}}{\omega^{2}} + \text{.}\text{.}\text{.} \right)$

This, by the way, can also be derived directly from the previous general
dispersion relation by integrating the dispersion integral by parts, and
truncating the expansion at second-order. The dispersion relation for
electrostatic waves is shown below.

> ![](media/image3.jpeg){width="4.749305555555556in"
> height="4.391666666666667in"}
>
> Example: Two-stream Instability -- fluid model

We consider the simple case of a uniform plasma in which the ions are
fixed, and the electrons have an average velocity $v_{o}$ with respect
to the ions. We assume the plasma is cold
$\left( T_{e} = T_{i} = 0 \right)$. The linearized equations of motion
are for the fluid quantities are

$\begin{array}{r}
m_{i}n_{0}\frac{\partial v_{i}}{\partial t} = \text{en}_{0}Ε \\
m_{e}n_{0}\left( \frac{\partial v_{e}}{\partial t} + v_{o}\frac{\partial v_{e}}{\partial x} \right) = - \text{en}_{0}Ε
\end{array}$

![](media/image4.jpeg){width="4.35in" height="4.018055555555556in"}

We take the Fourier transform of these equations, and look for
electrostatic solutions, namely those where the electric field is purely
longitudinal with respect to the drift

$\begin{array}{r}
 - {i\omega m}_{i}n_{o}v_{i} = \text{en}_{o}Ε \\
m_{e}n_{o}\left( - i\omega + \text{ikv}_{o} \right)v_{e} = - \text{en}_{o}Ε
\end{array}$

We can readily solve for the first-order velocities

> $\begin{array}{r}
> v_{i} = \frac{\text{ie}}{m_{i}\omega}Ε \\
> v_{e} = - \frac{\text{ie}}{m_{e}}\frac{Ε}{\omega - \text{kv}_{0}}
> \end{array}$

The continuity equation for ions is given by

$n_{i} = \frac{k}{\omega}n_{0}v_{i} = \frac{\text{ien}_{0}k}{m_{i}\omega^{2}}Ε$

And the continuity equation for electrons is

$\left( - i\omega + \text{ikv}_{0} \right)n_{e} + \text{ikn}_{0}v_{e} = 0$

And we use the good old Poisson equation to relate the field to the
perturbed density

$\text{ik}\varepsilon_{0}Ε = e\left( \text{ien}_{0}kΕ \right)\left( \frac{1}{m_{i}\omega^{2}} + \frac{1}{m_{e}\left\lbrack \omega - \text{kv}_{0} \right\rbrack^{2}} \right)$

Eliminating E leads to the following dispersion relation

$1 = \omega_{{\rho e}^{2}}\left( \frac{m_{e}m_{i}}{\omega^{2}} + \frac{1}{\left\lbrack \omega - \text{kv}_{0} \right\rbrack^{2}} \right)$

This is a fourth-order polynomial in $\omega$ and is fairly tedious to
solve exactly. But a rough idea of the behavior of the roots can be
found by examining just the RHS of the above dispersion relation. A
typical result is shown below.

> ![](media/image5.jpeg){width="5.05in" height="3.0930555555555554in"}
>
> For this case, there are four real roots for $\omega$and no wave
> growth. But if a slightly smaller
> $\frac{\text{kv}_{0}}{\omega_{\rho}}$ is chosen, then the following
> picture can result.

![](media/image6.jpeg){width="5.149305555555555in"
height="2.986111111111111in"}

In this case two of the roots must be complex (conjugates), which means
one of the modes is an unstable, growing wave. This is the classical
[two-stream instability]{.underline}.

Physically, the plasma frequency of electrons in the moving frame is
Doppler-shifted by the motion. It can then appear in the lab-frame as
the ion plasma frequency, and then coupling of the free-energy
associated with the drifting motion of the electrons to the ions can
occur.

Another way to interpret the instability is to realize that the waves
generated have effectively [negative energy]{.underline}. If the energy
of the drifting electrons alone is

$Ε_{0} = 12\text{mn}_{0}v_{0^{2}}$

then the energy of the perturbed system is given by

$Ε_{1} = 12m\left( n_{0} + n_{1} \right)\left( v_{0} + v_{1} \right)^{2}$

Due to the phase relation forced between n and v by the continuity
relation, the cross-product actually can reduce the overall energy of
the system relative to the unperturbed state. This is a prerequisite for
instability.

$\begin{array}{r}
\text{ΔΕ} = Ε_{1} - Ε_{0} \simeq v_{0^{2}}\frac{\text{kn}_{0}v_{e}}{\omega - \text{kv}_{0}} \\

\end{array}$

> $\Rightarrow \text{ΔΕ} < 0$ when $\text{kv}_{0} > \omega$

Example:

Two-stream instability in a storage ring: trapped particle modes

In a nonneutral stored beam, particles of opposite charge can become
trapped in the potential well of the beam.

![](media/image7.jpeg){width="3.3020833333333335in" height="2.36875in"}

For instance, in a stored anti-proton beam, ions can accumulate in the
negative potential well. The trapped particles undergo bounce
oscillations which can couple to the natural transverse (betatron) modes
of the beam itself.

The Vlasov equations for the beam and trapped ion populations are

$\begin{array}{r}
\frac{\partial\psi_{b}}{\partial t} + {\overset{˙}{y}}_{b}\frac{\partial\psi_{b}}{\partial y_{b}} + f_{b}\frac{\partial\psi_{b}}{\partial{\overset{˙}{y}}_{b}} + w_{0}\frac{\partial\psi_{b}}{\partial\theta} = 0 \\
\frac{\partial\psi_{i}}{\partial t} + {\overset{˙}{y}}_{i}\frac{\partial\psi i}{\partial y_{i}} + F_{i}\frac{\partial\psi_{i}}{\partial{\overset{˙}{y}}_{i}} = 0
\end{array}$

where the coupling forces are given by

$\begin{array}{r}
F_{b} = - \omega_{b}^{2}\left( y_{b} - \overline{y_{i}} \right) + F_{\omega} \\
F_{i} = - \omega_{i}^{2}\left( y_{i} - \overline{y_{b}} \right)
\end{array}$

> $\omega_{b}$ is the particle's betatron frequency and $\omega_{i}$ is
> the bounce frequency of the trapped ions. $F_{\omega}$ is the wake
> force given in terms of a wall impedance

$F_{\omega} = i\frac{\text{Ne}^{2}\omega_{o}}{{\gamma m}_{b}C}z(\Omega)$

Assuming small amplitude oscillations, we linearize this set by taking

$\begin{array}{r}
\psi_{b} = \psi_{\text{bo}} + \psi_{\text{bi}}e^{i(n\theta - \Omega t)} \\
\psi_{i} = \psi_{\text{io}} + \psi_{\frac{i}{}e^{i(n\theta - \Omega t)}}
\end{array}$

The linearized equations become

$\begin{array}{r}
 - i\left( \Omega - {n\omega}_{o} \right)\psi_{\frac{b}{} - w_{b}\frac{\partial\psi_{b1}}{\partial\varphi_{b}} + \left( \frac{\text{iNe}^{2}\omega_{0}}{{\gamma m}_{b}\text{CZ}(\Omega)}Y_{b} + \omega_{i^{2}}Y_{i} \right)\frac{\partial\psi_{\text{bo}}}{\partial{\overset{˙}{y}}_{b}} = 0} - i\Omega\psi_{\frac{i}{} - \omega_{i}\partial\frac{\psi_{\frac{i}{}}}{\partial\varphi_{i}} + w_{i^{2}}Y_{b}\frac{\partial\psi_{i0}}{\partial{\overset{˙}{y}}_{i}} = 0} \\

\end{array}$

which leads to the following dispersion relation

$1 = i\frac{\text{Ne}^{2}\omega_{o}}{\gamma_{\text{mb}}C}\int_{}^{}\frac{g(p)\text{dp}}{{\omega b}^{2} - \left( \Omega - {n\omega}_{0} \right)^{2}}\left( Z\lbrack\Omega\rbrack + Z_{i}\lbrack\Omega\rbrack \right)$

> $g(p) =$ form factor

Z~i~ is the effective impedance presented by the trapped ions

$Z_{i}(\Omega) = - i\frac{\gamma_{\text{mb}}C{\hat{\omega}}_{b}^{2}}{\omega_{0}\text{Ne}^{2}}\int_{}^{}\frac{\psi^{'}\left( \omega_{i^{2}} + \omega_{b^{2}} \right)\omega_{i}\left| y_{i} \right|^{2}\text{dy}_{i}}{\Omega^{2} - \omega_{i^{2}}}$

The ion impedance shows a maximum when the mode frequency Ω is near the
bounce frequency. This collective mode becomes unstable when the bounce
frequency overlaps with the beam's betatron frequency (or its
harmonics).

Beam Distributions

The distribution appropriate for a relativistic beam may be found in a
similar fashion. In the beam frame, the particle motion is given by the
familiar Hamiltonian

$H = \frac{m}{2}\left( v_{x}^{2} + v_{y}^{2} + v_{z}^{2} \right) + q\varphi(x,y,z)$

so that the time-independent distribution must be of the form

$f(H) = f_{0}\text{exp}\left( - \frac{H}{\text{kT}} \right)$

Using our relativistic scaling, we can write this in the lab frame as

$f_{\text{lab}} = A\text{exp}\left\{ - \frac{\gamma\text{mv}_{}^{2} + \gamma^{3}{m\Delta v}_{\text{//}}^{2}}{2\text{kT}_{b\gamma}} - \frac{q\varphi}{\text{kT}_{b\gamma}} \right\}$

which in the general case can be written

$\begin{array}{r}
f_{\text{lab}} = A\text{exp}\left\{ - \frac{\gamma\text{mv}_{}^{2}}{2\text{kT}_{}^{\text{lab}}} - \frac{q\varphi}{\text{kT}_{{}^{\text{lab}}}} \right\} \\
\text{exp}\left\{ - \frac{\gamma^{3}{m\Delta v}_{\text{//}}^{2}}{2\text{kT}_{\text{//}}^{\text{lab}}} - \frac{q\varphi}{\text{kT}_{\text{//}}^{\text{lab}}} \right\}
\end{array}$

Let us assume the potential is made up of an externally-applied
[linear]{.underline} part

$\varphi_{e} = 12{\gamma m}_{0}v_{0}^{2}k_{0}^{2}r^{2}$

where $k_{0}$ is the betatron wavenumber and $v_{0}$ is the mean
transverse velocity, and $\varphi_{s}$ is the self-potential.

The transverse distribution is then given by

$f_{} = f_{0}\text{exp}\left\{ - \frac{\gamma\text{mv}_{0}^{2}\left( x^{} + y^{} + k_{0}^{2}r^{2} \right)}{2\text{kT}_{}} - \frac{{q\varphi}_{s}}{\text{kT}_{}\gamma^{2}} \right\}$

Integrating with respect to $v_{}$ gives the [Boltzmann
relation]{.underline} for the transverse distribution

$n(r) = n_{o}\text{exp}\left\{ - \frac{\gamma\text{mv}_{0}^{2}k_{o}^{2}r^{r}}{2\text{kT}_{}} - \frac{{q\varphi}_{s}}{\text{kT}_{}\gamma^{2}} \right\}$

Taking the limit as $T_{} \rightarrow 0$ yields a constant density with
the constraint that

$\frac{\gamma\text{mv}_{0}^{2}k_{0}}{2} = \frac{q^{2}n_{0}}{4\varepsilon_{0}\gamma^{2}}$

This is the usual situation in a space-charge-dominated beam.

In the case where $T_{} \rightarrow \infty$, the space charge term is
negligible and the distribution is given by

$n(r) = n_{0}\text{exp}\left( - \frac{\gamma\text{mv}_{0}^{2}r^{2}}{2\text{kT}_{}} \right)$

For intermediate cases, the profile must be found numerically using
Poisson's equation to relate n and φ~s~.

![](media/image8.jpeg){width="5.875in" height="3.0375in"}

Ref: Reiser

Method of Multiple Scales

We now turn to a consideration of a formal solution of the full Vlasov
equation. Because it is a nonlinear equation, a common method is to
solve by means of a perturbation expansion around the linear solution.
But this approach can run into problems.

As an example, consider the simple equation of the [anharmonic
oscillator]{.underline}.

$\frac{\text{dx}^{2}}{\text{dt}^{2}} + \omega_{0}^{2}x = - {\alpha x}^{2} - {\beta x}^{3}$

Assume a solution in the form of an expansion

$x(t) = x^{(o)} + x^{(1)} + x^{(2)} + \text{.}\text{.}\text{.}$

where each $x^{(j)}$ is assumed to be of higher order than the previous
one, namely

$\frac{x^{(l + 1)}}{x^{(l)}}\text{\textasciitilde}a\text{<<}1$

Substituting into the nonlinear differential equation and equating each
order [separately]{.underline} to zero gives

$\begin{array}{r}
{\overset{¨}{x}}^{(1)} + \omega_{0}^{2}x^{(1)} = 0 \\
{\overset{¨}{x}}^{(2)} + \omega_{0}^{2}x^{(2)} = - \alpha\left\lbrack x^{(1)} \right\rbrack^{2} \\
{\overset{¨}{x}}^{(3)} + \omega_{0}^{2}x^{(3)} = - 2{\alpha x}^{(1)}x^{(2)} - \beta\left\lbrack x^{(1)} \right\rbrack^{3} \\
 \vdots 
\end{array}$

To first-order, we have a simple harmonic oscillator. The second-order
[particular]{.underline} solution is

$x^{(2)}(t) = x_{0^{(2)}} + x_{{2\omega_{o}}^{(2)}}\text{cos}2\omega_{o}t$

which shows a shift in the equilibrium position of the oscillator. At
third-order we have

$${\overset{¨}{x}}^{(3)} + \omega_{0}^{2}x^{(3)} = \left( \frac{5}{6}\frac{\alpha^{2}a^{3}}{\omega_{0}^{2}} - \frac{3}{4}{\beta a}^{3} \right)\text{cos}\omega_{0}t + \left( - \frac{\alpha^{2}a^{3}}{6\omega_{0}^{2}} - \frac{1}{4}{\beta a}^{3} \right)\text{cos}3\omega_{0}t$$

The critical new feature is that lower-order harmonics can beat together
to produce a [secular]{.underline} growth at the third harmonic

$x_{\omega_{0}}^{(3)} = \left( \frac{5}{6}\frac{{\alpha a}^{3}}{\omega_{0}^{2}} - \frac{3}{4}{\beta a}^{3} \right)t\text{sin}\omega_{0}t$

The unlimited growth is a spurious result of our perturbation method. To
find the correct behavior valid for all times (to third-order in
amplitude) we must use an expansion procedure that [also permits the
oscillation frequency to vary with amplitude.]{.underline}

$\begin{array}{r}
x(\Omega t) = x^{(1)} + x^{(2)} + \text{.}\text{.}\text{.} \\
\Omega = \omega^{(o)} + \omega^{(1)} + \omega^{(2)} + \text{.}\text{.}\text{.}
\end{array}$

Following the same procedure as before, we find the additional
conditions on the frequency terms at each order

$\begin{array}{r}
\omega^{(o)} = \omega_{o} \\
\omega^{(1)} = 0
\end{array}$

At third-order, we have

$$2\omega_{0}\omega^{2}\text{cos}\Omega t + \omega_{0}^{2}{\overset{¨}{x}}^{(3)} + \omega_{0}^{2}x^{(3)} = f_{\omega_{0}}^{(3)}\text{cos}\Omega t + f_{3\omega_{0}}^{(3)}\text{cos}3\Omega t$$

Equating like coefficients of $\text{cos}n\Omega t$ terms, we get for
the frequency shift

$\omega^{(2)} = \frac{\omega_{0}}{2}\left( - \frac{5}{6}\frac{\alpha^{2}}{\omega_{0}^{4}} + \frac{3}{4}\frac{\beta}{\omega_{0}^{2}} \right)a^{2}$

And for the third-order amplitude term

$x^{(3)} = a\left( \frac{\alpha^{2}a^{2}}{\text{48}\omega_{0}^{4}} + \frac{{\beta a}^{2}}{\text{32}\omega_{0}^{2}} \right)\text{cos}3\Omega t$

The resulting solution

$x(t) = \sum_{j = 1}^{3}x^{(j)}(t)$

is now valid for all time to order $a^{4}$. This is the essence of the
[multiple-time-scale expansion]{.underline} method.

Therefore, one way to solve the Vlasov equation is to divide the motion
into two different time scales: a slow time-scale over which the average
quantities evolve, and a fast time-scale that describes fluctuations
about the mean. This looks similar to the linearization process, but it
still includes features of the non-linearity of the equation.

To illustrate the analysis, we choose the case of a coasting beam in a
storage ring and consider only the longitudinal degree of freedom. The
canonical variables which take the roles of x and v in this case are θ
and ε. The resulting Vlasov Equation is given by

$\frac{\partial f}{\partial t} + \overset{˙}{\theta}\frac{\partial f}{\partial\theta} + U\frac{\partial f}{\partial\varepsilon} = 0$

We note that the force term is represented by a voltage U. This can be
an external voltage, or one generated by the wall currents (wakes).
Recalling the discussion about the longitudinal oscillations, we note
that the impedance of the ring plays the role of a complex mass, and
that the wakefield is related to the impedance by the well-known
definition

$\begin{array}{r}
U = - ZΙ \\
Ι = {e\omega}_{o}\int_{- \infty}^{\infty}{\text{fd}\varepsilon}
\end{array}$

We now Fourier analyze the Vlasov equation, which owing to the
periodicity of the system, becomes a Fourier series

$f = \sum_{m}^{}f_{m}e^{\text{im}\theta},U = \sum_{m}^{}U_{m}e^{\text{im}\theta}$

The resulting Vlasov Equation now takes the form

$\frac{\partial f_{m}}{\partial t} + \text{im}\left( \omega_{o} + k_{o} \in \right)f_{m} + \frac{{e\omega}_{o}}{2\pi}\sum_{n + k = m}^{}u_{n}\frac{\partial f_{k}}{\partial \in} = 0$

We now choose λ as an expansion parameter, defined to be the ratio of
the time scales $\lambda = \tau_{o}\tau_{1}$

$\begin{array}{r}
\frac{\partial}{\partial t} = \frac{\partial}{\partial\tau_{o}} + \lambda\frac{\partial}{\partial\tau_{1}} + \text{.}\text{.}\text{.} \\
f_{m} = {\lambda f}_{m^{(1)}} + \lambda^{2}f_{m^{(2)}} + \text{.}\text{.}\text{.} \\
U_{m} = {\lambda U}_{m^{(1)}} + \lambda^{2}U_{m^{(2)}} + \text{.}\text{.}\text{.}
\end{array}$

The perturbations are considered to be small, but large enough to
maintain next-order corrections to the linear dynamics.

Using the above assumptions, we substitute into the Vlasov equation,
treating each time scale as an independent variable and find the
following equations by equating each order of λ separately to zero.

First-order:

$\frac{\partial f_{m^{(1)}}}{\partial\tau_{o}} + \text{im}\left( \omega_{o} + k_{o} \in \right)f_{m^{(1)}} + \frac{{e\omega}_{o}}{2\pi}U_{m^{(1)}}\frac{\partial f_{o}}{\partial\varepsilon} = 0$

This is the same as the linearized Vlasov equation we have already
considered. It leads to the dispersion relation for longitudinal
oscillations quoted previously.

$1 = \frac{i\left( {e\omega}_{o} \right)^{2}}{(2\pi)^{2}}Ζ_{m}N\int_{- \infty}^{\infty}\frac{\frac{\partial f_{o}}{\partial\varepsilon}d\varepsilon}{\Omega - m\left( \omega_{o} + k_{o}\varepsilon \right)}$

Second-order:

$$\begin{array}{r}
\frac{\partial f_{m^{(2)}}}{\partial\tau_{0}} + \text{im}\left( \omega_{0} + k_{0}\varepsilon \right)f_{m^{(2)}} + \frac{{e\omega}_{0}}{2\pi}U_{m^{(2)}}\frac{\partial f_{0}}{\partial\varepsilon} = \frac{\text{ie}\omega_{0}}{2\pi}\frac{\frac{\partial f_{0}}{\partial\varepsilon}}{\omega - m\left( \omega_{0} + k_{0}\varepsilon \right)}e^{- i\left( \omega - {m\omega}_{o} \right)\tau_{o}}\frac{\partial U_{m^{(1)}}}{\partial\tau_{1}} \\
 + \frac{i\left( {e\omega}_{0} \right)^{2}}{(2\pi)^{2}}\sum_{m = n + k}^{}{U_{m^{(1)}}U_{k^{(1)}}}e^{- i\left( \omega_{n} + \omega_{k} \right)\tau_{o}}\frac{\partial}{\in}\left( \frac{\frac{\partial f_{0}}{\partial\varepsilon}}{\omega_{k} - k\left\lbrack \omega_{0} + k_{0}\varepsilon \right\rbrack} \right)
\end{array}$$

The second-order equation brings the first true nonlinear corrections to
the dynamics as it involves products of first-order quantities.

For the purpose of this discussion, we have assumed that the equilibrium
distribution function $f_{0}$ is [unchanged]{.underline} by the
high-frequency oscillations, and have dropped
$\frac{\partial f_{0}}{\partial\tau_{0}}$,$\frac{\partial f_{0}}{\partial\tau_{i}}$
from the second-order equation.

The left-hand side represents the linear response at a given frequency,
similar to the first-order result. However, the product term on the RHS
introduces frequency mixing components and allows for the coupling of
multiple frequencies, a characteristically nonlinear process.

We note that when $\omega_{m} = \omega_{n} + \omega_{k}$ there is
resonant driving of the equation which produces so-called [secular
growth]{.underline}. The multiple-time-scale method allows a degree of
freedom to remove this secularity. (why does it matter?)

Setting the secular terms = 0 leads to the following equation for the
slow change of the second-order amplitude

$\frac{\partial U_{m}}{\partial\tau_{1}} = \frac{{e\omega}_{o}}{2\pi}n\sum_{n}^{}U_{n}U_{m - n}\Im(\Omega)$

nonlinear susceptibility

The above formalism is the basis for the study of the wave-wave
interaction in plasmas, and is equally applicable to a variety of
possible couplings in both unbunched and bunched beams.

We will examine this behavior in the future lectures in more detail.

Example: Two-=stream Instability -- kinetic treatment

The mathematical formalism follows directly from that associated with
longitudinal oscillations. The relevant distribution function in this
case has two parts, shown below

![](media/image9.jpeg){width="3.65in" height="2.457638888888889in"}

$\begin{array}{r}
f_{e} = \frac{n_{o}}{\left( \frac{2\pi\text{kT}_{e}}{m_{e}} \right)^{32}}\text{exp}\left\{ \frac{- m_{e}V_{}^{2}}{2\text{kT}_{e}} - \frac{m_{e}\left( V_{z} - V_{0} \right)^{2}}{2\text{kT}_{e}} \right\} \\
f_{i} = \frac{n_{o}}{\left( \frac{2\pi\text{kT}_{i}}{m_{i}} \right)^{32}}\text{exp}\left\{ \frac{- m_{i}V_{}^{2}}{2\text{kT}_{i}} - \frac{m_{i}V_{z}^{2}}{2\text{kT}_{i}} \right\}
\end{array}$

Thus the dispersion relation is given by

$0 = 1 + \frac{\omega_{\text{pe}}^{2}}{k^{2}}\int_{}^{}{\text{dv}\frac{k\frac{\partial f_{e}}{\partial v}}{\omega - \text{kv}}} + \frac{\omega_{\text{pi}}^{2}}{k^{2}}\int_{}^{}{\text{dv}\frac{k\frac{\partial f_{i}}{\partial v}}{\omega - \text{kv}}}$

If we restrict our consideration to the case where the phase velocities
are well removed from the thermal velocities, then integrating the
integrals by parts (and throwing away higher-order terms) results in

$1 - \frac{\omega_{\text{pe}}^{2}}{\left( \omega - \text{kv}_{o} \right)^{2}} - \frac{\omega_{\text{pi}}^{2}}{\omega^{2}} = 0$

which is the same dispersion relation derived from the fluid theory.

There are many variants on this theme, but the basic concept is the
same:

(1) solve for the perturbed density using the continuity equation
    (zeroth moment of the Vlasov Equation)

(2) solve for the perturbed fields in terms of the perturbed density
    using the force balance equation (first moment of the Vlasov
    Equation)

(3) relate the perturbed density to the perturbed fields with one (or
    all) of the Maxwell's equations

Summary of Lecture #2

We have derived the [Vlasov]{.underline} equation which governs the
behavior of all systems composed of distributions of particles. (the
version including two-body interactions between the particles is known
as the [Boltzmann]{.underline} equation.)

We have derived a set of fluid equations, which are the velocity moments
of f, that describe the behavior of the plasma. These are appropriate
when the internal structure of the distribution is unimportant.

We have found dispersion relations for longitudinal oscillations using
the zeroth and first-order moments of the Vlasov Equation, in
conjunction with a field equation (e.g. Poisson's equation).

We have used the [method of multiple scales]{.underline} to construct
both the solution of the linearized Vlasov equation, and the next-order
nonlinear corrections. These will be used in a number of applications
describing important nonlinear effects in plasmas (and beams).

References:

1). N.A. Krall and A.W. Trivilpiece, Principles of Plasma Physics,
McGraw-Hill, 1973

2). J.L. Shohet, The Plasma State, Academic press, 1971

3). R.C. Davidson, Methods in Nonlinear Plasma Theory, Academic Press,
1972