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.