# Simulating Intense Ion Beams for Inertial Fusion Energy\*

## Alex Friedman, LLNL and LBNL

## Science area description

The Heavy Ion Fusion (HIF) program's goal is the development of the body
of knowledge needed for Inertial Fusion Energy (IFE) to realize its
promise. The intense ion beams that will drive HIF targets are
*nonneutral plasmas* and exhibit collective, nonlinear dynamics which
must be understood using the kinetic models of plasma physics. This beam
physics is both rich and subtle: a wide range in spatial and temporal
scales is involved (see figure 1), and effects associated with both
instabilities and non-ideal processes must be understood. Ion beams have
a "long memory," and initialization of a beam at mid-system with an
idealized particle distribution introduces uncertainties; thus, it will
be crucial to develop, and to extensively use, an integrated and
detailed "source-to-target" HIF beam simulation capability. We begin
with an overview of major issues.

![](media/image1.pct){width="3.8895833333333334in"
height="2.134027777777778in"}

Figure 1. Timescales in driver and fusion chamber. Spatial scale lengths
range from electron gyroradii in magnets \~ 0.01 mm, to beam Debye
length \~ 1 mm, to beam radius \~ 1 cm, to machine length \~ km\'s.

[Long-term evolution of space-charge-dominated beams]{.underline}: In
the driver, the array of beams is accelerated by inductive electric
fields, and is confined by applied "focusing" fields. The beams dynamics
are "space charge dominated," that is, governed by a balance between the
applied fields and space charge forces. This contrasts with the more
usual situation in high-energy accelerators, where thermal pressure
typically dominates over space charge forces. The beam dynamics is
collisionless and Liouvillean, that is, the phase space density remains
constant along particle orbits. As a result, "emittance growth"
(dilution of the phase space) takes place through complicated
distortions driven by collective processes, imperfect applied fields,
image fields from nearby conductors and inter-beam forces. Such dilution
must be kept minimal, because of the necessity to focus the beams
ultimately onto a small (few mm) focal spot on the fusion target.
Simulations must capture the effects of small influences which act over
long distances. In addition, collective beam modes, and interactions
with the external environment which can drive resistive-wall
instabilities, must be understood and minimized. Other challenges
include the need to accurately simulate time-dependent
space-charge-limited emission from curved surfaces, a singular problem.
This area is computationally challenging because of the need for an
efficient but detailed description of the applied fields, and the needs
for good statistics and mesh resolution: 10^7^-10^8^ particles, \~10^5^
steps, more than 100 cells in each transverse direction, and thousands
of cells in the longitudinal direction, for a simulation with a
co-moving mesh that "treadmills" with the beam.

[Beam halo generation]{.underline}: In modern particle accelerators, the
confining fields are dominated by an "alternating-gradient" transverse
quadrupole moment from a sequence of electric or magnetic lenses. Thus,
the confining fields are non-steady in the beam frame, complicating
analysis; no exact nonsingular equilibria are known, and it may be that
none exist. Oscillations of the beam "core" can parametrically pump
particles into an outlying, or "halo," population. For focusability and
also to avoid the adverse effects of ions impinging on walls, beam halo
must be kept minimal. Here, PIC methods have been used, but emerging
continuum-Vlasov and nonlinear-perturbative methods may offer
advantages.

[Multispecies effects in driver]{.underline}: collective beam
interactions with "stray" electrons in the accelerator and transport
lines, including electron generation and trapping within the beam, must
be understood quantitatively. This area is computationally challenging
because of the ratio between the fast time scale for electron motion and
the slow time scale for electron build-up within the beam; the need to
efficiently gather/scatter and communicate multi-species information for
ionization and surface-physics processes; and the needs for efficient
dynamic load balancing and perhaps an adaptive mesh. Electron scales in
the driver will not be resolved in full end-to-end simulations in the
near future, so coupled disparate-timescale simulations capturing
electron dynamics and/or "subscale" modeling will be employed.

[Beam interactions with fusion chamber environment]{.underline}: 3-D
simulations of the propagation of the cluster of beams through the final
focusing optics, and onward through the fusion chamber's environment of
gas and plasma, are required in order to provide a realistically
complete model of the target illumination. The beam and background
plasma dynamics include: multibeam effects; return current formation and
dynamics (streaming instabilities); imperfect neutralization; beam
stripping; emittance growth; and photo-ionization of the beam ions and
background gas. Of these effects, many of the greatest uncertainties and
computational challenges are associated with multiple-beam interactions
near the target, and these will be one important focus of research
efforts. Another important focus of research efforts in the fusion
chamber will be on collective instabilities, such as resistive hose,
filamentation and two-stream modes.

In the chamber as in the driver, it is appropriate to employ multiple
methods, and we will use PIC, hybrid PIC-fluid, and nonlinear
perturbative ("delta-f") methods. The chamber calculations must allow
exploration of various propagation modes, e.g. "neutralized-ballistic,"
"assisted-pinch," etc. The challenges include the need for complex
physics models, "outgoing-wave boundary conditions," an implicit hybrid
model for the dense-plasma scenarios, and of order 10^7^-10^8^
simulation particles. It will also be valuable to employ multiple
models, so as to compare, *e.g.*, implicit electromagnetic (EM) methods
(which can stably under-resolve fast time scales not essential to the
physics) with explicit EM methods and with magneto-inductive ("Darwin")
methods that eliminate light waves from the description.

## Key computational tools

The HIF program has developed tools to explore these physics areas,
primarily (but not exclusively) using particle-in-cell (PIC) methods.
Plans involve adaptation of exiting codes to run optimally on computers
that use a hybrid of shared and distributed memory, tight coupling
between those tools, production of new and improved numerical
algorithms, e.g., averaging techniques that allow larger time-steps, and
development of improved physics models. Some of this work already relies
heavily on modern scripting techniques for code steering, and advanced
data visualization is playing an increasing role. In all areas,
benchmarking with theory, with experiments, and among codes will
continue to be essential.

The codes to be improved, coupled, and employed are listed in Table 1.
All are well-positioned to move quickly to the new hardware platform.
WARP and LSP are fairly large and complex codes offering many options;
they are, more accurately, code frameworks. BEST and BPIC are smaller
codes that are attractive test-beds for new methods, in addition to
being useful tools in their own right. The majority of the required
simulations are between one and two orders-of-magnitude beyond current
practice; we anticipate typical run times of order a day. Here we
characterize the codes by their methods, and by their regimes of
applicability.

> Follow particles (plasma particle-in-cell method)
>
> WARP (driver): 3-D (or r,z or x,y) ES, detailed lattice
>
> LSP (chamber & driver): 3-D or (r,z) implicit (or explicit) EM or ES,
> hybrid (kinetic/fluid)
>
> BPIC (chamber): 3-D or (r,z) EM, moving grid, outgoing waves
>
> Follow particles and perturbation to distribution function (δƒ)
>
> BEST (chamber & driver): 3-D EM, Darwin, or ES, offers reduced noise
>
> Evolve distribution function (ƒ) on a grid
>
> WARP-SLV (driver): 2-D (x,y,p~x~,p~y~) ES Semi-Lagrangian Vlasov
> solver
>
> Evolve moments of distribution function
>
> WARP-CIRCE & WARP-HERMES (driver): transverse moments, longitudinal
> Lagrangian fluid

Table 1. Classes of codes for HIF beams, codes in use or being
developed, and domains of applicability

WARP offers 3D and transverse-slice 2-1/2D geometries, and is used
extensively throughout the Heavy Ion Fusion program for studies of beams
in the accelerator, pulse-compression line, and final focusing system.
WARP runs in parallel on NERSC's Cray T3E-900 and IBM SP, using message
passing. Good scaling has been obtained using up to 256 processors on
problems of intermediate size. The code is written in a coarse-grained
object-oriented superset of Fortran, and runs under the control of the
Python scripting language. A prototype continuum Vlasov package, SLV,
was implemented in the WARP code and is currently running on simple
axisymmetric beam physics problems. Moment based models CIRCE and
HERMES, useful for rapid scoping and synthesis, are also implemented
within the WARP framework.

LSP offers (*r*,*z*) and 3D geometries, implicit or explicit EM or ES
PIC and fluid models, a multi-block mesh which allows simulation of
non-rectangular (*e.g.*, L-shaped) regions, and domain decomposition
designed for multilevel memory access. Its implicit hybrid model enables
simulation of dense-plasma scenarios. LSP has extensive gas and surface
interaction physics models; it already offers secondary emission,
kinetic neutrals, ionization, scatter and neutral recycling, and has
achieved good scaling using up to 256 processors on problems of
intermediate size. LSP is written in C using an object-oriented style.

The BPIC chamber-propagation code offers an explicit 3-D or (r,z)
electromagnetic PIC model, and uses a time-evolving mesh and uses novel
methods for "advecting away" the field errors associated with the
evolving cell boundaries. It is written in Fortran 95.

BEST offers nonlinear-perturbative ("delta-f") simulation in 3D geometry
and has been parallelized using a combination of MPI and OpenMP and
two-dimensional domain decomposition suitable for a supercomputer
equipped with both shared and distributed memory. The code was designed
to elucidate mode structures by minimizing discrete-particle noise,
employs a new Darwin (magnetoinductive) model algorithm, and, to
compensate for the mass ratio (about 250,000) of the heavy ions to the
electrons, uses a newly-developed adiabatic pushing and deposition
algorithm. Good scaling has been obtained using up to 512 processors. It
is written in Fortran 95.

## Representative simulations

Existing codes are used for a wide variety of simulations, as
illustrated in Figure 2. We do not attempt to describe the relevant
physics in any detail, but merely show the range of ongoing activities.

![](media/image2.pct){width="6.498611111111111in"
height="4.091666666666667in"}

Figure 2. Representative output from HIF beam simulations: (a) WARP3d
simulation of space-charge-limited emission off a curved surface, and
acceleration in a 3-D structure, including subgrid-scale placement of
conductor boundaries (cut-cell method); (b) WARPxy study of beam
emittance versus time in an imperfectly-aligned beamline, for five
different intervals between applications of steering.; (c) WARP3d study
of longitudinal waves on beam, driven unstable by impedance of
accelerating structures; (d) accelerating waveforms for a possible
future experimental accelerator, for use in WARP3d simulations; (e) BEST
simulation of unstable electron-ion two-stream mode in a beamline; (f)
semi-Lagrangian Vlasov simulation of beam halo generation due to
anharmonic focusing fields, using prototype model in WARP-SLV; (g)
distorted beam phase space in final focusing, as simulated using WARPxy;
(h) BPIC simulation of beam transit through fusion chamber environment
and onto the target.

## Approach

A concept for source-to-target simulation is shown in Figure 3. In this
scenario, the beam is simulated from the source through the final
focusing optic using WARP3d, and the particle and field data are then
transferred into LSP (or BPIC) where the simulation is carried through
to the fusion target. At that point the particle data is used to
generate "ray" information for the ion beam source in the target
simulation code. Meanwhile, LSP is used to study electron effects in the
driver, especially sources and trapping in the beam; for this to be
accurate, it is necessary to understand beam halo quantitatively, and
for this the marker-following capabilities of BEST and/or the Vlasov
solver in WARP are employed in coupled side calculations. BEST is also
used to study beam instabilities in detail, using parameters transferred
from WARP and LSP.

![](media/image3.pct){width="4.2in" height="2.6666666666666665in"}

Figure 3. Depiction of source-to-target simulation strategy (see text)

The major developments required include: (i) optimization of codes for
efficiency on the new computer architecture; (ii) development of new and
improved numerical algorithms (*e.g.,* for larger timesteps and Vlasov
solution); and (iii) development of improved physics models (*e.g*., for
multibeam, converging beam, self-magnetic, atomic physics, and module
impedance effects) that will be made practical by the terascale
capability. The codes will be linked using scripting tools (especially
Python) for intercommunication and code steering, "workspace" tools for
heterogeneous computations, and self-describing data files (*e.g.,*
NetCDF). The "data glut" associated with saving information from the
many processors will be addressed by incorporating optimized parallel
I/O capabilities. The challenges of visualizing a time-dependent 6D
phase space will be addressed through the use of volume and isosurface
rendering, coupled with projection and range selection along the
non-visualized coordinate directions; animation will also be further
developed and employed. The simulations will entail self-consistent
field descriptions requiring interprocessor communication, but will
employ optimized domain decomposition and dynamic load balancing so as
to be scalable on terascale architectures.

## Relationships with other scientific disciplines

This research area involves nonlinear dynamics, self-consistent fields,
large-scale parallel computations, massive data handling, interactive
and script-driven code steering, and visualization of a time-dependent
multidimensional phase space. These aspects appear in many emerging
applications of terascale computing, and considerable cross
fertilization with other areas can be anticipated.

Other accelerator applications are moving toward higher beam
intensities, and the knowledge gained via this research into very strong
space-charge regimes will be relevant to a wide variety of applications.
We anticipate long-term benefits to such efforts as the Spallation
Neutron Source, the Very Large Hadron Collider, the Next Linear
Collider, the Accelerator for Transmutation of Waste, the Muon Collider,
and for application to Boron Neutron Capture Therapy.

HIF researchers are collaborating with the NERSC computational science
group in the integration of Adaptive Mesh Refinement (AMR) techniques
with the Heavy Ion Fusion PIC simulation code WARP3d. That group
initially developed the AMR method for application to combustion and
fluid flow studies. We anticipate that the method will be useful in
simulation studies of heavy ion beams in several contexts: mesh
refinement around the beam in a PIC code; around internal conducting
structures to capture subtle but important field details; and around key
phase-space structures in a continuum Vlasov calculation in 4D, 5D, and
ultimately 6D, where straightforward methods would require a very large
mesh.

[Acknowledgements:]{.underline} This document incorporates input from
others in the Heavy Ion Fusion Virtual National Laboratory (a
partnership between LBNL, LLNL, and PPPL) and its collaborators:
R. C. Davidson, PPPL; J. J. Barnard and R. H. Cohen, LLNL; J-L. Vay,
LBNL; and D. R. Welch, Mission Research Corporation.

\*This work was performed under the auspices of the U.S. Department of
Energy by the University of California, Lawrence Livermore National
Laboratory under Contract Number W-7405-Eng-48.