%
%
% This paper is a LaTeX file. It does not contain the figures.
% To get the figures send an email request to morrison@hagar.ph.utexas.edu
% and ask for a copy of the IFS report
%
%
%\newcommand{\Real} {\Bbb R}
\newcommand{\Integer} {\Bbb Z}
\newcommand {\Natural} {\Bbb N}
\newcommand {\Rational} {\Bbb Q}
\newcommand{\Complex} {\Bbb C}
\newcommand{\Torus} {\Bbb T}
\newcommand{\eqq}[1]{eq.~(\ref{#1})} %Physica D style
\newcommand{\Eq}[1]{Equation~(\ref{#1})} %Physica D style
\newcommand{\eqs}[1]{eqs.~(\ref{#1})} %Physica D style
\renewcommand{\S}[1]{Sec.~\ref{#1}} %Physica D style
\newcommand{\fig}[1]{fig.~(\ref{#1})} %Physica D style
\newcommand{\Section}[1]{\setcounter{equation}{0}\section{#1}}
\newtheorem{theorem}{Theorem}
\newtheorem{corollary} {Corollary}
\newcommand \half {\frac{1}{2}}
\newcommand \Dw {\Delta \omega}
\documentstyle[12pt]{article}
\hfuzz 5pt
\vfuzz 5pt
\renewcommand{\jot}{12pt}
\def\emode{{$\et_e$~mode}\ }
\def\imode{{$\et_i$~mode}\ }
\def\mfp{{\rm mfp}}
\def\smax{{\mathop{\scriptstyle\max}}}
\def\smin{{\mathop{\scriptstyle\min}}}
\def\pmb#1{\setbox0=\hbox{$#1$}
\kern-.025em\copy0\kern-\wd0
\kern.05em\copy0\kern-\wd0
\kern-.025em\raise.0433em\box0 }
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\def\proof{{\parskip7pt\noindent{\sl Proof}.\ \ }}
\def\move{\hspace*{20pt}}
\def\const{{\rm const.}}
\def\kev{{\rm keV}}
\def\mev{{\rm MeV}}
\def\cm{{\rm cm}\!\!}
\def\ibar{{\mbox{$\;\iota\!\!$}}}
\def\etal{{\it et~al.}\,}
\def\lambar{{\mathchar'26\mskip-9mu\lambda}}
\def\EB{{\bfE\times\bfB}}
\def\mhd{{\rm MHD}}
\def\smallin{\hbox{{\small\rm I}\kern-.2em\hbox{{\small\rm N}}}}
\def\IN{\hbox{{\rm I}\kern-.2em\hbox{{\rm N}}}}
\def\IR{\hbox{{\rm I}\kern-.2em\hbox{{\rm R}}}}
\def\IP{\hbox{{\rm I}\kern-.2em\hbox{{\rm P}}}}
\def\Re{\mathop{\rm Re}\nolimits}
\def\Im{\mathop{\rm Im}\nolimits}
\def\crit{{\rm crit}}
\def\ext{{\rm ext}}
\def\sgn{{\rm sgn}\,}
\def\eff{{\rm eff}}
\def\cosh{{\rm cosh}}
\def\Ampere{{\rm Amp\'ere}}
\def\Alfven{{\rm Alfv\'en}\ }
\def\Alfvenic{{\rm Alfv\'enic}\ }
\def\sech{{\rm sech}}
\def\rec{\rule{.8ex}{.9ex}}
\def\bs{{\scriptstyle\pmb\ast}}
\def\beginab{\begin{abstract}}
\def\endab{\end{abstract}}
\def\beginenum{\begin{enumerate}}
\def\endenum{\end{enumerate}}
\def\begindoc{\begin{document}}
\def\enddoc{\end{document}}
\def\bq{\begin{equation}}
\def\eq{\end{equation}}
\def\bqy{\begin{eqnarray}}
\def\eqy{\end{eqnarray}}
\def\bqyn{\begin{eqnarray*}}
\def\eqyn{\end{eqnarray*}}
\def\bc{\begin{center}}
\def\ec{\end{center}}
\def\bfll{\begin{flushleft}}
\def\efll{\end{flushleft}}
\def\bflr{\begin{flushright}}
\def\eflr{\end{flushright}}
\def\bigskip{\vspace{\bigskipamount}}
\def\medskip{\vspace{\medskipamount}}
\def\smallskip{\vspace{\smallskipamount}}
\def\dfr{\displaystyle\frac}
\def\dsp{\displaystyle}
\def\bskip{\baselineskip}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CALIGRAPHY LETTERS (SCRIPT!) (upper case only)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\def\cala{{\cal A}}
\def\calb{{\cal B}}
\def\calc{{\cal C}}
\def\cald{{\cal D}}
\def\cale{{\cal E}}
\def\calf{{\cal F}}
\def\calg{{\cal G}}
\def\calh{{\cal H}}
\def\cali{{\cal I}}
\def\calj{{\cal J}}
\def\calk{{\cal K}}
\def\call{{\cal L}}
\def\calm{{\cal M}}
\def\caln{{\cal N}}
\def\calo{{\cal O}}
\def\calp{{\cal P}}
\def\calq{{\cal Q}}
\def\calr{{\cal R}}
\def\cals{{\cal S}}
\def\calt{{\cal T}}
\def\calu{{\cal U}}
\def\calv{{\cal V}}
\def\calw{{\cal W}}
\def\calx{{\cal X}}
\def\caly{{\cal Y}}
\def\calz{{\cal Z}}
\arraycolsep 1.5pt
\tabcolsep 1.5pt
\def\m@th{\mathsurround=0pt}
\def\n@space{\nulldelimiterspace=0pt \m@th}%1
\def\biggg#1{{\mbox{$\left#1\vbox to 20.5pt{}\right.\n@space$}}}%2
\def\Biggg#1{{\mbox{$\left#1\vbox to 23.5pt{}\right.\n@space$}}}%3
\def\Bigggg#1{{\mbox{$\left#1\vbox to 40pt{}\right.\n@space$}}}%4
\def\Biggggg#1{{\mbox{$\left#1\vbox to 50pt{}\right.\n@space$}}}%5
\def\Bigggggg#1{{\mbox{$\left#1\vbox to 60pt{}\right.\n@space$}}}%6
\def\Biggggggg#1{{\mbox{$\left#1\vbox to 80pt{}\right.\n@space$}}}%7
\def\noal#1{\noalign{\smallskip\noindent\rm#1\smallskip}}
\def\simsim#1{\begin{array}{c}\si\\[-7pt]\approx\end{array}}
\def\wh#1{{\widehat{#1}}}
\def\wt#1{{\widetilde{#1}}}
\def\pp#1{\frac{\p}{\p #1}}
\def\dd#1{\frac{d}{d #1}}
\def\simle{{\mathop{\stackrel{\sim}{\scriptstyle <}}\nolimits}}
\def\simgr{{\mathop{\stackrel{\sim}{\scriptstyle >}}\nolimits}}
\def\gele{{\mathop{\stackrel{>}{\scriptstyle <}}\nolimits}}
\def\lege{{\mathop{\stackrel{<}{\scriptstyle >}}\nolimits}}
\def\dotedot{\mathop{=}\limits^\cdot_{^{\scriptstyle\cdot}}}
%\def\ded{\mathop{\stackrel{\doteq}{\ .\ }}\nolimits}
\def\chix{{{\raise.5ex\hbox{$\chi$}}}}
\def\nms{\mathsurround=0pt}
\def\gtapprox{\mathrel{\mathpalette\overapprox>}} % greater than or approx.
\def\ltapprox{\mathrel{\mathpalette\overapprox<}} % less than or approx.
\def\overapprox#1#2{\lower 2pt\vbox{\baselineskip 0pt \lineskip - 1pt
\ialign{$\nms#1\hfil##\hfil$\crcr#2\crcr\approx\crcr}}}
\def\gtsim{\mathrel{\mathpalette\oversim>}} % greater than or sim.
\def\ltsim{\mathrel{\mathpalette\oversim<}} % less than or sim.
\def\oversim#1#2{\lower 2pt\vbox{\baselineskip 0pt \lineskip 1pt
\ialign{$\nms#1\hfil##\hfil$\crcr#2\crcr\sim\crcr}}}
\def\gtequal{\mathrel{\mathpalette\overequal>}} % greater than or equal.
\def\ltequal{\mathrel{\mathpalette\overequal<}} % less than or equal.
\def\overequal#1#2{\lower 2pt\vbox{\baselineskip 0pt \lineskip 1pt
\ialign{$\nms#1\hfil##\hfil$\crcr#2\crcr=\crcr}}}
\def\leftrightarrowfill{$\nms \mathord\leftarrow \mkern-6mu
\cleaders\hbox{$\mkern-2mu \mathord- \mkern-2mu$}\hfill
\mkern-6mu \mathord\rightarrow$}
\def\overdoublearrow#1{\vbox{\ialign{##\crcr
\leftrightarrowfill \crcr\noalign{\kern-1pt\nointerlineskip}
$\hfil\displaystyle{#1}\hfil$\crcr}}}
\def\mapright#1{\smash{\mathop{\longrightarrow}\limits^{#1}}}
\def\maprightsuper#1{\smash{\mathop{\longrightarrow}\limits^{#1}}}
\def\mapleftsuper#1{\smash{\mathop{\longleftarrow}\limits^{#1}}}
\def\maprightsub#1{\smash{\mathop{\longrightarrow}\limits_{#1}}}
\def\mapleftsub#1{\smash{\mathop{\longleftarrow}\limits_{#1}}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ABBREVIATED GREEK
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\def\p{\partial}
\def\ti{\tilde}
\def\varep{\varepsilon}
\def\p{\partial}
\def\al{\alpha}
\def\Be{\Beta}
\def\be{\beta}
\def\de{\delta}
\def\De{\Delta}
\def\ep{\epsilon}
\def\et{\eta}
\def\ga{\gamma}
\def\Ga{\Gamma}
\def\io{\iota}
\def\ka{\kappa}
\def\la{\lambda}
\def\La{\Lambda}
\def\na{\nabla}
\def\om{\omega}
\def\Om{\Omega}
\def\ph{\phi}
\def\Ph{\Phi}
\def\ps{\psi}
\def\Ps{\Psi}
\def\rh{\rho}
\def\si{\sigma}
\def\Si{\Sigma}
\def\ta{\tau}
\def\th{\theta}
\def\Th{\Theta}
\def\Up{\Upsilon}
\def\ze{\zeta}
%
%THIS NEXT MACRO REPLACES THE TEX \bar COMMAND
%\def\baroverletter#1{\setbox1=\hbox{$#1$}
% \dimen1=\ht1
% \advance\dimen1 by 1.25pt % was 1pt
% \dimen2=\wd1
% \advance\dimen2 by -2pt % was -1pt
% \rlap{\hspace{1pt}\rule[\dimen1] % was .5pt
% {\dimen2}{.35pt}}\box1} % was .25pt
%
\def\baroverletter#1{\setbox1=\hbox{$#1$}
\dimen1=\ht1
\advance\dimen1 by 1pt
\dimen2=\wd1
\advance\dimen2 by -1pt
\rlap{\hspace{.5pt}\rule[\dimen1]
{\dimen2}{.35pt}}\box1} % was .25pt chgd. 2/11/89 whm
%
%%%%
\def\physifs{Department of Physics and\\ Institute for Fusion Studies\\ The
University of Texas at Austin\\ Austin, Texas\ \ 78712}
%\input boldchar
%\textwidth 6.5in
%\oddsidemargin 0in \evensidemargin 0in
%\topmargin -.5in
%\textheight 8.5in
\begindoc
\bibliographystyle{unsrt}
\title{Self-Consistent Chaos in the Beam-Plasma Instability}
\author{
{\bf J.L. Tennyson}
\thanks{Posthumous. Prepared by JDM and PJM.}
\\Stanford Linear Accelerator Center\\Stanford University \\Palo Alto, CA
94305\\
{\bf J.D. Meiss}
\\Program in Applied Mathematics\\ Box 526, University of Colorado\\
Boulder CO 80309\\
{\bf P.J. Morrison}
\\ \physifs\\}
\date{January 21, 1993}
\baselineskip 12 pt %set single spacing for title
\maketitle
\newpage
\begin{abstract}
\bskip24pt
\baselineskip 14 pt %set spacing for abstract
The effect of self-consistency on Hamiltonian systems with a large number of
degrees-of-freedom is investigated for the beam-plasma instability using
the single-wave model of
O'Neil, Winfrey, and Malmberg. The single-wave model is reviewed and
then rederived within the Hamiltonian context, which leads naturally to
canonical
action-angle variables. Simulations are performed with a large ($10^4$)
number of beam
particles interacting with the single wave. It is observed that the system
relaxes into a
time asymptotic periodic state where only a few collective degrees are
active; namely, a clump of
trapped particles oscillating in a modulated wave, within a uniform chaotic
sea with
oscillating phase space boundaries. Thus self-consistency is seen to
effectively {\sl
reduce} the number of degrees-of-freedom. A simple low degree-of-freedom
model is derived
that treats the clump as a single {\it macroparticle}, interacting with the
wave and
chaotic sea. The uniform chaotic sea is modeled by a fluid waterbag, where
the waterbag
boundaries correspond approximately to invariant tori. This low
degree-of-freedom
model is seen to compare well with the simulation.
\end{abstract}
\setlength{\baselineskip}{16pt}
\newpage
\section{Introduction} \label{sec:introduction}
Chaotic motion in Hamiltonian systems is common whenever there are more than one
degree-of-freedom \cite{MacKay-Meiss}. Often, the systems studied are low
dimensional
approximations of many degree-of-freedom systems. In some cases, such as
planetary dynamics,
highly accurate descriptions can be given with only a few
degrees-of-freedom. However,
there are many situations, such as galactic dynamics, where the number of
degrees-of-freedom is essentially infinite. Generally, one expects such
systems to
exhibit greater chaos when the dimension increases; furthermore, some
effects such as ``Arnold diffusion" are possible only with a larger number
of degrees-of-freedom.
Even in high dimensional cases it is often of interest to study a low
dimensional approximation, for example to study the motion of a single star
in a given galactic
gravitational potential---this was the motivation for the classic study of
H\'enon and
Heiles (see references in \cite{MacKay-Meiss}). Such an approximation is not
{\it self-consistent}. Other well studied examples of this type include the
motion of
charged particles in electromagnetic fields, where the fields produced by
the particles are
ignored; the motion of tracer particles in a fluid, where the influence of
these particles on
the fluid velocity field is ignored (the passive advection problem); and
the so-called
``kinematic" dynamo, where a velocity field can intensify a magnetic field,
but the
back-reaction of the field is ignored.
There has been little work on the effect of self-consistency. In this paper
we show
how it is possible in a system with a large number of degrees-of-freedom for the
inclusion of self-consistency to result in dynamics with ``effectively" few
degrees-of-freedom.
The model considered here treats self-consistency in the simplest possible
way. A large
number of noninteracting particles are influenced by a time dependent
potential. The
Hamiltonian for each particle has one and one half degrees-of-freedom, and
so the
motion can be chaotic. However, each of the particles is charged and
therefore contributes to
the potential---this is the self-consistent effect. In the model, we assume
that the
field has a single degree-of-freedom; i.e., that the spatial dependence of
the field is
given. Thus each particle experiences forces due to the others, but only to
the extent
that the other particles contribute to the single mode of the field. This is in
contrast to the fully self-consistent n-body dynamics, where each particle
interacts
directly with every other particle. This latter case is considerably more
difficult.
Models similar the one described above may be appropriate for many physical
situations; for example, a galaxy with a predominantly axisymmetric
gravitational potential that is perturbed by a small number of modes, say
those corresponding to spiral density waves. Each star contributes to
these modes, and also has a possibly chaotic motion in the corresponding
field. Similar effects
occur for planetary rings, beam-beam interactions in accelerators, tearing
mode-plasma
interactions in tokamaks, etc.
The specific problem we consider is the beam-plasma instability. The
formulation is due to
O'Neil, Winfrey and Malmberg (hereafter referred to as {\it OWM})
\cite{OWM}, see also \cite{Shevchenko}.
Here a beam of charged particles moves in a background neutral plasma. The
system is unstable to
the formation of electrostatic plasma waves. Following \cite{OWM} we
suppose that the
most unstable of these waves predominates, and neglect the remainder of the
modes. OWM
showed that the wave grows in amplitude until it traps the beam particles.
It then
saturates and begins to oscillate in amplitude as the beam particles slosh
in the wave
potential. In a later study Mynick and Kaufman \cite{Mynick-Kaufman}
developed a simple
model to account for these oscillations---they assumed that the beam could
be represented
by a rigid bar in phase space. When the beam is trapped in the wave, the
rigid bar begins
to rotate. Mynick and Kaufman computed the frequency shift and amplitude
oscillations of
the wave as it interacts with this rigid bar. Smith and Pereira
\cite{Smith-Pereira}
noted that since the amplitude of the plasma wave oscillates, the beam
particles can
experience chaotic motion. They studied the motion of a test particle in a
model of this
oscillating field and showed that much of the test particle phase space is
indeed
chaotic. However, there is an island in the phase space where the motion is
regular; they
noted that some fraction of the beam particles in the numerical experiments
of OWM
should find themselves in the correct region of phase space to be trapped
in this
oscillating island. Later Adam, Laval and Mendonca \cite{Adam} studied a
model in which a
single bunch of particles, what we will call a {\it macroparticle},
interacts with the plasma wave. As we will show below using the Hamiltonian
formulation, this
two degree-of-freedom system is integrable due to the existence of a phase
symmetry. In
\cite{Adam} it was shown that the macroparticle system has solutions which
correspond to
periodic oscillations of the bunch in the wave.
In \S{sec:model} we review the derivation of the OWM model, deriving the
Hamiltonian
formulation of \cite{Mynick-Kaufman} directly from the Vlasov-Poisson
equations.
\S{sec:simulations} discusses numerical solutions of the OWM equations with
up to $10^5$ beam
particles. As was observed in \cite{Shevchenko}, the simulations show that
the linear instability saturates and oscillates about a finite amplitude.
We observed at least $100$ periods of these oscillations; as far as we can
determine, the oscillations persist and the system becomes asymptotically
periodic. We study the motion of a single test particle in this periodic
potential, showing that a substantial portion of the original beam is
indeed trapped in a
stable island in the test particle phase space. However more than half the
beam finds itself in
the chaotic region of phase space, and spreads more or less uniformly over
this region. The
upper and lower boundaries of this ``chaotic sea" are formed from invariant
tori of the test particle
system.
In \S{sec:sea} we construct a four degree-of-freedom model. One
degree-of-freedom is the
wave, the second corresponds to the trapped bunch of beam particles---the
macroparticle, and the
last two describe the chaotic sea. These two degrees-of-freedom describe
the oscillations of the
boundary of the chaotic sea and are
derived from the ``waterbag" approximation. A waterbag consists of a
constant phase space
density between two moving boundaries. In our case the simulations show
that the phase space
density of the chaotic particles is indeed nearly constant and the
boundaries of the chaotic zone are
formed from invariant surfaces well outside the oscillating separatrix of
the wave. We model these
boundaries with sinusoidal curves, an assumption consistent with that of the single mode in
the potential. Finally, the frequency shift of the trapped particle
oscillations due to the
chaotic sea is computed analytically from this model.
\section{Single Wave Model} \label{sec:model}
O'Neil, Winfrey, and Malmberg (OWM) \cite{OWM} introduced the single
wave model to describe the growth and saturation of the weak beam-plasma
instability. In this section we briefly review the approximations
made in their derivation, rederive the model within the Hamiltonian context and
obtain the Hamiltonian form previously presented by Mynick and Kaufman
\cite{Mynick-Kaufman}, discuss linear instability, and finally consider a
special case
where only a single beam particle is included.
\subsection{Derivation} \label{sec:derivation}
To obtain the single wave model, the response of the plasma and beam to
the fields are treated separately. We consider only the one dimensional,
collisionless,
nonrelativistic, electrostatic case. The total electron density
\begin{equation}
\label{density}
n(x,t) = n_p(x,t)\ +\ n_b(x,t)\ =
\ n_p(x,t)\ +\ \sum\limits_{ j = 1}^{N} \delta (x-x_j(t))
\end{equation}
is a sum of contributions from the plasma and from a set of $N$ beam
particles at positions $x_j(t)$. Beam particles
evolve according to the fully nonlinear force equation
\begin{equation}\label{newton}
m \ddot{x}_j = - e E (x_j,t)\,,
\end{equation}
where the electron charge is $-e$.
A basic assumption of the model is that the phase velocity of the resulting
instability is much larger than the velocities of particles in the
background plasma:
the plasma responds non-resonantly, and trapping effects of plasma particles
in the wave can be neglected. This implies that the plasma responds
approximately linearly to the wave so that we can define a linear hermetian
dielectric operator $\hat{\epsilon}$ such that
\begin{equation}
\label{dielectric}
4 \pi e n_p(x,t) = (1\ -\ \hat{\epsilon }) \varphi '' (x,t)\,,
\end{equation}
where $\varphi$ is the electrostatic potential, $E = -\varphi'$.
Substituting this into Poisson's equation for $\varphi$ yields
\begin{equation}\label{poisson}
\hat{\epsilon} \varphi = 4 \pi e n_b(x,t)\,.
\end{equation}
For a homogeneous plasma, $\hat{\epsilon}$ is a convolution operator in
both $x$ and $t$, and Poisson's equation is most easily treated by Fourier
transform. In this representation the dielectric becomes a real function,
$\epsilon(k,\omega)$. For a weak beam, the right hand side of
\eqq{poisson} is ``small" implying that the electrostatic response is
dominated by the zeros, $\epsilon(k,\omega_0(k)) = 0$, of the dielectric.
It is a reasonable approximation to expand $\epsilon$ about one
such zero retaining only the first derivative of $\epsilon$ with
respect to $\omega$:
\begin{equation}
\epsilon(k,\omega)\ \approx \ \epsilon (k,\omega_0)\ + \\
\left.{\frac{\partial
\epsilon}{\partial \omega }}\right| _{\omega = \omega_0} (\omega-\omega_0) =
\epsilon '(\omega -{\omega }_0) \,.
\end{equation}
For example, for a cold plasma $\epsilon = 1 - \omega_p^2 / \omega^2$,
and $\partial \epsilon / \partial \omega |_ {\omega = \omega_0}
\equiv \epsilon ' = 2/ \omega_p$.
Transforming back to the time domain and using \eqq{poisson} then gives
\begin{equation}
\label{phi}
\dot E_k + i\omega_0 E_k = \frac{4\pi e}{k L \epsilon'} \sum\limits_{j=1}^{N}
e^{-ik{x}_{j}(t)}\,,
\end{equation}
where we have used the Fourier transform of the beam density of \eqq{density}
\begin{equation}
n_b(k,t) = \frac{1}{L} \int_{0}^{L}dx\ e^{-ikx}{n}_{b}(x,t)\ =
\frac{1}{L} \sum\limits_{j = 1}^{N} e^{-ik{x}_{j}(t)}\,,
\end{equation}
with L defining the periodicity length of the system.
At this point we assume that the electrostatic field has only one spatial
Fourier component, i.e., that there is a single wave. It is expected that the
single wave approximation is appropriate when the width of the unstable
spectrum in $k$-space is relatively narrow in units of $2\pi/L$. In this
case, if $k$ represents the most unstable mode, the amplitude of all
other Fourier components will be exponentially smaller than the single
wave during the linear growth stage. Of course, some time after nonlinear
saturation of the instability, the other modes will become important. As
is well known, the width of the unstable spectrum depends on the small
parameter $(n_b/n_p)^{1/3}$, so that the single wave model will be most
appropriate in the weak beam case.
Under the single wave assumption \eqq{newton} yields
\bqy
\label{hamilton}
\dot x_j&=&p_j/m \nonumber\\
\dot p_j&=&-e\left(E_ke^{ikx_j}+E_{-k}e^{-ikx_j}\right)\,.
\eqy
\Eq{hamilton} together with \eqq{phi} are the closed dynamical system
that governs the interaction of a single wave with the beam particles.
\subsection{Hamiltonian Structure and Derivation} \label{sec: hamiltonian}
Now consider the derivation of the equations of motion, \eqs{phi} and
(\ref{hamilton}), within the Hamiltonian context. The derivation proceeds
by first
considering the kinematics, i.e.~the dynamical variables used to
describe the state of the system, and then the dynamics,
obtained by finding the appropriate Hamiltonian.
We begin the first part by supposing that the electrons are described by
specifying their phase space coordinates $(x_j,p_j)$, where $j=1,2,...M$.
The first
$N(0, \psi =\pi)$ are less intuitive. These exist
only if $P> 3 N_m\left(\frac{N_m}{4N}\right)^{1/3}$. They correspond to
a particle perched on the top of the potential well. The lower momentum
particle is unstable, while the larger momentum particle is actually
stable. This can be seen in \fig{fig:singleparticle} which shows contours of
constant $H$ in the $(p,\psi)$ phase space.
Small oscillations about the stable fixed point are such that as the
particle starts to fall off the potential hill, the wave phase shifts so as
to catch the particle.
The final ``fixed point" is the degenerate case $(p = P, \psi =
\rm{arbitrary})$ for which the wave amplitude is zero. This corresponds
to the unstable equilibrium from which the beam plasma instability
results, as discussed in the previous section.
Small oscillations about these fixed points, with $(J_0=P-p_0$ are
governed by the linearized Hamiltonian
\begin{eqnarray}\label{linearHam}
\delta^2 H &=& \frac{\delta p^2}{2M_e}
+ N_m \left(\frac{J_0}{N}\right)^\half \delta \psi^2
\,; \nonumber\\
M_e &\equiv& N_m \frac{2 J_0}{2 J_0 - p_0} \,.
\end{eqnarray}
Here $M_e$ is the effective mass due to the self-consistent coupling.
Note that since $p_0<0$, the effective mass is smaller than that of a
particle in a fixed potential well. This gives a bounce frequency of
\begin{equation}
\omega_B^2 = 2 \left(\frac{J_0}{N}\right)^\half \left( 1-\frac{p_0}{2J_0}
\right)\,.
\end{equation}
The first factor is just the bounce frequency of a particle in a fixed well,
normalized in accord with \eqq{nondimensional}. The latter term provides
an increase in the frequency, as also discussed in \cite{Mynick-Kaufman}.
To compare these with the simulations in the next section, we set the total
momentum, $P$, to
zero, corresponding to the initial conditions $J(0) = p(0) = 0$.
In this case $J(\tau) = -p(\tau)$ and there is only one fixed point
\bq
\Omega \equiv \frac{p_0}{N_m} = - \left(\frac{N_m}{N}\right)^{1/3} \,,
\quad \psi = 0 \,.
\eq
>From the simulations we will find $N_m \simeq .4N$, and thus that
\begin{eqnarray} \label{zero_p}
\Omega &=& -0.74 \nonumber \\
|\Phi| &=& 0.29 \nonumber \\
\omega_B &=& 0.94\,.
\end{eqnarray}
\section{Simulations} \label{sec:simulations}
Simulations of the model of \eqs{motion} were carried out for several
initial conditions and a number of different integration algorithms. The results
shown below were obtained with a symplectic, leap frog method.
\subsection{OWM Model} \label{sec:numerics-OWM}
In the simulations the particles are initialized as a cold beam with
uniform spacing $\xi_j =
2\pi j/N \ , j = 0,...N$, and zero momentum $p_j = 0$. The initial wave
amplitude is set to be
very small. In \fig{fig:amplitude}, the wave amplitude is shown as a
function of time for
\eqq{motion} with $N=10^4$. For small $\tau$ the wave amplitude
grows exponentially in time at the rate predicted by \eqq{growth_rate} and
at the phase
velocity $v_{\phi} = {\cal R}(e^{2\pi i/3}) = -0.5$. %
%
% Jeff uses physical dimensions. The scaling between Jeff's frequencies and
ours is
% \omega_Jeff = 0.12599 \omega
%
As the wave grows, the beam experiences a growing sinusodal perturbation,
and as can be seen in
the density plot of \fig{fig:density}, the beam density also varies
sinusoidaly. Near $\tau = 16$ the
beam curls over as the particles begin to oscillate in the wave.
Consequently the amplitude of the
wave reaches a maximum.
At this stage the beam density, $n_b(\xi,\tau)$, develops cusps at the positions
where the beam curls over, see \fig{fig:density}. Note that though there
are many spatial
harmonics in the beam density, the single wave model does not allow the
development of similar
harmonics in the potential. These would lead to the growth of
other waves and undoubtedly greatly change the subsequent behavior of the
system.
None-the-less, the subsequent development of the OWM dynamics is quite
interesting. As the
beam particles begin to oscillate in the wave, their oscillation
frequencies depend upon
their energy, just as for a single particle in a fixed potential. Thus as
the beam begins to
rotate about the potential minimum, those particles closer to the center
have larger
oscillation frequencies than those near the ``separatrix."
If the wave amplitude were fixed, one would see phase mixing of the
particles (visualized
as an ever tighter spiral in the particle phase space), and the
oscillations in the
particle total energy would damp away---this is the mechanism of Landau
damping in a
large amplitude wave discussed by O'Neil \cite{O'Neil}.
However since $v_\phi \approx -0.5$ and the beam is
initialized at $v = 0$, when the beam particles oscillate in the wave, their
net momentum also oscillates. Hence, because of the conservation of total
momentum,
\eqq{conservation}, the wave momentum, $J$, must also oscillate as well.
Therefore the wave
amplitude is not fixed and each beam particle experiences an oscillating
potential. As
is well known, the phase space for a single beam particle in a given
oscillating potential has
chaotic zones, especially for the regions near the oscillating separatrix.
In the beam
particle phase space, \fig{fig:phase1}, we show the phase space positions
of each of the beam
particles at a fixed time. Note the two distinct regions: one chaotic in
appearance with a nearly
uniform density, and the other a more coherent cluster of particles. In the
cluster one still sees
evidence of the initial beam, though it has wrapped around a number of times.
In the simulations, which were carried out up to $\tau \approx 650$, the
oscillations in the potential persisted, and indeed, as can be seen in
\fig{fig:amplitude} the oscillations become increasingly periodic with
time. Furthermore
as we varied $N$ up to $10^5$ and the changed integration accuracy, we
noted that these
oscillations became more periodic and constant in amplitude as the number
of particles increased
and as the accuracy improved. Thus we believe that the asymptotic state is
a periodic one.
Meanwhile, the particle phase space exhibits quite complicated behavior.
About 60\% of the
particles---those with relatively large energies in the wave
frame---experience chaotic
motion, and spread out roughly uniformly in a region of phase space whose
average width
is $\Dw = 4.7$. We called the chaotic particles the ``stochastic sea" in
1983---a more modern terminology is the ``chaotic sea". The remaining 40\%
of the particles
are clustered tightly together in a clump which has a width of order
$\pi/2$ and a density
about 6 times larger than the background, recall \fig{fig:density}. These
particles continue
to oscillate coherently in the wave, as seen in the sequence of phase
space pictures,
\fig{fig:phaseseq}. Also shown in the figures is the instantaneous
separatrix of the wave---that is
the separatrix that a single beam particle would see in a fixed potential
at the current
amplitude. Note that wave amplitude and the momentum of the cluster
oscillate $180^{\circ}$ out
of phase. This clump of particles is treated as a single particle, the
macroparticle, in the model of
\S{sec:sea}.
In addition to the uniform cold beam, several different initial conditions
have been partially
studied. For example we started particles on a circular ring, $\xi^2 + p^2
= r^2$, and also
considered a bunched beam with particles placed uniformly in a limited
range, $ -1 <\xi < 1$ with
$p \ne 0$. In these cases some subset of the particles remained coherently bunched as a
macroparticle, and the system did not appear to settle into an equilibrium.
Cold beams with nonzero momenta also lead to oscillations as was shown in
\cite{Shevchenko}, though the amplitude of the .We have not investigated
this in any generality and do not know the extent of the initial conditions
which will give rise to a periodic final state.
\subsection{Test Particle} \label{sec:test-particle}
To investigate further dynamics of the beam particles,
consider the ``test particle" motion of a single particle in a given oscillating
potential. These are obtained from the nonself-consistent, one and one half
degree-of-freedom Hamiltonian
\bq \label{testH}
H_t(p, \xi,\tau) = \half p^2 - 2 \left(\frac{J(\tau)}{N}\right)^\half
\cos{(\xi-\theta(\tau))}\,,
\eq
where $J$ and $\theta$ are considered to be given periodic functions of
$\tau$. This
system has been studied in the same context by Smith and Pereira
\cite{Smith-Pereira}, who
used a simple model for $J(\tau)$ and $\theta(\tau)$.
Here we determine $J$ and $\theta$ numerically, from the simulations of
\S{sec:numerics-OWM},
building these functions from an average over a number of periods of the
oscillations.
A stroboscopic plot of the test particle dynamics is shown in
\fig{fig:test_particle}
for several different values of $\theta$. The dots represent the
trajectories of a
number of different test particles. As was also noted in
\cite{Smith-Pereira}, there is a
prominent stable island in the test particle phase space which oscillates
exactly out of
phase with the potential; much of the rest of the phase space is chaotic.
Also shown in
\fig{fig:test_particle} are the corresponding plots of the beam particle phase
space---here of course each point in the plots is the position of one of
the $10,000$
beam particles. Note that the macroparticle clump sits, as near as can be
ascertained,
at the position of the test particle island. This verifies an assertion in
\cite{Smith-Pereira}, where it was merely noted that some fraction of the
beam particles
initially stretched across the position of the test particle island.
\section{Chaotic Sea Model} \label{sec:sea}
As we have seen from the simulations, the asymptotic state of the cold
beam initial
condition, evolved under the OWM Hamiltonian, appears to be almost exactly
periodic.
Approximately 40\% of the initial beam forms a the clump of particles that
oscillates in the
potential well formed by the wave. The remaining particles spread out
approximately uniformly
in phase space between two oscillating boundaries, which are approximately
invariant
surfaces of the test particle Hamiltonian $H_t$.
In this section we use these ideas to obtain a reduced Hamiltonian model of four
degrees-of-freedom, which approximately describes the full 10,001
degree-of-freedom system.
In the model, as noted above, we assume that the clump of regularly
oscillating particles is
localized enough so that all these particles can be treated as one
located at $(\xi, p)$. This macroparticle contains $N_m$ particles
and hence a mass $mN_m$ and charge $-eN_m$. The approximation that the $N_m$
particles can be treated as a single particle ignores any
internal degrees-of-freedom of the cluster (of which there are $N_m-1$). It
is clear that a
correction to the model would include, at the least, some sort of
rotational degree-of-freedom;
this would be similar to the analysis of Mynick and Kaufman
\cite{Mynick-Kaufman}
who assumed that the cluster of particles formed a ``bar" in phase space
(an assumption that is
reasonable only for the first few trapping oscillations). In our case the
particles in the
macroparticle could be thought of as rotating on the invariant surfaces of
the test
particle Hamiltonian, \eqq{testH}, about the central periodic orbit. This
behavior could be
modelled by the addition of an oscillator degree-of-freedom to the system;
however, we neglect
this motion here.
Much more interesting is the treatment of the chaotic particles. As noted in the
\S{sec:numerics-OWM}, the phase space density of these particles appears to
be nearly
constant up to fairly sharp boundaries in velocity space. We assume that
these particles can
be treated as a continuum with a constant phase space density $f_c$,
between the boundaries
$v_+(x,t)$ and $v_-(x,t)$; thus,
\begin{equation}\label{ssdensity}
n_c (x,t) = \int_{v_-}^{v_+} f_c dv = f_c (v_+ - v_-)
\end{equation}
is the density of the chaotic sea. The total number of such particles in
the length $L$ will
be denoted by $N_c = N-N_m$. Particles in the chaotic sea evolve according
to \eqq{newton},
and hence $f_c$ evolves according to the Vlasov equation. As is well known,
and easy to see, the Vlasov equation in this situation can be exactly
reduced to two equations for
the evolution of the boundaries \cite{Davidson}. These equations are called
the waterbag
equations:
\begin{eqnarray}\label{waterbag}
{\frac{\partial v_+}{\partial t}}\ +\ v_+{\frac{\partial v_+}{\partial x}}
&=& -eE
\nonumber\\
{\frac{\partial v_-}{\partial t}} \ +\ v_-{\frac{\partial v_-}{\partial
x}} &=& -eE\,.
\end{eqnarray}
Following the philosophy of the derivation of the OWM model, where only a
single wavenumber
in
$\varphi$ is kept, we retain only one wavenumber in $v_\pm(x)$:
\begin{equation}
v_\pm = v_\pm^0 + \tilde{v}_\pm e^{ikx} + \tilde{v}^*_{\pm} e^{-ikx} \,.
\end{equation}
The equations of motion then become
\begin{equation}
\left( \frac{\partial}{\partial t}+ ikv_\pm^0 \right) \tilde{v}_\pm = -e
E_k \,,
\end{equation}
where $v_\pm^0$ are the mean velocities (which do not change in time according
to \eqq{waterbag}). Note that this implies that the boundary of the chaotic
sea, an invariant surface of $H_t$, is sinusoidal in shape. We see from
\fig{fig:phaseseq} that this is reasonable since the boundary deviation is
small and
it is not too close to the ``separatrix" of the wave.
The evolution of the field is obtained from Poisson's equation, as
previously, except that the density of the chaotic sea, as given by
\eqq{ssdensity}, must be
included. Thus \eqq{phi} becomes
\begin{equation}
\dot{E}_k +i\omega_0 E_k = \frac{4\pi e}{k\epsilon'}
\left[
f_c (\tilde{v}_+ - \tilde{v}_-) + \frac{N_m}{L} e^{-ikx_m}\,.
\right]
\end{equation}
These equations are nondimensionalized, using \eqq{nondimensional}; in addition
defining
\bq
\omega_\pm \equiv \frac{kv_\pm^0 - \omega_0}{\omega_b} \,, \quad
V_\pm \equiv \frac{k}{\omega_b}\tilde{v}_\pm e^{i\omega_0 t} \,.
\eq
In terms of these variables the equations of motion become
\begin{eqnarray}
\dot{\Phi} &=& i \frac{N_c}{N \Dw} (V_{+} - V_-) +i\frac{N_m}{N}
e^{-i \xi}
\,,\nonumber\\
\ddot{\xi} &=& i \Ph e^{i \xi } - i\Ph^* e^{-i \xi }\,,\nonumber\\
\dot{V}_\pm &=& -i \omega_\pm V_\pm + i \Ph\,,
\end{eqnarray}
where $ \Dw = \omega_+ - \omega_- $ is the average, nondimensional width of
the chaotic sea.
This set of equations is also a Hamiltonian system, with the wave
action-amplitude variables
defined in \eqq{action}, and the new action-amplitude variables for the
chaotic sea defined by
\begin{equation}
V_+ = \left(\frac{J_+ \Dw}{N_c}\right)^\half e^{-i\theta_+} \,, \quad
V_- = \left(\frac{J_- \Dw}{N_c}\right)^\half e^{+i\theta_-} \,,
\end{equation}
which results in the Poisson bracket relations
\begin{equation}
\left[ V_\pm^* , V_\pm \right] = \pm i \frac{\Dw}{N_c} \,.
\end{equation}
The Hamiltonian takes the form
\bqy
H\ =\ {\frac{{p}^{2}}{2{N}_m}}&+&
\frac{N_c \omega_+}{\Dw} {\left|{V_+}\right|}^2 -
\frac{N_c \omega_-}{\Dw} {\left|{V_-}\right|}^2 \nonumber \\
&-&\left\{\frac{N_c}{\Delta \omega} \Ph^*(V_+ - V_-) + N_m \Ph^* e^{-i\xi
} + \rm{c.c.}\right\}\,.
\eqy
In terms of canonical variables this becomes
\begin{eqnarray}\label{ssHam}
H = \frac{p^2}{2N_m} &+& \omega_+ J_+ - \omega_- J_- \nonumber \\
&-& 2 \alpha \Big[\sqrt{J J_+} \cos{(\theta - \theta_+)} -
\sqrt{J J_-} \cos{(\theta +
\theta_-)}\Big] \nonumber \\
&-& 2 \beta \sqrt{J} \cos{(\theta-\xi)} \,.
\end{eqnarray}
where the coefficients are given by
\bq
\alpha = \left( \frac{N_c}{N \Dw} \right)^{1/2} \,, \quad \beta =
\frac{N_m}{\sqrt{N}} \,.
\eq
The first three terms in the Hamiltonian represent the total particle
kinetic energy as a
sum of macroparticle energy and harmonic terms for the oscillations of
the chaotic sea boundary. The last three terms give the electrostatic
interaction energy.
Thus we have reduced the 10,001 degree-of-freedom Hamiltonian to one of
four degrees-of-freedom, which describes well the motion in the periodic
final state of the
simulations, {\sl provided} the three parameters $\omega_+, \omega_-$ and
$N_m$ are given.
\subsection{Symmetry, Equilibrium, and Stability} \label{sec:symmetry}
The chaotic sea Hamiltonian conserves the total momentum, given by
\begin{equation} \label{conservation2}
P = J + p + J_+ - J_- + \frac{N_c}{2}(\omega_+ + \omega_-)\,.
\end{equation}
Here the first two terms represent the momentum in the wave and
macroparticle, and the
last three are the contributions of the chaotic sea. These latter terms
include the momentum
in the oscillations of the waterbag boundary, $J_+ -J_-$, and the mean
momentum in the
chaotic sea represented by the last term in \eqq{conservation2}. This
equation can also be
derived from \eqq{conservation} by splitting off the contribution of the
particle momentum for
the $N_c$ particles in the chaotic sea, treating them as a Vlasov fluid in
the waterbag.
As far as we know, there are no other conserved quantities besides
the energy. Thus the system has three effective degrees-of-freedom, and
should exhibit the
full complexity of such systems, including chaotic motion and Arnold
diffusion.
The stable equilibrium of \eqq{ssHam} corresponding to the macroparticle
sitting in the potential
minimum is given by the equations
\begin{eqnarray} \label{zero_p_mod}
\Omega \equiv \frac{p}{N_m} &=& \dot{\xi} = \dot{\theta} = \dot{\theta}_+
= - \dot{\theta}_-
\nonumber\\
\Omega &=& \alpha^2 \left( \frac{1}{\Omega - \omega_+} - \frac{1}{\Omega -
\omega_-} \right)
- \frac{\beta}{\sqrt{J}} \nonumber \\
J_\pm &=& \frac{\alpha^2}{(\Omega - \omega_\pm)^2} J \,.
\end{eqnarray}
>From the simulation results of \fig{fig:phaseseq}, we found $ N_m \approx
0.4 N$, $\omega_+
\approx 1.9$, and $\omega_- \approx -2.8$. This gives $\alpha = .36$,
$\beta = 40$, and $\Dw =
4.7$. For these parameters we can solve \eqq{zero_p_mod} to obtain
\begin{eqnarray}\label{ssparams}
\Omega &=& -0.66 \, , \quad |\Phi| = 0.73 \nonumber \\
|V_+| &=& 0.29 \, , \quad |V_-| = 0.35 \,.
\end{eqnarray}
On the other hand, using \fig{fig:phaseseq} to determine the average phase
velocity of the wave in
the simulations we obtain $\Omega \approx -0.67$. The average value of
$|\Phi|$ from
\fig{fig:amplitude} is $0.75$. Both of these values
certainly agree with \eqq{ssparams} more closely than they do with the
single particle equilibrium
\eqq{zero_p}. In fact, if we take the measured value of $\Omega$ as given,
then
\eqq{zero_p_mod} gives the macroparticle fraction as 42\% and $|\Phi|$
becomes $0.74$, both of
which are consistent with our estimates from the numerics.
To compute the frequency of small oscillations about the equilibrium, we
first eliminate the action
$J$ using the conservation law (\ref{conservation2}), defining phases $\psi
= \theta - \xi$, and
$\psi_\pm = \theta_\pm \mp \theta$, conjugate to the momenta $p$ and
$J_\pm$ respectively.
Upon
linearization the Hamiltonian becomes
\bq
\delta^2 H = \frac{1}{2} \delta {\bf p} \cdot {\bf M}^{-1} \cdot \delta
{\bf p}
+ \frac{1}{2} \delta {\mbox{\boldmath $\psi$}} \cdot {\bf K} \cdot \delta
{\mbox{\boldmath $\psi$}}\,,
\eq
where $\delta {\bf p} = (\delta p, \delta J_+, \delta J_-)$ and $\delta
{\mbox{\boldmath $\psi$}} =
(\delta \psi, \delta \psi_+, \delta \psi_-)$ are the deviations from
equilibrium. The matrix $\bf{M}$,
the effective mass matrix, is positive definite. The matrix $\bf{K}$, the
effective spring constant
matrix, turns out to be diagonal. In order that the three terms in $\bf{K}$ be positive, we must
assume that $\psi =\psi_+= 0$, while $\psi_- = \pi$. This is consistent
with the fact that the
lower boundary of the chaotic sea is observed to have a $180^{\circ}$
phase lag with respect to
the
upper boundary.
The frequencies of small oscillation are given by the square roots of
eigenvalues of the matrix ${\bf
K} {\bf M}^{-1}$. For the parameters of the simulation, the mass matrix is
diagonal to a good
approximation. The element $M_{11}^{-1}$ turns out to be identical to
$1/M_e$ of
\eqq{linearHam}; neglecting terms of order $J_\pm /J$, the other diagonal
elements are
\bq
M_{22}^{-1} = - \frac{1}{2} \frac{\Omega - \omega_+}{J_+} \,, \quad
M_{33}^{-1} = \frac{1}{2} \frac{\Omega - \omega_-}{J_-}
\eq
The matrix ${\bf K}$ is
\bq
{\bf K}= dia(2 \beta \sqrt{J} , \, 2 \alpha \sqrt{J J_+} , \, 2 \alpha
\sqrt{J J_+})\,.
\eq
Using the values obtained before for the equilibrium, we determined the
eigenvalues
numerically from the \underline{full} matrix. The three oscillation
frequencies are
\bq
\omega_B = (1.33, \, 2.61, \, 2.11)
\eq
The first of these corresponds to the large oscillations observed in the
simulations. The
eigenfunction of this mode corresponds primarily to the $\delta p$
degree-of-freedom that
describes oscillation of the macroparticle. From \fig{fig:amplitude} we
measure the frequency of
oscillations of the asymptotic state, obtaining $\omega_B = 1.28$---again
a good agreement
with the calculated value. We have no evidence for the excitation of the
other two normal modes;
however, it might be possible to determine these through careful simulation.
Thus we conclude that the chaotic sea model agrees with the simulations
extremely well,
and much better than the single particle calculation, \eqq{zero_p}.
\section{Conclusions}
We have studied, the long time dynamics of the electrostatic interaction of many
particles with a plasma wave. The wave arises from an instability (the
beam-plasma
instability) of the initial state corresponding to a cold beam of
particles. In the
simulations, the asymptotic state corresponds to a periodically oscillating
wave amplitude
together with a trapped clump of particles. About $42\%$ of the particles
are trapped
by approximate invariant surface within the oscillating wave, while the
remaining particles
move chaotically---becoming successively trapped and detrapped.
We modeled this motion by a four degree-of-freedom Hamiltonian system,
where one degree-of-
freedom
corresponds to the clump, two to the chaotic sea, and one to the wave.
This model quantitatively captures the asymptotic state of the effectively
infinite degree-of-freedom system.
One would like to speculate that there are other physical systems for which
the effect of
self-consistency would be similar. For example in the case of galactic dynamics,
the self-consistent propagation of a density wave would lead to chaotic
dynamics for many
stars, but there could be a set of stars which are coherently interacting
with the wave.
As usual a number of open questions remain:
\begin{itemize}
\item The periodic asymptotic state appears to arise for a number of
different initial
conditions in the OWM model. What is the ``basin" of initial conditions
which lead to such a
state? One could, for example, consider a warm beam initial state, or
widely different initial
states such as those discussed at the end of \S{sec:numerics-OWM}.
\item Is there a way of self-consistently calculating the quantities $N_m$,
and $\omega_\pm$,
which we took from the simulations, for the chaotic sea model?
\item For the single particle Hamiltonian, there is a second stable
equilibrium corresponding
to the particle at rest at the maximum of the potential. Is there a
periodic state of the many
particle system near this equilibrium?
\item In the OWM model, only a single Fourier harmonic of the electrostatic
field is kept.
What is the effect of adding additional harmonics?
\item The chaotic sea Hamiltonian should exhibit the full array of
possible Hamiltonian motions. Are there many particle states that correspond to these
motions?
\end{itemize}
\vskip .5in
\centerline{ACKNOWLEDGMENTS}
\par
Work on this project began in 1982 when all the authors were at the
Institute for Fusion Studies. Due to subsequent dispersal of the authors,
the paper was not completed until now. Support was provided by the
U. S. Dept. of Energy Contract No. DE-FG05-80ET-53088
to the University of Texas at Austin.
\pagebreak
\section*{Figure Captions}
\begin{enumerate}
\item Contours of $H$, \eqq{macroham}, in the single particle phase space
$(p,\ps)$. Units of
$p$ (and $P$) are $N_m\left(\frac{N_m}{N}\right)^{1/3}$. In the upper
figure P = 0 and there is
one fixed point; in the lower figure P = 2 and there are three fixed points.
\label{fig:singleparticle}
\item Plot of $|\Ph(\tau)|$, the normalized wave amplitude, for $N=10,000$
particles initialized as a
cold beam.
\label{fig:amplitude}
\item Plot of the beam density as a function of position. The sinusoidal
distortion of the density due
to the growing wave is shown in (a) at $\tau = 12.6$. By $\tau = 69.3$, in
(b), about $10$
bounces have occurred and the macroparticle has formed, represented here by
the sharp peak
around $\xi = 0$. The remaining, chaotic sea particles have a sinusoidal
density variation.
\label{fig:density}
\item Plot of the beam particle phase space at $\tau = 641$ showing a well
defined macroparticle
and chaotic sea.
\label{fig:phase1}
\item Sequence of beam phase space plots over one bounce period. Note that
the macroparticle
bounces coherently in the wave, and the wave amplitude and chaotic sea
boundaries also oscillate
periodically. \label{fig:phaseseq}
\item Phase space plots comparing the full simulation of
\S{sec:simulations} with the dynamics of
a test particle in a given time-dependent potential $\Ph(\tau)$ as
determined by the simulation.
Shown are several test particle initial conditions at three different times
during a cycle.
\label{fig:test_particle}
\end{enumerate}
%\pagebreak
\begin{thebibliography}{99}
\bibitem{MacKay-Meiss}
R.S. MacKay and J.D. Meiss, {\it Hamiltonian Dynamical Systems: a Reprint
Selection},
(Adam-Hilgar, 1987).
\bibitem{OWM}
T.M. O'Neil, J.H. Winfrey, and J.H. Malmberg, Nonlinear interaction of a
small cold beam
and a plasma, {\it Phys. Fluids} {\bf 14} (1971), 1204-1212.
\bibitem{Shevchenko}
N. G. Matsiborko, I. N. Onishchenko, V. D. Shapiro and V. I. Shevchenko, On
Nonlinear Theory of Instability of a Monoenergetic Electron Beam,
\newblock {\it Plasma Phys.} {\bf 14} (1972), 591-600.
\bibitem{Mynick-Kaufman}
H.E. Mynick and A.N. Kaufman, Soluble theory of nonlinear beam-plasma
interaction,
{\it Phys. Fluids} {\bf 21} (1978), 653-663.
\bibitem{Smith-Pereira}
G.R. Smith and N.R. Pereira, Phase-locked particle motion in a large
amplituded plasma wave,
{\it Phys. Fluids} {\bf 21} (1978), 2253-2262.
\bibitem{Adam}
J.C. Adam, G. Laval, and I. Mendonca, Time-dependent nonlinear Langmuir waves,
{\it Phys. Fluids} {\bf 24} (1981), 260-267.
\bibitem{Morrison}
P.J. Morrison, The Maxwell-Vlasov equations as a continuous Hamiltonian system,
{Phys. Lett.}
{\bf A80} (1980), 383-386.
\bibitem{Morrison-Pfirsch}
P.J. Morrison and D. Pfirsch, Dielectric energy versus plasma energy, and
Hamiltonian
action-angle variables for the Vlasov equation,
{\it Phys. Fluids} {\bf B4} (1992), 3038-3057.
\bibitem{Shadwick-Morrison}
B. Shadwick and P.J. Morrison, Neutral modes and negative energy modes,
preprint, Institute for
Fusion Studies Report, IFSR\# 546 (1993) University of Texas at Austin.
\bibitem{Kueny}
C. Kueny, Ph.d thesis, Nonlinear instability and chaos in plasma wave-wave
interactions, University of Texas at Austin, 1993.
\bibitem{O'Neil}
T.M. O'Neil, Collsionless damping of nonlinear plasma oscillations,
{\it Phys. Fluids} {\bf 8} (1965), 2255-2262.
\bibitem{Davidson}
R.C. Davidson, {\it Methods of Nonlinear Plasma Theory}, (Academic, New
York, 1972), pp. 45-47.
\end{thebibliography}
\end{document}