MODELING OF ACCELERATOR SYSTEMS AND EXPERIMENTAL VERIFICATION
OF QUARTER-WAVE RESONATOR STEERING
By
Carla Benatti

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
Physics - Doctor of Philosophy
2014

ABSTRACT
MODELING OF ACCELERATOR SYSTEMS AND EXPERIMENTAL
VERIFICATION OF QUARTER-WAVE RESONATOR STEERING
By
Carla Benatti
Increasingly complicated accelerator systems depend more and more on computing power
and computer simulations for their operation as progress in the field has led to cutting-edge
advances that require finer control and better understanding to achieve optimal performance. Greater ambitions coupled with the technical complexity of today’s state-of-the-art
accelerators necessitate corresponding advances in available accelerator modeling resources.
Modeling is a critical component of any field of physics, accelerator physics being no
exception. It is extremely important to not only understand the basic underlying physics
principles but to implement this understanding through the development of relevant modeling
tools that provide the ability to investigate and study various complex effects. Moreover,
these tools can lead to new insight and applications that facilitate control room operations
and enable advances in the field that would not otherwise be possible.
The ability to accurately model accelerator systems aids in the successful operation of
machines designed specifically to deliver beams to experiments across a wide variety of
fields, ranging from material science research to nuclear astrophysics. One such accelerator
discussed throughout this work is the ReA facility at the National Superconducting Cyclotron Laboratory (NSCL) which re-accelerates rare isotope beams for nuclear astrophysics
experiments.
A major component of the ReA facility, as well as the future Facility for Rare Isotope
Beams (FRIB) among other accelerators, is the Quarter-Wave Resonator (QWR), a coaxial

accelerating cavity convenient for efficient acceleration of low-velocity particles. This device is
very important to model accurately as it operates in the critical low-velocity region where the
beam’s acceleration gains are proportionally larger than they are through the later stages of
acceleration. Compounding this matter, QWRs defocus the beam, and are also asymmetric
with respect to the beam pipe, which has the potential to induce steering on the beam. These
additional complications make this a significant device to study in order to optimize the
accelerator’s overall performance.
The NSCL and ReA, along with FRIB, are first introduced to provide background and
motivate the central modeling objectives presented throughout this work. In the next chapter,
underlying beam physics principles are then discussed, as they form the basis from which
modeling methods are derived. The modeling methods presented include multi-particle
tracking and beam envelope matrix transport. The following chapter investigates modeling
elements in more detail, including quadrupoles, solenoids, and coaxial accelerating cavities.
Assemblies of accelerator elements, or lattices, have been modeled as well, and a method for
modeling multiple charge state transport using linear matrix methods is also given.
Finally, an experiment studying beam steering induced by QWR resonators is presented, the
first systematic experimental investigation of this effect. As mentioned earlier, characterization
of this steering on beam properties is important for accurate modeling of the beam transport
through the linac. The measurement technique devised at ReA investigates the effect’s
dependence on the beam’s vertical offset within the cavity, the cavity amplitude, and the
beam energy upon entrance into the cavity. The results from this experiment agree well
with the analytical predictions based on geometrical parameters calculated from on-axis field
profiles. The incorporation of this effect into modeling codes has the potential to speed up
complex accelerator operations and tuning procedures in systems using QWRs.

ACKNOWLEDGMENTS

To begin with, I could not have accomplished any of this without the guidance of my advisor,
Mike Syphers. His instruction and patience were truly essential factors in making me the
scientist that I am today. If there is one person I wish to emulate in my career, it is Mike. I
can only hope to one day be half the mentor to others that he has been for me.
I also received indispensable advice and guidance from Paul Chu, for which I am deeply
grateful. I cannot thank him enough for taking me on as a student who was in need of
much direction and guidance, which he provided. He has instilled in me many interests in
computing and programming which I plan to pursue throughout my career.
The remaining members of my guidance committee, Sekhar Chivukula, Felix Marti, and
Stuart Tessmer, also deserve my thanks for their advice and support throughout the years.
In addition, the kind support from Walter Wittmer, Daniela Leitner, Alberto Rodriguez, and
Matt Hodek has been amazing. They have all helped me in their own ways throughout my
time as a student at MSU, enabling me to grow and develop in a supportive environment.
I’m also extremely grateful for Marc Doleans and Oliver Kester, who started me off on my
adventure in accelerator physics.
The measurements I took at ReA would not have been possible without the help from its
capable operators: Shannon Krause, Randy Rencsok, Sam Nash, Ken Galbraith, Eric Mollon,
and Jay Paquette. Larry Tobos kept the LB source behaving throughout the measurements,
for which I’m sincerely grateful. I also thank Alberto Facco for his consultation during the
calculations preceding these measurements. Xiaoyu Wu, Qiang Zhao, and Mauricio Portillo
helped me tremendously with modeling, especially with IMPACT and COSY, and for that I
am very thankful.
iv

I am grateful for the assistance from Jack Ottarson, Ben Arend, Don Lawton, and Erin
O’Donnell, which led to many of the figures and diagrams in this thesis. I also have the
physics graduate secretary, Debbie Barratt, to thank for all of her help on the administrative
side. I also have Eduard Pozdeyev and Winter Zheng to thank for their contributions with
the solenoid parameterization and cavity field profile data respectively.
My fellow accelerator physics graduate students, especially Dan Alt, Jeremiah Holzbauer,
Kritsada Kittimanapun, Ling-Ying Lingy, and Daniel Winklehner, along with Post-Doc Bojan
Durickovic, have also helped me tremendously throughout my time at MSU. I would like to
extend my gratitude to them for all they’ve done in helping me throughout my graduate
career at MSU.
I’d also like to sincerely thank my family and friends for their continuous support. My
parents, Liz and Henry Benatti, especially have always supported me in every endeavor,
and my brother Dominic keeps me grounded by continuing to beat me in video games.
My soon-to-be-family Cheryl Barquest and Jack Robinson, and Jim Barquest and Linda
Slaughter, have also supported me every step of the way. My friends Mandie Gehring, Kate
Mikels, Emma Morrison, Alexis Prasov, and Kathy Walsh (among many others) all know the
struggles I’ve been through and have been there to help me throughout them all.
Finally, above all, I’d like to acknowledge my Bradley, who believed in me when I didn’t
believe in myself. He is the main reason I survived graduate school, and I can never thank
him enough for everything he has done for me. Thank you for inspiring me, thank you for
believing in me, thank you for supporting me. I am forever grateful.

v

TABLE OF CONTENTS

LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii
LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

ix

Chapter 1 Introduction . . . . . . . . . . .
1.1 NSCL . . . . . . . . . . . . . . . . . . .
1.1.1 ECR Ion Source . . . . . . . . . .
1.1.2 Cyclotrons . . . . . . . . . . . . .
1.1.3 RIB Production . . . . . . . . . .
1.1.4 Fragment Separator . . . . . . . .
1.1.5 Experimental Stations . . . . . .
1.1.6 ReAcceleration Motivation . . . .
1.1.7 ReA . . . . . . . . . . . . . . . .
1.2 Future Facilities . . . . . . . . . . . . . .
1.2.1 ReA Upgrades . . . . . . . . . . .
1.2.2 FRIB . . . . . . . . . . . . . . . .
1.3 Modeling Considerations . . . . . . . . .
1.3.1 Modeling Methods . . . . . . . .
1.3.2 Multiple Charge State Transport
1.3.3 Quarter Wave Resonator Steering

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

1
3
5
6
8
9
10
11
13
15
15
16
18
18
20
21

Chapter 2 Review of Beam Physics . . . . . . . . . . . . . . . . . . . . . . . .
2.1 Transverse Motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Longitudinal Motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

22
24
29

Chapter 3 Modeling Elements . . . . . . . .
3.1 Quadrupoles . . . . . . . . . . . . . . . . .
3.2 Solenoids . . . . . . . . . . . . . . . . . . .
3.2.1 Hard-Edge Model Parameterization
3.2.2 Solenoid with Corrector . . . . . .
3.2.3 Solenoid Error Study . . . . . . . .
3.3 Accelerating Cavities . . . . . . . . . . . .
3.3.1 Electromagnetic Fields . . . . . . .
3.3.2 Energy Gain . . . . . . . . . . . . .
3.3.2.1 TTF Calculation . . . . .
3.3.2.2 Cavity Geometry . . . . .
3.3.3 Focusing Effects . . . . . . . . . . .
3.3.4 Benchmark Considerations . . . . .
3.3.4.1 TTF Computation . . . .
3.3.4.2 Thin-lens Modeling . . . .

35
36
42
47
52
56
58
59
61
61
64
65
67
67
68

vi

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

3.4

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

70
70
72
73
73
75
84

Chapter 4 Quarter-Wave Resonator Steering
4.1 Steering Predictions . . . . . . . . . . . . . .
4.1.1 QWR Choice . . . . . . . . . . . . .
4.1.2 Vertical Steering . . . . . . . . . . .
4.1.3 ReA Calculation . . . . . . . . . . .
4.2 Experiment . . . . . . . . . . . . . . . . . .
4.2.1 Procedure . . . . . . . . . . . . . . .
4.2.1.1 Vertical Offsets . . . . . . .
4.2.1.2 Profile Measurement . . . .
4.2.2 Beam Axis Definition . . . . . . . . .
4.2.3 Results and Analysis . . . . . . . . .
4.3 Conclusions and Outlook . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

87
88
88
91
94
97
99
100
102
106
107
112

3.5

Lattice Examples . . . . . . . . .
3.4.1 ReA Linac . . . . . . . . .
3.4.2 FRIB Linac Segment 1 . .
Multiple Charge State Transport
3.5.1 Background . . . . . . . .
3.5.2 Methodology . . . . . . .
3.5.3 Benchmark . . . . . . . .

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

Chapter 5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
5.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
5.2 Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116
APPENDICES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
Appendix A Additional Cavity Measurement . . . . . . . . . . . . . . . . . . . . 118
Appendix B L092 Drive Calibration . . . . . . . . . . . . . . . . . . . . . . . . . 120
BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124

vii

LIST OF TABLES

Table 3.1

Analysis Results of Quad Scan Measurement Data (2-24-11) . . . . .

Table 4.1

Summary of cavity steering data at β = 0.036, 0.043, 0.046, and 0.051. 111

viii

41

LIST OF FIGURES

Figure 1.1

The National Superconducting Cyclotron Laboratory (NSCL) in 2012.

3

Figure 1.2

The NSCL layout highlighting areas of production, fast beams, gas
stopping, stopped beams in the low energy area, and reaccelerated
beams from ReA. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

ECR schematic, including solenoid and hexapole coils. Microwave
power, ωRF , creates the resonance zone in the center, and ion extraction
is indicated to the right. . . . . . . . . . . . . . . . . . . . . . . . .

5

Cyclotron Schematic, depicting simplified beam path in red, alternating charge on “dees” with beam’s relative position, and magnetic field
into the page. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

ISOL production yields at ISOL(DE), the Isotope mass Separator
On-Line facility at CERN [6], compared with PF production yields
at the NSCL [7]. Note that there are bands of isotopes that are
inaccessible in the ISOL production method due to its chemically
dependent nature, while projectile fragmentation does not have this
limitation. However, ISOL can produce a very intense beams with a
higher yield than PF for certain elements. . . . . . . . . . . . . . . .

8

A1900 schematic, including four dipole magnets with maximum magnetic rigidity of 6 Tm. Production target, intermediate image, and
focal plane locations also indicated. . . . . . . . . . . . . . . . . . .

9

ReA layout including the accelerating linac on an elevated platform
and the experimental stations on the ground level. The elevated deck
design facilitates continuous access for laboratory personnel to other
experimental areas. . . . . . . . . . . . . . . . . . . . . . . . . . . .

13

A proposed layout for Energy Upgrade Phase I: ReA6. Note that additional single-cavity cryomodules have been added downstream from
the original linac for longitudinal bunching, as well as an additional
β = 0.085 cryomodule for further acceleration. New experimental
detectors will also be commissioned with the construction of ReA6. .

16

Facility for Rare Isotopes (FRIB) layout. . . . . . . . . . . . . . . .

17

Figure 1.3

Figure 1.4

Figure 1.5

Figure 1.6

Figure 1.7

Figure 1.8

Figure 1.9

ix

Figure 2.1

Coordinate System . . . . . . . . . . . . . . . . . . . . . . . . . . .

23

Figure 2.2

Phase space ellipse with negative α (defocusing). . . . . . . . . . . .

26

Figure 2.3

Phase space transformation of a converging beam along a drift space.
After achieving a waist, the converging beam becomes diverging. . .

27

Longitudinal motion energy gain is proportional to cos(φ). Points
of zero energy gain (rebunching mode), stable acceleration (negative
phase), and unstable region (positive phase) indicated. . . . . . . . .

30

Basic FODO cell consisting of focusing (QF) and defocusing (QD)
quadrupoles, separated by drift spaces of length L. In this example,
the FODO cell begins at the center of the first focusing quadrupole,
ending in the center of the next focusing quadrupole. . . . . . . . .

36

Left (a): Equipotential surfaces of quadrupole fields plotted for several different voltages. Right (b): The electrodes inside an actual
electrostatic quadrupole are shaped along equipotential surfaces in
order to produce a pure quadrupole field. . . . . . . . . . . . . . . .

37

Left (a): Quad Scan Measurements performed before RFQ at ReA.
Solid lines indicate measurement data, while dashed lines indicate
SVD analysis. Right (b): Analyzed viewer data for L057QB=−5.5 kV
with overlay of beam position and RMS beam size in green. Beam
profiles in red are shown above and to the left of the beam spot. . .

41

Equivalent hard-edge model for LEBT 0.6 T solenoid located before RFQ. The measured longitudinal field profile is shown with
Bmax =0.6926 T, while the equivalent hard-edge model consists of a
drift section (B =0 T), a constant field section (B =0.6010 T), followed
by a second drift section (B =0 T). . . . . . . . . . . . . . . . . . . .

43

Figure 2.4

Figure 3.1

Figure 3.2

Figure 3.3

Figure 3.4

Figure 3.5

The ReA solenoid L090SN scaled to 1T for Bρ = 0.5 T m. The
solution that minimized Mtot is found to be at Le =0.24 m, Be =0.8 T. 48

x

Figure 3.6

Figure 3.7

√
XAL simulation of transverse beamsizes,
x2 and
y 2 , σ11x
√
and σ11y as defined in Equation 3.2, for the ReA Linac without
the third cryomodule using the solenoid parameterization defined
in Equation 3.25. Vertical dashed lines indicate diagnostic locations.
Note that the beam is focused in both x and y at the solenoid locations,
and the beam sizes at the third and fourth diagnostic locations, at
∼6 m and ∼11 m respectively, are of similar magnitude, which was
observed experimentally. . . . . . . . . . . . . . . . . . . . . . . . .

50

Numerical method results of Be vs. Le shown calculated at various
Bρ and B0 values that cover the relevant ReA ranges. The red regions
indicate areas with high divergence between the hard-edge and fieldprofile models, while the blue regions indicate areas where the two
models agree. Note that solutions are present (blue regions exist) for
low beam rigidities, Bρ, only at low solenoid fields, B0 . The solution
range increases to higher fields at higher rigidities. . . . . . . . . . .

51

Figure 3.8

Percent difference between the “full-field” and “thin-lens” models for
beamsizes calculated at the solenoid exit, z = 0.2 m, as a function of
relative corrector field strength, By , to the main solenoid field, Bz = 2 T. 55

Figure 3.9

FRIB Linac Segment 1 (XAL) with random Gaussian error added to
each solenoid. Gaussian errors of about (a)0.5% and (b)1% have been
used — note how the RMS of the relative difference in σx,y grows by
approximately the percentage of the Gaussian error introduced to the
system. It is an important result that this growth levels out without
increasing exponentially. . . . . . . . . . . . . . . . . . . . . . . . .

57

Cut-away image of the ReA β = 0.041 QWR cavity type. The two
accelerating gaps that the beam travels through towards the bottom
of the cavity are shown in detail. . . . . . . . . . . . . . . . . . . . .

58

ReA β = 0.041 Cavity Fields, including the accelerating electric field,
Ez , and the transverse fields Bx and Ey . Note that Ey is very small,
and so therefore has been multiplied by ten for the purposes of this
graph. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

60

Figure 3.12

TTF of β = 0.041 and β = 0.085 cavities . . . . . . . . . . . . . . .

63

Figure 3.13

QWR and HWR Cavities [29]. The QWR cavities operate at 80.5 MHz,
while the HWR cavities operate at 322 MHz. . . . . . . . . . . . . .

64

Electric field lines in an RF gap. . . . . . . . . . . . . . . . . . . . .

65

Figure 3.10

Figure 3.11

Figure 3.14

xi

Figure 3.15

Figure 3.16

Figure 3.17

Figure 3.18

Figure 3.19

Figure 3.20

Figure 3.21

Figure 3.22

Figure 3.23

ReA T T F for β = 0.041 and β = 0.085 cavities, with operating
regimes highlighted. The β = 0.041 cavities are located in the first
two cryomodules, CK01 and CL01, while the β = 0.085 cavities are
located in the third cyromodule, CM01. . . . . . . . . . . . . . . . .

68

Modeling two gap acceleration in XAL as compared with IMPACT
for a β = 0.041 accelerating cavity. . . . . . . . . . . . . . . . . . . .

69

Benchmarking results of the ReA linac through its second cryomodule
show good agreement between the matrix modeling code XAL, and
the multi-particle tracking code IMPACT. Plots show the results for
transverse β functions (a) βx and (b) βy , longitudinal β function (c)
βz , and finally the kinetic energy per nucleon, (d) Ek . . . . . . . . .

71

XAL Cavity Benchmark with IMPACT for a U 33+ beam through
FRIB’s Linac Segment 1 to the Stripper Foil, including modeling of
β = 0.041 and β = 0.085 cavities. . . . . . . . . . . . . . . . . . . .

72

Results from ANL multiple charge state transport: pepper-pot image
of the combined bismuth beam with charge states 20+ and 21+ (on
the left) and superimposed images of 209 Bi20+ (blue) and 209 Bi21+
(red) beams transported individually (on the right) [50]. . . . . . . .

74

Multiple Charge State Energy Gain: ∆φ = 0 on left, ∆φ = ωRF ∆t on
right. As discussed in the text, modeling the separate charge states
with ∆φ = 0 yields the unphysical result where each charge state
arrives at the accelerating phase of the reference charge state, resulting
in higher charge states receiving higher energy gains and lower charge
states receiving lower energy gains than the reference charge state.
Instead, each individual charge state’s arrival time at each cavity must
be calculated in order to correctly model the multiple charge state’s
energy gain, as seen here on the right. . . . . . . . . . . . . . . . . .

77

FRIB Linac Segment 1 as simulated by XAL for U 33+ and U 34+
charge states, as well as the result of the two charge states combined.

81

End of FRIB Linac Segment 1 chicane as simulated by XAL for U33+
and U34+ charge states, as well as the result of the two charge states
combined. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

82

Schematic representation of three beamlets with means µi , standard
deviations σi , and intensities Ni . . . . . . . . . . . . . . . . . . . .

83

xii

Figure 3.24

Results of XAL and IMPACT calculations for FRIB FS1 to Stripper.
The IMPACT results are shown with dashed lines, while the XAL
results are shown with the solid lines. . . . . . . . . . . . . . . . . .

85

Figure 3.25

XAL Schema for adding multiple charge state run capabilities. . . .

86

Figure 4.1

FRIB Quarter-Waver Resonator (QWR) and Half-Wave Resonator
(HWR) cavity types shown with optimal β and operational frequencies.
The FRIB HWRs are designed for four times the frequency of the
QWRs, resulting in their length being about half as long. . . . . . .

89

Coaxial Cavity Fields: the beam is accelerated at the open end
(bottom) where the electric field is maximized, while the maximum
magnetic field circulates at the closed end (top). . . . . . . . . . . .

90

The calculated ReA β = 0.041 cavity D-factors from Equation 4.6,
DBx , DEy , and Drf , with Ea = 16 MV/m and φ = −90 ◦ along with
the sum of the three terms. . . . . . . . . . . . . . . . . . . . . . .

94

Calculation of ∆y vs. β for the ReA (a) β = 0.041 and (b) β = 0.085
cavity types, plotted for centered and ±1 mm vertical offset scenarios.
A typical accelerating cavity setting is used: φ = −20◦ and Ea =
16 MV/m. The beam species is H2 with q/A = 0.5. The x-range for
each plot is from β = 0.036, corresponding to the RFQ energy of
0.6 MeV/u, to β = 0.16, corresponding to the highest beam energy
for H2 at ReA, 12 MeV/u. . . . . . . . . . . . . . . . . . . . . . . .

95

Left (a): Plot of ReA β = 0.041 calculation ∆y vs. β for offsets dy =
−6, −3, 0, 3, and 6 mm. Right (b): ∆y vs. cavity amplitude Ea for
offsets dy = −6, −3, 0, 3, and 6 mm for beam energy Ek = 0.6 MeV/u,
corresponding to β = 0.036. . . . . . . . . . . . . . . . . . . . . . . .

96

ReA Layout from 2013, including EBIT charge breeder, Q/A section,
commissioning source, room temperature RFQ and superconducting
linac including two cryomodules containing β = 0.041 accelerating
cavities. Note the drift tube after the second cryomodule is the
location for the third cryomodule containing β = 0.085 accelerating
cavities, to be installed. . . . . . . . . . . . . . . . . . . . . . . . . .

98

Figure 4.2

Figure 4.3

Figure 4.4

Figure 4.5

Figure 4.6

xiii

Figure 4.7

Schematic of ReA Linac as commissioned in 2013, before the installation of the third cryomodule. Square elements represent superconducting solenoids, while circular elements represent quarter-wave
resonators. The drift space between the diagnostics after the second
cryomodule is ∼5.2 m. . . . . . . . . . . . . . . . . . . . . . . . . . .

99

Figure 4.8

Beam position at cavity midpoint as projected from measurements
at two downstream locations. Error bars are the RMS beam size.
Note the improvement of the vertical separation and consistency of
horizontal position achieved with the 10-29 data set which utilized
the response matrix measurement technique. . . . . . . . . . . . . . 100

Figure 4.9

45◦ Beam Profile Diagnostic setup. In the top left, the plate is shown
with three slits and one aperture in the order that they intercept the
beam during the profile scans: (1) the horizontal slit measuring the
vertical beam profile, (2) the 45◦ slit for a measure of the coupling between the x and y planes, (3) the vertical slit measuring the horizontal
beam profile, and (4) the 5 mm aperture for beam tuning. . . . . . . 103

Figure 4.10

Plots of beam current vs. slit position for raw data from 45◦ slit scan
taken at (a) L092 and (b) L110, and the Gaussian fit results for beam
position and RMS beam size. Note that the 45◦ plate at L110 does
not have the 45◦ slit, which accounts for the missing peak in (b). . . 104

Figure 4.11

45◦ Beam Profile Diagnostic Plate shown with distances from Equation 4.12, dref and doffset , labeled. . . . . . . . . . . . . . . . . . . . 105

Figure 4.12

Schematic diagram of relevant measurement quantities depicting the
difference between the diagnostic axis and the cavity center axis. . . 106

Figure 4.13

Plots of ∆y vs. Ea for beam energies (a) Ek = 0.6 MeV/u (β =
0.036), (b) Ek = 0.85 MeV/u (β = 0.043), (c) Ek = 1.08 MeV/u
(β = 0.046), and (d) Ek = 1.24 MeV/u (β = 0.051). Lines correspond
to the predicted response calculated using cavity design data and
experimental conditions. . . . . . . . . . . . . . . . . . . . . . . . . 108

Figure 4.14

Left (a): Plot of ∆y /Ea vs. y at the midpoint of the test cavity,
including each measurement data set. Note the agreement in the
slopes for the different energies. Right (b): Plot of ∆y /Ea vs. β
including each measurement data set. . . . . . . . . . . . . . . . . . 110

xiv

Figure 5.1

Left (a): Raw experimental data of ∆y vs. Ea at Ek = 0.6 MeV/u
(β = 0.036). Right (b): The same data set with a constant vertical
offset, δ =1.29 mm, applied during the analysis. . . . . . . . . . . . . 115

Figure A.1

Final measurement of L091 data set on 10-31. The Ea = 0.5 MV/m
L092 data point had all peaks shifted by the same amount. The dashed
lines indicate where this overall shift has been manually accounted
for in the analysis. L110 data indicates it was an issue with the
L092 drive, not a sudden change in the beam, for instance from an
instability in the LB source. . . . . . . . . . . . . . . . . . . . . . . 119

Figure A.2

Plot of ∆y vs. Ea for beam energies Ek = 1.08 MeV/u (β = 0.046).
Lines correspond to the predicted response calculated using cavity
design data and experimental conditions. . . . . . . . . . . . . . . . 120

xv

Chapter 1
Introduction
Advanced accelerator systems today require sophisticated computer simulations to develop a
complete understanding of their operation and functionality due to the increasing complexity
of the involved systems. Real-time control and modeling, more now than ever before, are
steadily becoming essential components in the operation of accelerators in order to enhance
performance.
Accelerator physics is concerned with the creation of beams of particles, their acceleration,
and subsequent transport to experimental detectors. During transport, beam manipulation
and advanced diagnostic techniques often demand accurate modeling in order to fully understand the accelerator and beam properties to ultimately achieve the necessary requirements
for state-of-the-art experiments.
The degree of detail incorporated into these accelerator models is primarily a function of
the particular needs of the application the model is serving. There are many different types of
accelerator facilities around the world that cover a wide variety of research goals, each having
its own unique modeling challenges. Synchrotron light sources utilize high energy electron
beams in a storage ring to study topics in materials science and condensed matter physics.
High energy proton synchrotrons or linear colliders are used to study subatomic particles in
high energy physics. Low to medium energy beams of heavy nuclei are important for nuclear
physics experiments, as in the case for Michigan State University’s National Superconducting
1

Cyclotron Laboratory (NSCL).
Accelerators in nuclear physics, such as the NSCL and the Facility for Rare Isotope
Beams (FRIB), alter the particle species and energy ranges used to study the vast nuclear
landscape. These constantly changing conditions, which accelerators must be designed to
accommodate, motivate the necessity of accurate modeling methods in real-time in order
to improve accelerator tuning procedures and minimize the operational down-time when
switching between beams.
Additionally, new devices need to be properly characterized in order to be incorporated
accurately into modeling codes. It is especially important to correctly characterize the effect
that elements have on steering of the beam when considering accelerator tuning. In the case
that an effect is observed but not yet incorporated into present modeling codes, detailed
systematic measurements should be made in order to verify its correct description. One prime
example of this is the Quarter-Wave Resonator (QWR) steering experiment performed at the
ReA facility at the NSCL, which will be presented in greater detail in the following chapters.
To begin with, though, this introductory chapter presents an overview of the accelerators
that will appear throughout this dissertation, focusing on the NSCL in Section 1.1. A
background of future facilities of interest is included in Section 1.2. Within a brief summary
of the different methods used for modeling these accelerators in Section 1.3, topics which
provide unique challenges to the modeling of modern low-energy heavy ion accelerators,
including QWR steering, are introduced.

2

1.1

NSCL

The National Superconducting Laboratory (NSCL) has been a highly successful laboratory
for doing world class nuclear research over the past three decades. By commissioning the new
reaccelerator, ReA, the lab continues to push the boundaries for nuclear research and expand
its capabilities to provide rare isotopes at relevant energy regimes for experimenters. In the
future, the Facility for Rare Isotope Beams (FRIB) will replace the coupled cyclotrons with a
high power driver linac to become the preeminent facility for nuclear research in the world.
The NSCL, as shown in 2012 in Figure 1.1, is located at Michigan State University (MSU).
The NSCL currently consists of a stable Electron Cyclotron Resonance (ECR) ion source,
SuSI (the Superconducting Source for Ions), two coupled superconducting cyclotrons, the
K500 and the K1200, and the A1900 fragment separator.

Figure 1.1: The National Superconducting Cyclotron Laboratory (NSCL) in 2012.

Experimental stations at the NSCL are categorized by the energy of the beams that
they use. High energy experiments take beam directly from the A1900, while low energy
experiments are delivered beam from the gas stopping stations.
In addition to the NSCL’s coupled superconducting cyclotrons which provide the primary
beam on target, another accelerator has also been built on site to reaccelerate rare isotopes,
3

ReA. The original design for ReA includes accelerating components capable of reaccelerating
beams of uranium up to 3 MeV/u. The gas stoppers also provide beam to ReA which delivers
re-accelerated beams to its dedicated experiments.
The components of the NSCL to be discussed in further detail are shown in the schematic
layout in Figure 1.2, including the ion sources, cyclotrons, fragment separator, experimental
stations, and the reaccelerator, ReA. The highlighted regions show the areas of fast beams,
gas stopping, stopped beams in the low energy area, and reaccelerated beams.

Figure 1.2: The NSCL layout highlighting areas of production, fast beams, gas stopping,
stopped beams in the low energy area, and reaccelerated beams from ReA.

The particle beam path begins in an ion source, continuing though the cyclotrons and the
fragment separator, the production area of Figure 1.2. It then proceeds to travel to either the
fast beam area for fast beam experiments, or the gas stopping area. After the gas stopping
area the beam is transported to either the stopped beam experiments or the reaccelerator,
ReA, where it can be reaccelerated to its dedicated experiments. The reacceleration approach
allows intermediate energies to be attained, as will be explained further in Sections 1.1.6
and 1.1.7.

4

1.1.1

ECR Ion Source

Stable beams are first produced by SuSI, a superconducting Electron Cyclotron Resonance
(ECR) source capable of achieving very intense and highly charged beams, for example more
than 380 eµA of Xe26+ [1]. The SI unit for current is the Ampere (A), which is defined as a
Coulomb per second (C/S). When describing beams from ECRs, this unit is typically cited
as an electrical amp (eA), which includes the charge state of the beam. The current of the
beam as measured on a Faraday Cup, for instance, is measured in electrical Amps (eA). If
the particle current is being discussed, this is defined as the electrical current divided by the
charge state.
ECRs produce these highly intense beams by resonantly exciting electrons in the plasma
ECR
with cyclotron frequency ωRF = qBm
, where q and m are the charge and mass of an electron,

BECR is the resonant magnetic field, and ωRF is the cyclotron frequency. A schematic of the
SuSI plasma chamber is shown in Figure 1.3.

Figure 1.3: ECR schematic, including solenoid and hexapole coils. Microwave power, ωRF ,
creates the resonance zone in the center, and ion extraction is indicated to the right.

5

The microwave power heats the electrons, which in turn ionize the neutral gas or metal
vapor atoms, increasing the charge state of the desired element. SuSI operates at a microwave
frequency of fµ = 18 GHz, corresponding to a resonant magnetic field of BECR = 0.64 T.
The plasma inside an ECR is confined axially by solenoid fields: in SuSI on the injection
side Binj ≈ 4BECR and on the extraction side Bext ≈ 2BECR . The plasma is confined radially
through a hexapole field, where Brad ≈ 2BECR . All magnets in SuSI are superconducting,
which reduces the power consumption and also allows for tunable radial magnetic fields, an
advantage over permanent hexapole magnets used in room temperature ECRs [2]. Once the
desired charge state is reached, the ions are extracted from the magnetic bottle by applying
a voltage to the extraction electrodes [3] and transported to the cyclotrons for acceleration.

1.1.2

Cyclotrons

The K500 was the first ever superconducting cyclotron. A schematic of a cyclotron is shown
in Figure 1.4. The name derives from the energy that a proton beam could potentially be
accelerated to, 500 MeV. Heavy ion beams have protons and neutrons in the isotope’s nucleus
which contribute to their mass. Because of this it is common to cite an energy per nucleon
instead of the total beam energy. The K500 is capable of accelerating heavy ion beams up
to 17 MeV per nucleon, or 17 MeV/u. The total beam energy is recovered by multiplying this
value with the total number of nucelons.
The beam from the ion source is injected into the center of the cyclotron, where it is
then accelerated by an alternating electric field between the “dees”. The first cyclotron
accelerators had two electrodes, and were so named “dees” for their “D”-shaped electrodes.
Both the K500 and the K1200 have three “dees”. To utilize the acceleration through these
electrodes multiple times, a static magnetic field is applied.
6

Figure 1.4: Cyclotron Schematic, depicting simplified beam path in red, alternating charge
on “dees” with beam’s relative position, and magnetic field into the page.
mv as their energy increases.
The ions spiral through this field, with increasing radius r = |q|B

Here m, q, and v, are the mass, charge, and velocity, of the ion beam respectively. B is
the magnetic field of the cyclotron. The accelerating voltage is applied to the “dees” at the
|q|B

cyclotron resonance frequency of the beam, f = 2πm which results in the beam experiencing
an accelerating electric field through each of the gaps between the electrodes. Once the
particles’ radius reaches the radius of the cyclotron, the accelerated beam is ejected.
The K500 was the NSCL’s main accelerator from 1982 to 1989, when it was replaced with
the K1200 cyclotron which operated until 1999, when the cyclotrons were coupled [4]. This
configuration, known as the coupled cyclotron facility (CCF), allows for much higher energies
to be achieved, because a stripper foil is utilized before injection into the K1200. This yields
a higher charge state, and in turn a more efficient acceleration is achieved through the K1200
acceleration. The K1200, the larger of the two cyclotrons, can then accelerate ions to a range
of energies, from 20 (MeV/u) to 200 MeV/u [5].

7

1.1.3

RIB Production

The fast beams of rare isotopes at the NSCL are produced through a process called Projectile
Fragmentation (PF). The stable high energy (∼200MeV/u) ion beam from the cyclotrons is
impinged on a thin light production target, such as Beryllium, to produce rare (radioactive)
isotope beams (RIBs). The isotopes capable of being produced at the NSCL’s CCF is shown
on the right in Figure 1.5.

Figure 1.5: ISOL production yields at ISOL(DE), the Isotope mass Separator On-Line facility
at CERN [6], compared with PF production yields at the NSCL [7]. Note that there are
bands of isotopes that are inaccessible in the ISOL production method due to its chemically
dependent nature, while projectile fragmentation does not have this limitation. However,
ISOL can produce a very intense beams with a higher yield than PF for certain elements.

The nuclear reactions inside the target foil at relatively high energies remove a significant
amount of nucleons to produce a shower of particles with less mass than that of the projectile
and target nuclei. A fraction of the beam reacts with the target, producing rare isotopes
lighter than the beam element in a statistical process. By using a heavy beam of an element
such as Uranium, many rare isotopes can be produced this way [8].
Another production method for rare isotope beams is called Isotope Selection On-Line, or
ISOL, which is a chemically dependent process in which a high energy beam is impinged on a
8

thick target where the rare isotopes diffuse out of the target and are then subsequently accelerated. This process can take on order of microseconds, so for rapidly decaying beams projectile
fragmentation is the necessary production method. Because of ISOL’s chemically-dependent
nature, it may also be necessary to produce certain beams with projectile fragmentation, as
seen in Figure 1.5.
Projectile fragmentation (PF) has the advantage that it is a chemically independent
process and avoids this limitation of ISOL production. However, as PF is a statistical process,
with RIBs farther from stability produced at slower rates, separation of the many different
reaction products is necessary in order to form a beam of one desired species.

1.1.4

Fragment Separator

The A1900 fragment separator at the NSCL, seen in Figure 1.6, uses several methods to
efficiently separate a desired isotope beam from the many reaction products created through
the PF production method [9]. Its name is derived from “A” for Analysis, and 1900 for the
maximum proton kinetic energy in MeV that can be bent by its four dipole magnets.

Figure 1.6: A1900 schematic, including four dipole magnets with maximum magnetic rigidity
of 6 Tm. Production target, intermediate image, and focal plane locations also indicated.
9

The first method for isolating the desired RIB utilized by the A1900 is magnetic separation.
The dipole magnets disperse the beam according to its momentum to charge ratio, defined as
the magnetic rigidity, Bρ = p/q. Since the velocities of the different species are fairly similar,
this amounts to essentially a selection by mass-to-charge ratio.
This reduces the number of species considerably, but there are still many RIBs with
the same A/Q ratio, so further separation is needed. Selection by proton number (Z) is
performed by placing a degrader at the intermediate image location (I2 in Figure 1.6) at the
center of the A1900. As energy loss in matter is Z-dependent, as given by the Bethe-Bloch
equation [10], the velocity of the different RIBs changes as a function of their charge alone,
which allows them to then be separated out by their magnetic rigidity, Bρ.

1.1.5

Experimental Stations

After particle fragment separation, there are several different experimental stations to which
the RIB may be delivered. Experiments requiring high energies receive the beam directly from
the A1900. These experiments use various detectors including the Modular Neutron Array
and Large Multi-Institutional Scintillator Array (MoNA-LISA) [11], the High Resolution
Array Detector (HiRA) [12], and the S800 Spectrograph [13] to study nuclear reactions and
properties such as half-life, cross-sections, and branching ratios which all require fast beams.
Other experiments probe nuclear properties that require low energy beams, such as
measuring the mass as done with the Low Energy Beam and Ion Trap (LEBIT) [14], or for
nuclear magnetic resonance (β-NMR) studies, as performed with the Beam Cooler and Laser
spectroscopy (BECOLA) [15]. The fast beams from the cyclotrons must therefore have their
energy reduced before being transported to these low energy experiments.
The majority of the beam’s energy is lost through a solid degrader, followed by further
10

beam cooling in a gas cell filled with an inert gas, such as Helium. Beam cooling in this
context refers to the reduction of beam energy due to collisions in the gas which essentially
stop the beam. Once this is accomplished, the beam can then be extracted from the gas cell.
This process reduces the energy spread of the beam and its area in phase space, or emittance,
to an acceptable level for the high-precision low-energy experiments [16]. A new type of gas
stopper, called the Cyclotron Gas Stopper, is currently in the process of being assembled
and commissioned. The goal for this new gas stopper is to have higher efficiencies for lighter
beams than the current linear gas cell [17].
Other nuclear experiments, particularly those of interest in nuclear astrophysics, require
higher energies than the stopped beam experiments. The desired range of energies, from
∼0.33-3 MeV

per nucleon, is of scientific interest because these are the typical energies that

exist in stars. This demand from the nuclear science community drove the decision to
design the ReA facility to re-accelerate heavy rare isotope beams to 3 MeV/u [18]. There
also remains interesting nuclear physics to probe at slightly higher energies at or above the
Coloumb barrier, for which upgrades to ReA6 and ReA12 have also been proposed for the
future.

1.1.6

ReAcceleration Motivation

Before the construction of the ReA facility, the NSCL was capable of providing fast
(∼200MeV/u) rare isotope beams as well as low energy beams (∼1eV/u) through its gas
stopping station.
These two energy regimes are necessary for studying different nuclear reactions and
different properties of the nucleus. For example, higher energy beams are required to produce
transfer reactions. Inverse kinematics, using a heavy beam on a light target, can then be
11

used to measure the B(E2) value which is important in nuclear physics.
High energy beams are also used to study the symmetry energy term in the nuclear
equation of state, as well as many other properties of the nucleus. The S800, GRETINA and
MONA are some of the experimental stations that use high energy beams at the NSCL. On
the other hand, low energy beams are necessary for other types of high precision measurements
that are performed at LEBIT and BECOLA which measure the mass and hyperfine structure
of the nucleus respectively.
With the addition of the ReA facility, the NSCL has added a new energy region at
∼3-12 MeV/u

and now is capable of providing beams in the low, mid-range, and high energy

range for nuclear physics experiments. This mid-range region of energy is of interest especially
for nuclear astrophysics reactions as well as some other experiments that are well suited in
this region. When the facility is upgraded with additional cryomodules (in planning) the
energy range will cover physics just above the Coulomb barrier, which is again especially of
interest to nuclear astrophysics.
Fast rare isotope beams are produced and slowed down and then reaccelerated to reach the
desired energy. This process is necessary in order to preserve the emittance, or phase-space
quality, of the beam. The beam cannot be slowed immediately to the desirable mid-range
energy, as the emittance and energy spread would grow to a point where the beam is no
longer suitable for experimental conditions.
It is therefore favorable to slow down the beam in a solid degrader where the beam loses
over 95% of its energy, followed by a gas cell where collisions with inert gas atoms reduce its
energy further. The ions will not neutralize completely during the collisions with the buffer
gas atoms if the buffer gas has a very high ionization potential — typically a noble gas is
chosen for this role as its filled outer electron shell make these gases extremely inert. In most
12

cases this leaves a singly-charged beam that can then be extracted and reaccelerated. The
ReA facility also includes an Electron Beam Ion Trap (EBIT) which increases the beam’s
charge state, or charge breeds the beam, for more efficient re-acceleration.

1.1.7

ReA

The ReA facility, shown in Figure 1.7, consists of an Electron Beam Ion Trap (EBIT), a
charge to mass or Q/A selection section, a normal conducting RFQ, a superconducting linac,
and transport beam lines that deliver the beam to the experiments [19]. The beam optics
designs were developed using COSY [20] and IMPACT [21].

Figure 1.7: ReA layout including the accelerating linac on an elevated platform and the
experimental stations on the ground level. The elevated deck design facilitates continuous
access for laboratory personnel to other experimental areas.

Gas-stopped beams are typically singly ionized. In order to achieve more efficient acceleration, these beams are injected into the ReA Electron Beam Ion Trap (EBIT) for charge
breeding. An EBIT uses a magnetic bottle to confine the ions that it collects, while a
high-density electron beam increases the charge state of the ions by stripping away outer-shell
electrons, similar to the mechanism of an ECR [22]. Once the desired charge state is reached,

13

the highly charged ions are ejected from the trap.
An achromatic mass separator consisting of an electrostatic sector and a magnetic
sector provides a mass resolving power of R∼100 and energy spread of δE = 2% [23]. A
small longitudinal emittance is attained by using a Multi-Harmonic Buncher (MHB) before
acceleration [24]. The MHB excites up to three harmonics, approximating a Sine wave, in
order to bunch the beam for acceptance into ReA’s accelerating elements.
The accelerating components of ReA include a normal conducting room temperature RadioFrequency Quadrupole (RFQ) [25], and three cryomodules with two types of superconducting
quarter-wave RF cavities operating at 80.5 MHz [26]. The maximum final energy achievable
in this first stage of ReA is 3 MeV/u for heavy nuclei such as uranium, and 6 MeV/u for ions
with A ≤ 50.
The ReA facility at the NSCL has the capability to provide low energy rare isotope beams
which are typically difficult to produce at Isotope Separation On-Line (ISOL) facilities. As
mentioned in Section 1.1.3, ISOL facilities produce RIBs at rest by impinging a high-energy
stable beam on a thick production target. The rare isotopes then diffuse out of the target
where they can subsequently be accelerated to the desired energy.
Due to the individual chemistry of each species and target ion source, however, some beams
are extremely difficult to produce at the necessary rate for experiments [27]. The ReA facility
can take advantage of the chemically-independent projectile fragmentation production method
at the NSCL to offer experimenters beams at energies that were previously unavailable.

14

1.2

Future Facilities

There are many exciting developments planned for the future of the NSCL. Following is
a brief summary of these proposals for the future facilities at MSU, including upgrades to
the existing reaccelerator and an entirely new primary driver linac to replace the coupled
cyclotrons with the Facility for Rare Isotope Beams.

1.2.1

ReA Upgrades

High-quality low-energy beams are extremely desirable in the experimental nuclear physics
community. With ReA able to achieve 3 MeV/u beams, many of these low-energy reactions
for nuclear astrophysics can be studied. These reaction measurements occur at astrophysical
energies, and ultimately motivated the construction of ReA. However, as alluded to earlier,
many experiments need slightly higher energies, still well below the energies of the fast beams
at production (∼100 MeV/u).
Some of these experiments include the study of transfer reactions of exotic nuclei, multistep Coulomb excitation, and nuclear shell structure evolution far from stability, which
all require reactions with very short-lived beams at energies around 6-12 MeV/u. Energy
upgrades to ReA have therefore been planned in stages, titled ReA6 and ReA12, which will
reaccelerate the rare isotope beams to 6 and 12 MeV/u respectively. The layout for ReA6 is
shown in Figure 1.8.
The energy upgrade for ReA6 will require four additional cryomodules, three singlecavity cryomodules for longitudinal matching and rebunching, and one cryomodule with 8
accelerating cavities [28]. This cryomodule will be identical in design to the cryomodules
designed for FRIB. Drift space after this cryomodule is left during this phase for the second

15

Figure 1.8: A proposed layout for Energy Upgrade Phase I: ReA6. Note that additional singlecavity cryomodules have been added downstream from the original linac for longitudinal
bunching, as well as an additional β = 0.085 cryomodule for further acceleration. New
experimental detectors will also be commissioned with the construction of ReA6.
energy upgrade to ReA12 which requires an additional two accelerating cryomodules.

1.2.2

FRIB

The Facility for Rare Isotope Beams, FRIB, will take the place of the coupled cyclotrons at
the NSCL as the primary beam accelerator. The plans for FRIB also include a new high
power target station. After the target, the A1900 fragment separator will be reconfigured to
deliver the secondary beam to existing experimental stations as well as additional detectors
designed specifically to use beams from FRIB [29]. The configuration can be seen in the
schematic layout of Figure 1.9.
The low-β (β = v/c) portion of FRIB linac utilizes Quarter-Wave Resonators (QWRs),
similar to the accelerating cavities of ReA. This fact reinforces the importance of understanding

16

Figure 1.9: Facility for Rare Isotopes (FRIB) layout.
the performance, and characterizing the behavior, of these QWR type cavities. This topic
which will be explored in more detail in later chapters.
FRIB will increase the beam power by a factor of ten, from the ∼40kW on target at the
NSCL currently to ∼400kW on target with the final FRIB design parameters. One way to
accomplish this is to increase the current, as beam power is a product of beam current and
electric potential.
State-of-the-art ECRs are able to provide the necessary increased current, but the possibility exists to dramatically increase the beam current by accelerating multiple charge states
at once, as discussed in Section 3.5. By accelerating multiple charge states simultaneously,
the total beam current can be increased by almost a factor of five [29].

17

Many upgrade paths have also been incorporated into the design of FRIB in order to
extend its longevity [30]. There is the possibility for an upgrade to include a 1 GeV ISOL
facility that would produce certain primary beams at higher intensities than can be achieved
with the chemically-independent projectile fragmentation method.
Space has also been left vacant in the accelerator tunnel for additional cryomodules to be
added for future energy upgrades. Integrating the design of future upgrades into the project
plan ensures FRIB will continue to be at the forefront of nuclear physics research for years to
come.

1.3

Modeling Considerations

There are unique challenges that accelerators, such as ReA and FRIB, face which require
both in-depth modeling for detailed studies and fast online modeling to facilitate beam
tuning. Some typical accelerator modeling methods are introduced in Section 1.3.1, followed
by the introduction of some of these unique challenges: multiple charge state transport in
Section 1.3.2 and Quarter-Wave Resonator (QWR) steering in Section 1.3.3.

1.3.1

Modeling Methods

There are different methods available for modeling accelerators, including multi-particle
tracking codes and beam envelope calculations using matrix based transport codes. The
appropriate modeling method depends on the beam properties being studied.
Tracking individual particles through magnetic and electric fields of elements of the
accelerator takes the result from those hundreds of thousands of particle runs together to
simulate the beam. This multi-particle tracking method is excellent for designing machines

18

as it considers higher order effects, and can take into consideration misalignments and other
errors and uncertainties.
Instead of tracking particles individually, the beam can be considered as a “macro-particle”
in which the calculated envelope is actually the RMS extent of particles that make up the
beam. By considering the beam statistically, its various statistical moments can be propagated,
some of which in turn are observables that can be measured with beam diagnostics. Each
element in the accelerator is represented by a first order matrix which acts upon the beam.
The matrices are multiplied together to get a model of the full accelerator.
Before choosing a modeling method, it is always important to first establish that the
relevant physics phenomena are being appropriately described. Only once the correct physics
is incorporated into the model can accurate results expect to be gained. To this end,
it is essential that the limitations of each modeling method is understood. Small-angle
approximations or linearity assumptions are often used to simplify complicated systems.
These descriptions are often valid when characterizing first-order effects, but care must be
taken for interpretations not to exceed the limits that each model presents.
Computational resources and restrictions are another factor to consider when choosing an
appropriate modeling method. Multi-particle tracking codes are able to provide a very detailed
description of the beam, but they also have a proportionally high demand computationally.
While envelope or single particle tracking may provide less detail, they are very fast methods
that can be used in on-line modeling applications. This is especially important when it is
necessary to model multiple different ion species through the accelerator, as with the NSCL
and FRIB.
Many beam properties can be studied with a simple linear model, but as mentioned
previously, some phenomena require more detailed modeling. For instance, higher order effects
19

such as field non-linearities and beam-beam effects are typically studied with multi-particle
tracking codes which can model 3-D fields, while matrix-based transport codes are sufficient
for studying linear motion.
The focus of this dissertation will be based in matrix-based transport modeling, as
presented in Chapter 2. Many elements can be studied with this method of modeling, and
benchmarking efforts confirm its accuracy when compared with multi-particle tracking codes,
as shown in Chapter 3. The modeling codes used throughout include XAL [31], IMPACT [32],
and COSY [33]. XAL is a linear matrix-modeling code, IMPACT is a particle tracking code,
and COSY is a higher-order matrix transport model.

1.3.2

Multiple Charge State Transport

The high intensities required on target in order to produce the desired beam rates of extremely
rare isotopes for nuclear experiments at facilities such as FRIB may be difficult to achieve
even using state-of-the-art ECR (electron cyclotron resonance) ion sources. The transport of
two or more charge states simultaneously provides a method by which the high intensities
required by FRIB may be attained; however, this is a great technical challenge to attempt,
and requires much study and simulation in order to ensure that the multiple charge state
transport scheme is practical and realistic in its objectives.
Multiple charge state beam transport has typically been modeled through multi-particle
tracking codes which can take a considerable time to run. In order to model multiple charge
state beams in real-time with an online model, a fast modeling method must be established
and benchmarked. This topic will be examined in more detail in Chapter 3 along with a brief
summary of the history of multiple charge state transport.

20

1.3.3

Quarter Wave Resonator Steering

Another challenge facing linac accelerators of heavy ions is to include the steering from
Quarter Wave Resonators into the modeling of these cavity types. The asymmetry of the
QWRs leads to a vertical kick of the beam, even when the beam enters the cavity centered
on-axis. Despite this steering effect, QWRs are the preferred cavity type for low-β acceleration,
as will be discussed further in Section 4.1.1.
With respect to accelerators, “steering” is typically thought of as originating from field
errors, or possibly misaligned beam line elements. In the case of QWRs used in the low-β
portion of the linacs of ReA and FRIB, the steering is inherent in the design of the cavity,
meaning even on-axis particles will be steered. Although design considerations can be made
to mitigate this effect, simpler cavity geometries may still be preferable for the advantages
they provide by way of cost and time saving in production. Moreover, these design alterations
may not entirely eliminate this effect for all beam conditions and velocities. Therefore a
predictive model to study this effect is desired.
Currently this effect is counteracted through the use of vertical and horizontal correctors
installed before and after the cavities. The steering had not yet been measured experimentally
to verify the theoretical predictions, or added to online models for tuning purposes. The first
measurement data of this effect is presented in Chapter 4, with concluding remarks given in
Chapter 5.

21

Chapter 2
Review of Beam Physics
Fundamental accelerator physics concepts are briefly introduced in this chapter for reference
in further discussions. Transverse and longitudinal beam physics principles are presented,
followed by a summary of the modeling methods used throughout this work.
In modeling the behavior of charged particles through an accelerator, the motion is
determined by the Lorentz force which the particles experience when passing through the
electromagnetic fields of the beam line elements. The equation of motion for a particle of
charge q with velocity v in fields E and B is therefore given by the following formula:

dp
=q E+v×B .
dt

(2.1)

The electric and magnetic fields must satisfy Maxwell’s equations, where boundary conditions
are given by the beam line composed of electromagnetic elements. It is often useful to
describe the motion of particle beams in terms of their transverse motion, perpendicular to
the direction of transport, and longitudinal motion, along the direction of transport.
In Figure 2.1, transverse motion corresponds to momentums px and py , while longitudinal
motion corresponds to pz . A moving frame of reference with respect to a reference particle
along the ideal trajectory, s, is used for a coordinate system.
These coordinates define the position of a particle in an accelerator as a six dimensional

22

yˆ

sˆ

zˆ

particle

xˆ

ref. particle
Figure 2.1: Coordinate System
vector as given in Equation 2.2, where primes denote derivatives with respect to s, and
δ = p∆p , the normalized difference in longitudinal momentum of the particle with respect to
ref
the reference particle.













U =









x 


x 



y 



y 



z 


δ

(2.2)

With the coordinates established, particle dynamics can now be explored. Independent
transverse and longitudinal motions are introduced in Sections 2.1 and 2.2 before considering
coupling between the transverse degrees of freedom, such as in the case of solenoids, discussed
in Section 3.2. A full six dimensional treatment is necessary in the case of elements that
couple the longitudinal and transverse planes, such as bending dipoles, and an additional
seventh dimension is introduced in Section 3.5. Throughout this work the dimensionality of
each problem will be considered and motivated on a case by case basis.
23

2.1

Transverse Motion

Solving for the transverse linear motion of particle beams yields a solution that takes the
form of a harmonic oscillator with its damping constant dependent on the position s in the
beam line, as derived in [34]. This is exactly the form of the solution to Hill’s equation, a
linear second-order differential equation that has been studied since the nineteenth century.

x + K(s)x = 0

(2.3)

Hill’s equation, given by Equation 2.3, essentially describes a harmonic oscillator with a
position-dependent spring constant, K(s). This motivates the choice to let x take the following
form, where ω and ψ are functions of s:

x(s) = A1 ω(s) cos(ψ(s)) + A2 ω(s) sin(ψ(s)).

(2.4)

Solving for x and plugging into equation 2.3 results in

(A1 C2 + A2 C1 ) cos(ψ) + (−A1 C1 + A2 C2 ) sin(ψ) = 0

(2.5)

where

C1 =

ωψ + 2ω ψ

C2 = ω − ω(ψ )2 + Kω.

(2.6)
(2.7)

In order for the formula to remain valid for all values of s, the coefficients in front of the
Cosine and Sine terms in Equation 2.4 must independently go to zero. Requiring nonzero A1
24

and A2 constants leads to the conclusion that the expressions C1 and C2 be zero. Solving for
C1 = 0 in Equation 2.6 results in the formula

ψ =

k
ω2

(2.8)

where k is a constant of integration.
We define β(s) = 1/ψ as the betatron function, along with α = − 12 β and γ = (1 + α2 )/β.
The choice for these definitions becomes evident when solving C2 = 0 in Equation 2.7. In
terms of β and γ, the solution is simplified to

β + 2Kβ − 2γ = 0.

(2.9)

Plugging in α = − 12 β , Equation 2.9 becomes

α = Kβ − γ.

Making a change of variables to express x in terms of constants

(2.10)

√

and δ, the equation of

motion in Equation 2.4 and its derivative can be written as

x=
x = −

√

β cos(ψ + δ)

β (α cos(ψ

+ δ) + sin(ψ + δ))

(2.11)
(2.12)

where the previously defined expressions for β and α have also been used. The choice for
the form

√

becomes clear when the results from Equations 2.11 and 2.12 are plugged into

the expression γx2 + 2αxx + βx 2 . With much simplification, it can be shown that this

25

combination results in Equation 2.13.

γx2 + 2αxx + βx 2 =

(2.13)

This is the general equation for an ellipse, which results in the physical implication that
in the case of periodic lattices, on subsequent revolutions the particle will trace out an ellipse
in transverse phase space (xx or yy ) of area A = π , with maximum extents xmax =
√
and xmax = β .

√

γ

Figure 2.2: Phase space ellipse with negative α (defocusing).

Figure 2.2 shows the ellipse in phase space. Solving Equation 2.13 for the intersections
at the x- and x - axes gives xint =

/γ and xint =

/β. The area, π , is known as the

emittance of the beam — the phase space area of the beam — and is typically given in units
of πmm·mrad. The conservation of phase space area in an inertial reference frame is known
as Liouville’s Theorem [35].
This parametrization of the beam with β, γ, α, and , was derived by Courant and Snyder
in their theory of strong-focusing synchrotrons [36]. These parameters may also be referred to

26

as Twiss parameters [37]. These parameters are important and convenient tools for describing
the state of the beam.
The evolution of the beam in transverse phase space through a drift space is shown in
Figure 2.3. A converging beam, with α > 0, reaches a beam waist, where the beam is upright
in phase space with α = 0. After this point it diverges, with α < 0. Here the beam is
spreading out in position, x or y, and a focusing element would be necessary to maintain
beam transport.

Figure 2.3: Phase space transformation of a converging beam along a drift space. After
achieving a waist, the converging beam becomes diverging.

Through certain elements, quadrupoles and drift spaces for example, the transverse phase
spaces, xx and yy , decouple, allowing the four dimensional transverse phase space to be
reduced to a two dimensional problem. The two-dimensional transformation of the transverse
coordinates defined in Equation 2.2 between U0 and U1 can be written as a 2×2 matrix, M ,
as shown in Equation 2.14. This matrix is referred to as a beam transport matrix.








 x0 
 x1 
U1 = M U0 ⇒ 
=M

x1
x0
T

(2.14)

T

By multiplying Equation 2.14 with its transpose, U1 U1 = M U0 (M U0 )T = M U0 U0 M T ,
27

it is clear that the 2×2 matrix defined in Equation 2.15,




 x 
UUT = 

x



x2

x x



xx 

2
xx x


=

(2.15)

can also be transformed by the transport M matrix. Defining the RMS of this 2×2 matrix,
as given in Equation 2.16, it follows that the σ matrix transforms as σ1 = M σ0 M T .


σ=

x2
xx


xx
x

2




From σ = U U T we find that U T σ −1 U = ✶. As xx

(2.16)

= x x , σ12 = σ21 , σ is a

symmetric matrix. This allows for simplification when multiplying out U T σ −1 U to result in
Equation 2.17.
σ22 x2 − 2σ12 xx + σ11 x 2 = det σ
Multiplying Equation 2.13 through by

(2.17)

and comparing term by term with Equation 2.17,

it is found that the two-dimensional σ matrix can be written in the form of Equation 2.18,




 β −α 
σ= 

−α γ

(2.18)

with σ transforming as σ1 = M σ0 M T . This 2-D σ matrix describes the beam in two
dimensions as a phase space ellipse seen in Figure 2.2. The full 6-D representation of the
beam is sometimes referred to as a hyperellipsoid, with the σ matrix elements therefore
representing the hyperellipsoid coefficients.
Another useful description of the beam utilizes the phase advance, ∆ψ, of the beam

28

between two points. By inserting initial and final conditions, x0 , x0 and x1 , x1 respectively,
defined by Equations 2.11 and 2.12 into Equation 2.14, after much derivation it can be shown
that for linear motion, the transport matrix M takes the following form:

M=



cos(∆ψ) + α0 sin(∆ψ)
β0 sin(∆ψ)
β1 



β0
1 ((α − α ) cos(∆ψ) − (1 + α α ) sin(∆ψ)) cos(∆ψ) − α sin(∆ψ)
0
1
0 1
1
β
1

(2.19)
where ∆ψ = ψ(s1 ) − ψ(s0 ) is the phase advance of the beam between U0 and U1 [38]. This
form is particularly useful in period lattices.
Until this point, only transverse beam motion has been discussed. Following is an
introduction of motion in the longitudinal plane, along the direction of beam propagation.

2.2

Longitudinal Motion

In the longitudinal plane, the reference particle is also known as the synchronous particle
because it is synchronized with the RF frequency of the accelerating fields. During acceleration,
this particle arrives at a given point with the exact design velocity for the accelerator, while
particles with slightly less energy arrive later, and those with slightly more energy arrive
earlier.
Stable acceleration phases are located before the peak of the RF amplitude, as labeled
in Figure 2.4. In this region, particles that arrive earlier (i.e. with more energy) than the
synchronous particle experience a smaller accelerating field, while particles that arrive later
(i.e. with less energy) than the synchronous particle experience a larger accelerating field. In
this way, particles are kept oscillating around a stable phase and longitudinal focusing, or
phase stability, is achieved.
29

Figure 2.4: Longitudinal motion energy gain is proportional to cos(φ). Points of zero energy
gain (rebunching mode), stable acceleration (negative phase), and unstable region (positive
phase) indicated.
On the opposite side of the crest of the RF, the situation is reversed. Now earlier particles
will receive more acceleration while later particles will receive less acceleration, creating a
region of instability and longitudinal defocusing.
A third point is labeled in Figure 2.4 at the x-intercept at φ = −90 ◦ , the “bunching”
phase. At this point, the reference particle does not gain any energy. Later particles gain
energy, and earlier particles actually lose energy (negative portion of the curve), resulting in a
longitudinal “bunching” effect for which this point is named. The x-intercept on the unstable
side of the curve would also correspond to a point where the reference particle does not gain
any energy. Here, though, unlike the stable “bunching” point, earlier particles would gain
energy and later particles would lose energy, spreading the beam out longitudinally. This
point is therefore not used in typical beam transport situations, but may be utilized in cases
where longitudinal phase space manipulation is desired [39].
30

Similarly to transverse motion, the longitudinal motion of the beam can also be linearized
in cases where higher-order effects can be neglected. Accelerator modeling codes have several
different choices of longitudinal coordinates, some of which will be presented here.
The coordinates used to describe longitudinal motion are often Energy and Time, as these
are canonical variables [40]. These variables are often referred to as the change in Energy
and change in Time, ∆W and ∆t, defined with respect to the reference particle. These
coordinates conveniently lend themselves as a natural basis for longitudinal motion as their
phase space emittance is conserved according to Liouville’s theorem. This basis is especially
convenient when describing bunch length, or energy spread, as the widths are directly evident
from the longitudinal phase space diagram.
Another longitudinal coordinate basis of canonical variables is z and pz . Conventionally
the coordinates will often be defined as z and z , or z and δ = ∆p/pref , as seen in Equation 2.2,
but it is important to remember that these bases are not strictly canonical variables.
One matrix modeling code, XAL, available as open source code through the Open XAL
project [41], uses z and z where z = 2∆p . Here γ0 refers to the relativistic factor, not
γ p
0 ref

to be confused with γ from the Courant-Snyder parameters defined previously. There is a
simple transformation between XAL’s coordinate definition and the z and δ basis, as shown
in Equation 2.20.










z




  1 0  z 
=


2
∆p/pref
0 γ0
z

(2.20)

It is easy to arrive at this transformation, as the same first coordinate, z, is used in
both definitions, and z is simply a function of δ: z = δ2 = 2∆p . More complicated
γ0
γ0 pref
transformations are necessary when the longitudinal coordinates are more dissimilar.
The modeling code COSY [33], for example, uses the more esoteric longitudinal coordinates
31

γ

0 , and ∆W/W is the relative change in energy. This basis
l and ∆W/W , where l = −∆tv0 (1+γ
)
0

is chosen to give a coordinate along the beam line, l, that accounts for the relativistic factor
as the longitudinal coordinate points in the direction of propagation. The transformation
between this definition and the z and δ basis is slightly more complicated, as shown in
Equation 2.21.




 z


 
=

∆p
pref



γ0 +1
γ0

0

1
γ0 (γ0 +1)
γ0
γ0 +1




l




∆W
W




(2.21)

This can be shown in the following manner. Beginning with the definitions for p and
W , p = βγm0 c and W = (γ − 1)m0 c2 , it can be shown that for a change in time, ∆t,
Equations 2.22 and 2.23 hold true, noting here that p, W , γ and β are all functions of time.

∆p = γ02 p

∆β
β0

∆W = γ0 β02 m0 c2

(2.22)

∆p
pref

(2.23)

The relative change in W , ∆W/W is found by dividing ∆W in Equation 2.23 by W =
(γ0 − 1)m0 c2 . This result can be simplified to yield the relationship shown in Equation 2.24.
∆p
γ0 ∆W
=
pref
γ0 + 1 W

(2.24)

This completes the derivation for the second row of matrix Equation 2.21. For the first
row, across the change in time, ∆t, z is defined by z = −∆t(v0 + ∆v). Substituting in the
c∆β
∆β
definition for ∆v from v = βc ⇒ ∆v
v0 = cβ0 = β0 leads to Equation 2.25.

z = −∆tv0 (1 +

∆β
) = −∆tv0
β
32

1 ∆p
1+ 2
γ0 pref

(2.25)

Plugging in the definition from Equation 2.24 for ∆p/pref as well as the COSY definition
of −∆tv0 =

l(1+γ0 )
γ0

into Equation 2.25 along with some rearranging of terms results in

Equation 2.26.

l=



γ0
z
1 + γ0
1+

1
1
∆W
γ0 (1+γ0 ) W



(2.26)

1 = 1 + x + . . . can be used to simplify this result, making the
The Taylor series 1−x
γ

∆W is small. This yields l ≈
1
∆W ), which can
1
0 z(1 −
assumption that γ (1+γ
1+γ
W
)
γ
(1+γ
)
0
0
0
0
0 W

be rearranged to form Equation 2.27.

z=

γ0 + 1
γ0

l+

1
γ0 (1 + γ0 )

∆W
W

(2.27)

This is exactly the first row as given in the matrix Equation 2.21, completing this
derivation. As previously stated, this matrix transforms the longitudinal coordinates between
the two bases, z and ∆p/pref , and l and ∆W/W .
If the transport matrix, M , is expressed in one basis, it also must also be transformed
in order to be used in the different basis. Defining T the transformation matrix of the
coordinates given in Equation 2.21, the following transformation is necessary to transform
the transport matrix.
MCOSY = T −1 Mδ T

(2.28)

In Equation 2.28 Mδ refers to a transport matrix expressed in terms of the (z, ∆p/pref ) basis, while MCOSY refers to the transport matrix in terms of the COSY longitudinal coordinates.
Again, the transformation matrix T is defined by Uδ = T UCOSY .
This transport matrix transformation definition can be shown in the following manner.
Beginning with Uδ2 = Mδ Uδ1 , plug in the definition Uδ = T UCOSY . This results in T UCOSY 2 =

33

Mδ T UCOSY 1 . Multiplying both sides by T −1 yields Equation 2.29.

UCOSY 2 = T −1 Mδ T UCOSY 1

(2.29)

This can be compared with the definition UCOSY 2 = MCOSY UCOSY 1 to arrive at the result
shown in Equation 2.28. In this way the transport matrices in different bases may be arrived
at in a straightforward manner once the transformation between coordinate bases is known.
When comparing modeling results across these different codes these appropriate transformations must be made in order to be consistent. It should also be noted that although some
bases are used as they are more convenient in terms of calculations, such as the case with the
momentum coordinates z and δ, and z and z , these coordinates are not canonical variables
and therefore their phase space emittance will not be conserved.
The next chapter will focus on the matrix modeling of several accelerator elements,
including quadrupoles, solenoids, and accelerating cavities. To reinforce the validity of these
models, benchmarking results with multi-particle tracking code IMPACT are also included.

34

Chapter 3
Modeling Elements
The ability to accurately model beam line elements is fundamental in the field of accelerator
physics. Once a model for each element is established, the combination of these individual
elements forms a model of the entire beam line. Developing further studies with these
individual models promotes a deeper understanding of the beam’s behavior through the
accelerator, and also produces opportunities to measure basic beam properties, such as
emittance.
Such is the case for quadrupole emittance measurements which use our understanding of
modeling quadrupoles in a beam line to manipulate the beam in such a way to derive the
beam’s emittance through a simple procedure. The quadrupole model and the results from
emittance measurements taken with the quadrupole doublets before the RFQ at ReA are
presented in the first section of this chapter.
Following is a detailed discussion of the investigation into two of the main elements that
make up the ReA linac, the solenoid and the Quarter-Wave RF cavity. Solenoidal fields
provide focusing, while RF cavities provide acceleration, or, in some cases, deceleration. Each
element is addressed in its own section, and results from modeling both elements together
are also included. Also discussed is the modeling of combined solenoid and corrector magnets
and studying the addition of errors, while further investigation of the RF cavities is left to
Chapter 4.
35

This chapter concludes with a brief history of multiple charge state transport, as this is
an important emerging technique for increasing beam intensity beyond the current limits of
state-of-the-art ion sources. Matrix-based codes have not traditionally incorporated multiple
charge state modeling functionality, motivating the matrix-based envelope method presented
later in the chapter. This method greatly reduces the calculation time required compared
with multi-particle tracking, allowing for on-line based modeling of multiple charge states.

3.1

Quadrupoles

Quadrupole magnets are important elements in transport beam lines–in fact, they make up
one of the most basic units for beam transport, forming a ”‘Focus-Drift-Defocus-Drift”’, or
FODO, lattice [38]. This common lattice assembly of alternating focusing and defocusing
quadrupoles allows for focusing to be accomplished in both transverse planes. The basic
FODO cell is shown in Figure 3.1.

Figure 3.1: Basic FODO cell consisting of focusing (QF) and defocusing (QD) quadrupoles,
separated by drift spaces of length L. In this example, the FODO cell begins at the center of
the first focusing quadrupole, ending in the center of the next focusing quadrupole.

The ability to utilize strong focusing and furthermore the ability to decouple the x
and y transverse planes make quadrupoles an essential component in beam transport lines.
Quadrupoles focus in one plane, and defocus in the other, without coupling the two planes.
36

+

−

−

+
(b)

(a)

Figure 3.2: Left (a): Equipotential surfaces of quadrupole fields plotted for several different
voltages. Right (b): The electrodes inside an actual electrostatic quadrupole are shaped
along equipotential surfaces in order to produce a pure quadrupole field.
Hyperbolic equipotential surfaces uniquely determine the quadrupole field, as shown
in Figure 3.2. Ideal quadrupole fields are linear near the beam axis. Quadrupoles can be
magnetic or electrostatic elements; magnetic quadrupole fields take the form of Bx = Gy,
By = Gx, and Bz = 0, where G is the focusing gradient. Electrostatic quadrupole fields
take the form of Ex = 2V x/a2 , Ey = −2V y/a2 , and Ez = 0, where V is the voltage on the
hyperbolic electrodes, and a is the minimum radial aperture. As both produce transverse
forces proportional to the distance from the central axis, they both have the same resulting
transport matrix form, given in Equation 3.1. This matrix represents the case of a quadrupole
of length L and focusing strength k, focusing in the xx plane.


sin(kL)
k

0
0
 cos(kL)


 −k sin(kL) cos(kL)
0
0

=

Mquad

sinh(kL)

0
0
cosh(kL)
k


0
0
k sinh(kL) cosh(kL)

37













(3.1)

Notice that the 2 × 2 submatrix in the focusing plane, xx , contains sines and cosines,
while the 2 × 2 submatrix in the defocusing plane, yy , contains hyperbolic sines and cosines.
A quadrupole focusing in the yy plane would have these submatrices reversed accordingly.
Taking the limit as kL → 0 of the focusing submatrix, the diagonal elements become
cos(kL) → 1. The off-diagonal elements become

sin(kL)
k

→ L and −k sin(kL) → −k 2 L. In

this way a thin-lens focusing approximation is recovered, with focal strength f1 = k 2 L.
As stated previously, the form of the quadrupole transport matrix, given in Equation 3.1,
does not depend on whether the quadrupole is magnetic or electrostatic. The focusing
strength, k, however, is uniquely defined by the quadrupole type. If the quadrupole is
G , where G[T/m] is the focusing gradient
magnetic, the focusing strength k is given by k 2 = Bρ

and Bρ[T m] is the beam’s magnetic rigidity. In the case of an electrostatic quadrupole,
Vquad 1
k2 = Q
A T a2 where Vquad [V] is the voltage on the electrode, T [eV/u] is the kinetic energy,

and a[m] is the radial aperture.
If by varying the voltage of a quadrupole a waist can be established at the point of a
diagnostic device downstream, it is possible to analyze the resulting data to make an indirect
measurement of the beam’s emittance. Although there are several methods of analysis, the
one highlighted here uses Singular Value Decomposition (SVD) to compute a pseudoinverse
of a known matrix to solve the least-squares problem presented.
The σ matrix for one degree of freedom, x or y, can be written in the following form.






x2

 σ11 σ12  
=

≡
σ
σ21 σ22
xx

38


xx
x

2




(3.2)

Recall that the σ matrix in this propagates from point 1 to point 2 in a beam line as
σ2 = M σ1 M T . Multiplying out these M and σ matrices and solving for σ11 at point 2 in
terms of the σ elements at point 1 yields Equation 3.3. The fact that σ12 = σ21 has been
used for simplification.

2 (σ ) + 2M M (σ ) + M 2 (σ )
(σ11 )2 = x2 = M11
11 1
11 12 12 1
12 22 1

(3.3)

In this equation, the left-hand side, x2 , is found by squaring the RMS measurement of
the beamsize (

x2 ). The M transport matrix elements can be calculated in order to solve

for the elements of σ1 which in turn determine the emittance.
By performing many measurements, a system of equations of the form shown in Equation 3.4 can be built. If this is written as (σ11 )2 = A(M )(σcolumn )1 , then A(M ) is an n × 3
real matrix, with n the number of measurements. This matrix is not suitable for inversion,
as it is not square. However, a factorization of A, the SVD, may be found which enables the
calculation of the pseudoinverse of A. The pseudoinverse, also known as the Moore-Penrose
pseudoinverse, is a generalized inverse matrix which can be used to find the least squares
solution for an over-constrained system of equations where a unique solution may not exist.








x2


1

..
.
x2 n









=






2

2 )
(M11
1

2(M11 )1 (M12 )1

2 )
(M12
1

..
.

..
.

..
.

2 )
2
(M11
n 2(M11 )n (M12 )n (M12 )n

39












σ11 


σ12 


σ22

1

(3.4)

The SVD factorization is given in Equation 3.5. U is an n × n unitary matrix (U ∗ U =
U U ∗ = In ), Σ is an n × 3 diagonal matrix of non-negative reals, and V ∗ is a 3 × 3 unitary
matrix, the conjugate transpose of V . The matrices U and V act as rotations, while Σ is like
a scaling transformation. These matrices are used to find A+ , given in Equation 3.6.

A = U ΣV ∗

(3.5)

A+ = V Σ+ U ∗

(3.6)

In Equation 3.6, Σ+ is the pseudoinverse of Σ, which can be found by replacing the
non-zero entries along the diagonal by their reciprocals, and taking the transpose. Once
A+ is solved for, it may be multiplied with the column vector of measured beamsizes to
determine the σ matrix elements column vector, (σcolumn )1 = A+ (σ11 )1 . Having solved for
σcolumn , the σ matrix is simply the symmectric 2x2 matrix formed from the elements of the
σcolumn vector through Equation 3.2, noting again that σ12 = σ21 . It is then simple to find
the emittance by taking the determinant of the σ matrix and solving through =

det(σ).

This is the method of analysis chosen to determine the emittance from quadrupole scan
data taken before the RFQ at ReA. One example scan is shown in Figure 3.3a. There may
have been errors in measuring the beamsize from the viewer plate, or instabilities in the beam
itself which cause the measured data to not be entirely parabolic. When using a parabolic
fit to the data, the resulting emittance was found to strongly depend on the choice of the
subset of data analyzed, whereas the SVD method gave very consistent emittance values
across the various subsets of the data. The results of the SVD analysis is summarized in
Table 3.1. Note that Equation 2.18 has been used to define , β, α, and γ in this table.
The final RMS-emittances as seen in Table 3.1 were found to be x ≈ 1.3πmm · mrad and
40

(a)

(b)

Figure 3.3: Left (a): Quad Scan Measurements performed before RFQ at ReA. Solid lines
indicate measurement data, while dashed lines indicate SVD analysis. Right (b): Analyzed
viewer data for L057QB=−5.5 kV with overlay of beam position and RMS beam size in green.
Beam profiles in red are shown above and to the left of the beam spot.
Table 3.1: Analysis Results of Quad Scan Measurement Data (2-24-11)
L057QB
xRMS [πmm · mrad]
βx [m/rad]
−αx
−αx
γx [rad/m]
y

xx

scan
1.31
2.43 1.74
1.74 1.66

L057QB
y RMS [πmm · mrad]
βy [m/rad]
−αy
−αy
γy [rad/m]

yy

scan
1.09
1.73 -0.67
-0.67 0.87

≈ 1.1πmm · mrad, which are consistent with predicted values. The computed values for

βx,y , αx,y , and γx,y , vary depending on the chosen fitting method, while the emittance stays
fairly constant as it is the determinant of the matrix. It is important to note that the x and
y Courant-Snyder parameters are at least of the same order of magnitude. The signs of αx
and αy are opposite from one another, indicating that the beam is defocusing in the x-plane
and focusing in the y-plane, consistent with transport through a drift following a quadrupole,
as with this Quad Scan measurement.
During this measurement, a He+ beam from the LB commissioning source was tuned
through the center of the quadrupole L057QB, the second quadrupole in a quadrupole doublet
that comes before the RFQ. After the centering is established, the quad’s voltage is varied

41

from 5.5 kV to −2.5 kV in steps of 0.5 kV. The beamsize was measured on a CaF scintillating
viewer plate located in a diagnostic box immediately after the quadrupole. One of the
analyzed viewer images can be seen in Figure 3.3b. The transfer matrices used to complete
the SVD analysis were computed in COSY, with fringe-fields on. The inclusion of fringe-fields
had less than a 1% effect on the end result.

3.2

Solenoids

Solenoid magnets are one of the focusing elements used in ReA and FRIB that need to be
modeled accurately. There are multiple models for focusing elements, including thick and
thin lens approximations, as well as particle tracking through the magnetic field. COSY has
several different solenoid models, including a tanh-based on-axis field model with optional
fringe fields and a thick element model with a theoretical on-axis field distribution. XAL and
many other matrix modeling codes use an equivalent hard-edge model of the solenoid.
This so-called equivalent hard-edge model method utilizes an effective solenoid length
and field to approximate the magnetic field profile experienced by the beam, resulting in
“hard” edges to the solenoid instead of the smooth increase of real fields. An example of this
hard-edge equivalent model can be seen in Figure 3.4.
The motion in a solenoid can be seen as a series of three consecutive transformations
on the beam, and motivates the method for computing the full transport matrix for a
hard-edge solenoid model. First, the beam experiences fringe fields upon its entrance into the
solenoid. Next, the beam travels through the body of the solenoid where it interacts with the
longitudinal component of the field, Bz . Finally, the beam exits the solenoid and encounters
fringe fields once again.

42

Figure 3.4: Equivalent hard-edge model for LEBT 0.6 T solenoid located before RFQ. The
measured longitudinal field profile is shown with Bmax =0.6926 T, while the equivalent
hard-edge model consists of a drift section (B =0 T), a constant field section (B =0.6010 T),
followed by a second drift section (B =0 T).
The entrance fringe fields of a solenoid produce an azimuthal velocity in the beam.
Similarly, the exit fringe fields produce an azimuthal velocity in the opposite direction,
essentially canceling the effect of the entrance fringe fields by the time the beam has exited
the solenoid.
This can be seen by examining the forces acting on the beam as it moves through a region
of changing longitudinal magnetic field, where dB
dz > 0. Beginning from the Lorentz force,
dp
dt

= qv × B, when Bθ = 0 the azimuthal component becomes
dpθ
= qvz Br .
dt

(3.7)

The radial magnetic field, Br , is related to the change in the longitudinal magnetic field,
dBz
dz ,

through Amp`ere’s Law. From

B · dl = µ0 Ienc in a current free region, along a circular

43

path of radius r within the solenoid fringe field region, we arrive at

B · ds = Bz πr2 + 2πr

Br dz = 0.

(3.8)

Solving for Br in terms of Bz yields

r dBz
r
.
Br dz = − Bz ⇒ Br = −
2
2 dz

(3.9)

As the beam is moving along the z-direction, we can rewrite the expression in Equation 3.7
as
dpθ dz
= qvz Br .
dz dt

(3.10)
dp

dBz
Plugging in the solution for Br and simplifying yields dzθ = − qr
2 dz , which becomes

∆pθ = − qr
2 ∆Bz . Here ∆Bz , defined across the “hard-edge” boundary, is simply ∆Bz = Bz .
Looking at the x-component, ∆px = −∆pθ sin θ, and using y = r sin θ we find

∆px =

qr
qBz
Bz sin θ =
y.
2
2

(3.11)

x
Rewriting this in terms of the angle ∆x = ∆p
pz , and using the definition of magnetic

rigidity, Bρ = pqz , we arrive at the equation

∆x =

qBz
Bz
y=
y.
2pz
2Bρ

(3.12)

qBz
Bz
x=−
x.
2pz
2Bρ

(3.13)

Similarly, it can be shown that

∆y = −

44

This derivation leads to the matrix for a step into(+), or out of(-), a constant longitudinal
magnetic field, as




 1


 0

Medge± = 

 0


∓k

0

0

0 


1 ±k 0 



0 1 0 


0 0 1

where k is given by
k=

Bz
2Bρ

(3.14)

and Bz is the longitudinal magnetic field of the solenoid. The transport matrices used when
describing motion through a solenoid are 4 × 4 matrices, as the xx and yy planes are coupled.
z
Note that when stepping into a magnetic field, the field is increasing such that dB
dz > 0.
z
While stepping out of a magnetic field, the field is decreasing such that dB
dz < 0 and this is

what changes the sign that appears in front of the two k terms in the Medge matrix.
A similar derivation for the rotation of a charged particle in a constant magnetic field is used
for describing the motion of the beam through the body of a solenoid. The transformation of
the (x, x , y, y ) vector through the body of a solenoid with a constant magnetic field B = Bz zˆ
is given by the following matrix:


sin(θ)
2k

0
 1


 0
cos(θ)
0

Mbody = 

 0 −(1−cos(θ))
1
2k


0
− sin(θ)
0
where θ = 2kL.

45

−(1−cos(θ))
2k

sin(θ)
sin(θ)
2k

cos(θ)













(3.15)

When the entrance, body, and exit matrices are multiplied together, they give the final
matrix of the hard-edge equivalent solenoid model.

Msol = Mexit × Mbody × Mentrance

(3.16)

where

Mentrance = Medge+

(3.17)

Mexit = Medge−

(3.18)

resulting in


cos(kL)2

sin(kL) cos(kL)
k

sin(kL)2
k



sin(kL) cos(kL)






2
2
 −k sin(kL) cos(kL)
cos(kL)
−k sin(kL)
sin(kL) cos(kL) 



Msol = 
2


−
sin(kL)
sin(kL)
cos(kL)
2
 − sin(kL) cos(kL)

cos(kL)
k
k




k sin(kL)2
− sin(kL) cos(kL) −k sin(kL) cos(kL)
cos(kL)2
the full hard-edge equivalent solenoid transport matrix.
As noted, the hard-edge model of a solenoid simplifies the smooth growth and decline of
the magnetic field on axis to a rectangular function which consists of a step up, constant
field, and step down approximation. The matrix representation of this model is made up of
matrices representing the entrance into the solenoid, the solenoid body, and the exit from the
solenoid. The correct effective length and field strength, however, must be chosen for the
hard edge solenoid model to be valid.

46

3.2.1

Hard-Edge Model Parameterization

It is important to determine appropriate solenoid parameters when using the hard-edge
model. These parameters include the effective length, Le , and the effective field strength,
which may be represented by the magnitude of the scaling coefficient between current and
magnetic field, Be = µe [T/A] × I[A]. In order to determine which combinations of these
two factors would be appropriate in accurately modeling a given solenoid field, a numerical
method for comparing the hard-edge model with an approximated “real-field” model, which
closely follows the actual solenoid field profile, has been developed.
This numerical method computes the error between the transport matrix that steps
through the solenoid field in many small sections, approximating the “real-field” profile, and
the matrix where there is only one step through the solenoid, the hard-edge model. This
error is defined as the sum of the relative difference in the matrix elements, as given in
Equation 3.19. This simple difference definition provides a method of evaluation to determine
the degree to which the two matrices agree, allowing them to be compared numerically. Here
M [i, j] represents the i-th row and j-th column of matrix M .

∆Mtot =
i

j

Mstep [i, j] − MLe Be [i, j] 2
Mstep [i, j]

(3.19)

This calculation is performed for a wide range of combinations of Le and Be , mapping
out an effective solution space. The result is dependent on both the scaling of the magnetic
field, as well as the beam’s rigidity. A sample 2 × 2 plot of this error as a function of effective
field strength and effective length produced by this method is shown in Figure 3.5.
In this case, the magnetic field profile has been scaled down to 1 T. The beam has
magnetic rigidity Bρ = 0.5 T m, which corresponds to a H2 beam at 3 MeV/u. In Figure 3.5,
47

Figure 3.5: The ReA solenoid L090SN scaled to 1T for Bρ = 0.5 T m. The solution that
minimized Mtot is found to be at Le =0.24 m, Be =0.8 T.
red represents high error as calculated from Equation 3.19, while blue represents a region of
minimal error. The optimum solution found is Le = 0.24 m, with Be = 0.8 T. The coefficient
µe can be computed from Equation 3.20, and comes out to 0.0883 T/A.

Be [T] =

Bmax
Be
× I [A] = µe × I [A]
×
I
Bscaled
profile

(3.20)

The difficulty with using this method to accurately calculate the correct effective length
and field strength for the hard-edge model is that it depends on the magnitude of the magnetic
field in the solenoid, as well as the magnetic rigidity of the beam. Unfortunately, this brute
force method of calculating the relevant solenoid parameters is not practical when many
different types of beams with various Q/A ratios and energies will be used at ReA. It can,
however, prove useful as a cross-check in determining which Le and Be combinations are
reasonable, as calculated by a more general solution which would be valid for all cases.
One natural method would be to calculate Le and Be directly from the solenoid field

48

distribution [42]. The effective length is defined intuitively, as given in Equation 3.21. The
B ), motivates
importance of the focal length depending on B 2 , not B ( f1 = k 2 L where k = 2Bρ

the determination of the effective field strength in Equation 3.21.

Le =

Bz dz
Be

Be2 =

Bz2 dz
Le

(3.21)

Multiplying through by the denominators, the integral quantities may be defined as given
in Equations 3.22 and 3.23.

I1 = Be Le =

Bz dz

(3.22)

I2 = Be2 Le =

Bz2 dz

(3.23)

I2
Le = 1 .
I2

(3.24)

This system of equations can be solved as

I
Be = 2
I1

The end result of this parameterization is given by Equation 3.25.

Le =

B dz
B 2 dz

2

µe =

B 2 dz
I B dz

(3.25)

This method has the advantage that is does not depend on the central field of the solenoid,
B0 , or the beam’s rigidity, Bρ. This method also produced reasonable beam responses through
the ReA linac when modeling the beam with the machine settings from the QWR steering
experiment on August 29th, 2013, as seen in Figure 3.6. The H2 beam is at 0.6 MeV/u, with
no cavities turned on, and the original initial conditions given by the RFQ output simulation

49

have been used. The vertical dashed lines indicate the location of the diagnostic boxes. The
beamsize was fairly equal at both the L092 and L110 locations, represented by the vertical
lines at ∼6 m and ∼11 m. This is reproduced well with the solenoid parameterization given
by this method.

√
√
Figure 3.6: XAL simulation of transverse beamsizes,
x2 and
y 2 , σ11x and σ11y
as defined in Equation 3.2, for the ReA Linac without the third cryomodule using the solenoid
parameterization defined in Equation 3.25. Vertical dashed lines indicate diagnostic locations.
Note that the beam is focused in both x and y at the solenoid locations, and the beam sizes
at the third and fourth diagnostic locations, at ∼6 m and ∼11 m respectively, are of similar
magnitude, which was observed experimentally.

It is important to recognize the regions where this method is valid. Because this method
does not give any information about the beam’s rigidity or the scaling of the solenoid, we
should be sure to check that we are in the region where the equivalent method is comparable
to the step by step transport matrix. If the fields are too large for the beam’s rigidity, as
shown in Figure 3.7, then the beam’s motion will be influenced by the steeply increasing
fields at edges of the solenoid. The linear transport matrix approximation is only valid where
these nonlinear effects can be neglected, so this should be taken into consideration.
50

Figure 3.7: Numerical method results of Be vs. Le shown calculated at various Bρ and
B0 values that cover the relevant ReA ranges. The red regions indicate areas with high
divergence between the hard-edge and field-profile models, while the blue regions indicate
areas where the two models agree. Note that solutions are present (blue regions exist) for low
beam rigidities, Bρ, only at low solenoid fields, B0 . The solution range increases to higher
fields at higher rigidities.
In Figure 3.7, Bρ ranges from 0.22 Tm to 1.00 Tm, while Be ranges from 1 T to 8 T. The
rows from top to bottom represent Bρ = 0.22 Tm, corresponding to H2 at 0.6 MeV/u, Bρ =
0.32 Tm, Bρ = 0.44 Tm corresponding to He at 0.6 MeV/u, Bρ = 0.50 Tm, corresponding to
H2 at 3 MeV/u, Bρ = 0.75 Tm, and Bρ = 1.00 Tm corresponding to He at 3 MeV/u. The
columns represent Be = 1 − 8 T in steps of 1 T. This covers the entire range of relevant beam
rigidities and solenoid fields available at ReA.
The blue regions in Figure 3.7 are where the hard-edge model converges to the field-profile
model. For beams of low rigidity (upper rows), this is only valid at the lowest field strengths
(left most columns). As the beam’s rigidity increases, the model is valid at higher and higher
fields. This is because higher energy beams are more difficult to influence, and therefore
the hard-edge model is still valid even at the stronger fields. It is important to ensure that

51

when the hard-edge model is utilized, the relevant parameterization corresponds with an
appropriate region as given by Figure 3.7. This is typically not an issue when modeling
machine settings because during operation the solenoid’s fields are optimally kept at the
lowest possible values, which means we will most likely always be in an acceptable region,
shown in blue in Figure 3.7.

3.2.2

Solenoid with Corrector

Each solenoid at ReA contains horizontal and vertical correctors which can steer a misaligned
beam back toward its design trajectory. In this section we will restrict our discussion to
horizontal correctors only (By ), as the treatment is analogous for vertical correctors (Bx ). In
order to model this configuration simply, we insert a thin horizontal corrector dipole element
into the center of the body of the hard edge solenoid model. This total “thin-lens” matrix
model will then consist of three matrices multiplied together: the second half of the solenoid
body, a thin-lens corrector, and the first half of the solenoid body, as seen in Equation 3.26.

L
L
Mthin-corr/sol = Mbody ( sol ) × Mthin-corr × Mbody ( sol )
2
2

(3.26)

Here we must deal with 6 × 6 matrices, as the thin dipole element couples the longitudinal
motion with the transverse. The M [1, 4; 1, 4] sub-matrix of Mbody is defined in Equation 3.15,
L

plugging in L = sol
2 . The M [5, 6; 5, 6] sub-matrix of Mbody is simply a drift matrix with
L
L = sol
2 , and all other elements are zero.

The terms in the matrix for a thin corrector are simply derived from the equations of
motion for a charged particle with longitudinal velocity z˙ in the corrector’s magnetic field:
By for a horizontal corrector, Bx for a vertical corrector. In the case of a horizontal corrector,
52

from the Lorentz force we have m¨
x = −qvBy = −q zB
˙ y . Simplifying with Bρ = pq = mqz˙
results in x¨ =

−By z˙ 2
Bρ .

d (z z)
Using dt
˙ = z˙ z˙ + z z¨ and recognizing that there is no acceleration
−B

−B

d (z z).
˙ Integrating both sides produces x˙ = Bρy z z.
(¨
z = 0), the equation becomes x¨ = Bρy dt
˙

Therefore the only non-zero off-diagonal element in the thin horizontal corrector matrix is
element [2; 6], how x˙ depends on z.
˙
Substituting z = L, the length of the corrector, in the above equation, Mthin-corr for a
thin horizontal corrector with field By becomes simply the identity matrix with element
By L
x L instead.
M [2; 6] = − Bρ
. A vertical corrector with field Bx has element M [4; 6] = BBρ

Multiplying the matrices out results in the formula for the “thin-lens” transfer matrix
z L . B is the
as given in Equation 3.27. Here we define the dimensionless quantity ξ = BBρ
z

magnetic field of the solenoid, and By is the magnetic field of the corrector.










Mthin-corr/sol = 










1

Bρ sin ξ
Bz

0

cos ξ

0 −

Bρ(1−cos ξ)
Bz

0

Bρ(1−cos ξ)
Bz

0

0

sin ξ

0

1

Bρ sin ξ
Bz

0

ξ
By L sin 2
−
2Bz
ξ
By L cos 2
−
2Bρ
ξ
By L 1−cos 2
2Bz
ξ
By L sin 2
2Bρ

0

− sin ξ

0

cos ξ

0

0

0

0

0

1

L

0

0

0

0

0

1






















(3.27)

The “full-field” matrix solution is obtained by considering the total combined fields of
the solenoid and vertical dipole corrector, B = By yˆ + Bz zˆ, across the entire length of the
element, as opposed to breaking up the interactions as in the “thin-lens” approximation. The
Lorentz force, q v × B = vy Bz − vz By xˆ + (−vx Bz ) yˆ + vx By zˆ, produces the resulting
equations of motion, given in Equations 3.28-3.30.

53

mγ x¨ = qBz y˙ − qBy z˙

(3.28)

mγ y¨ =

−qBz x˙

(3.29)

mγ z¨ =

qBy x˙

(3.30)

Solving these equations of motion yields the “full-field” matrix solution as shown in
Equation 3.31. Here the quantities Bα =












Mfull-field = 











1

Bρ sin η
Bα

0

0

cos η

0

0 −
0
0
0

BρBz (1−cos η)
2
Bα
Bz sin η
Bα

BρBy (1−cos η)
2
Bα
By sin η
Bα

1
0
0
0

α L are defined.
By2 + Bz2 , and η = BBρ

BρBz (1−cos η)
2
Bα
Bz sin η
Bα
Bρ sin η
Bz2 L−
2
Bα
L−
2
Bα
2
B (1−cos η)
1− z 2
Bα
Bρ sin η
By Bz L− B
α
2
Bα
By Bz (1−cos η)
2
Bα

0
0
0
0
1
0

BρBy (1−cos η)
−
2
Bα
B sin η
− yB
α
Bρ sin η
By Bz L−
2
Bα
2
Bα
By Bz (1−cos η)
2
Bα
Bρ sin η
By2 L−
2
Bα
L−
2
Bα
2
B (1−cos η)
1− y 2
Bα
























(3.31)

These two solutions, the “full-field” and the “thin-lens”, agree when we make the assumption that the corrector field is much weaker than the solenoid field, By

Bz , such

that By /Bz → 0. If we can confidently determine that this assumption is valid in the case
of Hydrogen beams, H2 , which require the weakest solenoid fields, then it will therefore be
valid for all beams.
The two matrices given in Equations 3.27 and 3.31 are applied to a σ matrix of parameters
typical of a beam exiting the ReA RFQ. The ratio of By to Bz is varied across both models,
and the beam sizes at the end of the solenoid models are then compared.
An approximately 1% effect on the beam size is found when By is approximately 25% the

54

Beamsize Percent Differences at Solenoid Exit
15
x
y
z

Percent Difference

10

5

0

−5

−10

−15
0

0.5

1

1.5

By/Bz

Figure 3.8: Percent difference between the “full-field” and “thin-lens” models for beamsizes
calculated at the solenoid exit, z = 0.2 m, as a function of relative corrector field strength,
By , to the main solenoid field, Bz = 2 T.
strength of Bz . In Figure 3.8, the percent difference between the two models, the “thin-lens”
and “full-field” solutions, of the beamsize at the solenoid exit is plotted as a function of
corrector field strength, By where the solenoid field has been set at Bz = 2 T. Typical
solenoid fields for H2 beams in the ReA linac are on the order of a few Tesla, while the
corrector fields only reach a few milliTesla. This is less than 0.1% of typical Bz values, and in
the region in Figure 3.8 where the percent difference is

1%. Therefore the error introduced

by the simplified model is completely negligible in this region, as we hoped. It is concluded
that the thin-lens model of a corrector at the center of the solenoid is therefore valid.

55

3.2.3

Solenoid Error Study

As with any kind of modeling, it is always important to remember that the physical elements
being modeled are subject to real-life conditions, with minor imperfections such as asymmetries
in the field profiles, misalignments of the elements, or fluctuating temperatures potentially
affecting their performance. It is therefore necessary to study what effect small errors will
introduce to the system being modeled. In particular, it is extremely important to understand
how small errors affect a system that contains many of the same element, as in these cases
small systematic errors can compound quickly and cause serious problems.
For example, FRIB uses a total of 87 superconducting solenoids for beam focusing through
its linac and folding segments [29]. It is therefore important to know the effect of errors in
the solenoid strength on the beamsize exiting the solenoid. To illustrate this technique, we
introduce small errors, di,k (on order of 1.0% and 0.5%), to the design solenoid field strengths
in FRIB Linac Segment 1, and calculate the resulting change in beamsize, ∆σ. The new
solenoid strengths are produced through the formula B1,k = B0,k × 1 + di,k , where k is
the solenoid index, and i is the run number through the XAL model. The error introduced,
di,k , is a Gaussian error di,k = µ + dRMS × randn; where µ = 0, dRMS is the maximum error
amplitude, and randn is a pseudorandom value generated by MATLAB drawn from the
standard normal distribution.
With the error introduced, there are fluctuations in the beam envelope about the original
case, as seen in Figure 3.9. The RMS of the relative difference in σ in each case approaches
the value of the error that was introduced to the solenoid fields. A 1% random error on the
solenoid fields produced ∆σσRMS ≈ 0.1 while A 0.5% error resulted in ∆σσrms ≈ 0.05. From the
plots of

∆σRMS
σ

on the right-hand side of Figure 3.9, one can see that the error has essentially

56

0.15

20

x ,y

0 0

x
y
0.1

RMS of ∆σ

σx [mm]

x,y

/σ

15

10

0.05

0

5
0

20

40

60

80

0

100

20

40

60

80

100

60

80

100

S [m]

S [m]

(a)
20

0.15
x
y

18
0 0

x ,y

/σ

12

0.1

x,y

14

RMS of ∆σ

σx [mm]

16

10
8

0.05

6
4

0
0

20

40

60

80

100

0

S [m]

20

40

S [m]

(b)
Figure 3.9: FRIB Linac Segment 1 (XAL) with random Gaussian error added to each solenoid.
Gaussian errors of about (a)0.5% and (b)1% have been used — note how the RMS of the
relative difference in σx,y grows by approximately the percentage of the Gaussian error
introduced to the system. It is an important result that this growth levels out without
increasing exponentially.
saturated after ∼10 m, and does not increase with the additional segments of solenoids.
This result is encouraging in that random errors on the solenoid fields produce an
acceptable error in the beam size that does not increase rapidly throughout the beam line.
These errors should be quite low as superconducting magnets can be kept stable at a precise
field. Even with small errors throughout the beam line, only a small change in the beamsize
is observed, not a compounding effect throughout the linac which would be problematic for
operations.

57

3.3

Accelerating Cavities

Accelerating cavities are arguably the most important beam line elements to be modeled
in a linac. ReA, as well as FRIB, uses Superconducting Radiofrequency (SRF) cavities to
accelerate the beam. This section focuses on the Quarter-Wave Resonator (QWR) type, seen
in Figure 3.10, as this is the type used during the experiment detailed in Chapter 4.

Figure 3.10: Cut-away image of the ReA β = 0.041 QWR cavity type. The two accelerating
gaps that the beam travels through towards the bottom of the cavity are shown in detail.

This section discusses the theoretical background of particle acceleration through QWRs
in order to understand how these structures are modeled. It is important to first be familiar
with the electromagnetic fields present in the cavity. These fields ultimately determine the
motion of particles through these accelerating structures, so they shall be discussed in more
detail here.

58

3.3.1

Electromagnetic Fields

In a vacuum, the wave equation is given by equation 3.32 where ∇2 is the Laplacian and
c is the speed of light. This equation is derived from Maxwell’s equations, and determines
the propagation of electromagnetic waves. The electric field, E, or magnetic field, B, can be
substituted for the wave function, u(r, t).
1 ∂2
∇2 − 2 2
c ∂t

u(r, t) = 0

(3.32)

By separating variables, u(r, t) = u(r)eiωt , a time-independent spatial equation for steady
state oscillations is achieved, known as the Helmholtz equation, given in 3.33. In this equation
the wave vector, k, is given by k = ω/c.

∇2 + k 2 u(r) = 0

(3.33)

In a hollow resonator, such as the case with accelerating cavities, the electric field satisfies
the boundary condition that the electric field vanishes at the surface [43, p. 310]. Using the
boundary conditions and the Helmholtz equation, a solution for the electromagnetic fields
inside the cavity can be found.
Due to the complex surface structures, analytical solutions are foregone in favor of
determining the solution numerically on a mesh. Dedicated commercial electromagnetic field
solvers are available for this purpose, such as CST Microwave Studio [44].
The fields present in the ReA β = 0.041 cavity are shown in Figure 3.11. These field
profiles are calculated along the z-axis of the 3D field maps [45]. The z-axis is centered on
the beam port of the cavity, and taken as the ideal trajectory of the beam. Figure 3.11 plots

59

the electric and magnetic fields vs. z, the longitudinal position in the cavity, with z = 0
positioned at the midpoint of the cavity.

ReA β=0.041 Cavity

15

Ez

10 × Ey

10

cBx

Field (MV/m)

5

0

−5

−10

−15
−100

−50

0

z (mm)

50

100

Figure 3.11: ReA β = 0.041 Cavity Fields, including the accelerating electric field, Ez , and
the transverse fields Bx and Ey . Note that Ey is very small, and so therefore has been
multiplied by ten for the purposes of this graph.

Ez is the electric field in the longitudinal direction which is responsible for the acceleration
of the beam through the cavity gaps. Ey is the electric field in the vertical direction, and
can be positive or negative depending upon the design of the cavity. From the symmetry of
ˆ which results in B = Bx (z)ˆ
the cavity, B = Bφ (r)φ,
x along the z-axis. For a right-handed
coordinate system with the z-axis to the right and y to the top of the page, the x-axis points
into the page. This corresponds to Bx pointing out of the page for negative z, and into the
page for positive z, as is seen by the sign of Bx going from negative to positive along the
z-axis in Figure 3.11.
In this instance, the second gap, located at positive z, is accelerating the beam. The
60

fields oscillate at the RF frequency of these cavities. From this we can see that at an earlier
time, the first gap, located at negative z, would have the positive accelerating Ez field.
The longitudinal electric field is asymmetric with respect to z, as is the magnetic field, Bx ,
while the shape of Ey is symmetric with respect to z. It is also important to note that the
magnitude of Ey has been multiplied by 10. For QWR fields, Ez

3.3.2

cBx

Ey .

Energy Gain

As previously stated, the longitudinal electric field is responsible for the acceleration of the
beam. The energy gain through a cavity is given by the following formula:

∆W = qEa L · T T F (β) cos φ,

(3.34)

where
Ea =

1 L/2
Ez (z) dz.
L −L/2

(3.35)

In Equation 3.34, q is the charge, L and φ are the cavity length and phase respectively.
Ea is the average accelerating electric field, given by equation 3.35 and T T F (β) is the transit
time factor. The transit time factor gives the relative acceleration for beams traveling at
different velocities, as it is a function of β.

3.3.2.1

TTF Calculation

The full definition of the transit time factor is given by equation 3.36.

T T F (β) =

L/2
E (z) cos (ωt(z)) dz
−L/2 z
L/2
E (z) dz
−L/2 z

61

(3.36)

Assuming that the change in velocity across the gap is small compared with the beam’s
incoming velocity, ∆v
v

1, the formula for the transit time factor can be simplified. This

assumption is valid for velocity gains ∆v
v

10%, which cover all cavities in ReA and FRIB.

The lowest velocity beams experience an energy gain increase across one accelerating gap
of ∼0.03 MeV/u out of ∼0.6 MeV/u. This 5% energy gain corresponds to an ∼2% velocity
gain, which allows us to elimate v, and thereby t, as a function of z in Equation 3.36.
ωz
Neglecting the change in velocity across the gap, we can rewrite ωt(z) = ωz
v = βc [46].

An analytical solution can then be calculated, which results in Equation 3.37.

TEz (β) = sin

πdEz
βλ

sin

πgEz
βλ
πgEz
βλ

= sin

πdEz
βλ

Tg (β, gEz )

(3.37)

This form is useful for calculations, and is valid in most practical cases. The gap-to-gap
distance, dEz , and the gap length, gEz , are defined by Equations 3.38 and 3.39 respectively.
Here α represents each index, while Fα stands in for each field, Ez , Ey , and Bx , respectively.

dα = 2

gα = 4

L/2
0

L/2
zFα (z) dz
0
L/2
Fα (z) dz
0

z − d2α Fα (z) dz
L/2
Fα (z) dz
0

(3.38)

(3.39)

Figure 3.12 shows calculated T T F (β) curves for ReA β = 0.041 and β = 0.085 QWRs.
At each curve’s optimal point, the beam has the exact velocity required for the energy gained
from the first gap to result in the beam passing between the two gaps in the time that it
takes the RF to switch 180◦ . The distance between the gaps is designed to be ∼βopt λ/2.
Faster than optimal velocity beams travel between the two gaps of the cavities in less time

62

Figure 3.12: TTF of β = 0.041 and β = 0.085 cavities
than it takes RF to fully oscillate 180◦ , and therefore experience slightly less acceleration
than an optimal velocity beam. In a similar manner, slower than optimal velocity beams
take a longer time between the two gaps, resulting in the RF oscillating more than 180◦ and
receiving less acceleration.
Following the curve down from the optimal point to slower and slower velocities to the left,
there even comes a point where the curve equals zero. This is at the point where the beam’s
velocity is such that the transit time between the two gaps results in equal acceleration in
the first gap and deceleration in the second gap, for zero total acceleration. The cavities are
optimally operated well away from this point.

63

3.3.2.2

Cavity Geometry

The cavity types that will be utilized at FRIB are called Quarter-Wave Resonators (QWRs)
and Half-Wave Resonators (HWRs), as seen in Figure 3.13. As discussed previously, each
cavity type is optimized for a certain velocity, at which point the beam experiences the
maximum kick at both gaps. This is a function of the distance between the two gaps and the
frequency of the RF.

Figure 3.13: QWR and HWR Cavities [29]. The QWR cavities operate at 80.5 MHz, while
the HWR cavities operate at 322 MHz.

The two types of QWR cavities in ReA have been optimized for β = 0.041 and β = 0.085.
The β = 0.085 cavities operate at the same frequency, 80.5 MHz, as the β = 0.041 cavities,
but the spacing between the gaps is larger, resulting in a higher optimal β. The half-wave
resonators (HWRs) of FRIB operate at a higher harmonic frequency, 322 MHz, and have
been optimized for β = 0.29 and β = 0.53.
The difference in gap spacing between the cavity types can clearly be seen in Figure 3.13.
64

The inner conductor, which comes down through the center of the cavity, is noticeably larger
in the β = 0.085 cavity from the β = 0.041 cavity, thereby optimizing the transit of a faster
beam. The inner conductor shields the beam from experiencing the cavity fields when the
RF is changing. The same is true for the HWR cavities.
If the beam were allowed to see the fields continuously, it would be accelerated and
decelerated the same amount, resulting in no net acceleration. The longitudinal electric field,
oscillating back and forth between pointing in the positive, and then negative, z-direction, is
shielded by the cavity geometry from the beam when it is pointing in the negative z-direction.
In this way, the beam can be accelerated in both gaps as it passes through the cavity.

3.3.3

Focusing Effects

As the beam passes through the cavity, recall that it experiences not only the accelerating
longitudinal electric field, Ez , but also transverse fields. Figure 3.14 depicts schematically
the electric field lines present in an RF gap. As the beam’s position deviates off the central
axis, the transverse components of the electric field are seen in addition to Ez .

Figure 3.14: Electric field lines in an RF gap.

65

The radial electric field, Er , has the opposite sign of Ez following from Maxwell’s equation
∇ · E = 0, or in cylindrical coordinates, 1r

∂(rEr )
∂r

z
+ ∂E
∂z = 0. This leads to the fact that when

the electric field is accelerating the beam, as in Figure 3.11 for positive z, Ez is decreasing
with respect to z toward the end of the gap, meaning Er must be positive. This is the case
of RF defocusing [46].
The transverse focusing, with the cavity modeled as a thin lens, can therefore be calculated
as having a strength kr as given in Equation 3.40.

kr =

¯
πqEz T T F (β)Lf
sin (−φ)
2
¯γ
cEr β¯

(3.40)

Note that this quantity is negative, indicating transverse defocusing. The longitudinal
focusing has the opposite sign from the transverse defocusing, as given in Equation 3.41.

kz = −2kr γ¯ 2

(3.41)

With the cavity changing the energy of the beam, these quantities have been calculated
at the midpoint of the cavity, using the beam’s energy at the midpoint of the cavity, W to
W , and then β¯ =
calculate γ¯ = 1 + E
r

1 − 12 . It is important that β¯ is not calculated as the
γ¯

average of βi and βf , as β is not linear and therefore this will not yield the correct result.
In addition to these focusing effects, quarter-wave resonator cavities are known to steer
the beam vertically due to fields present from the asymmetry of this cavity type with respect
to the beampipe. This effect is very important to study, as it has been previously calculated
but not measured. This will be discussed in more detail in Chapter 4, along with presentation
of measurement data taken at ReA investigating this effect.

66

3.3.4

Benchmark Considerations

Multi-particle tracking through simulated fields is one method of beam line modeling that is
used, while another is a thin-lens adaptation for matrix modeling. Multi-particle tracking
is crucial for studying higher-order effects, while matrix modeling is utilized in many cases
where speed of calculation is a factor, such as with on-line modeling, or studies with large
variable spaces to explore.
The energy gain from an RF cavity is the one of the first things to investigate when
comparing cavity models in multi-particle tracking codes with cavity models in envelope
matrix models. As previously discussed, the energy gain imparted to a beam from an RF
cavity is a function of the beam’s charge to mass ratio, and the cavity’s voltage, phase and
geometry, as given by Equation 3.34. The geometry of the cavity determines its transit time
factor, or T T F , which is a function of the cavity’s efficiency with respect to beam velocity.

3.3.4.1

TTF Computation

Many multi-particle tracking codes, such as IMPACT, look up T T F values in a table,
interpolating for points in between. Other codes, such as matrix-based XAL, fit a high-order
polynomial to the region of interest, and use this polynomial and its derivative to calculate
the T T F to use. Seemingly small differences between these two methods can result in large
differences in the overall energy gain calculated, so it is important to study this carefully in
order to make sure both codes are using the same values when benchmarking the codes.
This is particularly important when the T T F is changing rapidly, as is the case when
entering the third cryomodule at ReA. In Figure 3.15, note that where the transition between
the the β = 0.041 and β = 0.085 cavities takes place, at β ≈ 0.05, the β = 0.085 T T F curve
is much steeper than the curve for the β = 0.041 T T F , which is fairly flat.
67

Figure 3.15: ReA T T F for β = 0.041 and β = 0.085 cavities, with operating regimes
highlighted. The β = 0.041 cavities are located in the first two cryomodules, CK01 and CL01,
while the β = 0.085 cavities are located in the third cyromodule, CM01.
The small difference in energy gain as calculated by the two models for one cavity may
not be a large effect,

1%, but when there are multiple cavities one after another, such as in

a linac, these errors compound and can produce increasingly inaccurate results.

3.3.4.2

Thin-lens Modeling

One final aspect to consider when benchmarking the energy gain through an RF cavity
between the two types of modeling methods involves the intrinsic difference between the
modeling methods. Multi-particle tracking follows the beam smoothly through the changing
cavity fields, resulting in a smooth energy gain profile, while the matrix representation of a
cavity models the energy gain as a thin-lens, with a single kick at the center of the cavity, or

68

Figure 3.16: Modeling two gap acceleration in XAL as compared with IMPACT for a β = 0.041
accelerating cavity.
with two kicks located at the cavity’s physical gap locations. Figure 3.16 shows the difference
in energy gain profiles between the smooth multi-particle tracking model and the two-gap
thin-lens matrix model.
When modeling the second gap kick, the relative phase between the two gaps has to be
calculated. Since the envelope for the reference particle is being modeled, it should arrive
at the design phase (nominally −20◦ for ReA). This is accomplished by running the matrix
model twice, first calculating the energy gain and correct phases for the cavities, and finally
using these correct phases to track the beam through the lattice.
For simplification, the one-kick model at the midpoint of the cavity is often used instead
as it produces nearly identical results. Moreover, the one-kick model reflects the fact that
the thin-lens model is used simply as an approximation to the real-field. When speed of
calculation is crucial, the simplest reliable model is desired, and so the thin-lens model is
appropriate in these cases.

69

3.4

Lattice Examples

With the models established for the linac beam line elements, solenoid and RF cavities, the
full accelerator lattice can now be modeled for various beams and machine settings. It is also
important to benchmark the results of multi-particle tracking, which uses full field profiles in
its modeling, with that of the simplified matrix method approach which only include linear
terms, and is therefore a linear approximation.
After many benchmarking considerations have been handled, as discussed in Section 3.3.4,
a reliable model can finally be established for the linac in question. This section will present
models for both the ReA and FRIB linacs which have been benchmarked in XAL, a matrix
modeling code, against the design code, IMPACT, a multi-particle tracking code.

3.4.1

ReA Linac

The solenoid focusing results shown in Section 3.2 in Figure 3.6 was for the case of an RFQ
energy beam,0.6 MeV/u, without any cavities turned on. In the case where the cavities are
on, we also need to make sure that the energy gain is modeled correctly, as well as the beam
size. Using the models previously discussed for both the solenoids and RF cavities, the ReA
linac can be accurately modeled.
Figure 3.17 shows the benchmarking modeling result for a H2 beam, Q/A = 0.5. The
nominal settings for this beam are modeled with the first rebunching cavity on and the
β = 0.041 cavities set at 16 MV/m and φ = −20 ◦ accelerating phase. In XAL, the solenoid
hard-edge effective lengths and field strengths are calculated through the parameterization
discussed in Section 3.2.1 using the real measured field profiles. The β = 0.041 cavities are
modeled using the single-kick “thin-lens” model.

70

6

6
IMPACT
XAL

5

5

4

4

β [m/rad]

3

3

y

βx [m/rad]

IMPACT
XAL

2

2

1

1

0

0
0

1

2

3

4

5

6

0

1

2

3

S [m]

(a)

5

6

(b)

60

1.5

IMPACT
XAL

IMPACT
XAL

1.4

50

1.3
1.2

Ek [MeV/u]

40

βz [deg/%]

4

S [m]

30

20

1.1
1
0.9
0.8
0.7

10

0.6
0.5

0
0

1

2

3

4

5

0

6

1

2

3

4

5

6

S [m]

S [m]

(d)

(c)

Figure 3.17: Benchmarking results of the ReA linac through its second cryomodule show
good agreement between the matrix modeling code XAL, and the multi-particle tracking code
IMPACT. Plots show the results for transverse β functions (a) βx and (b) βy , longitudinal β
function (c) βz , and finally the kinetic energy per nucleon, (d) Ek .
The benchmarking result of modeling the ReA linac in Figure 3.17 shows XAL’s calculation
in red as compared with the design lattice model run from IMPACT in blue. Through thick
elements, such as the solenoids, IMPACT outputs results at the beginning and end of the
elements, while XAL tracks through the thick elements. This is the reason for the difference
between the display of the two curves. A good agreement is found for transverse and
longitudinal motion, as well as energy gain through the cavities. Slight discrepancies in the
longitudinal plane are expected due to the inherent differences in the modeling methods.

71

3.4.2

FRIB Linac Segment 1

The benchmark results for FRIB’s Linac Segment 1 between IMPACT and XAL using a
single kick at the center of the cavity, the “thin-lens” model, is shown in Figure 3.18. Note
that in this example, the phases for the matching cavities for the transition between the
β = 0.041 cavities and the β = 0.085 cavities had to be adjusted in order to account for the
difference in T T F values calculated, as discussed earlier in Section 3.3.4.

Figure 3.18: XAL Cavity Benchmark with IMPACT for a U 33+ beam through FRIB’s Linac
Segment 1 to the Stripper Foil, including modeling of β = 0.041 and β = 0.085 cavities.

FRIB’s Linac Segment 1 (LS1) consists of 16 β = 0.041 cavities, 96 β = 0.085 cavities, and
44 superconducting solenoids. The energy for the nominal case of a U 33+ beam through LS1
increases from 0.3 MeV/u to 16.6 MeV/u. FRIB’s Folding Segment 1 (FS1) is also modeled,
which brings the beam through a 180◦ bending section after a stripper foil which increases
the beam’s charge state from 33 − 34+ to 76 − 80+ for more efficient subsequent acceleration.

72

The nominal design for FRIB’s LS1 and FS1 actually transports both U 33+ and U 34+ to
satisfy the increased beam intensity requirements of modern experiments. The concept and
challenges associated with multiple charge state transport are presented in the next section.

3.5

Multiple Charge State Transport

There are many challenges associated with multiple charge state transport, but this method of
beam delivery is essential in order to attain the high currents required by modern experiments.
A brief summary of the history of multiple charge state transport is presented in this section,
followed by the methodology and benchmarking process of matrix modeling methods. Multiparticle tracking results of multiple charge state simulations in the FRIB Linac Segment 1
beam line (LS1) have been used as the basis for the benchmarking process. As stated before,
it is important to model multiple charge states with a matrix based model in the case of an
on-line model, where the speed and performance of the calculation becomes an issue. For
this purpose, the slower multi-particle tracking simulations are no longer a viable option.

3.5.1

Background

One of the primary performance criteria for FRIB’s driver linac beam is to deliver beams with
energy greater than 200 MeV/u with high beam currents of ∼400 eµA. During the design
of FRIB, modern ECR ion sources were able to produce enough current in a single charge
state to satisfy these requirements for all ion species up to Xe [29]. However, even the best
sources, such as VENUS at Lawrence Berkeley National Lab (LBL), were only able to deliver
∼200 eµA

of U 33+ or U 34+ [47] as of 2008. Therefore, ion beams heavier than Xe required

design of simultaneous transport of at least two charge states in order to achieve the power

73

requirements of FRIB. Improvements in beam development have resulted in state-of-the-art
ECR ion sources producing ∼400 eµA of U 33+ or U 34+ which meets the requirement for
single charge state transport [48], however multiple charge state transport is still an important
technique that may be used for increasing the delivered beam current on target.
Experimenters at Argonne National Lab (ANL) in 2007 set up ECR-LEBT beam line
sufficient for the study of transport and merging of multiple charge state Bi beams [49]. The
setup consisted of an ECR ion source, a 100 keV platform and a Low Energy Beam Transport
section, designed as an injector for a high-intensity heavy-ion linac, such as FRIB.
In 2008, they successfully demonstrated the separation and recombination of 209 Bi20+
and 209 Bi21+ [50]. Figure 3.19 shows the simultaneous transport of both charge states on the
left, as well as a superimposed image of the transport of each charge state individually. This
is an important proof-of-concept that two charge states can be transported simultaneously to
effectively double the beam intensity.

Figure 3.19: Results from ANL multiple charge state transport: pepper-pot image of the
combined bismuth beam with charge states 20+ and 21+ (on the left) and superimposed
images of 209 Bi20+ (blue) and 209 Bi21+ (red) beams transported individually (on the
right) [50].

Multiple charge beam dynamics have been studied further in acceleration through an
RFQ [51] and through an ion linac [52]. The verification of multiple charge state beams

74

from an ECR ion source doubling the intensity of a high quality heavy-ion beam is highly
encouraging, and motivates this technique as an option for FRIB to utilize for its heavy-ion
beams.

3.5.2

Methodology

The following is a description of a method that can be used to model multiple charge state
beams with a matrix envelope calculation instead of multi-particle tracking. The principle
lies in tracking individual charge states, followed by a step where these individual charge
states are combined in order to get the total extent of the beam. This approach is valid
as long as beam-beam effects, such as space charge or interaction between the individual
charge states can be neglected. This is appropriate for the cases we will consider, and, in
fact, multi-particle tracking also ignores these effects. The space charge potential or other
higher order terms can be accounted for, but we constrain this section’s discussion to the
methodology of transporting multiple charge states in regions where these effects can be
neglected.
To begin with, we track the reference charge state. In the case of the FRIB driver linac,
for U238 this would be the 33+ charge state. This step is performed like any other matrix
model simulation of a single charge state beam. We then track each additional charge state
with respect to this reference charge state orbit, as established in the first run.
In order to do this, we need to carefully consider elements that change the beam’s energy,
such as the accelerating cavities. In order to handle the case of modeling the cavities correctly
for multiple charge states, it is necessary to calculate ∆φ for each different charge state at each
cavity with respect to the reference charge state. As the energy gain (given in Equation 3.34)
depends on the charge, q, each additional charge state will have slightly different energies
75

than the reference charge state. This means that they will arrive at slightly different times,
and therefore slightly different phases, at each cavity. The equation for calculating the time
of flight for each charge state is given in Equation 3.42.

ti =

d
c 1 − 12

(3.42)

γi

Charge states gain different amounts of energy based on their different charge and also
based on their different arrival times to the cavity gaps which results in them experiencing
different cavity phases. As higher charge states gain more energy through the first gaps, they
arrive earlier at the subsequent gaps which means they experience a lower phase and therefore
gain less energy than the lower charge states. This continues until they have roughly the
same energy at which point the process repeats itself, with higher charge states gaining more
energy than the lower charge states. This results in a longitudinal simple harmonic motion of
the charge states about the nominal charge state.
We can then calculate the change in arrival time of the additional charge states from the
reference charge state, ∆ti = ti − t0 , which gives the change in RF phase that the charge
state will experience in the cavity, ∆φi = ωRF ∆t. Using this formula, we get the result for
the energy gain that the different charge states experience in Equation 3.43.

∆Ei =

Qi
eET L cos (φ + ∆φi )
A

(3.43)

This formula gives the total energy gain for each charge state, with the reference charge
state experiencing an energy gain of ∆E0 , calculated using the reference charge state Q0
and at the reference phase such that ∆φ0 = 0. We might expect at first glance that because

76

the higher charge states will get a larger energy gain due to their higher Q, that they would
receive more energy than the reference charge state overall, while the the opposite would
be true for the lower charge states. We can model this with FRIB Linac Segment 1 (LS1),
where the reference beam is U 33+ . Additional charge states are shown for ∆Q = ±1, ±2. For
practical purposes FRIB LS1 has only been designed to accelerate two charge states, Q = 33+
and 34+, but the goal here is to show the characteristics of the energy gain of multiple charge
states simultaneously, so we have included more charge states to illustrate this effect.
This result is shown in Figure 3.20 on the left, where ∆φ = 0. However, if the higher
charge states arrive earlier to the next cavity because of a larger energy gain from the
initial cavity, they consequently see a lower cavity phase (which would be rising to the ideal
accelerating phase, for example −20◦ ) and get a lower energy gain from the second cavity.

Figure 3.20: Multiple Charge State Energy Gain: ∆φ = 0 on left, ∆φ = ωRF ∆t on right. As
discussed in the text, modeling the separate charge states with ∆φ = 0 yields the unphysical
result where each charge state arrives at the accelerating phase of the reference charge
state, resulting in higher charge states receiving higher energy gains and lower charge states
receiving lower energy gains than the reference charge state. Instead, each individual charge
state’s arrival time at each cavity must be calculated in order to correctly model the multiple
charge state’s energy gain, as seen here on the right.

The opposite happens for lower charge states than the reference charge state. This bunches
the beams longitudinally, and keeps the energies of the additional charge states oscillating
77

about the reference charge state’s energy. The same stability condition that works for one
charge state, which is why accelerating cavities are operated at a phase off from the crest, as
discussed in Section 2.2, operates in the same way for multiple charge states. We see that
when we include ∆φ in the calculation in order to account for the difference in arrival times
of the different charge states, we get the correct result as seen in the right of Figure 3.20.
Another element that has to be treated carefully are bending dipoles, which change the
beam’s reference orbit. In order to calculate the separation of charge states through bend
elements, we need to account for the different bend radii of the different charge states.
It is clear that the different charge states will each have its own bend radius from the fact
p

that each charge state has a different rigidity, Bρi = qi . In order to calculate the position
i
of each charge state’s beamlet with respect to the reference charge state’s orbit, we must
shift the beam centroid of the addition charge state’s calculation off by the amount that is
calculated to occur in theory. This way when plotting the results, each charge state will have
the same reference orbit, instead of setting its own trajectory through the the dipole as zero.
As previously stated, charge states spread apart in dipoles because of their different
rigidities. In order to get the correct dispersion, each additional charge state is given an
offset to its average momentum in the bend axis coordinate plane which is calculated from
the difference in rigidities. We can accomplish this by modifying an element in XAL’s 7 × 7
sigma matrix, defined by Equation 3.44 [53].

78




 x2


 xx



 xy



σ =
 xy


 xz



 xz


x

xx

xy

xy

xz

xz

x

x2

xy

xy

xz

xz

x

xy

y2

yy

yz

yz

y

xy

yy

y2

yz

yz

y

xz

yz

yz

z2

zz

z

xz

yz

yz

zz

z2

z

x

y

y

z

z

1
























(3.44)

The seventh row and column in XAL’s sigma matrix refer to the respective centroids
of the beam which transform linearly with the other beam moments, x2 , xx , etc. The
beam sigma matrix, σ, given in Equation 3.44, transforms according to Equation 3.45, where
M is the transfer matrix. The transfer matrix in its general form is Equation 3.46.

σ1 = M σ 0 M T

(3.45)


 (x|x)


 x |x



 (y|x)



=
 y |x
M



 (z|x)



 z |x


( x |x)


x|x

(x|y)

x|y

(x|z)

x|z

x |x

x |y

x |y

x |z

x |z

y|x

(y|y)

y|y

(y|z)

y|z

y |x

y |y

y |y

y |z

y |z

z|x

(z|y)

z|y

(z|z)

z|z

z |x

z |y

z |y

z |z

z |z

x |x

( y |y)

y |y

( z |z)

z |z

(x| x ) 



x| x



(y| y ) 



 (3.46)
y| y



(z| z ) 



z| z



1

The notation in Equation 3.46 is adopted from [53]. For example, x|x

simply refers to

the relevant component of M which gives the dependence of x on x through the element

79

being modeled by matrix M . These terms depend on the fields present inside the element,
and can be calculated as shown earlier in this chapter.
If we consider the case of a horizontally bending dipole, where we expect the x-position
of the beam to spread out, we can look at the reference charge state’s average position in x,
x , as a function of the matrix elements for a dipole, and the vector representing the beam
centroid, the seventh column in the σ matrix. This equation is given in 3.47.

x 1 = ( x |x) x 0 +

x |x

x 0 + ( x |y) y 0 +

+ ( x |z) z 0 +

x |z

x |y

y 0

(3.47)

z 0 + ( x |x)

For our charge states, x 0 = x 0 = 0. The dipole transfer matrix also has ( x |y) =
x |y

= ( x |z) = ( x |x) = 0. This leads to the only nonzero element contributing from

the dipole to be ( x |z). The formula for this term for an additional charge state is given
10p[GeV/c]

in Equation 3.48. Here the definition of magnetic rigidity, Bρ [T · m] = 2.9979Q

[34], has

been used to solve for ρ. In Equation 3.48 all quantities are given for the additional charge
state, such that θ = Lρ = eBL
pc . The only contributing element from the σ matrix is then z 0 .
This formula is given in Equation 3.49. Here, “ref” indices indicate the reference charge state,
while all other quantities are given for the additional charge state being calculated.

( x |z) = ρ (1 − cos θ) =

1 ∆Bρ
1
z 0= 2
= 2
γ0 Bρref
γ0

10p0
(1 − cos θ)
2.9979Q0 B


γ02 − 1



Q0

Qref
2
γref

−1

−

(3.48)

2 −1
γref



Qref



(3.49)

The full simulation for FRIB Linac Segment 1 (LS1) using XAL for U33+ and U34+ is
given in Figure 3.21. A close up of the end of FS1, which has a chicane consisting of four
80

dipoles, bending 5◦ each, is shown in Figure 3.22. This chicane section is included in the beam
line design to eliminate a direct line-of-sight of the beam path between the stripper foil and
the linac. This precaution protects the accelerating cavities from any stray particles emitting
from the collisions at the stripper foil, which would decrease the cavity’s performance if they
happened to hit the cavity’s surface.

Figure 3.21: FRIB Linac Segment 1 as simulated by XAL for U 33+ and U 34+ charge states,
as well as the result of the two charge states combined.
Notice that in the bottom graph of Figure 3.22 for U33+ , the reference charge state,
x = 0. This is the reference orbit from which the deviation for U34+ is calculated in the red
plot. We can see that this beamlet spreads out by almost 1.5 mm from the reference charge
state, and returns to ∼0 mm by the end of the chicane. We can also note that the beamsize,
σx remains small for each charge state, but the combination grows as they spread apart in
the first dipole, and come back together through the rest of the chicane. The averages of
these charge states as plotted is the sum of the distributions of the individual charge states
81

Figure 3.22: End of FRIB Linac Segment 1 chicane as simulated by XAL for U33+ and U34+
charge states, as well as the result of the two charge states combined.
which is derived in the following manner.
Each charge state can be modeled in each phase space (x, y, z) as a Gaussian distribution
with its average equaling the centroid location, µ, and a standard deviation of σ. These
Gaussian distributions for each charge state can be combined to determine the mean and RMS
of the full distribution. Equation 3.50 is thus the total distribution function. A schematic
representation of the distributions being added is shown in Figure 3.23.
Ni − 12
√ e
ftot (x) =
σ
2π
i
i

x−µi 2
σi

(3.50)

The standard deviation of this distribution is the square root of the normalized variance
of the function ftot , defined by Equation 3.51, where µtot is the mean of the distribution.
The mean of the total distribution is the weighted average of the individual gaussian peaks,

82

Figure 3.23: Schematic representation of three beamlets with means µi , standard deviations
σi , and intensities Ni
as given in Equation 3.52.

2
σtot

=

x2

tot

=

∞
2
−∞ (x − µtot ) ftot (x)}, dx
∞
−∞ ftot (x) dx

µtot = x tot =

i Ni µi
i Ni

(3.51)

(3.52)

Plugging Equation 3.50 into Equation 3.51 and making the substitution u = x − µtot ,
2 is simplified to Equation 3.53. Note here that the fact
with du = dx, the formula for σtot

that integration of a Gaussian distribution over all space being normalized to one has been
used to simplify the denominator.

2 =
σtot

1
i Ni

∞
Ni
√
u2 + 2u (µi − µtot ) + (µi − µtot )2 du
σ
2π
−∞
i
i

(3.53)

∞ n −ax
Performing the integrals using −∞
x e
dx = n+1 results in Equation 3.54. Simplifya
2(n!)

ing the terms, completing the square, and taking the square root yields the final result for

83

the standard deviation in Equation 3.55.

2 =
σtot

1
i Ni

Ni
√
σ
2π
i
i

√

π (1!) 2−1 2σi2

3/2

1
√
+ 2 (µi − µtot ) π (0!) 2−1/2 2σi2

√
+ (µi − µtot )2 σi 2π
(3.54)

σtot =

i Ni (σi

+ (µi − µtot ))2
i Ni

(3.55)

This equation for σtot along with Equation 3.52 for µtot define the distribution of the
charge states. By using these equations, individual charge states can be run through a matrix
modeling code and the resulting individual distributions can be combined to define the total
combined distribution of the charge states. Again, this means of describing the beam is
intended to be used for determining the centroid motion of the combined charge states as
well as their relative distribution size as a means to compare the matrix modeling of multiple
charge state transport to results from multi-particle tracking codes.
Now that this methodology has been implemented in XAL, in order to validate this method
the results must be benchmarked with a multi-particle tracking code, such as IMPACT. The
results from this benchmarking procedure are presented in the following section.

3.5.3

Benchmark

The FRIB Folding Segment 1 (FS1) to the Stripper Foil is an ideal beam line to use as a
benchmark for this method, as it contains both matching cavities and bending dipoles where
the separation of charge states can be seen. The results are shown in Figure 3.24.
This method has been shown to agree with the results from the multi-particle tracking
84

Figure 3.24: Results of XAL and IMPACT calculations for FRIB FS1 to Stripper. The
IMPACT results are shown with dashed lines, while the XAL results are shown with the solid
lines.
code, IMPACT. The XAL run and plotting only takes a few seconds to complete, compared
with IMPACT that can take minutes to run. This makes XAL a feasible option for an on-line
model. With this benchmark complete, the only part left to implement would be to include
this new multiple charge state run capability into an established modeling code such as XAL.
A schema is given in Figure 3.25 which represents how this would be easy to realize. The
left-hand side of the diagram in Figure 3.25 depicts how a single charge state XAL model run
is accomplished presently. The XAL QS1 model run populates the XAL object from which
applications draw the relevant beam modeling parameters.
Multiple charge state XAL model runs may easily be implemented by wrapping the
summing and < z > routines around individual XAL runs of different charge states. This
is depicted in the right-hand side of the diagram in Figure 3.25. In this way, the XAL QS1
model run is used to calculate < z > for the subsequent model runs, XAL QS2 ... XAL
QSN. The final result from each model run is given as input to the summing routine, which

85

Figure 3.25: XAL Schema for adding multiple charge state run capabilities.
populates a single XAL object, mirroring the single charge state XAL model run method.
Consequently, applications would be able to run multi-charge state modeling in XAL
with the advantage that there is no need for major source code editing, as the layer between
applications and the XAL object remains unaltered. Although this method should be simple
to implement, this development has been postponed until the OpenXAL collaboration is able
to use this functionality, in order to reduce duplication of effort as the code develops.

86

Chapter 4
Quarter-Wave Resonator Steering
The focus of this chapter is on the behavior of the low-β accelerating cavities utilized at ReA,
Quarter-Wave Resonators (QWRs). Understanding and characterizing the performance of
these cavity types is crucial for the most efficient operation of not only ReA, but other future
linacs that will use QWRs, such as FRIB.
One attribute specific to QWRs is the unique asymmetric geometry with respect to the
beampipe, unlike other accelerating cavity types which are symmetric with respect to the
beampipe. This asymmetry, typically in the vertical direction, leads to a deflection of the
beam in the direction of the asymmetry. Calculations of the predicted QWR steering effect
as derived from electromagnetic field profiles have been performed, but they had not been
experimentally verified. In order to accomplish this, calculations for the ReA β = 0.041
cavity type have been compared with the results from measurements performed at ReA
characterizing this effect.
The analytic steering predictions are presented in Section 4.1, along with the calculation
for the ReA β = 0.041 and β = 0.085 QWR steering prediction in Section 4.1.3. A detailed
description of the experiments performed at ReA with the β = 0.041 QWR cavity type is
presented in Section 4.2, concluding with analysis and results from these experiments.

87

4.1

Steering Predictions

An analytic approach to the behavior of QWR cavities with respect to vertical steering has
been derived by Facco and Zvyagintsev [54]. They arrive at a formula for the vertical steering
of the beam which consists of three terms.
The first two terms arise from the asymmetry of the QWR structure, depending on the
x-component of the magnetic field, and the y-component of the electric field. The third term,
proportional to the offset of the beam from the central cavity axis, arises from RF defocusing,
as introduced in Section 3.3.3.
Before considering this formula, the choice of using QWR type accelerating cavities is
briefly discussed in Section 4.1.1 to motivate the reasoning behind using an element that
is asymmetric to the beampipe, and therefore produces steering of the beam due to this
asymmetry. Following is a review of the derivation of the vertical steering formula, presented
in detail in Section 4.1.2. A final discussion of the relevant ReA calculations is presented in
Section 4.1.3.

4.1.1

QWR Choice

Quarter-Wave Resonators are the accelerating cavity of choice for low-β beams in linacs
such as ReA and FRIB. Low-β beams require longer RF wavelengths, and therefore lower
frequencies, due to the required spacing between the accelerating gaps.
For optimal acceleration, the beam travels the distance between the gaps in the time that
it takes the RF phase to switch by ∆φ = π (π-mode) [46]. This leads to ∆φ = ωRF ∆t = π,
where ωRF = 2πc/λ, and ∆t = l/v = l/ (cβ). Solving for l, the spacing between the RF gaps,
results in l = βλ/2.

88

To be mechanically stable with a reasonable gap-to-gap distance, l, the RF wavelength, λ,
must be increased when β is low. This results in a lower resonant RF frequency, f , recalling
that λ = c/f . In order to avoid impractically large structures when requiring cavities to
operate at lower frequencies, QWR type cavities are used, as its length is ∼λ/4, as its name
suggests, as opposed to the HWR length of ∼λ/2, as seen in Figure 4.1.

Figure 4.1: FRIB Quarter-Waver Resonator (QWR) and Half-Wave Resonator (HWR) cavity
types shown with optimal β and operational frequencies. The FRIB HWRs are designed for
four times the frequency of the QWRs, resulting in their length being about half as long.

The geometry of the QWR cavity type sets up a standing wave of length λ/4 excited
at the cavity’s resonant frequency. The shape that the resonator must take to achieve this
quarter wave length standing wave results in its asymmetry with respect to the beampipe.
The QWR cavity can be thought of as a coaxial cable of length λ/4 that is shorted at one
end — the circulating magnetic field is maximum at the closed end, while the radial electric
field that accelerates the beam is maximum at the open end. This can be seen in Figure 4.2.
Any asymmetries, whether part of the design or produced during the manufacturing
89

Figure 4.2: Coaxial Cavity Fields: the beam is accelerated at the open end (bottom) where
the electric field is maximized, while the maximum magnetic field circulates at the closed
end (top).
process, of elements with respect to the beam pipe can introduce aberrations to the beam.
These asymmetries introduce fields on-axis that can deflect the beam off its desired trajectory.
This becomes an issue with QWR cavities because of the inherent asymmetry of their design.
There are two direct correction methods to address this issue, either through the design
of the beam port shaping to counteract the field asymmetry, or through the displacement of
the QWR axis during installation to offset the majority of the effect [55]. These methods
may eliminate most of the steering introduced by the asymmetry, however this can be
more complicated mechanically and lead to higher costs in production. Avoiding these
complications with the simplified “straight” beam port design geometries and installing the
QWRs on-axis is often the most realistic and feasible approach. Therefore it is important in
this case to study and understand the degree of the steering in order to be able to analyze
this effect through modeling and aim to mitigate it through the use of corrective steerers.

90

4.1.2

Vertical Steering

The vertical steering induced through a QWR can be derived as a function of the average
accelerating field, the cavity phase, the beam velocity and the vertical offset of the cavity.
The fields that contribute to the change in vertical momentum that are present in a QWR
cavity are Ey and Bx , from the Lorentz force, as seen in Equation 4.1.

∆py = q

Ey (z, t) + βcBx (z, t) dt

(4.1)

The fields Ey and Bx are calculated on axis (x = y = 0), and provide the lowest order
steering calculation for particles entering the cavity near the design trajectory along the
central axis. Recall that these field profiles have been shown in Figure 3.11.
By using the paraxial approximation (pz ≈ p), the vertical steering effect is then given by
∆p

∆y = p y , where ∆py is the change in transverse momentum imparted by the cavity, and p
is the longitudinal momentum. By plugging in p = βγAm0 c and performing the integrals
while assuming only small changes in β through the accelerating gap, such that z ≈ z0 + βct,
we arrive at Equation 4.2.

∆y (β) =

qeEa L sin φ
× DBx (β) + DEy (β)
Am0 c2 TEz (β0 )

(4.2)

Recall that TEz has been defined as the transit time factor in Equation 3.37. The first
term in Equation 4.2 arises from the fact that there is a nonzero x-component of the magnetic
field. The formula for DBx is given in Equation 4.3.

DBx (β) =

GBx
Tg (β, gBx ) sin
βγ

91

πdBx
βλ

(4.3)

The second term in Equation 4.2 is due to the y-component of the electric field. The formula
for DEy is given in Equation 4.4.

DEy (β) =

GEy
Tg (β, gEy ) cos
β 2γ

πdEy
βλ

(4.4)

This vertical electric field is much weaker than the transverse magnetic field, as mentioned
earlier, which results in the strength of this term being an order of magnitude smaller. We
also note this term is also inversely proportional to the square of β, so for low β beams,
such as the case where we are interested, this term does become stronger. The geometrical
constants, GBx and GEy , are given by Equations 4.5.
L/2

c
Bx (z) dz
GBx = 0
L/2
Ez (z) dz
0

GEy =

L/2
Ey (z) dz
0
L/2
Ez (z) dz
0

(4.5)

Tg (β, g) has been defined in Equation 3.37. The gap-to-gap distances, dEz , dEy , and dBx ,
are given by Equation 3.38. The gap lengths, gEz , gEy , and gBx , are given by Equation 3.39.
These terms arise from fields due to the asymmetry of the cavity. These QWR on-axis
fields were presented in Section 3.3.1. It is the result of these two terms that there is a net
vertical steering for QWR cavities, even when the beam is centered in the aperture of the
cavity. Once the beam is not centered in the cavity, however, these are not the only factors
behind beam steering. There is also a contribution due to the offset of the beam entering into
the cavity. This contribution is present due to an effect called RF defocusing, as previously
presented in Section 3.3.3.
This radial transverse defocusing effect depends on r = rˆ
r = xˆ
x + y yˆ. Restricting our
discussion to beam steering in the vertical yˆ-direction, the RF defocusing term depends

92

simply on y, the vertical deviation of the beam from the beam axis. When added to the first
two terms introduced due to the QWR cavity asymmetry, the full equation for the change in
vertical angle as a function of β as presented in [54] is given in Equation 4.6.

∆y (β) =

qeEa L sin φ
× DBx (β) + DEy (β) + Drf (β, y)
Am0 c2 TEz (β0 )

(4.6)

The term due to RF defocusing is given in Equation 4.7.

Drf (β, y) = −

ω

Tg (β, gEz ) sin
2β 3 γ 3 c

πdEz
βλ

y

(4.7)

Note again that it depends linearly with y, and has a negative sign, due to the fact that Er
and Ez have opposite signs, as stated in Section 3.3.3. This term is typically larger than DBx
and DEy at low velocities, but decreases at higher velocities. This can be seen in Figure 4.3
where each D-factor term is plotted across a typical β range calculated at ReA β = 0.041
cavity accelerating parameters of Ea = 16 MV/m and φ = −90 ◦ . The sum of the three
terms is also plotted.
With all three terms introduced, there are some observations to be made. The magnitude
of the terms has been presented, and it is important to recognize that when measuring the
effect in an experiment, all terms will be present at once. Another important note is that the
three terms that contribute to Equation 4.6 do not all have the same sign. The last term has
a minus sign, which indicates that there is a “sweet-spot” where the contribution from each
term will cancel each other out. This is predicted to be at a slightly positive vertical offset.
This “sweet-spot” as predicted by theory may allow for the determination of the alignment
of the cavity center with respect to the diagnostics center with a beam-based measurement.

93

Figure 4.3: The calculated ReA β = 0.041 cavity D-factors from Equation 4.6, DBx , DEy ,
and Drf , with Ea = 16 MV/m and φ = −90 ◦ along with the sum of the three terms.

4.1.3

ReA Calculation

A calculation of Equation 4.6 has been performed and the results are shown here. As
mentioned previously, the 3D fields from the ReA cavities are calculated with CST Microwave
Studio. The on-axis fields are used to calculate the values for the vertical steering result. The
β = 0.041 cavity’s on-axis fields used for the ReA calculation are shown in Figure 3.11. The
FRIB cavity fields are similar, and result in comparable, but not identical, calculations. Here
we will focus solely on the calculations for the ReA cavities, to be compared with experiments.
A sample calculation of the change in vertical deflection from Equation 4.6 is presented
in Figure 4.4. Typical accelerating parameters of Ea = 16 MV/m and φ = −20◦ have been
used for this calculation. The two graphs depict the result for the ReA β = 0.041 QWR
fields, along with the calculations performed for the ReA β = 0.085 QWR fields. Each graph
94

includes a vertical dashed line at the optimal beta of the cavity. The central theory line, in
blue, plots the dy = 0 mm case, while the other lines above and below show the result of a
±1 mm vertical beam offset from the ideal trajectory through the center of the cavity.
2

β=0.041

1.5

β=0.085

y= 1 mm
y= 0 mm
y=−1 mm

1.5

1

1

0.5

0.5

∆y’ (mrad)

∆y’ (mrad)

2

y= 1 mm
y= 0 mm
y=−1 mm

0
−0.5

0
−0.5

−1

−1

−1.5

−1.5

−2

−2
0.04

0.06

0.08

0.1

β (v/c)

0.12

0.14

0.16

0.04

(a)

0.06

0.08

0.1

β (v/c)

0.12

0.14

0.16

(b)

Figure 4.4: Calculation of ∆y vs. β for the ReA (a) β = 0.041 and (b) β = 0.085 cavity
types, plotted for centered and ±1 mm vertical offset scenarios. A typical accelerating cavity
setting is used: φ = −20◦ and Ea = 16 MV/m. The beam species is H2 with q/A = 0.5. The
x-range for each plot is from β = 0.036, corresponding to the RFQ energy of 0.6 MeV/u, to
β = 0.16, corresponding to the highest beam energy for H2 at ReA, 12 MeV/u.

This calculation was performed at Ea = 16 MV/m; however, many of the ReA cavities
can be operated up to fields of Ea = ∼30 MV/m. This will almost double the ∆y as shown
in Figure 4.4. Another factor to consider is that the steering will increase when the cavities
are operated in rebunching mode at φ = −90◦ , instead of accelerating at φ = −20◦ , but this
is typically done at much lower fields which then reduces the effect. Each different case will
change the overall magnitude, but the trends of the curves remain the same. Keeping this
in mind, we can use the case shown in Figure 4.4 to look more into this steering effect as it
relates to the ReA QWR cavities specifically.
It is important to note that the regions in which the β = 0.041 and β = 0.085 cavities
operate differ, and the calculated deflections are much larger in the lower β regions where
95

β = 0.041 cavities are used. This is due to the effect being inversely proportional to β, with
slower beams being deflected more. This is important to note, because it means that the
effect is greatest with the lower velocity, lightest beams.
In these plots, the same offset produces a much larger effect at the β = 0.041 cavity’s
optimal β than it does at the β = 0.085 cavity’s optimal β. From the β = 0.041 plot in
Figure 4.4, it is noted that at β = 0.041, the difference between the ±1 mm cases is ∼2 mrad,
while the same difference at β = 0.085 in the β = 0.041 plot in Figure 4.4 only produces a
difference of ∼0.5 mrad. This indicates that the β = 0.041 cavities are much more sensitive
to an offset from the beam line center, which we also hope to measure.
Plotted in Figure 4.5a is the predicted change in vertical beam angle normalized by the
cavity field for a ReA β = 0.041 cavity operating at −90 ◦ phase vs. β across a range of
energies available at ReA.
2

0

y= 6 mm
y= 3 mm
y= 0 mm
y=−3 mm
y=−6 mm

6
4

∆y’ (mrad)

1

∆ y’/Ea (mrad/(MV/m))

8

y = 6 mm
y = 3 mm
y = 0 mm
y = −3 mm
y = −6 mm

2
0
−2
−4

−1

−6
−8

−2
0.036

0.038

0.04

0.042

0.044

β (v/c)

0.046

0.048

0.05

0

0.052

0.5

1

1.5

2

E (MV/m)

2.5

3

3.5

4

a

(a)

(b)

Figure 4.5: Left (a): Plot of ReA β = 0.041 calculation ∆y vs. β for offsets dy = −6, −3, 0, 3,
and 6 mm. Right (b): ∆y vs. cavity amplitude Ea for offsets dy = −6, −3, 0, 3, and 6 mm
for beam energy Ek = 0.6 MeV/u, corresponding to β = 0.036.

From the plot in Figure 4.5a, the simple observation can be made that a negative kick

96

is generally expected from the cavity for negative offsets. Similarly, a positive kick from
the cavity is expected for positive offsets. Recall that the central theory line, plotted for
dy = 0 mm, however, will not equal zero, as confirmed in the graph. From the trend of
the curves, note that lower beam energies, corresponding to lower β, result in a larger net
displacement. Again, comparing the plots it is also clear that larger offsets will result in
larger displacement.
For the calculations in Figure 4.5a, the vertical kick has been normalized by Ea . Varying
the cavity amplitude for a fixed energy instead results in the trends shown in Figure 4.5b.
For a fixed energy of 0.6 MeV/u, corresponding to β = 0.036, Figure 4.5b plots ∆y vs. Ea .
Looking at the dy = ±3 mm plots at Ea = 4 MV/m, the expected ∆y is 3 mrad and −4 mrad
respectively. This demonstrates that the deflection is not symmetric. Adding a slightly
positive offset dy to the beam entering the cavity can minimize the deflection of the beam due
to the QWR steering. This is the “sweet-spot” indicated in section 4.1.2. This theoretical
offset to eliminate steering in the QWR cavities at ReA for 0.6 MeV/u beams is ∼0.2 mm.

4.2

Experiment

The vertical steering can be studied with a carefully designed experiment at ReA, which
exclusively uses QWRs in its linac. The layout of ReA during this experiment is depicted in
Figure 4.6. The deflection can be maximized by using a low mass beam at low energy and
operating one of ReA’s relatively low beta cavities in rebunching mode. The deflection can
be scaled down by either increasing the beam mass, increasing the beam energy, lowering the
cavity’s amplitude or operating the cavity at an accelerating phase instead of rebunching,
where theory predicts the amplitude of the steering is greatest. For stability, the test cavity

97

is left at a fixed phase of -90◦ . As energy is one of the variables under investigation in this
experiment, it is also beneficial to measure the deflection at the -90◦ phase setting where the
energy of the beam is unchanged.

Figure 4.6: ReA Layout from 2013, including EBIT charge breeder, Q/A section, commissioning source, room temperature RFQ and superconducting linac including two cryomodules
containing β = 0.041 accelerating cavities. Note the drift tube after the second cryomodule
is the location for the third cryomodule containing β = 0.085 accelerating cavities, to be
installed.
The cavity’s amplitude is varied to compare its response with the linear relationship as
predicted by Equation 4.6. The offset of the beam going into the cavity is also varied to see
the dependence of the third term in Equation 4.6. These results are then compared with the
calculations shown in Figure 4.5b. Lastly, this measurement is performed at multiple beam
energies in order to compare the experimental measurement with the theoretical predictions
as seen in Figure 4.5a.
98

4.2.1

Procedure

The procedure for this measurement begins with tuning the linac to the desired energy with
a parallel beam through the long drift tube following the linac. The four measurement
energies are Ek =0.60 MeV/u, 0.85 MeV/u, 1.08 MeV/u, and 1.24 MeV/u. The last cavity is
off during this phase. The parallel beam is achieved using the two horizontal and vertical
steerers located before the last cavity in the second cryomodule, labeled in Figure 4.7 as the
test cavity. Each solenoid in the linac is equipped with horizontal and vertical steerers. The
two cavities before the test cavity remain off during each measurement point, resulting in a
drift space between the two solenoids utilized to establish a parallel beam.
The beam diagnostic locations at which the beam profiles are measured are also labeled
in Figure 4.7. These two beam diagnostic locations will be referred to as L092 and L110
respectively. The energy is measured using the dipole-slit procedure: a dipole magnet is
located directly after L110 (not shown). Adjusting the calibrated dipole field to achieve a
beam centered at the collimator slits determines the energy of the beam. The centroid at the
two locations after the linac is measured with the 45◦ slit scanners to establish the baseline
beam trajectory.

Solenoids

Test Cavity

Diagnostic Locations

Cavities

Figure 4.7: Schematic of ReA Linac as commissioned in 2013, before the installation of
the third cryomodule. Square elements represent superconducting solenoids, while circular
elements represent quarter-wave resonators. The drift space between the diagnostics after
the second cryomodule is ∼5.2 m.

99

The test cavity is then turned on in rebunching mode and the slit scanners measure
the beam centroid position at the two downstream locations. The cavity amplitude is then
changed multiple times, measuring the centroid positions each time. After these measurements,
the cavity is turned off and a new parallel beam is achieved with a vertical offset. The same
procedure is repeated for multiple different offsets. The method for establishing a vertical
offset using the two correctors before the test cavity is described in following section.

4.2.1.1

Vertical Offsets

Figure 4.8 is a compilation of the data taken in two different measurements, at 1.3 MeV/u
and 1.24 MeV/u beam energy. The plot shows the extrapolated position of the beam, ym
from measuring the position at L092, y1 and the angle y from measuring the position at
L110, y2 . The different offsets were produced by changing the corrector settings in front of
the test cavity can, as seen with the variation of the points in the y-direction.
8

05−07 Ek=1.30MeV/u

Beam Y at Cavity Midpoint (mm)

10−29 Ek=1.24MeV/u
4

0

−4

−8
−8

−4

0

4

Beam X at Cavity Midpoint (mm)

8

Figure 4.8: Beam position at cavity midpoint as projected from measurements at two
downstream locations. Error bars are the RMS beam size. Note the improvement of the
vertical separation and consistency of horizontal position achieved with the 10-29 data set
which utilized the response matrix measurement technique.

100

The procedure during the first measurement in May was to steer the beam up by
introducing small vertical steps in the slits, and tuning on the beam transmission. This stepby-step method resulted in a large shift in the x-direction, since the horizontal slit (measuring
the vertical position) does not constrain the beam in the x-direction. It is preferable for this
experiment to isolate the vertical motion of the beam.
A new method was therefore devised to determine the beam’s response from each individual
corrector and use the resulting response matrix to achieve a purely vertical displacement
of the beam. The difference between the beam’s position at the two downstream locations,
−→
∆x given in Equation 4.8, was measured as a function of the change in the four individual
−→
corrector settings, ∆I.




 ∆x1 




 ∆y1 

−→ = 


∆x


 ∆x2 




∆y2





 ∆IH1



−→ =  ∆IV 1

∆I

 ∆IH2


∆IV 2












(4.8)

By changing each corrector individually and recording the response of the beam, the matrix
M in Equation 4.9 is mapped out element by element.
−→
−→
∆x =M ∆I
−→
−→
∆I =M −1 ∆x

101

(4.9)

Inverting M and plugging in the desired shift in beam position, given in Equation 4.10,
established a ratio of the corrector settings to use in order to offset the beam vertically by a
desired amount.





 0 




 ∆y 


−→ =


∆x


 0 




∆y

(4.10)

This procedure was more robust than the previous method, and successfully prevented the
beam from moving in the x-direction, as can be seen in the October data in Figure 4.8. With
this response matrix approach a large vertical shift of the beam with very little movement in
the horizontal direction can be achieved, as desired.

4.2.1.2

Profile Measurement

Beam profile scans at both diagnostic locations discussed previously, L092 and L110, determine
the centroid and RMS beam size. These measurements are performed with 45◦ slit scanners.
The 45◦ slit scanner plate is shown in top left of Figure 4.9. This plate intercepts the beam at
a 45◦ angle, allowing the horizontal slit to measure the vertical beam profile, and the vertical
slit to measure the horizontal beam profile. A 5 mm collimation hole is also drilled in the
end of the plate for beam tuning. Sample scans are shown in Figure 4.10.
The scan is performed by fully inserting the drive (drive position=−100 mm) and moving
the drive in 1 mm steps until it no longer intercepts the beam (drive position=40 mm). With
this direction of running the scan, the peak on the far left of the plot corresponds to the
vertical profile of the beam. The next peak corresponds to the horizontal profile of the beam,
while the peak centered at drive position = 0 mm corresponds to the 5 mm hole. The drive

102

Figure 4.9: 45◦ Beam Profile Diagnostic setup. In the top left, the plate is shown with three
slits and one aperture in the order that they intercept the beam during the profile scans: (1)
the horizontal slit measuring the vertical beam profile, (2) the 45◦ slit for a measure of the
coupling between the x and y planes, (3) the vertical slit measuring the horizontal beam
profile, and (4) the 5 mm aperture for beam tuning.
at L092 has an additional central slit perpendicular to the motion of the drive to measure
the 45◦ profile of the beam, which accounts for the central peak in the L092 drive scan that
is absent from the L110 scan.
In Figures 4.10a and 4.10b, appropriate regions have been used for fitting a Gaussian curve
to each peak. A low threshold level is subtracted from the raw data before the fitting analysis
is performed to reduce background noise fluctuations. The resulting Gaussian curves fit to
each peak is shown in the dashed red lines. The fit parameters from the Gaussian curve also
give the centroid and RMS beam size for the horizontal and vertical profiles. Equation 4.11
gives the equation for a Gaussian function.
103

L092_dy_0mm_L091_0MVpm.dat

450
400

400

350

300

250
200

200
150

100

100

50

50
0
−120

0
−100

−80

−60

−40

Raw Drive Position (mm)

−20

0

20

Yc = −0.25 mm
Y = 0.39 mm
rms

250

150

−120

Xc = 0.07 mm
X = 1.39 mm
rms

350
Yc = 0.95 mm
Yrms = 0.46 mm

Current (pA)

Current (pA)

300

L110_dy_0mm_L091_0MVpm.dat

450

X = −0.78 mm
c
Xrms = 0.60 mm

−100

−80

−60

−40

−20

Raw Drive Position (mm)

0

20

(b)

(a)

Figure 4.10: Plots of beam current vs. slit position for raw data from 45◦ slit scan taken at
(a) L092 and (b) L110, and the Gaussian fit results for beam position and RMS beam size.
Note that the 45◦ plate at L110 does not have the 45◦ slit, which accounts for the missing
peak in (b).

y=

−(x−µ)2
Ae 2σ2

(4.11)

The raw drive position, as plotted in Figure 4.10, is along the 45◦ axis, so a simple transformation has to be made to extract the absolute x and y positions and beam sizes. The raw
position, given by µ in Equation 4.11,

xc , yc =

µx,y − doffset ∓ dref
√
2

(4.12)

is centered and then the relative distance to the slit position is either added or subtracted to
√
get the position of the slit measuring either y or x. This value must then be divided by 2
to get the distance along the perpendicular axis instead of along the 45◦ angle. Similarly,
in order to get the RMS beam size, xRMS or yRMS , σ as determined by the Gaussian fit in
√
Equation 4.11 must be divided by 2, as given in Equation 4.13.
104

σx,y
xRMS , yRMS = √
2

(4.13)

These transformations are given in Equations 4.12 and 4.13, with a schematic representation
of the distances dref and doffset shown in Figure 4.11.

Figure 4.11: 45◦ Beam Profile Diagnostic Plate shown with distances from Equation 4.12,
dref and doffset , labeled.
In the scan in Figure 4.10a, the centroid position of the beam at L092 is about (−0.8 mm,
1.0 mm) with RMS size (0.6 mm, 0.5 mm). After the long drift space, the beam profile is
measured again at L110, shown in the scan in Figure 4.10b. The centroid position of the
beam has changed to about (0.1 mm, −0.3 mm) with RMS size of about (1.4 mm, 0.4 mm).
The change in x-position of ∼0.9 mm along a ∼5.2 m drift corresponds to x ≈ 0.17 mrad. In
y the change in position of ∼−1.3 mm corresponds to y ≈ −0.25 mrad.
All measurement data sets include these initial measurements taken with the test cavity
turned off. During the analysis, all measurement points will subtract these small initial angles
in order to measure changes in angle, instead of absolute angle. Multiple calibration points
are also measured through the run time to ensure the stability of these initial conditions.
105

4.2.2

Beam Axis Definition

The measurement of the beam centroid position at both diagnostic locations enables the
computation of the angle of the beam from the change in position between them, as they are
separated by a drift space. This is schematically represented in Figure 4.12. The difference
of the measurements at the first and second diagnostics determine the beam’s deflection,
y −y
y = 2d 1.

Figure 4.12: Schematic diagram of relevant measurement quantities depicting the difference
between the diagnostic axis and the cavity center axis.

The beam enters the cavity at a certain unknown vertical offset, with the cavity off. The
parallel beam is measured at the two diagnostic locations, but as noted in the previous
section, the beam may have an initial angle measured with the cavity off (not shown). This is
subtracted off from each measurement with the cavity turned on in order to look at relative
changes, ∆y . Figure 4.12 shows a simplified perspective where the incoming angle is zero.
Each beam measurement, y1 and y2 , are measured with respect to the diagnostic axis,
with zero defined as the center of the beampipe as defined by the alignment of the diagnostic
elements. Once y is calculated from y1 and y2 , the position of the beam at the center of
the test cavity, ym , can be extrapolated. This quantity is also calculated with respect to the
106

diagnostic axis.
In Equation 4.6 for the QWR steering prediction, vertical offsets y are measured with
respect to the cavity center. There exists a small unknown vertical displacement of the cavity
center with respect to the diagnostic axis, δ. The relationship y = ym + δ, is true for all
measurements, including across the different energy measurement, as δ is a constant.
A chi-squared minimization of the entire data set calculates the optimal value of this
quantity to be δ =1.29 mm. This result is used to convert the raw measured values, ym , to
displacements with respect to the cavity center, y. In this way, all the data results shown in
the following section along with the calculated predictions for a given offset y are consistent,
defined with respect to the cavity center.

4.2.3

Results and Analysis

Figure 4.13 shows the results of the measurement data for beam energies 0.6 MeV/u,
0.85 MeV/u, 1.08 MeV/u, and 1.24 MeV/u. The plots show ∆y vs. Ea (cavity amplitude setting) for various offsets, which were produced by changing the corrector settings in
front of the cavity. With the cavity off, the beam positions at the two downstream locations
are measured to give the baseline angle, which is then subtracted from each measurement
point to get a change in angle, ∆y . This is why all the data sets overlap at the origin. This
baseline angle is also used along with the position measured at the first measurement location
to calculate the beam’s position at the cavity center. These calculated values are given in the
legends of the figures, as well as in the ∆y column of Table 4.1, and are used to calculate the
prediction for each data set.

107

E =0.60MeV/u β=0.036

y = 4.65 mm
y = 2.41 mm
y = 0.81 mm
y = −1.29 mm
y = −3.28 mm

1

∆y’ (mrad)

0

k

2

y = 5.07 mm
y = 3.24 mm
y = 0.04 mm
y = −0.92 mm
y = −2.94 mm
y = −5.42 mm

1

∆y’ (mrad)

E =0.85MeV/u β=0.043

k

2

−1

0

−1

−2

−2
0

0.5

1

1.5

Ea (MV/m)

2

2.5

3

0

0.5

1

(a)

Ek=1.24MeV/u β=0.051

k

2

y = 4.42 mm
y = 3.07 mm
y = 1.12 mm
y = −1.16 mm
y = −3.03 mm

2.5

3

0

y = 6.75 mm
y = 3.37 mm
y = 0.09 mm
y = −2.57 mm
y = −5.90 mm

1

∆y’ (mrad)

1

∆y’ (mrad)

2

(b)

E =1.08MeV/u β=0.046
2

1.5

Ea (MV/m)

−1

0

−1

−2
0

0.5

1

1.5

Ea (MV/m)

2

2.5

−2

3

0

(c)

0.5

1

1.5

Ea (MV/m)

2

2.5

3

(d)

Figure 4.13: Plots of ∆y vs. Ea for beam energies (a) Ek = 0.6 MeV/u (β = 0.036), (b)
Ek = 0.85 MeV/u (β = 0.043), (c) Ek = 1.08 MeV/u (β = 0.046), and (d) Ek = 1.24 MeV/u
(β = 0.051). Lines correspond to the predicted response calculated using cavity design data
and experimental conditions.

108

The vertical error bars are determined from the statistical error inherent in fitting the
beam’s profile with a Gaussian. A larger RMS size of the beam at the measurement locations,
or features of the profile such as a double peak, will increase the uncertainty to which the
centroid of the beam is known at each location.
These errors in the beam’s position are then added in quadrature and divided by the
distance between the two measurement points in order to obtain the error on the angle. The
horizontal errors are determined from the cavity’s phase and amplitude stability measurements,
yielding an RMS value of 0.005 MV/m.
When comparing similar vertical offsets from the cavity center between the plots, as
expected, a larger kick was experienced by the lower energy beams. The cavity amplitude
also could not be increased as much for the lower energy beams as it could for the higher
energy beams, due to beam losses.
The setup at ReA limited the maximum possible measured ∆y to ∼2 mrad which corresponds to a ∼10 mm shift over the 5.2 m drift between measurement points. Although the
slit drives can measure a range of ∼±15 mm, the beam may have already been displaced by
up to ∼5 mm from the center of the cavity, resulting in a maximum measurable displacement
of ∼10 mm in that same direction.
An important observation at this point that each data set is linear with varying the cavity
amplitude, as we expect from Equation 4.6. It is also encouraging that the different offsets at
the cavity fit fairly well within the error bars to its respective theoretical prediction.
The result from plotting ∆y normalized by the amplitude of the test cavity setting
for each point as a function of beam position at the midpoint of the test cavity is shown
in Figure 4.14a. Note that the slopes of the relative energy measurements are ordered as
expected, and agree reasonably well with theory.
109

2
β = 0.036
β = 0.043
β = 0.051
β = 0.046

∆y’/Ea (mrad/(MV/m))

1

0

−1

−2
−8

−6

−4

−2

0

2

Cavity Center y (mm)

4

6

8

(a)
2

y = 6 mm
y = 3 mm
y = 0 mm
y = −3 mm
y = −6 mm

∆y’/E a (mrad/(MV/m))

1

0

−1

−2
0.036

0.038

0.04

0.042

0.044

β (v/c)

0.046

0.048

0.05

0.052

(b)
Figure 4.14: Left (a): Plot of ∆y /Ea vs. y at the midpoint of the test cavity, including each
measurement data set. Note the agreement in the slopes for the different energies. Right (b):
Plot of ∆y /Ea vs. β including each measurement data set.

110

Another way to visualize the data is as a function of β, as shown in Figure 4.14b. Here
the data points indicate an average of a measured data set (therefore plotting one point per
offset). Open circles are the corresponding theory values for the zero-data measured offset at
the cavity. Dashed lines are equally spaced theory lines for reference. Vertical error bars here
are statistical, from the uncertainty in measurement of centroid from the Gaussian fitting of
the beam profiles. The horizontal error bars on β correspond to a ±1% uncertainty on the
energy measurement.
The full data set of values plotted in Figure 4.14b is available in Table 4.1. Note that
numbers have been rounded to the hundredths place after all calculations.
Table 4.1: Summary of cavity steering data at β = 0.036, 0.043, 0.046, and 0.051.

β
0.036

0.043

0.046

0.051

∆y
∆y /Ea Meas. ∆y /Ea Calc.
(mm) (mm/MV/m) (mm/MV/m)
5.07
3.24
0.04
-0.92
-2.94
-5.42
4.49
2.26
0.66
-1.44
-3.44
4.26
2.92
0.97
-1.31
-3.18
6.60
3.22
-0.06
-2.73
-6.05

1.54
0.91
-0.09
-0.39
-1.02
-1.53
0.77
0.37
-0.01
-0.33
-0.76
0.72
0.35
0.11
-0.25
-0.55
0.75
0.34
0.05
-0.30
-0.68

±
±
±
±
±
±
±
±
±
±
±
±
±
±
±
±
±
±
±
±
±

0.16
0.11
0.10
0.10
0.11
0.19
0.16
0.18
0.19
0.16
0.14
0.06
0.07
0.06
0.06
0.07
0.10
0.09
0.09
0.09
0.09

111

Diff

Diff/Error

1.40 0.14
0.87 0.03
-0.05 -0.04
-0.33 -0.06
-0.91 -0.12
-1.62 0.09
0.74 0.03
0.35 0.01
0.08 -0.08
-0.29 -0.03
-0.64 -0.13
0.57 0.15
0.38 -0.02
0.10 0.01
-0.22 -0.04
-0.48 -0.07
0.64 0.12
0.30 0.04
-0.03 0.08
-0.30 0.00
-0.63 -0.05

0.89
0.31
-0.45
-0.62
-1.09
0.50
0.20
0.07
-0.44
-0.22
-0.87
2.39
-0.37
0.12
-0.58
-1.13
1.17
0.45
0.86
0.04
-0.55

The four sections of this table, defined by the β column, correspond to the four measurement sets shown in Figure 4.13. Each row corresponds to a set of data taken with a distinct
combination of corrector settings. The offset produced by these settings when Ea =0 MV/m
is given in the ∆y column. Note that these are the same values used to create the lines shown
in Figure 4.13.
The ∆y /Ea Measurement column gives the average of all the data points for this given
offset (after normalization by Ea ). The stated error is the root of the sum of the squares of
the errors for that measurement set normalized by the number of measurement points.
The ∆y /Ea Calc. column gives the result of the calculation using the offset from the ∆y
column, while the Diff column gives the difference between the Measurement and Calculation
columns. The Diff/Error column can be interpreted as a measure of how far off the data
value is from the theory value with respect to the given error.

4.3

Conclusions and Outlook

A set of measurements has successfully been completed at ReA to study the vertical kick of the
accelerating cavities, quarter-wave resonators (QWRs), in the linac. This study investigated
the effect of the beam’s vertical offset within the cavity, the cavity amplitude, and the beam’s
energy on the vertical steering. This is the first systematic measurement of QWR steering
and its dependence on these various factors.
The basic procedure involves measuring the beam centroid at two downstream locations,
from which the beam’s vertical trajectory can be calculated. The beam’s position at the
center of the cavity is then used to calculate the predicted deflection of the beam, which
can then be compared with the measurements. The prediction indicates that increasing the

112

cavity accelerating field, Ea , as well as increasing the offset of the beam from the central axis
will result in a larger kick, ∆y .
In practice this has been known, but this is the first time that a systematic measurement
of QWR beam steering has been performed. Having found a good agreement between the
measured data and prediction, this effect can now confidently be included into future beam
line modeling of ReA, and will be of relevance to other linacs using QWRs for low energy
acceleration, such as the Facility for Rare Isotope Beams (FRIB).
These initial measurements have focused on the last β = 0.041 cavity of ReA, but in
principle the established procedure could be performed with other cavities at ReA as well. In
the future this method may prove to be useful in obtaining a beam-based measurement of
the relative alignment of the cavities with respect to one another.
The only requirement for this method is that all elements after the test cavity be off in
order to produce a drift space through which the beam’s outgoing angle may be measured.
This becomes more difficult with earlier cavities in the cryomodule, as the drift space to the
first diagnostic increases.
It may be possible to use corrector fields after the test cavity to help transport the beam,
as these weak fields induce a constant deflection on the beam, and the calculation is concerned
with changes in deflection, rather than the absolute deflection. These corrector fields would
have to be taken into account when calculating the beam’s position in the cavity.
Overall, a very promising method now exists for determining the deflection of the beam
from the QWR accelerating cavities, which can be used to inform future experiments and
modeling efforts. This experimental verification of the steering predictions is an important
step in fully understanding and characterizing the performance of the low-β accelerating
cavities at ReA and FRIB.
113

Chapter 5
Conclusions
The various topics covered through this work are united under a central theme of modeling
accelerator systems. From simulations of error propagation through an ideal design lattice,
to experimental measurements of beam steering due to asymmetric accelerating cavities, the
overall focus remains on accurately modeling accelerators in order to provide understanding,
inform decisions, and ultimately improve performance.

5.1

Summary

There are many methods used for modeling accelerators. Multi-particle tracking and matrix
transport are two that have been utilized throughout this work. The modeling of beam line
elements including quadrupoles, solenoids, and coaxial accelerating cavities have been shown
in detail, including considerations determined through the modeling studies.
The continuous development of online modeling of accelerator systems is enormously
important for the future of accelerator physics. This crucial component is what will allow
heavy ion accelerators of the future to be successful in running many different types of
beam, including multiple charge states simultaneously, and minimize tuning times to increase
efficiency.
This effort requires a method for modeling multiple charge states through envelope

114

calculations, which has been proposed and benchmarked with multi-particle tracking results.
A series of experiments has been conducted at ReA which measured the vertical quarter
wave resonator steering. The results have been analyzed and compared with theory predictions.
The measured data from these experiments consists of changes in the vertical deflection,
δy as a function of cavity accelerating field, Ea , for various energies. The cavity is set to
-90◦ phase.
The raw data from the measurements initially does not correspond with the steering
predictions, as shown in Figure 5.1a. This is due to the raw measurement data being taken
with respect to an axis defined by the beam diagnostics, whereas the steering predictions are
made with respect to the central cavity axis, as discussed in Section 4.2.2.

y = 5.07 mm
y = 3.24 mm
y = 0.04 mm
y = −0.92 mm
y = −2.94 mm
y = −5.42 mm

1

∆y’ (mrad)

0

k

2

y = 6.36 mm
y = 4.53 mm
y = 1.34 mm
y = 0.37 mm
y = −1.65 mm
y = −4.12 mm

1

∆y’ (mrad)

E =0.60MeV/u β=0.036

Ek=0.60MeV/u β=0.036

2

−1

0

−1

−2

−2
0

0.5

1

1.5

Ea (MV/m)

2

2.5

3

0

(a)

0.5

1

1.5

Ea (MV/m)

2

2.5

3

(b)

Figure 5.1: Left (a): Raw experimental data of ∆y vs. Ea at Ek = 0.6 MeV/u (β = 0.036).
Right (b): The same data set with a constant vertical offset, δ =1.29 mm, applied during the
analysis.

Introducing a single constant vertical offset, found to be δ =1.29 mm through a leastsquares fitting, lines up the experimental data with the predicted calculations, as seen in
Figure 5.1b. The same constant δ works to line up all the data sets regardless of energy, not

115

just the shown example at Ek = 0.6 MeV/u, as it arises from the experimental setup and is
therefore constant throughout all measurements.
This is a powerful statement for the verification of the steering predictions. As this is the
first time such a systematic measurement of the QWR vertical steering has been performed,
it is important to note that across all energies measured, and for different tunes that include
varied offsets through the test cavity, a single constant δ may be introduced to the analysis
to produce a near perfect match of the measured data with the steering predictions.

5.2

Outlook

There is much work that can be continued into the future. Accurate online modeling of
accelerator systems is an increasingly important objective to achieve for machines such as
ReA and FRIB that require flexible tuning routines and algorithms. Other challenges such as
incorporating multiple charge state tracking and QWR steering effects into an online model
must be handled sufficiently in the hopes of attaining the best possible working model.
The possibility of incorporating multiple charge state modeling into XAL is an idea for
expanding XAL’s future online model capabilities. Being able to track multiple charge states
in an online model will be especially important for machines such as FRIB.
In the future it is entirely conceivable that quarter wave resonator steering modeling could
be included in XAL’s online model. This effect, especially now that it has been measured
and confirmed experimentally, is critical to introduce into modeling codes if a model of the
centroid motion through QWRs and corrector magnets is to be achieved. This is especially
significant for modeling objectives such as orbit correction and scaling of linac settings.

116

APPENDICES

117

Appendix A
Additional Cavity Measurement
On October 30th, 2013, during the preliminary tuning for the Ek = 1.08 MeV/u measurement
using the L091 test cavity, it was noted that by chance the final tune resulted in the last
solenoid operating at a very low current compared with its nominal value. This solenoid,
L090, had a setting of ∼4 A, while the other solenoids were tuned to values in the range of
10-20 A.
As discussed in Section 4.3, the only additional requirement necessary to use this same
measurement technique for the cavities prior to the L091 test cavity is that there is a drift
space between the cavity and the diagnostics. As the L090 solenoid field was already minimal,
further tuning was pursued until this solenoid and its correctors could be turned off completely.
This is the tune that was used to measure the L091 test cavity at Ek = 1.08 MeV/u on
October 31st, 2013.
After this measurement was completed with the L091 cavity, the same measurement
procedure was performed with the L089 cavity. However, on the last data point of the L091
data set there was an abrupt change in the L092 drive reading, as seen in Figure A.1a. There
was no corresponding shift in the L110 drive at the same setting, which was taken immediately
after the L092 measurement at Ea = 0.5 MV/m, indicating that this was an issue with the
L092 drive instead of a shift of the beam.
Indeed, the L092 drive failed after the final Ea =1.5 MV/m measurement point of the
118

(b)

(a)

Figure A.1: Final measurement of L091 data set on 10-31. The Ea = 0.5 MV/m L092 data
point had all peaks shifted by the same amount. The dashed lines indicate where this overall
shift has been manually accounted for in the analysis. L110 data indicates it was an issue
with the L092 drive, not a sudden change in the beam, for instance from an instability in the
LB source.
L089 data set, ending this data set measurement prematurely. This nearly complete data set
is shown in Figure A.2 with two missing measurement points in the y = 4.21 mm data set.
The hope was to complete this measurement set with the L089 cavity in addition to the
L091 measurement and to extract a relative offset of the cavities with respect to each other,
as previously proposed in Chapter 4. By using the procedure introduced in Section 4.2.2 to
find δL091 and δL089 , the difference ∆ = δL091 − δL089 , would be equal to the relative vertical
offset of the cavities with respect to one another.
Using the raw data from both measurement to calculate δL091 and δL089 , it is found that
∆ = 0.97 mm. With this being a positive value, δL091 > δL089 , therefore the L089 cavity is
dy = 0.97 mm above the L091 cavity. This value is suspect, however, due to the change in
the L092 drive at the end of the L091 data set. When this shift of 0.83 mm is accounted for,
the difference becomes ∆ = 0.14 mm. This would indicate that L089 is dy = 0.14 mm above
L091. There is uncertainty about this result because of the malfunctioning L092 drive during
the measurement. It would therefore be beneficial to take the L091 and L089 (and possibly

119

L089 E =1.08MeV/u β=0.046
k

2

y = 4.21mm
y = 1.99mm
y = −0.27mm
y = −1.59mm
y = −3.81mm

∆y’ (mrad)

1

0

−1

−2
0

0.5

1

1.5

Ea (MV/m)

2

2.5

3

Figure A.2: Plot of ∆y vs. Ea for beam energies Ek = 1.08 MeV/u (β = 0.046). Lines
correspond to the predicted response calculated using cavity design data and experimental
conditions.
even L088) data again without interruption from the L092 drive.

120

Appendix B
L092 Drive Calibration
The results and analysis shown in the following section are for the QWR steering measurements
done at beam energies of 1.3 MeV/u for which β = 0.053, and 0.6 MeV/u corresponding to
β = 0.036.
The first measurement in the Quarter-Wave resonator steering series of experiments took
place on May 6-7, 2013. The second measurement took place on August 29, 2013. The
measurement in May was performed at 1.3 MeV/u, while the August measurement used the
RFQ energy of 0.6 MeV/u.
During the initial run in May, the slit drive at L092 would not maintain its calibration,
and therefore needed to be recalibrated after each slit scan. During preparations for the
August experiment, the same slit drive at L092 became inoperable, and a substitute drive
had to be installed in its place. It was aligned by eye, with official alignment to take place
at a later date. Due to these variables, it was necessary to introduce different calibration
offsets for each measurement set in the analysis of the L092 data. The L110 data from
both measurement sets used the original calibration during the analysis, as the L110 drive
functioned properly during both measurements.
Changing the dref constant moves the x and y data in opposite directions, as seen through
Equation 4.12. By the same equation, we see that changing the doffset constant results in the
x and y data moving in the same direction. By adjusting these two constants, it is possible
121

to fit the May data from L092 to the theory predictions. We could not measure the actual
offset of the drive during the May experiment, as the drive failed and had to be removed from
the beam line, but the true offset of the drive during the August experiment was measured
and compared with the offset value that was used during the analysis.
The L110 calibration constants used in this analysis are doffset = −59.4233 mm and
dref = 29.731 mm from Equation 4.12. The L092 calibration constants for the May experiment
(1.3 MeV/u) are doffset = −59.4233 mm and dref = 33.5 mm. The L092 calibration constants
for the August experiment (0.6 MeV/u) are doffset = −56.83 mm and dref = 29.731 mm. This
analysis was performed with only the May and August data. The procedure for producing
the doffset values is equivalent to the method for determining δ in Section 4.2.2.
It appears that in May, the L092 slit drive may not have been moving along an exact 45◦
trajectory. Another possibility is that it may have been offset in the direction perpendicular
to the drive motion. During the measurement in May, the drive would continuously slip out
of calibration. When it was finally removed from the beam line in August, it was observed
that the drive had completely stripped the threads of the bolt holding it in place. This would
be consistent with the drive not functioning properly, as we suspected.
In order to match the data to theory, the dref for L092 in May had to be increased from
its original value, which was identical to the L110 calibration. We can interpret this to mean
that the beam intercepted the slits along a path that was longer than the path along the
central 45◦ angle, either being due to a change in the angle of the drive or a misalignment
along the perpendicular 45◦ direction. The original doffset value was able to be used to match
the data to theory, which indicates that the drive’s zero position was properly aligned at the
5 mm hole, with no overall offset along the 45◦ axis.
For the August data, the drive had to be installed quickly without being properly aligned,
122

as it was needed immediately. This was acceptable since we knew any misalignment that
might occur would be along the 45◦ axis, which could then be compensated for in the
post-processing of the data. As it turns out, the doffset value had changed by ∼3 mm. This
corresponds to an offset along the 45◦ axis, as predicted would occur since it was aligned by
eye. The dref constant returned to its original value, implying that this drive was properly
functioning along the 45◦ axis, and intercepting the beam as expected.
The interpretation of the May measurements with respect to how well it agrees with
theory predictions needs to take into consideration that these offsets were included in the
analysis of its L092 data. Even though the reasons for these offsets as described in this
section are well understood, it was determined that the best approach would be to redo this
energy measurement and avoid using this data set in further analysis, which was the course
of action taken.

123

BIBLIOGRAPHY

124

BIBLIOGRAPHY
[1] LT Sun, J Brandon, D Cole, M Doleans, D Morris, G Machicoane, E Pozdeyev, T Ropponen, and L Tobos. Intense beam production with SuSI. ECRIS, 10:4, 2010.
[2] PA Z´avodszky, B Arend, D Cole, J DeKamp, G Machicoane, F Marti, P Miller, J Moskalik,
J Ottarson, J Vincent, and et al. Design of SuSI–superconducting source for ions at
NSCL/MSU–II. the conventional parts. Nuclear Instruments and Methods in Physics
Research Section B: Beam Interactions with Materials and Atoms, 241(1):959–964, 2005.
[3] G Machicoane, C Benatti, D Cole, M Doleans, OK Kester, F Marti, X Wu, PA Zavodszky,
and C Zhang. Performance investigation of the nscl 18 ghz superconducting ecr ion
source susi. In in Proceedings of the 2009 Particle Accelerator Conference (Vancouver,
Canada 2009), 2009.
[4] RC York, HG Blosser, TL Grimm, F Marti, J Vincent, X Wu, and A Zeller. The
NSCL coupled cyclotron project-overview and status. In Proceedings of the Fifteenth
International Conference on Cyclotrons and their Applications, pages 687–691, 1998.
[5] DJ Morrissey. The coupled cyclotron project at the NSCL. Nuclear Physics A, 616(12):45–55, 1997. Radioactive Nuclear Beams.
[6] The ISOLDE radioactive ion beam facility. http://isolde.web.cern.ch/, December 2011.
[7] The NSCL at MSU. http://nscl.msu.edu/, March 2014.
[8] DJ Morrissey and BM Sherrill. Radioactive nuclear beam facilities based on projectile
fragmentation. Philosophical Transactions A, 356(1744):1985–2006, 1998.
[9] DJ Morrissey, BM Sherrill, M Steiner, A Stolz, and I Wiedenhoever. Commissioning
the A1900 projectile fragment separator. Nuclear Instruments and Methods in Physics
Research Section B: Beam Interactions with Materials and Atoms, 204(0):90–96, 2003.
14th International Conference on Electromagnetic Isotope Separators and Techniques
Related to their Applications.
[10] J Beringer and et al. (Particle Data Group). Review of Particle Physics. Journal of
Physics G, 33:258, 2006.

125

[11] K Purtell, K Rethman, A Haagsma, J Finck, J Smith, and J Snyder. Performance
comparison of MoNA and LISA neutron detectors. Bulletin of the American Physical
Society, 55, 2010.
[12] MS Wallace, MA Famiano, MJ van Goethem, AM Rogers, WG Lynch, J Clifford,
F Delaunay, J Lee, S Labostov, M Mocko, and et al. The high resolution array (HiRA)
for rare isotope beam experiments. Nuclear Instruments and Methods in Physics Research
Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 583(2):302–
312, 2007.
[13] D Bazin, JA Caggiano, BM Sherrill, J Yurkon, and A Zeller. The S800 spectrograph.
Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions
with Materials and Atoms, 204:629–633, 2003.
[14] S Schwarz, G Bollen, D Lawton, P Lofy, DJ Morrissey, J Ottarson, R Ringle, P Schury,
T Sun, V Varentsov, and et al. The low-energy-beam and ion-trap facility at NSCL/MSU.
Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions
with Materials and Atoms, 204:507–511, 2003.
[15] K Minamisono, BR Barquest, G Bollen, PF Mantica, DJ Morrissey, R Ringle, and
S Schwarz. Beam cooling and laser spectroscopy (BECOLA) project at NSCL. Bulletin
of the American Physical Society, 54, 2009.
[16] G Bollen, DJ Morrissey, and S Schwarz. A study of gas-stopping of intense energetic
rare isotope beams. Nuclear Instruments and Methods in Physics Research Section A:
Accelerators, Spectrometers, Detectors and Associated Equipment, 550(1):27–38, 2005.
[17] G Bollen, C Campbell, S Chouhan, C Guenaut, D Lawton, F Marti, DJ Morrissey,
J Ottarson, G Pang, S Schwarz, and et al. Manipulation of rare isotope beams–from
high to low energies. Nuclear Instruments and Methods in Physics Research Section B:
Beam Interactions with Materials and Atoms, 266(19):4442–4448, 2008.
[18] X Wu, G Bollen, M Doleans, TL Grimm, W Hartung, F Marti, S Schwarz, RC York,
and Q Zhao. MSU re-accelerator-the re-acceleration of low energy RIBs at the NSCL.
In Proceedings of the 13th Workshop on RF Superconductivity, Beijing, China, page 28,
2007.
[19] D Leitner, J Popielarski, A Lapierre, F Montes, G Perdikakis, S Schwarz, W Wittmer,
and X Wu. Status of the ReAccelerator facility ReA for rare isotopes. Bulletin of the
American Physical Society, 56, 2011.

126

[20] M Berz and K Makino. COSY INFINITY 9.0. MSU Report MSUHEP-060803, pages
1–62, 2006.
[21] J Qiang, RD Ryne, S Habib, and V Decyk. An object-oriented parallel particle-in-cell
code for beam dynamics simulation in linear accelerators. Journal of Computational
Physics, 163(2):434–451, 2000.
[22] S Schwarz, G Bollen, M Doleans, O Kester, K Kittimanapun, M Kostin, A Lapierre,
M Portillo, C Wilson, JRC L´opez-Urrutia, and J Dilling. Status of the EBIT in the
ReA3 reaccelerator at NSCL. Hyperfine Interactions, 199(1-3):261–267, 2011.
[23] M Portillo, G Bollen, S Chouhan, O Kester, G Machicoane, J Ottarson, S Schwarz, and
A Zeller. An achromatic mass separator design for ions from the EBIT charge breeder
at the NSCL. PAC2009, Vancouver, May, 2009.
[24] Oliver Kester, D Bazin, C Benatti, J Bierwagen, G Bollen, S Bricker, S Chouhan,
C Compton, K Davidson, J DeLauter, and et al. The MSU/NSCL re-accelerator ReA3.
In International Conference on Radio Frequence Superconductivity, pages 57–61, 2009.
[25] J Schmidt, J Maus, N Mueller, A Schempp, O Kester, and J Haeuser. Tuning of the
4-Rod-RFQ for MSU. 2010.
[26] W Hartung, J Bierwagen, S Bricker, C Compton, J DeLauter, P Glennon, M Hodek,
M Johnson, F Marti, P Miller, and et al. Superconducting quarter-wave resonator cavity
and cryomodule development for a heavy ion re-accelerator. In Proceedings of XXIV
International Linear Accelerator Conference, Victoria, BC, 2008.
[27] O Kamigaito. World-wide efforts on rare isotope and radioactive ion beams. volume 10,
2010.
[28] X Wu, B Arend, C Compton, A Facco, M Johnson, D Lawton, D Leitner, F Montes,
S Nash, J Ottarson, and et al. The design and commissioning of the accelerator system
of the rare isotope reaccelerator–rea3 at michigan state university. High-Intensity
High-Brightness Hadron Beams (HB), Beijing, China, 2012.
[29] Conceptual Design Report for the Facility for Rare Isotope Beams (FRIB) at Michigan
State University, 2010.
[30] R York, S Assadi, G Bollen, T Glasmacher, W Hartung, MJ Johnson, G Machicoane,
F Marti, E Pozdeyev, M Syphers, and et al. Technical challenges in design and construction of FRIB. FROBN2, PAC, 11:2561, 2011.
127

[31] T Pelaia II. Open XAL status report 2013. 2012.
[32] J Qiang. IMPACT-T user document version 1.6. Lawrence Berkeley National Laboratory:
LBNL Paper LBNL-62326 http://www.escholarship.org/uc/item/2sw7g9c3, 2010.
[33] M Berz and K Makino. COSY INFINITY version 8.1. user’s guide and reference manual.
Department of Physics and Astronomy MSUHEP-20704, Michigan State University,
2002.
[34] DA Edwards and MJ Syphers. Transverse Linear Motion, pages 57–98. Wiley-VCH
Verlag GmbH, 2008.
[35] LN Hand and JD Finch. Analytical Mechanics. Cambridge University Press, 1998.
[36] ED Courant and HS Snyder. Theory of the alternating-gradient synchrotron. Annals of
Physics, 3(1):1–48, 1958.
[37] DC Carey. The Optics of Charged Particle Beams. Second Printing with Corrections.
Harwood Academic Publishers, Reading, UK, 1992.
[38] H Wiedemann. Particle Accelerator Physics. Springer, 2007.
[39] MG Minty and F Zimmermann. Longitudinal phase space manipulation. In Measurement
and Control of Charged Particle Beams, pages 175–209. Springer, 2003.
[40] H Goldstein, C Poole, and J Safko. Classical mechanics. Classical mechanics (3rd ed.),
1, 2002.
[41] Open XAL. http://xaldev.sourceforge.net/, 2014.
[42] Eduard Pozdeyev. Private Communication.
[43] LD Landau and EM Lifshitz. Electrodynamics of continuous media. Elsevier et al.,
Amsterdam et al., 2008.
[44] CST-computer simulation technology, 2013.
[45] Winter Zheng. Private Communication.

128

[46] TP Wangler. Principles of RF linear accelerators. Wiley Beam Phys. Accel. Technol.
Wiley, New York, NY, 1998.
[47] D Leitner, ML Galloway, TJ Loew, CM Lyneis, IC Rodriguez, and DS Todd. High
intensity production of high and medium charge state uranium and other heavy ion
beams with venusa). Review of Scientific Instruments, 79(2):02C710, 2008.
[48] G Machicoane, D Cole, D Leitner, D Neben, and L Tobos. Metallic beam development
for the facility for rare isotope beama). Review of Scientific Instruments, 85(2):02A957,
2014.
[49] SA Kondrashev, A Barcikowski, B Mustapha, PN Ostroumov, RH Scott, SI Sharamentov,
and NV Vinogradov. Experimental results on multi-charge-state LEBT approach. In
Particle Accelerator Conference, 2007. PAC. IEEE, pages 1658–1660. IEEE, 2007.
[50] PN Ostroumov, A Barcikowski, SA Kondrashev, RH B Mustapha, NE SI Sharamentov,
and et al. Extraction from ecr and recombination of multiple-multiplecharge state
heavy-ion beams in lebt. LINAC, 2008.
[51] J Tamura, T Hattori, N Hayashizaki, T Ishibashi, T Ito, T Kanesue, H Kashiwagi, and
M Okamura. Multiple charge state ion beam acceleration with an RFQ LINAC. Venice,
ITALY, 2009. HIAT09.
[52] PN Ostroumov and KW Shepard. Multiple-charge beam dynamics in an ion linac.
Physical review special topics-accelerators and beams, 3(3):030101, 2000.
[53] C Allen. Theory and technique of beam envelope simulation: Simulation of bunched
particle beams with ellipsoidal symmetry and linear space charge forces. LANL Report
No. LA-UR-02-4979, August 2002.
[54] A Facco and V Zvyagintsev. Beam steering in superconducting quarter-wave resonators: An analytical approach. Physical Review Special Topics Accelerators and
Beams, 14(7):070101, July 2011.
[55] PN Ostroumov and KW Shepard. Correction of beam-steering effects in low-velocity
superconducting quarter-wave cavities. Physical Review Special Topics-Accelerators and
Beams, 4(11):110101, 2001.

129