RESONANCE COMPENSATION STUDIES AT THE FNAL RECYCLER RING

By

Cristhian Gonzalez-Ortiz

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

Physics—Doctor of Philosophy

2024

ABSTRACT

Circular accelerators, such as the ones in the Fermilab Accelerator Complex, are indispensable tools

driving scientific discovery by propelling particles to high energies and high intensities. In the

context of the Proton Improvement Plan II (PIP-II) at Fermilab, the Recycler Ring (RR) is confronted

with the challenge of high-intensity operation. The space charge tune shifts at these intensities

necessitate the mitigation of third-order resonances. By measuring the Resonance Driving Terms

(RDTs) for each resonance line, their strengths and damaging effects can be quantified.

The study is focused on the utilization of dedicated normal and skew sextupoles for the com-

pensation of these resonance lines. By employing the response matrix method to measure the

Resonance Driving Terms (RDTs) relative to the currents in these sextupoles, a methodical ap-

proach to simultaneously compensate individual and multiple third-order resonances has been

developed and assessed. The assessment and verification incorporate mapping out beam losses

through dynamic and static tune scans, along with the utilization of the Ion Profile Monitor (IPM)

system for beam size measurements.

Furthermore, optimization algorithms have been explored and incorporated into the resonance

compensation scheme. Parallel, to the work done at Fermilab, the study extends its scope by

performing related experiments at the CERN Proton Synchrotron Booster (PSB). For these experi-

ments, multiple resonance lines in the tune diagram were compensated by optimizing the currents

in the compensation elements with the aid of advanced optimization algorithms.

High-intensity operation of the Recycler Ring involves understanding the interplay between the

space charge tune spread, resonance lines, and specialized subsystems such as transverse dampers.

The following study also delves into experiments that try to discern each of these phenomena, while

highlighting challenges still ahead. Overall, these studies indicate the feasibility of compensating

for multiple third-order resonances, contributing to the objectives of PIP-II for the reliable delivery

of a 1.2 MW proton beam to experimental facilities.

ACKNOWLEDGEMENTS

I would like to acknowledge Dr. Robert Ainsworth for taking those couple of nutmegs when we

played footy together. Thank you, Rob, for making me a better scientist and teaching me so many

valuable life lessons. Thank you for always having my back.

I would like to acknowledge the instrumental role of Dr. Peter Ostroumov in these years of my

Ph.D. at Michigan State University. Thank you, Peter, for pushing me to be a better physicist and

scientist, every day. Thank you for believing in me, even in the moments I didn’t believe in myself.

Thank you to the members of my doctoral committee, Steve Lund, Yue Hao, Paul Gueye

and Kendall Mahn. I truly value every comment and every word of encouragement at all of our

meetings. We did not meet frequently, but your guidance was crucial through this process. Also

thank you to whole support network at MSU and FRIB—these are institutions where you will never

feel alone and always have someone support you when needed.

I would like to acknowledge the whole MI/RR department at Fermilab. Thank you, Kyle, Adam,

Nino, Osama, Betiay, Denton, Marty, Mike, Phil, MaryKate, Matilda and Meiqin, and whoever

I might have left out. Thank you for having the patience and your unconditional support while I

pushed through with my experiments. You guys always had a smile even when I asked the most

basic questions. Special acknowledgement to the operations team at Fermilab for always being

open to questions and discussion. Thank you for having the patience to withstand me when I would

continuously trip the permit on beam losses. I would also like to acknowledge the hospitality and

awesomeness of our CERN collaborators, specially Foteini and Tirsi.

Luisa. Tú has sido el soporte y la inspiración estos últimos años para poder estar a puertas de

cruzar la meta. Gracias por siempre estar ahí y no rendirte nunca con lo nuestro. Vamos a lograr

grandes cosas juntos. Ya casi adoptamos un perrito, manguitos szah.

Gracias a toda mi familia, incluyendo a esos amigos que considero mi familia, por apoyarme en

este camino. Desde Bogotá hasta llegar a Chicago pasando por Lansing. Nunca pensé en llegar tan

lejos de la loma. No hubiera podido terminar este doctorado sin el apoyo de todos ustedes. Gracias

por estar ahí en cada momento de mi vida. Thank you, Sjuju, the great celestial whale overseer.

iii

LIST OF ABBREVIATIONS .

.

.

.

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

vi

TABLE OF CONTENTS

CHAPTER 1

INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . .

.
1.1 Circular Accelerators and Storage Rings . . . . . . . . . . . . . . . . . . . . .
.
1.2 Fermilab .
.
.
1.3 Outline

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

. .
.
.

. .
.
.

.
.

.
.

.
.

1
2
3
6

CHAPTER 2

7
BEAM DYNAMICS IN RINGS . . . . . . . . . . . . . . . . . . . . . .
2.1
7
Introductory Accelerator Physics . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Lie Maps in Accelerator Physics . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.3 One-turn Map and Normal Form . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.4 Resonances in Circular Accelerators . . . . . . . . . . . . . . . . . . . . . . . 18
2.5 Resonance Driving Terms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.6 Amplitude-Dependent Tune Shift . . . . . . . . . . . . . . . . . . . . . . . . . 24
. 28
2.7 Space Charge Tune Shift

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

CHAPTER 3

. .

THE FNAL RECYCLER RING . . . . . . . . . . . . . . . . . . . . . . 35
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.1
.
Introduction . .
3.2 General Description .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.3 Tune Diagram and Resonances . . . . . . . . . . . . . . . . . . . . . . . . . . 42
. . . . . . . . . . . . . . . . . . . . . . . . 45
3.4 High Intensity and Tune Footprint
3.5 Diagnostic Devices .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
3.6 Sextupoles for Resonance Compensation . . . . . . . . . . . . . . . . . . . . . 51

.

CHAPTER 4

COMPENSATION OF THIRD-ORDER RESONANCES AT LOW
INTENSITIES .

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

. 54
4.1 Global RDTs and Lattice Model
. . . . . . . . . . . . . . . . . . . . . . . . . 54
. 57
4.2 Measurement of Third Order RDTs . . . . . . . . . . . . . . . . . . . . . . .
4.3 Compensation of RDTs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
. 86
4.4 Optimization of Compensation Currents . . . . . . . . . . . . . . . . . . . .
. 89
4.5 Experimental Verification of Compensation . . . . . . . . . . . . . . . . . .
4.6 Additional Sextupoles for Resonance Compensation . . . . . . . . . . . . . . . 100

CHAPTER 5

RESONANCE COMPENSATION STUDIES AT THE CERN
PROTON SYNCHROTRON BOOSTER . . . . . . . . . . . . . . . . . 108
5.1 General Description .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
5.2 Tune Diagram and Operation . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
. 113
5.3 Optimization Algorithms for Resonance Compensation . . . . . . . . . . . .
. 123
5.4 Experimental Verification of Compensation . . . . . . . . . . . . . . . . . .

CHAPTER 6

. . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.1 Space Charge Tune Shift
6.2 Measurement of Tune Spread . . . . . . . . . . . . . . . . . . . . . . . . . .
Intensity-Dependent Effects and Non-Gaussian Beam Profiles . . . . . . . . .
6.3

HIGH INTENSITY STUDIES . . . . . . . . . . . . . . . . . . . . . . . 125
. 125
. 126
. 129

iv

6.4 Static Tune Scans at Different Intensities . . . . . . . . . . . . . . . . . . . . . 131
6.5 Effect of Transverse Dampers . . . . . . . . . . . . . . . . . . . . . . . . . . . 139

CHAPTER 7

7.1 Conclusions . .
7.2 Future Work .

CONCLUSIONS AND FUTURE WORK . . . . . . . . . . . . . . . .

. 142
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142
. 144
.

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

. .
.

. .

BIBLIOGRAPHY . .

. .

.

.

.

.

.
.

.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148

v

LIST OF ABBREVIATIONS

Michigan State University

Fermilab National Accelerator Laboratory

Recycler Ring

Main Injector

Resonance Driving Terms

Neutrinos at the Main Injector

CERN Proton Synchrotron Booster

MSU

FNAL

RR

MI

RDTs

NuMI

PSB

ACNET Accelerator Control Network

TbT

NAFF

GFCs

Turn-by-Turn Data

Numerical Analysis of Fundamental Frequencies

Generating Function Coefficients

vi

CHAPTER 1

INTRODUCTION

Particle accelerators are the workhorses for modern scientific discoveries. Experimental nuclear and

particle physics research benefits greatly from the progress of accelerator physics and technology.

Accelerator physics is a rich field of applied physics living on the intersection of electromagnetism,

solid-state and atomic physics, nonlinear mechanics, plasma physics, and quantum mechanics, just

to name a few [1]. Furthermore, the design and operation of modern accelerator projects require

costly enterprises of scientists, engineers, operators, and politicians coming together under one

metaphorical roof. Everyone coming together to perform “megascience” [2].

The scientific principle of particle accelerators involves accelerating, steering, and/or storing

charged particles through electromagnetic manipulations. These manipulations occur through a

plethora of devices and components that can control electromagnetic fields, e.g., magnets and

electrical cavities. The group of particles subject to this electromagnetic handling is called "the

beam". The field of beam dynamics studies the interaction between the beam and the steering

devices, as well as the Coulomb interactions between the beam itself—this is known as space

charge physics. One can make an additional distinction when these steering devices are configured

circularly or linearly. This is the distinction between circular accelerators and linear accelerators

(Linacs).

Furthermore, one can categorize particle accelerators by the type of elementary particles com-

posing the beam and how close to the speed of light they travel. The first category refers to the

distinction between hadrons and leptons—particles that interact or do not interact through the

strong force, respectively [3]. For example, protons and heavy ions are considered hadrons, while

electrons and muons are considered leptons. The second category refers to whether particles in the

machine travel at high or low energy. An example of a low-energy hadron machine is the heavy ion

Linac at FRIB (Facility for Rare Isotope Beams) [4]. An example of a high-energy lepton machine

was the Stanford Linear Accelerator located at SLAC National Accelerator Laboratory [5]. The

two most famous high-energy hadron machines in history are the Tevatron [6], which operated

1

at Fermilab, and the LHC (Large Hadron Collider), operating at CERN [7]. Furthermore, there

are accelerator projects that encompass several categories, such as the future EIC (Electron-Ion

Collider) being built at Brookhaven National Laboratory [8], which will use an electron ring—a

lepton machine—and a heavy-ion circular accelerator—a hadron machine—to probe new physics.

A plethora of accelerator projects worldwide are either operational, under commissioning, or being

designed.

The following thesis will explore the beam dynamics of a circular machine. The research results

will be applied to the Fermilab Recycler Ring, which stores high-energy protons.

1.1 Circular Accelerators and Storage Rings

As will become more apparent in Ch. 2, a particle accelerator can be thought of as a composition

of accelerator-themed LEGO® bricks [9]. Each elemental LEGO® brick can be thought of as an

accelerator component performing some mapping on the charged particles entering it. As it turns

out, one can assemble these LEGO® bricks circularly to give rise to circular accelerators. The

assembly of these blocks in a particular shape gives rise to what is known as the lattice of the

accelerator.

The acceleration part of these structures comes from elements inside the lattice that introduce

some electromotive force in the longitudinal direction. The most common example for these blocks

is radio-frequency (RF) cavities [1], with super-conducting RF cavities also as an established

technology [10]. Particles that go through these elements gain energy on every pass. In the case

where there are no acceleration blocks on these structures, a storage ring arises. Nevertheless,

storage rings can also have RF cavities just for longitudinal beam manipulation but no overall

acceleration—such is the case of the Recycler Ring.

Circular accelerators are unique since particles have to pass thousands or even millions of turns

through the same LEGO® blocks. This multipass feature gives birth to exciting and complex

dynamics inside these machines. One of these phenomena is called betatron resonances. The most

simple lattice of a high-energy machine is composed of focusing and steering blocks, which are

dipole and quadrupole magnets, i.e., the lowest-order multipole magnets. For reasons that will

2

become apparent in Ch. 2, these elements, in combination with free drift spaces, represent linear

blocks. For the simplest circular machine, these linear LEGO® bricks are assembled to create

a linear lattice. One can design a linear lattice with stable particle orbits around the accelerator.

Nevertheless, accounted or unaccounted elements, described by linear or nonlinear blocks, around

the machine can perturb the stable orbits. Ultimately, the effect of these perturbations can add up

coherently over many turns to push the beam out of the acceptance of the lattice, i.e., the particles

hit the enclosing vacuum pipe and create beam loss. This whole process is known as a betatron

resonance in a circular accelerator. Chapter 2 provides a mathematical description of this process.

The following thesis describes an effort to mitigate the deleterious effect of these resonances in

the Recycler Ring. After dipole and quadrupole, the third order of multipole magnetic fields is the

sextupole component. Therefore, sextupole fields around the lattice are the source of third-order

betatron resonances. Specifically, this thesis explores mitigation techniques to these third-order

resonances, mainly in the Fermilab Recycler Ring (see Ch. 3, Ch. 4 and Ch. 6), but also with some

experiments done at the CERN Proton Synchrotron Booster (see Ch. 5).

1.2 Fermilab

The best introduction to Fermilab is to cite an excerpt from Ref. [2]:

[...] A passenger peers through the window of an airplane. As his plane flies into

Chicago’s O’Hare Field from the west, he notices a large ring on the ground below

(see Fig. 1.1). Near it, he sees a towering white structure, a group of colorful smaller

buildings, an expanse of forest, open fields, and lakes.

"What is that ring?" he asks his neighbor.

"Fermilab," she replies. "It’s a physics laboratory. The government supports research

there into what the universe is made of."

"Why the ring?"

"It’s the four-mile-round main ring of a machine called the Tevatron. It turns protons

into tools for looking inside the atomic nucleus. Huge magnets steer the protons

3

around the ring, while high voltages accelerate them. [...]" (pp. 1)

Figure 1.1 Aerial view of the Fermi National Accelerator Laboratory (FNAL) located in Batavia,
IL, USA [11].

The Fermi National Accelerator Laboratory (FNAL), better known as Fermilab, has a long

and rich history of designing, building, and operating high-energy particle accelerators. Ever

since the founding director of Fermilab, Robert R. Wilson, envisioned the 400 GeV Main Ring

back in 1967, Fermilab has been at the forefront of accelerator physics [2, 6, 12]. The most

famous accelerator project hosted by Fermilab has been the Tevatron, a proton-antiproton circular

collider with a circumference of around 6.28 km. This machine took protons and antiprotons from

smaller machines still in operation or repurposed as of 2024, e.g., the Recycler Ring. The Tevatron

operated up until 2011, leaving an indelible legacy in the field of high energy and accelerator

physics. Nostalgia aside, Fermilab still hosts a deluge of particle physics experiments connected to

its main accelerator complex.

4

Figure 1.2 summarizes the current layout of the Fermilab Accelerator Complex. As of 2024, the

Fermilab Accelerator Complex is composed of an 𝐻− source that connects to a linear accelerator,

accelerating the ions to an energy of 400 MeV. This linear accelerator feeds to the first circular

machine—the Booster—where protons are achieved and accelerated to an energy of 8 GeV. After

the Booster, the protons are transported to the Recycler Ring (RR), which is the second circular

machine. The Recycler Ring stores and stacks protons to increase the beam intensity delivered

to the Main Injector (MI). This last circular accelerator is where protons are accelerated from an

energy of 8 GeV to 120 GeV. Once at this energy, the protons are transported to the Neutrinos

at the Main Injector (NuMI) experiment to create the world’s most intense neutrino beam [13].

Nevertheless, throughout the chain of accelerators, the beam is also delivered to many other

experiments. Therefore, the facility has several modes of operation depending on the online

experiments. A more detailed and technical study of the current Fermilab Accelerator Complex,

focusing on the Recycler Ring, is given in Ch. 3.

Figure 1.2 Schematic layout of the Fermilab Accelerator Complex as of 2024. Original plot
provided by R. Ainsworth.

5

1.3 Outline

The following thesis will explore the compensation of third-order resonances in the Fermilab

Recycler Ring. Chapter 1 introduces the motivation behind this thesis work. Chapter 2 summarizes

single particle dynamics with exponential Lie operators’ help and introduces a relevant concept of

collective beam dynamics: the space charge tune shift. This theoretical overview gives a segue

into this thesis’ Ch. 3, where the Recycler Ring is introduced and described in detail. This chapter

motivates the compensation of third-order resonances under the framework of the current and

future operation of the RR. With the basic physics concepts and the description of the machine

put in place, Ch. 4 describes in full detail the scheme and experiments developed to compensate

for third-order resonances at low intensities. Before moving to explore the Recycler Ring at high

intensities, Ch. 5 provides an interlude to show a series of experiments done at the CERN PS

Booster. These experiments explore the use of advanced optimization algorithms to compensate

for multiple resonance lines simultaneously. Coming back to Fermilab, Ch. 6 showcases the

studies and experiments done at high intensities in the RR to understand the interplay between the

compensation of resonance lines and space charge effects. Finally, Ch. 7 brings down the curtain

by providing some general conclusions and future work stemming from this thesis.

6

CHAPTER 2

BEAM DYNAMICS IN RINGS

2.1

Introductory Accelerator Physics

The basic building blocks of a particle accelerator are the elements that steer and focus the

beam in the transverse direction. This is done by utilizing and manipulating the Lorentz force (cid:174)𝐹𝐿

by means of the electromagnetic fields (cid:174)𝐸 and (cid:174)𝐵 to act on some particle with charge 𝑞 and velocity
(cid:17). While a handful of electromagnetic devices can do this, the most

(cid:16) (cid:174)𝐸 + (cid:174)𝑣 × (cid:174)𝐵

(cid:174)𝑣, i.e., (cid:174)𝐹𝐿 = 𝑞

prominent ones in high-energy accelerators are magnets with no electric field, (cid:174)𝐸 = 0. In particular,

there are pure dipole magnets for steering and pure quadrupole magnets for focusing. Nevertheless,

some machines such as the Recycler Ring have combined function magnets with both types of

magnets—and even higher order magnets—embedded in one to steer and focus the particles. The

previous information assumes that every magnet can be expanded and described as a decomposition

of magnetic multipoles, where dipole and quadrupole components are the lowest order terms of the

expansion. Therefore, using the Beth representation [1], the multipole expansion for an arbitrary

multipole magnet reads:

1
[𝐵𝜌]

(cid:0)𝐵𝑦 (𝑥, 𝑦) + 𝑖𝐵𝑥 (𝑥, 𝑦)(cid:1) = −

∞
∑︁

(𝑏𝑛 + 𝑖𝑎𝑛) (𝑥 + 𝑖𝑦)𝑛 ,

1
𝜌

𝑛=0
where 𝑏𝑛 and 𝑎𝑛 are the multipole coefficients defined by

𝑏𝑛 =

1
𝐵0𝑛!

𝜕𝑛𝐵𝑦
𝜕𝑥𝑛

(cid:12)
(cid:12)
(cid:12)
(cid:12)𝑥=𝑦=0

,

𝑎𝑛 =

1
𝐵0𝑛!

𝜕𝑛𝐵𝑥
𝜕𝑥𝑛

(cid:12)
(cid:12)
(cid:12)
(cid:12)𝑥=𝑦=0

.

(2.1)

(2.2)

For Eqs. 2.1 and 2.2, 𝑥 and 𝑦 represent the Cartesian coordinates, the product [𝐵𝜌] represents

the magnetic rigidity of the beam with 𝐵0 being the main dipole field and 𝜌 is the bending radius.

At the same time, 𝐵𝑥 and 𝐵𝑦 are the transverse magnetic fields in the magnets. Specifically, the

coefficients 𝑏𝑛 and 𝑎𝑛 represent the multipole coefficient of the magnet with the dipole coefficient

defined as 𝑏0 = 1, such that 𝐵0𝑏0 = − [𝐵𝜌] /𝜌. Following the expansion, the term 𝑎0 is known as

the dipole roll coefficient, 𝑏1 as the quadrupole coefficient, 𝑎1 as the skew quadrupole coefficient,

𝑏2 as the sextupole coefficient, 𝑎2 for the skew sextupole coefficient, 𝑏3 for the octupole coefficient,

𝑎3 for the skew octupole coefficient, and so forth. This expansion follows the U.S. convention.

7

The most basic circular accelerator of circumference 𝐶 is composed of LEGO® blocks chosen

from a pile of dipoles, quadrupoles, and free drift spaces. Equations 2.1 and 2.2 describe each

of these elements, i.e., for the quadrupole case, the only non-zero coefficient is 𝑏1. The particles

inside this ring have a longitudinal relativistic velocity of 𝑣 = 𝛽𝐿𝑐, and therefore have a revolution

frequency of 𝑓𝑟𝑒𝑣 = 𝛽𝐿𝑐/𝐶. The Hamiltonian of a single particle with position coordinates 𝑥, 𝑦 and

momentum coordinates 𝑝𝑥, 𝑝𝑦 traversing through such a system at an independent time coordinate

𝑠 is:

𝐻 =

(cid:16)

1
2

𝐾𝑥 (𝑠)𝑥2 + 𝐾𝑦 (𝑠)𝑦2 + 𝑝2

𝑥 + 𝑝2
𝑦

(cid:17)

,

where 𝐾𝑥 (𝑠), 𝐾𝑦 (𝑠) are the effective focusing functions and are defined as:

𝐾𝑥 (𝑠) =

1
𝜌2 −

𝑏1(𝑠)
𝜌

, 𝐾𝑦 (𝑠) =

𝑏1(𝑠)
𝜌

(2.3)

(2.4)

assuming the definition of 𝑏1(𝑠) in 2.2 extends to describe the distribution of the horizontal

quadrupole coefficient around the ring. For the case where 𝜌 ≫ 1 (high-energy limit), the function

𝐾𝑦 (𝑠) = −𝐾𝑥 (𝑠), i.e., horizontally focusing quadrupoles will have a defocusing effect in the vertical

direction and vice versa.

In classic accelerator references, such as Refs. [1, 14, 15], the equations of motion derived from

the Hamiltonian in Eq. 2.3 are known as Hill’s equations. The usual accelerator-physics method

will define linear mappings bringing some initial state vector (cid:174)𝑋0 =

to solve this type of equation is to introduce transfer matrices for each type of linear element. This
(cid:17) to a final vector
(cid:17) using a symplectic matrix 𝑀, i.e., (cid:174)𝑋 𝑓 = 𝑀 (cid:174)𝑋0. Table 2.1 shows the 4D
transfer matrices for common linear elements found in accelerators. It is worth noting that these

𝑓 , 𝑦 𝑓 , 𝑦′
𝑓

0, 𝑦0, 𝑦′
0

(cid:174)𝑋 𝑓 =

𝑥 𝑓 , 𝑥′

𝑥0, 𝑥′

(cid:16)

(cid:16)

4D matrices can include coupling elements that couple the 𝑥 − 𝑦 plane. For this case, the off-block

coefficients of the matrices would be non-zero. An example of this last case is shown in Table 2.1

through the thin skew quadrupole case. Nevertheless, the starting point of this work is to consider

a circular accelerator built from linear non-coupling elements, while other coupling or nonlinear

elements are considered perturbative.

In Sec. 2.2, Lie operators extend these mappings to the

nonlinear regime.

8

Table 2.1 Transfer matrices for common accelerator elements in the high-energy regime and paraxial
approximation.

Element

Drift space
of length 𝐿

Dipole
of bending radius 𝜌,
length ℓ and
orbiting angle 𝜃 = ℓ/𝜌

Thin quadrupole
of focal length

𝑓 = lim
ℓ→0

1
𝑘ℓ

> 0

4D Transfer Matrix

𝑀 =

1 𝐿 0 0


0 1 0 0


0 0 1 L


0 0 0 1












𝑀 =

cos 𝜃
− 1
𝜌 sin 𝜃
0
0










𝜌 sin 𝜃 0 0
0 0
cos 𝜃
1 ℓ
0
0 1
0










𝑀 =

1
−1/ 𝑓
0
0










0
0
0
1
0
1
0 1/ 𝑓

0


0


0


1



Thick quadrupole
of strength 𝑘 > 0
and length ℓ

𝑀 =

√

𝑘ℓ

cos

√

𝑘ℓ

sin

1
√
𝑘

√

−

√

𝑘ℓ

𝑘 sin

√

𝑘ℓ

cos

0

0

0

0

0

0

0

0

cosh √︁|𝑘 |ℓ

sinh √︁|𝑘 |ℓ

1√
|𝑘 |

√︁|𝑘 | sinh √︁|𝑘 |ℓ

cosh √︁|𝑘 |ℓ



































Thin skew quadrupole
of focal length

𝑓𝑠 = lim
ℓ→0

1
𝑘 𝑠ℓ

> 0

𝑀 =

1
0
0

0
0
0


1 1/ 𝑓𝑠 0


0
1
0


1
0
1/ 𝑓𝑠 0












Like stacking LEGO® bricks, one can stack these transfer matrices to calculate the total mapping

through a stack of elements. The total linear mapping of a consecution of accelerator elements is

just the matrix multiplication of the corresponding transfer matrices, i.e., for a lattice with 𝑁 linear

9

blocks, the total transfer matrix reads:

𝑀𝑇 𝑜𝑡𝑎𝑙 = 𝑀𝑁 𝑀𝑁−1...𝑀1.

(2.5)

One can introduce the Courant-Snyder (CS) parametrization and Twiss parametrization to

understand linear lattices further. One can parametrize any linear element or stack of linear

elements functions by the Twiss parameters 𝛽𝑢 (𝑠), 𝛼𝑢 (𝑠) and 𝛾𝑢 (𝑠) and a phase advance defined
as 𝜙𝑢 (𝑠) = ∫ 𝑠

𝑑𝑠/𝛽𝑢 (𝑠), where 𝑢 stands for either the 𝑥 or 𝑦 plane. This results in the following

0

transfer matrix from location 𝑠 = 0 to an arbitrary location 𝑠:

𝑀 (𝑠) =

where 𝑀𝑢 is a square matrix defined as:

𝑀𝑥 (𝑠)

02×2

02×2 𝑀𝑦 (𝑠)









,









√︃ 𝛽𝑢 (𝑠)

𝛽𝑢 (0) (cos 𝜙𝑢 (𝑠) + 𝛼𝑢 (0) sin 𝜙𝑢 (𝑠))

𝛼𝑢 (0)−𝛼𝑢 (𝑠)
𝛽𝑢 (𝑠) 𝛽𝑢 (0)

cos 𝜙𝑢 (𝑠) −

1+𝛼𝑢 (0)𝛼𝑢 (𝑠)
𝛽𝑢 (𝑠) 𝛽𝑢 (0)

sin 𝜙𝑢 (𝑠)

√︁𝛽𝑢 (𝑠) 𝛽𝑢 (0) sin 𝜙𝑢 (𝑠)

√︃ 𝛽𝑢 (0)

𝛽𝑢 (𝑠) (cos 𝜙𝑢 (𝑠) − 𝛼𝑢 (𝑠) sin 𝜙𝑢 (𝑠))

𝑀𝑢 =









(2.6)

,








(2.7)

and 02×2 is a square matrix of size 2 × 2 with zeros in all of its entries. As a consequence of the

Twiss parametrization, one can express the position coordinates at an arbitrary 𝑠 along a linear

lattice as:

𝑢(𝑠) = √︁2𝐽𝑢 𝛽𝑢 (𝑠) cos (𝜙𝑢 (𝑠) + 𝜙𝑢0) ,

(2.8)

where 𝐽𝑢 is a constant of motion known as the action and 𝜙𝑢0 is an arbitrary phase constant.

Equation 2.8 describes what is known as betatron oscillations. Additionally, the action 𝐽𝑢 can be

calculated as

2𝐽𝑢 = 𝛾𝑢𝑢2 + 2𝛼𝑢𝑢𝑢′ + 𝛽𝑢𝑢′2,

(2.9)

for any 𝑠 using the property 𝛽𝑢𝛾𝑢 = 1+𝛼2

𝑢. Therefore, the particles traversing through the lattice will

exhibit oscillatory motion, with amplitude dictated by the 𝛽𝑢 (𝑠) function and frequency dictated

by the 𝜙𝑢 (𝑠) function.

Until now, one can apply this formalism to a linear or circular accelerator. No assumption

has been made on how the stacking of the LEGO® bricks in the lattice affects this formalism.

10

Nevertheless, in a circular accelerator of circumference 𝐶, there is an additional constraint on the

Twiss functions given that periodic boundary conditions exist, i.e., 𝛽(𝑠) = 𝛽(𝑠 + 𝐶), and similarly

for the other Twiss functions. Another critical parameter for a circular lattice is the tune 𝑄𝑢, defined

as the total phase advance over one turn of the machine, i.e.,

2𝜋𝑄𝑢 = 𝜙𝑢 (𝐶) =

∫ 𝑠+𝐶

𝑠

𝑑𝑠
𝛽𝑢 (𝑠)

.

(2.10)

With the definition of the tune 𝑄𝑢, Eq. 2.7 can be rewritten in order to calculate the one-turn

matrix 𝑀 (𝐶) of the circular accelerator, such that:

𝑀 (𝐶) =

𝑀𝑥 (𝐶)

02×2

02×2 𝑀𝑦 (𝐶)









,









where each plane will have its one-turn transfer matrix reading:

𝑀𝑢 (𝐶) =









cos 2𝜋𝑄𝑢 + 𝛼𝑢 (𝐶) sin 2𝜋𝑄𝑢

𝛽𝑢 (𝐶) sin 2𝜋𝑄𝑢

−𝛾𝑢 (𝐶) sin 2𝜋𝑄𝑢

cos 2𝜋𝑄𝑢 − 𝛼𝑢 (𝐶) sin 2𝜋𝑄𝑢

.









(2.11)

(2.12)

Ultimately, Eq. 2.12 can be used to calculate the particle’s state vector after 𝑁 turns, where the

total transfer matrix 𝑀𝑢 (𝑁𝐶) will be 𝑀𝑢 (𝑁𝐶) = 𝑀𝑢 (𝐶) 𝑁 . If one tracks a particle with some initial

conditions for enough turns and its geometrical coordinates (𝑥, 𝑥′, 𝑦, 𝑦′) are recorded through some

diagnostic device every turn, the resulting trajectory would lie on an ellipse—the phase space

ellipse. The astute reader would have already identified that Eq. 2.9 hinted at this fact, given that

this is just the implicit definition for an ellipse. In essence, the Twiss parametrization will give

the geometry to describe the ellipse fully. Figure 2.1 illustrates these statements by associating the

phase space ellipse geometry with the Twiss functions.

11

Figure 2.1 Phase space ellipse in geometrical coordinates with Twiss parametrization and its
counterpart transformation in normalized phase space.

Figure 2.1 also hints that one can do a linear transformation 𝑇𝑢 (𝑠) to transform the phase space

ellipse into a circle. This transformation is known as a Floquet transformation. Therefore, a

change of coordinates from geometrical coordinates (𝑢, 𝑢′) to normalized coordinates ( ˆ𝑢, ˆ𝑝𝑢) can

be achieved through the following linear transformation:









ˆ𝑢

ˆ𝑝𝑢









=

1
√
𝛽𝑢
𝛼𝑢√
𝛽𝑢









0
√

𝛽𝑢









𝑢





𝑢′












= √︁2𝐽𝑢

cos 𝜙𝑢 (𝑠)

− sin 𝜙𝑢 (𝑠)









.









Ultimately, one can express the one-turn Hamiltonian for circular accelerators as:

𝐻0 = 2𝜋𝑄𝑥 𝐽𝑥 + 2𝜋𝑄 𝑦𝐽𝑦,

(2.13)

(2.14)

which is simpler and more succinct than the Hamiltonian in Eq. 2.3. Conclusively, one can map

all the dynamics of a linear-circular accelerator with its intricacies to rotations on a simple circle.

The rotation matrix 𝑅(𝑠) summarizes the dynamics in normalized phase space, which will depend

on the lattice itself and is analogous to the transfer matrices 𝑀 (𝑠) in geometrical space. Therefore,

12

the following commutative diagram summarizes the linear dynamics:

𝑥, 𝑥′

𝑦, 𝑦′

(cid:169)
(cid:173)
(cid:173)
(cid:171)

(cid:170)
(cid:174)
(cid:174)
(cid:172)0

𝑇 (𝑠)

𝐽𝑥, 𝜙𝑥

𝐽𝑦, 𝜙𝑦

(cid:169)
(cid:173)
(cid:173)
(cid:171)

(cid:170)
(cid:174)
(cid:174)
(cid:172)0

𝑀 (𝑠)

𝑅(𝑠)

𝑥, 𝑥′

𝑦, 𝑦′

(cid:169)
(cid:173)
(cid:173)
(cid:171)

(cid:170)
(cid:174)
(cid:174)
(cid:172) 𝑓

𝑇 (𝑠)

𝐽𝑥, 𝜙𝑥

𝐽𝑦, 𝜙𝑦

(cid:169)
(cid:173)
(cid:173)
(cid:171)

.

(cid:170)
(cid:174)
(cid:174)
(cid:172) 𝑓

(2.15)

Lie operators are introduced in the following sections to generalize to non-linear mappings.

For a linear-circular accelerator, the last section has shown that starting from Hill’s equation, linear

transformations can be applied to a complex machine such as an accelerator in order to end up with

a simple mathematical equation such as the one described in Eq. 2.14. Put, nonlinear elements

will distort the phase space circle of Fig. 2.1, destroying the system’s linearity. The premise here is

that non-linear elements in circular accelerators are inevitable, and they will come from anywhere

and everywhere in the lattice, either from accounted or unaccounted sources. Therefore, one can

use higher-tier mathematical tools to describe non-linear dynamics in a circular accelerator. Linear

matrices will only get us so far.

2.2 Lie Maps in Accelerator Physics

The most basic element of a particle accelerator can be considered a LEGO® brick acting

as a black box transformation for a single particle. This black box takes some single charged
(cid:17), as defined in a Frenet-Serret coordinate
particle with initial transverse coordinates (cid:16)
(cid:17). For simplicity, this analysis
system, and maps them to some final coordinates (cid:16)

0, 𝑦0, 𝑦′
0

𝑥0, 𝑥′

𝑥 𝑓 , 𝑥′

𝑓 , 𝑦 𝑓 , 𝑦′
𝑓

will not consider any longitudinal effect, but it can be easily incorporated. By gathering the initial
(cid:17), and doing the same for the final coordinates,

coordinates into a vector, i.e., (cid:174)𝑋0 =

𝑥0, 𝑥′

(cid:16)

0, 𝑦0, 𝑦′
0

i.e., (cid:174)𝑋 𝑓 =

(cid:16)

𝑥 𝑓 , 𝑥′

𝑓 , 𝑦 𝑓 , 𝑦′
𝑓

(cid:17), one can define the mapping M that relates both vectors, such that:

(cid:174)𝑋 𝑓 = M (cid:174)𝑋0.

13

(2.16)

Unlike the previous section, M need not be a linear mapping. For a charged particle inside some

accelerator element (described using Hamiltonian dynamics), Poisson brackets and exponential Lie

operators can help describe this mapping M [9, 14, 16–18].

Let (cid:174)𝑋 = (𝑞1, 𝑝1, . . . , 𝑞𝑛, 𝑝𝑛) be a 2n dimensional vector, made from 𝑛 pairs of canonical

coordinates (𝑞𝑖, 𝑝𝑖) that make up the 2𝑛 dimensional phase space. Moreover, let two arbitrary
(cid:17) be functions of (cid:174)𝑋 and 𝑠, where 𝑠 plays the role of the independent

functions 𝑓

(cid:17) and 𝑔

(cid:16) (cid:174)𝑋; 𝑠

(cid:16) (cid:174)𝑋; 𝑠

"time" coordinate. The Poisson brackets [•, •] can be defined as:

[ 𝑓 , 𝑔] =

𝑛
∑︁

𝑖=1

𝜕 𝑓
𝜕𝑞𝑖

𝜕𝑔
𝜕 𝑝𝑖

−

𝜕 𝑓
𝜕 𝑝𝑖

𝜕𝑔
𝜕𝑞𝑖

.

(2.17)

Using this definition, one can explicitly write the Poisson bracket definition for a four-dimensional

phase space described by state vector (cid:174)𝑋 = (𝑥, 𝑥′, 𝑦, 𝑦′). This definition reads:

[ 𝑓 , 𝑔] =

𝜕 𝑓
𝜕𝑥

𝜕𝑔
𝜕𝑥′ −

𝜕 𝑓
𝜕𝑥′

𝜕𝑔
𝜕𝑥

+

𝜕 𝑓
𝜕𝑦

𝜕𝑔
𝜕𝑦′ −

𝜕 𝑓
𝜕𝑦′

𝜕𝑔
𝜕𝑦

.

(2.18)

The Lie operator : 𝑓 : acts on some function 𝑔 and is the adjoint operator of the Poisson bracket

operator. Its definition reads:

: 𝑓 : 𝑔 = [ 𝑓 , 𝑔] .

(2.19)

This specific : • : notation allows a compact notation to define the exponential Lie operator. The

exponential Lie operator of an arbitrary function 𝑓 is defined as

𝑒: 𝑓 :• =

∞
∑︁

𝑘=0

1
𝑘!

(: 𝑓 :) 𝑘 • .

(2.20)

For a Hamiltonian system, the mapping of coordinates from (cid:174)𝑋0 to (cid:174)𝑋 𝑓 follows the expression:

(cid:174)𝑋 𝑓 = 𝑒−ℓ:𝐻: (cid:174)𝑋

(cid:12)
(cid:12)
(cid:12)
(cid:12) (cid:174)𝑋= (cid:174)𝑋0

,

(2.21)

which is known as a Lie Map [16].

In this case, ℓ corresponds to the integration length of

the independent coordinate. For example, for a particle traversing a magnet with length 𝐿, the

14

integration length is ℓ = 𝐿. When looking at the one-turn map, the integration length corresponds

to the circumference 𝐶 of the accelerator over an effective Hamiltonian 𝐻𝑒 𝑓 𝑓 . Furthermore, the

integration length ℓ would be the phase advance 𝜇 if working with action-angle variables.

2.3 One-turn Map and Normal Form

LEGO® bricks can be stacked together to create complex structures. Analogously, one can

assemble accelerator elements to create complex ring-shaped structures such as circular acceler-

ators. In such structures, particles will experience the same one-turn mapping over thousands or

millions of turns. The one-turn map M1 of a circular accelerator is the composition (◦) of mappings

from every LEGO® element in the ring. Choosing an arbitrary initial point at 𝑠 = 0 and going

around the ring, the one-turn map describes the transformation of coordinates after one turn, i.e.,

(cid:174)𝑋𝑁=1 = M1 (cid:174)𝑋0. This map composition reads:

M1 = 𝑀𝑁+1 ◦ 𝑒:ℎ𝑁 : ◦ · · · ◦ 𝑒:ℎ2: ◦ 𝑀2 ◦ 𝑒:ℎ1: ◦ 𝑀1 = 𝑀𝑁+1𝑒:ℎ𝑁 : . . . 𝑒:ℎ2:𝑀2𝑒:ℎ1:𝑀1,

(2.22)

where 𝑀𝑖 is the matrix representation of a linear mapping that does not couple the 𝑥 − 𝑦 plane, e.g.,

drift space mapping or quadrupole mapping. On the other hand, the map 𝑒:ℎ𝑖: represents any linear

or non-linear mapping that can be found around the machine and can be considered a perturbation

to the ideal lattice, including coupling elements, e.g., skew quadrupoles, higher order multipole

elements. Figure 2.2 illustrates the procedure to build the one-turn map for a circular accelerator.

15

Figure 2.2 Diagram of an arbitrary circular accelerator to illustrate the one-turn map.

Through the use of the Baker-Campbell-Hausdorff formula [19], Eq. 2.22 can be collapsed to

the expression

M1 = 𝑒−𝐶:𝐻𝑒 𝑓 𝑓 :,

(2.23)

where 𝐶 is the circumference of the ring, and 𝐻𝑒 𝑓 𝑓 is the effective Hamiltonian of the machine over

one turn. As mentioned earlier, for most cases, it is of interest to look at the perturbations to the

linear uncoupled dynamics of the design lattice. With this in mind, Eq. 2.23 can be rewritten as:

M1 = 𝑒:ℎ:𝑅,

(2.24)

where 𝑅 is a rotation matrix encoding the linear uncoupled dynamics of the ideal lattice. On the

other hand, the term 𝑒:ℎ: encodes the perturbations to this ideal situation. It is worth pointing out

that for the case ℎ = 0, the traditional Courant-Snyder variables are recovered.

One can write the Courant-Snyder variables ( ˆ𝑥, ˆ𝑝𝑥, ˆ𝑦, ˆ𝑝𝑦) or normalized phase space coordinates

16

for a linear uncoupled case as:

ˆ𝑢 = √︁2𝐽𝑢 cos (cid:0)𝜙𝑢 + 𝜙𝑢0

(cid:1) ;

ˆ𝑝𝑢 = −√︁2𝐽𝑢 sin (cid:0)𝜙𝑢 + 𝜙𝑢0

(cid:1) ,

(2.25)

(2.26)

where 𝑢 can stand either for the 𝑥 or 𝑦 coordinate, 𝐽𝑢 and 𝜙𝑢 correspond to the action-angle variables

and 𝜙𝑢0 corresponds to the initial phase. For the case where perturbations exist, i.e., ℎ ≠ 0, the

action 𝐽𝑢 is not constant anymore and will be a function of 𝜙𝑢.

The Normal Form formalism is introduced at this point in order to find action-angle coordinates

𝐼𝑢 and 𝜓𝑢, such that the motion depends on 𝜓𝑢 at a constant 𝐼𝑢, with some initial phase 𝜓𝑢0. These

are known as non-linear action-angle variables. The variables 𝐼𝑢 and 𝜓𝑢 are calculated from the

transformation 𝑒−:𝐹: acting on 𝐽𝑢 and 𝜙𝑢. The whole point is to find these variables that allow the

Hamiltonian to be only amplitude-dependent. The following commutative diagram can summarize

these Normal Form gymnastics:

𝐽𝑥, 𝜙𝑥

𝐽𝑦, 𝜙𝑦

(cid:169)
(cid:173)
(cid:173)
(cid:171)

(cid:170)
(cid:174)
(cid:174)
(cid:172)0

𝑒−:𝐹:

𝐼𝑥, 𝜓𝑥

𝐼𝑦, 𝜓𝑦

(cid:169)
(cid:173)
(cid:173)
(cid:171)

(cid:170)
(cid:174)
(cid:174)
(cid:172)0

𝑒:ℎ ( 𝐽𝑢 , 𝜙𝑢 ):𝑅

𝑒:𝐻 (𝐼𝑢 ):

𝐽𝑥, 𝜙𝑥

𝐽𝑦, 𝜙𝑦

(cid:169)
(cid:173)
(cid:173)
(cid:171)

(cid:170)
(cid:174)
(cid:174)
(cid:172) 𝑓

𝑒−:𝐹:

𝐼𝑥, 𝜓𝑥

𝐼𝑦, 𝜓𝑦

(cid:169)
(cid:173)
(cid:173)
(cid:171)

.

(cid:170)
(cid:174)
(cid:174)
(cid:172) 𝑓

(2.27)

Without loss of generality, the generating function 𝐹 can be written as a Fourier expansion over

the objective space (𝐼𝑥, 𝜓𝑥, 𝐼𝑦, 𝜓𝑦) such that:

𝐹 =

∑︁

𝑗 𝑘𝑙𝑚

𝑓 𝑗 𝑘𝑙𝑚 (2𝐼𝑥)

𝑗+𝑘
2

(cid:0)2𝐼𝑦(cid:1) 𝑙+𝑚

2 𝑒𝑖[( 𝑗−𝑘)(𝜓𝑥+𝜓𝑥0)+(𝑙−𝑚)(𝜓𝑦+𝜓𝑦0)] .

(2.28)

Similarly, the argument of the Lie operator 𝑒:ℎ: from Eq. 2.24 can be expanded as:

ℎ =

∑︁

𝑗 𝑘𝑙𝑚

ℎ 𝑗 𝑘𝑙𝑚 (2𝐽𝑥)

𝑗+𝑘
2

(cid:0)2𝐽𝑦(cid:1) 𝑙+𝑚

2 𝑒𝑖[( 𝑗−𝑘)(𝜙𝑥+𝜙𝑥0)+(𝑙−𝑚)(𝜙𝑦+𝜙𝑦0)] .

(2.29)

17

For Eqs. 2.28 and 2.29, the integer indices 𝑗, 𝑘, 𝑙, 𝑚 run from 0 to infinity and correspond to the

four degrees of freedom for transverse phase space.

The terms 𝑓 𝑗 𝑘𝑙𝑚 are known as generating function coefficients. The terms ℎ 𝑗 𝑘𝑙𝑚 are known as

Hamiltonian coefficients or resonance driving terms (RDTs). Section 2.5 describes the non-linear

dynamics of accelerators using RDTs. The generating function coefficients 𝑓 𝑗 𝑘𝑙𝑚 can be related to

the Hamiltonian resonance driving terms ℎ 𝑗 𝑘𝑙𝑚 through the following relation [17, 20]:

𝑓 𝑗 𝑘𝑙𝑚 =

ℎ 𝑗 𝑘𝑙𝑚
1 − 𝑒2𝜋𝑖[( 𝑗−𝑘)𝑄 𝑥+(𝑙−𝑚)𝑄 𝑦]

,

(2.30)

where 𝑄𝑥 and 𝑄 𝑦 represent the transverse uncoupled and unperturbed tunes of the accelerator. The

transverse tunes of a circular accelerator are defined as the phase advances in each plane over one

turn, in units of 2𝜋, i.e., 𝑄𝑢 = 𝜙𝑢 (𝑠 = 𝐶)/2𝜋.

The terms ℎ 𝑗 𝑘𝑙𝑚 are generally defined by the order in which they enter the one-turn normal form

Hamiltonian [20]. With the assumption of thin elements around the ring, the general expression to

define RDTs reads:

ℎ 𝑗 𝑘𝑙𝑚 = Ξ 𝑗 𝑘𝑙𝑚

∑︁

𝑖

𝐿𝑖𝐾𝑛−1,𝑖 𝛽

𝑗+𝑘
2

𝑥,𝑖 𝛽

𝑙+𝑚
2

𝑦,𝑖 𝑒𝑖[( 𝑗−𝑘)𝜙𝑥,𝑖+(𝑙−𝑚)𝜙𝑦,𝑖] ,

where Ξ 𝑗 𝑘𝑙𝑚 is just a constant defined as:

Ξ 𝑗 𝑘𝑙𝑚 = −

1
2𝑛

1
𝑛!

(cid:18) 𝑛
𝑙 + 𝑚

(cid:19) (cid:18) 𝑗 + 𝑘
𝑗

(cid:19) (cid:18)𝑙 + 𝑚
𝑙

(cid:19)

.

(2.31)

(2.32)

For Eqs. 2.31 and 2.32, 𝑛 = 𝑗 + 𝑘 + 𝑙 + 𝑚 represents the order of the RDT. The sum over 𝑖 is done

over all multipoles of order 𝑛 and length 𝐿𝑖 that either have a normal component 𝐾𝑛−1,𝑖 = 𝑏𝑛−1,𝑖/𝜌
if 𝑙 + 𝑚 is even, or a skew component 𝐾𝑛−1,𝑖 = 𝐾 (𝑠)
𝑛−1,𝑖 = 𝑎𝑛−1,𝑖/𝜌 if 𝑙 + 𝑚 is odd, remembering 𝜌 is
the bending radius as it appeared on Eq. 2.1. The notation 𝐾𝑛−1,𝑖 is to keep up with the MAD-X

convention for naming normal multipole coefficients [21]. The symbols for 𝛽𝑥,𝑖, 𝛽𝑦,𝑖, 𝜙𝑥,𝑖 and 𝜙𝑦,𝑖

represent the unperturbed beta functions and phase advances in each plane, respectively.

2.4 Resonances in Circular Accelerators

Equation 2.30 diverges when the denominator tends to zero. Specifically, this defines a condition

reading:

( 𝑗 − 𝑘) 𝑄𝑥 + (𝑙 − 𝑚) 𝑄 𝑦 = 𝑝,

(2.33)

18

where 𝑝 can be any integer. Equation 2.33 defines resonance lines in tune space of order 𝑛 =

𝑗 + 𝑘 + 𝑙 + 𝑚. Suppose the accelerator is tuned to operate on top of these resonances. In that case,

the perturbations will add up coherently turn-to-turn and kick the resonant particles out of their

original trajectory. Operating close to or on top of a resonance line is generally harmful as particles

will be lost. This is especially true for lower-order resonances, i.e., for 𝑛 < 4. Generally, the higher

the resonance order, the weaker it is [15]. This thesis focuses on third-order resonances, i.e., 𝑛 = 3,

and how to mitigate their deleterious effect.

Figure 2.3 shows the tune diagram with resonance lines, as defined by Eq. 2.33, drawn up to

the fifth order. The integer part of both tunes is chosen to include the actual area of operation of

the Recycler Ring. The integer part of the tune carries information about the distribution of the

driving terms around the accelerator. These driving terms can be systematic terms found in every

cell of the accelerator or non-systematic terms appearing randomly around the accelerator [1].

These definitions give rise to the distinction between systematic and non-systematic resonances.

In general, a systematic resonance will be driven much stronger than a non-systematic one.

In

particular, Ch. 3 describes the operation, resonance lines and tune diagram for the Recycler Ring

in more detail. Typically, one designs the operation point of a circular accelerator to be clear of any

resonance line and as far away as possible from integer (𝑛 = 1) and half-integer (𝑛 = 2) resonances.

Nevertheless, in reality, two concepts complicate things. The first relates that resonance lines are

not infinitely thin and have some stop bandwidth. The second one concerns that particles will not

have localized tunes at high intensities but rather a distribution of tunes with some tune spread, i.e.,

a tune footprint. Section 2.7 looks at the space charge tune shift. Ultimately, choosing the operation

point on Fig. 2.3 involves localizing a resonance-free region where the intensity-dependent tune

footprint can be placed.

19

Figure 2.3 Tune diagram with resonance lines up to fifth order, enclosing the operation point of the
Recycler Ring.

It is worth stopping here and asking what is the driving force behind each of these resonance

lines. Classic accelerator references such as Refs. [1, 14, 15] will derive Eq. 2.33 by perturbing

Hill’s equation with different magnetic multipole orders. A closer look into each perturbation term

reveals that half-integer resonances are caused by quadrupole-like terms around the accelerator,

third-order resonances by sextupole-like terms, and fourth-order resonances by octupole-like terms,

following the pattern to higher orders. Nevertheless, the story is complicated when one considers

that pure multipole magnets can feed down or up to other order terms if installation misalignments

exist, e.g., a misaligned sextupole feed to skew quadrupole-like terms.

Figure 2.4 zooms into the region of interest for the Recycler Ring operation in the tune diagram,

as shown in Fig. 2.3. As mentioned before, the operation point of an accelerator in the tune diagram

is not a singular point but a footprint. While the lattice can be tuned to a specific nominal point,

20

the particles will interact with other particles through the Coulomb force. This Coulomb force acts

as a defocusing force, and it depends on the particles position in the bunch. Consequently, each

particle will feel a different tune shift depending on its position within the bunch of particles. This

effect is called the incoherent space charge tune shift, and it will be the largest for particles in the

core of the bunch, i.e., the beam core. At low particle intensities, such as the one used to produce

Fig. 2.4, the tune spread of the particles in the bunch is small enough to approximate the physics

to single-particle dynamics. For beams with low particle intensities and a minor tune spread, such

as the one depicted in Fig. 2.4, operating clear from any low-order resonance lines is not generally

a problem. Nevertheless, at high intensities, the situation changes.

Figure 2.4 Approximate operational tune footprint at low intensities calculated with PySCRDT
[22], i.e., 1e10 particles per bunch.

Figure 2.4 plots all resonance lines up to fourth order in this region of interest. The half-integer

line 𝑄𝑥 − 𝑄 𝑦 = 1, known as a difference coupling resonance, is usually driven by solenoidal- and

21

skew-quadrupole-like fields in the lattice. Sextupole-like fields drive the third-order lines 3𝑄𝑥 = 76

and 𝑄𝑥 + 2𝑄 𝑦 = 74. Skew sextupole-like terms drive the other third-order lines 3𝑄 𝑦 = 73 and

2𝑄𝑥 + 𝑄 𝑦 = 75. Moreover, finally, octupole-like fields drive the fourth order lines −𝑄𝑥 + 3𝑄 𝑦 = 48

and 3𝑄𝑥 − 𝑄 𝑦 = 52. This assumes a rectangular multipole expansion notation of the magnetic field,

such as the one presented in Eq. 2.1.

2.5 Resonance Driving Terms

The RDTs ℎ 𝑗 𝑘𝑙𝑚 are related to the strength of the resonance ( 𝑗 − 𝑘) 𝑄𝑥 + (𝑙 − 𝑚) 𝑄 𝑦. Therefore,

controlling and measuring these RDTs is of special interest to accelerator physics. The following

section explains how to get to a useful expression that can be used in order to measure the ℎ 𝑗 𝑘𝑙𝑚

terms through Fourier expansions.

The whole point of introducing the Normal Form coordinates (𝐼𝑢, 𝜓𝑢) through the transformation

𝑒−:𝐹: as defined in Eq. 2.28 is to transfer complicated non-linear dynamics to simple dynamics

that lie on a circle—𝐼𝑢 and (cid:164)𝜓𝑢 are constant. When this happens, a set of canonical coordinates

(cid:16)

(cid:174)𝜁 =

𝜁 +
𝑥 , 𝜁 −

𝑥 , 𝜁 +

𝑦 , 𝜁 −
𝑦

(cid:17) can be defined as:

𝑢 = √︁2𝐼𝑢𝑒∓𝑖(𝜓𝑢+𝜓𝑢0),
𝜁 ±

(2.34)

always keeping in mind that 𝐼𝑢 is a constant of motion and 𝜓𝑢0 is a constant initial phase set by the

initial conditions. The Poisson brackets for a pair of these quantities are:

(cid:2)𝜁 +

𝑢 , 𝜁 −
𝑢

(cid:3)

=

𝜓𝑢,𝐽𝑢

𝜕𝜁 +
𝑢
𝜕𝜓𝑢

𝜕𝜁 −
𝑢
𝜕𝐽𝑢

−

𝜕𝜁 +
𝑢
𝜕𝐽𝑢

𝜕𝜁 −
𝑢
𝜕𝜓𝑢

= −2𝑖,

(2.35)

for the same plane 𝑢 and using a reduced form of Eq. 2.18.

In this notation, the subindices

from [•, •]𝜓𝑢,𝐽𝑢 refer to the variables used to calculate the Poisson brackets. Using Eq. 2.35, the

following useful property can be derived:

(cid:104)

𝜁 +
𝑥

𝑗 𝜁 −
𝑥

𝑘 𝜁 +
𝑦

𝑙 𝜁 −
𝑦

𝑚, 𝜁 −
𝑥

(cid:105)

𝜓𝑥,𝐽𝑥

(cid:16)

𝜁 +
𝑦

𝑙 𝜁 −
𝑦

=

𝑚(cid:17) (cid:104)

𝜁 +
𝑥

𝑗 𝜁 −
𝑥

𝑘 , 𝜁 −
𝑥

(cid:105)

𝜓𝑥,𝐽𝑥

= −2𝑖 𝑗 𝜁 +
𝑥

𝑗−1𝜁 −
𝑥

𝑘 𝜁 +
𝑦

𝑙 𝜁 −
𝑦

𝑚,

(2.36)

where the last step is achieved using the Leibnitz rule for Poisson brackets, i.e., [ 𝑓 𝑔, ℎ] = [ 𝑓 , ℎ]𝑔 +

𝑓 [𝑔, ℎ].

22

On the other hand, going back to the Courant-Snyder phase space, a set of coordinates known
(cid:17) can be defined. Similarly to Eq. 2.34, the resonance basis

as a resonance basis (cid:174)ℎ =

𝑥 , ℎ−
ℎ+

𝑥 , ℎ+

𝑦 , ℎ−
𝑦

(cid:16)

reads:

𝑢 = ˆ𝑢 ± ˆ𝑝𝑢 = √︁2𝐽𝑢𝑒∓𝑖(𝜙𝑢+𝜙𝑢0),
ℎ±

(2.37)

always keeping in mind that in the Courant-Snyder phase space, the action 𝐽𝑢 is a function of the

phase 𝜙𝑢, i.e., 𝐽𝑢 = 𝐽𝑢 (𝜙𝑢) and is not constant. The initial phase 𝜙𝑢0 is again a constant set by the

initial conditions.

The basis grouped in (cid:174)ℎ and the one grouped in (cid:174)𝜁 are related by the transformation:

(cid:174)ℎ = 𝑒:𝐹

(cid:17): (cid:174)𝜁,
(cid:16) (cid:174)𝜁

(2.38)

where 𝐹 ( (cid:174)𝜁) is the generating function written in terms of the basis (cid:174)𝜁. The inverse transformation

to Eq. 2.38 reads:

(cid:174)𝜁 = 𝑒−:𝐹

(cid:17): (cid:174)ℎ.
(cid:16) (cid:174)𝜁

(2.39)

Writing out the generating function 𝐹 ( (cid:174)𝜁) in a general polynomial form, this functions reads:

(cid:17)
(cid:16) (cid:174)𝜁

𝐹

=

∑︁

𝑗 𝑘𝑙𝑚

𝑓 𝑗 𝑘𝑙𝑚 (cid:0)𝜁 +
𝑥

(cid:1) 𝑗 (cid:0)𝜁 −
𝑥

(cid:1) 𝑘 (cid:16)

(cid:17) 𝑙 (cid:16)

(cid:17) 𝑚

.

𝜁 −
𝑦

𝜁 +
𝑦

(2.40)

Inserting the definitions in Eq. 2.34 into Eq. 2.40, the proposed definition in Eq. 2.28 can be

recovered.

Expanding Eq. 2.38 by using the exponential Lie operator definition from Eq. 2.20 reads:

(cid:174)ℎ = (cid:174)𝜁 +

(cid:104)

𝐹

(cid:17)
(cid:16) (cid:174)𝜁

(cid:105)
, (cid:174)𝜁

+

(cid:104)

(cid:104)

𝐹

(cid:105) (cid:105)

𝐹, (cid:174)𝜁

+ . . . ,

1
2

(2.41)

where this expression was truncated to second order in the Poisson brackets. By taking only the

first two terms of the expansion and introducing the expression from Eq. 2.40, one can find an

approximated expression for ℎ−

𝑥 which reads:

𝑥 ≈ 𝜁 −
ℎ−

𝑥 +

(cid:104)

𝐹

(cid:17)
(cid:16) (cid:174)𝜁

(cid:105)

, 𝜁 −
𝑥

= 𝜁 −

𝑥 +

∑︁

𝑗 𝑘𝑙𝑚

(cid:104)

𝑓 𝑗 𝑘𝑙𝑚

𝜁 +
𝑥

𝑗 𝜁 −
𝑥

𝑘 𝜁 +
𝑦

𝑙 𝜁 −
𝑦

𝑚, 𝜁 −
𝑥

(cid:105)

,

(2.42)

23

At this point, the usefulness of Eq. 2.36 comes into play. Introducing the explicit result from Eq.

2.36 into Eq. 2.42 yields the following expression:

𝑥 ≈ 𝜁 −
ℎ−

𝑥 − 2𝑖 ∑︁

𝑗 𝑓 𝑗 𝑘𝑙𝑚 𝜁 +
𝑥

𝑗−1𝜁 −
𝑥

𝑘 𝜁 +
𝑦

𝑙 𝜁 −
𝑦

𝑚.

𝑗 𝑘𝑙𝑚

(2.43)

Manipulating this expression further, the definition for 𝜁𝑢 as described in Eq. 2.34 can be introduced

into Eq. 2.43. This yields:

𝑥 (𝑁) = √︁2𝐼𝑥𝑒𝑖(𝜓𝑥+𝜓𝑥0)
ℎ−

− 2𝑖 ∑︁

𝑗 𝑓 𝑗 𝑘𝑙𝑚 (2𝐼𝑥)

𝑗+𝑘−1
2

(cid:0)2𝐼𝑦(cid:1) 𝑙+𝑚

2 𝑒𝑖[(1− 𝑗+𝑘)(𝜓𝑥+𝜓𝑥0)+(𝑚−𝑙)(𝜓𝑦+𝜓𝑦0)] .

(2.44)

𝑗 𝑘𝑙𝑚

At this point, Eq. 2.44 is starting to look like a useful Fourier expansion. Ultimately, the diagnostic

data from a circular accelerator will come from a diagnostic device triggered every turn, i.e., turn-

by-turn data. For that reason, it will be useful to rewrite Eq. 2.44 in terms of the 𝑁 number of

turns of particles in the accelerator. The expression relating the phase advances to the turn number

reads:

𝜓𝑢 = 2𝜋𝑄𝑢 𝑁,

(2.45)

where 2𝜋𝑄𝑢 is the respective phase advance over one turn of the accelerator, i.e., the tune of the

circular accelerator.

Therefore, one can build the resonance basis by getting the quantity ℎ±

𝑢 = ˆ𝑢 ± ˆ𝑝𝑢 in terms of the

number of turns 𝑁 and using Eq. 2.45. Specifically, for ℎ−

𝑥 this reads:

𝑥 (𝑁) = √︁2𝐼𝑥𝑒𝑖(2𝜋𝑄 𝑥 𝑁+𝜓𝑥0)
ℎ−

− 2𝑖 ∑︁

𝑗 𝑓 𝑗 𝑘𝑙𝑚 (2𝐼𝑥)

𝑗+𝑘−1
2

(cid:0)2𝐼𝑦(cid:1) 𝑙+𝑚

2 𝑒𝑖[(1− 𝑗+𝑘)(2𝜋𝑄 𝑥 𝑁+𝜓𝑥0)+(𝑚−𝑙)(2𝜋𝑄 𝑦 𝑁+𝜓𝑦0)] ,

(2.46)

𝑗 𝑘𝑙𝑚

where 𝑄𝑥 and 𝑄 𝑦 are the horizontal and vertical uncoupled tunes. Note that one can extend this

analysis to calculate the other elements in (cid:174)ℎ.

2.6 Amplitude-Dependent Tune Shift

The RDT formalism allows one to calculate an important quantity in accelerator physics called

the amplitude-dependent tune shift. The Hamiltonian for a single particle in a linear-circular lattice

24

with perturbation elements reads:

𝐻 (𝑥, 𝑦, 𝑠) = 𝐻0(𝐽𝑥, 𝐽𝑦) + 𝐻1(𝐽𝑥, 𝜙𝑥, 𝐽𝑦, 𝜙𝑦),

(2.47)

where 𝐻0 is the unperturbed linear Hamiltonian with tunes 𝑄𝑥 and 𝑄 𝑦, and 𝐻1 is the perturbation

Hamiltonian stemming from linear and non-linear unaccounted blocks in the lattice.

From Secs. 2.1 and 2.3, the expressions for 𝐻0 and 𝐻1 have been explicitly written in Eqs. 2.14

and 2.29 in terms of 𝐽𝑥, 𝜙𝑥, 𝐽𝑦, 𝜙𝑦, therefore the sum of both expression reads:

𝐻0 + 𝐻1 = 2𝜋𝑄𝑥 𝐽𝑥 + 2𝜋𝑄 𝑦𝐽𝑦 +

∑︁

𝑗 𝑘𝑙𝑚

ℎ 𝑗 𝑘𝑙𝑚 (2𝐽𝑥)

𝑗+𝑘
2

(cid:0)2𝐽𝑦(cid:1) 𝑙+𝑚

2 𝑒𝑖[( 𝑗−𝑘)(𝜙𝑥+𝜙𝑥0)+(𝑙−𝑚)(𝜙𝑦+𝜙𝑦0)] . (2.48)

Nevertheless, it is essential to remember that 𝐻1 is the compilation of all perturbations after one

turn to the linear Hamiltonian 𝐻0 and is therefore perturbative.

Consequently, for Eq. 2.48, the independent time variable is the number of turns 𝑁. For this

case, the equations of motion taking 𝑁 as the number of turns around the circular accelerator are

just:

and

𝜕𝐽𝑢
𝜕𝑁

= −

𝜕𝐻
𝜕𝜙𝑢

= −

𝜕𝐻1
𝜕𝜙𝑢

,

𝜕𝜙𝑢
𝜕𝑁

=

𝜕𝐻
𝜕𝐽𝑢

= 2𝜋𝑄𝑢 +

𝜕𝐻1
𝜕𝐽𝑢

.

(2.49)

(2.50)

This last term 𝜕𝐻1/𝜕𝐽𝑢 in Eq. 2.50 will define the amplitude-dependent tune shift. If the whole

lattice were linear, the betatron oscillations would gain a phase change of 2𝜋𝑄𝑢 every turn, i.e.,

Δ𝜙(𝑁) = 2𝜋𝑄𝑢 𝑁. Nevertheless, given that this new term 𝐻1 perturbs the dynamics in the

accelerator, the phase change will now depend on the amplitude 𝐽𝑢 of the betatron oscillations of

the single particle. Therefore, given that this effect acts on each particle, it is incoherent. The

effective result from this new term is to detune the circular lattice from its original tune 𝑄𝑢 for each

particle.

Explicitly calculating the expression from Eq. 2.50 for 𝐽𝑥 and 𝜙𝑥 using the Hamiltonian in Eq.

2.47 yields

𝜕𝜙𝑥
𝜕𝑁

= 2𝜋𝑄𝑥 +

∑︁

𝑗 𝑘𝑙𝑚

ℎ 𝑗 𝑘𝑙𝑚 ( 𝑗 + 𝑘) (2𝐽𝑥)

𝑗+𝑘

2 −1 (cid:0)2𝐽𝑦(cid:1) 𝑙+𝑚

2 𝑒𝑖[( 𝑗−𝑘)(𝜙𝑥+𝜙𝑥0)+(𝑙−𝑚)(𝜙𝑦+𝜙𝑦0)] .

(2.51)

25

In particular, looking at the average limit where particles have undergone many turns around the

circular accelerator 𝑁 → ∞ is interesting. This is done to wash out any oscillatory behavior in

𝜕𝐻1/𝜕𝐽𝑢. Therefore, the following quantity is of interest:

lim
𝑁→∞

(cid:29)

(cid:28) 𝜕𝜙𝑥
𝜕𝑁

𝑁

= lim
𝑁→∞

1
𝑁

∫ 𝑁

0

𝑑𝑁′ 𝜕𝜙𝑥
𝜕𝑁′

,

(2.52)

which is just the definition for the average of 𝜕𝜙𝑥/𝜕𝑁 over 𝑁 turns, for many turns. Explicitly

calculating this quantity gives:

(cid:29)

(cid:28) 𝜕𝜙𝑥
𝜕𝑁

𝑁

= 2𝜋𝑄𝑥

ℎ 𝑗 𝑘𝑙𝑚 ( 𝑗 + 𝑘) lim
𝑁→∞

1
𝑁

∫ 𝑁

0

𝑑𝑁′ (2𝐽𝑥)

lim
𝑁→∞

+

∑︁

𝑗 𝑘𝑙𝑚

𝑗+𝑘

2 −1 (cid:0)2𝐽𝑦(cid:1) 𝑙+𝑚

2 𝑒𝑖[( 𝑗−𝑘)(𝜙𝑥+𝜙𝑥0)+(𝑙−𝑚)(𝜙𝑦+𝜙𝑦0)] .

(2.53)

In general, it is known that 𝐽𝑢, 𝜙𝑢 depend on the number of turns, i.e., 𝐽𝑢, 𝜙𝑢 = 𝐽𝑢 (𝑁), 𝜙𝑢 (𝑁).

Nevertheless, Eq. 2.53 can be approximated by inserting the unperturbed solution of 𝐻0, which

means that 𝐽𝑢 is constant and 𝜙𝑢 = 2𝜋𝑄𝑢 𝑁. With this in mind and assuming the constants 𝜙𝑢0 = 0

without loss of generality, Eq. 2.53 reduces to

lim
𝑁→∞

(cid:29)

(cid:28) 𝜕𝜙𝑥
𝜕𝑁

𝑁

= 2𝜋𝑄𝑥

+

∑︁

𝑗 𝑘𝑙𝑚

ℎ 𝑗 𝑘𝑙𝑚 ( 𝑗 + 𝑘) (2𝐽𝑥)

𝑗+𝑘

2 −1 (cid:0)2𝐽𝑦(cid:1) 𝑙+𝑚

2

lim
𝑁→∞

1
𝑁

∫ 𝑁

0

𝑑𝑁′ 𝑒2𝜋𝑖[( 𝑗−𝑘)𝑄 𝑥+(𝑙−𝑚)𝑄 𝑦] 𝑁 ′ .

(2.54)

A closer look into the integral in Eq. 2.54 reveals that this integral in the limit where 𝑁 → ∞ is

just a delta function reading 𝛿 (cid:0)( 𝑗 − 𝑘)𝑄𝑥 + (𝑙 − 𝑚)𝑄 𝑦(cid:1). As a reminder, the RDT approximation

breaks down if Eq. 2.33 holds. Therefore, the delta function can only be zero if 𝑗 = 𝑘 and 𝑙 = 𝑚.

Thus, this delta function effectively becomes two Kronecker deltas—𝛿 𝑗 𝑘 and 𝛿𝑙𝑚. Inserting this

into Eq. 2.54 yields:

lim
𝑁→∞

(cid:29)

(cid:28) 𝜕𝜙𝑥
𝜕𝑁

𝑁

= 2𝜋𝑄𝑥 +

∑︁

𝑗 𝑘𝑙𝑚

ℎ 𝑗 𝑘𝑙𝑚 ( 𝑗 + 𝑘) (2𝐽𝑥)

𝑗+𝑘

2 −1 (cid:0)2𝐽𝑦(cid:1) 𝑙+𝑚

2 𝛿 𝑗 𝑘 𝛿𝑙𝑚.

Reducing this expression further with the properties of Kronecker deltas reads:

lim
𝑁→∞

(cid:29)

(cid:28) 𝜕𝜙𝑥
𝜕𝑁

𝑁

= 2𝜋𝑄𝑥 + 2 ∑︁

ℎ 𝑗 𝑗𝑙𝑙 𝑗 (2𝐽𝑥) 𝑗−1 (cid:0)2𝐽𝑦(cid:1) 𝑙 .

𝑗𝑙

26

(2.55)

(2.56)

Equation 2.56 states that the terms where 𝑗 = 𝑘 and 𝑙 = 𝑚 define the constant detuning terms in

the accelerator. As a consequence of Eq. 2.31, for this case, the RDTs will be real numbers, i.e.,

ℎ 𝑗 𝑘𝑙𝑚 ∈ R. Therefore, the constant detuning terms will come from even orders of 𝑛 = 𝑗 + 𝑘 + 𝑙 + 𝑚.

It is worth remembering that this is a first approximation given that the assumption is that the

dynamics are mainly governed by 𝐻0. Higher-order approximations would involve recursively

solving Eqs. 2.49 and 2.50 as a Taylor expansion of the actions 𝐽𝑢 [23].

For example, Eq. 2.56 can be used to calculate the detuning due to horizontal quadrupole

errors 𝑛 = 2. For this, the calculation would only involve calculating ℎ1100 from Eq. 2.31, the only

surviving term. Therefore, the detuning in 𝑥 due to quadrupole errors would read

lim
𝑁→∞

(cid:29)

(cid:28) 𝜕𝜙𝑥
𝜕𝑁

𝑁

= 2𝜋𝑄𝑥 + 2ℎ1100 = 2𝜋𝑄𝑥 −

1
2

∑︁

𝑖

𝐿𝑖𝐾1,𝑖 𝛽𝑥,𝑖 = 2𝜋𝑄𝑥 −

1
2

∑︁

𝐿𝑖

𝑖

𝑏1,𝑖
𝜌

𝛽𝑥,𝑖,

(2.57)

where the sum over 𝑖 goes around all the quadrupole errors in the ring. Equation 2.57 gives a

well-known result in accelerator physics [1].

Another example is calculating the detuning terms due to octupole components around the ring,

i.e., 𝑛 = 4. In this case, the calculation yields

lim
𝑁→∞

(cid:29)

(cid:28) 𝜕𝜙𝑥
𝜕𝑁

𝑁

= 2𝜋𝑄𝑥 + 4ℎ1111𝐽𝑦 + 8ℎ2200𝐽𝑥,

(2.58)

where if ℎ1100 and ℎ2200 are explicitly calculated from Eq. 2.31, this gives:

lim
𝑁→∞

(cid:29)

(cid:28) 𝜕𝜙𝑥
𝜕𝑁

𝑁

= 2𝜋𝑄𝑥 +

𝐽𝑦
4

∑︁

𝑖

𝐿𝑖𝐾3,𝑖 𝛽𝑥,𝑖 𝛽𝑦,𝑖 +

𝐽𝑥
8

∑︁

𝑖

𝐿𝑖𝐾3,𝑖 𝛽2

𝑥,𝑖.

(2.59)

It is essential to clarify that the sum over 𝑖 goes around all the octupole perturbations in the ring.

To summarize, the amplitude-dependent tune shift Δ (2𝜋𝑄𝑢) is an important quantity that will

change the dynamics of particles inside a circular accelerator. In order to calculate this incoherent

effect to first order in perturbation theory, the quantity:

Δ (2𝜋𝑄𝑢) = 2𝜋𝑄𝑢 − lim
𝑁→∞

(cid:29)

(cid:28) 𝜕𝜙𝑢
𝜕𝑁

𝑁

= lim
𝑁→∞

(cid:29)

(cid:28) 𝜕𝐻1
𝜕𝐽𝑢

,

𝑁

(2.60)

needs to be calculated.

27

2.7 Space Charge Tune Shift

Up to this point, the beam dynamics of high-energy particle accelerators has been explained

in terms of single-particle dynamics. Up until now, there are a couple of implicit assumptions:

(a) particles do not interact with each other, and (b) the basic blocks composing the lattice have

been idealized to create the fields but without supplying any electromagnetic boundary conditions.

Nevertheless, to have a model closer to reality, one has to consider the interaction between particles

through the Coulomb force. Furthermore, particles also interact with the electromagnetic properties

of the basic elements, ultimately creating unwanted electromagnetic wake fields. This latter

phenomenon opens another branch of accelerator physics that studies collective beam instabilities

[24]. However, the scope of this thesis is only interested in the first bullet point regarding direct

particle-particle interactions through the Coulomb force—widely known as space charge physics.

It is worth specifying that this restricts the analysis to a single bunch of particles.

As mentioned in Sec. 2.4, the Coulomb force will act as a defocusing and detuning force on

each individual particle.

In order to explain this statement using a Hamiltonian formalism, the

starting point needs to be the single-particle Hamiltonian that includes the Coulomb potential from

the charge distribution in the bunch [25]. The expression for this system reads:

𝐻 (𝑥, 𝑦) = 𝐻0(𝑥, 𝑦) + 𝐻1(𝑥, 𝑦) + Ψ(𝑥, 𝑦, ˜𝑧),

(2.61)

keeping in mind that 𝑥 and 𝑦 are interchangeable for their respective action-phase variables

𝐽𝑥, 𝜙𝑥, 𝐽𝑦, 𝜙𝑦. For a bunched beam, the variable ˜𝑧 = 𝑧 − 𝑧𝑏 is introduced in order to represent

the longitudinal distance from the center of the bunch, always keeping in mind that the reference

system moves with the bunch as described by a Frenet-Serret coordinate system. The center of

the bunch with coordinates (0, 0, 𝑧𝑏) is determined from the centroid of the longitudinal charge

distribution or line density. This longitudinal line density will dictate how much detuning a given

particle experiences. Synchrotron motion—longitudinal oscillations of the particles inside the

bunch—will also cause fluctuations in the detuning force experienced by a given particle [1].

The one-turn space-charge potential Ψ, is the average potential from all the space charge kicks

around the accelerator for one turn. Nevertheless, the functional form should be similar to Eq.

28

2.31, to apply the same methods of Sec. 2.6. Therefore, this reads:

Ψ(𝐽𝑥, 𝜙𝑥, 𝐽𝑦, 𝜙𝑦, ˜𝑧) =

∑︁

𝑗 𝑘𝑙𝑚

𝐺 𝑗 𝑘𝑙𝑚 (2𝐽𝑥)

𝑗+𝑘
2

(cid:0)2𝐽𝑦(cid:1) 𝑙+𝑚

2 𝑒𝑖[( 𝑗−𝑘)(𝜙𝑥+𝜙𝑥0)+(𝑙−𝑚)(𝜙𝑦+𝜙𝑦0)] .

(2.62)

Analogous to Eq. 2.31, the terms 𝐺 𝑗 𝑘𝑙𝑚 are named the global space-charge resonance driving terms

(GSCRDTs). One can calculate these terms as the one-turn average of the instantaneous driving

terms of order 𝑙, 𝑗, 𝑘 and 𝑚, which have an explicit 𝑠 dependence. Therefore, the definition for the

GSCRDTs 𝐺 𝑗 𝑘𝑙𝑚 reads:

𝐺 𝑗 𝑘𝑙𝑚 =

1
𝐶

∫ 𝑠0+𝐶

𝑠0

˜𝑉𝑗 𝑘𝑙𝑚 (𝑠)𝑒𝑖[( 𝑗−𝑘)(𝜙𝑥 (𝑠)+𝜙𝑥0)+(𝑙−𝑚)(𝜙𝑦 (𝑠)+𝜙𝑦0)] 𝑑𝑠,

(2.63)

where the terms ˜𝑉𝑗 𝑘𝑙𝑚 (𝑠) are the instantaneous driving terms of the expansion for the self-field

potential 𝜓 of the bunch at a point 𝑠 of the accelerator, i.e.,

𝜓(𝐽𝑢, 𝜙𝑢, ˜𝑧, 𝑠) = 𝑓 ( ˜𝑧) +

∑︁

𝑗 𝑘𝑙𝑚

˜𝑉𝑗 𝑘𝑙𝑚 (𝑠) (2𝐽𝑥)

𝑗+𝑘
2

(cid:0)2𝐽𝑦(cid:1) 𝑙+𝑚

2 𝑒𝑖[( 𝑗−𝑘)(𝜙𝑥+𝜙𝑥0)+(𝑙−𝑚)(𝜙𝑦+𝜙𝑦0)] .

(2.64)

The term 𝑓 ( ˜𝑧) will hold any longitudinal dependence of the potential, but it is generally

irrelevant to the space charge detuning calculation. This approach is similar to the one taken by

Ref. [26]. This assumes the self-field potential has been calculated and a Floquet transformation

has been performed using Eq. 2.8. The self-potential 𝜓(𝑥, 𝑦, ˜𝑧, 𝑠) is determined self-consistently

from the 3D Poisson equation with a particle number density of the bunch 𝑛𝑏 (𝑥, 𝑦, ˜𝑧, 𝑠) reading:

(cid:18) 𝜕
𝜕𝑥

+

𝜕
𝜕𝑦

(cid:19)

𝜕
𝜕 ˜𝑧

+

𝜓(𝑥, 𝑦, ˜𝑧, 𝑠) = −

2𝜋𝐾𝑏
𝑁𝑏

𝑛𝑏 (𝑥, 𝑦, ˜𝑧, 𝑠),

where 𝐾𝑏 is a dimensionless parameter known as the self-field perveance defined as

𝐾𝑏 =

2𝑁𝑏𝑒2
𝑏
𝛾3
𝐿𝑚𝑏 𝛽2
𝐿𝑐2

,

(2.65)

(2.66)

with 𝑁𝑏 being the total number of particles in the bunch defined as 𝑁𝑏 = ∫ 𝑑𝑥 𝑑𝑦 𝑑 ˜𝑧 𝑛𝑏 (𝑥, 𝑦, ˜𝑧, 𝑠),

𝑒𝑏 is the charge of one beam particle, 𝑚𝑏 is the rest mass of one beam particle, 𝛾𝐿 and 𝛽𝐿 being the

relativistic longitudinal factors of a beam with total energy 𝐸𝑇 = 𝛾𝐿𝑚𝑏𝑐2, and 𝑐 being the speed of

light.

29

One can find the solution to Eq. 2.65 in any book of mathematical methods for physics, see Refs.

[27, 28]. The solution to this equation involves using Green’s function and finding the convolution

with the particle number density 𝑛𝑏 (𝑥, 𝑦, ˜𝑧, 𝑠) = 𝑛𝑏 ((cid:174)𝑟, 𝑠). This reads:

𝜓(𝑥, 𝑦, ˜𝑧, 𝑠) = −

2𝜋𝐾𝑏
𝑁𝑏

∫

R3

𝑑(cid:174)𝑟′ 𝑛𝑏 ((cid:174)𝑟, 𝑠)
|(cid:174)𝑟 − (cid:174)𝑟′|

.

(2.67)

In section 2.6, it was interesting to look at how particles at different amplitudes undergo different

detuning due to 𝐻1 given their 𝐽𝑢 coordinates. The same will be valid for the effect of Ψ on the

single-particle dynamics. The quantity of interest will be the detuning due to Coulomb forces—

the space-charge tune shift. Nevertheless, it is essential to note that the space charge will be a

continuous force around the accelerator. Therefore, the terms 𝐺 𝑗 𝑘𝑙𝑚 in order to enter the one-turn

Hamiltonian in Eq. 2.61 will have to be an average of the space-charge force around the ring for

a single particle—as defined in Eq. 2.63. Ultimately, the space-charge tune shift will also be an

incoherent quantity.

For now, let 𝐻1 = 0, but let space-charge dictate the functional form of Ψ. For this case, and in

analogy to Eq. 2.50, the space charge tune shift will be dictated by:

𝜕𝜙𝑢
𝜕𝑁

=

𝜕𝐻
𝜕𝐽𝑢

= 2𝜋𝑄𝑢 +

𝜕Ψ
𝜕𝐽𝑢

.

(2.68)

The non-trivial part of this calculation is figuring 𝜕Ψ/𝜕𝐽𝑢 for a specific bunch distribution after

calculating the integral in Eq. 2.67.

The first case to analyze is when the bunch has a uniform charge density enclosed in a 3D

ellipsoid. The charge bunch distribution reads:

𝑛𝑏 (𝑥, 𝑦, ˜𝑧, 𝑠) =





ˆ𝑛𝑏 (𝑠),

if 0 ≤ 𝑥2
𝑎2 (𝑠)

0,

+ ˜𝑧2
𝑐2 (𝑠)

≤ 1

,

+

𝑦2
𝑏2 (𝑠)
if else.

(2.69)

where 𝑎(𝑠), 𝑏(𝑠) and 𝑐(𝑠) are the length of the semi-axes of the ellipsoid around the accelerator

and can be calculated from the beam-envelope equations [25]. Additionally, the constant particle

number density over an ellipsoid is defined as ˆ𝑛𝑏 = 𝑁𝑏/𝑉𝑒, where 𝑉𝑒 is the volume of an ellipsoid

reading 𝑉𝑒 = (4/3)𝜋𝑎𝑏𝑐. Inputting this distribution into Eq. 2.67 gives the expression for the

30

potential. This involves solving a 3D integral inside the ellipsoid. This reads:

𝜓(𝑥, 𝑦, ˜𝑧, 𝑠) = −

2𝜋𝐾𝑏 ˆ𝑛𝑏
𝑁𝑏

∫

S

𝑑 (cid:174)𝑟′

1
|(cid:174)𝑟 − (cid:174)𝑟′|

,

(2.70)

where S is the region enclosed by the ellipsoid. The integral has been studied with ellipsoidal

coordinates, and the solution can be found in Refs. [25, 29]. In these references, it is shown that

the solution can be expressed as a quadratic function of 𝑥, 𝑦 and ˜𝑧 such that:

𝜓(𝑥, 𝑦, ˜𝑧, 𝑠) = −

𝜋𝐾𝑏 ˆ𝑛𝑏𝑎𝑏𝑐
𝑁𝑏

(cid:2)−𝐴𝑥2 − 𝐵𝑦2 − ˜𝐶 ˜𝑧2 + 𝐷(cid:3) ,

(2.71)

such that 𝐴, 𝐵, ˜𝐶 and 𝐷 are calculated as:

𝐴(𝑠) =

∫ ∞

0

𝑑𝜍
(cid:0)𝑎2 + 𝜍(cid:1) √︃(cid:0)𝑎2 + 𝜍(cid:1) (cid:0)𝑏2 + 𝜍(cid:1) (cid:0)𝑐2 + 𝜍(cid:1)

𝐵(𝑠) =

∫ ∞

0

˜𝐶 (𝑠) =

∫ ∞

0

𝐷 (𝑠) =

𝑑𝜍
(cid:0)𝑏2 + 𝜍(cid:1) √︃(cid:0)𝑎2 + 𝜍(cid:1) (cid:0)𝑏2 + 𝜍(cid:1) (cid:0)𝑐2 + 𝜍(cid:1)
𝑑𝜍
(cid:0)𝑐2 + 𝜍(cid:1) √︃(cid:0)𝑎2 + 𝜍(cid:1) (cid:0)𝑏2 + 𝜍(cid:1) (cid:0)𝑐2 + 𝜍(cid:1)
∫ ∞
𝑑𝜍
√︃(cid:0)𝑎2 + 𝜍(cid:1) (cid:0)𝑏2 + 𝜍(cid:1) (cid:0)𝑐2 + 𝜍(cid:1)

0

.

,

,

,

(2.72)

(2.73)

(2.74)

(2.75)

Furthermore, Eq. 2.71 can be rewritten using Eq. 2.8 to highlight the dependence on 𝐽𝑢, 𝜙𝑢

explicitly. This reads:

𝜓(𝐽𝑢, 𝜙𝑢, ˜𝑧, 𝑠) = −

(cid:20)

𝜋𝐾𝑏 ˆ𝑛𝑏𝑎𝑏𝑐
𝑁𝑏

−𝐴 (2𝛽𝑥 𝐽𝑥)

−𝐵 (cid:0)2𝛽𝑦𝐽𝑦(cid:1)

(cid:19)

(cid:18) 𝑒2𝑖(𝜙𝑥+𝜙𝑥0) + 𝑒−2𝑖(𝜙𝑥+𝜙𝑥0) + 2
4
(cid:18) 𝑒2𝑖(𝜙𝑦+𝜙𝑦0) + 𝑒−2𝑖(𝜙𝑦+𝜙𝑦0) + 2
4

(cid:19)

(cid:21)

− ˜𝐶 ˜𝑧2 + 𝐷

.

(2.76)

Equation 2.76 can be cross-examined using Eq. 2.64 in order to retrieve the driving terms ˜𝑉𝑗 𝑘𝑙𝑚 and

𝑓 ( ˜𝑧). Therefore, for a uniform distribution, the following expansion for the self-potential holds:

𝜓(𝐽𝑢, 𝜙𝑢, ˜𝑧, 𝑠) = 𝑓0( ˜𝑧) + ˜𝑉2000𝐽𝑥𝑒2𝑖(𝜙𝑥+𝜙𝑥0) + ˜𝑉0200𝐽𝑥𝑒−2𝑖(𝜙𝑥+𝜙𝑥0) + ˜𝑉1100𝐽𝑥

+ ˜𝑉0020𝐽𝑦𝑒2𝑖(𝜙𝑦+𝜙𝑦0) + ˜𝑉0002𝐽𝑦𝑒−2𝑖(𝜙𝑦+𝜙𝑦0) + ˜𝑉0011𝐽𝑦.

(2.77)

31

The instantaneous driving terms ˜𝑉𝑗 𝑘𝑙𝑚 from Eq. 2.77 are summarized in Table 2.2. Table 2.2 also
shows how to calculate the SCRDTs 𝐺 𝑗 𝑘𝑙𝑚 from averaging out the ˜𝑉𝑗 𝑘𝑙𝑚 (𝑠) all around the circular

accelerator as an example of using Eq. 2.63. One can calculate the integrals for the 𝐺 𝑗 𝑘𝑙𝑚 with

previous knowledge of the unperturbed beta functions and phase advances around the ring. The

envelope equations, as described in Ref. [25], should also be solved in order to get 𝑎(𝑠), 𝑏(𝑠) and

𝑐(𝑠) around the lattice, ultimately, aiding in the evaluation of 𝐴(𝑠), 𝐵(𝑠) and 𝐶 (𝑠).

Table 2.2 Instantaneous driving terms and GSCRDTs for a uniform beam.
˜𝑉𝑗 𝑘𝑙𝑚

˜𝑉𝑗 𝑘𝑙𝑚 (𝑠)

𝐺 𝑗 𝑘𝑙𝑚

𝐺 𝑗 𝑘𝑙𝑚

˜𝑉2000

3
8 𝐾𝑏 𝐴(𝑠) 𝛽𝑥 (𝑠) 𝐺2000

3𝐾𝑏
8𝐶

∫ 𝑠0+𝐶

𝑠0

𝐴(𝑠) 𝛽𝑥 (𝑠)𝑒2𝑖(𝜙𝑥 (𝑠)+𝜙𝑥0) 𝑑𝑠

˜𝑉0200

3
8 𝐾𝑏 𝐴(𝑠) 𝛽𝑥 (𝑠) 𝐺0200

3𝐾𝑏
8𝐶

∫ 𝑠0+𝐶

𝑠0

𝐴(𝑠) 𝛽𝑥 (𝑠)𝑒−2𝑖(𝜙𝑥 (𝑠)+𝜙𝑥0) 𝑑𝑠

˜𝑉1100

3
4 𝐾𝑏 𝐴(𝑠) 𝛽𝑥 (𝑠) 𝐺1100

3𝐾𝑏
4𝐶

∫ 𝑠0+𝐶

𝑠0

𝐴(𝑠) 𝛽𝑥 (𝑠) 𝑑𝑠

˜𝑉0020

3
8 𝐾𝑏 𝐵(𝑠) 𝛽𝑦 (𝑠) 𝐺0020

3𝐾𝑏
8𝐶

∫ 𝑠0+𝐶

𝑠0

𝐵(𝑠) 𝛽𝑦 (𝑠)𝑒2𝑖(𝜙𝑦 (𝑠)+𝜙𝑦0) 𝑑𝑠

˜𝑉0002

3
8 𝐾𝑏 𝐵(𝑠) 𝛽𝑦 (𝑠) 𝐺0002

3𝐾𝑏
8𝐶

∫ 𝑠0+𝐶

𝑠0

𝐵(𝑠) 𝛽𝑦 (𝑠)𝑒−2𝑖(𝜙𝑦 (𝑠)+𝜙𝑦0) 𝑑𝑠

˜𝑉0011

3
4 𝐾𝑏 𝐵(𝑠) 𝛽𝑦 (𝑠) 𝐺0011

3𝐾𝑏
4𝐶

∫ 𝑠0+𝐶

𝑠0

𝐵(𝑠) 𝛽𝑦 (𝑠) 𝑑𝑠

Just like in Sec. 2.6, looking at the terms where 𝑗 = 𝑘 and 𝑙 = 𝑚 is interesting. This calculation

will give nonlinear detuning terms due to a uniform charge distribution. Therefore, similar to Eq.

32

2.56, the nonlinear detuning terms due to space charge for the horizontal plane are:

lim
𝑁→∞

(cid:29)

(cid:28) 𝜕𝜙𝑥
𝜕𝑁

𝑁

= 2𝜋𝑄𝑥 + 2 ∑︁

𝐺 𝑗 𝑗𝑙𝑙 𝑗 (2𝐽𝑥) 𝑗−1 (cid:0)2𝐽𝑦(cid:1) 𝑙 .

(2.78)

𝑗𝑙

Therefore, it is interesting to look at the terms 𝐺1100 and 𝐺0011, which will cause constant detuning

of the unperturbed tunes.

One can do the previous analysis for a uniform beam, i.e., the particle number density is constant

inside the bunch. Nevertheless, this is an idealized case for a circular accelerator. A more realistic

approach would be to consider a Gaussian distribution. This multivariate distribution, assuming

no correlation, reads:

𝑛𝑏 (𝑥, 𝑦, ˜𝑧, 𝑠) =

1
(2𝜋)3/2 𝜎𝑥𝜎𝑦𝜎˜𝑧

𝑒

− 𝑥2
2𝜎2
𝑥

−

𝑦2
2𝜎2
𝑦

− ˜𝑧2
2𝜎2
˜𝑧

(2.79)

This approach is the approach taken in Refs. [22, 26]. Specifically, the Python library developed at

CERN called PySCRDT [22] allows for the calculation of GSCRDTs and the tune footprint for a

Gaussian beam distribution. As mentioned before in the discussion of Fig. 2.4, the tune footprint

depends on the intensity of the beam interchangeable with the beam space charge perveance

parameter 𝐾𝑏. Figure 2.5 shows the space charge tune footprint for the Recycler Ring for typical

beam parameters used in studies. One can calculate the footprint using PySCRDT. Nevertheless,

the approach behind PySCRDT is similar to the example provided with a uniform beam distribution.

One special feature of the Gaussian beam is that the space charge potential is its largest at the core of

the beam. Hence, the detuning is largest for particles close to the middle of the bunch. In contrast,

the particles away from the core of the beam will not be subject to the full space charge potential

and, therefore, will oscillate close to the unperturbed nominal tune of 𝑄𝑢. Figure 2.5 illustrates this

concept by drawing a color map at different actions, close to the beam core and 3𝜎𝑢 away from the

center of the bunch—location where particles oscillate close to nominal tunes.

33

Figure 2.5 Tune footprint for a Gaussian beam in the Recycler Ring at an intensity of 5e10 particles
per bunch created with PySCRDT.

Up until now, the assumption 𝐻1 = 0 has held for this section. Nevertheless, there will be

an interplay between the nonlinear Hamiltonian and the space charge potential in high-intensity

particle accelerators. Ultimately, the terms ℎ 𝑗 𝑘𝑙𝑚 and 𝐺 𝑗 𝑘𝑙𝑚 will dictate the perturbation to the

linear Hamiltonian by their effective detuning. Therefore, the following Hamiltonian will describe

the perturbation due to nonlinear elements in the circular lattice and due to space charge forces:

𝐻 (𝐽𝑥, 𝜙𝑥, 𝐽𝑦, 𝜙𝑦, ˜𝑧) = 2𝜋𝑄𝑥 𝐽𝑥 + 2𝜋𝑄 𝑦𝐽𝑦+

(cid:0)ℎ 𝑗 𝑘𝑙𝑚 + 𝐺 𝑗 𝑘𝑙𝑚(cid:1) (2𝐽𝑥)

𝑗+𝑘
2

(cid:0)2𝐽𝑦(cid:1) 𝑙+𝑚

2 𝑒𝑖[( 𝑗−𝑘)(𝜙𝑥+𝜙𝑥0)+(𝑙−𝑚)(𝜙𝑦+𝜙𝑦0)] .

(2.80)

∑︁

𝑗 𝑘𝑙𝑚

34

CHAPTER 3

THE FNAL RECYCLER RING

3.1

Introduction

The Fermilab Recycler Ring (RR) is one of the circular accelerators located in the Fermilab

Accelerator Complex. It was originally designed to store and accumulate antiprotons that remained

from a Tevatron event [30]. The recycling of antiprotons was deemed ineffective and was never

operationally implemented [31]. Since 2011, the RR has been repurposed as a pre-injector to the

Main Injector (MI) by storing and accumulating protons [32].

It is worth pointing out that the

MI and the RR share the same tunnel, with a circumference of 3.319 km (2.062 mi). The work

done for this thesis focuses on the Recycler Ring. The following chapter starts by giving a general

description of the operation and physics of the Recycler Ring. The following sections introduce

and motivate the compensation of third-order resonances for high-intensity operation.

The MI/RR complex is fed protons by the Proton Source, which by itself consists of the Pre-

Accelerator, the Linear Accelerator (Linac), and the Booster. The Pre-Accelerator systems provide

𝐻− ions to the Linac, where they are accelerated to an energy of 400 MeV. After this, the beam

is injected into the Booster Ring. The Booster is a rapid-cycling synchrotron operating at a 15

Hz repetition rate. During this injection process, the 𝐻− beam passes through a carbon stripping

foil and incorporates into the circulating proton beam. The Booster ramps the energy up from

400 MeV to 8 GeV. This 8 GeV proton beam can either go to the Booster Neutrino Experiments

or get injected into the Recycler Ring. Once in RR, the beam has two possible destinations: 1)

high-energy neutrino experiments through MI or 2) Muon Campus. For the latter, the proton beam

gets rebunched from 53 MHz to 2.5 MHz and transported to the Muon Campus. The proton beam

gets slip-stacked for high-energy neutrino experiments, doubling the intensity injected into the

Main Injector. Once in MI, the beam is accelerated to 120 GeV and sent to the NuMI (Neutrinos at

the Main Injector) beam facility [13, 31, 32]. Figure 3.1 describes the current accelerator complex,

including the experimental beamlines that feed neutrino, muon, and fixed target experiments.

35

Figure 3.1 Current operational layout of the Fermilab Accelerator Complex as of 2024. The original
plot was provided by R. Ainsworth, but modified for this document.

The Proton Improvement Plan II (PIP-II) is the first step in establishing the Fermilab Accelerator

Complex as a multi-MW proton facility [33]. The near-future objective is to deliver a 1.2 MW

proton beam to the Deep Underground Neutrino Experiment (DUNE) through the Long-Baseline

Neutrino Facility (LBNF) [34], still in construction.

In order to meet this goal, Fermilab is

planning several upgrades in the accelerator complex, including a new 800 MeV superconducting

linear accelerator. The plan for the layout of the Fermilab Accelerator Complex is shown in Fig.

3.2. With minimal upgrades to the Main Injector and Recycler Ring but a substantial overhaul of

the Booster Ring, this will allow for a 50% increase in particles per pulse intensity. Table 3.1 also

specifies some upgrades that will happen for the PIP-II era. Some examples include an increase

of the particle per bunch intensity, a shortening of the Main Injector acceleration ramp, and an

increase in the Booster ramping rate. As the Recycler Ring starts to deal with higher intensities

from the PIP-II upgrade, it is important to mitigate the effects of space charge, as discussed in Secs.

2.4 and 2.7. Particles along the bunch will experience space charge forces leading to detuning

in their betatron frequencies. Given the incoherent nature of this process, this leads to the beam

having a larger tune spread in the tune diagram and having particles operate on top of resonances.

36

Figure 3.2 Future layout of the Fermilab Accelerator Complex for the Proton Improvement Plan II
(PIP-II). The original plot was provided by R. Ainsworth but modified for this document.

3.2 General Description

The RR is a permanent magnet storage ring operating at a fixed momentum of 8.835 GeV/c,

equivalent to an energy of 8 GeV. The basic cell structures of this machine are FODO (Focusing

Quadrupole - Drift - Defocusing Quadrupole - Drift) cells. During its conception, the need for a

quick and non-expensive design spurred the idea of combining quadrupole and dipole magnets into

one combined function magnet. Figure 3.3 shows these combined function magnets as the green-

covered magnets on the top ring—the Recycler Ring. In order to further reduce costs during its

construction, these magnets were chosen to be permanent magnets made out of a strontium ferrite

[30]. One advantage of having permanent magnets is that there is no need for power supplies,

cooling systems, or power distribution cables. Consequently, these types of magnets are very stable

against time and temperature. Nevertheless, the magnetic field of such magnets does degrade over

time. Reference [32] shows how the fields in RR-type magnets can degrade around 1% after 20

years. Ultimately, this slightly changes the nominal energy in the machine.

37

Figure 3.3 Picture of the Main Injector (blue and red magnets in the bottom) and the Recycler Ring
(green magnets up top) tunnel.

Figure 3.3 reminds the reader that the Main Injector and the Recycler Ring share the same

tunnel. This tunnel is divided into six sections with labels: 100, 200, 300, 400, 500 and 600.

Injection from the Booster Ring into the Recycler takes place just after the beginning of section

100. Specifically, this happens at a Lambertson magnet labeled as LAM102. This case is just

an example of how every RR or MI element will be labeled according to its position in one of

these sections. Figure 3.4 shows a schematic of the Recycler Ring section with some labels for

critical subsystems inside them. In particular, it shows the location close to 232, where the beam

is transferred from the Recycler Ring to the Main Injector. Figure 3.4 also presents the location in

section 500 where the beam is transferred from RR to the Muon Campus, as well as the location

after 401 where the beam is dumped towards the abort line [32, 35].

As highlighted in Fig. 3.4, a critical subsystem is the tune trombone. The RR has two of

them, one in sections 601-609 and another in 301-309. A tune or phase trombone is a linear insert

38

composed of quadrupoles introducing local phase advances. These quadrupoles are powered to add

a controlled phase advance while leaving the Twiss parameters unchanged at the end of the insert

and matched to the start of the trombone [30, 36]. Ultimately, these subsystems introduce local tune

changes that allow the control of both tunes in a ±0.5 range, i.e., Δ𝑄𝑢 = ±0.5. The Recycler Ring

has accelerator applications through ACNET that allow the control of these tune trombones and set

the tunes to some desired values. Furthermore, this application allows tune ramps and continuous

linear changes in the tunes. This feature is crucial for building dynamic loss maps—an important

tool to probe resonance lines as will be described in Sec. 4.5.1.

Figure 3.4 Schematic layout of the Recycler Ring and its corresponding sections. Original plot
provided by R. Ainsworth, first published on Ref. [32].

The Recycler Ring has 104 FODO cells distributed into three main groups around its 3319.4-

meter circumference [30]. The first group of cells comprises two permanent magnet quadrupoles

arranged so that each cell has a phase advance of 90◦. There are 18 cells of this first group

throughout the Recycler. The second group comprises a cell with two combined function magnets

39

to simultaneously bend and focus the beam. This type of cell also has a phase advance of 90◦, with

54 cells of this type in the RR. Ultimately, the configuration of these cells will dictate the Twiss

parameters around the ring. Figure 3.5 shows the beta functions 𝛽𝑢 (𝑠) around the Recycler Ring

tuned to some particular set of tunes. Some horizontal BPMs are marked in the upper horizontal

axis to note the corresponding section in the tunnel.

Before describing the third group, it is worth clarifying to the reader that the transverse motion

inside any circular accelerator is dictated by:

𝑢(𝑠) = 𝑢 𝛽 (𝑠) + 𝐷𝑢 (𝑠)𝛿,

(3.1)

where 𝑢 is either the 𝑥 or 𝑦 plane, 𝑢 𝛽 (𝑠) is the betatron motion as described by Eq. 2.8, 𝐷𝑢 (𝑠)

is known as the dispersion function and 𝛿 = Δ𝑝/𝑝0 is the fractional momentum deviation of an

individual particle with respect to the reference longitudinal momentum 𝑝0. In Ch. 2, specifically

in Sec. 2.2, an assumption was made that the transverse coordinates were not going to be influenced

by any longitudinal property, focusing on on-momentum particles, i.e., 𝛿 = 0. Nevertheless, for

this section, it is relevant to introduce the dispersion function 𝐷𝑢 (𝑠) to describe the Recycler Ring

in all its glory and detail. The dispersion function depends on the distribution of quadrupole and

dipole components throughout the ring [1]. Even so, the scope of this work regarding resonance

compensation will not consider any longitudinal-transverse coupling.

With all this said, the third group comprises special dispersion suppressor cells, 𝐷𝑥 (𝑠) = 0 for

these regions. Each one has two combined function magnets and globally has a betatron phase

advance of again 90◦. This cell type will allow dispersion-free regions necessary for the Recycler’s

injection, extraction, and other particular subsystems. There are 32 cells of this type. Figure 3.6

shows the dispersion functions for the Recycler Ring. Thanks to this third type of FODO cells, the

horizontal dispersion function 𝐷𝑥 (𝑠) has regions of close-to-zero values. The vertical dispersion

is negligible, i.e., 𝐷 𝑦 (𝑠) ≈ 0.

40

Figure 3.5 Beta functions for the Recycler Ring lattice tuned to 𝑄𝑥 = 25.44 and 𝑄 𝑦 = 24.39. Lattice
functions are calculated from the lattice file using SYNERGIA.

Figure 3.6 Dispersion functions for the Recycler Ring lattice tuned to 𝑄𝑥 = 25.44 and 𝑄 𝑦 = 24.39.
Lattice functions are calculated from the lattice file using SYNERGIA.

The ultimate role of the Recycler Ring is to smooth the way for high-intensity proton beam

injection into the Main Injector. In particular, for the high-intensity beam delivered to the neutrino

experiments, the Recycler performs a beam manipulation known as slip-stacking. Although this

longitudinal manipulation is out of the main scope of this thesis, it is relevant to briefly explain

this procedure since this is the Recycler Ring’s main feature. Slip-stacking is the process through

which a pair of bunches are manipulated in longitudinal phase space to merge them into one higher-

intensity bunch [35]. In particular, the RR utilizes two RF stations to decelerate and accelerate

trains of bunches. The Recycler can fit 7 Booster Rings. Hence, 6+6 beam batches are injected

41

and slip-stacked, while 1 batch space is left as a gap for the kickers to fire. Each batch is composed

of 81 bunches following the harmonic number of the Recycler. A more detailed explanation of

the Recycler Ring’s slip-stacking procedure can be found in Refs. [32, 35, 37]. Ultimately, this

slip-stacking process is relevant because it doubles the particles per bunch intensity (ppb), which

has significant consequences for the transverse dynamics, as explained in Sec. 2.7. Furthermore,

Table 3.1 summarizes some important parameters for the RR slip-stacking and the general operation

of the neutrino-bound beam.

Table 3.1 Typical Recycler Ring properties for beam sent to NuMI, with some PIP-II nominal
parameters.

Parameter

Value

Circumference
Momentum
Revolution Period
Revolution Frequency
RF Frequency
RF Voltage
Harmonic Number
Synchrotron Tune
Slip Factor
Superperiodicity
Horizontal Tune
Vertical Tune
Horizontal Chromaticity
Vertical Chromaticity
95% Normalized Emittance
95% Longitudinal Emittance
Intensity

MI Ramp Time

Booster Frequency

3319.4
8.835
11.1
90.1
52.8
80
588
0.0028
-8.6 × 10−3
2
25.43
24.445
-6
-7
15
0.08
5 × 1010
8 × 1010 (PIP-II)
1.2
1.133
1.067
15
20 (PIP-II)

Unit

m
GeV/c
𝜇s
kHz
MHz
kV

𝜋 mm mrad
eV s
ppb
ppb
s
s
s
Hz
Hz

3.3 Tune Diagram and Resonances

Table 3.1 specifies nominal horizontal and vertical tunes for the operation of the Recycler Ring.

While this pair of tunes can be switched around during operation to reduce losses, they will always

42

stay clear of the surrounding transverse betatron resonances, as introduced in Sec. 2.4. Figure 3.7

plots the resonance lines relevant to the Recycler Ring operation, with their corresponding order

in the multipole expansion. Specifically, this study looks at normal sextupole lines 3𝑄𝑥 = 76 and

𝑄𝑥 + 2𝑄 𝑦 = 74, plus skew sextupole lines 3𝑄 𝑦 = 73 and 2𝑄𝑥 + 𝑄 𝑦 = 75.

Figure 3.7 Portion of the tune diagram enclosing the operational tunes of the Recycler Ring.

Starting from the lowest order resonance line, Fig. 3.7 shows the second-order linear coupling

line 𝑄𝑥 − 𝑄 𝑦 = 1. This resonance line comes from the skew-quadrupole component accumulated

all around the ring. While operating on this resonance line will not lead to chaotic motion, the

43

dynamics are characterized by action exchange between the betatron modes [1]. This exchange

can cause emittance growth, eventually leading to beam loss. Due to this fact, the Recycler Ring

has two families of skew quadrupoles, which correct this resonance line. Consequently, operating

close to this resonance line is feasible with this correction.

The third-order resonance lines are the main focus of this thesis. They are driven by normal-

and skew-sextupole components distributed around the ring.

In particular, it is known that the

combined function magnets are the main drivers of these resonance lines. Even though these

magnets were designed to be quadrupole and dipole only, they can still have sextupole components

from construction errors. Nevertheless, since they are permanent magnets, little can be done about

the magnets themselves. As a result, the following work explores a compensation technique to

bring down the lattice RDTs that drive these resonance lines. A summary of the RDTs, their

corresponding resonance line, and the frequency position of their spectral lines are summarized in

Table 3.2. The spectral lines are a concept critical to the measurement of RDTs, and it is introduced

in the next chapter. The compensation of these RDTs involves the usage of additional compensation

sextupoles. Section 3.6 describes these elements. Chapter 4 goes into a deep dive explaining this

approach from a theoretical and experimental point of view.

Table 3.2 Corresponding RDTs and location of spectral lines for each resonance line.

Resonance Line
3𝑄𝑥 = 76
𝑄𝑥 + 2𝑄 𝑦 = 74
3𝑄 𝑦 = 73
2𝑄𝑥 + 𝑄 𝑦 = 75

Source

Normal Sextupole
Normal Sextupole
Skew Sextupole
Skew Sextupole

RDT Hor. Spect. Vert. Spect.
ℎ3000
ℎ1020
ℎ0030
ℎ2010

(-2,0)
(0,-2)
-
(-1,-1)

-
(-1,-1)
(0,-2)
(-2,0)

Finally, there are also fourth-order resonance lines close to the nominal operation of the Recycler.

Octupole components around the ring drive these lines. In practice, −𝑄𝑥 +3𝑄 𝑦 = 48 and 3𝑄𝑥 −𝑄 𝑦 =

52 are not particularly concerning to the operation of the RR. Their corresponding lattice RDTs

are not large enough to consider them dangerous during the RR cycle, i.e., the beam is stored for

approximately 1 second. Nevertheless, these lines are interesting because fourth-order resonance

44

lines are susceptible to being amplified by the space charge potential through the SCRDTs, as

described in Sec. 2.7.

3.4 High Intensity and Tune Footprint

In order to achieve the PIP-II beam power objective, the Recycler will be required to store and

accumulate 50% more beam than current operations [33]. As explained in Sec. 2.7, larger beam

intensities lead to larger space charge tune spread. Figure 3.8 shows how, for the typical operation

to the neutrino experiments under future PIP-II specifications, there is a space charge tune spread of

around 0.1 in both planes. At nominal tunes, it is clear that for high-intensity operation, the particles

in the core of the beam will start to operate on top of third-order resonances. One feature further

complicates this picture: particles will undergo synchrotron oscillations within the diamond-shaped

region. Particles will go from the core of the beam to its outskirts, and its incoherent space charge

tune shift will fluctuate accordingly. Therefore, particles will cross in and out of these third-order

resonances at a synchrotron period. Ultimately, Fig. 3.8 just tries to take a snapshot of the region

of these excursions.

45

Figure 3.8 Approximate operational tune footprint at high intensities, i.e., 1e11 particles per bunch.

Figure 3.8 and the latter discussion motivate the need for mitigating the deleterious effects of

3𝑄𝑥 = 76 and 𝑄𝑥 + 2𝑄 𝑦 = 74, and weaker effects from 3𝑄 𝑦 = 73 and 2𝑄𝑥 + 𝑄 𝑦 = 75. If one

reduces the RDT of each resonance line, then their strength becomes less, losses in those regions

are reduced, and ultimately there will be a larger area in the tune diagram for operation. This

mechanism is desirable for PIP-II. Nevertheless, Ch. 4 describes compensating one resonance line

might make another one worse and vice versa. Here lies the difficulty of this approach. Additionally,

these third-order resonance lines must be adequately characterized at low intensities before going

to PIP-II intensities. Experimentally, everything starts to blend together in the tune diagram at high

intensities due to the space charge tune shift. This amalgam of phenomena motivates the fact to

have a distinction between Ch. 4 and 6.

46

3.5 Diagnostic Devices

The Recycler Ring has two primary diagnostic devices that fall under the scope of this work:

the Beam Position Monitors (BPMs) and the Ion Profile Monitors (IPMs). Although both systems

quantify the properties of the beam, each one gives different information about the beam distribution.

In particular, BPMs probe the first moment of the transverse distribution, i.e., the bunch’s mean

position or beam centroid. Depending on the BPM system, this measurement is done in one

plane or simultaneously in both planes. On the other hand, IPMs go one order higher and give

information about the second-order moment of the transverse beam distribution—information about

the variance and spread of the bunch. Specifically for the Recycler case, there is one IPM for the

horizontal direction and another for the vertical case.

There are 208 BPMs in the Recycler Ring in total. Specifically, there are 104 horizontal BPMs

and 104 vertical BPMs. Each class is oriented in order to measure the corresponding plane. The

BPMs are labeled according to their position around the ring as described in Sec. 3.2 and in Fig.

3.4. One BPM consists of two parallel pick-up electrodes that produce an electric signal once the

beam passes. The beam position is determined from the relative amplitude between the signals of

the opposing channels [38]—also known as the ( 𝐴 − 𝐵)/( 𝐴 + 𝐵) signal. This signal is digitized and

calibrated to include the scaling factors and offsets, ultimately, in order to represent the transverse

beam position. The digitized signal from all 208 BPMs is interfaced with ACNET for accelerator

applications. The resulting data is digitized every turn, which is called turn-by-turn (TbT) data.

Figure 3.9 shows an example of this TbT BPM data for a beam pinged in the horizontal direction

and recorded for 2048 turns.

47

Figure 3.9 BPM turn-by-turn data for an arbitrary kick at horizontal BPM R:HP620.

One important example of using TbT data is using it to perform tune measurements. As men-

tioned in Sec. 2.1, specifically in Eq. 2.8, the motion inside the Recycler Ring will exhibit betatron

oscillations. In particular, the tune frequency dictates the main harmonic of these oscillations. Say a

particular set of BPM TbT data has been recorded for 𝑁 number of turns. Therefore, a Fast-Fourier

Transform (FFT) will help uncover and measure FFT accuracy—𝜎𝑄 ≈ 1/𝑁—the main frequency

of these oscillations. Figure 3.10 shows how the prominent peak of the FFT can be identified to

measure the horizontal tune 𝑄𝑥 of the circular accelerator. Furthermore, Ch. 4 will explore how a

more involved Fourier Transform algorithm, such as NAFF (Numerical Analysis of Fundamental

Frequencies) [39], will help analyze the spectrum of TbT data and measure higher order harmonics,

ultimately, leading to the measurement of RDTs in the RR.

48

Figure 3.10 Fast Fourier Transform amplitude for the turn-by-turn data presented in Fig. 3.9.

The second type of diagnostics relevant to this work are the Ion Profile Monitors (IPMs). As

mentioned before, this device measures the beam size in the Recycler Ring. The RR has two IPM

systems: one for the horizontal plane and another one for the vertical. One unique characteristic

of IPMs relies on their non-destructive nature. The working principle of this system is to collect

ions created in the interaction between the proton beam and residual molecules in the vacuum.

Micro-channel plates (MCP) detect the secondary ions and are biased at an arbitrary high voltage

[40]. The counts at each MCP are logged in and digitized to be compatible with ACNET. Similar

to the BPM systems, an application was developed to control and use the IPM system.

49

Figure 3.11 Reconstructed beam profile for the horizontal plane assuming Gaussian distributions
along 1024 turns for an arbitrary beam in the Recycler Ring.

Figure 3.11 shows the reconstruction of the horizontal beam profile from the raw data of the

IPM micro-channel plates. At every turn, the raw count data from each of the 96 MCP is fit to

Gaussian beam distribution, resulting in some variance 𝜎. The color map shows the evolution of

this beam profile along 1024 turns. The IPM system can record up to 65000 turns. Nevertheless,

it batches the data into a 1024 × 96 array. In particular, Fig. 3.11 shows how, for these arbitrary

experimental conditions, the size of the beam grows with the number of turns—this is pointing at

emittance growth. Other distributions different from a Gaussian can be used to fit the raw data to

extract more information about the skewness or tail population of the beam, e.g., a q-Gaussian or

a Skew-Gaussian distribution. In particular, Ch. 6 shows how to use q-Gaussian distributions to

quantify beam tail population.

50

3.6 Sextupoles for Resonance Compensation

The ideal linear lattice of a circular accelerator would be composed only of dipoles and

quadrupoles—assuming the beam is composed of mono-energetic particles. Nevertheless, this

is a strong assumption, given that a bunch of particles will always have an energy spread. Particles

orbiting at different energies will feel different quadrupole fields, resulting in betatron detuning.

This wording is the definition of chromaticity, and therefore, its mathematical definition reads:

𝐶𝑢 =

𝑑 (Δ𝑄𝑢)
𝑑𝛿

.

(3.2)

Every real ring will have a non-zero chromaticity. Sextupoles are crucial to control and correct

chromaticity [1]. Therefore, they are indispensable in any lattice. Nevertheless, as explained

in Ch. 2, introducing non-linear elements inevitably leads to driving betatron resonances, e.g.,

sextupoles will drive third-order resonances. In this sense, and dramatically speaking, sextupoles

are double-edged swords. They are essential for chromaticity correction but deleterious by inducing

resonances in the accelerator.

The Recycler Ring has 76 sextupoles, including normal and skew sextupoles. Additionally, as

described in Sec. 3.2, each gradient magnet has end shims on both ends and is used to correct

the field shape [41]. These end shims can also induce a sextupole-like field on the beam. It is

also worth reminding the reader that the gradient magnets also have sextupole components in their

permanent fields. All of these build up to the fact that the chromaticity of the Recycler Ring can be

controlled in a range of -2 to -20 in both planes, i.e., 𝑄𝑢 ∈ (−2, −20).

Table 3.3 List of sextupoles used to compensate third-order resonances in the Recycler Ring.

Name

Type

SC220(a/b) Normal Sextupole
SC222(a/b) Normal Sextupole
SC319(a/b) Normal Sextupole
SC321(a/b) Normal Sextupole
Skew Sextupole
SS222(a/b)
Skew Sextupole
SS319(a/b)
Skew Sextupole
SS321(a/b)
Skew Sextupole
SS323(a/b)

51

In particular, out of all the sextupoles in the Recycler Ring, this thesis work is interested in

the ones listed in Table 3.3. These elements have been placed in the Recycler Ring with the

sole objective of mitigating the harmful effect of the third-order resonances [37]. The elements

listed in Table 3.3 are in regions of close-to-zero dispersion, 𝐷𝑥 ≈ 0, to cancel out any additional

chromaticity due to their inclusion. In principle, normal sextupoles are used to correct 3𝑄𝑥 = 76

and 𝑄𝑥 + 2𝑄 𝑦 = 74, while skew sextupoles are used to correct 3𝑄 𝑦 = 73 and 2𝑄𝑥 + 𝑄 𝑦 = 75.

Nevertheless, normal sextupoles can have small tilts around the magnetic axis, introducing a

skew-sextupole field and vice versa.

Figure 3.12 shows an example of these sextupoles installed in the Recycler Ring. As can be

seen, two yellow-coated elements appear in the top ring. The same power supply powers both

elements. Hence, the notation in Table 3.3 alludes to a and b elements, e.g., normal sextupole

SC220(a/b). The numbers in the name allude to their location in the RR/MI tunnel, as illustrated

in Fig. 3.4. Given that each pair of elements, a and b, are connected to the same power supply,

these two elements effectively become one knob. The currents for the power supply are controlled

through ACNET. This concept is important for the reader to remember, given that in Ch. 4, the

currents of the sextupoles listed in Table 3.3 will be the main knobs turned in order to mitigate the

strength of these resonances.

52

Figure 3.12 Picture of compensation sextupoles (yellow magnets on top) installed in the Recycler
Ring. Picture provided by Dr. Robert Ainsworth.

53

CHAPTER 4

COMPENSATION OF THIRD-ORDER RESONANCES AT LOW INTENSITIES

4.1 Global RDTs and Lattice Model

The following chapter explores how to mitigate the effect of third-order resonances from the

Recycler Ring by minimizing the Resonance Driving Terms (RDTs) that drive each resonance. The

resonances in question are introduced in Figs. 3.7 and 3.8, and are summarized in Table 4.1. The

RDTs for these third-order resonance lines can be calculated from Eqs. 2.31 and 2.32. Table 4.1

shows the explicit expression for each third-order RDT relevant to this work. The sum over 𝑖 goes

through each element of the lattice beam line and asks if it has some sextupole component in its

definition—it can be normal 𝐾2,𝑖 or skew 𝐾 (𝑠)
2,𝑖

sextupole component. If it has this multipole, it will

add it to the RDT sum by weighting it with the beta functions 𝛽𝑢,𝑖 and phase advances 𝜙𝑢,𝑖 from the

linear approximation at those particular locations. Ultimately, the ℎ 𝑗 𝑘𝑙𝑚 RDT will be a complex

number whose amplitude |ℎ 𝑗 𝑘𝑙𝑚 | should be minimized.

Table 4.1 Corresponding RDTs for each resonance line.

Resonance Line

RDT Expression

3𝑄𝑥 = 76

ℎ3000 = −

𝑄𝑥 + 2𝑄 𝑦 = 74

ℎ1020 = −

1
16

∑︁

𝐾2,𝑖 𝐿𝑖 𝛽

3
2

𝑥,𝑖𝑒3𝑖𝜙𝑥,𝑖

∑︁

𝑖
𝐾2,𝑖 𝐿𝑖 𝛽

1
2

𝑥,𝑖 𝛽𝑦,𝑖𝑒𝑖[𝜙𝑥,𝑖+2𝜙𝑦,𝑖]

3𝑄 𝑦 = 73

ℎ0030 = −

2𝑄𝑥 + 𝑄 𝑦 = 75

ℎ2010 = −

1
16

∑︁

𝑖

∑︁

𝐾 (𝑠)

2,𝑖 𝐿𝑖 𝛽

3
2

𝑦,𝑖𝑒3𝑖𝜙𝑦,𝑖

𝑖
𝐾 (𝑠)
2,𝑖 𝐿𝑖 𝛽𝑥,𝑖 𝛽

1
2

𝑦,𝑖𝑒𝑖[2𝜙𝑥,𝑖+𝜙𝑦,𝑖]

1
48

𝑖

1
48

Figure 4.1 shows a visual representation for calculating the ℎ3000 RDT. It shows the amplitude

of the complex cumulative sum as it goes around the ring (thick, solid orange line). This plot

also shows the amplitude of each contribution for every 𝑖-th element in the lattice with a sextupole

component (thin purple line). This quantity visualizes where and how the sextupole component

is distributed around the ring. Ultimately, after doing this sum around the ring, the final result is

54

a complex number with amplitude and phase corresponding to the ℎ3000 RDT, as calculated from

some arbitrary location in the lattice. The amplitude of the ℎ3000 term is plotted in Fig. 4.1 with a

red dashed line. Figure 4.2 shows a similar exercise for the ℎ1020 term. These calculations are done

with a lattice model with a list of components and magnet coefficients, which should be very close

to what is inside the tunnel. The particular RR model used was the RR2020V0922FLAT lattice

model, provided by R. Ainsworth and M. Xiao.

A question that promptly arises is: does the arbitrary starting position for the sum of Eq. 2.31

change the RDT result? The short answer is yes; the RDT will change depending on the initial

position for the sum. Reference [18], specifically in its Ch. 5, goes into depth on how to correlate

the RDT calculated from a starting point 𝑠1 to one measured at starting point 𝑠2. The difference

in this case relates to the amount of multipole component between both calculation points, e.g.,

the amount of elements that have sextupole component between an 𝑠1 and 𝑠2 observation point.

Nevertheless, given that there is an infinite amount of 𝑠1 and 𝑠2 observation points and only so much

real state in this thesis, the plots are for an arbitrary observation point in the lattice. Additionally,

given that the sextupole components are evenly distributed around the ring, the RDT values will

not oscillate much.

As mentioned in Ch. 3, the Recycler Ring comprises permanent gradient magnets. From

looking at Figs. 4.1 and 4.2, one can see that the sextupole component is evenly distributed around

some sections of the ring. If one were to plot the distribution of permanent gradient on these plots,

their location would coincide with the peaks of the individual contributions to the RDT of Figs. 4.1

and 4.2. Therefore, the sources that drive the Recycler Ring’s normal sextupole resonances come

from the permanent magnets themselves—this is known as a systematic-driven resonance instead

of a random-error-driven resonance. Highly periodic and linear machines use the fluctuations

between RDT measurements from BPMs to locate any sextupole errors in the lattice and try to fix

them [18]. Nevertheless, this is not the case for the Recycler, given its low superperiodicity of 2

and its uniform sextupole component distribution.

55

Figure 4.1 Distribution of the ℎ3000 term around the ring with individual contributions from each
relevant element and the cumulative sum from an arbitrary starting point.

Figure 4.2 Distribution of the ℎ1020 term around the ring with individual contributions from each
relevant element and the cumulative sum from an arbitrary starting point.

56

4.2 Measurement of Third Order RDTs

Calculating the theoretical RDTs from the lattice model is a matter of calculating the sums

outlined in Table 4.1. Nevertheless, measuring the third-order RDTs requires following a long and

involved recipe. This recipe is based on previous work from Refs. [18, 20], but with many original

steps specifically developed for the Recycler Ring. In summary, these are the steps used in order to

measure RDTs at the Recycler:

1. Kick the beam with dipole kickers and save turn-by-turn BPM data for offline analysis.

2. Go through data and estimate momentum coordinates with previously calculated transfer

matrices from the lattice model.

3. Estimate Twiss parameters and normalized coordinates ( ˆ𝑢, ˆ𝑝𝑢) at every BPM location.

4. Create resonance basis ℎ±

𝑢 (see Eq. 2.37).

5. Get spectrum of resonance basis using NAFF (Numerical Analysis of Fundamental Frequen-

cies) through the SUSSIX software [42].

6. Identify resonance lines from the spectrum corresponding to the RDT of interest.

7. Calculate the RDT at each BPM location from the equivalence relation of spectral lines and

RDT expansion (see Eq. 2.46).

The following subsections will explore the steps outlined in the previous list in more detail.

4.2.1 Dipole Kick and BPM data

The first step towards measuring RDTs is to kick (or ping) the beam in either transverse

direction(s) to excite betatron oscillations. Betatron oscillations are the natural oscillations of

particles around their equilibrium orbit in a circular accelerator. Horizontal and vertical dipole

kickers induce a kick in the beam. In particular, the devices used for this were the kicker devices with

ACNET names R:K4XXX and R:KVXXX, horizontally and vertically, respectively. In principle,

these are the Recycler abort kickers—devices used to send the beam to the abort line in Recycler.

57

Nevertheless, these pingers’ high-voltage settings and timings can be changed to give a small kick

to the beam—a ping that is small enough not to steer the beam to the abort line but large enough to

excite betatron oscillations in one or both transverse directions. The beam was kicked exclusively

in the horizontal direction in order to measure purely horizontal RDTs, e.g., ℎ3000, solely in the

vertical direction for vertical RDTs, e.g., ℎ0030, or pinged in both directions for all of them, including

coupling RDTs, e.g., ℎ1020 and ℎ2010.

Once the kickers are set correctly, the next step is to take BPM data. As mentioned in Sec.

3.5, ACNET applications allow gathering and saving BPM data for offline analysis. The BPM data

from all 208 BPMs (104 horizontal and 104 vertical) can be saved in one file. Figure 4.3 shows

an instance of kicking the beam in the horizontal direction and recording beam centroid data for

2048 turns at an arbitrary BPM. The amount of turns recorded is also a customizable quantity. The

ping shown in Fig. 4.3 happens early in the cycle, at around 50 turns. The next hundred turns hold

information about the betatron oscillations. This data is used to extract the tunes 𝑄𝑢 and RDTs of

the machine.

Figure 4.3 BPM data of an arbitrary horizontal kick in the beam at horizontal BPM R:HP620.

58

The physics of a kicked beam is well explained in Refs. [43, 44]. Once the beam is pinged,

the centroid response will oscillate at betatron tunes 𝑄𝑢. Nonlinearities and chromaticity of the

machine will dictate the envelope of these oscillations. How fast this envelope decays is a measure

of the decoherence of the beam. Decoherence of a particle beam in accelerator physics refers to

the process by which a beam that initially has particles oscillating in phase—–meaning they have

similar amplitudes and frequencies—–gradually becomes out of phase over time. This decoherence

effect results in a spread in the particles’ positions and momenta, leading to a more diffuse beam.

This decoherence is caused by nonlinearities in the machine and the transverse chromaticities that

will detune the beam out of coherence, as explained by Eqs. 2.51 and 3.2. Ultimately, this process

will diffuse and maim the signal recorded by the BPMs.

Therefore, when the BPM data for the kicked beam is taken, special care must be taken to sustain

coherent oscillations. This is done by manipulating the machine’s chromaticities using specialized

sextupoles. Table 3.1 showed the nominal chromaticities at which the Recycler Ring operates. For

these RDT measurements, the vertical chromaticity was changed between the range of -7 and -3 to

find sustained vertical oscillations. For the horizontal case, a chromaticity of -5 would usually be

enough for 1 mm oscillations. The other system that affects decoherence is the transverse dampers.

As per its name, these devices dampen any oscillation in the beam. Therefore, they were turned off

for these particular studies. Considering these factors, one could get consequential oscillations in

both transverse planes as a first step to measure Resonance Driving Terms.

4.2.2 Estimation of Momentum Coordinate

At the most fundamental level, BPM data holds only information about the centroid position

of the beam. Nevertheless, in a particle accelerator, it is of interest to look at the whole phase

space picture ( ˆ𝑢, ˆ𝑝𝑢)—including the momentum coordinate. Therefore, explaining how the mo-

mentum coordinate is calculated from the TbT position from every BPM is relevant. The approach

used involves a model-based perspective. Therefore, the lattice model plays a crucial role in the

momentum estimation.

The approach to estimating the momentum coordinate involves solving a least-squares problem.

59

This approach is developed in Ref. [45]. The first step is to calculate the model’s transfer matrices

(linear approximation) from one fixed point in the accelerator to all the horizontal and vertical BPM

locations—in total, there should be 208 transfer matrices. For these studies, the fixed point chosen

was the starting point of the turn count, which is the location of the vertical BPM R:VP601. For

example, the following paragraphs show how to calculate the horizontal phase space coordinates,

including the relative momentum deviation 𝛿. In particular, the main objective is to calculate the
(cid:17)
𝑋0 matrix with dimensions 3 × 2048, which corresponds to the phase space coordinates (cid:16)

(cid:174)𝑥0, (cid:174)𝑥′

0, (cid:174)𝛿0

at the location of R:VP601. The 𝑋0 array will have the following definition:

𝑋0 =

(cid:174)𝑥0

(cid:174)𝑥′
0
(cid:174)𝛿0

(cid:169)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:171)

(cid:170)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:172)

=

(cid:169)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:171)

(cid:2)𝑥′

[𝑥0(𝑁 = 1), 𝑥0(𝑁 = 2), ..., 𝑥0(𝑁 = 2048)]
0(𝑁 = 2048)(cid:3)
0(𝑁 = 2), ..., 𝑥′
[𝛿0(𝑁 = 1), 𝛿0(𝑁 = 2), ..., 𝛿0(𝑁 = 2048)]

0(𝑁 = 1), 𝑥′

,

(cid:170)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:172)

(4.1)

and holds the information over the 2048 turns. The least squares problem is defined as the solution

to the following system:

𝐴𝑋0 = 𝐵,

(4.2)

where 𝐴 is the matrix made from the horizontal coefficients of the model’s transfer matrices, and

it reads:

𝐴 =

(𝑀11

𝑀12

(𝑀11

𝑀12

𝑀13) 𝐵𝑃𝑀 (𝑖−10)

𝑀13) 𝐵𝑃𝑀 (𝑖)

...

...

(𝑀11

𝑀12

𝑀13) 𝐵𝑃𝑀 (𝑖+10)

(cid:169)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:171)

(cid:170)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:172)

.

(4.3)

The notation (...)𝐵𝑃𝑀 ( 𝑗) means that all the matrix coefficients inside the parenthesis are indexed by

the BPM(j) and should be copied from that particular transfer matrix correlating the fixed point to

BPM(j). For this case, the BPM(i) corresponds to the particular BPM location where one wants

to calculate the phase space coordinates 𝑋𝐵𝑃𝑀 (𝑖), after calculating 𝑋0. One can note that only the

ten upstream BPMs and ten downstream BPMs of BPM(i) are included in this calculation. This

60

number is easily customizable in this estimation but should not be too large, i.e., no larger than 30,

to preserve some sense of locality.

One can define the 𝐵 matrix from the BPM observations, and its transpose is just the BPM

responses stacked horizontally. This notation reads explicitly:

...

(cid:69)

...

(cid:68)

𝑥𝐵𝑃𝑀 (𝑖−10)
...

...

𝑥𝐵𝑃𝑀 (𝑖)
...

...

(cid:69)

...

(cid:68)

𝑥𝐵𝑃𝑀 (𝑖+10)
...

.

(cid:69)

(cid:170)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:172)

𝐵𝑇 =

(cid:68)

(cid:169)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:171)

(4.4)

For Eq. 4.4, the triangular bracket notation ⟨...⟩ is used to specify that the BPM data inside the

brackets has already been averaged out—the oscillations recorded in the BPM data is centered

around 0. Again, to use least-squares approximation, this calculation should take the 10 BPMs

upstream and the 10 BPMs downstream of the BPM(i), whose phase space coordinates are being

calculated.

The least-squares solution ˆ𝑋0 to this problem is given by:

ˆ𝑋0 = ( 𝐴𝑇 𝐴)−1 𝐴𝑇 𝐵.

(4.5)

Once ˆ𝑋0 is calculated from the data of 10 BPMs upstream and downstream of BPM(i). The phase
space coordinates ˆ𝑋𝐵𝑃𝑀 (𝑖) at BPM(i) can be calculated from:

ˆ𝑋𝐵𝑃𝑀 (𝑖) =

(cid:174)𝑥

(cid:174)𝑥′

(cid:174)𝛿

(cid:169)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:171)

(cid:170)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:172)𝐵𝑃𝑀 (𝑖)

=

(cid:169)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:171)

[𝑥(𝑁 = 1), ..., 𝑥(𝑁 = 2048)]

[𝑥′(𝑁 = 1), ..., 𝑥′(𝑁 = 2048)]

[𝛿(𝑁 = 1), ..., 𝛿(𝑁 = 2048)]

(cid:170)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:172)

= [𝑀11 𝑀12 𝑀13] 𝐵𝑃𝑀 (𝑖)

ˆ𝑋0.

(4.6)

The hat notation ˆ𝑋 symbolizes data estimated based on the least-squares solution to Eq. 4.5.

Similar to the horizontal case, for the vertical case, the phase space coordinates at each BPM

location can be estimated by first estimating the phase space coordinates 𝑌0 at an arbitrary location

in the lattice—vertical BPM R:VP601 for this case. At this location, the 𝑌0 array will be defined

as:

𝑌0 = (cid:169)
(cid:173)
(cid:173)
(cid:171)

(cid:174)𝑦0

(cid:174)𝑦′
0

(cid:170)
(cid:174)
(cid:174)
(cid:172)

= (cid:169)
(cid:173)
(cid:173)
(cid:171)

[𝑦0(𝑁 = 1), 𝑦0(𝑁 = 2), ..., 𝑦0(𝑁 = 2048)]
0(𝑁 = 2048)(cid:3)

0(𝑁 = 2), ..., 𝑦′

0(𝑁 = 1), 𝑦′

(cid:2)𝑦′

.

(cid:170)
(cid:174)
(cid:174)
(cid:172)

(4.7)

61

The least-squares estimate ˆ𝑌0 for the array in Eq. 4.7 from the recorded BPM data will be given by:

ˆ𝑌0 = ( 𝐴𝑇

𝑦 𝐴𝑦)−1 𝐴𝑇

𝑦 𝐵𝑦,

where similar to the horizontal case, the matrix 𝐴𝑦 is defined by:

𝐴𝑦 =

(𝑀21

(𝑀21

𝑀22) 𝐵𝑃𝑀 (𝑖−10)
...

𝑀22) 𝐵𝑃𝑀 (𝑖)
...

(𝑀21

𝑀22) 𝐵𝑃𝑀 (𝑖+10)

(cid:169)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:171)

,

(cid:170)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:172)

and 𝐵𝑦 is defined by

...

(cid:69)

...

(cid:68)

𝑦𝐵𝑃𝑀 (𝑖−10)
...

...

𝑦𝐵𝑃𝑀 (𝑖)
...

...

(cid:69)

...

(cid:68)

𝑦𝐵𝑃𝑀 (𝑖+10)
...

.

(cid:69)

(cid:170)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:172)

𝐵𝑇

𝑦 =

(cid:68)

(cid:169)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:171)

(4.8)

(4.9)

(4.10)

Once one calculates ˆ𝑌0, it can be transferred to the location of BPM(i) utilizing the transfer matrices.
This is the way of calculating ˆ𝑌𝐵𝑃𝑀 (𝑖), which reads

ˆ𝑌𝐵𝑃𝑀 (𝑖) = (cid:169)
(cid:173)
(cid:173)
(cid:171)

(cid:174)𝑦

(cid:174)𝑦′

(cid:170)
(cid:174)
(cid:174)
(cid:172)𝐵𝑃𝑀 (𝑖)

= (cid:169)
(cid:173)
(cid:173)
(cid:171)

[𝑦(𝑁 = 1), ..., 𝑦(𝑁 = 2048)]

[𝑦′(𝑁 = 1), ..., 𝑦′(𝑁 = 2048)]

(cid:170)
(cid:174)
(cid:174)
(cid:172)

= [𝑀21 𝑀22] 𝐵𝑃𝑀 (𝑖)

ˆ𝑌0.

(4.11)

It is worth pointing out that for the vertical case, one picks out the appropriate elements of the

transfer matrices, i.e., 𝑀21 and 𝑀22 instead of the horizontal coefficients 𝑀11, 𝑀12 and 𝑀13. The

other thing to note is that any momentum dependence is dropped for the vertical case, given that the

vertical dispersion is negligible in the Recycler Ring, as shown in Fig. 3.6. The previous procedure

of calculating ˆ𝑈𝐵𝑃𝑀 (𝑖) is done and saved for each of the 104 horizontal and 104 vertical BPMs—
ˆ𝑋𝐵𝑃𝑀 (𝑖) or ˆ𝑌𝐵𝑃𝑀 (𝑖), accordingly. Figure 4.4 shows an application of this momentum reconstruction
technique for a horizontal BPM R:HP620 and its vertical neighbor R:VP621.

With model-based approaches, it is important to be confident that the model is as close to

the real accelerator as possible. The beta-beating measures how well the model’s beta functions

62

describe the beta functions from the real-world accelerator.

In particular, M. Xiao has shown

that the beta-beating along the Recycler is below 10% [46]—an acceptable quantity for modern

accelerators. Therefore, this proves that the model used is reliable up to some significance level.

The beta-beating quantity is ultimately limited by how linear the accelerator is and any ripple noise

from the power supplies feeding the quadrupole and dipoles.

Figure 4.4 Phase space coordinates reconstruction for two neighboring BPMs—one horizontal
R:HP620 and one vertical R:VP621—windowed for 300 turns.

4.2.3 Twiss Parameters and Normalized Phase Space

The next step to measure RDTs, once the phase space coordinates have been reconstructed for

every BPM, is to build the normal phase space coordinates ( ˆ𝑢, ˆ𝑝𝑢). This is done by using Eqs. 2.9

and 2.13, and the information provided in Fig. 2.1. In order to build the normalized phase space,

the Twiss parameters first have to be estimated. This is done by performing a least-squares fit to

estimate the ellipse parameters of reconstructed phase space, such as the one shown in Fig. 4.4.

Therefore, an ellipse is fit to the data, and after that, the Twiss parameters are retrieved from the fit,

including the centroid action 2𝜋⟨𝐽𝑢⟩ = 𝜀𝑥 (see Eq. 2.9).

Figure 4.5 shows an example of this procedure. The left plot shows reconstructed phase space

data with the best ellipse fit. The parameters indicated on the inset provide detailed characteristics

of this particular fit: (a) 𝛽𝑥 (beta function) describes the spatial spread of the beam, (b) 𝛼𝑥 (alpha

63

function) relates to the angle the beam particles make with the reference orbit, and it is a measure

of the beam divergence, and (c) 𝜀𝑥 (centroid emittance) quantifies the area of the phase space that

the beam centroid occupies. In particular, the centroid action can be calculated from 2𝜋⟨𝐽𝑢⟩ = 𝜀𝑥.

The right plot of Fig. 4.5 shows the reconstructed normalized phase space using the Twiss

parameters, as estimated from the ellipse fit. This normalization process involves scaling by the

square root of the Twiss parameter 𝛽𝑢 at that particular location, per Floquet’s transformation as

defined in Eq. 2.13. This procedure is done for the reconstructed data of every BPM, and, finally

( ˆ𝑢, ˆ𝑝𝑢) is recorded.

Figure 4.5 Reconstructed phase space data with the best ellipse fit (left plot). The right plot shows
the reconstructed normalized phase space using Eq. 2.13 with the Twiss parameters, as estimated
from the ellipse fit. This data corresponds to the horizontal data shown in Fig. 4.4 for the R:HP620
BPM.

4.2.4 Resonance Basis and Spectral Decomposition

Once one reconstructs the normalized phase space coordinates for every BPM, the next step is

to build the resonance basis ℎ±

𝑢 (𝑁), as defined in Eq. 2.37, i.e., ℎ±

𝑢 = ˆ𝑢 ± ˆ𝑝𝑢. This resonance basis

is a function of the number of turns 𝑁 that have elapsed. As shown in Eq. 2.46, this quantity can

have a spectral decomposition as a Fourier series. In general, this spectral decomposition in the

64

horizontal dimension reads:

ℎ−
𝑥 (𝑁) = ˆ𝑥 − 𝑖 ˆ𝑝𝑥 =

and in the vertical dimension:

ℎ−
𝑦 (𝑁) = ˆ𝑦 − 𝑖 ˆ𝑝𝑦 =

∑︁

𝑗 𝑘𝑙𝑚

∑︁

𝑗 𝑘𝑙𝑚

𝐻𝑆𝐿 𝑗 𝑘𝑙𝑚𝑒2𝜋𝑖𝑁 [(1− 𝑗+𝑘)𝑄 𝑥+(𝑚−𝑙)𝑄 𝑦] ,

(4.12)

𝑉 𝑆𝐿 𝑗 𝑘𝑙𝑚𝑒2𝜋𝑖𝑁 [(𝑘− 𝑗)𝑄 𝑥+(1−𝑙+𝑚)𝑄 𝑦] .

(4.13)

These sums run over the ( 𝑗, 𝑘, 𝑙, 𝑚) indices. Theoretically, they run to infinity, but one can truncate

them. The 𝐻𝑆𝐿 𝑗 𝑘𝑙𝑚 term corresponds to the complex amplitude defining the horizontal spectral

line at frequency location (1 − 𝑗 + 𝑘)𝑄𝑥 + (𝑚 − 𝑙)𝑄 𝑦. The 𝑉 𝑆𝐿 𝑗 𝑘𝑙𝑚 term corresponds to the complex

amplitude defining the vertical spectral line at location (𝑘 − 𝑗)𝑄𝑥 + (1 − 𝑙 + 𝑚)𝑄 𝑦. Moreover, as

always, 𝑄𝑢 corresponds to the betatron tunes. These definitions help to understand the experimental

data in order to measure RDTs.

Nevertheless, from the theoretical point of view, the spectral decomposition can be calculated

with the Lie algebra gymnastics shown in Sec. 2.5. These spectral decompositions read for the

horizontal plane:

𝑥 (𝑁) = √︁2𝐼𝑥𝑒𝑖(2𝜋𝑄 𝑥 𝑁+𝜓𝑥0)
ℎ−

− 2𝑖 ∑︁

𝑗 𝑓 𝑗 𝑘𝑙𝑚 (2𝐼𝑥)

𝑗+𝑘−1
2

(cid:0)2𝐼𝑦(cid:1) 𝑙+𝑚

2 𝑒𝑖[(1− 𝑗+𝑘)(2𝜋𝑄 𝑥 𝑁+𝜓𝑥0)+(𝑚−𝑙)(2𝜋𝑄 𝑦 𝑁+𝜓𝑦0)] ,

(4.14)

𝑗 𝑘𝑙𝑚

and for the vertical case:

𝑦 (𝑁) = √︁2𝐼𝑦𝑒𝑖(2𝜋𝑄 𝑦 𝑁+𝜓𝑦0)
ℎ−

− 2𝑖 ∑︁

𝑙 𝑓 𝑗 𝑘𝑙𝑚 (2𝐼𝑥)

𝑗+𝑘
2

(cid:0)2𝐼𝑦(cid:1) 𝑙+𝑚−1

2

𝑒𝑖[(𝑘− 𝑗)(2𝜋𝑄 𝑥 𝑁+𝜓𝑥0)+(1−𝑙+𝑚)(2𝜋𝑄 𝑦 𝑁+𝜓𝑦0)] .

(4.15)

𝑗 𝑘𝑙𝑚

The RDT calculation exploits the equivalence between Eqs. 4.12 and 4.13 to Eqs. 4.14 and 4.15.

Therefore, the generating function coefficients ( 𝑓 𝑗 𝑘𝑙𝑚) can be related to the horizontal and spectral

line coefficients (𝐻𝑆𝐿 𝑗 𝑘𝑙𝑚 and 𝑉 𝑆𝐿 𝑗 𝑘𝑙𝑚). Ultimately, the 𝑓 𝑗 𝑘𝑙𝑚 terms can be related to Resonance

Driving Terms ℎ 𝑗 𝑘𝑙𝑚 through Eq. 2.30.

A quick comparison between Eqs. 4.12 and 4.14 allow the building of an equivalence table

between the generating function coefficients (GFCs) and the spectral lines. Table 4.2 is the result

65

of this. One can find the original table in Refs. [20, 42]. Table 4.2 shows how to calculate the

GFCs from its corresponding spectral line. In particular, the amplitude and phase of a spectral line

will be given by |𝑈𝑆𝐿 𝑗 𝑘𝑙𝑚 | and arg(𝑈𝑆𝐿 𝑗 𝑘𝑙𝑚), while its location in its corresponding frequency

space will be given by 𝑄(𝑈𝑆𝐿 𝑗 𝑘𝑙𝑚). It is worth pointing out that in this notation, 𝑈 can be either

the horizontal or vertical plane.

Table 4.2 Equivalence table between the generating function coefficients and the spectral lines.

Generating Function
Coefficient

Spectral Line

Amplitude

(cid:12)
(cid:12) 𝑓 𝑗 𝑘𝑙𝑚

(cid:12)
(cid:12)

(cid:12)
(cid:12)𝐻𝑆𝐿 𝑗 𝑘𝑙𝑚

(cid:12)
(cid:12) = 2 𝑗 (2𝐼𝑥)

𝑗+𝑘−1
2

(cid:0)2𝐼𝑦(cid:1) 𝑙+𝑚

2

(cid:12)
(cid:12) 𝑓 𝑗 𝑘𝑙𝑚

(cid:12)
(cid:12)

(cid:12)
(cid:12)𝑉 𝑆𝐿 𝑗 𝑘𝑙𝑚

(cid:12)
(cid:12) = 2 𝑙 (2𝐼𝑥)

𝑗+𝑘
2

(cid:0)2𝐼𝑦(cid:1) 𝑙+𝑚−1

2

(cid:12)
(cid:12) 𝑓 𝑗 𝑘𝑙𝑚

(cid:12)
(cid:12)

Phase

𝜙 𝑗 𝑘𝑙𝑚 = arg ( 𝑓 𝑗 𝑘𝑙𝑚)

arg (𝐻𝑆𝐿 𝑗 𝑘𝑙𝑚) = 𝜙 𝑗 𝑘𝑙𝑚 + 𝜓𝑥0 − 𝜋
2

arg (𝑉 𝑆𝐿 𝑗 𝑘𝑙𝑚) = 𝜙 𝑗 𝑘𝑙𝑚 + 𝜓𝑦0 − 𝜋
2

Spectral
Harmonic

N/A

𝑄 (cid:0)𝐻𝑆𝐿 𝑗 𝑘𝑙𝑚(cid:1) = (1 − 𝑗 + 𝑘)𝑄𝑥 + (𝑚 − 𝑙)𝑄 𝑦

𝑄 (cid:0)𝑉 𝑆𝐿 𝑗 𝑘𝑙𝑚(cid:1) = (𝑘 − 𝑗)𝑄𝑥 + (1 − 𝑙 + 𝑚)𝑄 𝑦

4.2.5 Resonance Basis Spectrum

As mentioned in the last section, the resonance basis spectral lines hold enough information

to reconstruct the generating function coefficients (GFCs) and, ultimately, the resonance driving

terms (RDTs). Therefore, it is interesting to have a program that reconstructs the spectral lines from

𝑢 (𝑁) data. This is exactly what SUSSIX [42] does. SUSSIX is a software developed at CERN,
ℎ±

and it utilizes the Numerical Analysis of Fundamental Frequencies (NAFF) algorithm to identify

and calculate the resonance lines up to an arbitrary order. SUSSIX can calculate spectral lines for

the horizontal, the vertical plane, or both planes.

Figure 4.6 shows an example of spectrum data calculated using SUSSIX. The horizontal axis

66

represents the frequency (or tune) of the spectral lines, while the vertical axis shows the amplitude

of the lines, i.e., |𝐻𝑆𝐿 𝑗 𝑘𝑙𝑚 |. All the lines that show up in Fig. 4.6 represent the harmonics present

in ℎ−

𝑥 (𝑁). The largest peak corresponds to the horizontal tune 𝑄𝑥. In this case, the second-largest

tune corresponds to the vertical tune 𝑄 𝑦, indicating some residual coupling. The lines 𝐻𝑆𝐿2010,

𝐻𝑆𝐿3000 and 𝐻𝑆𝐿1020 are also marked with their location at the plot. One can use these lines

to calculate relevant GFCs and corresponding RDTs. Other unmarked spectral lines illustrate

how higher-order harmonics and RDTs might be present in the oscillations. It is especially worth

highlighting how several RDTs can feed to the line at 𝑄𝑥 = 0. Figure 4.7 shows this exercise

done but for the vertical resonance basis ℎ−

𝑦 (𝑁). In this case, the largest peak corresponds to the

vertical tune 𝑄 𝑦, while the second-largest corresponds to the horizontal tune 𝑄𝑥. The lines 𝑉 𝑆𝐿2010,

𝑉 𝑆𝐿3000, and 𝑉 𝑆𝐿1020 also show up, but at different locations than their horizontal counterparts,

as expected. The amplitudes of these spectral lines are at least three orders of magnitude less than

the main harmonic of the oscillations.

Figure 4.6 Spectral lines of ℎ−
the 104 Horizontal BPMs. The spectrum for all BPMs is superimposed in this plot.

𝑥 calculated with SUSSIX [42]. The ℎ−

𝑥 signal was reconstructed for

67

Figure 4.7 Spectral lines of ℎ−
the 104 Vertical BPMs. The spectrum for all BPMs is superimposed in this plot.

𝑦 calculated with SUSSIX [42]. The ℎ−

𝑦 signal was reconstructed for

4.2.6 Spectral Lines and RDT calculation

The spectral lines 𝐻𝑆𝐿3000, 𝐻𝑆𝐿2010, and 𝐻𝑆𝐿1020 identified in Fig. 4.6 are the necessary

lines to calculate the ℎ3000, ℎ2010 and ℎ1020 RDTs. In combination with the 𝑉 𝑆𝐿0030 spectral line

and its corresponding ℎ0030 RDT, these are all the spectral lines needed to calculate the RDTs for

each resonance line shown in Fig. 3.7 and specified in Table 3.2. Table 4.3 shows the explicit

expressions for each resonance line RDT as a function of its GFC. Once one calculates the 𝑓 𝑗 𝑘𝑙𝑚

terms using the expressions in Table 4.2, the RDT is just one multiplication away.

Table 4.3 Corresponding RDTs from GFCs for each resonance line.

Resonance Line

RDT

3𝑄𝑥 = 76

ℎ3000 = 𝑓3000 (cid:0)1 − 𝑒6𝜋𝑖𝑄 𝑥 (cid:1)

𝑄𝑥 + 2𝑄 𝑦 = 74

ℎ1020 = 𝑓1020

(cid:16)1 − 𝑒2𝜋𝑖[𝑄 𝑥+2𝑄 𝑦] (cid:17)

3𝑄 𝑦 = 73

ℎ0030 = 𝑓0030 (cid:0)1 − 𝑒6𝜋𝑖𝑄 𝑦 (cid:1)

2𝑄𝑥 + 𝑄 𝑦 = 75

ℎ2010 = 𝑓2010

(cid:16)1 − 𝑒2𝜋𝑖[2𝑄 𝑥+𝑄 𝑦] (cid:17)

68

4.2.7 Third Order RDTs at every BPM location

Once one transforms the data from spectral lines to GFCs 𝑓 𝑗 𝑘𝑙𝑚, it is just a matter of using

the expressions in Table 4.3 to calculate the RDTs ℎ 𝑗 𝑘𝑙𝑚. Figure 4.8 shows several measurements

of the ℎ3000 term for the Recycler Ring from one particular BPM, i.e., R:HP620. This figure is a

polar plot to indicate the amplitude and phase of this complex quantity. The left plot displays a set

of 36 vectors representing individual measurements of the ℎ3000 term at a specific Beam Position

Monitor (BPM) located at R:HP620. From this plot, one can see some spread in the data from each

measurement. This spread is due to changes in the beam from Booster for every injection and to

noise in BPM data. The right plot shows the average vector of the measurements, with its spread

indicated by the shaded area. This spread signifies the variance in the data from measurement to

measurement.

Figure 4.8 Polar plot for the measurement of ℎ3000 for 36 measurements at horizontal BPM R:HP620.

For every BPM, measurements like those presented in Fig. 4.8 are done. In the case of the

69

ℎ3000, one can use the 104 horizontal BPMs. For the case of the ℎ0030 term, the 104 vertical

BPMs are used. Figure 4.9 shows a plot of the amplitude of the measured ℎ3000 as a function

of the BPM used. Each label (e.g., R:HP100, R:HP106,..., R:HP638) corresponds to a different

BPM in the accelerator. The labels are ordered according to their physical positions along the

Recycler Ring. Each point represents the amplitude of ℎ3000 measured at the corresponding BPM.

The error bars on each data point indicate the uncertainty in the measurement, coming from the

statistical variance as shown in Fig. 4.8. The shaded background of Fig. 4.9 visually emphasizes

the spread and distribution of the data points around the average. While there is a global spread,

the local variations of ℎ3000 from BPM to BPM correspond to the amount of sextupole between

them. Nevertheless, the large discontinuity at the 600 region—before R:HP602—comes from the

momentum-estimation method, which calculates the coordinates at R:VP601. Following all the

previous steps in this recipe will allow us to measure the RDTs at the Recycler Ring.

Figure 4.9 Measurements of ℎ3000 at all horizontal BPMs around the ring. The average (dashed
line) corresponds to the average ℎ3000 from all the measurements.

4.3 Compensation of RDTs

4.3.1 Theoretical Approach from Lattice Model

Once there is a way to measure and calculate RDTs, the next step is to find a way to cancel out

the RDTs. As mentioned in Sec. 3.6, four normal sextupoles and four skew sextupoles have been

installed previously for this particular purpose. They have also been included in the lattice model.

This section focuses on canceling the RDTs from a lattice model perspective—just from theoretical

70

considerations. One can calculate the RDTs from the lattice model as shown in Sec. 4.1.

The first step is to understand how to cancel out a single RDT—a single resonance line. For

example, to compensate the 3𝑄𝑥 line, the ℎ3000 would need to be canceled.

In order to cancel

this term, one needs to set the appropriate 𝐾2 coefficients to the sextupoles. This procedure is

equivalent to setting the appropriate currents in the sextupoles. The appropriate 𝐾2 coefficients are

calculated from the following equation:

𝐾 (𝑠𝑐220)

−|ℎ3000| sin 𝜓3000

−|ℎ3000| cos 𝜓3000













 (𝐵𝑎𝑟𝑒)
where the ℎ3000 RDT is a complex quantity with a real and imaginary part, i.e., ℎ3000 = |ℎ3000| 𝑒𝑖𝜓3000 =

𝐾 (𝑠𝑐321)

𝐾 (𝑠𝑐319)

𝐾 (𝑠𝑐222)











































(4.16)

= 𝑴

0

0

2

2

2

2

,

|ℎ3000| cos 𝜓3000 + 𝑖|ℎ3000| sin 𝜓3000. The 𝑴 term corresponds to the RDT response matrix. Fi-
nally, the 𝐾 (𝑖)

2 corresponds to the 𝑖-th compensation sextupole coefficient. The response matrix for
compensating only 3𝑄𝑥 and arbitrarily setting two sextupole currents to the negative of each other

reads:

𝑴 =















𝑀11

𝑀21

1

0

𝑀12

𝑀22

1

0

𝑀13

𝑀23

0

1

𝑀14

𝑀24

0

1















,

(4.17)

where the theoretical matrix coefficients are calculated from the derivatives of the expressions in

Table 4.1 concerning 𝐾 (𝑖)

2 , and they read:

𝑀11 =

𝑀12 =

𝑀13 =

𝑀14 =

1
48

1
48
1
48
1
48

𝑙𝑠𝑐 𝛽3/2

𝑥(𝑠𝑐220)

cos 3𝜙𝑥(𝑠𝑐220),

𝑙𝑠𝑐 𝛽3/2

𝑥(𝑠𝑐222)

cos 3𝜙𝑥(𝑠𝑐222),

𝑙𝑠𝑐 𝛽3/2

𝑥(𝑠𝑐319)

cos 3𝜙𝑥(𝑠𝑐319),

𝑙𝑠𝑐 𝛽3/2

𝑥(𝑠𝑐321)

cos 3𝜙𝑥(𝑠𝑐321),

71

(4.18)

(4.19)

(4.20)

(4.21)

𝑀21 =

𝑀22 =

𝑀23 =

𝑀24 =

1
48
1
48
1
48
1
48

𝑙𝑠𝑐 𝛽3/2

𝑥(𝑠𝑐220)

sin 3𝜙𝑥(𝑠𝑐220),

𝑙𝑠𝑐 𝛽3/2

𝑥(𝑠𝑐222)

sin 3𝜙𝑥(𝑠𝑐222),

𝑙𝑠𝑐 𝛽3/2

𝑥(𝑠𝑐319)

sin 3𝜙𝑥(𝑠𝑐319),

𝑙𝑠𝑐 𝛽3/2

𝑥(𝑠𝑐321)

sin 3𝜙𝑥(𝑠𝑐321),

(4.22)

(4.23)

(4.24)

(4.25)

where each beta function 𝛽𝑢(𝑖) and phase evaluation 𝜙𝑢(𝑖) in the 𝑢 plane—horizontal or vertical—is

taken at the location of the 𝑖-th compensation element. Every compensation magnet has the same

length of 𝑙𝑠𝑐.

Once one calculates the response matrix, the solution can be found by calculating its inverse

and operating on the bare machine RDT vector. This RDT vector is a known quantity that can be

calculated from the lattice model. The solution for compensation thus reads:

𝐾 (𝑠𝑐220)

2

2

𝐾 (𝑠𝑐222)

−|ℎ3000| cos 𝜓3000













 (𝐵𝑎𝑟𝑒)
Table 4.4 shows this equation’s solution for arbitrary lattice tunes.













 (𝐶𝑜𝑚 𝑝)

−|ℎ3000| sin 𝜓3000

𝐾 (𝑠𝑐321)

𝐾 (𝑠𝑐319)

= 𝑴−1





























0

0

2

2

.

(4.26)

It is worth pointing out

that the RDT response matrix, as defined in Eq. 4.17, should have an inverse. Nevertheless, one

can relax the problem constraints by dropping the two last rows of Eq. 4.17. For this case, the

corresponding equation to Eq. 4.26 would involve calculating the pseudo-inverse using a least-

squares approach. In this case, every sextupole should have different values for the 𝐾2 coefficients

that compensate ℎ3000.

72

Table 4.4 𝐾2 coefficients that cancel out the ℎ3000 term. These values were calculated for a Recycler
lattice tuned to 𝑄𝑥 = 25.4422 and 𝑄 𝑦 = 24.3915.

Sextupole 𝐾2 [m−3]

SC220

-0.594702

SC222

0.594702

SC319

0.930019

SC321

-0.930019

Figure 4.10 shows the calculation of the ℎ3000 RDT once one sets the 𝐾2 coefficients in the

sextupoles to the values in Table 4.4. For this case, as opposed to the bare machine calculation from

Fig. 4.1, the cumulative sum of ℎ3000 goes to 0 at the end of the sum. This compensation means

that one has canceled the global RDT. Furthermore, if Fig. 4.1 is compared to Fig. 4.10, one can

see the big spikes in the middle of the lattice that corresponds to the individual contributions from

the compensation sextupoles. These kicks ultimately bring the ℎ3000 term to zero.

73

Figure 4.10 Distribution of the ℎ3000 term around the ring with individual contributions from each
relevant element and the cumulative sum from an arbitrary starting point when correction elements
are set to compensate 3𝑄𝑥 = 76, i.e., ℎ3000 = 0.

On the other hand, there is the ℎ1020 term. This one also gets affected if normal sextupole

components are changed in the lattice, as explained by the expressions in Table 4.1. Therefore,

looking at how ℎ1020 changes as ℎ3000 is compensated is interesting. Figure 4.11 shows the

calculation for this term as one compensates 3𝑄𝑥. Again, four additional spikes show up in the

individual contribution plot exactly where the compensation sextupoles are placed. If one compares

this plot to Fig. 4.2, it can also be seen that the total sum of the ℎ1020 term increases from 1.739

to 1.988 m−1/2. This shows that by compensating the ℎ3000 RDT, the ℎ1020 increases. This effect

means that compensating one resonance line might worsen other lines.

74

Figure 4.11 Distribution of the ℎ1020 term around the ring with individual contributions from each
relevant element and the cumulative sum from an arbitrary starting point when correction elements
are set to compensate 3𝑄𝑥 = 76, i.e., ℎ3000 = 0.

An effort can be made to extend Eq. 4.16 to include the cancellation of ℎ1020. The simultaneous

compensation of the ℎ3000 and the ℎ1020 term is given by the following equation:

−|ℎ3000| cos 𝜓3000

−|ℎ3000| sin 𝜓3000

−|ℎ1020| cos 𝜓1020

−|ℎ1020| sin 𝜓1020



























 (𝐵𝑎𝑟𝑒)

= 𝑴

𝐾 (𝑠𝑐220)

2

𝐾 (𝑠𝑐222)

2

𝐾 (𝑠𝑐319)

2

𝐾 (𝑠𝑐321)

2





























.

For this case, one needs to modify the response matrix 𝑴. This new response matrix reads:

𝑴 =

𝑀11

𝑀21

𝑀31

𝑀41















𝑀13

𝑀23

𝑀33

𝑀43

,

𝑀14

𝑀24

𝑀34

𝑀44















𝑀12

𝑀22

𝑀32

𝑀42

75

(4.27)

(4.28)

where the first two rows are defined the same as before by Eqs. 4.19-4.25. Similarly, the coefficients

in the last two rows are defined by:

𝑀31 =

𝑀32 =

𝑀33 =

𝑀34 =

𝑀41 =

𝑀42 =

𝑀43 =

𝑀44 =

1
16

1
16
1
16
1
16
1
16
1
16
1
16
1
16

𝑙𝑠𝑐 𝛽1/2

𝑥(𝑠𝑐220)

𝛽𝑦(𝑠𝑐220) cos (cid:0)𝜙𝑥(𝑠𝑐220) + 2𝜙𝑦(𝑠𝑐220)

(cid:1) ,

𝑙𝑠𝑐 𝛽1/2

𝑥(𝑠𝑐222)

𝑙𝑠𝑐 𝛽1/2

𝑥(𝑠𝑐319)

𝑙𝑠𝑐 𝛽1/2

𝑥(𝑠𝑐321)

𝑙𝑠𝑐 𝛽1/2

𝑥(𝑠𝑐220)

𝑙𝑠𝑐 𝛽1/2

𝑥(𝑠𝑐222)

𝑙𝑠𝑐 𝛽1/2

𝑥(𝑠𝑐319)

𝑙𝑠𝑐 𝛽1/2

𝑥(𝑠𝑐321)

𝛽𝑦(𝑠𝑐222) cos (cid:0)𝜙𝑥(𝑠𝑐222) + 2𝜙𝑦(𝑠𝑐222)

(cid:1) ,

𝛽𝑦(𝑠𝑐319) cos (cid:0)𝜙𝑥(𝑠𝑐319) + 2𝜙𝑦(𝑠𝑐319)

(cid:1) ,

𝛽𝑦(𝑠𝑐321) cos (cid:0)𝜙𝑥(𝑠𝑐321) + 2𝜙𝑦(𝑠𝑐321)

(cid:1) ,

𝛽𝑦(𝑠𝑐220) sin (cid:0)𝜙𝑥(𝑠𝑐220) + 2𝜙𝑦(𝑠𝑐220)

(cid:1) ,

𝛽𝑦(𝑠𝑐222) sin (cid:0)𝜙𝑥(𝑠𝑐222) + 2𝜙𝑦(𝑠𝑐222)

(cid:1) ,

𝛽𝑦(𝑠𝑐319) sin (cid:0)𝜙𝑥(𝑠𝑐319) + 2𝜙𝑦(𝑠𝑐319)

(cid:1) ,

𝛽𝑦(𝑠𝑐321) sin (cid:0)𝜙𝑥(𝑠𝑐321) + 2𝜙𝑦(𝑠𝑐321)

(cid:1) ,

(4.29)

(4.30)

(4.31)

(4.32)

(4.33)

(4.34)

(4.35)

(4.36)

as calculated from the derivative of ℎ1020 with respect to 𝐾 (𝑖)

2 , with the explicit expression in Table

4.1.

The solution to Eq. 4.27 is just given by the inverse of the response matrix applied to the bare

machine RDT vector. This perspective is equivalent to the approach in Eq. 4.26. The result of this

approach is the calculation of the 𝐾2 coefficients in the compensation sextupoles that cancel out

both the ℎ3000 and the ℎ1020 term. The solution to Eq. 4.27 is given in Table 4.5 for an arbitrarily

tuned lattice. One thing to note is that these values are almost one order of magnitude larger than

those in Table 4.4.

76

Table 4.5 𝐾2 coefficients that cancel out the ℎ3000 and ℎ1020 term. These values were calculated for
a Recycler lattice tuned to 𝑄𝑥 = 25.4422 and 𝑄 𝑦 = 24.3915.

Sextupole 𝐾2 [m−3]

SC220

6.039888

SC222

4.826964

SC319

-4.529674

SC321

-6.167827

Once one inputs the magnet coefficients from Table 4.5 into the lattice model, the RDT calcu-

lations can be done to verify that the terms effectively go to 0. Figures 4.12 and 4.13 show these

calculations. Moreover, one can see how the RDT sums add up to 0. Nevertheless, it can be seen

from these plots how the spikes are now very large compared to the previous cases and compared

to the individual contributions from around the lattice. This effect means the sextupole component

will be localized in the compensation sextupoles’ locations. This is an unwanted effect, given that

these large sextupole kicks will overfocus the beam in these locations. In particular, the sextupole

component will no longer be perturbative in these locations.

77

Figure 4.12 Distribution of the ℎ3000 term around the ring with individual contributions from each
relevant element and the cumulative sum from an arbitrary starting point. This sum is done with
existing sextupoles powered at the correct currents to simultaneously cancel out ℎ3000 and ℎ1020.

The results from solving Eq. 4.27, as specified in Table 4.5, indicate that the 𝐾2 coefficients

are very large compared to the case where there was only ℎ3000 compensation. Experimentally, the

currents in the sextupoles that compensate both ℎ3000 and ℎ1020 are too large and exceed the range

of operational limits. While there are no limits for 𝐾2 in the sextupoles from the lattice model,

the solutions are very large. Suppose one were to translate the operational current limits of the

sextupoles to the lattice model 𝐾2 coefficients. In that case, one can show that the solutions lie out

of this range, i.e., the maximum 𝐾2 coefficients are around 1 m−3. Therefore, this solution does not

work for the operational beam. Section 4.6 shows how to circumvent this by adding two additional

sextupoles in the lattice at optimized locations to effectively bring down the compensation 𝐾2

coefficients—the compensation currents.

78

Figure 4.13 Distribution of the ℎ1020 term around the ring with individual contributions from each
relevant element and the cumulative sum from an arbitrary starting point. This distribution is with
existing sextupoles powered at the correct currents to simultaneously cancel out ℎ3000 and ℎ1020.

While the previous section focused only on the normal sextupole RDTs, all of this can easily be

extended to the skew sextupole RDTs, i.e., ℎ0030 and ℎ2010. All of this assumes the lattice model has

a good approximation of the skew sextupole distribution around the ring. For the Recycler Ring,

the skew sextupole resonances are not as strong as the normal sextupole resonance 3𝑄𝑥 = 76. This

resonance dominates the losses in this region. Therefore, special attention has been paid to the

3𝑄𝑥 = 76 resonance and its ℎ3000 resonance driving term. Nevertheless, using the skew sextupole

correctors, this procedure works for 3𝑄 𝑦 = 73 and 2𝑄𝑥 + 𝑄 𝑦 = 75 and its corresponding RDTs.

4.3.2 Response Matrix Approach

The next step is translating the procedure from the previous section to the real machine. This step

is done by measuring the response matrix. As mentioned before, for resonance compensation there

are four dedicated normal sextupoles with currents that can be set to (𝐼𝑠𝑐220, 𝐼𝑠𝑐222, 𝐼𝑠𝑐319, 𝐼𝑠𝑐321)

and four dedicated skew sextupoles with currents that can be set to (𝐼𝑠𝑠323, 𝐼𝑠𝑠323, 𝐼𝑠𝑠319, 𝐼𝑠𝑠321). As

79

shown in the previous section, one RDT can be canceled out with the right kick from the correction

elements, which means the resonances are corrected to the first order.

Nevertheless, compensating for one resonance line might worsen other resonances. This is

why, for simultaneous compensation, compensation currents will vary depending on the subsets

of resonances to compensate. In principle, the currents 𝐼𝑥 needed in each correction element in

order to cancel out the four bare machine RDTs are given by the solution to this linear system of

equations:

−|ℎ3000| cos(𝜓3000)

−|ℎ3000| sin(𝜓3000)

−|ℎ1020| cos(𝜓1020)

−|ℎ1020| sin(𝜓1020)

−|ℎ0030| cos(𝜓3000)

−|ℎ0030| sin(𝜓3000)

−|ℎ2010| cos(𝜓1020)

−|ℎ2010| sin(𝜓1020)





















































 (𝐵𝑎𝑟𝑒)

= 𝑴

𝐼𝑠𝑐220

𝐼𝑠𝑐222

𝐼𝑠𝑐319

𝐼𝑠𝑐321

𝐼𝑠𝑠223

𝐼𝑠𝑠323

𝐼𝑠𝑠319

𝐼𝑠𝑠321























































(4.37)

where 𝑀𝑖 𝑗 is the response matrix for the RDTs with respect to the currents and includes any roll

that can happen for the correction sextupoles. It is worth emphasizing that one can calculate the

response matrix from measurements instead of the model for this approach. This response matrix

𝑀𝑖 𝑗 can be calculated by scanning the currents in each correction element and looking at the

response from the real and imaginary part of the RDTs, i.e., ℎ 𝑗 𝑘𝑙𝑚 = |ℎ 𝑗 𝑘𝑙𝑚 |𝑒𝑖𝜓 𝑗 𝑘𝑙𝑚.

In reality, there are limitations to solving Eq. 4.37. First, all the RDTs (ℎ 𝑗 𝑘𝑙𝑚) may not be

accessible for measurement, given that they may not show up as a spectral line. Another limitation

is that the solution for the currents may be outside the maximum limits for the correction elements.

One can also try canceling a subset of RDTs from Eq. 4.37, including only one RDT. For

example, in order to compensate 3𝑄𝑥 = 76, the system of equations to be solved is:

80

−|ℎ3000| sin(𝜓3000)


−|ℎ3000| cos(𝜓3000)













0

0













 (𝐵𝑎𝑟𝑒)

=















𝑀11 𝑀12 𝑀13 𝑀14

𝑀21 𝑀22 𝑀23 𝑀24

1

0

1

0

0

1

0

1















𝐼𝑠𝑐220

𝐼𝑠𝑐222

𝐼𝑠𝑐319

𝐼𝑠𝑐321





























,

(4.38)

which is analogous to Eq. 4.16, only now the sextupole currents are the knobs. Nevertheless, in

this case, the RDT is measured from the procedure shown in Sec. 4.2 for the bare machine—no

compensation sextupoles turned on. In order to measure the response matrix 𝑴, the ℎ3000 RDT

has to be measured as a function of the compensation currents in the sextupoles. The procedure

summarizing all of this is the following:

1. Measure the bare machine RDT vector for several shots, e.g., measure ℎ3000 for 20 instances.

2. Choose a compensation sextupole for the first scan.

3. Set the compensation sextupole to the lower limit of the scan range, e.g., -10 Amps.

4. Measure the ℎ3000 RDT (or arbitrary RDT) using following the recipe in Sec. 4.2.

5. Increase the current in the compensation sextupole by some step size, e.g., by 2 Amps, and

repeat the RDT measurement.

6. Repeat steps 2-4 until the upper limit of the scan range, e.g., 10 Amps.

7. For every measurement, decompose the RDT into its real and imaginary part, e.g., ℎ3000 =

|ℎ3000| cos 𝜓3000 + 𝑖|ℎ3000| sin 𝜓3000.

8. Plot the real and imaginary part of the RDT as a function of sextupole current and perform a

linear fit. The slopes of these fits are the corresponding 𝑀𝑖 𝑗 response matrix coefficients.

9. Repeat steps 2-8 for the other compensation sextupoles until the response matrix is fully

calculated. This procedure can be done for a subset of correctors, reducing the dimension of

the response matrix.

81

10. Once the response matrix is specified, calculate its inverse (or pseudo-inverse) and operate

on to the bare machine RDT vector in order to get the compensation currents on every BPM

(see Fig. 4.16).

Figures 4.14 and 4.15 show an example of performing two of these scans for the ℎ3000 term. In

Fig. 4.14, it can be seen how the relationship between the real part and the imaginary part follows a

linear relation within the shaded region. The fit for these plots was done only for this region, given

that prior empirical data showed that the currents should not exceed the absolute value of 7 Amps.

One can extract the response matrix coefficients from the slope values provided by the fit.

Figure 4.14 Real and imaginary part of ℎ3000 as a function of current fed to compensation sextupole
SC222. Measurement is shown for R:HP126 data.

Figures 4.14 and 4.15 are taken at two different BPMs, i.e., horizontal BPMs R:HP126 and

R:HP612. Comparing both plots shows that the RDT sensitivity of ℎ3000 changes depending on the

sextupole used. Furthermore, at high currents, it can be seen that some non-linear behavior starts

to kick in. This behavior was present in almost all BPMs. Nevertheless, this only happened at high

currents far from the region of interest, which was capped at an absolute value of 6 Amps.

After performing this process, the result is an array of compensating currents predicted by every

BPM. Figure 4.16 shows an example of these results. The bars at each BPM represent the predicted

compensation current for a particular sextupole. While each BPM will give a local approximation

82

of what currents minimize the ℎ3000 RDT, the sextupoles can only take one input. One must

average all these currents to get a global compensation current that minimizes the RDT at every

location. Nevertheless, given some uncertainty on the bare machine RDTs and the response matrix

coefficients, this needs to be included in developing this global solution.

Figure 4.15 Real and imaginary part of ℎ3000 as a function of current fed to compensation sextupole
SC319. Measurement is shown for R:HP612 data.

In order to build a final setting for the compensation currents, not all BPMs are necessarily used.

For this prediction, the BPMs that showed the most linear fit from the RDT scan as a function of

current were used. One can use the R-squared statistics and the reduced 𝜒2 to quantify this linearity.

Once the best BPMs were selected, the compensation currents were selected. This selection is why

not all BPMs are plotted in Fig. 4.16. One can build a prediction function assuming that the

currents predictions shown in Fig. 4.16 follow a bi-Gaussian distribution centered at the mean

values with a variance equal to the uncertainty in the error bar. The prediction function is just the

sum of all the bi-Gaussian distributions. A contour plot example of such a prediction function is

shown in Fig. 4.17. Ultimately, this prediction function will specify a region in the currents’ space

of each sextupole where the best setting lies.

83

Figure 4.16 Compensation currents calculated from BPMs that showed the best R-square statistic
from the RDT scans.

84

Figure 4.17 Prediction function built from the best R-squared statistics BPMs to predict the best
global currents for ℎ3000 compensation.

The prediction function in Fig. 4.17 can be used to predict an optimal setting for the currents

in the compensation sextupoles. Specifically, the centroid of this prediction distribution can be

calculated and used as the optimum setting. This centroid corresponds to the cyan cross in Fig. 4.17.

Later in this chapter, a whole section is dedicated to verifying the compensation of resonances; see

Sec. 4.5. Nevertheless, a quick way to verify resonance compensation is to use a dynamic tune ramp

that crosses the 3𝑄𝑥 = 76. Uncorrected for, around 95% of the beam is lost once this resonance

is crossed. One can use the interval where the prediction function thinks it is best to operate and

scan the parameter space (sextupole currents) while recording beam transmission across the 3𝑄𝑥

resonance. The fuchsia X in Fig. 4.17 shows the setting with the best transmission—the best

transmission now corresponds to less than 5% losses. The centroid of the prediction function is

around 2.4% (relative distance) away from this experimental value. This difference shows that the

predicted value from the response matrix approach does indeed reduce the losses from the 3𝑄𝑥 = 76

line—meaning the ℎ3000 RDT is minimized in this region.

85

4.4 Optimization of Compensation Currents

The last section showed how to predict the optimum setting for the sextupoles by measuring

the response matrix. Nevertheless, this procedure takes a few hours of study time and might not

be accessible for routine operations. Therefore, one can perform this procedure once or twice per

run. It is interesting to explore optimization algorithms that can help predict the optimum setting

that minimizes losses from the resonance lines.

Figure 4.18 Nelder-Mead optimization from an instance where the response matrix prediction was
used as the initial point. The contour plot shows the transmission through 3𝑄𝑥 as a function of the
currents.

Figure 4.18 shows an example of using a Simplex (Nelder-Mead) optimization procedure to find

the optimum setting in the compensation sextupoles. In the background of this plot, the contour

plot corresponds to the experimentally measured transmission through 3𝑄𝑥. This Nelder-Mead

instance was started with the currents predicted from the response matrix approach. After 45

iterations, the optimization algorithm converged to the value specified in the legend of Fig. 4.18.

This value is fairly close to the value predicted by the response matrix approach. The maximum

transmission from the contour plot (red cross in Fig. 4.18) agrees up to some significance with

the currents from the optimization procedure and the response matrix approach. They are all fairly

86

close within a region. Therefore, by taking the response-matrix predicted currents as an initial

point, one can fine-tune this configuration using this optimization algorithm.

In principle, any

numerical optimization algorithm should converge to the answer given that the underlying function

is not especially difficult—in the context of numerical optimization.

Figure 4.19 Multiobjective Nelder-Mead optimization of the transmission through 3𝑄𝑥 and 3𝑄 𝑦.
The objective function to be minimized was defined as an equally weighted sum of both minus
transmissions—corresponding to the fractional losses through the resonance lines.

Figure 4.19 takes this further by performing multiobjective optimization on two resonance

lines, i.e., 3𝑄𝑥 and 3𝑄 𝑦. The eight compensation sextupoles (see Table 3.3), including normal and

skew, were used for this optimization. While, in principle, powering up the normal compensation

sextupoles should not affect the 3𝑄 𝑦 line, it was observed that the losses due to this line increased

when 3𝑄𝑥 was corrected. This observation motivated the multiobjective optimization of both

lines. Figure 4.19 shows how the other one got worse by trying to correct one line. Ultimately,

87

the best configuration would lie on the Pareto front of these functions. The y-axis in Fig. 4.19

represents the fractional beam loss when crossing the resonance lines. The objective function

minimized by a Nelder-Mead algorithm was defined as an equally weighted sum of both fractional

losses. Ultimately, Fig. 4.19 shows how the algorithm slowly drifts toward the best solution to this

problem.

Figure 4.20 Bayesian optimization of the compensation currents for 3𝑄𝑥 transmission. The contour
plot shows the Gaussian Process Posterior function trained by the sampled data (orange dots).

One can use more advanced optimization algorithms to fine-tune the best optimization currents

in the compensation sextupoles. Figure 4.20 shows an instance where Bayesian optimization was

used for this purpose. A summary of Bayesian optimization can be found in Ref.

[47].

It is

worth clarifying that the optimum currents from Fig. 4.20 differ from the ones in 4.18 given

that the first plot was generated from data taken in the 2022-2023 run and the other one from the

2021-2022 run. The tuning of the machine changes from year to year—even from week to week.

The background contour plot of Fig. 4.20 corresponds to the Gaussian Process posterior trained

from the Bayesian-guided transmission sampling through 3𝑄𝑥. Ultimately, the cluster of points

around the maximum shows how the algorithm converged and found the best configuration to feed

88

the sextupoles. These currents are close to the ones shown in Fig. 4.17. Chapter 5 shows some

studies done at the Proton-Synchrotron Booster (PSB) at CERN, where this optimization approach

is fully embraced to compensate resonance lines.

4.5 Experimental Verification of Compensation

4.5.1 Dynamic Loss Maps

One effective method for visualizing resonance compensation involves constructing dynamic

loss maps. In order to generate these loss maps, specialized quadrupoles responsible for controlling

the tune of the Recycler are gradually adjusted to map out the desired tune area. These are

quadrupoles inside the tune trombone region. The beam loss rate is measured using DCCT devices

across the specified region throughout this process. This procedure is done in the horizontal and

vertical directions. One generates the initial horizontal scan by maintaining a constant vertical tune

while implementing a horizontal tune ramp ranging from 𝑄𝑥 = 25.47 to 𝑄𝑥 = 25.31. Subsequently,

the vertical tune, initially set as constant at 𝑄 𝑦 = 24.47, is adjusted incrementally to 𝑄 𝑦 = 24.31 in

steps of 0.005, with intensity data recorded at each step. Conversely, the roles are reversed for the

vertical scan: the horizontal tune remains constant while a vertical tune ramp progresses from 𝑄 𝑦 =

24.47 to 𝑄 𝑦 = 24.31. Then, the constant horizontal tune is varied from 𝑄𝑥 = 25.47 to 𝑄 𝑦 = 25.31

in steps of 0.005. The resulting intensity data from both scans can be differentiated, normalized by

the instantaneous intensity, and interpolated within a two-dimensional grid to construct plots akin

to those depicted in Figs. 4.21 and 4.22. Figure 4.21 demonstrates the initial machine scan without

compensation. If plotted alongside the theoretical positions of the lines as in Fig. 4.22, the beam

loss bands align with the resonance lines.

Figure 4.22 illustrates the correspondence between the loss patterns and the theoretical positions

of resonance lines. A slight deviation exists between the set and actual tunes due to calibration

adjustments from the tune trombone program. However, despite this variance, the resonance line

configuration within the loss pattern facilitates the visualization of each resonance’s strengths.

Specifically, a higher normalized loss at a particular tune location indicates a stronger Resonance

Driving Term (RDT) for the corresponding resonance line.

Inside Figs. 4.21 and 4.22, third,

89

fourth, and even traces of fifth-order resonance lines are discernible, with third-order resonance

lines exhibiting the greatest prominence.

Figure 4.21 Dynamic loss map from ramping the tunes with an interval of Δ𝑄𝑢 = 0.005 in both
directions. The directions of the scan are from left to right and top to bottom. The results are
superimposed in this plot.

Figure 4.23 depicts dynamic loss maps representing various configurations of the compensation

sextupoles. Specifically, Fig. 4.23a illustrates the loss map for the bare machine, where no

compensation sextupoles are activated, while Fig. 4.23b and Fig. 4.23c demonstrate compensation

for a single resonance line each. For 3𝑄𝑥 compensation, one can set the four normal sextupoles

to the calculated compensation currents using the RDT response matrix method. Moreover, a

comparison between Fig. 4.23b and Fig. 4.23a indicates a reduction of normalized losses by two

orders of magnitude at the 3𝑄𝑥 line with compensation. This observation also holds for the 3𝑄 𝑦

90

compensation, as shown in Fig. 4.23c.

Figure 4.22 Dynamic loss map with the corresponding lines from Fig. 3.7 drawn on top.

Figures 4.23d, 4.23e, and 4.23f showcase the optimal configurations of compensation sextupoles

designed to address multiple resonance lines simultaneously.

It is important to note that while

attempting to compensate for one or multiple resonance lines, there is a possibility that other

resonance lines may strengthen. This effect is evident in the explicit case depicted in Fig. 4.23f,

where compensating for 3𝑄 𝑦 and 𝑄𝑥 + 2𝑄 𝑦 leads to the amplification of the 2𝑄𝑥 + 𝑄 𝑦 resonance.

Such occurrences pose a limitation when aiming to compensate for more than two resonance lines,

as the compensation currents tend to increase. There exists a constraint on the currents supplied

to the compensation sextupoles. For instance, the required currents exceed the current limits when

compensating both normal sextupole lines, 3𝑄𝑥 and 𝑄𝑥 + 2𝑄 𝑦. Ongoing efforts are focused on

reducing the compensation currents in this specific scenario. Section 4.6 summarizes some of these

91

efforts, where additional sextupoles have been installed in order to decrease the currents that cancel

out both the ℎ3000 and ℎ1020 RDTs.

Another notable detail is evident in Figs. 4.23a-4.23f is the presence of white areas within the

loss maps, indicating regions with insufficient beam particles to map out the losses accurately. In

certain configurations of the compensation sextupoles, the combined weakening of the third-order

resonance lines occurs in a manner that leaves some beam remaining beyond these lines. One can

argue that conducting two additional scans, injected from the left and bottom, could effectively map

out these inaccessible regions. Such an enhancement could be considered a future upgrade to these

dynamic loss maps. Ultimately, all the plots presented in Fig. 4.23 demonstrate various potential

configurations that open up regions of tune space for utilization during operations, enabling the

accommodation of high-intensity beams.

92

(a) Bare machine.

(b) 3𝑄 𝑥 Compensation.

(c) 3𝑄 𝑦 Compensation.

(d) 3𝑄 𝑥 and 3𝑄 𝑦 Compensation.

(e) 3𝑄 𝑥 and 2𝑄 𝑥 + 𝑄 𝑦 Compensation.

(f) 3𝑄 𝑦 and 𝑄 𝑥 + 2𝑄 𝑦 Compensation.

Figure 4.23 Dynamic loss maps for several configurations of compensation sextupoles.

93

4.5.2 Static Tune Scans

Another method for visualizing resonance compensation involves using static tune scans. While

the loss maps detailed earlier illustrate the dynamic crossing of resonances, an alternative method

entails setting the tune to a specific value and assessing the beam survival ratio alongside the beam

size over a defined time interval. This time interval was around 0.8 seconds, which is equivalent to

approximately 72800 turns. The beam survival ratio was calculated from the beam intensity data

captured from the DC Current Transformer (DCCT). The beam size is quantified using the Ion

Profile Monitor System (IPM) introduced in Sec. 3.5.

Figure 4.24 Static tune scan for the bare machine with comparisons between experimental data
and SYNERGIA simulations at a low intensity of 0.5e10 particles per bunch (2 Booster Turns of
equivalent intensity).

Figure 4.24 shows a static tune scan only including beam survival ratio data. This plot includes

94

data from experiments and data from SYNERGIA simulations done with no space charge. No

space-charge simulator was needed because these experiments were done at low intensities. The

simulations were done for 72800 turns, and the experimental data was partitioned to only this

equivalent time window. Figure 4.24 shows the resonance stop band corresponding to the 3𝑄𝑥

resonance. Effectively, no beam survives when parked on top of the resonance line. There is a good

agreement between the stop bandwidth of the resonance from the simulations and the experimental

data. This observation hints that the model roughly captures the sextupole component around the

ring, even if imperfect. Nevertheless, left of the resonance, the experimental data deviates from the

model, given that the Recycler Ring is not optimized to run in this region.

On the other hand, Fig. 4.25 shows a similar comparison between static tune scans done with

SYNERGIA simulations and experimental results. The simulations were done by calculating the 𝐾2

coefficients that compensate 3𝑄𝑥 and then launching the simulations for 72800 turns. The intensity

and beam emittance data were recorded every 100 turns for these simulations. Experimentally, the

compensation currents were calculated from the response matrix method, and losses on top of 3𝑄𝑥

were optimized with a Nelder-Mead algorithm as mentioned in Sec. 4.4. The first thing to note is

that the beam survival ratio increases close to the resonance line with compensation—for both the

simulations and the experiment. Nevertheless, in this case, the simulation data does not agree well

with the experiments. In principle, the compensation done experimentally is better than that done

with the lattice model. Therefore, something is missing from the compensation algorithm applied

to the lattice model. The lattice model is imperfect and does not necessarily capture all the physics

in the real accelerator. More work must be done to match the simulations and experiments of static

tune scans with resonance line compensation.

95

Figure 4.25 Plot for static tune scan for the machine with 3𝑄𝑥 compensation, including a comparison
between experimental data and SYNERGIA simulations.

Additional experiments performing static tune scans were done. Figures 4.26 and 4.27 are

additional experiments performed at different dates than the ones presented in Figs. 4.24 and 4.25.

An accelerator is a machine that drifts in time; every day, it gets tuned and manipulated. From day

to day, different fields can be sampled from different elements. This is why the data can sometimes

change from experiment to experiment, especially if it was taken six months apart—such is the

case with these datasets.

Nevertheless, these additional experiments allow a better understanding of the accelerator and

improve the resonance compensation algorithms. Figures 4.26 and 4.27 show static tune scans

with beam size calculated from IPM data. Figure 4.26 shows the case for when the bare machine

was sampled—compensation sextupoles were off. As the working tune point gets closer to the 3𝑄𝑥

96

resonance line, Fig. 4.26 shows how the beam survival ratio goes down, and the beam size goes

up. The losses increase due to the beam size growth coming from the resonance. Eventually, it hits

the beam pipe, and losses appear. When operating on top of the resonance with no compensation,

almost no beam survives, and the IPM data can not be used.

On the other hand, Fig. 4.27 shows what happens with compensation turned on. Once the

compensation currents are set into the sextupoles, the beam survival ratio across the 3𝑄𝑥 resonance

line increases. Up to the point where approximately 95% of the beam survives while operating

on top of the resonance line—compare this to the uncompensated case where no beam survives.

Nevertheless, this region has some beam size growth, but not as large as the uncompensated case.

There is still sufficient beam that can be used to probe the beam size. With compensation at low

intensity, the beam size growth is suppressed such that 95% of the beam can now survive for 0.8

seconds on top of this third-order resonance line.

Figure 4.26 Static tune scan with beam survival ratio and IPM data box plots for the bare machine
at 2 Booster Turns of equivalent intensity.

It is worth specifying that the plots shown in Figs. 4.26 and 4.27 are box plots of the sigmas

of Gaussian-fitted IPM data. Box plots are used to represent distributions. The beam size data

is a distribution depending on the number of turns. Ultimately, the beam size can grow with the

97

number of turns, as shown in Fig 4.28 and in Sec. 6.4. Therefore, a simple scatter plot of the mean

with some error bars would be incorrect to illustrate IPM data—this type of plot can only be used

with centrally distributed data. This is the reason why box plots have been used in these plots.

Figure 4.27 Static tune scan with beam survival ratio and IPM data box plots for the machine with
3𝑄𝑥 compensation at 2 Booster Turns of equivalent intensity.

The 95% normalized emittance 𝜀𝑁,95% of the beam in the 𝑢 plane—either horizontal or vertical—

is the phase space area that encapsulates 95% of the particles in a beam traveling close to the speed

of light. The definition for this quantity is:

𝜀𝑛,𝑢,95% = 6𝛽𝐿𝛾𝐿𝜀𝑢,

(4.39)

where (𝛽𝐿, 𝛾𝐿) are the longitudinal relativistic factors. The geometrical RMS emittance 𝜀𝑢 was

calculated from:

𝜀𝑢 =

𝜎2
𝑢
𝛽𝑢

−

𝐷2
𝑢
𝛽𝑢

(cid:18) 𝛿 𝑝
𝑝0

(cid:19) 2

,

(4.40)

where 𝛽𝑥 = 27.10 meters, 𝐷𝑥 = 0.32 meters and 𝛿 𝑝/𝑝0 = 5𝑒 − 4 were the parameters estimated for

the horizontal IPM. The data from the horizontal IPM can be used to calculate the beam size 𝜎𝑥,

as shown in Figs. 4.26 and 4.27. All these values were extracted from the Ref. [40] information.

Figure 4.28 shows the 95% normalized horizontal emittance calculation from the IPM data

presented in Fig. 4.27. In this case, each point from each scatter plot is indexed by its decimated

98

turn number. The decimated turn number is the data point’s position in the IPM return array. Even

though the IPM can analyze up to 65000 turns, it can only output an array of 1024 decimated

turns. The data analyzed in these plots was for 65000 turns decimated into 1024 data points. As

mentioned, this data can not be represented with an error bar plot because it does not necessarily

follow a normal distribution. This can be seen in the emittance evolution against the turn number

illustrated in Fig. 4.28.

Figure 4.28 Static tune scan with beam survival ratio and emittance data calculated from IPM data
for the machine with 3𝑄𝑥 compensation at 2 Booster Turns of equivalent intensity.

Just as the beam size grows in Fig. 4.27 while on top of the resonance line, this translates into

normalized emittance growth shown in Fig. 4.28. Once the beam is injected close to the resonance

line at a normalized emittance of around 11 𝜋 mm mrad, it grows with the number of turns. In

particular, the Recycler Ring has an acceptance of 40 𝜋 mm mrad. This means that particles will

start getting lost close to this number. This can be seen for the values of around 𝑄𝑥 = 25.344 where

the emittance grows up to values larger than 40 𝜋 mm mrad, and, indeed, there is some loss in the

99

beam survival ratio plot (red dashed lines). Ultimately, this section has shown how using DCCT

and IPM data allows for the characterization of resonance lines and their compensation.

4.6 Additional Sextupoles for Resonance Compensation

As mentioned before in the discussion of Figs. 4.10, 4.11 and 4.23, the currents in the

compensation sextupoles that cancel out the ℎ3000 and ℎ1020 RDTs are too large. This section

explores the idea of introducing two additional normal compensation sextupoles to bring down the

𝐾2 coefficients that cancel out both of these RDTs, effectively bringing down the currents in all the

sextupoles. Therefore, this section looks at optimizing the location for these two new sextupoles.

The introduction of two new sextupoles implies that the RDT equation from Eq. 4.16 is modified,

and this reads:

−|ℎ3000| cos 𝜓3000

−|ℎ3000| sin 𝜓3000

−|ℎ1020| cos 𝜓1020

−|ℎ1020| sin 𝜓1020



























 (𝐵𝑎𝑟𝑒)

= 𝑴𝑨

𝐾 (𝑠𝑐220)

2

𝐾 (𝑠𝑐222)

2

𝐾 (𝑠𝑐319)

2

𝐾 (𝑠𝑐321)

2
𝐾 (1)
2
𝐾 (2)
2











































.

(4.41)

In this case, the response matrix gets appended to two new columns with the RDT sensitivity of

the two new sextupoles. This reads:

𝑴𝑨 =

(cid:16)

(cid:17)

𝑴















cos 3𝜙𝑥(1)

1

48 𝑙𝑠𝑐 𝛽3/2
𝑥(1)
48 𝑙𝑠𝑐 𝛽3/2
𝑥(1)

1

cos 3𝜙𝑥(2)

1

48 𝑙𝑠𝑐 𝛽3/2
𝑥(2)
48 𝑙𝑠𝑐 𝛽3/2
𝑥(2)

1

sin 3𝜙𝑥(1)
𝛽𝑦(1) cos (cid:0)𝜙𝑥(1) + 2𝜙𝑦(1)
𝛽𝑦(1) sin (cid:0)𝜙𝑥(1) + 2𝜙𝑦(1)

1

16 𝑙𝑠𝑐 𝛽1/2
𝑥(1)
16 𝑙𝑠𝑐 𝛽1/2
𝑥(1)

1

(cid:1)

(cid:1)

1

16 𝑙𝑠𝑐 𝛽1/2
𝑥(2)
16 𝑙𝑠𝑐 𝛽1/2
𝑥(2)

1

sin 3𝜙𝑥(2)
𝛽𝑦(2) cos (cid:0)𝜙𝑥(2) + 2𝜙𝑦(2)
𝛽𝑦(2) sin (cid:0)𝜙𝑥(2) + 2𝜙𝑦(2)

,














(4.42)

(cid:1)

(cid:1)

where the beta functions and phases are evaluated at a lattice location that has still yet to be

determined. The original matrix 𝑴 is defined from Eq. 4.27.

100

The first step was to minimize the function

𝑓 (cid:0)𝛽𝑥(𝑖), 𝛽𝑦(𝑖), 𝜙𝑥(𝑖), 𝜙𝑦(𝑖)

(cid:1) = max

𝐾 (𝑠𝑐220)

2

𝐾 (𝑠𝑐222)

2

𝐾 (𝑠𝑐319)

2

𝐾 (𝑠𝑐321)

2
𝐾 (1)
2
𝐾 (2)
2






















(cid:169)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:173)
(cid:171)

,






















(cid:170)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:174)
(cid:172)

(4.43)

where the max operator takes the absolute maximum of the vector. Therefore, it would correspond

to the maximum 𝐾2 coefficient in the vector. The minimization procedure was done using a

Nelder-Mead algorithm. Ultimately, the result would be a set of 𝛽𝑥(𝑖), 𝛽𝑦(𝑖), 𝜓𝑥(𝑖), 𝜓𝑦(𝑖), where 𝑖

runs from 1 to 2 (representing the two additional sextupoles), and this specifies the location where

the maximum 𝐾2 coefficient from all the sextupoles gets minimized. An instance of Nelder-Mead

was launched 10000 times; each time, the starting point would be chosen randomly to ensure most

of the parameter space was explored. All the results could be plotted in histograms to create

constraints for locations depending on the beta and phase advance. Figure 4.29 shows an example

of such histograms. This plot shows the relations and constraints the new sextupoles must comply

with to minimize the maximum 𝐾2 coefficient. For example, Fig. 4.29 shows that the horizontal

beta functions of both sextupoles should have a ratio of around 1. Additionally, it shows that the

horizontal phase advances should lie close to thirds of 𝜋, e.g., 𝜋/3, 2𝜋/3, 𝜋, etc. . . .

101

Figure 4.29 2D Histogram of optimum lattice locations where the compensating 𝐾2 coefficients
with the two new sextupoles are minimized. Red dashed lines correspond to thirds of 𝜋.

With the previous constraints already set up, the next step is to look into all the possible pairs of

drift spaces in the lattice where two sextupoles would fit. Therefore, all the possible combinations

of two available drift spaces longer than 1 meter and horizontal dispersion less than 0.1 meters were

explored. This last constraint came from the fact that these new sextupoles should not introduce any

additional chromaticity into the Recycler. In total, the number of total possible combinations was

39340. From all of these possibilities, the first filter analyzed only the combinations that complied

with the constraints imposed by minimizing Eq. 4.43. With this first filter, the number of pair

candidates was narrowed to 776. In these 776 pairs of locations, the sextupoles were introduced,

and the resulting maximum 𝐾2 coefficient was calculated. Figure 4.30 shows these candidates

with their corresponding metric for this particular constraint superimposed onto the 2D histogram

of Fig. 4.29. The color bar of Fig. 4.30 labels the candidates by their particular maximum 𝐾2

coefficient, as defined in Eq. 4.43.

102

Figure 4.30 2D Histogram of optimum lattice locations where the compensating 𝐾2 coefficients
with the two new sextupoles are minimized with available drift spaces candidates scattered on top
of it.

With the number of candidates narrowed down, the next step is to select the idle pair of locations

to introduce these new sextupoles. Figure 4.31 was used as a tool for this purpose. The left vertical

axis shows the dispersion curve (golden solid continuous line) of the Recycler Ring to ensure that

the new candidates are indeed in areas of low dispersion. The red crosses and the blue exes represent

the pair of drift space candidates where the sextupoles could go. The right vertical axis specifies

the maximum 𝐾2 coefficients when the sextupoles are introduced in these locations. The minimum

of this quantity is given by the pair of sextupoles introduced in the region 620. By observing Fig.

4.31 and making sure the locations in 620 were viable in terms of space and logistics, the new

locations for two new compensation sextupoles were chosen.

103

Figure 4.31 Dispersion function of the Recycler Ring with possible new locations to introduce a
pair of sextupole magnets that cancel out ℎ3000 and ℎ1020 simultaneously. The right y-axis shows
the maximum compensation 𝐾2 coefficient needed with these new candidates.

Once the new locations are pinned down, the compensation of ℎ3000 and ℎ1020 can be verified.

Table 4.6 shows the 𝐾2 coefficients that need to go into the sextupoles for 3𝑄𝑥 = 76 and 𝑄𝑥 + 2𝑄 𝑦 =

74 compensation. These values can be compared to those in Table 4.5, which are one order of

magnitude less. Furthermore, all the values in Table 4.6 are less than 1.0 m−3, the equivalent

current limit for the 𝐾2 coefficients. Ultimately, these two new sextupoles, SC620a and SC620b,

have brought down the compensating sextupole coefficients by one order of magnitude. Their

installation would enable a viable option to compensate both resonance lines simultaneously.

104

Table 4.6 𝐾2 coefficients that cancel out the ℎ3000 and ℎ1020 term when the two new 620 sextupoles
are introduced. These values were calculated for a Recycler lattice tuned to 𝑄𝑥 = 25.4422 and
𝑄 𝑦 = 24.3915.

Sextupole 𝐾2 [m−3]

SC220

-0.263835

SC222

0.463129

SC319

0.480173

SC321

-0.407482

SC620a

0.581187

SC620b

0.511139

As it was previously done in Sec. 4.1, the ℎ3000 and ℎ1020 RDTs can be calculated after

introducing the values of 𝐾2 coefficients specified in Table 4.6. Figures 4.32 and 4.33 show the

RDT calculations with these settings. Both plots show how the global RDTs go down to 0 at the end

of the sum along the lattice. Nevertheless, more noteworthy is the fact that the RDT calculations

do not show huge variations in the cumulative sum, as they did in Figs. 4.12 and 4.13. This means

that these two additional sextupoles will delocalize the amount of sextupole that must be introduced

into the ring to cancel out both RDTs. Therefore, the RDT variation around the ring will be smaller

with these new settings than with the four original compensation sextupoles.

As of Spring 2024, these new 620 sextupoles have been installed in the Recycler Ring. Nev-

ertheless, they have not been commissioned nor tested for operations. The next step would be to

verify that both lines can be compensated with these six normal compensation sextupoles. The

final performance of these new sextupoles will be published in a future publication.

105

Figure 4.32 Distribution of the ℎ3000 term around the ring with individual contributions from each
relevant element and the cumulative sum from an arbitrary starting point. This calculation includes
the new 620 sextupoles powered at the correct currents to simultaneously cancel out ℎ3000 and
ℎ1020.

106

Figure 4.33 Distribution of the ℎ1020 term around the ring with individual contributions from each
relevant element and the cumulative sum from an arbitrary starting point. This calculation includes
the new 620 sextupoles powered at the correct currents to simultaneously cancel out ℎ3000 and
ℎ1020.

107

CHAPTER 5

RESONANCE COMPENSATION STUDIES AT THE CERN
PROTON SYNCHROTRON BOOSTER

5.1 General Description

The Proton Synchrotron Booster (PSB) is the first circular accelerator in the CERN accelerator

complex that ultimately leads to the LHC. Figure 5.1 shows the entire chain of accelerators at

CERN, feeding a variety of physics experiments [48]. Following the successful implementation

of the LHC Injectors Upgrade (LIU) [49], the PSB receives 𝐻− ion beam from the Linac4 at an

energy of 160 MeV. Interestingly, the PSB is not just one ring but four identical synchrotron rings

stacked on each other. This design counteracts the space charge effects, which are the largest

in low-energy machines. Once the ion beam enters the PSB rings, the electrons are stripped off

through a charge-exchange process with a carbon foil, and a proton beam is achieved [50]. The

proton beam is then accelerated from an energy of 160 MeV to 2 GeV. Each ring extracts one

bunch and is injected in different buckets into the Proton Synchrotron (PS). This description is true

for LHC-type beams. Nevertheless, the PSB can also feed protons to other customers such as its

highest-intensity user—ISOLDE (Isotope mass Separator On-Line facility) [51].

One ring of the PS Booster has a total circumference of 157.08 meters. Multiplying this

quantity by the four rings, one gets a length of 628.32 meters, the exact circumference of the next

accelerator, the Proton Synchrotron (PS). The PSB has a superperiodicity of 16, meaning it has 16

identical fundamental cells. Each cell has a length of 9.82 meters, housing a sequence of bending

magnet, focusing quadrupole, defocusing quadrupole, focusing quadrupole, and bending magnet

[52, 53]. Between these main components are drift spaces used for RF insertions, diagnostic

devices, injection and extraction devices, and additional multipole corrector magnets. At injection

energy, 160 MeV, protons have a revolution period of 1.01 𝜇𝑠, while at the extraction energy of 2.0

GeV, they have a revolution period of 0.553 𝜇𝑠. The beam stays in the PS Booster for around 530

microseconds, where CERN and its accelerators are synchronized with a 1.2-second basic period.

108

Figure 5.1 A graphic overview of all accelerators in operation at CERN as of 2022. Original
image taken from Ref. [48]. This file is licensed under the Creative Commons Attribution 4.0
International license.

5.2 Tune Diagram and Operation

Figure 5.2 illustrates the tune diagram dynamics that the LHC-type beam undergoes at the PS

Booster [51, 54, 55]. As mentioned before, the beam gets injected at an energy of 160 MeV. At this

low energy, the tune footprint is large enough that the spread can reach up to 0.5, i.e., Δ𝑄𝑢 ≈ −0.5.

The nominal injection tunes are around 𝑄𝑥 = 4.40 and 𝑄 𝑦 = 4.45, in order to accommodate the

footprint between the integer resonance lines 𝑄𝑢 = 4.0 and the half-integer line 2𝑄 𝑦 = 9. As the

beam is accelerated, the quadrupoles are ramped up to match the increasing beam rigidity, but,

additionally, a tune ramp is introduced to move the shrinking footprint to a less resonance-populated

area in the tune diagram. The nominal extraction tunes are around 𝑄𝑥 = 4.17 and 𝑄 𝑦 = 4.23. At

extraction, the beam tune footprint has shrunk by a factor of (𝛾3

𝐿), as explained by the beam
perveance definition in Eq. 2.66. At extraction, the footprint is smaller than 0.05, i.e., |Δ𝑄𝑢 | ≲ 0.05.

𝐿 𝛽2

109

Figure 5.2 Operational tune footprint for PSB beam at injection (cool color map) and footprint after
the beam has been accelerated to 2 GeV (warm color map). During acceleration, there is a tune
ramp illustrated with the fuchsia arrow.

As illustrated in Fig. 5.2 with the fuchsia arrow, the nominal tune ramp crosses several resonance

lines. Uncorrected, these resonance lines will lead to beam loss during the tune ramping. These

include four third-order resonance lines and four fourth-order lines. It is worth reminding the reader

that the third-order lines are excited by sextupole-like components in the ring, while fourth-order

lines are excited by octupole-like fields. The third order resonances include two normal sextupole

lines, 3𝑄𝑥 = 13 and 𝑄𝑥 + 2𝑄 𝑦 = 13, and two skew sextupole lines, 3𝑄 𝑦 = 13 and 2𝑄𝑥 + 𝑄 𝑦 = 13.

For the octupole case, these include two normal octupole lines, 4𝑄𝑥 = 17 and 𝑄𝑥 + 3𝑄 𝑦 = 17,

two skew octupole lines, 4𝑄 𝑦 = 17 and 3𝑄𝑥 + 𝑄 𝑦 = 17, and the octupole coupling sum resonance,

2𝑄𝑥 +2𝑄 𝑦 = 17. Figure 5.3 shows all these resonance lines summarized in one tune diagram. All of

these resonance lines have different strengths in each Booster ring. It is worth pointing out that the

coupling resonance 𝑄𝑥 − 𝑄 𝑦 = 0 is not strongly excited in the low-intensity operation of the PSB.

At higher intensities, it has been demonstrated that this resonance line is driven in fourth-order by

110

space charge [53].

Figure 5.3 Portion of the tune diagram enclosing the operational tunes of the PS Booster with
relevant resonance lines labeled.

Like the plots shown in Fig. 4.21 and Fig. 4.22, loss maps can be used to visualize the strength

of the resonance lines in the CERN PS Booster. Figure 5.4 shows loss maps for each of the four

rings in the PSB. These plots are created slightly differently from the ones in the Recycler Ring.

The plots shown in Fig. 5.4 are an average from four loss maps. One where losses are mapped from

(a) left to right, i.e., fixing 𝑄 𝑦 and ramping from 𝑄𝑥 ≈ 4.49 to 𝑄𝑥 ≈ 4.15, (b) another one from

111

right to left, (c) one from top to bottom, i.e., fixing 𝑄𝑥 and ramping from 𝑄 𝑦 ≈ 4.49 to 𝑄 𝑦 ≈ 4.15,

(d) and one from bottom to top. Therefore, these plots show an average from mapping the losses in

four directions.

Figure 5.4 Dynamic loss maps for the bare machine of the four rings (R1, R2, R3, and R4) in the
PS Booster. The plots are an average of scanning in 4 directions. Plot provided by F. Asvesta.

Figure 5.4 shows how resonances are excited differently for each ring at the Proton-Synchrotron

Booster. For example, third-order lines are more excited for Ring 2 than Ring 1. This feature hints

that every ring will have different values for the corrector magnets used for resonance compensation.

The fourth-order order skew resonances are not observed in this loss map. The work in the following

sections looks to calculate the currents of the corrector magnets used for resonance compensation.

This approach involves using two optimization procedures to find the values that clear out the losses

from the plots in Fig. 5.4.

112

5.3 Optimization Algorithms for Resonance Compensation

The whole objective of the following work is to minimize the losses during the PS Booster’s

operational cycle, as explained in Fig. 5.2. The losses during the cycle come from particles

falling on top of third-order and fourth-order resonance lines, as identified in Fig. 5.3. In order

to decrease the strength of these resonances, sextupoles and octupoles can be used to control their

Resonance Driving Terms (RDTs). Nevertheless, measuring and minimizing the RDTs may not

always be ideal. For this work, the main observable to optimize was the beam loss from crossing

the resonance lines, which is related to the amplitude of the RDTs. Figure 5.5 shows an illustration

of the experimental setup used in order to measure this beam loss.

Figure 5.5 Experimental setup of the tune diagram dynamics for optimizing resonance compensation
used in the PS Booster.

The experimental setup introduced in the PS Booster involved several steps to find the optimal

113

compensation currents. First, a low-brightness beam was injected into every ring at an energy of

160 MeV. Energy is not ramped up for this configuration, and the machine stays at a flat bottom.

Second, a tune ramp was programmed into the quadrupoles in the ring in order to go from initial

tunes of 𝑄𝑥 = 4.40 and 𝑄 𝑦 = 4.45 to a final setting of 𝑄𝑥 = 4.17 and 𝑄 𝑦 = 4.23. This particular

setting was introduced to mimic the operational tune ramp of the LHC-type beam. The start of this

tune ramp occurs within 𝑡0 = 300 ms of the start of the cycle and ends at 𝑡 𝑓 = 600 ms, i.e., all of

this occurs within a 300 ms time window. During this time window, the beam loss is measured by

comparing the beam current at the end of the window to the initial value of the beam current. The

currents fed to the corrector magnets remain constant during this measurement. Figures 5.5 and

5.6 summarize this experimental setup.

While monitoring the beam loss, the corrector magnets used for compensation are varied every

cycle according to the optimization algorithm. Table 5.1 summarizes the 11 elements used for

resonance compensation for this work. Out of these elements, there are 6 normal sextupoles, 3

skew sextupoles, and 2 normal octupoles. Figure 5.6 shows an example of the power cycle in these

magnets for one optimization step. Before the tune ramp, most actors show currents at or near zero.

As preparation for the tune ramp, the magnets are powered to the set values—per the optimizer

calculation. During the tune ramp, they are set to a constant value and powered off once the cycle is

finished. These magnets were varied for each ring, and each ring had its independent optimization

run. Theoretically, in order to fully correct eight resonance lines, one needs at least 16 correctors.

Nevertheless, this work aimed to find a solution to this over-constrained problem through advanced

optimization algorithms.

The two optimization algorithms used were Bayesian Optimization and BOBYQA (Bound

Optimization BY Quadratic Approximation). In order to implement these algorithms, the special

application GeOFF (Generic Optimization Frontend and Framework) was used [56]. This graph-

ical application is designed to facilitate numerical optimization through various algorithms and

reinforcement learning on CERN accelerators. It incorporates programmable interfaces that can

be used to specify the hyperparameters of the optimization algorithms.

114

Table 5.1 List of elements (optimization actors) in the PS Booster at CERN used for resonance
compensation optimization as present in all four PS Booster rings.

Actor
1
2
3
4
5
6
7
8
9
10
11

Type
Name
Normal Sextupole
XN04L1
Normal Sextupole
XN06L1
XN09L1
Normal Sextupole
XN012L1 Normal Sextupole
XN0311L1 Normal Sextupole
XN0816L1 Normal Sextupole
ON0311L1 Normal Octupole
ON0816L1 Normal Octupole
Skew Sextupole
XSK2L4
Skew Sextupole
XSK4L1
Skew Sextupole
XSK6L4

Figure 5.6 Waveform for the currents fed to the correctors (actors) as a function of time in the study
cycle for one particular iteration of an arbitrary optimization instance.

115

5.3.1 Bayesian Optimization

Bayesian optimization is a well-established technique used in machine learning and optimization

tasks to efficiently search for the optimal solution within a parameter space, particularly when the

objective function is expensive or time-consuming to evaluate—this is the case for an accelerator.

It combines probabilistic modeling with the principles of Bayesian inference to iteratively update a

model of the objective function based on observed data, gradually refining its understanding of the

parameter space [47]. By balancing exploration (searching for promising regions) and exploitation

(leveraging known information to identify the best areas), Bayesian optimization aims to find the

global optimum while minimizing the number of evaluations required. This method is particularly

useful in hyperparameter tuning for complex models, where traditional grid search or random

search approaches may be impractical due to computational costs.

At the heart of Bayesian optimization are Gaussian processes. A Gaussian process (GP) is

a statistical model representing a collection of random variables, any finite subset with a joint

Gaussian distribution [47]. GPs are defined by a mean function and a covariance function, which

capture prior beliefs about the underlying function being modeled and the correlations between

different points in the input space. These special types of statistical models are key to estimating

the underlying function that wants to be optimized. Through Bayesian inference, GPs can be

updated with observed data, yielding posterior distributions that can inform predictions and lead to

an optimized sampling for the Bayesian optimization algorithm.

Figure 5.7 shows an example of the evolution of the objective function during a Bayesian

optimization procedure. In this case, the objective function is the normalized beam loss after the

tune ramp illustrated by Fig. 5.5. It can be seen how the Bayesian optimizer finds solutions that

effectively cancel out the beam loss from crossing the resonances. Nevertheless, given that this

optimizer is built to find a global minimum, it will keep sampling other regions to ensure the best

solution is not a local minimum. The color map of Fig. 5.7 shows how for some of these cases,

the optimizer prioritizes exploration and drifts to some unknown region where the losses are high.

For these cases, the underlying Gaussian process will learn that there is no worth in exploring these

116

regions. Ultimately, the configuration with the least relative beam loss—the best configuration—is

saved and kept as the optimum solution.

Figure 5.7 Normalized beam current plots for the Bayesian Optimization method done at Ring 1.

The bottom plot of Fig. 5.8 shows the explicit steps of each actor versus the number of

iterations. Additionally, the top plot shows the trend of the objective function (relative beam loss)

as the number of iterations increases. There are significant oscillations in the relative beam loss

values, especially at the beginning, but a general trend towards minimization as the algorithm

progresses through iterations. The early iterations reflect the exploration phase. BO is sampling

points that give a broad understanding of the objective function’s landscape. As iterations progress,

there is a trend toward certain regions in the parameter space. This trend indicates a shift from

exploration to exploitation, where the algorithm samples more from areas it believes to be near the

optimum. The narrowing of actor current variability suggests a reduction in uncertainty about the

location of the minimum beam loss. The GP model is becoming more informed and better trained.

At the end of the optimization instance, such as the one shown in Fig. 5.8, the configuration that

117

gave the smallest relative beam loss is saved. GeOFF sets the default configuration of the correctors

to these best values.

The corrector values for the best configurations are important to discuss from Fig. 5.8. It can

be seen from Fig. 5.8 that some correctors land on the limit values, e.g., the limits for the normal

sextupoles (XNO magnets) were [-50,50]. This behavior is especially apparent for the octupole

correctors (ONO magnets), which have limits from -80 to 80, e.g., the magnet ONO816L1 is maxed

out. Nevertheless, this was expected given that only two octupole correctors were used to correct

four fourth-order lines, whereas one would need eight correctors to cancel out all the fourth-order

RDTs fully. It is important to remind the reader that to implement these types of algorithms with

several actors efficiently, the currents need to be normalized between [0,1] to ensure that all the

data is on the same scale and improve the model performance. This procedure is done by GeOFF

internally.

Figure 5.8 Summary for Bayesian optimization of resonance compensation applied to Ring 2 in the
CERN PSB.

While Fig. 5.8 shows a relatively smooth Bayesian optimization instance for Ring 2, Figs. 5.9

and 5.10 show a different story for Ring 3. This particular optimization instance failed to find a

configuration yielding a relative beam loss of less than 5% after 200 iterations—unacceptable for

operation. This feature can be seen on the top plot of Fig. 5.9, which follows the objective function

118

across iterations.

In order to find an acceptable configuration, the BO instance was performed

again, but in this case, the initial Gaussian process was trained with the data from the previous

incomplete run. The additional 350 iterations from this second instance are shown in Fig. 5.10.

This second instance yielded configurations that led to acceptable values for the relative beam loss.

Possible explanations for this hurdle include inadequate exploration of the parameter space, poor

choice of the kernel for the Gaussian Process, poor choice of hyperparameters for this BO, and

noise in the measurements.

Figure 5.9 Summary for the first instance of Bayesian optimization towards resonance compensation
applied to Ring 3 in the CERN PSB.

119

Figure 5.10 Summary for the second instance of Bayesian optimization of resonance compensation
applied to Ring 3 in the CERN PSB. This second optimization instance was trained with the data
from the first instance.

5.3.2 BOBYQA (Bound Optimization BY Quadratic Approximation)

BOBYQA, which stands for Bound Optimization BY Quadratic Approximation, is an opti-

mization algorithm commonly used for solving constrained optimization problems [57]. Unlike

gradient-based methods, BOBYQA belongs to the class of derivative-free optimization algorithms,

making it suitable for scenarios where the objective function is not differentiable or computationally

expensive to evaluate. The algorithm iteratively builds a quadratic approximation of the objective

function within a trust region, a bounded area around the current solution. By iteratively updating

the quadratic model and moving towards the predicted optimum within the trust region, BOBYQA

efficiently explores the search space while minimizing the number of function evaluations.

It

employs a bound constraint strategy to ensure that the search remains within specified bounds. In

this case, BOBYQA was used to solve the same optimization task from the last section but using a

special type of beam at the PSB, the BCMS (Bunch Compression Merging and Splitting) beam.

The BCMS configuration is special given that it is produced with a smaller intensity (compared

to the LHC-type beam from the previous section), and, hence, it has a smaller tune spread. For

this case, the injection tunes were changed from those shown in Fig. 5.5. The injection tunes

120

were chosen to avoid one normal sextupole and one skew-sextupole resonance. This different

configuration changes the optimum compensation currents compared to the Bayesian optimization

results.

Similar to the exercise done in Fig. 5.7, Fig. 5.11 shows the relative beam loss profiles for an

instance of the BOBYQA optimization done on Ring 1. As the iterations progress (transitioning

from blue to yellow), the spread in beam current values decreases. This decrease suggests that

the algorithm is refining its search based on the feedback from the objective function and moving

towards regions of the parameter space that yield better optimization results. The later iterations,

indicated by yellow lines, show beam current trajectories closer together. This convergence pattern

signifies that the algorithm has likely identified a promising region in the parameter space and

is now fine-tuning the parameters.

It minimizes the variance between iterations, suggesting an

approach toward a stable solution.

Figure 5.11 Normalized beam current plots for the BOBYQA method done at Ring 1.

Similar to the exercise done for the Bayesian optimization, Fig. 5.12 shows the evolution of

121

the objective function and actor currents for a BOBYQA optimization instance on a BCMS beam.

The initial values were set from the Bayesian optimization results for each ring, as discussed in

the last section. The trend of the objective function (top plot) shows a general decrease in relative

beam loss as the number of iterations increases, indicating that the optimization successfully finds

parameters that result in less beam loss. The bottom graph illustrates the variation of the actor

currents in arbitrary units across iterations. A unique line pattern and color distinguish each actor.

Some correctors are kept at 0, given that they were unavailable for Ring 3. Implementing BOBYQA

through GeOFF ensures the algorithm does not fall to a local minimum by bumping the settings to

a new configuration once it falls to a stable one. That can be seen from the sudden oscillations in

the objective function plot from Fig. 5.12.

Figure 5.12 Summary for BOBYQA optimization of resonance compensation applied to Ring 3 in
the CERN PSB for BCMS (Bunch Compression Merging and Splitting) beam.

Comparing this BOBYQA instance with the Bayesian optimization results, one can see that the

exploration phase is much shorter and less broad for the BOBYQA instance. Similarly, the variation

steps from BOBYQA in the actor currents are smaller than BO. This comparison is true because of

the underlying mechanism used to estimate the objective function and its uncertainty. BOBYQA

uses deterministic local quadratic models, while Bayesian Optimization uses probabilistic models

that measure uncertainty. BOBYQA is generally considered a local search technique, meaning it

122

may be more prone to finding a local minimum than a global minimum. Hence, the initial point for

BOBYQA holds more importance. Bayesian optimization is generally used for global optimiza-

tion and is particularly strong in high-dimensional spaces or when function evaluations are very

expensive. In summary, while both algorithms are powerful tools for optimization in accelerator

physics problems, the choice between BOBYQA and Bayesian optimization will depend on the

specific characteristics of the problem at hand, including the nature of the objective function, the

presence of constraints, the dimensionality of the problem, practical considerations like the avail-

ability of computational and experimental resources, and, ultimately, availability of study/machine

development time.

5.4 Experimental Verification of Compensation

The whole point of the optimization algorithms explained in the last sections was to reduce the

losses that show up in the loss map from Fig. 5.4. Figure 5.13 shows a new loss map with the best

configurations found for each ring using the Bayesian optimization procedure on LHC-type beam.

When comparing both loss maps, it is clear that with these new configurations, the loss maps have

been cleared out of losses in the region of interest. The immediate losses are decreased by nearly one

order of magnitude in the region occupied by the tune ramp. In particular, these new configurations

largely suppress the third-order resonances that dominated the losses in Fig 5.4. Nevertheless, some

resonance lines are still visible in the loss maps in Fig. 5.13, e.g., 𝑄𝑥 − 2𝑄 𝑦 = −4. Given that these

lines are not in the region of PSB operation, they are not of particular concern.

This chapter has shown a different approach to resonance compensation from the one presented

in Ch. 4. Chapter 4 showed a physics-informed approach by minimizing the Recycler’s RDTs

and, ultimately, leading to the reduction of beam loss. On the other hand, the previous sections

of this chapter showed an optimization-based approach, which can be considered a brute-force

line of action. In this approach, all the actors are thrown into an optimization algorithm, which

finds a numerical solution that minimizes the objective function. There is no physics involved, just

numerical optimization. In particular, BOBYQA’s deterministic approach can quickly determine a

solution when the accelerator operates near a resonance compensation optimum. With its proba-

123

bilistic nature, Bayesian Optimization is better suited for exploring unknown or poorly understood

operational regimes. All the previous approaches have been proven to work and yield satisfactory

results in the context of resonance compensation. The choice between optimization-based and

physics-informed approaches—–or a combination thereof—–depends on the specific context of the

particle accelerator’s operational goals, the available data, and the accuracy of the model.

Figure 5.13 Dynamic loss maps for the four rings in the PS Booster with the best configuration from
the Bayesian optimization of the resonance compensation. The plots are an average of scanning in
4 directions. Plot provided by F. Asvesta.

124

CHAPTER 6

HIGH INTENSITY STUDIES

6.1 Space Charge Tune Shift

For this chapter, the focus goes back to the Fermilab Recycler Ring. All the experiments and

measurements done in Ch. 4 were done at low intensities, i.e., less than 1e10 particles per bunch

(ppb). Nevertheless, current operations and future operations under PIP-II objectives are done at

higher intensities, i.e., 5e10 ppb for current operations and 8e10 ppb for PIP-II. Therefore, it is

relevant to explore resonance compensation at higher intensities.

An important parameter for high-intensity operation is the space charge tune shift. As mentioned

in Ch. 2, the space charge tune shift is incoherent, meaning every particle will feel a different

magnitude of this effect. Nevertheless, as shown in Figs. 2.5 and 3.8, this incoherent quantity

will be contained within a necktie distribution. This necktie profile will define the space charge

tune spread. When particles are in the beam core, they will feel the largest space charge potential,

leading to the largest detuning and defining the edge of the necktie distribution.

Equation 2.78 showed a general way to calculate the space charge tune shift for particles at

different amplitudes 𝐽𝑢. Nevertheless, this calculation is an involved process where the envelope

equation has to be solved around the ring simultaneously as the Poisson equation, Eq. 2.65.

Different approximations can be made in order to make rapid estimates of the space charge tune

shift. Furthermore, a more important quantity is the maximum tune shift of the beam core

particles. This tune shift will ultimately define the space charge tune spread. While this quantity

can be calculated using PySCRDT [26], a cruder estimate can also be calculated. This estimate

involves using a smooth-lattice approximation and circular beam pipe approximation, as used in

Ref. [58].

With these simplifications, the space charge tune spread can be found as follows:

Δ𝑄𝑢,𝑠𝑐 =

−3𝑁𝑏𝑟0𝑅𝑆

4𝜎𝑧 𝛽𝐿𝛾2

𝐿𝜀𝑛,𝑢,95%

,

(6.1)

where 𝑁𝑏 is the number of protons per bunch, 𝑟0 = 1.535×10−18 meters is the classical radius of the

125

proton, 𝜎𝑧 = 0.5726 meters is the bunch length, 𝑅 = 𝐶/(2𝜋) is the radius of the Recycler Ring with

a circumference of 𝐶 = 3319.4 meters, 𝑆 = 1.596 is a geometrical factor of the bunch, 𝜀𝑛,𝑢,95% is

the 95% normalized emittance in the 𝑢 plane, and (𝛽𝐿, 𝛾𝐿) are the longitudinal relativistic factors.

The following section explores a method to measure this tune spread.

6.2 Measurement of Tune Spread

The 3𝑄𝑥 = 76 line can be used to estimate the tune spread of the beam in the Recycler Ring.

As seen before, when a dynamic tune ramp is set, and 3𝑄𝑥 is crossed, 95% of the beam is lost.

By monitoring the tune ramp set into the tune trombone and assuming a linear tune ramp, one can

correlate the tune instance when losses first appear with the tune spread. The tune distance between

the losses first appearing and the 3𝑄𝑥 line is a crude estimate for the tune spread. This estimate is

crude because particles in the beam core—particles with the largest tune shift—will start feeling

the resonance and drift to higher amplitudes. The migration of particles from the beam core to the

tails means that the tune shift will be smaller. There is a time delay (tune distance delay) between

particles touching the resonance and particles being lost to the aperture. While this delay exists

given the large strength of the 3𝑄𝑥 line, just hitting the resonance head-on can still give a crude

estimate for the value of the tune spread and its behavior.

Figure 6.1 shows slices of loss maps at different intensities used to measure this tune spread.

The top slice corresponds to an intensity of 2 Booster Turns (BTs), i.e., approximately 0.6 × 1010

particles per bunch for this experiment. The bottom slice corresponds to an intensity of 17 BTs,

approximately 6.0×1010 ppb, close to the maximum intensity Booster can provide in a single batch.

The important thing to note here is that as the intensity increases, the losses happen earlier in the

tune ramp. This behavior can be correlated to the space charge tune spread using the location of

the 3𝑄𝑥 line, which corresponds to where all the beam is lost. For example, for the bottom slice,

the losses start at 𝑄𝑥 = 25.36, while all beam is lost at 𝑄𝑥 = 25.34. This feature means the space

charge tune spread can be approximated to 0.02. Ultimately, Fig. 6.1 shows beams with space

charge tune spreads that range from 0.005 to 0.02.

126

Figure 6.1 Dynamic loss strips at a fixed vertical set tune of 𝑄 𝑦 = 24.45 and on a range from
𝑄𝑥 = 25.39 to 𝑄𝑥 = 25.32 for different intensities.

127

Another thing to note from Fig. 6.1 is how the space charge tune spread tends to stabilize to a

value of around 0.02 for intensities higher than 12 BTs or 4.4 × 1010. This behavior is because the

equilibrium emittance grows exponentially at higher intensities. Specifically, Ref. [40] shows how

the beam emittance grows with intensity in the Recycler Ring. It looks at the first turn emittance

coming from the Booster Ring, which will dictate the initial tune shift of the Recycler Ring. This

study, combined with calculations of the space charge tune shift provided in Ref. [58], explains this

behavior. The emittance grows exponentially as the beam intensity increases, leading to a saturated

tune shift value. Figure 6.2 summarizes this behavior by plotting the measured tune spread from

losses, the calculated tune spread from the first-turn abort emittance, and, finally, the first-turn abort

95% emittance on the right vertical axis.

The maximum tune spread shown in Figs. 6.1 and 6.2 is around 0.02, which only compares to

1/5 of the PIP-II projection. It is worth pointing out that these experiments are done for a single

batch with no slip stacking. Therefore, with slip stacking and smaller beam emittance coming from

Booster, space charge tune shifts in the RR of around 0.1 will be possible in the PIP-II era. For

now, Fig. 6.2 shows how the tune spreads are not large enough to consider space-charge-dominated

losses in the Recycler Ring relevant. All current losses are emittance-dominated, meaning the

beam from Booster is hitting on the acceptance of the Recycler.

128

Figure 6.2 Tune spread measurements compared to tune spread calculations from horizontal emit-
tances measured with multiwire data on first-turn abort beam.

6.3

Intensity-Dependent Effects and Non-Gaussian Beam Profiles

As mentioned in Sec. 3.5, the Ion Profile Monitor system is useful for extracting information

about the transverse beam distribution. In particular, fitting a Gaussian distribution and extracting

the scaling parameter or sigma 𝜎𝑢 can help characterize the transverse beam size. Such as it was

done in Figs. 4.26 and 4.27, the beam size growth can be correlated to the beam loss in the machine.

Nevertheless, one can go further and look into deviations from Gaussian distributions. Reference

[59] looks into an application and explanation for these non-Gaussian profiles.

In particular, looking at how the beam tails populate with incoming particles from the beam core

is interesting. This population happens when the space charge detuning is large enough for particles

in the beam core to start touching resonance lines and, therefore, migrate to larger amplitudes, i.e.,

129

the beam tails. In order to quantify this effect, one can use q-Gaussian distributions [59, 60]. The

definition for the probability density function 𝑓𝑞 (𝑥) of this distribution for a given 𝑞 reads:

𝑓𝑞 (𝑥) =

1

𝜎𝐶𝑞

𝑒𝑞

√
2

(cid:18)

−

𝑥2
2𝜎2

(cid:19)

,

(6.2)

where 𝜎 is the scale parameter for the distribution and corresponds to the usual standard deviation

of a Gaussian 𝜎 when 𝑞 = 1. Additionally, 𝐶𝑞 is the normalization constant. The auxiliary function

𝑒𝑞 (𝑥) is defined as:

𝑒𝑞 (𝑥) =






exp (𝑥),

if 𝑞 = 1

(1 + (1 − 𝑞) 𝑥)

1
1−𝑞 ,

if 𝑞 ≠ 1 and (1 + (1 − 𝑞) 𝑥) > 0

,

(6.3)

0

if 𝑞 ≠ 1 and (1 + (1 − 𝑞) 𝑥) ≤ 0

where the usual Gaussian distribution is recovered when 𝑞 = 1. The distribution is said to have

heavy tails when 𝑞 > 1 and light tails when 𝑞 < 1. The usual Gaussian distribution is recovered

when 𝑞 = 1, and the usual Cauchy-Lorentzian distribution is recovered when 𝑞 = 2. Figure 6.3

compares the probability density function for several values of 𝑞. It shows that tails can become

more populated for larger values of 𝑞 or under-populated for smaller values of 𝑞.

Figure 6.3 Several q-Gaussian distributions for different q parameters normalized to unit amplitude
and with a scale parameter of 𝜎 = 1. The left plot uses linear scaling, and the right plot uses
logarithmic scaling for the y-axis.

130

6.4 Static Tune Scans at Different Intensities

Section 4.5.2 showed how static tune scans are useful to verify if the resonance compensation

scheme benefits the Recycler Ring. Nevertheless, these studies in Ch. 4 were done at low intensities.

It is interesting to go to higher intensities and perform these static tune scans with a larger spread.

Figure 6.4 shows a static tune scan crossing the 3𝑄𝑥 line with no compensation at an intensity

of approximately 4.5e10 ppb (14 BTs) and a set vertical tune of 𝑄 𝑦 = 24.44. The beam survival

ratio drops dramatically after the horizontal tune of 𝑄𝑥 = 25.37, as the beam starts hitting the

resonance line. At the same time, the beam size grows exponentially. If this plot is compared to the

low-intensity case from Fig. 4.26, one can see that the losses start earlier at high intensities because

the space charge tune spread is larger for this case. The beam size growth is also enhanced at high

intensities. With large tune spreads, the beam survival ratio plots get more messy, given that the

beam can operate simultaneously on top of multiple resonance lines and there is no stable beam.

This process can be seen for horizontal tunes lower than 25.35 in Fig. 6.4. This fact is aggravated

by the transverse dampers being on for this particular experiment, but they have been fine-tuned to

operate in another tune region.

Figure 6.4 Static tune scan at an equivalent intensity of 14 Booster turns or approximately 4.5e10
ppb with no 3𝑄𝑥 compensation and transverse dampers on.

131

In Ch. 4, it was shown that when the 3𝑄𝑥 compensation was introduced, this helped improve

the beam survival ratio close to and on top of the resonance line. Figures 6.5, 6.6 and 6.7 show

what happens to the beam survival ratio and the beam size at different intensities during these static

tune scans. It is worth pointing out that the transverse dampers were turned off for these particular

measurements. A deeper discussion into the effect of transverse dampers is done in Sec. 6.5.

Additionally, the color map of the scatter plot utilized for these figures represents the decimated

turn number of the data point. Blue means that the data point happened earlier in the cycle—closer

to injection—and yellow means it happened later in the cycle. A decimated number of turns of

1024 corresponds to 65000 real turns in the machine.

Figure 6.5 Static tune scan at an equivalent intensity of 2 Booster turns or approximately 0.5e10
ppb.

The data for Fig. 6.5 has already been shown in Figs. 4.27 and 4.28. Nevertheless, it is

interesting to reference it again as the base case to compare against the high-intensity cases. In

particular, Fig. 6.5 shows the beam size and beam survival ratio as a function of horizontal tune

132

and decimated turn number. When compensated for, more than 95% of the beam survives on top

of the 3𝑄𝑥 resonance for 0.8 seconds. Nevertheless, beam size growth occurs when operating close

to a set tune of 𝑄𝑥 = 25.345, where the 3𝑄𝑥 = 76 resonance lies. This beam size growth does

not happen immediately but after many turns, as seen by the yellow color on large beam sizes.

Nevertheless, out of this particular region, there is no noticeable growth in beam size for other tune

values.

Figure 6.6 Static tune scan at an equivalent intensity of 8 Booster turns or approximately 3e10 ppb.

Figure 6.6 shows a static tune scan done at a higher intensity than the one in Fig. 6.5, i.e.,

3e10 ppb, approximately six times larger. From Fig. 6.2, the corresponding space charge tune

spread corresponds to approximately 0.015. Therefore, the beam size growth and losses should

start earlier in the scan. This process can be seen in Fig. 6.6 compared to the low-intensity case.

The other feature is the beam size growth at every tune value for the scan. The beam is injected and

undergoes some mechanism that leads to beam size growth. This mechanism is related to space

133

charge, given that it is not present at low intensities. Furthermore, it can be seen that the beam size

grows as a function of tune before leading to a big dip in the beam survival ratio. At this point, the

beam size is large enough to hit the aperture of the Recycler. Ultimately, it is interesting to see the

interplay between space charge, beam size growth, and the physical aperture of the machine.

Figure 6.7 Static tune scan at an equivalent intensity of 14 Booster turns or approximately 4.5e10
ppb.

Going almost to double the previous intensity, Fig. 6.7 shows a static tune scan done for intensity

of 5e10 ppb (particles per bunch). While the intensity is almost doubled, the tune spread only jumps

to 0.02. This is because the emittance grows exponentially, and the tune spread saturates, as shown

in Fig. 6.2. This small change in tune shift explains why Figs. 6.6 and 6.7 exhibit very similar

behaviors. Nevertheless, a closer look shows that the injection beam size is slightly higher for this

higher intensity, as expected from the emittance data shown in Fig. 6.2. Furthermore, the beam

survival ratio also decreases slightly.

134

Figure 6.8 IPM data and multi-wire first turn abort data for 𝑄𝑥 = 25.370 at different intensities.

Figure 6.8 compares emittance evolution data at these different intensities, as calculated from

Eq. 4.40. A set tune of 𝑄𝑥 = 25.37 was used for comparison. On this plot, first-turn-abort

emittances are also plotted as fuchsia triangles. Ideally, these emittances calculated from multiwire

data should correspond to those calculated from the IPM data. At low intensities, the emittances

from the IPM data are centered around the multiwire value. At higher intensities, the IPM data for

the initial number of turns coincides with the first-turn-abort data. Nevertheless, once injected, the

beam emittance grows as it circulates many turns around the ring. This growth can ultimately lead

to beam loss. Ultimately, the IPM system is a useful tool to characterize the effects of space charge.

The multi-wires could not survive such high intensities for many turns.

135

Figure 6.9 Static tune scan with beam survival ratio (BSR), average q-factor for the q-Gaussian fits,
and normalized emittance at an equivalent intensity of 2 Booster turns or approximately 0.5e10
ppb.

Going further with IPM data, one can use q-Gaussian distributions to fit raw data from the

beam profile. Using Eq. 6.2 to fit IPM data, the 𝑞 parameter can be extracted to characterize the

beam tails. Figures 6.9, 6.10, and 6.11 show a plot of the average 𝑞 values at each tune for the

different intensities. These 𝑞 values should quantify how populated the tails are. A priori, one

would not expect the tails to be underpopulated —all 𝑞 values should be above 1, 𝑞 > 1. Figures

6.9, 6.10, and 6.11 also plot the beam survival ratio and horizontal normalized emittance calculated

from IPM data. Analyzing these quantities together helps one understand the loss mechanism that

reduces the beam survival ratio.

For low intensities, it has been shown that the beam survival ratio improves when the 3𝑄𝑥

compensation is introduced. Figure 6.9 looks at how the small losses develop through the 𝑞 factor.

Close to 𝑄𝑥 = 25.345, the beam survival ratio drops slightly, while the normalized emittance and

the average 𝑞-factor increase.

In particular, the emittance increases after tens of thousands of

turns in the Recycler. While the sextupole term that feeds the 3𝑄𝑥 line has been canceled, higher

136

order terms still feed the resonance line, causing emittance growth and beam loss in this region.

Furthermore, the beam tails are populated in this region—the average 𝑞 factor increases. Therefore,

a combination of emittance growth and tail population from the resonance line explains the beam

losses in this region. It is also worth pointing out that the average 𝑞-factor for these fits is slightly

above 1.0, around 1.2, meaning that the raw data from the IPM system is not necessarily Gaussian.

This feature is not present in Figs. 6.10 and 6.11, where 𝑞 values before the resonance are centered

around 1.0. It is worth reminding the reader that the IPM system is designed to operate at higher

intensities.

Figure 6.10 Static tune scan with beam survival ratio (BSR), average q-factor for the q-Gaussian
fits, and normalized emittance at an equivalent intensity of 8 Booster turns or approximately 3e10
ppb.

According to Ref.

[30], the transverse acceptance of the Recycler Ring is 40 𝜋 mm mrad.

Taking a closer look into the emittance plot of Fig. 6.9, it can be seen that when the normalized

emittance starts hitting this value, that is where losses start to appear. Figure 6.10 shows how at

injection—for the first couple of initial turns of 3e10 ppb beam—the emittance is well below 40 𝜋

mm mrad. Nevertheless, as the beam circulates the Recycler Ring, the emittance grows close to the

acceptance of the machine. Given that this is the 95% emittance, there is still some beam outside

137

this phase space region, particularly in the tails of the beam. This beam-tail population is why

the beam survival ratio is not necessarily close to 1.0, even far from the resonance. This feature

can also be seen at a higher intensity of 4.5e10 ppb beam through Fig. 6.11. While the injection

emittance is slightly higher, the beam still grows to the machine’s acceptance limit. While the beam

emittance grows, the space charge tune spread in the Recycler is decreasing from its original value.

Close to the remnants of the resonance, there is a small dip in the beam survival ratio caused by

the increase in beam emittance. As long as the emittance comes close to the acceptance of the

machine, there will be losses from beam tail particles hitting the machine’s aperture.

Figure 6.11 Static tune scan with beam survival ratio (BSR), average q-factor for the q-Gaussian fits,
and normalized emittance at an equivalent intensity of 14 Booster turns or approximately 4.5e10
ppb.

The other important feature from the plots shown in Figs. 6.10 and 6.11 is the average 𝑞-factor

behavior. Close to 𝑄𝑥 = 25.35, the average 𝑞-factor increases. This increase coincides with a dip

in the beam survival ratio. Given that the beam is already on the acceptance limit, any increase in

the beam tail population will inevitably lead to beam loss. This behavior is what Figs. 6.10 and

6.11 are essentially showing. A point is reached where there is no more space in the beam pipe to

138

accommodate more beam tails. Therefore, any tail population mechanism will lead to beam loss.

After this point, if the beam is wide enough, there are no beam tails, given that the beam pipe itself

is collimating them, and everything is just the beam core. This is what happens for tunes left of the

resonance line. Essentially, while emittance growth exists to the limits of the machine, any beam

tail population mechanism will lead to beam loss.

6.5 Effect of Transverse Dampers

As hinted, the transverse dampers play a role in the quality of resonance compensation. The

transverse dampers at the Recycler Ring are used to suppress the coupled bunch beam instability

from the slip-stacked beam [61, 62]. This feedback system consists of a BPM pickup, input

filters, digital signal processing, output amplifiers, and stripline kickers—one set for the horizontal

direction and another for the vertical. The general knobs to tune the damper system include timing

delays and gain settings. For the experiments shown previously, the dampers were completely

turned off. Nevertheless, some experiments were done with the damper system turned on to

observe its effect. Figures 6.12 and 6.13 show a comparison between static tune scans done at low

intensity (0.5e10 ppb) with dampers on and off.

From just looking at Fig. 6.12, one could conclude that the resonance compensation does not

fully work because only 60% of the beam survives close to the resonance. Nevertheless, these

losses disappear when the dampers are turned off for Fig. 6.13. When the dampers are turned on,

they increase beam size when the horizontal tunes are close to the resonance lines. This beam size

growth leads to beam loss. The dampers seem to impact the quality of the resonance compensation

directly. Nevertheless, the gain and phase of these dampers are not optimized to run too far away

from the operational point of the Recycler Ring, i.e., 𝑄𝑥 = 25.40 and 𝑄 𝑦 = 24.44. Once the set

tune strays too far, the dampers will do more harm than good. This behavior occurs because this

feedback system will amplify any beam size growth. Further quantifiable studies should be done

to understand the interplay between dampers and resonance lines.

139

Figure 6.12 Static tune scan with beam survival ratio and IPM data box plots with 3𝑄𝑥 compensation,
transverse dampers ON and 2 Booster Turns of equivalent intensity or approximately 0.5e10 ppb.

Figure 6.13 Static tune scan with beam survival ratio and IPM data box plots with 3𝑄𝑥 compensation,
transverse dampers OFF and 2 Booster Turns of equivalent intensity or approximately 0.5e10 ppb.

140

Figure 6.14 looks at how the beam grows when the transverse dampers are activated. It can be

seen that once the horizontal tune is set to the left of 𝑄𝑥 = 25.36, a beam size and emittance growth

as a function of the number of turns starts to appear. This feature was absent in its dampers-off

counterpart, as shown in Figs. 6.5 and 4.28. Figure 6.14 also shows how tails generate before the

beam is lost to the beam pipe. Nevertheless, there seems to be a tune shift with respect to where the

losses are now generated. It is also worth noting that the raw data from the IPM follows this trend

from Fig. 4.28, where the average q-factor is around 1.2, meaning the data is not truly Gaussian at

low intensities.

Figure 6.14 Static tune scan with beam survival ratio (BSR), average q-factor for the q-Gaussian
fits, and normalized emittance with 3𝑄𝑥 compensation ON, transverse dampers ON and 2 Booster
Turns of equivalent intensity or approximately 0.5e10 ppb.

141

CHAPTER 7

CONCLUSIONS AND FUTURE WORK

The work shown throughout this document can also be found in Refs. [63–65]. This section reports

on the conclusions and future work stemming from this thesis.

7.1 Conclusions

7.1.1 RDTs and Resonance Compensation in the Recycler Ring

The measurement and subsequent cancellation of RDTs were investigated in depth in Ch. 4. It

was shown that third-order RDTs in the Recycler Ring could be measured and, afterward, canceled

with the compensation sextupoles in place. Furthermore, optimization procedures have been

developed to have this compensation with physics-informed input and/or numerical optimization.

The Recycler Ring now has a resonance compensation method for when large tune spreads are

present. Depending on the operational tunes, any setting in the compensation sextupoles from

Fig. 4.23 can be enabled to cope with space-charge-induced losses. This advancement is a new

operational feature that can be enabled on demand and will be especially useful for the PIP-II era.

7.1.2 Physics-Informed vs. Optimization-Based Compensation

Chapters 4 and 5 showed two different approaches to the same problem of resonance com-

pensation. The first one performed at the Fermilab Recycler Ring showed how to implement the

response matrix approach, which can be traced to some underlying physics principles. On the

other hand, Ch. 5 shows how to implement numerical optimization algorithms for compensating

multiple resonance lines at the CERN PS Booster. Nevertheless, this last approach had no direct

physics input. The distinction between a physics-informed approach and a numerical-optimization

approach summarizes the differences between both chapters. The numerical-optimization approach

is a brute-force one that can be carried out without any physics input—just by minimizing the losses

in the resonance lines with the compensation sextupoles. Nevertheless, the physics-informed com-

pensation can significantly enhance numerical optimization by constraining a window in parameter

space where the solution lies. The best approach to resonance compensation is to use a hybrid

scheme where the response matrix method sets the limits in parameter space for subsequent and

142

faster numerical optimization.

7.1.3 High-Intensity Resonance Compensation

For current operations of the Fermilab Accelerator Complex, the third-order resonance com-

pensation scheme developed in this work will not have an effect on reducing losses in the Recycler

Ring. The current tune spreads in the RR need to be larger for the third-order resonances to be

especially harmful to operations. The beam coming out of Booster grows exponentially with the

beam intensity, as shown in Fig. 6.2. Thus, saturating the tune spread at values close to 0.02—not

particularly large. The losses in the Recycler Ring are purely emittance-dominated for current

operations. They are not related to space charge effects. Nonetheless, the resonance compensation

scheme still benefits operations at high intensities when close to the third-order resonances. The

operational tune should remain as far as possible from these lines.

7.1.4 Transverse Dampers and High-Intensity Compensation

Section 6.5 showed how the transverse dampers have an effect on the resonance compensation

scheme. Figure 7.1 reinforces this statement by showing the beam survival ratio at a high intensity

of 4.5e10 ppb without compensation and dampers on (bare machine), with compensation and

dampers on, and with compensation and dampers off. Ultimately, the dampers degrade the quality

of the compensation. In particular, when turned off, the overall beam survival ratio increases. The

configuration of the dampers—gain and phase of feedback—has been optimized away from the

resonance lines. The damper settings will create beam size growth when the tunes are set close to

the third-order resonance lines. They are not optimized to run in this region and will cause harmful

effects on the beam, even with the resonance compensation scheme enabled.

143

Figure 7.1 Transverse damper effect on resonance compensation for bare machine + dampers ON,
3𝑄𝑥 compensated + dampers ON, and 3𝑄𝑥 compensated + dampers OFF configuration. The red
band corresponds to the 3𝑄𝑥 stop band at an intensity of 0.5 ppb with dampers ON.

7.2 Future Work

7.2.1 Verification of Newly-Installed Sextupoles

Section 4.3 explained the motivation behind installing additional sextupoles that would bring

down the currents needed to compensate 3𝑄𝑥 = 76 and 𝑄𝑥 + 2𝑄 𝑦 = 74. Subsequently, Sec. 4.6

explained the procedure used to pin down the new locations for two new compensation sextupoles.

These sextupoles have been installed in the Recycler. There is still future work to be done

regarding the commissioning and connection of the new 620 sextupoles. The sextupoles and their

power supplies need to be interfaced with ACNET. Furthermore, the resonance compensation

enhancement still needs to be verified by performing an RDT scan. Specifically, the response

matrix coefficients corresponding to these new sextupoles must be measured and calculated. All of

this, following the procedure outlined in Secs. 4.2 and 4.3. Figure 7.2 shows a picture of the newly

installed sextupoles.

144

Figure 7.2 Newly installed compensation sextupoles in the 620 section.

7.2.2 Resonances and Transverse Dampers at High Intensities

It was shown in Sec. 6.5 that the transverse dampers in the Recycler Ring can interact with

betatron resonance lines depending on their configuration.

It would be interesting to perform

experiments with the configuration of these dampers to characterize these effects fully, i.e., change

the gain and phase of the dampers and study their effect on the strength of the resonance lines.

Furthermore, it would also be interesting to use the dampers as anti-dampers to sustain betatron

oscillations. With this configuration, one could perform tune measurements and evaluate this

alternative method to measure RDTs. This tool would be similar to installing an AC dipole in the

Recycler.

145

7.2.3 Effect of MI Ramp on RDTs

The ultimate objective of this work is to characterize this resonance compensation scheme

fully and incorporate it into high-intensity operations. Nevertheless, one additional factor has

also been identified as playing a role in this effort. This factor stems from the fact that the Main

Injector and the Recycler Ring share the same tunnel. It has been shown that the MI acceleration

ramp changes the beam dynamics inside the RR. In particular, it introduces orbit distortions and

tune shifts depending on the ramp position [66]. There is an ongoing effort to characterize any

higher-order magnetic effect from the MI to RR, e.g., any sextupole term introduced by the MI

acceleration ramp. The first results have shown that the compensation currents change depending

on the location of the study event with respect to the MI acceleration ramp. Therefore, in the future,

the resonance compensation described hereinabove should be modified to accommodate this effect

to be fully operational.

7.2.4 Limits of Resonance Compensation

Chapter 4, specifically Sec. 4.5.1, shows how individual and different pairs of resonance lines

can be compensated in the Recycler Ring. Nevertheless, there is still a question of whether all four

third-order resonance lines can be compensated simultaneously. While the previous thesis does not

explicitly answer this question, it does contribute to a better understanding of the resonance lines.

The source of these third-order resonance lines is a systematic third-order term in the permanent

combined function magnets distributed around the ring. The fact that these are systematic does not

allow for a straightforward compensation of the four resonance lines. Whenever more resonance

lines try to be corrected, the currents in the compensation sextupoles increase. Therefore, one

would need more compensation sextupoles around the ring at specific locations to bring down these

currents, just as it was done in Sec. 4.6. Ultimately, the approach to follow would be to try to

minimize the RDT fluctuations around the ring with the least amount of correctors. Such is the

approach in Ref. [67]. Another thing to consider for future theoretical studies is that introducing

too much sextupole in the ring might take the machine closer to an anharmonic regime. In this

anharmonic regime, the resonances will not necessarily be straight lines but rather polynomial

146

curves resulting from motion in a 4D coupled phase space [68, 69]. These are considerations for

future exploration of the limits of resonance compensation driven by systematic terms around the

ring.

7.2.5 Space Charge RDTs

Equation 2.80 shows how the lattice elements and the space charge potential drive betatron

resonances.

It would be interesting to explore further how to measure space charge resonance

driving terms (SCRDTs) and their effect on the operation of the Recycler Ring. Without taking

into account any collective instabilities, the ultimate limit of the Recycler Ring will be dictated

by how strong these third-order lines become in the space-charge-dominated regime. The current

future of the Main Injector and Recycler Ring lies in the high-intensity beam. As the space charge

potential grows for the Fermilab beams, it is important to quantify all effects from high-intensity

operation—being SCRDTs one of them.

147

BIBLIOGRAPHY

[1] S Y Lee. Accelerator Physics.

3rd. WORLD SCIENTIFIC, 2011. doi: 10.1142/
8335. eprint: https://www.worldscientific.com/doi/pdf/10.1142/8335. url: https://www.
worldscientific.com/doi/abs/10.1142/8335.

[2] Lillian Hoddeson, Adrienne W. Kolb, and Catherine Westfall. Fermilab: Physics, the Fron-
isbn: 9780226346236.
tier, and Megascience. University of Chicago Press, Dec. 2008.
doi: 10.7208/chicago/9780226346250.001.0001. url: https://doi.org/10.7208/chicago/
9780226346250.001.0001.

[3] David Griffiths.

“Historical Introduction to the Elementary Particles”.

Introduc-
tion to Elementary Particles.
John Wiley and Sons, Ltd, 1987. Chap. 1, pp. 11–53.
isbn: 9783527618460. doi: https://doi.org/10.1002/9783527618460.ch1. url: https:
//onlinelibrary.wiley.com/doi/abs/10.1002/9783527618460.ch1.

In:

[4] Michigan State University. Facility for Rare Isotope Beams at Michigan State University.

https://frib.msu.edu/index.php.

[5] SLAC National Accelerator Laboratory. The Stanford Linear Accelerator Center at SLAC.

https://www.slac.stanford.edu/gen/grad/GradHandbook/slac.html.

[6] S. Holmes et al. “Accelerator physics at the Tevatron Collider: Introduction”. In: Accelerator
physics at the Tevatron Collider. Ed. by Valery Lebedev and Vladimir Shiltsev. 2014, pp. 1–
28. doi: 10.1007/978-1-4939-0885-1_1.

[7] CERN. The Large Hadron Collider (LHC). https://home.cern/science/accelerators/large-

hadron-collider.

[8] Ferdinand Willeke and J. Beebe-Wang. Electron Ion Collider Conceptual Design Report
2021. Tech. rep. BNL-221006-2021-FORE. Brookhaven National Laboratory, Feb. 2021.
doi: 10.2172/1765663. url: https://www.osti.gov/biblio/1765663.

[9] Étienne Forest. Beam Dynamics: A New Attitude and Framework. Vol. 8. The Physics
and Technology of Particle and Photon Beams. Amsterdam, The Netherlands: Hardwood
Academic / CRC Press, 1998. isbn: 978-90-5702-574-7.

[10] E Kako, P Pierini, and C.E. Reece. https://cerncourier.com/a/teslas-high-gradient-march/.

Dec. 2020.

[11] Reidar Hahn. Fermilab. File: https://upload.wikimedia.org/wikipedia/commons/3/3f/Fermilab.jpg.

2003. url: https://en.wikipedia.org/wiki/File:Fermilab.jpg.

[12] Swapan Chattopadhyay and Joseph David Lykken. Fermilab at 50. WORLD SCIEN-
TIFIC, 2017. doi: 10.1142/10637. eprint: https://www.worldscientific.com/doi/pdf/10.

148

1142/10637. url: https://www.worldscientific.com/doi/abs/10.1142/10637.

[13] P. Adamson et al. “The NuMI neutrino beam”. In: Nuclear Instruments and Methods in
Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equip-
ment 806 (2016), pp. 279–306. issn: 0168-9002. doi: https://doi.org/10.1016/j.nima.2015.
08.063. url: https://www.sciencedirect.com/science/article/pii/S016890021501027X.

[14] Andrzej Wolski. Beam Dynamics in High Energy Particle Accelerators. 2nd. WORLD
SCIENTIFIC, 2023. doi: 10.1142/13333. url: https://www.worldscientific.com/doi/abs/
10.1142/13333.

[15] Helmut Wiedemann. Resonances. Cham: Springer International Publishing, 2015, pp. 539–
isbn: 978-3-319-18317-6. doi: 10.1007/978-3-319-18317-6_16. url: https://doi.

564.
org/10.1007/978-3-319-18317-6_16.

[16] T. Satogata. Lecture Notes on Poisson Brackets and Lie Operators. https://toddsatogata.net/

2011-USPAS/lieAlgebras.pdf. Jan. 2008.

[17] P. Urschutz. “Measurement and Compensation of Betatron Resonances at the CERN PS

Booster Synchrotron”. PhD thesis. Vienna, Austria, 2004.

[18] R. Tomas Garcia. “Direct Measurement of Resonance Driving Terms in the Super Proton
Synchrotron (SPS) of CERN using Beam Position Monitors”. PhD thesis. Valencia, Spain,
Jan. 2003.

[19] Rüdiger Achilles and Andrea Bonfiglioli. “The early proofs of the theorem of Campbell,
Baker, Hausdorff, and Dynkin”.
In: Archive for History of Exact Sciences 66.3 (May
2012), pp. 295–358. doi: 10.1007/s00407-012-0095-8. url: https://doi.org/10.1007/
s00407-012-0095-8.

[20] R Bartolini and F Schmidt. “Normal form via tracking or beam data”. In: Part. Accel. 59

(1998), pp. 93–106. url: https://cds.cern.ch/record/333077.

[21] MAD-X: Methodical Accelerator Design. http://madx.web.cern.ch/madx.

[22] F. Asvesta. PySCRDT. 2019.

[23] Kouichi Soutome and Hitoshi Tanaka. “Higher-order formulas of amplitude-dependent
In: Phys. Rev. Accel.
tune shift caused by a sextupole magnetic field distribution”.
Beams 20 (6 June 2017), p. 064001. doi: 10.1103/PhysRevAccelBeams.20.064001. url:
https://link.aps.org/doi/10.1103/PhysRevAccelBeams.20.064001.

[24] A. W. Chao. Physics of collective beam instabilities in high-energy accelerators. 1993.

isbn: 978-0-471-55184-3.

149

[25] Ronald C Davidson and Hong Qin. Physics of Intense Charged Particle Beams in High En-
ergy Accelerators. IMPERIAL COLLEGE PRESS CO-PUBLISHED WITH WORLD SCI-
ENTIFIC PUBLISHING CO, 2001. doi: 10.1142/p250. url: https://www.worldscientific.
com/doi/abs/10.1142/p250.

[26] F. Asvesta and H. Bartosik. Resonance Driving Terms From Space Charge Potential. Tech.
rep. CERN-ACC-NOTE-2019-0046. https://cds.cern.ch/record/2696190/files/CERN-ACC-
NOTE-2019-0046.pdf. 2019.

[27] George B. Arfken, Hans J. Weber, and Frank E. Harris. “Chapter 10 - Green’s Functions”. In:
Mathematical Methods for Physicists (Seventh Edition). Ed. by George B. Arfken, Hans J.
Weber, and Frank E. Harris. Seventh Edition. Boston: Academic Press, 2013, pp. 447–467.
isbn: 978-0-12-384654-9. doi: https://doi.org/10.1016/B978-0-12-384654-9.00010-4.

[28] Gabriel Tellez Acosta. Metodos Matemeticos. Ediciones Uniandes, 2022. isbn: 9789587981827.

[29] Oliver Dimon Kellog. “Foundations of potential theory”. In: Springer, 1967, pp. 192–194.

[30] G. Jackson. The Fermilab recycler ring technical design report. Revision 1.2. Tech. rep.
FNAL-TM-1991ON: DE97051388; BR: KAHEP. Fermilab, Nov. 1996. doi: 10.2172/
16029. url: https://www.osti.gov/biblio/16029.

[31] S. Nagaitsev. Fermilab Antiproton source, Recycler ring and Main Injector. Tech. rep.
FERMILAB-FN-0957-AD. Fermilab, Mar. 2013. doi: 10.2172/1127965. url: https:
//www.osti.gov/biblio/1127965.

[32] R. Ainsworth et al. “High intensity operation using proton stacking in the Fermilab Recycler
to deliver 700 kW of 120 GeV proton beam”. In: Phys. Rev. Accel. Beams 23 (12 Dec.
2020), p. 121002. doi: 10.1103/PhysRevAccelBeams.23.121002. url: https://link.aps.org/
doi/10.1103/PhysRevAccelBeams.23.121002.

[33] M. Ball et al. The PIP-II Conceptual Design Report. Tech. rep. FERMILAB-TM-2649-AD-
APC1516858. Fermilab, Mar. 2017. doi: 10.2172/1346823. url: https://www.osti.gov/
biblio/1346823.

[34] Fermi National Accelerator Laboratory. Deep Underground Neutrino Experiment. url:

https://www.dunescience.org/.

[35] Fermilab Accelerator Division — Operations Department. Concepts Rookie Book.

[36] M. Xiao. “Phase trombone program migration for the recycler ring at Fermilab.” In: Pro-
ceedings of the IEEE Particle Accelerator Conference. Vol. 2005. Proceedings of the
Particle Accelerator Conference, PAC 2005. Fermilab, 2005, pp. 3135-3137 –3137.

150

[37] R. Ainsworth et al. “High intensity space charge effects on slip stacked beam in the Fermi-
lab Recycler”. In: Phys. Rev. Accel. Beams 22 (2 Feb. 2019), p. 020404. doi: 10.1103/
PhysRevAccelBeams.22.020404. url: https://link.aps.org/doi/10.1103/PhysRevAccelBeams.
22.020404.

[38] R. Webber et al. “Fermilab Recycler Ring BPM Upgrade Based on Digital Receiver Technol-
ogy”. In: AIP Conference Proceedings 732.1 (Nov. 2004), pp. 190–200. issn: 0094-243X.
doi: 10.1063/1.1831147. url: https://doi.org/10.1063/1.1831147.

[39]

J. Laskar. “The chaotic motion of the solar system: A numerical estimate of the size of
issn: 0019-1035. doi: https:
the chaotic zones”. In: Icarus 88.2 (1990), pp. 266–291.
//doi.org/10.1016/0019-1035(90)90084-M. url: https://www.sciencedirect.com/science/
article/pii/001910359090084M.

[40] B. Babacan et al. The Ionization Profile Monitors in the Recycler Ring. 2023. arXiv:

2307.07017.

[41] M. Xiao. Recycler Chromaticities and End Shims for NOvA at Fermilab. 2013. arXiv:

1301.6666.

[42] SUSSIX: A Computer Code for Frequency Analysis of Nonlinear Betatron Motion.

[43] R. E. Meller et al. “Decoherence of Kicked Beams”. In: (May 1987).

[44] S. Y. Lee. “Decoherence of the Kicked Beams II”. in: (Feb. 1991).

[45] Ming-Jen Yang. Lattice function measurement with TBT BPM data. Tech. rep. Fermilab,

1996. url: https://lss.fnal.gov/archive/test-tm/1000/fermilab-tm-1922.pdf.

[46] Robert Ainsworth et al. “Improvements to the Recycler/Main Injector to Deliver 850 kW+”.
In: JACoW (Oct. 2022). doi: 10.18429/JACoW-NAPAC2022-WEYE3. url: https://www.
osti.gov/biblio/1909888.

[47] Roman Garnett. Bayesian Optimization. Cambridge University Press, 2023.

[48] Fabienne Landua. A graphic overview of all accelerators in operation at CERN.. url:

https://cds.cern.ch/record/2813716?ln=en.

[49] H Damerau et al. LHC Injectors Upgrade, Technical Design Report. 2014. doi: 10.17181/

CERN.7NHR.6HGC. url: https://cds.cern.ch/record/1976692.

[50] Wim Weterings et al. “The new injection region of the CERN PS Booster”. In: (2019),
WEPMP039. doi: 10.18429/JACoW-IPAC2019-WEPMP039. url: https://cds.cern.ch/
record/2693918.

151

[51] Foteini Asvesta et al. “Resonance Compensation for High Intensity and High Brightness
Beams in the CERN PSB”. in: JACoW HB 2021 (2022), pp. 40–45. doi: 10.18429/
JACoW-HB2021-MOP06. url: https://cds.cern.ch/record/2841816.

[52] Tirsi Prebibaj. “Optimization of the High Brightness Beam Performance of the CERN PSB
with H- Injection”. PhD thesis. Frankfurt am Main, Germany: 2023-06-22, June 2023.

[53] Foteini Asvesta. “Space charge and lattice driven resonances at the CERN injectors”. PhD

thesis. Natl. Tech. U., Athens, 2021.

IPAC’22 (Bangkok, Thailand).

[54] F. Asvesta et al. “High Intensity Studies in the CERN Proton Synchrotron Booster”.

In:
International Particle Accelerator Conference 13.
Proc.
JACoW Publishing, Geneva, Switzerland, July 2022, WEPOTK011, pp. 2056–2059. isbn:
978-3-95450-227-1. doi: 10.18429/JACoW-IPAC2022-WEPOTK011. url: https://jacow.
org/ipac2022/papers/wepotk011.pdf.

[55] Simon Albright et al. “New Longitudinal Beam Production Methods in the CERN Proton
Synchrotron Booster.” In: JACoW IPAC ’21 (2021), pp. 4130–4133. doi: 10.18429/
JACoW-IPAC2021-THPAB183. url: https://cds.cern.ch/record/2811646.

[56] Generic Optimization Frontend and Framework (GeOFF). https://gitlab.cern.ch/geoff/geoff-

app.

[57] M. J. D. Powell. “The BOBYQA algorithm for bound constrained optimization without

derivatives”. In: 2009. url: https://api.semanticscholar.org/CorpusID:2488733.

[58] S.Y. Zhang. CALCULATION OF INCOHERENT SPACE CHARGE TUNE SPREAD.
tech. rep. Brookhaven National Laboratory, 1996. url: https://technotes.bnl.gov/PDF?
publicationId=30778.

[59] S. Papadopoulou et al. “Impact of non-Gaussian beam profiles in the performance of hadron
colliders”. In: Phys. Rev. Accel. Beams 23 (10 Oct. 2020), p. 101004. doi: 10.1103/
PhysRevAccelBeams.23.101004. url: https://link.aps.org/doi/10.1103/PhysRevAccelBeams.
23.101004.

[60] Wikipedia contributors. Q-Gaussian distribution — Wikipedia, The Free Encyclopedia.
[Online; accessed 19-March-2024]. 2023. url: https://en.wikipedia.org/w/index.php?title=
Q-Gaussian_distribution&oldid=1177637416.

[61] L. R. Prost et al. Transverse Instabilities in the FERMILAB Recycler. 2012. arXiv: 1201.

4168.

[62] Nathan Eddy et al. “Transverse Damper Using Diodes for Slip Stacking in the Fermilab
Recycler”. In: 6th International Beam Instrumentation Conference. 2018, TUPCF21. doi:
10.18429/JACoW-IBIC2017-TUPCF21.

152

[63] C.E. Gonzalez-Ortiz, R. Ainsworth, and P.N. Ostroumov. “Third-Order Resonance Com-
IPAC’22 (Bangkok, Thailand).
In-
pensation at the FNAL Recycler Ring”.
JACoW Publishing, Geneva, Switzer-
ternational Particle Accelerator Conference 13.
land, July 2022, MOPOST050, pp. 195–198.
isbn: 978-3-95450-227-1. doi: 10.18429/
JACoW-IPAC2022-MOPOST050. url: https://jacow.org/ipac2022/papers/mopost050.pdf.

In: Proc.

[64] C.E. Gonzalez-Ortiz, R. Ainsworth, and P.N. Ostroumov. “Simultaneous Compensation
IPAC’23 (Venice).
of Third-Order Resonances at the FNAL Recycler Ring”.
IPAC’23 - 14th International Particle Accelerator Conference 14.
JACoW Publishing,
Geneva, Switzerland, May 2023, pp. 3328–3331. isbn: 978-3-95450-231-8. doi: doi:10.
18429/jacow-ipac2023-wepl112. url: https://www.ipac23.org/preproc/pdf/WEPL112.pdf.

In: Proc.

[65] C.E. Gonzalez-Ortiz, R. Ainsworth, and P.N. Ostroumov. “Compensation of Third-order
Resonances in the High Intensity Regime”. In: Proc. HB’23 (Geneva). HB’23. JACoW
Publishing, Geneva, Switzerland, Oct. 2023. url: https://hb2023.vrws.de/papers/wea2i1.
pdf.

[66] N. Chelidze, R. Ainsworth, and K.J. Hazelwood. “The Effect of the Main Injector Ramp
on the Recycler”. In: Proc. 5th Int. Particle Accel. Conf. (NAPAC’22) (Albuquerque,
NM, USA). International Particle Accelerator Conference 5. JACoW Publishing, Geneva,
Switzerland, Oct. 2022, MOPA19, pp. 90–92. isbn: 978-3-95450-232-5. doi: 10.18429/
JACoW-NAPAC2022-MOPA19. url: https://jacow.org/napac2022/papers/mopa19.pdf.

[67] Bingfeng Wei et al. “Minimizing the fluctuation of resonance driving terms in dynamic aper-
ture optimization”. In: Phys. Rev. Accel. Beams 26 (8 Aug. 2023), p. 084001. doi: 10.1103/
PhysRevAccelBeams.26.084001. url: https://link.aps.org/doi/10.1103/PhysRevAccelBeams.
26.084001.

[68] H. Bartosik, G. Franchetti, and F. Schmidt. “Observation of fixed lines induced by a nonlinear

resonance in the CERN Super Proton Synchrotron”. In: Nature Physics (2024).

[69] Ezio Todesco. “Analysis of resonant structures of 4D symplectic mappings, using normal
forms”. In: Phys. Rev. E 50.6 (1994), pp. 4298–4301. doi: 10.1103/PhysRevE.50.R4298.
url: https://cds.cern.ch/record/262010.

153