% This is a copy of the generic AMS article.tex, adapted for use with
% the Fields fic-l document class. [mjd,14-Mar-1995]
%
% The changes were:
%
% ---Changed amsproc to fic-l in the documentclass statement
%
% ---Removed current address and email addresses.
%
% ---Changed "authors" to "advisors" in the dedication
%
% ---Changed the bibliography
%
% ---Removed examples of extra list types
%
% ---Removed the table
\documentclass{fic-l}
\usepackage{graphicx}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\theoremstyle{definition}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{example}[theorem]{Example}
\newtheorem{xca}[theorem]{Exercise}
\theoremstyle{remark}
\newtheorem{remark}[theorem]{Remark}
\numberwithin{equation}{section}
% Absolute value notation
\newcommand{\abs}[1]{\lvert#1\rvert}
\begin{document}
\title{Introduction to Gyrokinetic Theory with Applications in Magnetic
Confinement Research in Plasma Physics}
% Information for first author
\author{W. M. Tang}
% Address of record for the research reported here
\address{Princeton Plasma Physics Laboratory\\ Princeton University\\
Princeton\\ New Jersey 08543}
\email{tang@pppl.gov}
\thanks{The author is grateful to Dr. Wei-li Lee for his clear
description of the basic schemes for particle-in-cell simulations and
to Dr. Greg Rewoldt for his invaluable assistance in assembling the
final version of this document.This work was supported by
U. S. Department of Energy contract DE-ACO2-76CH03073.}
% Information for second author
%\author{Author Two}
%\address{Mathematical Research Section\\ School of Mathematical
%Sciences\\ Australian National University\\ Canberra ACT 2601,
%Australia}
% General info
%\subjclass{Primary 54C40, 14E20; Secondary 46E25, 20C20}
\date{\today}
%\dedicatory{This paper is dedicated to our advisors.}
\begin{abstract}
The present lecture provides an introduction to the subject of
gyrokinetic theory with applications in the area of magnetic
confinement research in plasma physics -- the research arena from
which this formalism was originally developed. It was presented as a
component of the ``Short Course in Kinetic Theory'' within the
Thematic Program in Partial Differential Equations held at the Fields
Institute for Research in Mathematical Science (24 March 2004). This
lecture also discusses the connection between the gyrokinetic
formalism and powerful modern numerical simulations. Indeed,
simulation, which provides a natural bridge between theory and
experiment, is an essential modern tool for understanding complex
plasma behavior. Advances in both particle and fluid simulations of
fine-scale turbulence and large-scale dynamics have been enabled by
two key factors: (i) innovative advances in analytic and computational
methods for developing reduced descriptions of physics phenomena
spanning widely disparate temporal and spatial scales; and (ii) access
to powerful new computational resources. Excellent progress has been
made in developing codes for which computer run-time and problem size
scale well with the number of processors on massively parallel
processors (MPP's). Examples include the effective usage of the full
power of multi-teraflop (multi-trillion floating point computations
per second) supercomputers to produce three-dimensional, general
geometry, nonlinear particle simulations which have accelerated
advances in understanding the nature of turbulence self-regulation by
zonal flows. Although the current results from advanced simulations
provide great encouragement for being able to include increasingly
realistic dynamics to enable deeper physics insights into plasmas in
both natural and laboratory environments, it should be kept in mind
that there will remain limitations to what can be practically
achieved.
\end{abstract}
\maketitle
\section{Introduction}
Plasmas are ionized gases which are often referred to as ``the fourth
state of matter'' and comprise over 99\% of the visible universe.
They are rich in complex, collective phenomena and encompass major
areas of research including plasma astrophysics and fusion energy
science. Fusion is the power source of the sun and other stars, which
occurs when forms of the lightest atom, hydrogen, combine to make
helium in a very hot (100 million degrees centigrade) plasma. The
development of fusion as a secure and reliable energy system that is
environmentally and economically sustainable is a truly formidable
scientific and technological challenge facing the world in the
twenty-first century. As such, progress toward this goal requires the
acquisition of the basic scientific understanding to enable the
innovations that are still needed for making fusion energy a practical
realization. In this as well as other areas facing major scientific
challenges, it is well recognized that research in plasma science
requires the accelerated development of computational tools and
techniques that aid the acquisition of the scientific understanding
needed to develop predictive models which can prove superior to
empirical scaling. This is made possible by the rapid advances in
high performance computing technology which will allow simulations of
increasingly complex phenomena with greater physics fidelity.
Accordingly, advanced computational codes, properly benchmarked with
theory and experiment, are now generally recognized to be a powerful
new tool for scientific discovery. In the key area of turbulent
transport of the plasma, the development of the gyrokinetic formalism
and its subsequent implementation in advanced simulations have enabled
major progress
In a magnetically-confined plasma, the interplay between the complex
trajectories of individual charged particles and the collective
effects arising from the long-range nature of electromagnetic forces
leads to a wide range of waves, oscillations, and instabilities
characterizing the medium. As a result, there is an enormous range of
temporal and spatial scales involved in plasmas of interest. As
illustrated in Figure 1, the relevant physics can span over ten
decades in time and space. Associated processes include the
turbulence-driven (``anomalous'') transport of energy and particles
across a confining magnetic field, the abrupt rearrangements
(disruptions) of the plasma caused by large-scale instabilities, and
the interactions involving the plasma particles with electromagnetic
waves and also with neutral atoms. Many of the these phenomena
involve short length and time scales (nanoseconds and microns) while
others occur on long time scales (seconds and minutes) and length
scales on the order of the device size (meters). Although the
fundamental laws that determine the behavior of plasmas, such as
Maxwell's equations and those of classical statistical mechanics, are
well known, obtaining their solution under realistic conditions is a
scientific problem of extraordinary complexity. Effective prediction
of the properties of energy-producing fusion plasma systems depends on
the successful integration of many complex phenomena spanning vast
ranges. This is a formidable challenge that can only be met with
advanced scientific computation properly cross-validated against
experiment and analytic theory.
\begin{figure}[tb]
\begin{center}
\includegraphics[width=4 in]{TangFig1.eps}
\end{center}
\caption{Huge ranges in spatial and temporal scales present major challenges to
plasma theory and simulation.}
\label{fig1}
\end{figure}
In magnetic confinement fusion experiments, the plasma interacts
directly with the ``confining'' electromagnetic fields, which can come
from an external source and/or from currents produced within the
plasma. This can lead to unstable behavior, where the plasma rapidly
rearranges itself and relaxes to a lower energy state. The resultant
thermodynamically favored state is incompatible with the conditions
needed for fusion systems, which require that more power output be
generated than it takes to keep the hot plasma well confined.
However, the hot plasma state is naturally subject to both large and
small-scale disturbances (``instabilities'') which provide the
mechanisms for lowering its energy state. It is therefore necessary to
first gain an understanding of these complex, collective phenomena,
and then to devise the means to control them. The larger-scale
(``macro'') instabilities can produce rapid topological changes in the
confining magnetic field resulting in a catastrophic loss of fusion
power density. Even if these instabilities can be controlled and/or
prevented, there can remain smaller-scale (``micro'') instabilities
which prevent efficient hot plasma confinement by causing the
turbulent transport of energy and particles. In order to make
progress on these issues, researchers in this field have effectively
developed the requisite mathematical formalism embodied by gyrokinetic
theory to deal with the complexity of the kinetic electromagnetic
behavior of magnetically- confined plasmas.
\section{Basic Ordering}
The general ordering appropriate for dealing with particle motion in
strongly magnetized plasmas involves the assumption that the Larmor
radius (or gyro-radius), $\rho$, of the particles is small compared to
the spatial variation of the electromagnetic fields ($L_B$); i.e.,
$L_B \equiv |B/\nabla B|$, with $\epsilon_B \equiv |\rho/L_B| \ll 1$,
where $\rho = v/\Omega$, $v =$ thermal velocity, and $\Omega =$
gyrofrequency. In addition, when addressing low frequency, long
parallel wavelength phenomena, the following ordering is usually
adopted: $\omega < \Omega$, $k_\perp \rho < 1$, $k_\parallel <
k_\perp$, with $\mathbf{k}$ being the wave number. This leads to a
more tractable or simplified version of the Boltzmann-Maxwell system
of kinetic equations in which the key components of the particle
motion are the fast gyro-motion plus the slow guiding center motion.
In particular, the ordering of the terms in the Boltzmann Equation can be
represented as follows:
\begin{eqnarray}
\mbox{} &&\frac{\partial \hat F}{\partial t} + \mathbf{v} \cdot \nabla \hat F +
\frac{Ze}{m} \left( \mbox{\boldmath$\varepsilon$} + \frac{1}{c} \mathbf{v}
\times \mathbf{B} \right)
\nabla_v \hat F = C(\hat F,\hat F) \label{eq1} \\
\mbox{} && \ \ \ \ \ \ \ \ \ \ \ \delta \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \delta \ \ \ \ \ \ \ \ \ 1 \ \
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \delta \ \ \ \ \ \ \ \ \ \ \ \ (F) \nonumber \\
\mbox{} && \delta\,\ \ \ \ \ \ \ \ \ 1 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ 1 \ \ \ \ \ \ \ \ \ 1 \ \
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \delta \ \ \ \ \ \ \ \ \ \ \ \ (f) \nonumber
\end{eqnarray}
where $\hat F = F + f$ and $\mbox{\boldmath$\varepsilon$} =
-\nabla(\Phi_0 + \Phi)$, with $F$ and $f$ being respectively the
equilibrium and perturbed distributions and $\Phi_0$ and $\Phi$ being
respectively the equilibrium and perturbed potentials governing the
electric field, $C(\hat F,\hat F)$ being the collision operator, and
$\delta$ designating the smallness ordering parameter.
Accordingly, the gyrokinetic equation governing the equilibrium is
\begin{eqnarray}
\mathbf{v}\cdot\nabla F^{(0)} - \mathbf{v}\cdot [\mu \nabla B +
v_\parallel (\nabla \hat n) \cdot \mathbf{v}_\perp ] \frac{1}{B}
\frac{\partial}{\partial \mu} F^{(0)} \\
-\frac{Ze}{m} \nabla \Phi_0 \cdot \mathbf{v}_\perp \frac{1}{B}
\frac{\partial}{\partial \mu} F^{(0)} - \Omega
\frac{\partial}{\partial \phi} F^{(1)} = C[F^{(0)},F^{(0)}] \nonumber
\end{eqnarray}
where $\hat n$ is a unit vector in the unperturbed magnetic field
direction, and averaging over the gyrophase $\phi$ gives
\begin{eqnarray}
v_\parallel \hat n \cdot \nabla F^{(0)} = C[F^{(0)},F^{(0)}]
\end{eqnarray}
which implies that $F^{(0)}$ is a Maxwellian.
The equation governing the perturbed distribution function is:
\begin{equation}
\left( \mathbf{v}_\perp \cdot \nabla - \Omega \frac{\partial}{\partial
\phi} \right) h = 0
\end{equation}
where $h = f^{(0)} + F_M Ze \Phi / T$. The ``guiding-center''
coordinate system is given by
\begin{eqnarray}
\psi' &=& \psi - \frac{1}{\Omega} \hat n \times \mathbf{v}_\perp \cdot
\nabla \psi \\
\chi' &=& \chi - \frac{1}{\Omega} \hat n \times \mathbf{v}_\perp \cdot
\nabla \chi \\
\phi' &=& \phi
\end{eqnarray}
where $\mathbf{B} = B \hat n = B_p \hat \chi + B_t \hat \zeta$, with
$B_p \hat \chi = \nabla \chi = \nabla \zeta \times \nabla \psi$. This
coordinate transformation changes the gyrokinetic equation to
$\partial h/\partial \phi' = 0$. Note that this transformation is
sufficient for magnetostatics, but a somewhat different ``gyro-center''
coordinate transformation is needed for electromagnetic
perturbations. This topic is addressed in H. Qin's lecture on ``A Short
Introduction to General Geometric Gyrokinetic Theory'' within this set
of Fields Institute Communications.
The gyro-average of the next-order equation leads to:
\begin{eqnarray}
&& \left[ \frac{\partial}{\partial t} + v_\parallel \hat n \cdot \nabla s
\frac{\partial}{\partial s} + \mathbf{v}_D \cdot \left( \nabla \psi
\frac{\partial}{\partial \psi} + \nabla \chi \frac{\partial}{\partial
\chi} \right) \right] h \\
&& = \frac{Ze}{T} F_M \frac{\partial}{\partial t} \langle \Phi
\rangle_{\phi'} + \frac{Ze}{m} \frac{1}{\Omega} B \left(
\frac{\partial}{\partial \psi} F_M \right) \left(
\frac{\partial}{\partial \chi} \langle \Phi \rangle_{\phi'} \right) + \langle
C(F_M, h) \rangle_{\phi'} \nonumber
\end{eqnarray}
where $s$ measures distance along the unperturbed magnetic field line
and
\begin{equation}
\langle A \rangle_{\phi'} \equiv \frac{1}{2\pi} \int_0^{2\pi} d\phi'\
A(\phi')
\end{equation}
The Poisson equation is
\begin{equation}
\nabla^2 \Phi = - 4\pi \sum_j Z_j e \int d^3v\ f_j^{(0)}
\label{eq10} \end{equation}
which reduces here to the quasi-neutrality condition
\begin{equation}
\sum_j Z_j e \int d^3v\ f_j^{(0)} = 0
\end{equation}
when $k^2 \lambda_{De}^2 \ll 1$, with $\lambda_{De}$ being the Debye
length.
The preceding equations are valid in the electrostatic limit. For the
electromagnetic generalization of this formalism, the magnetic vector
potential must be taken into account along with Ampere's Law from
Maxwell's Equations. The governing gyro-kinetic equation for
perturbations then becomes:
\begin{eqnarray}
&& \left[ \frac{\partial}{\partial t} + v_\parallel \hat n \cdot \nabla s
\frac{\partial}{\partial s} + \mathbf{v}_D \cdot \left( \nabla \psi
\frac{\partial}{\partial \psi} + \nabla \chi \frac{\partial}{\partial
\chi} \right) \right] h \\
&& = \frac{Ze}{T} F_M \frac{\partial}{\partial t} \langle \Psi
\rangle_{\phi'} + \frac{Ze}{m} \frac{1}{\Omega} B \left(
\frac{\partial}{\partial \psi} F_M \right) \left(
\frac{\partial}{\partial \chi} \langle \Psi \rangle_{\phi'} \right) + \langle
C(F_M, h) \rangle_{\phi'} \nonumber
\end{eqnarray}
where
\begin{equation}
\Psi = J_0\left( \frac{k_\perp v_\perp}{\Omega} \right) \left( \Phi -
\frac{v_\parallel}{c} A_\parallel \right) +
\frac{J_1(k_\perp v_\perp/\Omega)}{k_\perp v_\perp/\Omega}
\frac{mv_\perp^2}{Ze} \frac{\delta B_\parallel}{B}
\end{equation}
Ampere's Law can be expressed as:
\begin{equation}
\nabla_\perp^2 A_\parallel = -\frac{4\pi}{c} \sum_j \int d^2v\ Z_j e
v_\parallel J_0(k_\perp v_\perp/\Omega) h_j
\end{equation}
and
\begin{equation}
\frac{\delta B_\parallel}{B} = -\frac{4\pi}{B^2} \sum_j \int d^2v\ m_j
v_\perp^2 \frac{J_1(k_\perp v_\perp/\Omega)}{k_\perp v_\perp/\Omega}
h_j \end{equation}
where $A_\parallel$ is the perturbed parallel magnetic vector
potential and $\delta B_\parallel$ is the perturbed parallel magnetic
field, with $J_0$ and $J_1$ being the familiar Bessel functions. The
nonlinear generalizations of these linearized equations are discussed
in H. Qin's lecture on ``A Short Introduction to General Geometric
Gyrokinetic Theory'' within this set of Fields Institute
Communications.
\section{Numerical Simulations: Particle-in-Cell (PIC) Approach}
The scientific challenges related to magnetically-confined plasmas can
be categorized into four areas: macroscopic stability, wave-particle
interactions, microturbulence and transport, and plasma-material
interactions. In addition, the integrated modeling of the physical
processes from all of these areas is needed to effectively (i) harvest
the physics knowledge from existing experiments and (ii) design future
devices. Because charged particles, momentum, and heat move more
rapidly along the magnetic field than across it, magnetic fusion
research has focused on magnetic traps in which the magnetic field
lines wrap back on themselves to cover a set of nested toroidal
surfaces called magnetic flux surfaces (because each surface encloses
a constant magnetic flux). Macroscopic stability is concerned with
large-scale spontaneous deformations of magnetic flux surfaces.
These major displacements or macroinstabilities are driven by the
large electrical currents flowing in the plasma and by the plasma
pressure. Wave-particle interactions deal with how particles and
plasma waves interact. Detailed calculations of particle motions in
background electromagnetic fields are needed to assess the application
of waves to heat the plasma as well as address the dynamics of
energetic particles resulting from intense auxiliary heating and/or
alpha particles from the fusion reactions. Microturbulence and the
associated transport come from fine-scale turbulence, driven by
inhomogeneities in the plasma density and temperature, which can cause
particles, momentum, and heat to leak across the flux surfaces from
the hot interior to be lost at the plasma edge. Plasma-material
interactions determine how high-temperature plasmas and material
surfaces can co-exist. Progress in the scientific understanding in
all of these areas contributes in an integrated sense to the
interpretation and future planning of fusion systems. This demands
significant advances in physics-based modeling capabilities -- a
formidable challenge which highlights the fact that advanced
scientific codes are a realistic measure of the state of understanding
of all natural and engineered systems.
As illustrated schematically in Figure 2, the path for developing
modern high performance computational codes as validated tools for
scientific discovery involves a multi-disciplinary collaborative
process. This begins with basic theoretical research laying the
foundations for the mathematical formulation of the physical phenomena
of interest observed in experiments. Computational scientists produce
the codes which solve these equations. In order to implement the best
possible algorithms which efficiently utilize modern high-performance
computers, the optimal approach is to do so in partnership with
applied mathematicians who provide the basic mathematical algorithms
and the computer scientists who provide the requisite computer systems
software. The next step is the critical scientific code validation
phase where the newly computed results are compared against
experimental/observational data. This is a major challenge involving
a hierarchy of benchmarking criteria which begin with cross-checks
against analytic theory, empirical trends, and suggestive ``pictorial''
levels of agreement. It then graduates to sensitivity studies, where
agreement is sought when key physical parameters are simultaneously
varied in the simulation and experiment/observation. At the next
level, richer physics validation is dependent on the availability of
advanced experimental diagnostics which can produce integrated
measurements of key physical quantities such as spectra, correlation
functions, heating rates, and other variables of interest. If the
simulation/experimental data comparisons are unsatisfactory at any of
these validation levels, the work flow moves back to: [i] the
theorists (in consultation with experimentalists) if the problem looks
to be with the mathematical model; and [ii] computational scientists
(in consultation with applied mathematicians and computer scientists)
if the problem appears to be with the computational method. Even when
the theory/experiment comparisons prove satisfactory, code performance
criteria for speed and efficiency could dictate another round in the
computational science box. If all criteria are met, then the new
``tool for scientific discovery'' can be effectively utilized for
interpreting experimental data, designing new experiments, and even
predicting new phenomena of interest. This cycle of development will
of course be repeated as new discoveries with associated modeling
challenges are encountered. In addition, it should be kept in mind
that the continuous development of a robust computational
infrastructure (including hardware, software, and networking) is
needed to enable capabilities which minimize time-to-solution for the
most challenging scientific problems. Of course, even with the
successful application of the most advanced hardware and software,
there will remain complex problems that defy numerical solution.
\begin{figure}[tb]
\begin{center}
\includegraphics[width=4 in]{TangFig2.eps}
\end{center}
\caption{Development path for high performance codes as validated tools for
scientific discovery.}
\label{fig2}
\end{figure}
It is clear that since any given plasma simulation can only address a
finite range of space and time scales, the associated domains have
both minimum and maximum limits on spatial and temporal resolution.
Accordingly, simulation models are commonly developed from simplified
sets of equations, or ``reduced equations,'' which are valid for only
limited ranges of time and space scales. Examples include
``gyrokinetic equations'' \cite{Rutherford68} for dealing with
turbulent transport problems and the ``MHD equations''
\cite{Freidberg87} for addressing the large-scale stability issues.
While the reduced equations have enabled progress in the past,
fundamental restrictions on their regions of validity have motivated
the drive for improvements. In actual laboratory or natural plasmas,
phenomena occurring on different time and space scales interact and
influence one another. Simulations with greater physics fidelity thus
demand increased simulation domains, which can only result from the
derivation and application of more general equations that are valid on
a wider range of space and time scales.
The most fundamental theoretical description of a plasma comes from
kinetic equations for the distribution function within a
six-dimensional phase-space of each particle species (plus time).
They are coupled to each other through self-consistent electric and
magnetic fields. Velocity moments of these kinetic equations produce
a hierarchy of fluid equations amenable to modeling. In general, the
simulation techniques used in plasma physics fall into two broad
categories: kinetic models and fluid models. The most mature kinetic
approach is the particle-in-cell method, pioneered by John Dawson and
others \cite{Dawson83}. This method involves integrating a (possibly
reduced) kinetic equation in time by advancing marker-particles along
a representative set of characteristics within the (possibly reduced)
phase space. It basically involves a Lagrangian formulation in which
the dynamics of an ensemble of gyro-averaged particles are tracked.
Simulation techniques such as ``finite sized particles''
\cite{Langdon70} (to reduce the ``noise'' due to discrete marker
particles), ``gyro-kinetics'' \cite{Lee83} (a reduction of the full kinetic
equation to a five-dimensional phase space which removes
high-frequency motion not important to turbulent transport), and
``delta-f'' \cite{Dimits93} (a prescription for further reducing the discrete
particle noise by separating the perturbed from the equilibrium part
of the distribution function before integrating the gyrokinetic
equation along the appropriate characteristics) have been developed
over the last 20 years. These advances have served to reduce the
requirements on the number of ``particles'' necessary to faithfully
represent the physics and contributed to dramatically increasing the
accuracy and realism of the particle-in-cell simulation technique. An
alternative approach in kinetic simulations is the Vlasov or
"continuum" method \cite{Knorr62}, which involve the direct solution of the
kinetic equation governing the distribution function (examples include
the Boltzmann and Gyrokinetic equations) on a fixed Eulerian grid in
both coordinate and velocity space. Progress in the development of
associated codes in recent years has also had a significant impact on
the ability to realistically simulate microturbulent transport
phenomena.
In order to carry out particle-in-cell simulations, the
starting point is the Boltzmann equation from by Eq \ref{eq1} -- a
nonlinear partial differential equation in Lagrangian coordinates.
For particle simulations, an equivalent form is given by
\begin{equation}
\frac{dF}{dt} = \frac{\partial F}{\partial t} + \mathbf{v} \cdot
\frac{\partial F}{\partial \mathbf{x}} + \left( \mathbf{E} +
\frac{1}{c} \mathbf{v} \times \mathbf{B} \right) \cdot \frac{\partial
F}{\partial \mathbf{v}} = C(F), \label{eq15}
\end{equation}
where the distribution $F$ is specified by the Klimontovich-Dupree
representation:
\begin{equation}
F = \sum_{j=1}^N \delta(\mathbf{x} - \mathbf{x}_j) \delta(\mathbf{v} -
\mathbf{v}_j), \label{eq16}
\end{equation}
with $\mathbf{x}_j$ and $\mathbf{v}_j$ being the phase-space positions
of the $j$-th particle and $N$ being the total number of particles in
the system.
Equation \ref{eq1} can be recovered from Eq \ref{eq15} by ensemble
averaging F via the introduction of finite-size particles in the
simulations \cite{Dawson83}. In particular, when dealing with
plasmas, it is clear (as shown in Figure 3) that the Coulomb potential
for finite- sized particles is modified by Debye shielding (see
Ref. \cite{Langdon70} and references cited therein). Short range
interactions are accordingly reduced dramatically because there are
equal numbers of electrons and ions within a Debye sphere.
\begin{figure}[tb]
\begin{center}
\includegraphics[width=4 in]{TangFig3.eps}
\end{center}
\caption{Debye-shielding modification of the Coulomb potential for a
finite-sized particle.}
\label{fig3}
\end{figure}
This leads to a great simplification of the expression for the
short-range force on the $i$-th particle due to the electric-field
generated by all of the other particles, which is generally given by:
\begin{equation}
\mathbf{F}_i = q_i \mathbf{E}(\mathbf{x}_i) = \sum_{j \neq i} q_i q_j
(\mathbf{x}_i - \mathbf{x}_j)/|\mathbf{x}_i - \mathbf{x}_j|^3 \label{eq17}
\end{equation}
Now, for short ranges where $|\mathbf{x}_i - \mathbf{x}_j| \leq 1$,
the denominator, $|\mathbf{x}_i - \mathbf{x}_j|^3$, in Eq. \ref{eq17} is
simply replaced by $\lambda_D^3$; i.e.,
\begin{equation}
\mathbf{F}_i = q_i \mathbf{E}(\mathbf{x}_i) = \sum_{j \neq i} q_i q_j
(\mathbf{x}_i - \mathbf{x}_j)/\lambda_D^3
\end{equation}
The point particles here are now effectively uniformly charged spheres
of Debye- length radius. Although collisional dynamics are eliminated
by this approximation, they can be recovered as ``subgrid'' phenomena
via Monte Carlo methods with collision operators that can account for
the scattering and diffusion of particles in velocity space (see
Ref. \cite{Dawson83} and references cited therein).
Equation \ref{eq15} can be solved by tracking the temporal change of the
phase space positions of the particles. The associated basic
equations of motion for the particles are specified by the ordinary
differential equations:
\begin{equation}
\frac{d\mathbf{x}_j}{dt} = \mathbf{v}_j,
\end{equation}
and
\begin{equation}
\frac{d\mathbf{v}_j}{dt} = \frac{q}{m} \left( \mathbf{E} + \frac{1}{c}
\mathbf{v}_j \times \mathbf{B} \right)_{\mathbf{x}_j}.
\end{equation}
Since the number density of species $\alpha$ is given by
\begin{equation}
n_\alpha(\mathbf{x}) = \int F_\alpha d\mathbf{v} = \sum_{j=1}^N
\delta(\mathbf{x} - \mathbf{x}_{\alpha j})
\end{equation}
the electrostatic potential $\phi$, as noted earlier in Eq \ref{eq10}
is governed by Poisson's Equation with $\mathbf{E} = -\nabla \phi$.
\begin{equation}
\nabla^2 \phi = - 4\pi \sum_\alpha q_\alpha \sum_{j=1}^N
\delta(\mathbf{x} - \mathbf{x}_{\alpha j})
\end{equation}
Note that this linear partial differential equation is in Eulerian
coordinates (lab frame). So, unlike the ``particle pushing''
following Eqs \ref{eq16} and \ref{eq17} in the $\mathbf{x}$ and
$\mathbf{v}$ phase space, the field calculation is carried out in the
lab frame. As described in Ref. \cite{Lee03}, the electromagnetic
generalization of this formalism requires inclusion of the gyrokinetic
form of Ampere's Law.
Major progress in the simulation of the gyrophase-averaged
Vlasov-Maxwell system of equations governing low frequency
microinstabilities followed the introduction of the gyro-kinetic
methodology by W. Lee \cite{Lee83}. This involved incorporating the ion
polarization density into Poisson's Equation, and, as illustrated in
Fig. 4, the effective separation of the particle gyro-motion from its
gyro-center motion. Basically, the actual spiral motion of a charged
particle in a magnetic field is modified into that of a rotating
charged ring subject to guiding center electric and magnetic drift
motion together with parallel acceleration.
\begin{figure}[tb]
\begin{center}
\includegraphics[width=4 in]{TangFig4.eps}
\end{center}
\caption{Spiral motion of a charged particle in a magnetic field
($\mathbf{B}$) is modified by gyrokinetic approximation into that of a
rotating charged ring subject to guiding center electric and magnetic
drift motion together with parallel acceleration.}
\label{fig4}
\end{figure}
As observed from computational results \cite{Lee83} and supported by
analytic studies \cite{Krommes86}, the noise level from the
gyro-kinetic PIC simulations was found to be dramatically reduced.
One interpretation of this property is that the Debye shielding
depicted in Fig. 3 is effectively replaced by the ``gyro-radius
shielding'' introduced by the presence of the ion polarization density
in the gyro-kinetic Poisson's Equation \cite{Lee03}. Along with the
``delta-f'' prescription for further reducing the discrete particle
noise via separation of the perturbed from the equilibrium part of the
distribution function \cite{Dimits93}, modern gyro-kinetic methods
have effectively speeded up computations by 3 to 6 orders of magnitude
in time steps and 2 to 3 orders of magnitude in spatial resolution.
The accuracy and realism of the associated simulations have
accordingly benefited from such advances.
\section{Some Results from Advanced Simulations}
Even if the larger scale macroscopic disturbances in a magnetically
confined plasma could be avoided, the inherent free energy (such as
the expansion free energy in the temperature and density gradients)
can still drive turbulent cross-field losses of heat, particles, and
momentum. In fact, such increased (``anomalously large'') transport
is experimentally observed to be significantly greater than levels
expected from the collisional relaxation of toroidally-confined
plasmas (``neoclassical theory''). This is particularly important for
fusion because the effective size (and therefore cost) of an ignition
experiment will be determined largely by the balance between fusion
self-heating and turbulent transport losses. The growth and
saturation of the associated drift-type microinstabilities
\cite{Tang77} have been extensively studied over the years because
understanding this turbulent plasma transport process is not only an
important practical problem but is generally recognized as a true
scientific grand challenge. With the advent of increasingly powerful
supercomputers, it is generally agreed that this problem is
particularly well-suited to be addressed by modern terascale MPP
computational resources.
Building on the continuous progress in this area, significantly
improved models with efficient grids aligned with the magnetic field
have now been developed to address realistic 3D (toroidal) geometry
with both global and local approaches \cite{Nevins01}. As noted
earlier, solution approaches include the particle-in-cell method,
which follows the gyro-averaged orbits of an ensemble of discrete
particles in a Lagrangian formulation, and the continuum (Vlasov)
method, which directly solves the gyrokinetic equation on a fixed
Eulerian grid in both coordinate and velocity space. With regard to
the geometry of the problems addressed, the ``flux tube'' codes can
concentrate on the fine-scale dynamical processes localized to an
annular region depicted in Figure 5. The associated coordinates can
be described as being extended along equilibrium magnetic field lines,
while being localized in the perpendicular directions. Global codes
have the more imposing multi-scale challenge of capturing the physics
both on the small scale of the fluctuations (microinstabilities) and
the large scale of the equilibrium profile variations. Improved
implementation of gyrokinetic particle-in-cell algorithms as well as
gyrokinetic continuum (Vlasov) approaches have been productively
advanced \cite{Waltz02}.
\begin{figure}[tb]
\begin{center}
\includegraphics[width=4 in]{TangFig5.eps}
\end{center}
\caption{Flux-tube simulation results of turbulence localized to an
annular region of a 3D toroidal plasma.}
\label{fig5}
\end{figure}
If reliably implemented, high resolution simulations of the
fundamental equations governing turbulent transport can provide a
cost-effective means to address key phenomena that would otherwise
require expensive empirical exploration of a huge parameter space.
The progress in capturing the ion dynamics has been impressive. For
example, studies of electrostatic turbulence suppression produced by
self-generated zonal flows within the plasma show that the suppression
of turbulence caused by prominent instabilities driven by ion
temperature gradients (ITGs) is produced by a shearing action which
destroys the finger-like density contours which promote increased
thermal transport in a 3D toroidal system \cite{Lin02a}. This dynamical
process is depicted by the sequences shown in Figure 6. The lower
panels, which show the nonlinear evolution of the turbulence in the
absence of flow, can be compared against the upper panel sequence
which illustrates the turbulence decorrelation caused by the
self-generated $\mathbf{E}\times\mathbf{B}$ flow. This is also a good
example of the effective use of powerful supercomputers [e.g., the 5
teraflop IBM-SP at the National Energy Research Supercomputing Center
(NERSC) in Berkeley, CA]. Typical global particle-in-cell simulations
\cite{Lin02b,Budny00} of this type have used one billion particles
with 125 million grid-points over 7000 time-steps to produce
significant physics results. In particular, large-scale simulations
have been carried out to explore some of the key consequences of
scaling up from present-day experimental devices (around 3 meters
radius for the largest existing machines) to those of reactor
dimensions (about 6 meters). As shown in Figure 7, transport driven
by electrostatic ITG turbulence in present scale devices can change
character in larger systems. This transition from Bohm-like scaling
to Larmor-orbit-dependent ``Gyro-Bohm'' scaling is a positive trend,
because simple empirical extrapolation of the smaller system findings
would be more pessimistic. Some experimental observations in a number
of representative present-day experiments indicate that the relative
level of turbulent heat loss increases with plasma size while the size
of these eddies remains the same \cite{Petty02}. However, exceptions
to this trend, where Gyro-Bohm-like scaling sensitive to plasma
rotation was observed, have also been reported in certain
high-confinement (``H-mode'') cases \cite{Hahm04}. Such experiments
on confinement scaling properties remain a challenging area of
investigation. Nevertheless, for the larger sized reactor-scale
plasmas of the future, the present simulations would suggest that the
relative level of turbulent heat loss from electrostatic turbulence
does not increase with size. The underlying causes for why such a
transition might occur around the 400 gyroradii range indicated by the
simulations have been explored and theoretical models based on the
spreading of turbulence have been proposed
\cite{Rosenbluth98}. Although this predicted trend is a very favorable
one, the fidelity of the analysis needs to be further examined by
investigating additional physics effects which might alter the present
predictions. The analysis of associated scientific issues will
naturally demand more comprehensive physics models within
microturbulence codes. In addition to addressing experimental
validation challenges, the interplay between analytic theory and
advanced simulations will be increasingly important. For example, in
addition to the turbulence spreading theory noted, progress in physics
understanding of the nonlinear processes associated with zonal flow
dynamics has resulted both from directions provided by analytic theory
as well as by simulation results which have inspired new analytic
models \cite{Chen00,Diamond01,Malkov01,Chen01}.
\begin{figure}[tb]
\begin{center}
\includegraphics[width=4 in]{TangFig6.eps}
\end{center}
\caption{Turbulence reduction via sheared plasma flow compared to case
with flow suppressed.}
\label{fig6}
\end{figure}
\begin{figure}[tb]
\begin{center}
\includegraphics[width=4 in]{TangFig7.eps}
\end{center}
\caption{Full torus particle-in-cell gyrokinetic simulations (GTC) of
turbulent transport scaling.}
\label{fig7}
\end{figure}
An important multi-scale challenge for particle-in-cell kinetic
simulations involves dealing with the realistic implementation of
complete electron (``non-adiabatic'') physics (including important
kinetic effects, such as trapping in equilibrium magnetic wells, drift
motions, and wave-particle resonances) and electromagnetic dynamics.
These effects have largely been incorporated into gyrokinetic flux
tube (local) codes \cite{Chen01}, and present capabilities in gyrokinetic
global codes for dealing with electrostatic perturbations have been
successfully extended to include non-adiabatic electrons \cite{Lin01}. Much
more challenging for the global simulations are the electromagnetic
perturbations, which can alter the stability properties of the
electrostatic modes and also generate separate instabilities
associated with deformations of magnetic surfaces. In fact, answering
the long standing question about what causes the ubiquitously-observed
anomalously large electron thermal transport is likely linked to the
ability to deal with magnetic perturbations. They can potentially
cause a great increase in electron heat flux either through transient
deformations of the magnetic field (``magnetic flutter'') or, more
plausibly, by producing an ergodic region in which the magnetic field
lines no longer rest on nested flux surfaces but wander instead
through a finite volume ``breaking'' the flux surfaces.
In general, significant challenges for gyrokinetic simulations remain
in extending present capabilities for dealing with electrostatic
perturbations to include magnetic perturbations in cases where they
are sufficiently large to alter the actual geometric properties of the
self-consistent magnetic field. In such circumstances,
microinstabilities can drive currents parallel to the equilibrium
magnetic field, which in turn produce magnetic perturbations in the
perpendicular direction. These kinetic electromagnetic waves can
modify the stability properties of the electrostatic modes or act as
separate instabilities, such as kinetic ballooning modes
\cite{Tang80}, which can alter the magnetic topology.
In order to effectively deal with the challenging scientific issues
highlighted here, the plasma science community must address advanced
code development tasks which are important for most areas of research.
The basic goal of enhancing the physics fidelity of the codes and
developing the significantly improved software to deal with highly
complex problems involves addressing: (i) multi-scale physics such as
kinetic electromagnetic dynamics which have been discussed in earlier
sections of this article; (ii) more efficient algorithms compatible
with evolving computational architectures; and (iii) scalability of
codes necessary for utilizing terascale platforms. As the
computational hardware advances to meet the demands from the largest,
most difficult problems, it is essential to also meet the continuing
challenge of improving the scientific applications software and the
associated algorithms. Innovative approaches, such as adaptive mesh
refinement for dealing with higher dimensionality phase-space
challenges, are expected to be actively explored.
With regard to efficiently implementing present-generation codes on
the most powerful MPP computers, it is encouraging that the plasma
science community has had success in developing codes for which
computer run-time and problem size scale well with the number of
processors. A good example of this trend is illustrated in Figure 8,
where the global microturbulence PIC code, GTC, has demonstrated
excellent scalability for more than 2000 processors on the IBM-SP
computer at the National Energy Research Supercomputer Center (NERSC).
This code is the representative from the Fusion Energy Science area
within the NERSC suite of demonstration/benchmark codes to evaluate
realistic performance on new advanced computational platforms. Active
collaboration on the world's most powerful supercomputer, the Earth
Simulator Computer (ESC) in Japan, has just commenced and involved the
recent porting of this code to the ESC site. The goal of this on-
going project is to evaluate the importance of the ESC's
vector-parallel architecture compared to the much more widespread
super-scalar MPP architecture. An active collaboration has also been
initiated with a complementary global particle simulation effort in
Japan \cite{Idomura03}. Results from these early benchmark runs were quite
impressive. Specifically, utilization of 64 ESC processors yielded
results which were not only more efficient (by about a factor of two)
but were more than 20\% faster than 1024 processors on the IBM-SP
supercomputer at NERSC. Efficiency in this context refers to measured
ability of a given code to achieve the theoretically-rated performance
level of the computer processors. In addition to the ESC, the new X1
vector supercomputer at the Oak Ridge National Laboratory has
commenced bench-marking activities involving the particle-in-cell
code, GTC \cite{Oliker03}. Overall, the practical goal here is to effectively
utilize the tools, technologies, and advanced hardware systems that
will help minimize the time-to-solution for the most challenging
computational plasma physics problems.
\begin{figure}[tb]
\begin{center}
\includegraphics[width=4 in]{TangFig8a.eps}
\end{center}
\caption{3D gyrokinetic global particle-in-cell codes have
demonstrated excellent scaling as the number of processors is
increased.}
\label{fig8}
\end{figure}
It should be emphasized that a natural consequence of the effective
utilization of supercomputers is the tremendous amount of data
generated, as illustrated in Figure 9. Terabytes of data are even now
generated at remote locations (e.g., where supercomputing centers are
located), presenting data management and data grid technology
challenges \cite{Klasky03}. The data must be efficiently analyzed to
compute derived quantities. New advanced visualization techniques are
needed to help identify key features in the data. There are also
significant programming and algorithmic challenges, which must be met
in order to enable computational capabilities for addressing more
complex scientific problems. These include multi- dimensional domain
decomposition in toroidal geometry and mixed distributed/shared memory
programming. Other problems include load balancing on computers with
large numbers of processors, optimization of fundamental
gather-scatter operation in particle-in-cell codes, and scalable
parallel input/output (I/O) operations for the petascale range of data
sets.
\begin{figure}[tb]
\begin{center}
\includegraphics[width=4 in]{TangFig9.eps}
\end{center}
\caption{Terabytes of data are now generated at remote locations, as for example
the heat potential shown here on 121 million grid points from a particle-in-cell
turbulence simulation.}
\label{fig9}
\end{figure}
As noted earlier, it should be kept in mind that even with access to
greatly improved computational hardware and software advances, there
will remain limitations to what can be practically achieved
\cite{Laughlin02}. Indeed, some of the most complex plasma phenomena
involving highly transient nonlinear behavior may defy mathematical
formulation and be beyond the reach of computational physics.
\section{Future Challenges and Directions}
The ``computational grand challenge'' nature of plasma physics in
general and fusion research in particular is a consequence of the fact
that in addition to dealing with vast ranges in space and time scales
which can span over ten decades, the relevant problems involve extreme
anisotropy, the interaction between large-scale fluid-like
(macroscopic) physics and fine-scale kinetic (microscopic) physics,
and the need to account for geometric detail. Moreover, the
requirement of causality (inability to parallelize over time) makes
this problem among the most challenging in computational physics.
There has been excellent progress during the past decade in
fundamental understanding of key individual phenomena in high
temperature plasmas. Modern magnetic fusion experiments are typically
not quiescent, but exhibit macroscopic motions that can affect their
performance, and in some cases can lead to catastrophic termination of
the discharge. Major advances have been achieved in the modeling of
such dynamics, which require an integration of fluid and kinetic
physics in complex magnetic geometry. Significant progress has also
been made in addressing the dynamics governing the breaking and
reconnection of magnetic field lines, which is a central scientific
issue for fusion energy as well as allied fields such as astrophysics
and space and solar physics. Another key topic, where there have been
exciting advances in understanding, is the degradation of confinement
of energy and particles in fusion plasmas caused by turbulence
associated with small spatial-scale plasma instabilities driven by
gradients in the plasma pressure. While progress has been impressive,
the detailed physics of the growth and saturation of these
instabilities, their impact on plasma confinement, and the knowledge
of how such turbulence might be controlled remain major scientific
challenges. Finally, the understanding of overall plasma performance
requires integrating all of these issues in a simulation that includes
interactions between phenomena which were previously studied as
essentially separate disciplines. To achieve the ultimate goal of
such integration, it becomes necessary to follow the evolution of the
global profiles of plasma temperature, density, current and magnetic
field on the energy-confinement time scale with the inclusion of
relevant physics on all important time scales. While this is a
formidable long term goal, research efforts are beginning to address
such cross-disciplinary studies and to increase the physics content of
existing integrated codes. This will necessitate the development of an
architecture for bringing together the disparate physics models,
combined with the algorithms and computational infrastructure that
enables the models to work together. The associated integration of
code modules, many of which are individually at the limits of
computational resources, will clearly require substantial increases in
computer power. With the rapid advances in high-end computing power,
researchers in plasma physics a well as other fields can expect to be
able to model such systems in far greater detail and complexity,
leading eventually to the ability to couple individual models into an
integrated capability which can enable improved understanding of an
entire system.
Overall, significant advances in the development of mathematical
methodology, such as gyro-kinetic theory, together with the rapid
progress in scientific computing capabilities have enabled key
contributions to all areas of plasma science. While formidable
challenges remain, the achievements and approaches highlighted in this
article hold great promise for improving scientific understanding of
experimental data, for stimulating new theoretical ideas, and for
helping produce innovations leading to the most attractive and viable
designs for future experimental facilities. This has helped transform
traditional research approaches and has provided a natural bridge for
fruitful collaborations between scientific disciplines leading to
mutual benefit to all areas. As a natural consequence, a stimulating
path forward is provided to address the challenge of attracting,
educating, and retaining the bright young talent essential for the
future health of the field.
%\bibliographystyle{plain}
%\bibliography{fic-art}
%\nocite{*}
\begin{thebibliography}{99}
%[1] Tang W M 2002 Phys. Plasmas 9 1856
%[2] Department of Energy Office of Science's "Scientific Discovery through
%Advanced Computing (SciDAC) Program"www.science.doe.gov/scidac/, 2001;
%Dunning T (private communication)
%[3] "An Assessment of the Department of Energy's Office of Fusion
%Energy SciencesProgram," National Research Council, Fusion Science
%Assessment Committee, 2001, Final Report, Washington, D.C.: National
%Academy Press
\bibitem{Rutherford68} Rutherford, P. H. and Frieman, E. A.,
Phys. Fluids \textbf{11} (1968), 569; Taylor, J. B. and Hastie, R. J.,
Plasma Phys. \textbf{10} (1968), 479; Catto, P. J., Tang, W. M. and
Baldwin, D. E., Plasma Phys. \textbf{23} (1981), 639;
Frieman, E. A. and Chen, L., Phys. Fluids \textbf{25} (1982), 502.
\bibitem{Freidberg87} Freidberg, J. P. \textit{Ideal
Magnetohydrodynamics}, Plenum Press, New York and London, 1987; White,
R. B. \textit{Theory of Toroidally Confined Plasmas} Imperial College
Press, London, 2001.
\bibitem{Dawson83} Dawson, J. M., Rev. Modern Phys. \textbf{55} (1983), 403.
\bibitem{Langdon70} Langdon, A. B. and Birdsall, C. K., Phys. Fluids
\textbf{13} (1970), 2115; Okuda, H. and Birdsall, C. K., Phys. Fluids
\textbf{13} (1970), 2123.
\bibitem{Lee83} Lee, W. W., Phys. Fluids \textbf{26} (1983), 556; Lee,
W. W., J. Comput. Phys. \textbf{72} (1987), 243.
\bibitem{Dimits93} Dimits, A. M. and Lee, W. W.,
J. Comput. Phys. \textbf{107} (1993), 309; Parker, S. E. and Lee,
W. W., Phys. Fluids B \textbf{5} (1993), 77.
\bibitem{Knorr62} Knorr, G., Nucl. Fusion \textbf{3} (1962), 1119;
Cheng, C. Z., J. Comput. Phys. \textbf{24} (1977), 348; Denavit,
J. and Kruer, W. L., Phys. Fluids \textbf{14} (1971), 1782.
\bibitem{Lee03} Lee, W. W. and Qin, H., Phys. Plasmas \textbf{10}
(2003), 3196; Lee, W. W., Computer Physics Communications (2004), (in press).
\bibitem {Krommes86} Krommes, J. A. \textit{et al.}, Phys. Fluids
\textbf{29} (1986), 2421.
\bibitem{Tang77} Tang, W. M., Nucl. Fusion \textbf{18} (1977), 1089;
Horton, W., Rev. Mod. Phys. \textbf{71} (1999), 735
\bibitem{Nevins01} DOE Fusion SciDAC Plasma Microturbulence Project,
http://fusion/gat/com/theory/pmp, (2001); Nevins, W. M. (private
communication).
\bibitem{Waltz02} Waltz, R. E., Candy, J. and Rosenbluth, M. N.,
Proc. 19th IAEA Fusion Energy Conf., Lyon, France, (2002), paper
TH/P1-19.
\bibitem{Lin02a} Lin, Z. \textit{et al.}, Phys. Rev. Lett. \textbf{88}
(2002), 195004.
\bibitem{Lin02b} Lin, Z. \textit{et al.}, Proc. 19th IAEA Fusion
Energy Conf., Lyon, France, (2002), paper TH/1-1.
\bibitem{Budny00} Budny, R. V., \textit{et al.}, Phys. Plasmas
\textbf{7} (2000), 5038; Rhodes, T. L., APS-DPP Invited Paper UI1-5,
Bull. Am. Phys. Soc. \textbf{46} (2001), 323.
\bibitem{Petty02} Petty, C. C. \textit{et al.}, Phys. Plasmas
\textbf{9} (2002), 128.
\bibitem{Hahm04} Hahm, T. S., Diamond, P. H., Lin, Z., Itoh, K. and
Itoh, S-I., Plasma Phys. Controlled Fusion \textbf{46} (2004), A323.
\bibitem{Rosenbluth98} Rosenbluth, M. N. and Hinton, F. L.,
Phys. Rev. Lett. \textbf{80} (1998), 724.
\bibitem{Chen00} Chen, L., Lin, Z. and White, R., Phys. Plasmas
\textbf{7} (2000), 3129.
\bibitem{Diamond01} Diamond, P. H. \textit{et al.}, Nucl. Fusion
\textbf{41} (2001), 1067.
\bibitem{Malkov01} Malkov, M. A., Diamond, P. H. and Rosenbluth,
M. N., Phys. Plasmas \textbf{8} (2001), 5073.
\bibitem{Chen01} Chen, Y. and Parker, S. E., Phys. Plasmas \textbf{8}
(2001), 2095; Chen, Y. and Parker, S. E.,
J. Comput. Phys. \textbf{189} (2003), 463.
\bibitem{Lin01} Lin, Z. and Chen, L., Phys. Plasmas \textbf{8} (2001),
1447; Lee, W. W. \textit{et al.}, Phys. Plasmas \textbf{8} (2001), 4435.
\bibitem{Tang80} Tang, W. M. \textit{et al.}, Nucl. Fusion \textbf{20}
(1980), 1439.
\bibitem{Idomura03} Idomura, Y. and Tokuda, S., Nucl. Fusion \textbf{43}
(2003), 234.
\bibitem{Oliker03} Oliker, L., Canning, A., Ethier, S. \textit{et
al.}, \textit{Evaluation of Cache-based Superscalar and Cacheless
Vector Architectures for Scientific Computations}, Proc. of the
Supercomputing (SC) 2003 Conference, Phoenix, Arizona, U. S.,
Nov. 2003, http://www.sc-conference.org/sc2003/paperpdfs/pap255.pdf;
Ethier, S. (private communication).
\bibitem{Klasky03} Klasky, S., Ethier, S., Lin, Z., \textit{et al.},
\textit{Grid-Based Parallel Data Streaming implemented for the
Gyrokinetic Toroidal Code}, Proc. of the Supercomputing (SC) 2003
Conference, Phoenix, Arizona, U. S., Nov. 2003,
http://www.sc-conference.org/sc2003/paperpdfs/pap207.pdf.
\bibitem{Laughlin02} See for example, Laughlin, R. B.
\textit{Physical Basis of Computability}, Computing in Science and
Engineering \textbf{4} (2002), 27.
\end{thebibliography}
\end{document}