High Performance Computing for Beam Physics Applications

                   Robert D. Ryne and Salman Habib

      Los Alamos National Laboratory, Los Alamos, NM 87544, USA

INTRODUCTION

Several countries are now involved in efforts aimed at utilizing
accelerator-driven technologies to solve problems of national and
international importance. These technologies have both economic and
environmental implications. The technologies include waste transmutation,
plutonium conversion, neutron production for materials science and
biological science research, neutron production for fusion materials
testing, fission energy production systems, and tritium production.
All of these projects require a high-intensity linear accelerator that
operates with extremely low beam loss. This presents a formidable
computational challenge: One must design and optimize over a kilometer of
complex accelerating structures while taking into account beam loss to
an accuracy of 10 parts per billion per meter. Such modeling is essential
if one is to have confidence that the accelerator will meet its beam loss
requirement, which ultimately affects system reliability, safety and cost.
At Los Alamos, we are developing a capability to model ultra-low loss
accelerators using the CM-5 at the Advanced Computing Laboratory. We are
developing PIC, Vlasov/Poisson, and Langevin/Fokker-Planck codes for this
purpose. With slight modification, we have also applied our codes to
modeling mesoscopic systems and astrophysical systems. In this paper,
we will first describe HPC activities in the accelerator community.
Then we will discuss the tools we have developed to model classical
and quantum evolution equations. Lastly we will describe how these
tools have been used to study beam halo in high current,
mismatched charged particle beams.

OVERVIEW OF HPC ACTIVITY IN THE ACCELERATOR COMMUNITY

To the best of our knowledge, the first use of parallel computing
in the accelerator community was performed by Paul Schoessow of
Argonne National Laboratory in 1989, who used a Thinking Machines CM-2
and an Alliant FX/8 to perform wake field calculations [1].
In the early 1990s, researchers at the Superconducting Supercollider
Laboratory (SSCL) used an iPSC/860 parallel computer to perform particle
orbit tracking [2]. Researchers at Lawrence Berkeley Laboratory have
done magnet calculations (based on the Biot-Savart law) using
a MasPar, a CM-5, and a cluster of Sun Sparc 10s with PVM [3]; this
work was originally done to study field quality of superconducting
magnets for the SSCL, but it is now being applied to model quadrupole
magnets of the interaction region of the SLAC B-factory. An example of
coarse-grained parallel computing is the use of PVM at Fermilab to
perform photoinjector design studies with the code Parmela [4].
A few researchers in the accelerator community are involved with
object oriented techniques; in particular, Leo Michelotti of Fermilab
has developed C++ classes for the design, analysis, and modeling of
accelerators and detectors [5]. The remainder of this paper will
describe an effort at Los Alamos National Laboratory to develop parallel
codes for modeling classical and quantum evolution equations, including
codes to model beam halo in next-generation linear accelerators. 

MASSIVELY PARALLEL SIMULATION OF CLASSICAL AND QUANTUM EVOLUTION EQUATIONS

(Note: This section has been modified for the ascii version of this paper.
The original discussion can be found in the TeX or PostScript version).

In order to model the evolution of intense charged particle beams
we utilize both particle simulations and direct Vlasov/Poisson solvers.
To solve the Vlasov/Poisson equation on the CM-5, we utilize a spectral
method coupled with split-operator symplectic integration algorithms [6].
In the context of symplectic integration algorithms, Yoshida showed
how to take an algorithm of order 2n and construct from it an
algorithm of order 2n+2 [7]. Using this approach, we have developed
second order and fourth order Vlasov/Poisson codes on the CM-5.
Besides Vlasov/Poisson codes, we have also developed particle simulation
codes for modeling intense beams. Our two-dimensional code uses area
weighting to deposit charge on the numerical grid, a convolution
technique to compute the space charge fields, and the symplectic
integration algorithms described above to advance particles.

Though this paper emphasizes modeling beam halo, it is worth
pointing out that this approach can be applied to several classical and
quantum systems. For example, the Vlasov code can be modified to study
gravitating systems. An analogous code can be written to model the
Schrodinger equation. In fact, a second order implementation of this
code turns out to be identical to the approach used by Feit and Fleck
to study a variety of quantum systems, including the vibrational energy
levels of triatomic molecules [8].  Using Yoshida's technique, it is in
principle possible to obtain high order Schrodinger codes. We have
developed fourth order codes on the CM-5 for modeling the Schrodinger
equation. Similar codes for evolving the density matrix and the Wigner
distribution function, both in the absence and presence of dissipation
and noise, are under development.

MODELING BEAM HALO IN ULTRA LOW LOSS ACCELERATORS

In recent years there has been a steady growth in activity worldwide
aimed at utilizing accelerator-driven schemes to solve a variety of
important problems. The technologies that are being pursued are:
(1) accelerator transmutation of waste (ATW), accelerator-based
conversion of plutonium (ABC) resulting from nuclear disarmament;
(3) accelerator-based production of spallation neutrons for materials
science and biological science research; (4) accelerator production of
14-MeV neutrons for use at fusion materials irradiation facilities (FMIF);
(5) accelerator-based energy production; and (6) accelerator production
of tritium. The fact that so many governments are pursuing these areas is
a testament to the importance of these emerging technologies. For example,
Japan is involved in ATW and FMIF; Russia is involved in ATW and ABC;
France is involved in ATW; in Europe, Carlo Rubbia has proposed a scheme
for accelerator-based fission energy systems (a scheme originally proposed
by Los Alamos physicist Charles Bowman); and the United States is involved
in efforts related to all the above technologies. The nation that succeeds
in developing these technologies will have the ability to tackle difficult
environmental problems while at the same time enhancing its international
economic competitiveness.

All of the above projects require a high intensity accelerator that
operates with extremely low beam loss (approximately ten parts per
billion per meter). It is now known that a major source of beam loss
is the formation of a very low density beam halo at a radius several
times the rms radius. Understanding the nature of this halo and finding
ways to minimize it is crucial to the future of the above mentioned
technologies. Ultimately, the goal of beam halo studies is to
reduce risk, reduce cost, improve reliability and improve efficiency
of high current linacs associated with these projects. 

Recently there has been much activity related to the so-called
core-halo model of beam halo evolution [9]. In this model the core of
the beam is assumed to oscillate because of an initial radial mismatch;
the fields inside the core are roughly proportional to r, while they
are inversely proportional to r outside the core. Thus, a halo particle
moving in and out of the core sees a strongly time-dependent nonlinear
field superposed on a linear external focusing field. A useful way to
analyze the resulting dynamics is to make a stroboscopic plot, in which
test particles are plotted once each cycle of the core oscillation.
Figure 1 shows a stroboscopic plot of 32 test particles in the core-halo
model [see halo515.eps]; points are plotted in (x,px)-space each time
the core radius reaches a minimum. The main features are: (1) a central
region that has an extent of slightly more than the core radius;
(2) a large amplitude region; (3) a period-2 resonant structure
(associated with the fixed points to the left and right of the central
region); and (4) a separatrix with an inner branch that encloses the
central region and an outer branch that separates the period-2 resonance
from the large amplitude region.  In this figure the separatrix is
actually a narrow chaotic band, and the outer edge of the chaotic band
has the approximate shape of a peanut.  This "peanut diagram" provides
a useful picture for describing halo formation: particles in an initially
well-defined core may diffuse toward the separatrix by a number of
mechanisms such as nonlinear instabilities and collisions; if they cross
the separatrix then they will be carried along its outer branch to large
amplitude. Features of the core-halo model agree well with detailed
particle simulations run on the CM-5 at the Advanced Computing Laboratory.
Figure 2 shows simulation results for an initially mismatched
Kapchinskij-Vladimirskij (KV) distribution in a constant focusing channel
[see peanut.eps]. Though the initial distribution has the property that
it is uniform in (x,px)-space, it is unstable for the parameters chosen,
and the resulting phase space of Figure 2 is highly nonuniform. The
simulation was run with 2 million particles. The first curve bounding
the chaotic band in Figure 1 is also shown in Figure 2 for comparison.
The CM-5 results show that, for this simple configuration (an axially
symmetric beam in a constant focusing channel), the maximum particle
amplitudes are in good agreement with the amplitude of the separatrix
in the core-halo model.

Particle simulations provide a means of modeling the Vlasov equation by
following the trajectories of a large number of interacting particles.
(Since the Vlasov equation does not include the effects of two-particle
collisions, the above statement assumes that collisions between
particles in the simulation have a negligible effect). In the direct
approach described earlier, the distribution function is specified on a
grid in phase space, and the values of the function on the grid are
stepped forward in time. At present we have developed a 2-dimensional
(x,px,y,py)-code that runs on the CM-5 [10]. An example output is shown
in Figure 3, which corresponds to an initially mismatched 4-dimensional
Gaussian beam in a quadrupole channel [see vlas235.eps]. The figure,
which shows the beam density, was obtained by integrating the 4-dimensional
distribution over the variables px and py. This simulation utilized
a (128 x 128 x 128 x 128) grid, for a total of 268 million grid points;
a simulation of 300 time steps on 512 nodes of the CM-5 required 6 GBytes
of memory and 2.5 hours of cpu time. At present, we are continuing to
develop both particle simulations and Vlasov/Poisson simulations so
that we can make comparisons of these two approaches.

CONCLUSIONS AND FUTURE WORK

We have successfully developed one-dimensional and two-dimensional
Vlasov/Poisson codes and particle simulation codes on the CM-5
for modeling intense charged particle beams. Our codes have helped
validate the core-halo model of beam halo propagation, and they will
be valuable for studying nonidealized systems that are not amenable
to purely analytic treatment. Also, the ability to perform simulations
with several million particles has enabled us to study fine structure
in beams with unprecedented resolution. Using nearly identical
techniques, we have developed Schrodinger equation codes. In the future,
we plan to use our quantum modeling capability to study mesoscopic
systems, including the interaction of atoms with intense laser fields.
We also intend to develop Langevin and Fokker-Plank codes to study
the effects of collisions in charged particle beams and astrophysical
systems. In the area of modeling beam halo in ultra-low loss
accelerators, we intend to begin working on three-dimensional codes
and the inclusion of nonideal effects such as beam misalignments
and magnet fringe fields. Through a combination of numerical and
analytical approaches, we are increasing our understanding of beam
halo so that we can systematically design, evaluate, and optimize
linacs for accelerator driven technologies.

ACKNOWLEDGEMENT

This research was performed in part using the resources located at the
Advanced Computing Laboratory of Los Alamos National Laboratory,
Los Alamos, NM 87545

REFERENCES
[1] P. Schoessow, Proceedings of the Conference on Computer Codes and
    the Linear Accelerator Community, Los Alamos, p. 377 (January 1990).
[2] B. Cole et al., Proceedings of the 1991 Particle Accelerator
    Conference, San Fransisco, CA, p. 204 (May 1991); G. Bourianoff et al.,
    Proceedings of the 1993 PAC, Washington, DC p. 128 (May 1993).
[3] S. Caspi and R. Hinkons, private communication.
[4] Eric Colby, private communication.
[5] L. Michelotti, Stability of Particle Motion in Storage Rings,
    AIP Conference Proceedings No. 292, Upton, NY, p. 488 (1992).
[6] E. Forest and R. Ruth, Physica D 43, p. 105 (1990).
[7] H. Yoshida, Phys. Lett. A 150, p. 262 (1990);
    E. Forest et al., Phys. Lett. A 158, p. 99 (1991).
[8] M. Feit and J. Fleck, Jr., J. Chem. Phys. 78(1), p. 301 (1983);
    M. Feit et al., J. Comp. Phys. 47, p. 412 (1982);
    M. Hermann and J. Fleck, Jr., Phys. Rev. A, 38(12), p. 38 (1988).
[9] J.O'Connell et al., Proceedings of the 1993 Particle Accelerator
    Conference, Washington, DC p. 3657 (May 1993);
    J. Lagniel, Nucl. Inst. Meth. A345, No. 1, p. 46 (1994);
    J. Lagniel, Nucl. Inst. Meth. A345, No. 3, p. 405 (1994);
    T. Wangler, Los Alamos preprint, LA-UR-94-1135 (1994);
    R. Gluckstern, to appear in Phys. Rev. Lett.
[10]S. Habib and R. Ryne, in preparation.

FIGURE CAPTIONS
Caption for Figure 1 (halo515.eps):
    Stroboscopic phase space plot based on the core-halo model with
    a tune depression of 0.5 and a mismatch of 1.5 (i.e. the beam
    oscillates with a maximum amplitude 1.5 times the matched value).
    The location of 32 test particles is plotted every time the beam
    envelope reaches a minimum, for a total of 1000 oscillations.

Caption for Figure 2 (peanut.eps):
    Beam phase space from a 2 million particle simulation on
    the CM-5 (65536 points are plotted). The outer peanut-shaped
    set of points were obtained from the core-halo model.
    The CM-5 results show that, for this configuration (an axially
    symmetric beam in a constant focusing channel), the core-halo
    model provides a good estimate of the maximum particle amplitudes.

Caption for Figure 3 (vlas235.eps):
    Output from a direct Vlasov/Poisson simulation. The 4-dimensional
    distribution function was integrated over px and py to obtain the
    beam density on a 128x128 grid.