‘I

but hIUI II ‘I I I

Michigan State
W University

 

This is to certify that the
dissertation entitled

OPTIMIZATION OF ACCELERATOR PARAMETERS USING
NORMAL FORM METHODS ON HIGH-ORDER TRANSFER
MAPS

presented by

PAVEL SNOPOK

has been accepted towards fulfillment
of the requirements for the

Ph.D. degree in Physics and Mathematics

 

 

Ravi/@— M

Major Professof§Signature
S I ‘l I o 7

Date

 

MSU is an affirmative-action, equal-opportunity employer

 

 

-.—.—.-.—<—-—-—--¢-I-I-O-n-l-l-O-O-l-I-I-u-I-.-o-l-u-0-l-l-t-I-I--Q-l-l-I-l-l-I-U-O-l-o-o-I—I-o-I-I-C-O-I-o-i-l-t-I-I-I-I-

 

OPTIMIZATION OF ACCELERATOR
PARAMETERS USING NORMAL FORM
METHODS ON HIGH-ORDER TRANSFER
MAPS

By

Pavel Snopok

A DISSERTATION

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

DOCTOR OF PHILOSOPHY

Department of Physics and Astronomy

2007

 

ABSTRACT

OPTIMIZATION OF ACCELERATOR
PARAMETERS USING NORMAL FORM
METHODS ON HIGH-ORDER TRANSFER
MAPS

By

Pavel Snopok

Methods of analysis of the dynamics of ensembles of charged particles in collider
rings are developed. The following problems are posed and solved using normal form

transformations and other methods of perturbative nonlinear dynamics:
0 Optimization of the Tevatron dynamics:
— Skew quadrupole correction of the dynamics of particles in the Tevatron
in the presence of the systematic skew quadrupole errors in dipoles;
- Calculation of the nonlinear tune shift with amplitude based on the results
of measurements and the linear lattice information;
0 Optimization of the Muon Collider storage ring:
— Computation and optimization of the dynamic aperture of the Muon Col-
lider 50x50 GeV storage ring using higher order correctors;

— 750x 750 GeV Muon Collider storage ring lattice design matching the Teva-

tron footprint.

 

The normal form coordinates have a very important advantage over the particle
optical coordinates: if the transformation can be carried out successfully (general
restrictions for that are not much stronger than the typical restrictions imposed on the
behavior of the particles in the accelerator) then the motion in the new coordinates has
a very clean representation allowing to extract more information about the dynamics
of particles, and they are very convenient for the purposes of visualization.

All the problem formulations include the derivation of the objective functions,
which are later used in the optimization process using various optimization algorithms.
Algorithms used to solve the problems are specific to collider rings, and applicable to
similar problems arising on other machines of the same type.

The details of the long-term behavior of the systems are studied to ensure the
their stability for the desired number of turns. The algorithm of the normal form
transformation is of great value for such problems as it gives much extra information
about the disturbing factors. In addition to the fact that the dynamics of particles
is represented in a way that is easy to understand, such important characteristics
as the strengths of the resonances and the tune shifts with amplitude and various
parameters of the system are calculated.

Each major section is supplied with the results of applying various numerical
optimization methods to the problems stated. The emphasis is made on the efficiency
comparison of various approaches and methods. The main simulation tool is the
arbitrary order code COSY INFINITY written by M. Berz, K. Makino, et al. at
Michigan State University. Also, the code MAD is utilized to design the 750x 750 GeV
Muon Collider storage ring baseline lattice.

The OptiM to COSY lattice converter is written specifically for the need of the

studies included into the dissertation, and tested on the Tevatron accelerator lattice.

Copyright by
Pavel Snopok

2007

 

i’I

 

Dedicated to my beloved Olga

 

ACKNOWLEDGMENTS

The acknowledgments page is really the opportunity to slow down for a moment,
to look back and observe the whole aeon of the graduate program studies, to recall
the names of the persons who helped me a great deal in reaching my goal in the form
of this dissertation piece. Without these people and their help the task of the Ph.D.
program would be next to impossible to accomplish.

For sure, the last four year were not the easiest for me, but definitely the most
interesting in the sense of new impressions, acquaintances, professional growth and
gaining invaluable experience. I am thankful to all the people who worked hard
on giving me the opportunity to pursue the graduate program at Michigan State
University: Dr. Martin Berz, my scientific adviser at MSU; Dr. Carol Johnstone,
my adviser at Fermi National Accelerator Laboratory; Dr. Dmitry Ovsyannikov and
Dr. Alexander Ovsyannikov, my scientific advisers at St.Petersburg State University,
both in undergraduate and graduate level programs; Dr. Viktor Yarba from Fermi
National Accelerator Laboratory.

I personally think it is very useful to be involved in a joint program, where the
participants can study at the university and work on their research at some national
laboratory. It gives you broader exposure and more areas to explore in addition to
fairly different working environments.

Graduate program at MSU/FNAL gave me a lot: the vision of the scientific
process as it is in the international scientific community; the opportunity to do my
research at Fermilab, one of the biggest accelerator physics centers in the world; the
multinational multicultural team work experience.

Besides that, I gained priceless experience of applying my knowledge and expertise

in Mathematics to the field of Beam Physics, to solving the theoretical and applied

vi

 

problems connected to the optimization of the parameters of the accelerating struc-
tures.

I hope during the course of my graduate studies I acquired the necessary experi-
ence to conduct independent research, scientific studies; to teach both undergraduate
and graduate level courses. This dissertation work should serve as a proof of my com-
petence. I also hope that the results derived in the dissertation will be interesting to
someone else, and not only me, and serve as a source of inspiration on how to improve
the methods of accelerator optimization proposed in the work.

I’m thankful to my fellow students from the Beam Theory group, who were always
willing to give a helping hand no matter what the situation was: Dr. Shashikant
Manikonda, Johannes Grote, Alex Wittig, Youn-Kyung Kim; and even more so to
my brother in arms Alexey Poklonskiy, with whom I went through the'graduate
studies side by side.

Dr. Kyoko Makino (MSU), Dr. David Neuffer (Fermilab), Dr. Kevin Paul
(Muons, Inc.), Dr. Michael Syphers (Fermilab) I thank for interesting discussions,
good advice and willingness to help at all stages of my research.

I cannot express how much I’m obliged for everything that happened in my life

 

and scientific career to my parents —— Wyacheslav Snopok] and Anna Snopok. I am

 

so sorry my father will not see me a Ph.D., but I’m sure he would be proud of me.
Endless thanks for love, constant support and boundless patience to my beloved
one —— Olga Safronova. I dedicate this work to her.
Unfortunately, I cannot go on mentioning people I’m thankful to, otherwise this
acknowledgment section will grow to the size of the dissertation itself. So, many
thanks to everyone taking part in my scientific and private life, I hope these people

will easily recognize themselves; a profound bow and my sincere gratitude.

vii

 

TABLE OF CONTENTS

LIST OF TABLES ................................ x
LIST OF FIGURES ............................... xi
1. Introduction .................................. 1
2. Equations of Motion, Fields and Potentials .............. 12
2.1 Equations of Motion of the Particles in Accelerators .......... 12
2.1.1 Curvilinear Coordinates ..................... 12

2.1.2 Equations of Motion in Curvilinear Coordinates ........ 13

2.2 Fields and Potentials ........................... 18
2.2.1 Fields and Potentials of the Magnetic Multipole (Stationary Case) 20

2.2.2 Multipole Expansion in Cylindrical Coordinates ........ 22

2.2.3 Multipole Expansion in Cartesian Coordinates ......... 24

3. Normal Form Algorithms ......................... 27
3.1 Nonlinear Normal Form Transformation ................ 27
3.2 Stable Symplectic Case, Nonlinear Tune Shift with Amplitude . . . . 31
3.3 Motivation for the Differential Algebra Approach ........... 36
3.4 Normal Form Transformation Result .................. 42

4. Dynamics Optimization in the Tevatron ................ 46
4.1 Skew Quadrupole Correction ....................... 46
4.1.1 Linear Problem .......................... 50

4.1.2 Nonlinear Problem ........................ 51

4.1.3 Time Preservation ......................... 52

4.1.4 Introduction to the Solution Methods .............. 53

4.1.5 Skew Quadrupole Error Sources ................. 55

4.1.6 Initial Data ............................ 56

4.1.7 Model Details and Tracking Results ............... 58

4.1.8 Skew Quadrupole Circuits Optimization Proposals ....... 59

4.1.9 Optimization Results ....................... 64

4.1.10 Summary ............................. 69

4.2 Calculation of Nonlinear Tune Shifts with Amplitude ......... 70
4.2.1 Objective Functions for the Study ................ 73

4.2.2 Calculation Results versus Measurement Results ........ 74

4.2.3 Sector Approximation, Uniform Particle Distribution ..... 78

4.2.4 Elliptical Beam, Uniform Distribution ............. 84

viii

4.2.5 Elliptical Beam, Normal or Arbitrary Distribution ....... 89

4.2.6 Numerical Experiment Results .................. 91
4.2.7 Summary ............................. 94
5. Muon Collider Storage Rings: Dynamics and Optimization . . . . 95
5.1 50x50 GeV Collider IR High-Order Correction ............. 95
5.1.1 Dynamic Aperture ........................ 97
5.1.2 Minimization of the Nonlinearities ................ 99
5.1.3 Minimization of the Resonances ................. 100
5.1.4 Dynamic Aperture as the Objective Function ......... 101
5.1.5 Advantages and Disadvantages of Different Optimization Ap-
proaches .............................. 102
5.1.6 Optimization Results ....................... 106
5.1.7 Summary ............................. 110
5.2 750x750 GeV Collider Storage Ring Lattice for the Tevatron Main Ringlll
5.2.1 Introduction ............................ 111
5.2.2 750x750 GeV Lattice ....................... 112
5.2.3 Final Focus Section ........................ 114
5.2.4 Chromatic Correction Section .................. 114
5.2.5 Arc Module ............................ 116
5.2.6 Advantages and Disadvantages of the New Design ....... 117
APPENDICES .................................. 120
A. Tevatron Lattice Converter from OptiM to COSY INFINITY . . 120
A2 Introduction ................................ 120
A3 Tools for Conversion ........................... 121
A.4 Source Code Details ........................... 121
A5 COSY Code Details in Connection to OptiM Code .......... 122
A6 Summary of Functionality of the Resulting
COSY Output ............................... 124
A.7 Code List of Element Wrappers ..................... 124
B. Graphical User Interface for the Tevatron Simulations in COSY IN-
FINITY .................................... 131
BIBLIOGRAPHY ................................ 136

ix

 

3.1
3.2

5.1
5.2

LIST OF TABLES

Transfer map in particle optical coordinates .............. 44
'Itansfer map in normal form coordinates ................ 45
Parameter comparison for various optimization approaches ...... 106
Parameter comparison for various storage rings lattices ........ 114

 

2.1

4.1
4.2
4.3

4.4

4.5

4.6

4.7

4.8

4.9

4.10

4.11
4.12

LIST OF FIGURES

Axes orientation for the cylindrical coordinates .............

Lattice scheme before any correction ..................

22

54

Beta functions comparison, top - OptiM32, bottom - COSY INFINITY 57

The (:17, a) plane and (y, b) plane phase portraits, only dipoles and pure
quadrupoles are active, and particles are launched along the x axis and
the y axis, respectively. The scales are 2.7- 10"3 m for the a: and y
axes and 4.0 - 10"3 for the a. and b axes. ................
The (as, a) plane and (y, b) plane phase portraits before the optimiza-
tion, particles are launched along the a: axis and the y axis, respec-
tively. The phase portraits include all sextupole and skew quadrupole
fields (correctors plus errors) in addition to quadrupoles and dipoles,
where 15% of the skew quadrupole errors have been removed in specific
dipoles. the scales are 2.7- 10"3 m for the a: and y axes and 4.0- 10‘3
for the a. and b axes. ...........................
Correction scheme 1: the skew quadrupole error is removed in dipoles
in each cell surrounding the defocusing quadrupole. ..........
Correction scheme II: the skew quadrupole errors in dipoles are fixed
in each even cell. .............................
The (11:,0.) plane and (y, b) plane phase portraits after the optimization
with 50% skew quadrupole errors in dipoles, errors are fixed according
to Scheme I: around each defocusing quadrupole; particles are launched
along the :1: axis and the y axis, respectively. The scales are 2.7- 10‘3
In for the :1: and y axes and 4.0 - 10‘3 for the a. and b axes. ......
The (:13, a) plane and (y, b) plane phase portraits after the optimization
according to Scheme II: the skew quadrupole errors are removed in
the two dipoles flanking each horizontally defocusing quadrupole, and
particles are launched along the x axis and the y axis, respectively.
The scales are 2.7- 10‘3 m for the :1: and y axes and 4.0 - 10‘3 for the
a and b axes. ...............................
Calculation order comparison and its influence on the phase portraits
in the (x,a) plane. The scales are 2.7- 10‘3 m for the :1: axis and
4.0- 10‘3 for the a axis. .........................
Measurement results: horizontal position of the center of mass over a
number of turns and its envelope ....................
Behavior of the particles in the normal form coordinates .......
Uniform particle distribution in the sector ...............

xi

61

63

65

67

68

 

4.13
4.14
4.15
4.16
4.17
4.18
4.19
4.20
4.21

4.22

5.1
5.2

5.3

5.4

5.5

5.6

5.7

5.8

5.9

5.10
5.11

B.1

FresnelS function graph .......................... 81

FresnelC function graph ......................... 81
Measurement samples ........................... 82
FresnelS(1.5:c) — FresnelS(:c) function graph .............. 83
Fl’esnelC(1.5:r) — FresnelC(a:) function graph .............. 84
R1 = d — p > 0 .............................. 85
R1 = d — p < 0 .............................. 86
The intersection of two circles ...................... 87
Particle distribution in the beam, the cross indicates the center of mass

of the beam ................................. 90
Calculation results and the comparison with the nonlinear model . . . 92
The grid of points to determine the dynamic aperture ......... 98

The comparison of the phase space portraits before and after the op-
timization using the first approach (minimization of individual high
order coefficients) in the (:r, a.) plane. Particles are launched along the
:1: axis with an increment of one a, a = 82 - 10"6 m. .......... 103
The comparison of the phase space portraits before and after the op—
timization using the second approach (minimization of the resonance
strengths) in the (:12, a) plane. Particles are launched along the :1: axis
with an increment of one a, a = 82 - 10—6 m ............... 105
Phase portraits after the optimization in the (as, a) and (y, b) planes.
Particles are launched along the :r and y axes, respectively, with an
increment of one a, a = 82 - 10‘6 m. .................. 107
The comparison of the phase portraits in the (x, a) plane before and
after the optimization. Particles are launched along the (:r, y) diagonal
with a step of one a. ........................... 108
The comparison of the phase portraits in the (y, b) plane before and
after the optimization. Particles are launched along the (x, y) diagonal

with a step of one a, a = 82 - 10‘6 m ................... 109
The octupole correctors layout ...................... 110
The baseline 50x 50 GeV lattice scheme compared to the 750x 750 GeV

lattice scheme ................................ 113
The final focus section beta functions and dispersion plots ...... 115
The chromaticity correction section beta functions and dispersion plots 116
The arc module beta functions and dispersion plots .......... 117
Tevatron simulation GUI ......................... 132

xii

 

1. INTRODUCTION

The problems of calculation and optimization of the dynamics of charged particles
in accelerators are studied by the field of physics called accelerator physics, or more
generally, beam physics. Beam physics deals with charged particle beams as well as
other beams in various electrophysical devices such as electron microscopes, lasers,
x-ray devices, and computer tomography devices, and also in the most general sense
—— with sets of nearby trajectories in a large phase space.

Despite the fact that the article by Courant and Snyder [21], where the details
of the linear dynamics of particles in the accelerating structures are discussed, was
published less than half a century ago, the theory of control and optimization of the
dynamics of particles in the accelerators nowadays has firm grounds as well as a large
number of classical studies and research methods.

However, this is not surprising, taking into account the pace of the growth of the
requirements and the constraints applied to the designed devices in the sense of the
final energy of the particles, the components manufacturing precision and the strength
of the steering and accelerating electromagnetic fields. For example, the accelerator
whose name will be seen again and again in this dissertation work — the Tevatron
— supplies a center of mass energy of 2 TeV at collision; particles typically stay in
the accelerating channel for 109—1010 turns, while the beam has the size of the order

013

of millimeters, which results in at least 1 of difference in scale. This fact, in turn,

 

results in the necessity of precise orbit calculations as well as accurate estimates of
the particle stability and engineering and alignment error limits.

Mathematical models accompany any accelerating structure during its whole life
span, starting from the instant when the accelerator is merely an idea, and to the last
particle acceleration cycle, or even after that, as in the case of the Tevatron accelera-
tor and one of the possible options for the Muon Collider storage ring discussed later
in this work. At this point, the question of the destiny of the Tevatron project is
unclear. Currently it is scheduled for the decommission by the end of 2009. However,
the question is what could replace the accelerator which was at the edge of scientific
and engineering achievements for many years. There is a tunnel of 6.28 km in cir-
cumference, and it could house a next generation accelerator, thus effectively cutting
down its cost. There are studies conducted toward this direction [49,50].

In this dissertation the lattice for a small Muon Collider storage ring with a center
of mass energy at collision of 100 GeV is discussed. Due to the small final energy of
the particles, the circumference of the ring measures 345 meters. At the same time,
a similar ring with a collision energy of about 1.5—2 TeV could be designed to match
the footprint of the Tevatron tunnel with minor changes. Until now there was no
explicit interest in this project, but recently this situation has changed dramatically.
Various Muon Collider lattices are studied at Fermi National Accelerator Laboratory
(FNAL or Fermilab) and a number of other laboratories around the world. These
studies include models that fit exactly or approximately the geometry of the Tevatron
tunnel to make use of its infrastructure. One variant of the lattice for such a storage
ring was proposed by the author and Dr. Carol Johnstone at the Muon Collabora-
tion Workshop at UCLA in January 2007 [49] and this proposal is included in the

dissertation work.

 

To support the lattice design process with adequate physical and mathematical
models and methods, a number of software packages are under constant development
by various research groups around the world [3,14,18,32]. In general these software
packages can be divided into two large categories, depending on the problems solved.
Some of the programs allow for the dynamics calculation of the individual particles of
the beam, and as such, they require solving of the systems of differential or integro—
differential equations for each different set of beam parameters. This approach is
justified if one needs to have as much information as possible about individual particles
of the beam or for taking into account the interaction of the particles in the beam
(space charge of the beam). On the other hand, one obvious disadvantage of such
an approach would be the need to recalculate the whole trajectory of each individual
particle, which results in an incredible volume of calculations in the case of a circular
accelerator, for which particles stay in the beam tube for millions of turns.

Another group of programs allows to calculate the overall effect of the accelerator
optics on the particles, e.g. in the form of the transfer map connecting the initial and
final coordinates of the particles, and then track an arbitrary distribution of particles
by applying the resulting transfer map to the coordinates of the particles in the
distribution. This gives a natural way to treat repetitive systems by merely applying
the transfer map over and over again for virtually infinite number of times to study
the behavior of the particles for a large number of turns, without spending much
computational effort in recalculating the dynamics of individual particles. COSY
INFINITY (COSY) [13,36], developed at Michigan State University by Dr. Martin
Berz, Dr. Kyoko Makino and others, falls into this category of programs. COSY is
one of the main analysis tools used for the studies presented in this work.

Nowadays it is hard to find a software package that falls exactly into one of the

categories described above. As a rule, programs combine the advantages of various
general approaches to beam dynamics calculation and optimization.

To make sure the calculations in the presented studies are correct, most of the
results have been cross-checked with other popular codes for beam dynamics opti-
mization, such as MAD [3] and OptiM [32] for aspects these codes can handle.

One of the most important criteria for the development of beam dynamics codes
is the running speed: how efficient the computation algorithms are, how fast one can
obtain the results and how much information should be sacrificed for speed consider-
ations, or how long it takes to obtain the results with the desired precision, or even
rigorous results if runtime is not a critical factor.

To increase the efficiency, COSY employs high-performance differential algebra
methods and Taylor model methods for rigorous computations [10—12,34].

In addition to that, COSY supports the transformation from particle optical co-
ordinates to a special set of coordinates called the normal form coordinates, which
could be obtained using a complicated order-by-order nonlinear transformation. The
normal form coordinates have a very important advantage over the particle optical
coordinates: if the transformation can be carried out successfully (general restrictions
for that are not much stronger than the typical restrictions imposed on the behav-
ior of the particles in the accelerator) then the motion in the new coordinates has a
very clean representation allowing to extract more information about the dynamics
of particles, and very convenient for the purposes of visualization. The algorithm of
the normal form transformation [12] is set forth in details in Chapter 3, and in the
subsequent chapters it is used multiple times to assist in solving various problems.

Logically the dissertation is divided into chapters and sections according to the

applied problems and subproblems stated and solved, as follows.

0 Optimization of the Tevatron dynamics:

— Skew quadrupole correction of the dynamics of particles in the Tevatron

in the presence of the systematic skew quadrupole errors in dipoles;

— Calculation of the nonlinear tune shift with amplitude based on the results

of measurements and the linear lattice information;
0 Optimization of the Muon Collider storage ring:

-— Computation and optimization of the dynamic aperture of the Muon Col-

lider 50x50 GeV storage ring using higher order correctors;

— A 750x750 GeV Muon Collider storage ring lattice design matching the

Tevatron footprint.

All the problem formulations include the derivation of the objective functions,
which are later used in the optimization process using various optimization algorithms.
Algorithms used to solve the problems are specific to collider rings, and applicable to
similar problems arising on other machines of the same type.

In linear accelerators, particles go down the accelerator channel only once, so the
most important aspect is to have as much information as possible about the steering
and accelerating fields and their fluctuations, about the particle interaction in the
beam and about the deviations of various parameters of the system. In other words,
one needs as much information about the linear and nonlinear components of the
dynamics, and it needs to be as precise as possible.

To a large extent the above statement is true for circular accelerators too, but
there is a significant difference: the particles in the beam pass the same elements of

the structure multiple times. So the most important task for this case is to study the

‘ ‘9'. VL‘I‘ .1“ “In“

!I-'V'A Li'- H m Z'OL_

details of the long-term behavior of the system and ensure its stability for the desired
number of turns [15]. This includes factors that could disturb the long-term stability.
The algorithm of the normal form transformation is of great value for such problems
as it gives a lot of extra information about the disturbing factors. In addition to the
fact that the dynamics of particles is represented in a way that is easy to understand,
such important characteristics as the strengths of the resonances and the tune shifts
with amplitude and various parameters of the system are calculated.

In Chapter 4, details of the optimization of the Tevatron accelerator parameters
are discussed. The Tevatron is a large and complicated accelerator complex with
lots of features complicating the process of establishing a mathematical model appro-
priately reflecting the state of the machine. Nonetheless, one can pick out various
independent families of the elements of the machine’s optics, and use as a figure of
merit certain characteristics controlled by each family or a group of families. In the
Tevatron at injection, the problem of the skew quadrupole error arises due to the dis-
placement (the so-called cold lift) of the superconducting coils in the main bending
dipole magnets, resulting in a skew quadrupole component increase in the magnetic
field. Two families of skew quadrupole correctors are used in the Tevatron to mini-
mize the cross-coupling between the horizontal and vertical planes. The problem of
finding optimized values for the correctors is stated and solved and the sum of the
absolute values of the transfer map terms responsible for the coupling between the
two planes is used as an objective function.

The skew quadrupole correction problem was studied by different scientists at Fer-
milab and the corresponding results can be found in [26,52,53]. At the same time, the
problem has never been stated before as that of minimization of the functions (4.1.7)

or (4.1.11) and (4.1.13). Also, the differential algebra approach has never been em-

 

ployed to attempt solving the problem. Certain correction schemes were proposed by
the author in collaboration with Dr. Carol Johnstone and Dr. Michael Syphers, and
the best scheme was chosen and implemented during the planned Tevatron shutdown
in August 2004 with minimal modifications due to additional manpower available to
fix more than 50% of the errors in dipoles, which was not the part of the initial prob-
lem statement [48]. Two main correction schemes together with the corresponding
numerical analysis and visual results are presented in Section 4.1 of Chapter 4.

The calculation of the nonlinear tune shift with amplitude based on the results
of measurements and the linear lattice information is discussed in Section 4.2 of
Chapter 4. Finding the nonlinear tune shift depending on the position of the particle
in the beam might not be an easy task, because the nonlinear component of the
dynamics is not known to the desired precision or because there are reasons to doubt
the correspondence between the model and the real machine.

At the same time, there is still a way to find the tune shift, if there is a set of
specific measurements and some extra information which is usually available, namely
that about the geometry of the beam (its size and particle distribution) and the linear
optics effect on the particles (in the form of the one-turn linear transfer matrix).

The method to solve this problem uses the potential of the normal form trans-
formation with its advantages to the full extent. In particular, it is known that in
the case that there are no resonances in the system, the dynamics in the normal
form coordinates is described by a rotation with the amplitude dependent frequency.
Besides, as the tune of the particle can be viewed as the limit of the total phase
advance divided by the number of turns when the number of turns goes to infinity,
the average tune for a large number of turns is the same in both sets of coordinates,

as the contribution of the normal form transformation and the inverse normal form

transformation becomes negligible.

Two different flavors of objective functions are introduced in Section 4.2.1, both
requiring the tracking of the beam center of mass r(N) for a thousand turns. Ex-
pressions for r(N) for various shapes of the beam and distributions of particles in the
beam are obtained in Sections 4.2.3 though 4.2.5.

To confirm that the proposed model for the nonlinear tune shift calculation is valid,
it is compared to both the numerical results for the nonlinear model of the Tevatron
and the independent approximate formula for the tune shift by Meller et al. [40].
Both comparisons lead to similar results: the proposed model shows a discrepancy
of less than 2%, which can be considered a very good result taking into account that
only the information on the one-turn linear transfer map and the geometry of the
beam has been used.

Chapter 5 addresses two problems connected to the Muon Collider storage ring
designs. The studies of the dynamics in muon colliders date back to 1960’s, when
Tikhonin, Budker [17, 54] and later Neuffer [42] proposed ideas of accelerating [1+ and
if particles. Muon accelerators form a distinctive class of accelerating structures
different from both proton and electron accelerators. The advantage that muons
(and electrons) have over protons is that they are truly elementary in the sense of the
Standard Model. Thus, when muons collide, they do not divide up the energy load
as does a particle that is made up of quarks. Muon—antimuon collisions are therefore
clean and the effective collision energy is about 10 times higher than that of the proton
beams with the same energy. Muons also have an advantage over electrons. Since
muons (and protons) are more massive than electrons, they produce less synchrotron
radiation, when in a circular path. There are also drawbacks in using muons: they

are unstable and decay rather quickly (7 = 2.2 psec). This imposes a restriction on

the length of the structure and the intensity of the acceleration. For example, in the
case of the 50 x 50 GeV storage ring, the number of turns for which particles stay in
the channel is not more than 1000, since otherwise the loss of luminosity due to the
decay is too significant. On the other hand, due to the relativistic effects of increase
of lifetime with velocity, this number stays approximately constant independent of
the final energy of muons, so for the 750 x 750 GeV storage ring the number of turns
is again of the order of 1000.

The problem of increasing the volume of the part of the phase space that stays
stable (not lost) after 1000 turns, and as a result, the number of particles being
stored and brought to a collision, represented by the luminosity, is stated and solved
in Section 5.1 of Chapter 5.

Basic principles and building blocks of the Muon Collider storage ring are discussed
in [5,7,27,28,43,55,56]. After the linear lattice and the chromaticity correction section
lattice have been designed to minimize Chromatic and geometric aberrations, higher
order aberrations must be optimized in order to enhance the system performance. The
interaction region is the primary target region for this problem. In general, the more
nonlinearities progress to the interaction point, the smaller the dynamic aperture and
the momentum acceptance are, which means smaller beam sizes, and hence, lower
luminosity. Various methods for solving the problem of this kind for the 50 x 50 GeV
storage ring are discussed. The baseline lattice is described in [28, 55]. Using the
existing model and deliberately adding certain high order multipole components to
the superconducting final focus magnets to compensate for the nonlinearities of the
lattice, one can achieve larger dynamic apertures, and hence better stability and
higher luminosity.

The main problem with the task of improving the dynamic aperture of the ac-

celerator channel is a large number of nonlinearities to control, while the number
of correctors is limited, and in general one wants to keep this number as small as
possible. One very effective approach for achieving the goal is to use the dynamic
aperture itself as a figure of merit and try to maximize it using optimization meth-
ods. The corresponding objective function is derived in Section 5.1.4, and it is given
by Eq. (5.1.3). Other approaches, such as the minimization of the most important
aberrations of higher orders or the minimization of the resonance strengths, show
less pronounced results. For each of the proposed approaches, the objective functions
are introduced in the form of Eqs. (5.1.1) and (5.1.2), correspondingly. The normal
form transformation plays an important role in two of the three approaches’ imple-
mentations, namely, the resonance strength minimization and the dynamic aperture
maximization. This transformation allows for the calculation of the values on the right
hand side of the equations and for representing the results of different approaches in
a uniform way (in one common coordinate system).

The second part of Chapter 5 addresses the issues with a totally different machine,
yet intimately connected to the 50x50 GeV Muon Collider storage ring lattice. It is
the problem of the design and optimization of the 750x750 GeV storage ring that
matches the footprint of the Tevatron accelerator. This lattice has all the same
building blocks as its 50 x50 GeV counterpart, and thus shares most of the advantages
of its clear and highly optimized design.

Each major section is supplied with the results of numerical experiments of ap—
plying various optimization methods to the problems stated. The emphasis is made
on the efficiency comparison of various approaches and methods.

Various optimization methods, such as the Nelder—Mead simplex method [41],

modified Newton method (Levenberg—Marquardt method [33,37]), and the simulated

10

annealing method [31] are used in the computations.

11

2. EQUATIONS OF MOTION, FIELDS

AND POTENTIALS

2.1 Equations Of Motion of the Particles in

Accelerators

2.1.1 Curvilinear Coordinates

Consider the motion of the particles in the electromagnetic field [8]:

d1?

dt

q(E'+27x 3'). (2.1.1)

Usually when studying the dynamics of particles, the time t plays the role of
the independent variable, and the positions 52’ and velocities 17' or momenta 15' serve
as phase space coordinates. In beam physics it is more convenient to use a system
coordinates relative to some reference particle. Besides, instead of the independent
variable t the arc length .9 along the reference orbit is used. In the new coordinates,

called particle optical coordinates, each point is determined by the state vector

( 1' l

y
2': I: W— t”) , (2.1.2)
a =pcc/P0

 

b = Py/PO
K6 = (E— Erma)

 

12

where .’L‘, y are the horizontal and vertical positions of the current particle with respect
to the reference particle, l is a space-like variable characterizing the deviation of the
current particle from the reference particle in the time Of flight (t — to) multiplied by
a constant k of the dimension of the velocity, a and b are the horizontal and vertical
reduced momenta, 6 is a relative difference between the total energy of the current
particle (E) and reference particle (E0), and 190 is the momentum of the reference
particle at the beginning.

By this definition of 2', the reference particle corresponds to Z = 0, hence the

vector 2' does indeed describe relative motion.

2.1.2 Equations of Motion in Curvilinear Coordinates

Consider the case when the reference orbit is allowed to bend in only one plane. Define
the momentary curvature of the reference orbit h(s) as a function of the arc length 3.
If the curvature is nonzero, the radius of curvature is then given by R(s) = 1/h(s).
Consider the bend angle that the reference orbit experiences as we move from position

so to position 3:

a = /a da = s: Rdé) = f: h(s)ds. (2.1.3)

0‘0 0

The equations of the dynamics in the Cartesian coordinates have the form (2.1.1),

or in the integral form
13(3) _.
as) = aso) + / mm
“30) (2.1.4)

3
= 13(30) + / F(§)t’d.§,
30
where t' = dt/ds, and F(s) depends on if and 15'.
The orientation of the locally attached particle optical coordinate system changes

from so to E, as it was rotated by the angle a from Eq. (2.1.3). So in the new local

13

coordinates
cosfsgohds 0 sinfsihds

§
mg) = 0_ 1 g - (13130) + / F(s)t’ds) . (2.1.5)
— sin fsso hds 0 cos [3% hds 80

Let M(so, 5) denote the matrix on the right hand side of Eq. (2.1.5). Differentiating
with respect to 5 and evaluating at so = s, 5 = .3 allows to conveniently obtain the

rate of change of the momentum p‘j:

153(5) = M(3,§)F(§)t’+M'(s,§)(fi(s)+/S§F(§)t’d§)]

0 0 h(s) =
= F’(s)t’ + o 0 0 5(3), (2.1.6)
——h(s) 0 0

where the matrix M is unity and the integral on the right hand side is zero under the
condition 5 = 3.

Hence, the first resulting term depends on the actual forces, and the factor t'
accounts for the fact that there was a change to s as an independent variable. The
second term is a pseudo force due to the fact that we are located in a rotating frame. If
there was an out—of-plane motion of the reference orbit allowed, the matrix M would
depend on two curvatures. Unfortunately, in this case an additional complication
arises from the fact that rotations around different axes do not generally commute.

Now consider the particle having a position a: yé 0. In this case

dl _a:+R

 

 

21-8- — R = 1+ has, (2.1.7)
dz _ dz dl _ da: _ dx/ds _ Pa:
ds — dl ds — dl<1+hx)_ dl/ds(1+hx)_(1+hx)p_s° (218)
Similarly for y one obtains
dy _ py
ds — (1 + hr)ps, (2.1.9)

14

and for the time of flight using the above expressions,

— = M32116- 9+ (ff-i)

2
= l<1+hx)\/px+py+ps =
1) p3

 

 

 

(1 + h$)p—ps, (2.1.10)

 

 

where p = \/pg, + p5 + p3 and v = \/ "0% + 123 + ”0:;2 are the magnitude of momentum
and velocity of the particle under consideration.

Therefore, the equations of motion in local coordinates with s as the independent
variable were obtained. The next step is to change to particle optical coordinates
(2.1.2).

To take into account relativistic effects, it is convenient to introduce the relativistic

measure
E - 63Wm. 31. 8)
mc2

 

, (2.1.11)

characterizing the ratio of kinetic energy to rest mass energy. It follows that

1
1—

<3qu

E: ’ymc2, E= mc2 +K,=K (—1)mc2,n=_=__,

1 2 277+772
_:\/ _(_1_+_n) = l/W’ (2.1.13)

£2..— T’l’flz 31: ,/n(2+n). (2.1.14)

me me C

and

 

Express the rate of change Of the particle optical variable 1 in terms Of particle

15

Optical quantities:

 

dz k p
, : —= I—t’ =—]_ h ——kt,
I d3 Mt 0) v( + @193 0
k
= u+ma5@£——-
POPS?) ’UO
kpo
= 1+hr——1+ m——
( )p0p3( )
1 k:
= [(1+hr) +U@_1]_, (2.1.15)
1'I'77019s ”0

where the following relations have been used:

1 1
v s 0 0

p p0
- = (1 +7ilm, — = (1 + Dolme-
v '00

The next step is to obtain pO/ps, which is on the right hand side Of the equation

 

 

 

(2.1.15):
2 s
12 p0 P_2__a2_b2 _ 77(2+71)m2 _ 2_b2
7’3 (/p2 - p2 - p2 1’12) "0(2 + "Olmb

where the expression for the momentum (2.1.14) is used.
Because it is always true that 17 [I if, it follows that 17/11 = 5/ p. This fact is useful

for the derivation of the equations for the reduced momenta:

1(a32e)=
d8 po’po’po

0 0 h(s) _,
_ 1 “ I P _
0 —h(s) 0 0 0
-.t’ “ ..
= zeE— + zeg x B(1+ hx)£—i + h (£0, Jim) =
190 v mm 190 290
E z' r B’
= —— (1+ 39) + 3 x ——(1 + MOE + h (50, JE), (2.1.17)
XeO k 130 XmO ps po :90

16

where the following notation has been introduced:

_m _mm
XmO "‘ —9 X60 _ °
26 28

The first term is transformed according to the following relation:

l’ 1 [’21

I I I I I O
—kt— -t +___ +—

I ( t0)”: 0 k ”0(1 k),

and in the second term the quantity p is reduced, while the numerator and the
denominator have an extra factor of po.

I
The expression (1 + I739) can be transformed using Eq. (2.1.15), so ultimately

1(22 a a) =
Po’Po’Po

d3
= _E_(1 + h,)_1+_na
xeo 1 + 770133
+ (st — BEBy, '12in -" aBs, aBy — be) (1+ h$)'—1—p'_0 +
p0 P0 XmO P3
+ “gm—135). (2.1.18)
190 190

The equation for the full energy deviation has a very simple form 6’ = 0, as long as
there is a conservation of energy. In RF cavities the equation has a more complicated
structure, but that is not the issue for the storage rings, where the energy of the
particles is constant.

Hence, the equations of motion in particle optical coordinates are

 

 

 

 

 

x’ = a(1+ hr)p—0-, (2.1.19)
P3
y’ = b(1+ moi-”2, (2.1.20)
173
1 k
l’ = [(1 + hm) + 77 E — 1] —, (2.1.21)
1+770 Ps '00
1 E B
a’ = (1+hz)[ +nmi+st p—°——y—]+h5’£, (2.1.22)
1 + 770193 Xe0 XmO Ps Xmo P0
1 E B B
b’ = (1+hx)] +nIfl—l—a 329+—$], (2.1.23)
1+ 770193 X60 XmO p3 XmO
6’ = 0, (2.1.24)

17

where n is determined by Eq. (2.1.11), and £131 — by Eq. (2.1.16).

2.2 Fields and Potentials

A large part of the problems of electrodynamics is described by the set of laws ex-
pressed by the system of Maxwell’s equations. In the international system of units

(SI), these equations have the following form:

(

 

‘7 - D = p (Coulomb’s law),
\7 - E = 0 (absence of magnetic monopoles),
{ .. _, BE _. (2.2.1)
V x E + I}? = 0 (Faraday’s law),
1 V x If -— 39-?- = j(Ampere’s law),

where D is the vector Of the electric flux density, p is the charge density, E is the
magnetic flux density, E is the vector of the electric field, If is the vector of the
magnetic field, and f is the electric current density.

The system of equations (2.2.1) is not closed, three more equations need to be
established: the so-called constitutive relations connecting D with E and If with E,
and one more relation connecting j with E. The first two of them have the following

general form:

D = 50E + I3,
(2.2.2)
B = MOH + M,
where 13, M are the vectors of electric and magnetic polarization of the medium;
60, [1.0 are the electric permittivity constant (dielectric constant of vacuum) and the
magnetic permeability constant of vacuum, respectively.

In the case of a homogeneous and isotropic medium, the equations (2.2.2) are

simplified, and besides, there is often one more relation, namely, Ohm’s law for the

18

linear medium:

D: d?
f: d?

Here 13,112 are proportional to the field applied, and hence, their influence can be
included in the form of new constants. Therefore, 5, [1. characterize the medium and
are referred to as the relative electric permittivity (dielectric constant) and the rel-
ative magnetic permeability, respectively. The constant a characterizes the electric
conductivity of the matter.

In many cases, it is of interest to study Maxwell’s equations in the absence of

matter at the point of interest. In this case
p=a J=6 @2Q

Thus, Maxwell’s equations take the simplified form:

625

||
5:

<1
m1
II

9

(2.2.5)

<1
X
PM
II
I
Q3
Q1

,__—5._
<1
X
m1
II
Q3] ‘39
FI-

and in the stationary case

<I1
U11
II

P

<1
Q31
ll

go

azm

<1
X
U11
ll
5:

<1
x
U31
II
o

19

2.2.1 Fields and Potentials of the Magnetic Multipole

(Stationary Case)

As a rule, the guiding magnets of the accelerating structure create a field which
changes slowly with time or does not change at all, and the beam of particles being
controlled by the magnets is inside a vacuum beam pipe. Therefore, one can use
Eqs. (2.2.6). There is one limitation however, in that this model does not take into
account the interaction of the particles. Nonetheless, in most cases this is not a
significant limitation, and if necessary, the particle interaction and beam—beam effects
can be considered separately.

In most cases only the magnetic component of the electromagnetic field is used
for focusing purposes, as the particles are running in the beam tube with the velocity
compared to that of light, and for this ultra-relativistic case a moderate magnetic
field of 1 T corresponds to a huge electric field of 3 - 108 V/ m, which is very hard to
create and control.

In various sources the multipole expansion for the transversal guide field described
below is treated differently [12,57]: expansion of the magnetic field E, expansion of
the magnetic vector potential A, expansion of the scalar potential V, using cylindrical
or Cartesian coordinates. The most appropriate approach is chosen for the specific
problem, but the physical meaning is always the same.

Consider the transversal dynamics of the particle. Denote z = :1: + z'y, where a:
is the horizontal displacement of the particle from the magnet axis, y is the vertical
displacement (the “verticality” and “horizontality” are relative to the plane of the

orbit of the reference particle). Assuming that the edge field contribution is negligibly

20

small, the second equation of (2.2.6) yields:
divB=v-B=0 => B=V7xxf=rot2fl (2.2.7)

where ff is the vector potential of the magnetic field. As the magnetic field is transver-
sal, the vector potential has only one component As in the direction of the motion of

the particles. The fourth equation of the system (2.2.6) yields that

d

rotE = 0 => B = —VV, (2.2.8)

where V is the scalar potential of the stationary transversal magnetic field.

Combining Eqs. (2.2.7) and (2.2.8), one obtains

 

 

6V 8A3 6V 8A3
Ba; — —E -— 8y and By — —5:'y- — — 811'; . (2.29)
Define a complex potential as a function of z = a: + 23/ in the following way:
A(z) = As(r, y) + z'V(:r,y) = As(z) + iV(z). (2.2.10)

Then Eqs. (2.2.9) are the Cauchy—Riemann conditions for 141(2). Therefore, the com-

plex potential is an analytic function and can be expanded in a power series
24(2) = Z ann; Kn = An + inn; (2.2.11)

where An, an are real constants.

The expansion (2.2.11) converges inside a circle [2| < Tc, where TC is the minimal
distance between the origin of the expansion and the iron yoke or the superconducting
coil or any other field source causing the non-homogeneity to appear on the right hand

side of Eq. (2.2.8). As a result, Eq. (2.2.9) does not hold anymore.

21

 

 

 

 

)
X
Figure 2.1. Axes orientation for the cylindrical coordinates

2.2.2 Multipole Expansion in Cylindrical Coordinates

In many cases it is convenient to express the field in cylindrical coordinates ('r, (0, s),

such that
{ = ”.0390 . (2.2.12)
y = rsm (p
In this case 2" = rne’mP = r"(cos 7190 + z'sin mp). Eqs. (2.2.10), (2.2.11) yield that
+oo
V(r, (,0) = :01.” cos mp + An sin mph”, (2.2.13)
11:0
+oo
A3(r, (p) = :01” cos nap — #n sin n1p)r".
n=0

Using Eqs. (2.2.8) and (2.2.13), one can get the expansion of the magnetic field

in cylindrical coordinates (see Fig. 2.1).

+00
1 8V
BSD = —;5<; = — n()\n cos ngo — an sin n<p)r"_1,
4:1 (2.2.14)
8V . _
Br 2 —5; = — "2:1an cos me + An smncp)r" 1.

22

Let us introduce the following notation: To is the reference radius for the multipole

expansion, usually chosen to be equal to the maximum deviation of the particles in

the beam from the reference orbit; Bmam is the main field component of the magnet

under consideration; bn is the so—called “normal” multipole coefficient, and an is a

“skew” multipole coefficient, defined as follows:

 

 

A
bn = — n n 7"ng and an = "Mn rn—l. (2.2.15)
Bmaz'n Bmain
The meaning Of the notation becomes clear from Eq. (2.2.20) below. Using this
notation,
+00 n—l
, r
B(p(r, (p) = Bmam :(bn cos up + an sm mp) (—) , (2.2.16)
n=1 TO
+00 n—I
, r
Br(r, (p) = Bmain Z(—an cos mp + bn sm mo) <—) , (2.2.17)
n=l 7‘0
+00 an bn 7' n
V , = —B - —— — ' — , 2.2.18
(r so) mama-0g n cos-212+ n sum) (0) ( >
+00 b a r "
As('r, (p) = —Bma,-n'ro 2(1 cos mp + il- sin mp) (—) , (2.2.19)
":1 n n m

where each value of 71 corresponds to the 2n—pole component of the field. For an ideal

2n-pole magnet, bn = 1, while all other coefficients vanish.

Consider
+00 1‘ 71-]
13¢ + 2'31. : Bmm-n Z (7) (bn(cos mp + i sin mp) -— z'an(cos mp + 2' sin mp)),
0
n=1
+00 7‘ n—l _
Btp + 131' = Bmain Z (E) “’71 - zIan)€m’p.

11:1

This representation allows to find the strength of the 2n-pole component:

(IgIln = (VB-12+ B920)" = Bmain ({6)n—1 V122,, 'I' 0%.

23

(2.2.20)

Eq. (2.2.20) yields that the magnitude of the 2n—pole component does not depend
on the azimuthal angle and scales with the (n — 1)-th power of 1'; an and bn are the
relative field contribution of the nth multipole component to the main field at the
reference radius r = To. This is the reason for the normalization given by Eq. (2.2.15).

Applications of the cylindrical coordinate representation include:

o Superconducting magnets design, where it is convenient to find the arrangement

of the current coils;

o Multipole component measurement with a coil rotating in the field, where the

induced voltage is proportional to \/ 32,, + (1%, while its phase is related to an/bn;

0 Problems where it is easy to see which multipole components are forbidden, if

specific symmetry properties of the field are assumed.

2.2.3 Multipole Expansion in Cartesian Coordinates

In Cartesian coordinates the expansion (2.2.11) assumes the form

~ +00 +00
A(z) = Z Isnz" = 201,, + ian)(:1: + 1y)". (2.2.21)
n=0 —0

Separation of real and imaginary parts of Eq. (2.2.21), yields the following using

24

the notation introduced in Eq. (2.2.15):

+oo
b a2
A3631?!) = Re 2 ann = —Bmain{blx 'I' aly + $1032 — :92) "I“ 71153] 'I'
1120 0 0
b a3
+ —%—(r3 — 3ry2) + ——2(3:c2y -— 313) +
3r0 3r0
b a
+ -i3(:1:4 6122312 + 314) + —§-(a:3y — ry3) + . . .}, (2.2.22)
41‘0 r0
+00 11 a2 2 2 b2
V($ay) : 1m 2: ”n2 = Bmain{a1$ _ bly ‘I' 5770.“: — y )_ airy 'I’
n=0
a3 b3
+ —§($3 3231/2) —2(3$2y - y3) +
31'0 r0
b
+ 3330174 — 6:1:2y2 + y4) - —%(z3y — 2311/3) + . . H} (2.2.23)
47‘0 r0

Similar to the case of cylindrical coordinates, the components of the magnetic field

can be obtained by virtue of Eq. (2.2.8):

6V b2 0-2 (13 b3
B$(x’y) = "'5— : Bmain{_a1 ‘I' —y _ —$ - 7(32 _ 312) ‘I" 721‘?! —
CD T0 To To To
+ iii-(933 31392) + :(32321/ - y3) + ] (2 2 24)
0 0
8V 0-2 b2 03 b3
By($1y) = —‘a—' : Bmain{b1+ —y —' —$ + 72133] + 7(332 — y2) +
11 1‘0 7‘0 To T0
(3
+ 2§(3$2y - 313) + —§(:c3 - 323/2) + . . H} (2.225)

Applications of the Cartesian coordinate representation:

0 If the dynamics of particles is described in Cartesian coordinates, it is easy to
identify the contribution of each individual multipole to the equation of motion.
For example, if all the “skew” coefficients an vanish, a flat beam (y E 0) remains
flat, which means there is no coupling Of the horizontal motion into the vertical

motion;

0 Each term in the representation (2.2.23) describes the shape of the magnet pole,

creating the field, as this contour is an equipotential surface.

25

Note that in the case of the electrostatic focusing elements, the derivation of the
formulas for fields and potentials will follow the above scheme due to the equivalence
of the first equation to the second and also that of the third equation to the fourth
equation of the system (2.2.6). On the other hand, as it was mentioned in the begin-
ning of the section, in most cases using the magnetic fields for the focusing purposes
is much more convenient.

Now, as the representations for the scalar magnetic potential V and the com-
ponents of the magnetic field E have been found in both cylindrical and Cartesian
coordinates, one can return to the equations of motion of the particles (2.1.19—2. 1.24),
and claim that all the quantities on the right hand side are known as functions of

fields and potentials: 17 = 77(V), E = (Br, By, B3), E = (Ex, Ey, E3).

26

3. NORMAL FORM ALGORITHMS

Normal form algorithms are used throughout this dissertation work, and the problems
considered in Chapters 4 and 5 make use of the algorithm either in the main opti-
mization step, or for better representation and visualization of the results obtained.
The main ideas of the nonlinear normal form transformation algorithm are presented
in [9,12]. The idea of this chapter is to introduce the notation used in the subsequent
chapters and place more emphasis on the steps of the algorithm which are not in the
book and might not be Obvious. All the systems considered in the dissertation are
symplectic, so the emphasis is made on the normal form transformation algorithm for

the systems satisfying the symplecticity condition.

3.1 Nonlinear Normal Form T1'ansformation

The idea of the transformation is to find such a nonlinear change of variables that
the motion in the new variables, up to a certain order, is circular with an amplitude-
dependent frequency.

Consider the nonlinear transfer map [8, 12,47] of a particle optical system

—0

5f =M(z",-,6), (3.1.1)

where Z" is the vector of 21) phase space coordinates, (f is the vector of parameters of

the system, 2'; is the initial state of the system, and if is its final state. Thus, M is

27

 

the function connecting the initial and final states of the system for some interval of
time.
The transformation to the normal form coordinates consists of a sequence Of non-

linear coordinate transformations of the form
IV: AoM oA‘l, (3.1.2)

where 0 denotes a composition of transfer maps.

Assume that the transformation matrix satisfies the following relation:
M(5, ") = (6, "). (3.1.3)

If this does not hold, one can find the fixed point and transform the map to this fixed
point, and therefore the new transfer map satisfies Eq. (3.1.3). To be able to do this
it should be assumed that none of the eigenvalues of the linear part of the map is 1,
which is always the case for stable multi-turn systems. Also, if the map has certain
parameters (the vector (f is not empty), M — I H must be invertible, where I H is the
block matrix with a square unity submatrix 21) x 21) in the upper-left corner and all
the other elements are zeros.

The first step of the transformation to the normal form coordinates is the diago-
nalization of the linear part of the map. This can only be done with the assumption
that all the 21) eigenvalues are different. For most of the circular accelerators this
does not represent a limitation.

All the eigenvalues can be written as complex conjugate pairs 736in , as in [12],
where for the stable symplectic systems Tj = 1 and 11,- is real Vj = W. In the

pi

basis of the complex conjugate vectors J corresponding to the complex conjugate

28

eigenvalues the linear part of the transfer map has the following form:

 

 

{ newl 0_ 0 0 \
0 rle—’“1 0 0
R: E E E 5 . (3.1.4)
0 0 rvefl‘v 0
K 0 0 0 rue-1““)

After the diagonal form of the linear transfer map is found, one proceeds to a
sequence of order-by-order transformations, each of which affects only one particular
order of nonlinearities. The aim of the transformation is to simplify the nonlinear part
of the transfer map as much as possible, ideally removing all the nonlinear elements
up to a certain order.

At the mth step one divides the transfer map into the linear part ’R. and the part
8m, containing all the nonlinearities. If the previous step was successful, the map Sm
contains only the terms of order m and higher.

The desired transformation has the form
Am = I + Tm, (3.1.5)

where Tm vanishes to order m — 1. As the linear part of Am is invertible, the map

Am itself is invertible [12], and the inverse to order m is
AT—nl =m I - Tm. (3.1.6)

After the transformation, one has up to order m
Am 0 M 0 A1711 =m
=m(1'+ Tm) o (R+ Sm) o (I — Tm) =m
=m(z+ T...) o (71 — 72 o rm + .9...) =m (3.1-7)
=mR—Ron+Sm+TmoR=m

29

 

The goal is to simplify the terms of order m of M, or in other words, it is necessary
to find Tm such that Sm = —-Cm, where Cm = (Tm o ’R — ’R. o Tm). If the range of
Cm is the full space, Sm can be removed entirely. However, there are circumstances
limiting the range of Cm.

Consider the representation of the unknown Tm in the form of the Taylor poly-
nomial of degree m in the eigenvector basis. Hereafter, the pairs 53’: Vj = f1? denote

the coordinates in this basis.
—. —._ [6+ ___ It‘— +
TrrIj=Z(TrrIjIk+,k)‘(Sf)1(31)1'°I-I(51I-)k (s. +v>’° (318)
- _ - ~_ 3+ _ r + ' '
T..,-=Z<7m,|k+,k >-<sr)1(s1)1 ...<s +>v (sink-

where (T,%|E+, It") are the coefficients of the corresponding orders in the expansion,
13+ = (3+,...,k;r), and E- = (k;,...,k;;).

Derive the equation expressing the coefficient Tm in terms of the known compo-
nents of Sm. It is true that Cm = (Tm o R — ’R o Tm), and also the form of ’R is given

by Eq. (3.1.4). Besides that, the goal is to achieve ij = —Smj. Hence,

—(s,$,j|E+,E-) = of m.|k‘+,k “-)= (3.1.9)

||
A

(7;, 071— 120%,}; |k+,E‘)=
21)
11.1, 1,)( war):
I:
v ++kf -- "+_"— :l:' - —- -_
(Hrlkl I) 6211“? k )"Tj'e 2H3) (T:j]k+,k ),

where A2j_1 = 738+le J , Agj = rje ” J and j = ,v are the eigenvalues of the linear

(n.
=<

part of M, or ultimately,

Si -0+ 4—
_. _, - .k ,k
(Ti.|k+,k‘) = ( Sm’I ) . (3.1.10)

k+ +kl_ ...-' ~_ - .
((Hl— 17.11 I) , ezp(k+—k ) _ 7'j . 8:13.113)

30

 

The algorithm efficiency, therefore, depends on the conditions for which the de-
nominator in the equation (3.1.10) is zero. In this case it is not possible to eliminate
all the elements of 8m. The next section addresses such conditions for the specific

case of rj and W-

3.2 Stable Symplectic Case, Nonlinear Tune Shift

with Amplitude

In the case of the stable symplectic system, all the rj’s are equal to 1, all the nj’s are

purely real, and the condition for (3.1.10) to become zero can be written as

17- (13+ — E—) = ipj mod 27r, (3.2.1)
where the sign on the right hand side corresponds to the member of the pair 173:.

Assume the equation [1' - fi = 0 mod 27r does not have any nontrivial solutions. In

this case, the only way Eq. (3.2.1) can hold is
if = 1.,— v1 74 j, k; = k].— 3: 1. (3.2.2)

Proposition 3.2.1: In the basis comprised of the complex conjugate vectors,

mj “3 the polynomials of Sis} - . . ,sj‘,s;‘ . The same holds true for the

pairs 8m -, Smj'

Proof. First of all, prove that 81:3” are the complex conjugate functions Of
3331‘, . . . , 53, s; . As after the diagonalization all the variables are split into in-
dependent pairs, it is sufficient to show that the proposition holds true for one such

pair.

According to Eq. (3.2.3), in the case of one pair of variables
21—1: (at-17*), A = (27+,r)’1, (3.24)

where the indices of 17+ and 1')” are left out. As 5— = v

_1_ a'I'Zb a—Zb — _‘_]_'_'— .
A _(c+id C_id)._(-u ,v ), (3.25)
i c—id —a+ib
=———— . 3.2.6
A 2(ad—bc)(—c—id a+ib) ( )
AOMOA—1=R+A0NOA—1, (3.2.7)

where N contains all the nonlinear components of the map M. Consider

N O A—1 ( 3+ ) : (N1((a + ib)s+ + (a — ib)s“; (c + id)s+ + (c — id)s") ) ‘
s“ N2((a + ib)s+ + (a — ib)s—; (c + id)s+ + (c - id)s—)
(3.2.8)

(a + ib)s+ + (a — ib)s_ = a(s+ + s‘) + ib(s+ — 3‘),
(3.2.9)
(0 + id)s+ + (c — id)s_ = c(s+ + s‘) + id(s+ — 3").

Before the coordinate transformation, all the coordinates of the vectors under
consideration were real. The matrix A—1 gives the transformation back to the old
coordinates. Hence the left hand side of Eq. (3.2.9) is purely real, which means
5+ + s— is purely real, while 3+ —— s" is purely imaginary.

As N1 and N2 are the polynomials with real coeflicients (Taylor polynomials of

the original map M, of which all the elements are real), and hence the values of these

polynomials after 3+ and s‘ are substituted, will be real:
N1((a + ib)s‘L + (a — ib)s"; (c + id)s+ + (c — id)s_) =
./\/‘1(a(s+ + 8—) + ib(s+ — sT);c(s+ + 8—) + id(s+ —— s—)) = 711 E R,
(3.2.10)
N2((a + ib)s+ + (a — ib)s_; (c + id)s+ + (c — id)s“) =

N2(a(s+ + s’) + ib(s+ — s’);c(s+ + s”) + id(s+ — s_)) = 172 E R.

32

Denote
1

____2(ad _ be) = ,,, (3.2.11)

then the ultimate form of the transformation to the new coordinates is

o o —1 3+ _ )1; i(c—id) —z'(a—ib) 771 _
A N A (3..) " (—i(c+id) i(a+ibl )(772)-

K.( (d+ic)n1 - (b+ia)n2 ) =
(d - M711 - (b - ia)‘I72

R (61771 - b1I2) + 2'(0771 - 0-772)
( (61771 - b172) - z'(6'01 - a772) ) ' (32.12)

All the expressions in the parentheses are real, and hence 872% are a pair of complex
conjugate polynomials of 3+, 3“.
To show that 77:]. = ng, it suffices to consider the denominators of (3.1.10) and

assure that they are complex conjugate as well.

In the symplectic case all rj = 1, and hence the expression for the denominator

reduces to
eifi. (12+ — IE") — ei‘W. (3.2.13)
By virtue of
hf” = 1:,- W75 j
, (3.2.14)
Q=gil

one can Obtain

13+ — is" = (0,...,0,:l:1,0,...,0),
with $1 at the position j. Hence,

eifi-(k+—k_) = eiiflj.

It is necessary to show that

 

aim-7 — 6+le = emu-7' — 6%“. (3.2.15)

33

Depending on the sign, there are two cases:

EB:

 

6W7 ._ 6W7 = a = 0 = 8—112 _ 6"”7';

 

eflflj — e+mJ = e+wJ — e_wJ.

In both cases Eq. (3.2.15) holds, which completes the proof of the proposition. I

According to the proposition above, in the symplectic stable case the elements of

the transfer map in the new coordinates si", 31-, . . . , 33’ , 5,? corresponding to the pairs
of the complex conjugate eigenvectors 171i . . .173: have the form
M“I = S—I— - f-(s+s— s+s—)
.7 1 1 - 1 v v
" J (3.2.16)
M]— = s; - fj(s'lfsl—,...,s:s,7)

The fact that s; in the first line and 3].— in the second line can be taken outside

the parentheses, and that f depends only on the products of the form 33's}: follows
from the condition (3.2.2):
+.
sj .
k+ = k'— +1 (5 l J (S—lk; = (8+)k;8+(8_) 7- : 3+(s+s_)kJ—-
.7 .7 ’ .7 .7 .7 .7 .7 .7 .7 .7 ’
s—
k+ k— k‘I‘ k'I' k‘I'
+ _ — _ + — _ + - — _ +
k —kj 1,(sJ)J(sJ-)J (3])J(sj)Jsj—s(sjsJ)J.

The conjugation in the second line of Eq. (3.2.16) follows from Proposition 3.2.1 and

the fact that 33-" = s].— as the coordinates in the complex conjugate vector basis.

Indeed, if the basis vectors are complex conjugate, the rows of the matrix A will be

'I'
.7

Eq. (3.2.12) the result in Eq. (3.2.16) follows for arbitrary 771, 772 6 IR.

complex conjugate as well, and hence s and s; are complex conjugate. Thus, from

34

 

It is not convenient to work with the pairs sji for the purpose of obtaining the
ultimate equation for M%, therefore the resulting map is expressed in terms of the

pairs 15%, such that

J J J . (3.2.17)
t].— = (s; — sz)/2z
Then
Sf :13; + th , (3.2.18)
3]- — t; — it]
33+ 3].— = (tj)2 + (1;)2,

A=1(1, 1); A—1=(1 ’,). (3.2.19)
2 —z i 1 —z
. :t .

.7 .7 I

+ + — + —
:t _ s--fj(slsl,...,s,,s,,) _
M7 _ (33-7(sfsl—,...,s;]'s;) ) —
1 1 1
= 2 ( —i i ) '
((tj+it;>f,-(<tr)2 +(t;)2 ,. .2332) +t..< -=>2))
(t; —z't;>7<(tt>2 +<t; )2 .(t+>2 +(t )2)
= _1_ 1 1 )aj ( (tj+itj— )(cos(<pj)+isin((pj)) ) =
2 —i i (ti—itj )(cos(cpj)—isin(<pj))

z ( 222-; -1...) ) (t2- ), .3...)

cos(<pj) t;

- . + 2 — 2 + 2 — 2
where fj = aj -eupJ((t1 ) +(t1 ) "”0” ) Ht”) ) and aj = const, as a symplectic mo—
tion is considered, and (p, = cpj((tif)2 + (t1—)2,...,(t.ff)2 + (5)2) depends on the
rotationally invariant quantity.
Therefore, in the new coordinates the dynamics of particles is described by a

rotation, the frequency of which depends only on the amplitudes (th-r)2 + (tj—)2 and

35

possibly the system parameters. The functions (p, for j = I717 give the tunes Of the
nonlinear motion. Indeed, for any repetitive system the tune of the particle is the
total polar angle advance divided by the number of turns in the limit of turn number
going to infinity if this limit exists. Consider an arbitrary particle of the beam. The
first change of the phase is due to the transformation to the new coordinates, and
after that each turn produces an equal polar angle 1,0]- which depends on the amplitude
and parameters of the particle. To obtain the final angle, it is necessary to perform
a transformation back to the original coordinates.

As the number of turns increases, the contribution of the initial and final po-
lar angles due to the transformation becomes increasingly insignificant, and in the
limit the tune results in only (pj. Hence, the limit exists and it can be computed
analytically as a by-product of the normal form transformation. This fact is a key
point to understanding the algorithm of the nonlinear tune shift calculation using the

measurement results, discussed in Section 4.2.

3.3 Motivation for the Differential Algebra

Approach

Consider one of the simplest possible nonlinear transfer maps, the sextupole kick
(thin lens approximation) superimposed on an idealized linear dynamics in a one-
dimensional phase space (7:, :r’), where :1: is the position of the particle, and 73’ is its
velocity [16, 58]:
( 7:1 > __ ( cos(27rn):1:o + sin(27r/I):c6
as’l — — sin(27r/I)a:o + cos(27rp):r[) + ks [cos(27rn):c0 + sin(27rp)a:()]2 ) ’
(3.3.1)

where k3 is the parameter, characterizing the sextupole strength.

36

 

Introduce the following change Of variables:

5: =2 :5 + 2411582 + 24123333, + 14220772
, (3.3.2)

515’ =2 73' + 311732 + 312737;, + 82203,)?
where the coefficients A11,A12,A22,Bn, 312, 322 are to be defined in such a way
that all the second order aberrations vanish. This is equivalent to the first step of
the nonlinear normal form transformation (3.1.5), described in Section 3.1. The “:2”
means that the partial derivatives of orders up to 2 of the right hand side and the
left hand side match at the origin.

As up to the second order it is true that

5:2 ___2 $2

535;, :2 $53, (3.3.3)

(all the other terms are of higher orders, and hence, omitted), it follows that

IE =2 53 - A1112 - 24125553, - 24220532
. (3.3.4)
as, :2 2’ — 311522 — Blzix’ — B2262?)2

Substituting these expressions into Eq. (3.3.1) yields in the new coordinates:

~ ~ ~ ~I ~I
( 9:1 — A1133? - 14123311131 — 14220131)2 ) =2

~I ~ ~ ~I ~I
$1 — 311$? - 312131171 - B22031)2

= ( COSQWXiO - A1155?) - 2412530536 - A22(536)2) )

- sin(27ru)(530 - «4115(2) - 2412520536 — A22If6)2) (3 3 5)
( + sin(27rn)(.i:6 — B1123 — 131250536 — B22636?) )

+ cos(27rp)(:i§6 — 8115133 - 312530536 ~ 822(i6)2)

0
+ 9
( k3 [cos(27r11):i:0 + sin(27Tl-‘)513()]2 )

37

or, carrying all the second order terms to the right hand side and utilizing Eqs. (3.3.3):

551 003(279”) (~ ~2 ~ ~l ~I 2
= — A :1: — A — A ) +
( 2’ ) 2 ( _Sin(2w) > $0 11 0 12230550 22(130)

) (if) - B112-ii — 312530516 - B22036?) +

0

+ ( k8 [cos(27rn):i:o + Sin(271#)ibl2 ) +

) (c032(27r,a)i(2) + sin(47ru):i:0i() + sin2(27r,a)(§:6)2) +
1 . ~2 ~ ~I 1 . ~I 2
——2- srn(47rn):1:0 + cos(47ru):1:0:1:0 + ~2— srn(47rp)(r0) +

A
+ ( B22 ) (sin2(27ra)i(2) — sin(47ru):i:0§:6 +cos2(21r,a)(:i:())2).
22

(3.3.6)

 

Gathering like terms in Eq. (3.3.6) yields a system of six equations allowing to

find all the second order coefficients of the transformation:

1
A11 0082(27r11) — §A12 sin(47r)u.) + A22 sin2(27r)u)-

(3.3.7)
— A11 cos(27rn) — Bu sin(27rn) = 0
A11 sin(47r,u) + A12 COS(47r/J.) — A22 Sin(47l'p.)—
(3.3.8)
— A12 COS(27T/.t) — B12 sin(27rp) = 0
1
A11 sin2(27r)u) + —A12 sin(47r/I) + A22 cos2(27rp.)—
2 (3.3.9)
-— A22 cos(27rn) - 822 sin(27ra) = 0
1
811 cosz(27rn) — 5312 sin(47rp) + 322 sin2(27r)u)+
(3.3.10)
+ A11 sin(27r)u) -— Bu cos(27rp) = —ks cos2(27rp)
Bu sin(47r,u) + 812 cos(47rn) — 322 sin(47ru)+
(3.3.11)

+ A12 sin(27r,u) — 312 cos(27ru) 2 —k3 sin(47ru)

38

1
811 sin2(27rn) + —B12 sin(47rn) + B22 cos2(27rn)+
2 (3.3.12)

+ A22 3mm”) — 822 cos(27ru) = —k, sin2(27r11)
The sum of Eqs. (3.3.7) and (3.3.9) gives:

A11 + A22 — COS(27Tp.)(A11 + A22) —- Sin(27ru)(B11 + 322) = 0. (3.3.13)
The sum of Eqs. (3.3.10) and (3.3.12) gives:
B11 + 322 - COS(2’IT/,l.)(B11 + 322) + Sin(27ru)(A11 + A22) = -k3. (3.3.14)
Eqs. (3.3.13) and (3.3.14) together yield the following:
1 — 008(27r1l) -— Sin(27ru) A11 + A22 _ 0 (3 3 15)
sin(27r/I) 1 — COS(27T/.l.) 811 + 322 -—ks I ' '
hence,
All + A22 _ 1 1 — cos(27rn) sin(27ru) 0 _
811 + 322 — 2(1 - cos(27rp)) — sin(27rn) 1 —- cos(27w) —k3

_ 1 —k3 sin(27rp) _ —%k3 cot(7rn)
_ 2(1 — cos(277n)) —ks(1 — cos(27r,u)) _ —%ks ’

 

 

(3.3.16)

The difference between Eqs. (3.3.7) and (3.3.9) gives:

c03(4r7T/4)(«‘111 - A22) - A12 sin(47m)-
. (3.3.17)
— (A11 - A22) COS-(2711‘) — (Bu - B22) Sin(27w) = 0-

The difference between Eqs. (3.3.10) and (3.3.12) gives:

c03(47r1u)(B11 - B22) - 312 8111(47Tul—
(3.3.18)

— (B11 — 322) COS(27l’/.l) + (A11 - A22) Sin(27r,u) = —k3 COS(47r,u).
Eqs. (3.3.8), (3.3.11), (3.3.17), and (3.3.18) together can be represented in the

matrix form:

Bu - 322 = 1 -M- A11 - A22
312 Sin(27W) A12
M- B11 — B22 + 8111(27w) A11 - A22 = —ks 693(47W)
312 A12 s1n(47r11)

39

, (3.3.19)

where

= cos(47r11) — cos(27ru) — sin(47rp)
M ( sin(47rp) cos(47ru) — cos(27rp) ) '

Eq. (3.3.19) yields

(M2 + sin2(27r/r) - I) ( AHA—1.21422 ) = —k3 sin(27r/.1) < :f:é:::; ) , (3320)

where I is the 2x2 unity matrix. Thus,

. A11--A:22 =_ sin 77 C05(47w)
R( ) .3 .2 ”(W ), (3....)

where

R __ 1 + cos(87ru) — cos(67rp) — cos(27rp) — sin(87m) + sin(67rp) + sin(27rp)
— sin(87rp) — sin(67rp) — sin(27rp) 1 + cos(87r/.t) — cos(67ru) — cos(27ru) °

det(R) = (1 + cos(87ru) — cos(67ru) — cos(27rp))2+

(sin(87w) — sin(67ru) — sin(27r/1))2 = 16 sin2(7r)u) sin2(37w).

Hence, one can invert the matrix R, effectively solving the system (3.3.21):

 

 

 

( A11 — A22 ) = _ ks sin(27rp) ~R—1 . ( cos(47ru) ) =
A12 16 sin2(7r11) sin2(37ru) sin(41ru)
ks cos W (3.3.22)
_ k3 cos(1r11) 2(cos(47r,u) — cos(27r)u)) = W
8 sin(7ru) sin2 (37m) 0 0 '
( 811 — 822 \ = _ k3 cos(7ru) ( cos(47ru) -— cos(27r)u) ) =
B12 / 2 sin(27r11) sin2(37ru) sin(47r,u)
(3.3.23)

 

1 l
1 cos(7rp = k cos(7ru) cos(27r;1) '
st (““3””) + sin 37m ) s sin(37r,u)

40

Ultimately,

 

L

 

 

k3 cos(7r/1) cos(27r11)
A11 = - .
23m(37ru)
A12 = 0
A _ _ks cos3(7ru)
22 _ sin(37rp)
ks
B = -—
11 2
k3 cos(7r)u) cos(27r11)
312 = .
sm(37m)
322 = 0

(3.3.24)

As one can see, solving the system (3.3.7) requires a lot of heavy algebraic and

trigonometric transformations. The most important part is that these transformations

cannot be generalized, and hence, the methods to solve the system will vary depending

on the transfer map under consideration.

If sin(31m) goes to zero, the coefficients A11, A22, B12 grow to infinity, and the

normal form algorithm fails. That is why it was assumed in Section 3.2 that the

equation [I - 7'13 = 0 mod 27r does not have any nontrivial solutions.

To verify that Eq. (3.3.24) is indeed the solution, substitute the resulting values

41

of the coefficients back into Eqs. (3.3.7)—(3.3.12):

1
A11cos2(27r11) - 5.412 sin(47rp) + A22 sin2(27rp) — A11 cos(27m) — B11sin(27rp) =

k3 cos3 (7m)
sin(37r11)

_ k3 cos(7rp) cos(27rp)
2 sin(37w)

 

0032(27r11) -— sin2(27ru)+

ks cos(7rp.) cos(27r/.L)
2 sin(37ru)

 

k
cos(27r,u) + 33 sin(27ru) =

k .
[cos3(27ru) + 2cos2(7w) sin2(27rp) — cos2(27ru)] + w =

__ _ k3 cos(7rp)
— 2 sin(37ru)
k3 cos(7r11)

)

= — 2sin(37rp [cos3(27ru) + cos(27rp) sin2(27r/1) + sin2(27rp.) - cos2(27r/1)] +

2
ks cos(7rp)
= — 2 sin(37r/r)
k3 cos(7ru)
= — 2 sin(37r)u)

k3 sin(27rp) _
2 _
k3 sin(27ru) _
2 _

[cos(27ru) — cos(47rp)] +

2 sin(37m) sin(7r11) + 0.

(3.3.25)
The other five equations of the system (3.3.7) can be verified using a similar technique.
It is already clear after one step that the manual computation of the normal form
transformation is not straightforward, and it does not imply any algorithm which can
be efficiently implemented and suitable for a large class of transfer maps.
On the other hand, the normal form algorithm in the differential algebra frame-
work is perfectly clear and very efficient to any order desired. It is implemented
in COSY INFINITY and used throughout the subsequent chapters to assist solving

various problems of the optimization of accelerator parameters.

3.4 Normal Form Transformation Result

To see the effect of the normal form transformation, let us consider a particular
transfer map and its representation after the normal form transformation.

The transfer map is a function, connecting the initial and final states of some

42

k3 sin(27ru)

deterministic dynamical system:

2‘ = 13(2), (3.4.1)

:2) = we. >. (3.42)

where F is continuous in 2' and continuous almost everywhere in the independent
variable t, the time, or s, the arc length, 2,- and 2f are the initial and final states of
the system, respectively, and 6 is a set of system parameters. Then M is a transfer
map of the system, connecting the initial state 2,- at t,- (or 3,) to the final state 2f at
t f (or sf). In the differential algebra framework one considers the equivalence classes
[M]n instead of M. Each class [M]n can be represented by the set of coefficient
of the Taylor expansion of M about the point z = 0. Depending on the order of
calculation (the highest order of the terms kept in the Taylor expansion), the number
of the coefficients can be very large, and also it grows fast as the calculation order is
increased. Table 3.1 shows the example of the human readable representation of the
transfer map in the particle optical coordinates before the normal form transformation
for a long dipole magnet.

Each row of the table describes one term of the Taylor expansion of the final coor-
dinates 23f and a f (divided by the horizontal line) in terms of the initial coordinates
2:,- and CL). The column labeled “EXPONENTS” describes the exponents of each of
the independent variables appearing in the respective term. The “ORDER” column
contains the total order of the term, that is, the sum of the exponents, and the first
column lists the double precision coefficient belonging to the respective term. For
example, the sixth row in the upper part of the table shows the Taylor expansion
coefficient of the (mlxa2) term, as the exponents are 1 for :1: and 2 for a. At the same

time, the sixth row in the lower part of the the table shows the coefficient of the

43

COEFFICIENT ORDER EXPONENTS
0.6234898018587336
0.7818314824680298
-.3056302334890786
.4874639560909119
.1173751344184454
-.3056302334890786
0.2437319780454559
-.4670491981129435E-01
0.1489837227175753
.1782158313658518
.9474825532788060E-01
.1736113620668300E-01
.9340983962258870E-01
.2234755840763630
.2546233897381215
16 0.1083071221753043
COEFFICIENT ORDER EXPONENTS
-.7818314824680298 1 1
0.6234898018587336 1 0
.3909157412340149 2 0
.1110223024625157E-15 4 1
4 0
5 3
5 2

1.;
O

toooxzoamtbwtou-LH
00

o—t H H 1—s H
op 00 [O H o
O I I O I

H
01
I
OOOOOOOOOOOOOOOO

mmmmpppppwwwwwp
OHMQDOI-‘MOJPOHOHMOH
mpwwwaHowMMI-tohs

 

.9772893530850374E-01
.2220446049250313E-15
.1110223024625157E-15

Nmmpwan—a
O

O

 

 

Table 3.1. Transfer map in particle optical coordinates

44

I COEFFICIENT ORDER EXPONENTS

 

1 0.6234898018587336 1 1 0 0
2 0.7818314824680299 1 0 1 0
I COEFFICIENT ORDER EXPONENTS
1 -.7818314824680299 1 1 0 0
2 0.6234898018587336 1 0 1 0

 

Table 3.2. Transfer map in normal form coordinates

(alx3a2) term, as those elements of the transfer map which are zero are omitted in
the representation.

It is clear from the representation of the transfer map, that the system under
consideration is nonlinear. At the same time, after the normal form transformation
is applied to this map, the resulting transfer map is linear (see Table 3.2). Thus, the
motion of particles in the normal form coOrdinates can be described by the linear

transfer matrix:

if _ 0.6234898018587336 0.7818314824680299 1:; (343)
t} —.7818314824680299 0.6234898018587336 ,3; ’ - -

where (tit?) and (ti, t?) are the initial and final normal form coordinates of the
particle, respectively.

The beauty of the differential algebraic approach to the normal form transfor-
mation is in that the algorithm is constructive, and it allows to obtain not only the
transfer map in the normal form coordinates, but also the transformation to that set

of coordinates and its inverse. Therefore,

23f _ A_10 0.6234898018587336 0.7818314824680299 0 A :17,-
af —.7818314824680299 0.6234898018587336 az- ’

(3.4.4)

where A is map of the transformation to the normal form coordinates:

(:)=A<:). was

45

4. DYNAMICS OPTIMIZATION IN

THE TEVATRON

4.1 Skew Quadrupole Correction

Consider the problem of the optimization of the dynamics of particles in the proton—
antiproton Tevatron collider. Increasing the luminosity reach of existing and future
colliders demands considered and precise optical design and predictability in opera-
tion. Driven by nonlinear fields, high-order beam dynamics is generally difficult to
control, calculate, and can severely limit a machine’s region of stable operation. An
approximately linear lattice is desirable for Operational simplicity and understanding,
and it also generally exhibits more robust, broader-range performance. Nonlinear
sources arising from field and alignment errors and required correction elements are
unavoidable. Successful management of the sources of nonlinearities, however, de-
pends on the linear lattice. Attributes of the linear lattice and relative locations
of sources generate interference, constructive or destructive, between the nonlinear
terms depending on their periodicity. In a highly effective linear lattice design, the
strongest nonlinear amplitudes can be mitigated passively by intelligently exploiting
periodicity, phase advance and optimal placement of nonlinear correctors. Such a

lattice enhances precision and predictability in the machine optics.

46

Passive cancellation, however, is generally not sufficient to address certain system-
atic or widespread field errors, and active correction in the form of added corrector
elements is usually required. The overall lattice approach must be evaluated not only
by its tolerance to errors, nonlinearities, and natural aberrations, but also by its po-
tential for active correction. Such correction may be “global” in the sense that an
error or aberration is corrected over one-turn optics. Global correction is not always
adequate to maintain sensitive collider optics. Immediate — or “local” correction —
of source terms, particularly if such terms propagate through the delicate optics of
the interaction regions, may be an additional requirement for stability and linearity.

The case addressed in this work is the state of the Tevatron collider as of August
2004, where a strong, systematic, skew quadrupole error was present in the operational
lattice as a result of a coil shift in the superconducting arc dipoles.

Technically the Tevatron is a superconducting synchrotron with six straight sec-
tions and six long arcs. Two of the straights are the well-known GDP and D0 detec-
tors, and four other straights are used for minor experiments, acceleration, injection
and extraction. The main part of each arc is occupied by the superconducting dipoles
designed in the late 1970’s to early 1980’s. Due to the extensive use, over time some
of the superconducting coils shifted, thus causing the skew quadrupole component to
appear and grow to the values that could not be corrected using only the “passive”
correctors available. Without introducing additional correctors or correcting the field
errors in dipoles it was not possible to keep and improve the quality of the beam, and
hence the luminosity.

The dynamics of the particles in the Tevatron can be described by a system of

differential equations of the form

2': 13(2), (4.1.1)

with

( x )

y
2': l: W _ t0) , (4.1.2)
a = Px/pO
b = Py/PO

 

 

\6 = (E—EOVEO/

where a: is the horizontal position of the particle relative to the reference particle
that goes through the centers of the magnets, y is the vertical relative position, l
characterizes the deviation in the time of flight with the coefficient k chosen in such
a way that I has the dimension of distance, a and b are the reduced momenta in
horizontal and vertical transversal directions respectively, 6 is the reduced difference
between the total energy of the particle under consideration (E) and that of the
reference particle (E0), and p0 is the initial momentum of the reference particle.

By this definition of 2', the reference particle has coordinates 2' = 0, and thus the
vector 2' indeed describes relative motion.

The general form of the vector 13“ on the right hand side can be found in Sec-
tion 2.1.2 of Chapter 2. As the particles move in the electromagnetic field, the system
(4.1.1) is Hamiltonian, as one can assume that the energy is constant at any given
time, as we consider the state of the Tevatron at injection. The right hand side of the
system possesses continuous partial derivatives with respect to the components of 2',
and almost everywhere continuous (excluding a finite number of points) with respect
to the independent variable 3 in some region U. Thus for any choice of the initial
coordinates in U the system (4.1.1) possesses a uniques solution and the behavior of
the system at s = 0 determines the dynamics of particles Vs > 0. The right hand
side of the system is a periodic function with a period 7' > 0, that is, F(s + T) = E(s)

and there is no such 0 < 7'2 < 7' that F(s + 72) = 15(3).

48

Assume that Eqs. (4.1.1) have been integrated to yield the function
as» = M33323»). (413)

where s,- and 3f are the initial and final values of the independent variable. The
function M describes how individual state space points “flow” as time progresses, and
in the theory of dynamical systems it is usually referred to as the flow of the differential
equation. Hereafter this function is referred to as the transfer map (describing how
state space points are transferred).

Let M = Msz‘fii'” be the transfer map of one full revolution, the so-called first
return map. Then one can trace the dynamics of particles after an arbitrary number
of turns by applying the operator M over and over again to the set of the initial data.

The optimization problem is to uncouple the dynamics of particles in the planes
(at, a) and (y, b). By definition the beam occupies a small subset of the phase space in
the neighborhood of the reference particle with zero coordinates. Hence we can employ
the differential algebra methods [12] and consider along with the transfer map M the
equivalence class [M]n 6 ”By, where n is the computation order, 2) is the number of
phase space variables, and nDy is the corresponding differential algebra. Differential
algebra methods are so efficient that the orders up to 10 can easily be handled by any
home PC. Also, as it was mentioned in Chapter 1, transfer map methods are more
efficient than the methods requiring integration of individual orbits for some tasks,
as one can calculate the transfer map (even with parameters) once and apply it over

and over again.

49

4.1.1 Linear Problem

Let the system (4.1.1) be linear, or alternatively, consider the linearization of the
full nonlinear system in the neighborhood of the origin: 5' = .Iac(I-7‘.)|O - 5. Then the
conditions that the transfer map must satisfy in order to decouple the motion in the

planes (2:, a) and (y, b) can be written as

 

3x0) _ ax(f) _ aa(f) _ aa(f) _

ayb) abb) ayb) 35(2) (4.1.4)
ay(f) _ ay(f) _ 31,0”) _ ab(f) _ 0

6x0) - ’ 6a“) _ ’ ax“) — ’ 6a“) _ ’

where Z = (:c,a,y,b,l,6)T, and the indices (2') and (f) denote the initial and final
conditions, respectively.

In the linear case the transfer map and its equivalence class [M]n can be repre-
sented in the form of the matrix for the coefficients connecting the initial and final

values:

(Z1|21) (21IZv) _.

2(8f)= E 5 (2(2)), (4-1-5)
(Zv|21) (ZUIZ‘U)

where (zilzj) is the coefficient of the dependence of the i-th component of the final
state of the particle on the j-th component of the initial state. Consider only the
transversal dynamics, that is, let 2' = (at, a, y, b)T. In this case the condition for the

decoupling of the horizontal and vertical planes has the form

(SUI-’17) (xla) 0 0

= (alz) (ala) o 0 ,3

(f) 0 0 (yly) (ylb) ((2)) (4.1.6)
0 0 (bly) (blb)

So, to achieve minimal interdependence between the two planes, it is necessary to

50

minimize the following objective function:

 

 

 

 

 

 

 

 

a$(f) 3550) aa(f) aa(f)
1: 5,73 m 5,73 + m +
ay(f) ay(f) ab(f) ab(f)
+ m m 52:75 Ea—(z—I = (4.1.7)

 

 

 

 

 

 

 

 

= |($ly)l + l(rv|b)l + |(aly)| + |(a|b)|+

+ |(y|$)l + Kyla)! + |(bl$)| + |(bla)|.
which depends on up to 8 control parameters, representing the strengths of the skew
quadrupole correctors in the Tevatron lattice, that is, J1 = J1(1‘Z), where the compo—
nents of 11' are the field strengths in the corresponding correctors. Both the compo-
nents of 2', and J1 depend continuously on 11‘ over some compact set for the parameters

17.

4. 1.2 Nonlinear Problem

If the transfer map M is nonlinear, then the condition for the transversal motion to

be decoupled assumes a more complicated form:
6i1+i2 (f) ai1+i2 (f)
(1) ,1 53’“ 0 2.2 =0 (') 2.1 y’“ 0 2.2 =0, k=1,2, (4.1.8)
2 Z Z Z
8(a) 3(a) 8(a) 8(32)

where E is given in the indexed form for notational convenience:

 

 

Z: (x1,x2,y1,y2)T = (ac,a, y, b)T, 2'1, i2 6 Z20; 2'1 + 22 6 171, n is the highest
order of the nonlinearities taken into account, indices (2') and (f) denote the initial
and final values, respectively.

In the nonlinear case the equivalence class cannot be represented in a matrix
form. Instead, the following notation is introduced to denote the coefficients of the

high order Taylor expansion

(2,-1sz . . . zgv), (4.1.9)

51

where v is the number of phase space variables, 3' is an index of the component of
20); 22:1 ik 6 I717, ik E Z20, n is the order of computations (the highest order in
the Taylor expansion). It is not hard to see that the linear matrix (4.1.5) follows this
notation. For the linear case n = 1.

)T

Using this notation for 2' = (2:1, 2:2, y1,y2 one obtains the following conditions:

(mlyily32)=0, (leyily32)=o, (yllxilx32)=0, (yzlxilx32)=0. (4.1.10)

for all possible combinations of 2'1, 2'2 6 Z20, such that i1 + 2'2 6 T,—n.
Thus, the nonlinear problem of uncoupling the dynamics in the horizontal and

vertical planes can be reduced to the minimization problem for the following objective

 

 

function:
2 n ail+i2$(f) 77' ajl +j2y(f)
J2 = _ . k . . + . - k . - =
E 33.3 W()(”) 2= WM)
2 n . . n - -
= Z Z: l($k|ylly;2)l+ Z |(yk|xilxg2)l , (4.1.11)
k=1 i1+i2=1 J'1+J'2=1

where, similar to the section 4.1.1, J2 depends on up to 8 control parameters, that
is, J2 = .1207), and the components of 17 are the strengths of the fields in the skew
quadrupole correctors. E and J1 depend continuously on 17 over some compact set for

the parameters 11'.

4. 1.3 Tune Preservation

In Sections 4.1.1 and 4.1.2, the statements of the two problems — linear and nonlinear
— for uncoupling of the transversal dynamics of the particles in the horizontal and
vertical planes were discussed. Changes made to the skew quadrupole content of the

lattice might result in the change of the tunes (eigenvalue phases of the transfer map).

52

Such a change is undesirable, and therefore along with one of the two problems of the
skew quadrupole correction one should state the problem of the preservation of the

tune. The constant part of the tune can be found using the following relations:
#3: — arccos 5 m + m ,

1 ay(f) ab(f)
[1y — arccos § Ell—('7 + m ,

is the vector of the phase space variables of the system (4.1.1)

(4.1.12)

where 2': (x, a, y, b)T

under consideration.

Hence, one should add one more objective function to the list of the functions to
1 3;;(f) aa(f)
um - arccos 5 5m)— + m +

l _ S l ay(f)+ab(f)
[“y me" 2 aye) aha) '

The strengths of the main quadrupoles are used as the control parameters for this par-

minimize:

J3:

 

(4.1.13)

 

ticular problem. As all the focusing quadrupoles and all the defocusing quadrupoles
are set to the same value of the strength throughout the ring, there are two indepen-

dent parameters to control.

4.1.4 Introduction to the Solution Methods

Next sections address the methods to solve the problems from Sections 4.1.1—4.1.3 and
corresponding algorithms, implemented during the work on the dissertation from 2003
to 2006. The computer implementation of the algorithms is done using Mathworks
Matlab R13 and COSY INFINITY versions 8.0—9.0. An additional utility program to

convert lattices [45] is written in PHP version 4.4.4 (PHP: Hypertext Preprocessor)

[1].

53

FODO 1
FODO 2
FODO 3
FODO 4
FODO 5
FODO 6
FODO 7
FODO 8
FODO 9 _
FODO 10 .
FODO 11 5
FODO 12 ,
FODO 13 : '

FODO 14 I
FODO 15

 

 

 

 

 

 

 

 

 

 

 

 

E 6*“ I 0:52
I 5'2 l D???

   
   

 

 

V SQC

    

 

Figure 4.1. Lattice scheme before any correction

According to the formulation of the problems in Sections 4.1.1—4.1.3, it is neces-
sary to minimize the function J1 or J2 depending on the kind of the problem under
consideration, and the function J3, the minimization of which guarantees that the
tune deviation from the initial value is minimal. These two problems can be solved
one after another, as the control parameters for J1 and J2 are the strengths of the
skew quadrupole correctors, the number of which varies from 6 to 8 for each sector,
while for J3 the control parameters are the strengths of the main quadrupoles —
focusing and defocusing.

It is desirable though to keep all the skew quadrupole correctors at the same
strength due to the features of power supplies. So, the number of various control
parameters goes down to 6 in all arcs. For an additional optimization step one can
use different sets of parameters. Straight sections also have four pairs of correctors
and if one wants to vary these pairwise, which gives 4 more control parameters.

Figure 4.1 shows a lattice scheme of the typical arc of the Tevtaron.

54

The scheme employs the following notation:
<FQ> — focusing quadrupole;
<DQ> — defocusing quadrupole;
<D*2> —— two sections of the main bending magnet;

<D*2 FIX> —-— two sections of the main bending magnet to be fixed during the

shutdown of August 2004;

<SQC> — skew quadrupole corrector, minimizing the influence of the skew

quadrupole errors in bending magnets.

4.1.5 Skew Quadrupole Error Sources

During the Tevatron operation since early 1980’s, due to the stress from the strong
magnetic fields, the superconducting coils of the 774 main bending dipoles shifted,
thus introducing a. significant skew quadrupole component in the expansion of the
magnetic field [57], which could not be effectively mitigated by the existing family of
skew quadrupole correctors. The status of the collider as of the beginning of August
2004 was as follows: about 15% of all the dipoles were fixed, and the other 85% still
had the skew quadrupole component in them. According to the plan of the shutdown
works, up to 50% of the dipoles (total) could be fixed. So, there was a need of an
effective pattern of the dipoles to fix to simplify the correction issue in such a way
that the existing family of correctors could uncouple the motion in the horizontal and
vertical planes. After the effective scheme was found the corrector strengths were
set to the values minimizing the functions J1, J2, J3 (one or more depending on the

nature of the problem).

55

4. 1.6 Initial Data

As all the calculations are performed in terms of the high-order Taylor expansions, it
is necessary to first obtain the representation of the transfer map of one full revolution
in terms of Taylor expansions. This is done using the code COSY INFINITY (COSY),
which implements multidimensional differential algebra. Among the conveniences of
such an approach is the fact that certain families of the optics elements can easily be
turned on and off in order to study particular effects and phenomena, their correlation
and the impact on the dynamics of particles in general.

The initial data — the description of the system lattice — were obtained from
one of the Fermilab employees, Valeri Lebedev, in the format of the code list in-
tended to run in the beam optics computation program OptiM [32]. To convert
the existing lattice to COSY internal format, a converter program was written (the
detailed description can be found in Appendix A). It is available on the web at
http://cosy.pa.msu.edu/converters/optiM2cosy/. PHP 4.4.4 [1] was the lan-
guage of choice for the converter. The main reasons for this choice were the simplicity
of the code writing and debugging and the fact that PHP is a web-oriented program-
ming language. As a result, the converter code can be used at virtually any web-server
supporting PHP. Currently the following set of optics elements is supported: dipoles,
quadrupoles, skew quadrupoles, sextupoles and skew sextupoles, thin lens approxima—
tion quadrupoles and sextupoles, solenoids, electrostatic separators, and octupoles.
The file generated by the converter can be used in COSY with no further change to
calculate the fifth order transfer map of one full revolution.

To ensure the conversion program works flawlessly, a number of tests were per-
formed. First of all, the linear transfer matrices were compared among COSY, OptiM

(Fig. 4.2), and also an independent implementation in COSY by Bela Erdelyi [22].

56

 

 

 

Figure 4.2. Beta functions comparison, top - OptiM32, bottom - COSY INFINITY

57

Also the values of the beta functions [57] and the dispersion functions were com-
pared. The quantitative result showed the discrepancy between the models within
one percent, which was most likely caused by the differences in the solenoid element

implementation in different codes [32,35].

4.1.7 Model Details and Tracking Results

Typically 11th to 15th order maps were required for the complete convergence of the
nonlinear effects (the phase portraits show no or little discrepancy when the order is
increased further), but lower orders (7th, for example) provided a quicker check and
served as the preliminary optimization results. The calculation order corresponds to
the degree of the Taylor polynomial which is substituted for the function on the right
hand side of the system (4.1.1).

Particles were launched at the CDF interaction point in steps of one sigma (for the
normalized emittance of 107r at injection, 0 = 1.2-10’4 m). Particles were tracked in
COSY by applying the transfer map repetitively for typically 10,000 turns. Only the
injection optics was being studied. It is important to note that the study is not, per
se, a dynamic aperture one for which particles are launched along phase space vectors
scaled to the linear injection ellipse and the transmitted transverse phase is mapped.
In a predominately linear lattice, tracking along a single vector in one plane of phase
space and then the other is sufficient to trace out the matched ellipse. Particles can be
launched along the :c or y axis, for example. The degradation of the linear motion as
evidenced by dissolution or distortion of the linear invariant ellipses is then studied.
Since the current study is directed toward optimizing the linear performance, this
is the approach used for tracking and the criterion for improvement. In some cases

other directions were also considered, for instance, the diagonals of the planes (513,0)

58

and (y, b).

The tracking results presented in this and subsequent sections are obtained for
10,000 turns with points plotted every 10th turn, and the scales are 2.7 - 10‘3 m
for the :1: and y axes and 4.0 - 10—3 for the a and b axes, if it is not mentioned
otherwise. The tracking is performed with a symplectification algorithm written by
Bela Erdélyi [23—25, 44] for the calculation order 7. All the particles are launched
either along the a: or y axis, which is explicitly mentioned in each figure caption. To
start comparing the impacts of different sets of the nonlinear elements, in Fig. 4.3
the phase portraits for the linear motion are shown. This includes only pure dipoles,
pure quadrupoles and no skew quadrupole errors in dipoles, and hence there is no
coupling between the (2:, a) and (y, b) planes.

According to the status of the Tevatron before August 2004, with 15% of the skew
quadrupole errors removed in selected dipoles, the otherwise unchanged lattice shows
significantly reduced stability (Fig. 4.4). Regular structures disappear and most of
the outer particles can be considered lost in just 10,000 turns, while the realistic
simulations require millions of turns. These phase portraits can be considered a

starting point of the study of different schemes of the skew quadrupole correction.

4.1.8 Skew Quadrupole Circuits Optimization Proposals

As the scheme suggests (Fig. 4.1), each of the Tevatron arcs in the model framework
has 15 FODO cells with the skew quadrupole correctors in every odd-numbered cell,
which means one corrector every two FODO cells. The “0” section of the FODO cell
contains four sections of the bending magnet, high order magnets, and various cor-
rectors. The skew quadrupole correctors are placed next to the horizontally-focusing

quadrupoles only.

59

 

 

(a) (11:, a) plane

 

 

 

(b) (y, b) plane

Figure 4.3. The (x, a) plane and (y, b) plane phase portraits, only dipoles and pure
quadrupoles are active, and particles are launched along the an axis and the y axis,
respectively. The scales are 2.7- 10‘3 m for the :z: and y axes and 4.0- 10‘3 for the a

and b axes.

60

 

(a) (:c, a) plane

 

(b) (y, b) plane

Figure 4.4. The (x,a) plane and (y, b) plane phase portraits before the optimiza:
tion, particles are launched along the cc axis and the y axis, respectively. The phase
portraits include all sextupole and skew quadrupole fields (correctors plus errors) in
addition to quadrupoles and dipoles, where 15% of the skew quadrupole errors have
been removed in specific dipoles. the scales are 2.7- 10‘3 m for the :c and y axes and
4.0- 10‘3 for the a and b axes.

61

It is clear from Fig. 4.4 that one family of skew quadrupole correctors is not
sufficient to correct all the skew quadrupole errors in dipoles along the ring. However,
if 50% of the coil shift errors in dipoles are fixed, one circuit of the skew quadrupole
correctors is capable of removing the coupling in the arcs. The problem here is to
discover both the optimal dipole pattern for error correction and the new strengths
for the skew quadrupole correctors.

An optimization with the corrector strengths set to different values is not practical.
All the correctors in each are have the same power supply, so it is more realistic to
use one strength for all the correctors arc-wise or even ring-wise.

The optimization process itself consists of two steps. First, the optimization of
each arc is performed using the skew quadrupole corrector strengths as control param-
eters. This optimization would be close to optimal if no skew quadrupole components
existed in the straight sections of the Tevatron, but there are skew quadrupole error
sources and correctors in the interaction regions. Because of these components and
the residual skew terms from the arcs since the arcs are not perfectly regular, the
one-turn transfer map has nonzero skew quadrupole terms which require correction
also. To remove this smaller, final stage of coupling a second step of the optimization
is required. In four of the six straight sections there exist eight skew quadrupole cor-
rectors and the strengths of these correctors are used to finish the skew-quadrupole
term cancellation in the one-turn map.

As a result of the optimization, two optimization schemes were selected out of six,
which differed in the dipole pattern used for correcting the skew quadrupole error.
The first scheme, proposed by the author and Dr. Carol Johnstone (Fermilab) [48],
attempts the elimination of skew quadrupole error source predominately in the ver-

tical plane by fixing the two dipoles flanking each vertically-focusing arc quadrupole.

62

Sector E scheme I

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

F000 1 EH2 Frxj
FODO 2 [an ax]
FODO 3 [Tan FIX]
FODO 4 I 0*: FIX]
FODO s [on mg
FODO 6 Lou FIX]
FODO 1 [m2 HZ]
FODO 8 [m2 FIX]
FODO 9 [on FIX]
FODO 1o [0*2 Fm
FODO 11 Lon F15]
FODO 12 . You mg
FODO 13 ,._,L'.. . WW F171
FODO 14 "‘ [on six]
FODO 1s ' [0*2 FIX]

 

 

Figure 4.5. Correction scheme 1: the skew quadrupole error is removed in dipoles in
each cell surrounding the defocusing quadrupole.

With one degree of freedom corrected, the remaining predominately horizontal sources
can be corrected with the existing single family of the skew quadrupole correctors.
The layout for Scheme I is shown in Fig. 4.5. One arc of the Tevatron is shown, and
the others are similar. The dipoles with the skew quadrupole error are marked with
the “D*2” symbol, and the fixed dipoles with the “D*2 FIX” symbol. The vertically-
focusing main quadrupole is marked with “DQ” as it is defocusing in the horizontal
plane.

The second approach proposed by Dr. Michael Syphers [51] is to correct the skew
quadrupole errors in each FODO cell missing a skew quadrupole corrector. This
scheme was further improved by removing specific correctors from the single family
to provide for more consistent correction in each are as a function of fractional phase
advance. The Tevatron has a periodicity of five, that is, if we consider the linear

dynamics, the phase advance of the particles after one cell is 72 degrees. Hence, after

63

five cells the particle has the same phase. Taking into account that all the skew
quadrupole correctors in each are have the same strength, it makes sense to have
the same number of correctors for the following groups of cells to avoid over / under-
correction: (1,6,11), (2,7,12), (3,8,13), (4,9,14) and (5,10,15). Looking at Fig. 4.1 (the
situation before the correction), the triplets (1,6,11), (3,8,13) and (5,10,15) have two
correctors, while the triplets (2,7,12) and (4,9,14) have one. For the Scheme II it is
proposed to turn off the correctors in the cells 1, 3, 15, so that every triplet of the
cells has one skew quadrupole corrector.

The underlying idea of the scheme is to correct the error locally at the source
and remove it from cells without local correction. This scheme gives improved per-
formance, particularly using one uniform corrector strength across all the arcs. The
layout for the scheme is given in Fig. 4.6. The skew quadrupole error sources reside
in odd cells with the skew quadrupole correctors marked with “SQC”. Removing part
of the correctors improved performance further, and these removed correctors are

marked with “SQC RMV”.

4.1.9 Optimization Results

The optimization itself is done in three stages. At first, the one-parameter minimiza-
tion of J1 or J2 is performed for each particular sector out of six (the transfer map
(4.1.3) of each individual arc is considered instead of the one-turn transfer map), as
all the sectors are somewhat different. For example, in the straight sections between
the sectors A and B, and C and D the two detectors CDF and DZero are located,
which imposes some restrictions onto the last cells of A and C, and the first cells of
B and D. This first optimization stage would be sufficient if there were no straight

sections connecting the arcs. The straights too have optics elements alternating the

64

Sector E scheme 11

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

FODO 1o
FODO 11
FODO 12
FODO 13
FODO 14
FODO 1s

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 4.6. Correction scheme II: the skew quadrupole errors in dipoles are fixed in
each even cell.

skew quadrupole component of the field as well as the skew quadrupole correctors.
So a second, “global” stage of the optimization is required with the same functions,
but with a different set of control parameters. Four pairs of correctors in the straight
sections or all of the available correctors (including these in the arcs) could be used
in the minimization procedure. Hence, the number of parameters varies from 4 to 10.

Various optimization methods and sequences of methods were used to solve the
problem. Among them are the simplex method [41], the modified Newton method [33],
and the simulated annealing method [31].

After the values for the skew quadrupole correctors minimizing J1 or J2 are found,
it is important to ensure that the change of parameters did not change the tune. For
the Tevatron #1: = 0.585, My = 0.575, and therefore all the values on the right hand
side of (4.1.13) are known, and one can search for the minimum of J3 using the

strengths of the main focusing and defocusing quadrupoles. They have the same

65

strengths for all the sectors, so there are two control parameters.

The phase portraits of the system after the multi-stage optimization according
to Scheme I are shown in Fig. 4.7. This scheme proved to be more stable under a
perturbation of the system parameters, so it could be more reliable from the practical
viewpoint, but it is less efficient in the sense of the horizontal and vertical plane
decoupling.

Fig. 4.8 shows the phase portraits of the system, optimized using the Scheme II
layout. From the qualitative point of view, this scheme application results in more
closed uncoupled trajectories: up to 200 in the horizontal plane and 12a in the vertical
plane, which is noticeably better than Scheme I. It also stays efficient even if all the
correctors in the arcs are set to one strength throughout the ring. Nonetheless, it
is worth noting that this scheme is sensitive to the small changes of the system
parameters.

Scheme II with minor changes was implemented during the planned Tevatron
shutdown in August 2004, and proved itself very efficient as a temporary solution
before all the errors in dipoles could be corrected.

As the main aim of the optimization was the decoupling of the dynamics in the
horizontal and vertical planes, to demonstrate the results, it is worth comparing the
linear parts of the transfer maps before and after the optimization, represented in the
matrix form for convenience.

The linear transfer map before the optimization is

—0.7553149 0.3637584 0.0663166 —0.0971958
—0.8640196 -0.9507699 —0.4305946 0.1421480
—0.0251393 —0.0621722 —0.8060247 0.3353297
—0.4966847 0.0614653 —0.7083212 —0.9862029

, (4.1.14)

66

 

(a) (1:, a) plane

 

(b) (y, b) plane

Figure 4.7. The (x, a) plane and (y, b) plane phase portraits after the optimization
with 50% skew quadrupole errors in dipoles, errors are fixed according to Scheme I:
around each defocusing quadrupole; particles are launched along the :1: axis and the
y axis, respectively. The scales are 2.7 - 10‘3 m for the x and y axes and 4.0- 10‘3
for the a and b axes.

67

 

 

(a) (z, a) plane

 

(b) (y, 1)) plane

Figure 4.8. The (x, a) plane and (y,b) plane phase portraits after the optimization
according to Scheme 11: the skew quadrupole errors are removed in the two dipoles
flanking each horizontally defocusing quadrupole, and particles are launched along
the :1: axis and the y axis, respectively. The scales are 2.7~ 10’3 m for the z and y
axes and 4.0 - 10‘3 for the a and b axes.

68

and after the optimizaton —

-0.8023857 0.3107970 0.0059011 —0.00254055

—0.7143330 —0.9696470 —0.005807 5 —0.00487758

0.0043365 —0.0022290 —0.8445417 0.28890307
-0.0077313 -0.0060658 —0.5409385 —0.9990799

(4.1.15)

Matrix elements responsible for the cross-coupling between the planes are marked
with bold font (four elements in the upper right corner and four elements in the lower
left corner). Eq. (4.1.15) shows that after the optimization these terms become up to

74 times smaller as compared to Eq. (4.1.14).

4.1.10 Summary

Studies conducted on the skew quadrupole correction case showed that if the skew
quadrupole errors are left unchanged, one family of the skew quadrupole correctors
is not enough to provide a sufficient correction, or the change of locations of the
correctors is in order, which is not a feasible option. At the same time, if 50% errors
are fixed, one corrector family can solve the dynamics decoupling problem successfully,
even with all the correctors in the arcs set to one particular strength.

All the calculations were carried out in COSY using the following two sets of

parameters:

0 draft optimization cycle: 104 turns, 7th calculation order — takes about 5—10

minutes on an average laptop: Intel Celeron 1.5 MHz and 256 Mb RAM;

0 final optimization cycle: 106 turns, 11th calculation order —— takes 8 to 12 hours

on the same machine.

It is worth noting that even the 7th order of calculation returns reliable results, allow-

ing to estimate the efficiency of the optimization approach before turning to a more

69

lengthy cycle. Figure 4.9 demonstrates that the phase portraits for the parameter

sets mentioned above are very similar, which justifies the use of the 7th order.

4.2 Calculation of Nonlinear Tune Shifts with

Amplitude

The tune of the system is one of the most important characteristics of the dynamics
of particles. For the linear system the tune stays constant, while in the nonlinear case
it might change, mainly depending on the position of the particle in the beam (the
so—called tune shift with amplitude), but also depending on other parameters of the
system.

Consider the problem of evaluation of the tune shift with amplitude in the non-
linear case using some extra information obtained by the specific kind of measure-
ments [46]. All the proposed methods have been tested on the Tevatron model [38]
and measurements [6], but that does not mean that the algorithm for finding the
tune shift with amplitude is not applicable to other machines. In fact, the algorithm
should stay valid for any other synchrotron, as long as one can proceed with a lin-
ear normal form transformation (see Section 3.1 for conditions). The normal form
transformation is the core of the method.

Suppose that one only has the information on the linear component of the dynam—
ics of the particles in the accelerator. Assume that there is some extra information
available: the size of the beam, the particle distribution type and also the results of
the special type of measurements of the beam position. The corrector is introduced
into the accelerator optics to kick the beam in the horizontal or vertical direction.

Once the strength of the corrector is known, the displacement of the center of the

70

 

(a) 7th order, 10000 turns

 

(b) 11th order, 1000000 turns

Figure 4.9. Calculation order comparison and its influence on the phase portraits in
the (x,a) plane. The scales are 2.7- 10"3 m for the 0: axis and 4.0- 10‘3 for the a
axis.

71

Average amplitude of the beam

 

 

0 2000 4000 6000 8000
Number of turns

Figure 4.10. Measurement results: horizontal position of the center of mass over a
number of turns and its envelope

beam can be found. After the corrector is turned on and off instantaneously, the
amplitude of the beam center of mass decreases due to the filamentation of the beam,
not the damping, as the motion is symplectic. The position of the center of mass
of the beam is then registered after each turn of the particles. One sample of the
measurement data for the horizontal position is shown in Fig. 4.10.

The normal form transformation yields that in the nonlinear case the tune can be

represented in the following form:
_ 2 4
u — uo+c1r +c2r +..., (4.2.1)

where no is a constant linear tune, c1, c2 are the coefficients of the higher order terms

in the expansion of the dependence of the tune a on the particle’s amplitude in the

72

 

normal form coordinates, where the amplitude is defined to be r = \/(t+)2 + (t‘)2
for the particle with normal form coordinates (t+, t").

As it is assumed that the linear component of the dynamics is known, the value
of #0 is also known for each pair of the conjugate variables describing the transversal
dynamics, while there is not enough information to find the coefficients in the ex-
pansion (4.2.1). The task is to find the connection between the number of turns N
required for the amplitude of the central particle after the kick 7‘(N) to fall to the
values comparable to the amplitude before the kick, which means that the absolute

value of the derivative 1" (N) of the amplitude falls below the certain value.

4.2.1 Objective Functions for the Study

The purpose of the study is to restore at least the cl coefficient in the expansion (4.2.1)
using the measurement results. If multiple measurements are available with various
initial kicks, it might be possible to restore the second coefficient in the expansion,
02.

Remark. N might be hard to find, so as an alternative one can use the following
standard trick: take the number of turns N1/2 required for the amplitude to fall to a
half of its initial value. This would be a good reference point for the fall-off curve.

Hence, the objective function assumes the form
J4(a, b) = |n(a, b) - N l, (4-2-2)

where {77 : 7"]77 < K, r’lm > h) Vm < 17}, T'lm is the rate of change of the beam center
of mass amplitude after k turns, and n is the initial value of the amplitude before the

kick;
J5(a,b) = I771/2(a,b) — N1/2I, (4-2-3)

73

where {711/2 : 7101/2 S %r|0}, rlnl/2 is the beam center of mass amplitude after 171/2
turns, and r|0 is the beam center of mass amplitude immediately after the kick (the
peak amplitude).

In both cases 77 = 0(7), and in turn, r = r(:rc, yc) where arc, ye are the coordinates
of the center of mass after the kick in the normal form coordinates. Depending on
the distribution of particles in the beam and the displacement of the beam due to the
kick, the formulas for are, ye will be different, and for certain cases they are derived
in Sections 4.2.3—4.2.5. The control parameters are the coefficients 01 and 02 in the
expansion of the tune (4.2.1).

All the numerical results of the subsequent subsections are based on the Tevatron
model available on the Fermi National Accelerator Laboratory website [38]. The
source code for the lattice is in the format of the MAD programming environment
[3], for which a converter to COSY INFINITY is readily available. It is currently
maintained by Dr. Kyoko Makino at Michigan State University [2], so there was no

need to reimplement the converter.

4.2.2 Calculation Results versus Measurement Results

After the normal form transformation all the particles follow circles with angular
velocity depending on the amplitude. This is the key fact allowing to establish a
connection between the nonlinear tune shift with amplitude and the behavior of the
beam.

Figure 4.11 schematically shows the positions of four particles with initial positions
chosen along some fixed direction (ray), after several turns. Particles cannot leave
their corresponding circles, but the rotation frequencies are different for different

radii. Assume that the outer particles move faster than the inner particles. In this

74

 

 

 

 

Figure 4.11. Behavior of the particles in the normal form coordinates

particular case the outer particle will leave the inner particle behind in the phase.
As a result of such a redistribution of particles the center of mass of the beam shifts
toward the origin of the coordinate system and then oscillates around it.

Hence, if the center of the initial distribution of particles is displaced from the
origin by means of an instantaneous transversal kick, as it is done to measure the
trajectory with the beam position monitors (BPM), then the amplitude of the center
of mass in the normal form coordinates decreases until it reaches a stable value.

As it is assumed that an accurate linear lattice description is available, one may
use the linear normal form transformation, for which the information on the linear
dynamics is sufficient, to get the information on the distribution of the beam in the
linear normal form coordinates after the kick. The linear normal form transformation

is discussed in great detail in [12], for a general transformation see also Chapter 3.

75

In such a coordinate system all the particles follow circles with the same angular
frequency. Hence, the linear transformation does not give enough information on the
nonlinear tune shift. At the same time the linear transformation is sufficient to get
an approximate initial distribution of the beam in the normal form coordinates.

2

It is a fact that the nonlinear tune shift with amplitude depends on 1' ,7‘4, and

so forth, where r is the amplitude of the particle in the normal form coordinates:

 

r = \/(t+)2 + (t‘)2, that is,
_ 2 4
u(r)—u0+clr +c2r +...,

As a rule, 017‘2 is the dominating term in the expansion. Accelerator designers try
to avoid high order nonlinearities, unless there is a specific need of them. Hence,
finding the coefficient c1 is the most important part. Later, if there are multiple
measurements available, the coefficient 62 could be attempted to be found as well.
In the nonlinear case all the particles of the beam follow concentric circles with

frequencies depending on the amplitude, so the transfer map has the form:

_ cos 2111(7) —sin 27m(r)
M _ ( sin 27r/1(r) cos 27rp('r) ) ' (4.2.4)

If the transfer map M is known, one can track the behavior of particles for arbi-
trary many turns. That, in turn, allows to find the number of turns corresponding
to the moment when the center of mass amplitude is back to its value before the
kick, N, or it is at the half of its value after the kick, N1/2. This establishes the
connection between a and N(N1/2). The number N(N1/2) can be found from the
measurements (Fig. 4.10). Hereafter, the common notation used is N, unless it is not
obvious, whether N or N1/2 is implied. Essentially these are the two different points
on the graph characterizing the same curve.

The general scheme for establishing a connection between c1 and N is as follows:

76

1. The outer particle of the beam having the amplitude R rotates with a frequency
MR), and hence, after 1/(11(R) — #0) turns this particle phase advance is 211
bigger (or smaller, depending on the sign of c1) than that of the particle close to
the origin; assume that the kick is weak enough and the beam is not displaced
too far from the origin, and in fact, in most cases the kick is such that the origin

is still inside the part of the phase space covered by the beam;

2. R depends on the strength of the kick and the initial particle distribution, the
value of R can be found as the maximum of the deviations of particle positions
from the origin after the linear normal form transformation, that is, all the

components to find R are known;

3. The value of c1 is not known, but one can always fix a certain c1 and using the
form of the transfer map (4.2.4) obtain the value of N as a function of c1 and

R;

4. Once the algorithm for finding N (cl) is established, it can be used multiple
times to obtain the correct value of the coefficient c1 for a known value of N
inferred from the measurements as discussed above; it is a typical one-parameter

optimization problem.

Hence, the problem under consideration has been reduced to establishing a de-
pendence of N on various values of cl and R. Depending on the initial distribution of
particles this can be a complicated task, which is not possible or not feasible to solve
analytically to get an explicit expression for c1 = c1(N). The next three subsections

are dedicated to solving this problem numerically for various distributions.

77

 

 

 

Figure 4.12. Uniform particle distribution in the sector
4.2.3 Sector Approximation, Uniform Particle Distribution

Let our initial distribution be uniform in the sector bounded by the two radii R1, R2
and two angles (p1, (p2 (Fig. 4.12). After N turns each particle of this distribution

will have the phase advance of
0N0“) = 27rN/J(T) = 27rA7010 + crr2 + W“)

(orders up to 4 are taken into account). Hence, the particle with radius R1 < 'r < R2
located on the front (back) line of the distribution will have a phase difference of
A0N('r) = 21rN (11(1") — 11(R1)) with respect to the inner particle.

To find the centroid of the resulting planar figure, bounded by two radii R1, R2

78

and two curves given by (p1 + 9N(r), (p2 + 6N(r), three integral formulas are used:

= // rdrdbl;
:1: — l//r2 cos6drd0'
C — S 1
1 2 ,
ya = E 1' sm 0drd0.

For the case under consideration:

R2 902 +0 N (r)
9:5,”) = l f 7'2 cos 6d6dr;
S R1 €p1+0N(r) (4 2 5)
(N) R2 902+0N (1‘) . .
yo = f 7'2 sin 6d6dr.
8 R1 <p1+0N(r)

(N) (N)

Hereafter, 3:0 and ye are the coordinates of the beam center of mass in the normal
form coordinate system (t+, t“).

Let us simplify the form of the last two integrals. Without loss of generality one
can assume —<,01 = 902 = (p (the angle can be changed as only the radius is the
quantity of interest). In addition to that, the coordinate 9 is changed to 2b + 0N(r).

Then one has d0 = dip, and the integrals transform to:

<p+0 N
:rgN) - R2/ 7‘2 cos 6d6dr— —
3 R1

(p+9N(r)
= R2]: 7‘2 cos (211+0N(r ))dv,/2dr; (4.2.6)
8' R1
<p+0 (r)
ygN) = R2/ N 7‘2 sin 0d6dr— —
E R1 <p+0N(r)
= R2/: 7'2 sin (z/J+0N(7‘))dzbd'r. (4.2.7)
5 R1

After integrating with respect to «[2 under the remaining integral one can use the

addition formulas:

cos(a + b) = cos a cos b — sin a sin b,
(4.2.8)
sin(a + b) = sin (1 cos b + cos a sin b,
79

-
.—

which ultimately results in

R
$52M = % f1; 7‘2 {sin(<p + 6N0» — Sin(—<p + 0N(7‘))} dr =

1 R2
= S/ r2 {singocosON + coscpsin 9N + sincpcos 6N — cosgosin 0N} dr =
R1

2 R2
. 2
= — srncp/ r cos 6N(r)dr.
5 R1

R
ygN) = %A12 7‘2 {“ 005(80 + 9N(7‘)) + cos(—<p + 6N(r))} dr =

1 1R2
7‘2 {— coscpcosBN + sincpsinQN + coscpcosQN + sincpsin 6N} dr =

122
= — sin 90/ 1‘2 sin 0N(r)dr.
S R1
(4.2.9)
These last integrals cannot be found analytically due to the polynomial nature of

the argument 6N. Even if one assumes 6 N oc 1‘2, the result is a complicated expression

given in terms of the Fresnel functions (Fig. 4.13, 4.14):

 

 

 

 

JR2
2 __ . 2 __ - 2 __
£21 1' cos 0N(r)dr — 47rNc1 (R2 sm(27rNc1R2) R1 sm(27rNclR1))
1 1/2 1/2 ,
_ 87r(Ncl)3/2 (FresnelS(2(Nc1) R2) — FresnelS(2(Nc1) R1)) ,
R 7‘ sm 6N('r)dr = -47rNc1 (Rgcos(27rNc1R2) — R1cos(27rNclR1)) —
1
1 1/2 1/2
_ 87r(Ncl)3/2 (FresnelC(2(Nc1) R2) — FresnelC(2(Nc1) R1)) ,

(4.2.10)
where FresnelS(a:) = fom cos(zzrtzfilt, and FresnelC(:c) 2 f6” sin(17z't2)dt [4]
The shape of the graph of the Fresnel functions explains the behavior of the beam
center of mass shown in Fig. 4.10. To see even more resemblance, some individual
BPM measurements are shown in Fig. 4.15.

This fact is easily explained by the properties of the Fresnel functions: both

80

FresnelC(x)

FresnelS(x)

 

0.6 -
0.4 -
0.2 -

-0.2 *

 

-
.
.
.
...........................................................................................................
.
........................................................................................................
.
o
.
.
x ..............................................................................................................
0 '
p ...........................................................................................................
c .
.
m ‘
o
.
.
h '
............................................................... ‘................................\. ............
m '
.
.

 

 

Figure 4.13. FfesnelS function graph

FresnelC(x)

 

06-.-
' I l .i
. 4 . .
. . .
. 4 o
. . n
. - .
p .......................................................................................................................... I
. . o
I . . .
. - .
. . .
o2. ............................................................................................................................... .
I . . .
. . .
. . .
.

0
n
n
v
n
.
o
.
o
o
u
c
1
n
.
4
r
a
u
a
u
u
u
o
u
'
n
u
e
n
o
o
-

 

Figure 4.14. FresnelC function graph

81

ooooooooooooooooo

 

Average amplitude of the beam

 

 

 

 

 

 

 

 

 

 

'30 2000 4000 6000 8000
Number of turns
.3
3x10 . g .
E 5
01
'8 2. ............................................................................................................................ .
d)
5
B 1 I ................................................................................................................... .
0 .
'5 ; l
g o *‘F l] .......... I
a x
g _1 .............................................................................................................. l
0 l
8’ l
g _2 . .......................................................................................................................... .
'30 2000 4000 6000 8000
Number of turns

Figure 4.15. Measurement samples

82

Properties of Fresnel Functions

l I r

 

 

FS(x)
- - - FS(1 .5x)
. —FS(1.5x)_—FS(x)

  
 

 

 

 

   

999
[Ohm

 

 

 

 

FS(x), FS(1 .5x),FS(1 .5x)-FS(x)

 

 

Figure 4.16. FresnelS(1.5:1:) — FresnelS(x) function graph

FresnelC(x) and FresnelS(x) oscillate around i as a: —+ 00 with slowly decreas-
ing amplitude. Hence, the difference of the two FresnelC or FresnelS functions
oscillates around zero, provided the arguments are proportional, which is the case
in Eqs. (4.2.10). The graphs of the functions HesnelS(1.53:) — FresnelS(x) and
FresnelC(1.5:c) — FresnelC(x) are shown in Fig. 4.16 and Fig. 4.17, respectively.

To calculate the integrals (4.2.6—4.2.7) in the general case, numerical integra-
tion methods should be employed. For this study the adaptive Simpson quadrature
method [39] is used. The main issue with the sector approximation is that it is not
sharp enough, and only works for the beams which are displaced by the transversal
kick in such a way that the whole beam is away from the origin. To handle the
situation with the beam that crosses the origin, and to be more precise with the

conclusions about the centroid, the attention should be paid to the exact shape and

83

Properties of Fresnel Functions

4

----- FC(x)
l' - - - FC(1.5x)
, — FC(1.5x)-_-FC(x)

 

   
   
   
  

 

 

 

 

999
”-505

O

 

 

I
P
M

 

 

 

,FC(1 .5x),FC(1 .5x)-FC(x)

 

 

Figure 4.17. FresnelC(1.5:1:) — FresnelC(x) function graph

position of the beam after the kick.

4.2.4 Elliptical Beam, Uniform Distribution

Let us assume that the particles in the beam are distributed uniformly (for simplicity,
we consider the general case later), and the beam has an elliptical shape. Then after
the transformation to the normal form coordinates the beam has the elliptical shape
again, and the axes of the transversal section of the beam are equal. Then the
boundary curve for the beam in the normal form coordinate pair is the circle, and
the parametric representation for it can be found in the form of the equations for two
half-circles: (r, 901(r)), (r, <p2(r)). Without loss of generality it can be assumed that
the resulting circle has its center on the horizontal axis, with the coordinates (61,0),

where d > 0 is known. Let p be the radius of the beam, then the beam lies between

84

 

 

 

 

 

Figure 4.18. R1 = d — p > 0

R1 = d — p and R2 = d + p, where it happens often that the radius R1 is less than
zero, which means that the origin (0,0) gets inside the beam (Fig. 4.18—4.19). Both
d and p parameters can be found by applying the linear normal form transformation
to the displaced beam boundaries. Below it will be shown that the two cases R1 > 0
and R1 < 0 can be treated in a uniform way. For the moment let us assume that
R1 > 0.
Similar to the previous section, the centroid of the beam has the coordinates
90(7‘)+9N(T)T
S=/ RQ/ dddr;
1R1 <p(r)+0N(r)
<P(T) +9N(T)
r(N)= 122/ r2 cos 9d0dr; (4.2.11)
—3 R1 «7(0) )+9N (T)

+9 N (7‘)
ycN - 32/—<P(7“) r2 sin Odddr;
S R1M(T)+9N(

the only difference being that (,0 is now a function of r. This double integral reduces

85

 

 

Figure 4.19. R1 = d - p < 0

to a simpler form, which is derived below, but the resulting single integral can only
be found using numerical methods. The integrand is not simplified further here, as
the non-uniform density beam case will be considered in the next subsection, which
only makes the integrand more complicated, thus not allowing any simplification of
the general form.

To simplify the integral expression, some additional information is needed on the
intersection of the two circles, as we are integrating along the arcs of a circle and the
boundary curve is also a circle. Let us consider the two circles: the first one centered
at the origin (0, 0) and having a radius r, and the second one centered at (d, 0) and
having a radius p (Fig. 4.20). This setup gives

$2 + y2 = r
(4.2.12)

(x-d)2+y2= p

86

 

 

 

 

Figure 4.20. The intersection of two circles

Solving the system for as one obtains

_r2—p2+d2
_ 2d ’

 

which yields

:1:_7‘2—p2+d2
'r_ 2dr

2 2 2
_ r -p +d
<p(r) — arccos( 2dr ).

9

Using the resulting expression for (p and Eq. (4.2.8), similar to what is done in

87

the previous subsection one gets

<p(‘r)+9N (1')
/RR 2/<p r2 cos 0d0dr— -
<p(1‘) +9N(7‘)

R2
2/ r2 sin<p(r )cosQN(r)dr =
31

 

 

 

and hence

 

 

 

 

 

 

= 2/Rl 7‘ sin arccos 2( 2dr ) cos 0N(7‘)dr =
R2 — p2 + d2 2
7.2
=2/R1 \/1—( 2dr ) cos 6N(r)dr,
R2 2 __ 2 2
)=/ 7'\/1—( 'r '0 2+d2) dr;
R1 2dr
R2 — 2 +d2 2
< (N 7.2 p - 4.2.14
:cc :8 R1 \/1—(: 2dr ) cost(r)dr, ( )
R _ 2 d2 2
(N)_ 2 2 /1_ p + . 6
\ yc —S R1 2dr ) sm N('r)dr.

 

In the special case of R1 = d — p < 0 with the layout corresponding to Fig. 4.19,

there are two choices:

0 For 0 < r < [R1|, the intersection of the two circles in (4.2.12) is purely imag-
inary, and hence the whole centour (1‘, go 6 [—7r, 7r)) belongs to the beam, and
one can assume for such 7‘ that go goes from -7r to 7r. This is the approach used

later for the non-uniform distribution;

As the set {(r,cp)] r E (0, [R1|),90 6 [—7r, 71)} is symmetrical with respect to
the origin, and the particles in the beam are distributed uniformly, the centroid
of this part of the beam is (0,0), so it can be omitted, and all of the above
formulas are applicable with no change if we put R1 = |d — p]. Unfortunately,

this approach does not work for the non-uniformly distributed beams.

88

4.2.5 Elliptical Beam, Normal or Arbitrary Distribution

Assume that the beam distribution is normal in both directions in every pair of
coordinates, and each two directions are independent. As the beam is round in
the normal form coordinates, the variances in both eigen—directions are the same
a = 09; = 0y, and hence the resulting density of the bivariate distribution is defined

by the formula

 

 

_ _ 2 2
f(z,y) = 27:02 exp( ((2: d) +y )) , (4.2.15)

as the mean values for the distribution are d and 0. Note that this formula is only
valid for the initial distribution, when 0N = O, and after N turns 0 N should be
subtracted from the value of the angle.
A typical particle distribution after various numbers of turns is shown in Fig. 4.21.
In the case of the non-uniform distribution the expressions for S, me, and ye are
essentially the same as in Eqs. (4.2.14), except that now the integrand is complicated
by the additional factor of f (r cos(6 — 6N), T sin(6 — 6N)). The componen “—6N” 18

introduced to always take the density of the initial normal distribution:
cp(r)+9N(r) .
NS=/RR 2/(p f(rcos(6 — 0N), rsin(6 — 0N))d9dr;
<p(r)+91v (T)

90(r)+9N(7“)
R2] 7‘ 2c—036f(7‘cos(9 0N), rsin(9— 0N))d0dr; (4.2.16)
:5 R1 <p(r)+0N(r)

<P(7')+9N(7')
—RRZ/(p r2 sin 6f(r cos(0— 6N), 'rsin(6— 6N))d0d'r.
<p(r)+0N(r)

In the case of negative R1, the part of the beam inside the circle of radius |R1|
cannot be omitted, and hence, the first of the two methods described at the end of the

previous subsection is used and the integration goes from 0 to R2 taking cp 6 [—7r, 7r)

for O < r < |R1|.

89

 

 

 

 

 

 

 

 

 

 

 

 

(b) after 1000 turns

 

 

 

 

 

 

 

 

 

(c) after 3000 turns

Figure 4.21. Particle distribution in the beam, the cross indicates the center of mass
of the beam.

90

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

5
x10_4 x1o‘4

(8) initial distribution (b) after 1000 turns

 

 

 

 

 

 

 

 

 

 

 

x110—4

(c) after 3000 turns

Figure 4.21. Particle distribution in the beam, the cross indicates the center of mass
of the beam.

90

4.2.6 Numerical Experiment Results

The calculation method described above allows one to find the dependence 7‘ =
r(N, cl, oz) for the elliptical beams with an arbitrary particle distribution, the only
requirement being that the initial distribution density function is known. Hence, for
each pair of values c1 and 02, one can find 77(c1, c2) or n1/2(c1, c2) introduced in the
beginning of the section as well as the corresponding values of the functions (4.2.2)
and (4.2.3). Having these data available and employing various optimization meth-
ods, one can find the correct values of c1 based on one particular measurement or
both coefficients cl and c2, provided that measurements for different kick strengths
are available.

The numerical results for the Tevatron correspond to the values obtained by track-
ing the nonlinear model of the machine. Consider, for example, the function J5, given
by Eq. (4.2.3) and the results of the measurements shown in Fig. 4.10. The number of
turns after which the amplitude of the center of mass falls down to a half of its value
varies depending on the BPM, one of the total of 115 reliable measurements. Tak-
ing the average over all the BPMs one obtains that N1/2 z 1000, or more precisely,
N1/2 6 [900,1100].

An optimization procedure results in the expected value of c1 = —2511 for the
initial beam amplitude after the kick of 7‘ = 0.24 - 10‘3. Taking into account that

[40 = 0.585, one gets
u z #0 + 617‘2 = 0.585 — 1.4463 . 10—4. (4.2.17)

To conceive how close the obtained value of cl is to the realistic value of the
tune shift with amplitude, a comparison was performed in COSY using the nonlinear

model of the Tevatron available at the official lattice page at Fermilab [38]. The model

91

 

 

 

 

 

 

 

 

 

 

x 10 ' r . .
—Calculated amplitude
a, - - -Modeled amplitude
13 i
1&2. ............................................................................................................................ .
E
an
E
.2
Tu . .............................................................................................................................. .
E 1
o
z
o 2 s ‘ M
0 1000 2000 3000 4000 5000
Number of turns

Figure 4.22. Calculation results and the comparison with the nonlinear model

reflects changes made to the accelerator as of November 2005, while the measurement
results are dated January 2006. Hence, there are strong reasons to assume that the
model distribution tracking should yield results comparable to the calculated value
of the tune shift, and therefore, the measurements. The COSY calculation shows
that the expected value of c1 for the nonlinear model should be -—2541, which means
the calculated value found by the optimization differs from the model value by not
more than 2%. At the same time, only the information about the distribution of the
particles in the beam, the size of the beam, and the linear dynamics was used to find
the nonlinear tune shift. Necessary additional information was extracted from the
measurements.

Figure 4.22 shows the graphs of the calculated amplitude with c1 = —2511 and

92

the model amplitude with cl = —2541. The slight difference between the graphs can
be explained not only by using different c1 ’8, but also by the fact that the fourth order
term c2r4 in the expansion of it has not been taken into account. Also the model is
just an approximation to the real machine’s optics. At the same time, the similarity
of the graphs allows to conclude that the model represents the real machine quite
accurately, at least for the low order nonlinearities affecting the tune shift (mainly
the sextupole content of the ring).

Also, the validity of the approach studied is perfectly supported by the indepen—
dent calculations done years ago. There is an estimate of the nonlinear tune shift by

R.Meller et al. [40], given by the following formula:

 

M R” (1’0 — ”A27
~ 1 (4.2.18)
47TN1/2,

where A is the amplitude of the center of mass of the beam, measured in 0 units
of the beam under consideration. This formula is derived for the beams with a
normal distribution of the particles, and it represents a good approximation when
the transversal kick is relatively weak (A is not too much greater then 1).

The value comparison of EA2 from Meller’s article to the value of cl'rz, obtained
by the calculation using the algorithm described above, gives the following results:

n = 7.96 - 10-5, A = 1.36,
)1 z #0 + 5A2 = 0.585 — 1.4723 - 10—4, (4.2.19)

that is, the difference between the values obtained using different approximations in

Eqs. (4.2.17) and (4.2.19) is less than 2%.

93

4.2.7 Summary

Correspondence is found between the first and the most important term in the ex-
pansion of the nonlinear tune shift with amplitude and the BPM measurement results
after the beam is kicked transversely. This correspondence can be found by minimiz-
ing the function (4.2.2) or (4.2.3). To be able to find the values of either one of the
two functions it is necessary to have the information on the behavior of the amplitude
7' of the center of mass. This amplitude can be found in the most general case using
the formulas (4.2.16) and 7‘ = W.

The COSY INFINITY’S normal form transformation algorithm and built-
in optimizers were employed to calculate the nonlinear tune shift with ampli-
tude along with the numerical integration algorithm and visualization tools of
Mathworks Matlab R13.

A method for the calculation of the nonlinear tune shift with amplitude was tested
on the Tevatron BPM measurement results and compared to the nonlinear model
calculations as well as the independent approximation formulas. In both cases the
discrepancy was less than 2%, which can be considered a very good result taking
into account that only the information on the one-turn linear transfer map and the
geometry of the beam has been used, while the lack of information on the nonlin-
ear behavior was compensated by a single BPM measurement with one particular
perturbation (kick) strength.

After the coefficient cl has been found, one might try finding the coefficient c2
if multiple measurement results are available. On the other hand, in the case of the

2

Tevatron, c27‘4 is 1—2 orders smaller than c1?" , so in this particular study there was

no attempt made to find c2.

94

5. MUON COLLIDER STORAGE
RINGS: DYNAMICS AND

OPTIMIZATION

5.1 50x50 GeV Collider IR High-Order

Correction

The progress in the development of the Interaction Regions (IRS) for advanced future
colliders, in particular various Muon Collider designs, requires a systematical high—
order dynamics study. The effects from the high-order beam dynamics become more
and more pronounced as the value of the fl-function at the Interaction Point decreases
to very small values, fl*=14 cm and fi*=4 cm for the 50x50 GeV Muon Collider and
Higgs Factory, respectively, and fi*=1 cm for the proposed 750x750 GeV Muon Col-
lider that matches the geometry of the Tevatron tunnel. This last lattice is described
in detail in Section 5.2. Simultaneous control of geometric and chromatic aberrations
is critical to the success of future machines and can only be achieved through the
deliberate addition of the nonlinear fields in the Interaction Region itself. Therefore,
there is a need of studying both the correction schemes and the unavoidable impact

of the high-order correctors — sextupoles, octupoles and even duodecapoles — located

95

in the Interaction Region close to the low-beta quadrupoles.

Muon colliders have certain advantages over high-energy proton and electron col-
liders. The size of a muon collider is significantly smaller in the circumference than
that of a proton collider since the rest energy of a muon is much less than that of a
proton. When compared to electrons, however, the larger mass represents an advan-
tage by reducing the synchrotron radiation and beamstrahlung, allowing for higher
magnetic confinement, and reducing the size of the collider relative to an electron
collider. On the other hand, muons are not stable and they decay relatively quickly
(7 = 2.2 us), so the acceleration must be rapid and the circumference of the ring must
be as small as possible to minimize the luminosity degradation due to the decay. Also
the muon decay is a source of secondary particles with relatively high energies —— elec-
trons and *7 particles. These, in turn, result in a powerful energy flow, which could
even destroy some important expensive components, especially the superconducting
magnets, or spoil the measurement results.

The muon accelerating complex consists of a number of stages: Muon Production,
Cooling, Acceleration, and finally the Collider Ring [5,7]. The collider ring circulates
beams of 12+ and ,u— traveling in the opposite directions. Considerable design effort
has been invested in a 50x50 GeV Muon Collider [28,55] and a 2X2 TeV one [27, 56].
To reach the luminosity goal of 1035 cm"2 8—1, the beta-function of the beam at the
interaction point should be extremely small for the 2 TeV Collider: [3* = 3 mm. The
collision spot is even smaller for the Linear Collider.

It is therefore useful to begin the study at the lower collision energy and less re-
strictive beam sizes. For the 50x50 GeV Higgs collider the value of the beta-function
at the interaction should be 6* = 4 cm for a momentum spread of dp/ p z 0.12%, rms.

Many issues associated with the design of such a collider are summarized in [55].

96

Specifically, the tune shift with amplitude is a very important figure of merit to
consider, especially the second order terms, which remain large in various proposed
schemes [28,55].

This study addresses the correction schemes and the effects of the high-order
correctors located in the IR itself proximate to the low-beta quadrupoles, or focusing
elements, to control the kinematic and geometric aberrations, the latter of which arise
from the use of sextupoles for the chromatic correction. The dynamic aperture of the
storage ring is the optimization criterion for the study. Finally, an increased dynamic

aperture is indicative of the potential further reductions of the value of 6*.

5. 1. 1 Dynamic Aperture

Dynamic aperture studies to date of the Muon Collider storage rings can be found
in [28]. In the following sections the dynamic aperture study and optimization are
considered more systematically. As the collider storage ring is a circular accelerator
where the muons stay for about 1000 turns, it is convenient to find the nonlinear
transfer map of one full revolution and then apply it repetitively to study the stability
of a particle distribution. According to [12], the study of the transfer map is sufficient
for characterizing the stability of the design optics. The main tool for the study is
COSY INFINITY [13] with its high order Taylor expansions of the flows of differential
equations and differential algebraic tools to determine flows of differential equations.
COSY allows such complicated objects as the nonlinear resonance strengths and the
nonlinear amplitude dependent tune shifts to be computed and optimized, or the
individual terms comprising these values to be addressed, all of which are important
for a careful demonstration of the optimization.

The dynamics of particles is described by a system of differential equations of the

97

Figure 5.1. The grid of points to determine the dynamic aperture

form (4.1.1) with the same restrictions as described in Section 4.1. The energy of the
particles does not change as the muons are already at the desired energy. The main
idea of the optimization is to minimize the particle losses for the maximal beam size.
If at some point the coordinate of the particle is greater than the beam pipe size, this
particle is considered lost. The number of the particles lost serves as a quantitative
criterion of the practical stability.

Hereafter, the coordinate (in the (2:,a) or (y, b) plane) of the particle farthest
from the origin (from the trajectory of the reference particle), which remains stable
for the 1000-turn storage cycle, is called the dynamic aperture. The dynamic aperture
can be found using the method illustrated in Fig. 5.1. A grid of particles with the
initial positions forming the concentric circles around the closed orbit (fixed point) is
launched and tracked for a thousand turns. The radii of the circles increase in fixed
steps until the trajectory is unstable, that is, until some of the outermost particles are
lost, numerically. The largest radius for which all the particles remain in the beam
pipe is considered to be the design dynamic aperture. The incremental step in the
radius is taken to be in units of one sigma (a) of the standard normal distribution of

the particles with the normalized rms emittance of 9077 mm mrad (a = 82 - 10'6 m).

98

The beam is considered to be round. From an optics standpoint, elliptical beams are
more manageable and less nonlinear than round beams in the design of Interaction
Regions. Using a [3* ratio of 1:4 of the horizontal size to the vertical size (factor of 2
in the relative beam sizes), however, causes a decrease in the luminosity of a factor of
2, which is unacceptable. Therefore, the condition of round beams at the Interaction
Point (IP) has been imposed in all current collider designs.

The optimization of the dynamic aperture by adding the high-order nonlinear
correctors such as octupoles, decapoles, or duodecapoles to the superconducting IR
quadrupoles and selecting the proper strengths for these correctors can be done in
various ways targeting different objective functions and nonlinear terms. The initial

dynamic aperture before any correction is introduced is 70 (see Fig. 5.2(a)).

5.1.2 Minimization of the Nonlinearities

The discussion about the impact of the nonlinearities on the dynamics, and their
correction in the Tevatron in the beginning of Section 4.1, stays valid for the Muon
Collider. By interacting in a complicated way, the nonlinear components of the dy-
namics may significantly reduce the stability of the machine and can even lead to
chaos. While in the linear case all the particles follow the closed orbits and stay
stable forever, so the dynamic aperture is always bigger than the beam pipe size, in
the nonlinear case this does not happen. Therefore, one of the ways to formulate a
problem of maximization of the dynamic aperture is this: minimize all the important
nonlinearities of at least lower orders (second and third) to make the dynamics as
close to linear as possible. In terms of the Taylor expansions of the transfer maps and

the corresponding equivalence classes [M]n, the objective function for the transversal

99

 

dynamics has the form:
3 . . . .
1. Z Z 2
J5 = E E [(zklzl1 222233z44) , (5.1.1)

where 2(1—1') = (z1 (21'), 22 (71'), 23 (11'), 24 (17))T is the vector of the phase space coordinates
of (4.1.1), 21' is the vector of the control parameters having the strengths of the high-
order correctors in the lattice as its components, (zklzil 2322:3214) is the aberration
of order 7:1 + 2'2 + 2'3 + 2'4, that is, the corresponding coefficient in the Taylor expansion
of the phase space variable zk. All the correctors are located in the neighborhood of
the interaction point, and can be effectively thought of as additional superconducting

coils in the superconducting quadrupoles of the final focus telescope. The number of

parameters may vary, but the general idea is to keep it as small as possible.

5.1.3 Minimization of the Resonances

Another systematic approach to the problem is to find a corrector scheme which
addresses the nonlinear resonances and/ or minimized tune shifts with amplitude.
Instead of minimizing specific nonlinear terms in the map, all the resonances or tune
shifts up to a particular order are minimized. The strength of all the resonances can be
found as a by-product of the application of the normal form transformation algorithm
described in Chapter 3. In this case the objective function can be represented as

4 2

J7 = Z Z (nglk’fl 6)]:

m=23=1

 

v + — s _.
k +k --o +_ — : .
(Hrll l)ez,u(k k )_,rjeiw] <K, ,
l=1

(5.1.2)
where n is some threshold characterizing the proximity to the resonance. J7 depends
on the control parameters 11' through the components Tnfj connected to the elements of
the transfer map 8% j after the corresponding stage of the normal form transformation.

The form of J7 is closely related to Eq. (3.1.10): the smaller the denominator is, the

100

harder it is to remove the unwanted Tnfj term. As the strengths {1' of the correctors
change, the form of the transfer map M, and hence, Sm ‘v’m = 2:5 change as well.
The dependence of J7 on the control parameters is continuous, as it follows from the

argument chain above.

5.1.4 Dynamic Aperture as the Objective Function

One more way to attack the problem of the storage ring performance optimization is

to use the dynamic aperture itself as the objective function:

J8 = {ma‘X(TI + T2)[Mk(xi7ai7yi7bi) = (13f,af,yf,bf), [xfl < 77: nyl < 777
i=1,N,k = 1,1000, V2222 + a? = 7'1, V3122 + b? = 72}, (5.1.3)

where k is the number of turns to keep particles confined by 77 in both coordinates
xf and yf, N is the number of particles in the grid which is used to determine the
dynamic aperture, and 2',- = (mi, 61,, yi, b,-)T is the phase space vector for the particle
lying in the plane (2:, a) on a circle of radius 71, and in the plane (y, b) on a circle of

radius 72. For the uniform distribution of particles along the circle one obtains
(“37107) = 7‘1 swim). 90j(N) = 277/N; (5.1.4)

(47,6) = r2 - e774”), MN) = 2.1/N, (5.1.5)

which corresponds to the distribution of particles shown in Fig. 5.1.
Again, the transfer map M depends on the strengths of the correctors 71', added
to the lattice for optimization purposes, and hence, J8 = J8(7'l). The problem is to

find a configuration of the corrector magnets which maximizes the value of J8.

101

5.1.5 Advantages and Disadvantages of Different

Optimization Approaches

The first approach to the optimization is to minimize some (or all, if possible)
nonlinear coefficients in the transfer map of one full revolution up to possibly higher
orders. The objective function J5 from Section 5.1.2 corresponds to this case.

Calculations show that this approach is not very effective due to the fact that the
number of the multipole correctors is very limited, while the number of the nonlinear
coefficients in the transfer map is large. For example, for the third order the map
of the storage ring lattice under consideration has 20 coefficients total, 10 in each
plane. At the same time, the nonlinear coefficients are always connected to each
other, so technically the number of parameters might be smaller than 20, but this
very connection is the source of the problems. By minimizing specific nonlinearities,
others are introduced of the same or higher order, and these often have an equally
destructive impact on the behavior of the system.

This approach has been attempted and improved the dynamic aperture from 70
to 100 (Fig. 5.2), which is not a very significant result compared to the other two
approaches.

The second approach (more systematic) is to find a correction scheme which
addresses nonlinear resonances and / or minimizes tune shifts with amplitude. Instead
of minimizing specific nonlinear terms in the map, all the resonances or tune shifts
up to a particular order are minimized. The objective function J7 from Section 5.1.3
corresponds to this approach. Unfortunately, the number of correctors available (three
locations at most) does not allow for the removal of all the desired terms of order three

or higher. Hence, the situation is similar to the previous approach: in minimizing

102

 

 

(a) before optimization

 

 

(b) after optimization

Figure 5.2. The comparison of the phase space portraits before and after the opti-
mization using the first approach (minimization of individual high order coefficients)
in the (x, a) plane. Particles are launched along the :1: axis with an increment of one

a, a = 82- 10"6 m.

103

 

one term, others increase and cause comparable stability problems. Still, among
different possible configurations there is one which requires three octupoles and one
duodecapole and increases the dynamic aperture from 70 to as high as 120. In
this configuration the following high order coefficients in the amplitude—dependent
tune shifts have been removed (2D simulation): (xlxz), (ylzr2), (xlyz), (yly2), (mlx4),
(yly4). There are no tune shifts of the third and fifth order, so the only terms up
to the fifth order that remain in the tune shifts with amplitude are (xly4), (ylz4),
(xlx2y2), (ylxzyz). The phase portraits in Fig. 5.3 show that the second approach
works significantly better than the first one.

The third approach — to take the dynamic aperture itself as the figure of merit
— has been found to be the most effective way to maximize the dynamic aperture,
and is at the same time the most straightforward. The objective function for this
approach is J3 from Section 5.1.4. Of course this approach leads to the increased
volume of calculations inside the optimization loop as one needs to recalculate the
dynamic aperture for each new value of the corrector strengths, which is more time
consuming than the calculation of the resonances or tune shifts because particles are
tracked for 1000 turns in each optimization step. This is not a serious limitation for a
7th order calculation though, and the 7th order is sufficient to determine the dynamic
aperture. The situation here is to a large extent similar to that of the Tevatron,
see Section 4.1.10. This approach is capable of finding a corrector configuration
that achieves a dynamic aperture of 130 using only two octupole correctors and no
duodecapoles at all. It is possible that 130 is the upper limit of the dynamic aperture
for the particular storage ring under consideration.

The results of the third approach are shown in Fig. 5.4. The two other approaches

have not been as successful and require further investigation. The problem with these

104

 

 

(a) before optimization

 

(b) after optimization

Figure 5.3. The comparison of the phase space portraits before and after the opti-
mization using the second approach (minimization of the resonance strengths) in the
(2:, a) plane. Particles are launched along the .7: axis with an increment of one a,
a = 82 ~ 10‘6 m.

105

 

nonlinear terms resonances dyn. aperture
minimization minimization max1m1zation

 

Dyn. aperture improvement 70 => 100 70 => 120 70 2. 130

Number of correctors 3 octupoles+ 3 octupoles+ 2 octupoles
2 duodecapoles 1 duodecapole

 

 

 

 

Table 5.1. Parameter comparison for various optimization approaches

approaches is the large number of values to minimize simultaneously and the absence
of a well established algorithm for the search of the connection between different high
order nonlinearities.

All the approaches attempted are summarized in Table 5.1.

5.1.6 Optimization Results

The tracking pictures showing the phase space portraits of the system after some
number of turns are, in general, an effective tool to serve as a qualitative estimate for
the dynamic aperture. To check that the third approach is efficient in general rather
than for a particular choice of the initial coordinates of the particles being tracked,
the Fig. 5.5—5.6 demonstrate the phase portrait comparison for the particles launched
along the diagonal (2:, y). Horn Fig. 5.5—5.6 it is clear that the dynamic aperture after
the optimization is significantly larger than before the optimization, even though it
is hard to estimate its exact value using phase portraits only. Therefore, for each

approach the quantitative results are given explicitly.

106

 

 

 

(b) (y, 5) plane

Figure 5.4. Phase portraits after the optimization in the (23,67) and (y, b) planes.
Particles are launched along the :17 and y axes, respectively, with an increment of one
a, a = 82 - 10"6 m.

107

F

 

(a) before the optimization

 

(b) after the optimization

Figure 5.5. The comparison of the phase portraits in the (2:, a) plane before and after
the optimization. Particles are launched along the (13,31) diagonal with a step of one
a.

108

 

(a) before the optimization

 

(b) after the optimization

Figure 5.6. The comparison of the phase portraits in the (y, b) plane before and after
the optimization. Particles are launched along the (:3, y) diagonal with a step of one
a, a = 82- 10-6 m.

109

 

 

 

 

 

1500.......I ............ , ................ ,. ......... . ....... ............ .............. -- B ....-q
3 g 2 5 i j 1

 
 
 

>~1000

'.

6,.

 

500.. ......... ,.-, .. .............. ; ................. . ., .......... . .....

 

 

 

0 20 40 60 80 100 120 140
Longitudinal coordlnate, z [m]

Figure 5.7. The octupole correctors layout
5.1.7 Summary

A variety of approaches aimed at increasing the dynamic aperture was studied, and
the one which produced the most effective way and produced the best results was
determined. The results of the optimization using the chosen method developed for
the study were presented and the dynamic aperture was found to increase almost by
a factor of 2 compared to the situation before the optimization: 7a to 130. Only two
octupole correctors were required in the proximity of the interaction point (at the
peaks of the beta functions in the a: and y planes). Though the aim was to maximize
the dynamic aperture, some of the higher order terms in the nonlinear resonance
matrix decreased naturally by an order or two in magnitude. Their reduction serves
as another strong, quantitative argument proving the efficiency of this optimization

approach. The next logical step in the study might be to add higher order multipoles

110

 

(decapoles, duodecapoles), but a preliminary study showed them to be ineffective in
the sense of improving the dynamic aperture further. For example, different config-
urations could be obtained with duodecapole and octupole fields which resulted in
particles being stable at 140, while at the same time the particles at 130 were lost. It
appears that the correction of geometric aberrations and kinematic effects are almost

completely represented by the third order terms.

5.2 750><7 50 GeV Collider Storage Ring Lattice

for the Tevatron Main Ring

A new lattice for the Muon Collider storage ring with design collision energy of
750x? 50 GeV is discussed. The important building blocks of the lattice are the Fi-
nal Focus Section, the Chromatic Correction Section, and the Arc Module, which
are described below in detail. These components of the collider have been designed
keeping in mind that the storage ring must approximately match the footprint of the
Tevatron Ring in order to take advantage of existing services and tunnels. The model
presented here relies heavily upon a previous, very efficient 50x50 GeV storage ring
lattice design. The current design value for 6* (beta function at the collision point)
is chosen to be 1 cm, which is less restrictive and has the advantage of lower chro-
maticities and longer bunch lengths (due to the hour-glass effect) than the previous

standard lattice with a 6* of 3 mm.

5.2. 1 Introduction

The idea is to design a lattice for the storage ring that fits or matches approximately

the footprint of the Tevatron Main Ring tunnel with all its bends and straights.

111

The alternative idea might be to make another tunnel, which is slightly larger in the
radius, around the currently available one to be able to use the current Tevatron
tunnel infrastructure. Taking into account the current status of the Tevatron project,
the Muon Collider might be a logical next step in utilizing the existing tunnel with its
infrastructure, thus saving a large amount of expenses connected to building a new

accelerator complex for muons.

5.2.2 750x750 GeV Lattice

Currently a 50% dipole packing fraction is used for the 750x750 GeV lattice, which
results in an arclength of 5.85 km (6283—0432 km of straights) in the dipole field of
5.3 T. At 750 GeV this field strength is reasonable, and in fact, the ultimate energy
might be increased to 1x1 TeV.

The 50x 50 GeV lattice [29,55] is used as a baseline, and its components are scaled
to handle the muons the the energy of 750 GeV. The layouts of both 50x50 GeV
and 750x 750 GeV lattices are shown in Fig. 5.8. The 50x50 GeV lattice is a highly
optimized one which, in turn, is based on the 2 X 2 TeV lattice [27], and therefore there
is a strong reason to assume the 750x750 GeV design shares most of the advantages
with these other lattices.

The 50x50 GeV ring has a roughly racetrack design with two circular arcs sepa-
rated by an experimental insertion on one side, and a utility insertion for injection,
extraction, and beam scraping on the other. The experimental insertion includes the
interaction region (IR) followed by a local chromatic correction section and a match-
ing section. The chromatic correction section is optimized to correct the ring’s linear
chromaticity, which is almost completely generated by the low beta quadrupoles in

the IR.

112

50x50 GeV 750x750 GeV

 

 

1 IR _l
' ‘ ARC ARC repeated
1R1
l 4
insertion
IR=IP+FF+CCS+MM

 

Figure 5.8. The baseline 50x50 GeV lattice scheme compared to the 750x750 GeV
lattice scheme.

The 750x 750 GeV lattice design uses the same building blocks as the 50x50 GeV
version, as it can be seen from Fig. 5.8: the final focusing section (FF), the chro-
maticity correction section (CCS), the matching module (MM), and the arc (ARC).
The difference in the layouts is dictated by the fact that the 750x750 GeV lattice
must match the Tevatron footprint. So the final focus section is placed in one of the
straight sections of the Tevatron ring, while the chromaticity correction section and
the matching section occupy parts of the Tevatron arc. The arcs of the original design
are repeated in the arcs of the Tevatron as many times as necessary.

There are two IRs in the proposed lattice design to compensate for the luminosity
loss due to the increased 6* as described in Section 5.2.6.

The main parameters of the 750x750 GeV lattice are summarized in Table 5.2
(second column). Similar parameters for the other lattices — 2x2 TeV (first column)

and 50x50 GeV (third column) — are shown for comparison.

113

 

 

 

4/1.5 TeV 1.5 TeV 100 GeV
(this design)
6* [mm] 3 10 40
l* (IP to quad) [m] 4 5.5 4.5
peak 6 [km] 145 35 1.4
IR quad aperture cm] 10 10 10
Poletip field [T 12 9 8
N(9)5% )[mm mrad] 8417r/3157r 13067r 21767r
Ap/p( 95%)é‘7o .01-.08 Z .018-.144 Z .036—.288
53,5112 + C -1500 -456 -53
a, IR + COS -2000 -645 -73
am 3.6 x 10—4 1.0 x 10-3 3.0 x 10—2
IR length [m] 1300 506 137
am -2.1 x 10—3 —-9.3 x 10-3 —9.5 x 10—2
Arc length [In] 187 70 31

 

 

 

 

Table 5.2. Parameter comparison for various storage rings lattices
5.2.3 Final Focus Section

The low beta function values at the IP are mainly produced by three strong su-
perconducting quadrupoles in the Final Focus Telescope with pole—tip fields of 9 T.
Because of significant, large-angle backgrounds from themuon decay, a background-
sweep dipole is included in the final focus telescope and placed near the IP to protect
the detector and the low-6 quadrupoles [30]. The bend starts at 35 meters, so the
FF section fits the Tevatron straight section footprint. The layout of the preliminary

Final Focus section design is shown in Fig. 5.9.

5.2.4 Chromatic Correction Section

A local chromatic correction of the Muon Collider interaction region is required to
achieve broad momentum acceptance. The Chromaticity Correction Section (CCS)
contains two pairs of sextupoles, one pair for each transverse plane, all placed at the

locations with high dispersion. The sextupoles of each pair are located at positions of

114

 

 

W

 

 

 

 

1P
Win32 version 8.51/15 01/02/07 00.01.09
2 40.---4-.-.-.....-1.2..
., 35 . fix 6 D in e
i- '7 31-0 q
,3 30. - :; ....... -0.9
r ‘ i: ; £08
A 25. ‘ ,‘I ' '
E , , _-0.7
n 20.3 5 , 73.;
15‘ . .......... 30,4
10.- : “~~.\ [~03
- : ‘ ...... 40.2
5 T ........... - 0.1
0.0 1 . .4 . . 4 . . . . . . ”0.0
0.0 40. 80. 120. 160.

Figure 5.9. The final focus section beta functions and dispersion plots

equal, high beta values in the plane (horizontal or vertical) whose chromaticity is to
be corrected, and very low beta waist in the other plane. Moreover, the two sextupoles
of each pair are separated by a betatron phase advance of near 77, and each sextupole
has a phase separation of (2n + 1)§ from the IP, where n is an integer. The result of
this arrangement is that the geometric aberrations of each sextupole is cancelled by
its companion while the chromaticity corrections add. The sextupoles of each pair
are centered about a minimum in the opposite plane (78min < 1), which provides a
chromatic correction with a minimal cross correlation between the planes. A further
advantage to locating the opposite planes minimum at the center of the sextupole, is
that this point is 35 away from, or “out of phase” with, the source of chromatic effects
in the final focus quadrupoles; that is, the plane not being chromatically corrected
is treated like the IP in terms of the phase to eliminate a second order chromatic

aberration generated by an “opposite-plane” sextupole. The repetitive symmetry

115

 

CCSMOD

 

 

 

 

 

 

 

 

Win32 version 8.51/15 01/02/07 00.01.09
«3000.-.-.fi--.-.+.-.-2.2..
é ~43; I5 ---D “-2.1 E-
m- 2500. - 1‘ 12.0 Q

- }1.9

2000. - - 1.8
1 31.7
1500. , _-I.6
. _-1.5
1000. - :1.4
. - 1.3
500. - E12
. , _-1.1
0.0 "-25 - . - . . - . - 1.0
0.0 10. 20. 30. 40. 50. 60. 70. 80.
s(m)

Figure 5.10. The chromaticity correction section beta functions and dispersion plots

and the fact that the transfer map of the section is unity implies that the important

aberration (xl66) vanishes as well. The layout of the CCS is shown in Fig. 5.10.

5.2.5 Arc Module

The Flexible Momentum Compaction module (Fig. 5.11) provides negative momen-
tum compaction values compensating for the positive momentum compaction gen-
erated by the Chromaticity Correction Section. Small beta functions are achieved
through the use of a doublet focusing structure which produces a low beta simultane-
ously in both planes. At the dual minimum, a strong focusing quadrupole is placed to
control the derivative of the dispersion with little impact on the beta functions. (The
center defocusing quadrupole is used only to clip the point of the highest dispersion.)
Ultimately a dispersion derivative can be generated which is negative enough to drive

the dispersion negative through the doublet and the intervening waist.

116

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

H1666

 

     
      

 

 

 

 

 

NORMAL_A
A .110. Win33'vefsi0'n 81.51/15: r ' 01/02/0/ 00.01.09 0.8 A
E 99 j- Bf". B! L)" E
:0. ' .13.. 3' :06 Q
88- ' - 0.4
77". ‘6' "-0.2
66'l 1'" -0.0
44. - .' , f ‘ '
33. - : T“
22. 16' _--0.6
-/ ,,,,, : --0.8
0.0' - . . . . - . - . - 31.0
0. 0 10. 20. 30. 40. 50. 60. 70. 80.

s (m)
Figure 5.11. The arc module beta functions and dispersion plots

5.2.6 Advantages and Disadvantages of the New Design

The proposed design has a lot of advantages. As the lattice is isochronous, a bunch
length change is prevented, which is very important for controlling the hour-glass
effect. 0* is chosen to be 1 cm, which has the advantage of lower chromaticities and
longer bunch lengths (due to the hour-glass effect), and also the apertures can be
chosen smaller than the 3 mm lattice ones. Smaller chromaticities lead to weaker
chromatic aberrations and larger momentum acceptance. All these facts contribute
to a larger dynamic aperture.

As for the disadvantages, the choice of larger 0* leads to an undesirable decrease

in luminosity. According to the formulas from [19,20]:
1
5;)

where L is the luminosity. For 0* = 3 mm the hour-glass reduction factor is n A =

Loc

117

0.76, while the disruption enhancement is fl) = 1.5. Overall, H D = 77AfD = 1.14.

For 038w = 1 cm one has 77A -+1,fD —+ 1, and hence

it
(£011) = 1.1497164” = 3.8.
Lnew eff Old

Therefore, the luminosity for 0* = 1 cm is 3.8 times smaller than for 0* = 3 mm.
However, the loss of luminosity can be compensated by the increased momentum
aperture and by using the 2 IRs in the ring (see the scheme in Fig. 5.8).

One other problem arises due to the fact that one is trying to match the lattice
to the existing geometry, which puts more constraints on the building blocks of the
lattice.

The last and the most important for now is the one of shielding the IP and securing
the temperature of the superconducting magnets from the undesired effect of electrons
and ’y rays produced by the muon decay. However, this problem is beyond the scope

of this dissertation work and will be considered separately.

118

APPENDICES

119

 

APPENDIX A

TEVATRON LATTICE CONVERTER

FROM OPTIM TO COSY INFINITY

In this Appendix the detailed description of the program [45] converting the Tevatron
machine lattice source file from OptiM [32] code to COSY IN FITITY [13,36] is given.

Only the transverse dynamics changing elements are considered.

A.2 Introduction

The Tevatron lattice description in the form of the code list for OptiM — the program
written at Fermi National Accelerator Laboratory — is based on the current magnets
layout with several corrections added for a model to correspond the experimental
measurements. To analyze the transversal dynamics of the particles in the Tevatron
using the program COSY INFINITY, the code list of the lattice description needs
to be converted to COSY lattice elements description format. The situation is being
complicated by the fact that OptiM and COSY use different units of measure: in
OptiM all the units are SGS, while in COSY most of the units are SI, so great care

about it should be taken during the conversion of the lattice.

120

A.3 Tools for Conversion

The converter program is written in PHP language version 4.4.4 [1]. This program—
ming environment has been chosen to provide the availability of fast online conversion
of the source code. As PHP is a web-oriented language, the resulting program can be
run on any web server having PHP support.

Moreover PHP has a variety of standard functions to work with strings, which is

very convenient as the converter program is mainly a text file parser.

A.4 Source Code Details

The source code in OptiM format has several sections to parse:

1. A section with different variables declaration and value assignment, the so—called

“math header”;

2. A section containing values for initial energy and mass of the particles under

consideration, emittances, betatron functions and other parameters;

3. A lattice description section, consisting of a set of references to the lattice
element names, as well as the number of periods, telling the system to repeat

the same elements set several times;

4. A section describing elements, including all the lattice element names and pa-
rameters; the number of mandatory and optional parameters varies with the

type of the element.

Each line of the source file can contain a comment, starting with the “#” sign. Such

a line is ignored at runtime.

121

A.5 COSY Code Details in Connection to OptiM

Code

To convert the source code to COSY, all the sections should be treated and parsed
separately. As the OptiM program requires no explicit variable declaration (all the
variables are available at the moment they are first referenced), while COSY requires

one, the first section of the OptiM program splits into two sections:
1. Section with variable declaration;
2. Section with value assignment.

These two sections go to different places of the destination COSY code, namely, one
before the main procedure call and the other one inside it.

Additional problem arises as COSY is case-insensitive. So the two variables “el”
and “E1”, treated as different by OptiM, will be treated as the same variable by
COSY. To solve this problem the idea of a cache has been utilized. If the case-
insensitive name of the variable coincides with the name existing in the cache, this
variable gets another unique name.

The section, describing various parameters characterizing the beam, is partially
converted to the set of calls to standard COSY initialization procedures, such as OV,
RP, etc. Currently, only two parameters are taken into account: the initial beam
energy and the particles mass.

Element names and their order in the lattice from the subsequent section of the
source file go into a separate array. The names themselves are not used in the final
COSY code, but they are preserved as comments to each lattice element function call

along with the number of the element:

122

QUADRUPDLE 81.4578 qu_3 0;{qH_3F #1948}
DRIFT O.17178;{oDQDUT #1949}

DRIFT O.168275;{oDPIN5LOC #1950}

DRIFT 0.381;{0DHQUADC #1951}

M_ELEMENT 2 o o;{mrsx #1952}

ILELEMENT 1 0 0;{mTQX #1953}

Each element name contains the information about the element type. The first letter
of the name is this type (e.g. “q” stands for “quadrupole”, “o” — for “drift”, “s” —
for “sextupole” and so on). This information is used by the converter to determine
the function corresponding to the element.

The set of parameters of each element is being taken from the last section of the
source OptiM file. For the moment only the limited set of elements is implemented
(according to the needs of the study), namely, those influencing the transverse motion
of the particles (no RF cavities). At the same time, the support for these elements
can easily be added, as COSY has a very wide range of standard elements. The only
problem one might experience is the problem of units of measure conversion, which
is unavoidable.

In OptiM optional parameter values can be assigned explicitly or not assigned at
all. For COSY output, the values not explicitly assigned are substituted with some
default values and treated correctly in each individual case.

Each section of the source code can contain lines with comments, these lines are
treated individually for each section and included into the destination file to maintain

the completeness of the information contained in the source file.

A.6 Summary of Functionality of the Resulting

123

COSY Output

Currently, the converter program is intended for use with a Tevatron model only
and convert only the transversal dynamics. Some of the elements absent in the
description of the Tevatron lattice are not implemented, e.g. such element types as
’transfer matrix’ element. This element type allows to substitute an element or a set
of elements with a predefined linear transfer map in OptiM.

The initial aim of the program was to just convert OptiM source code to COSY
code, which can be run right away and return the transfer map of the Tevatron for
one turn (the only thing to specify is the calculation order for COSY). But with minor
changes this program can do much more. As it advances all the elements of the given
lattice one by one, it is easy to add some more special function calls for a predefined
set of types of elements. Another opportunity is to edit a code wrapper for some
element type, and have this change in all subsequent conversions of the source OptiM
code.

As the converter program is written in PHP, it is easy to perform updates to the
COSY code if the OptiM description is changed.

OptiM is aimed to assist with the linear optics design of particle accelerators,
while COSY INFINITY is commonly used for high-order analysis of the lattice, fast

high-order transfer map retrieval and paritcle tracking.

A.7 Code List of Element Wrappers

Each of the implemented structure element types has its own designated procedure
—— a wrapper preparing standard COSY procedure call. Code lists of all available

wrappers with parameters and unit conversions are summarized in this section.

124

 

7

1. Drift —— free space, element type ‘0’ or ‘i

PROCEDURE DRIFT L;
VARIABLE LO 1;
LO:=L;

DL LO;

endprocedure;

The procedure DRIFT has only one input parameter: the length of the drift.
The conversion of the units from [cm] to [m] is done by the converter program,
so technically the procedure DRIFT does nothing, but the wrapper DRIFT
might still be needed in case of extending functionality of the procedure. E.g.
in one of the variations the summation of the lengths of individual elements has

been added to each wrapper, so the procedure DRIFT changed:

PROCEDURE DRIFT L;

VARIABLE L0 1;

LO:=L;

DL LO;
current_length:=current_1ength+LO;

endprocedure;

2. Dipole — element of type ‘b’ or "d’:

PROCEDURE DIPOLE L B G ANGLEI ANGLE2 TILT;
VARIABLE BO 1;
VARIABLE LO 1;

VARIABLE GO 1;

125

VARIABLE TILTO 1;

VARIABLE NT 1;

VARIABLE KN 1 1;

VA

RIABLE Si 1 1;

VARIABLE $2 1 1;

VA

RIABLE R 1;

VARIABLE PHI 1;

NT:

LO:

IF

EL

=1;

=L/100;{in meters}

(B=O);

DL LO;

SEIF TRUE;

BO:=ABS(B)/10;{in Tesla}

GO:=G*10;{in T/m}
R:=CONS(PO)/(CLIGHT*1e-6)/BO;
KN(1):=GO/R/BO;
81(1):=ANGLE1*DEGRAD;{rad}
82(1):=ANGLE2*DEGRAD;{rad}
PHI:=LO/R/DEGRAD;

TILTO:=TILT;

ROTATE TILTO;{rotate coordinate frame}
MCLK LO PHI 0.05 KN SI 82 NT;{Standard COSY call}

ROTATE -TILTO;{rotate back}

ENDIF;

126

ENDPROCEDURE;

This procedure takes six parameters: the length of the element in [cm], the
bending magnetic field in [19G], the field gradient [kG/cm], the entrance and exit
angles in [deg] and the roll angle in [deg]. All the necessary unit conversions are
done inside the wrapper that ultimately calls a standard COSY function MCLK.
The procedure ROTATE allows to work with “skew” elements. It rotates the

coordinate frame through the given angle:

procedure rotate phi;
variable M 3000 4;
variable i 1;
M(1):=COS(PHI*DEGRAD)*MAP(1)+SIN(PHI*DEGRAD)*MAP(3);
M(3):=-SIN(PHI*DEGRAD)*MAP(1)+COS(PHI*DEGRAD)*MAP(3);
M(2):=COS(PHI*DEGRAD)*MAP(2)+SIN(PHI*DEGRAD)*MAP(4);
M(4):=-SIN(PHI*DEGRAD)*MAP(2)+COS(PHI*DEGRAD)*MAP(4);
LOOP i 1 4;

MAP(i):=M(i);
ENDLOOP;

endprocedure;

,

. Long quadrupole — element of type ‘q :

PROCEDURE QUADRUPDLE L G TILT;
VARIABLE L0 1;
VARIABLE GO 1;

VARIABLE TILTO 1;

127

 

LO:=L/100;{in meters}
IF (G=O);

DL LO;{if there is no field - just put drift}
ELSEIF (i=1);

GO:=G*10;{in Tesla/m}

TILTO:=TILT;

ROTATE TILTO;

MQ LO GO*0.05 0.05;

ROTATE -TILTO;

ENDIF;

 

ENDPROCEDURE;

This procedure takes three parameters: the length in [cm], the field gradient in

[kG/cm] and the roll angle in [deg].

. Multipole kick -— element of type ’m’:

PROCEDURE M_ELEMENT m B TILT;
VARIABLE BO 1;

VARIABLE TILTO 1;

IF (B#O);{if the field is non-zero}
IF (m=1);{quadrupole}
TILTO:=TILT;
ROTATE TILTO;

BO:=B/10;{in Tesla}

128

MQ 16-6 (BO/1E-6)*0.05 0.05;
DL -1e-6;
ROTATE -TILTO;

ENDIF;

IF (m=2);{sextupole}
TILTO:=TILT;
ROTATE TILTO;
BO:=B*10;{in Tesla/m}
MH 16-6 (BO/1E-6)*0.05*0.05 0.05;
DL -1e-6;
ROTATE -TILTO;

ENDIF;

IF (m=3);{octupole}
TILTO:=TILT;
ROTATE TILTO;
BO:=B*1000;{in Tesla/m}
M0 19-6 (BO/1E-6)*0.05*0.05*0.05 0.05;
DL -1e-6;
ROTATE -TILTO;

ENDIF;

ENDIF;

ENDPROCEDURE;

This is a rather complicated procedure compared to all others. It is designed
to facilitate various kinds of multipole kicks —— short multipoles that can be

treated as ‘zero—length’ elements. The procedure takes three parameters: the

129

number of multipole harmonic (1 — quadrupole, 2 — sextupole, etc.), the integral

of the multipole gradient in [kG/cmm”1] and the roll angle in [deg].

. Solenoid —— element type ’c’:

PROCEDURE SOLENOID L B A;
VARIABLE L0 1;

VARIABLE BO 1;

 

VARIABLE A0 1;

 

LO:=L/100;{in meters}

BO:=B/10;{in T}

 

IF (BO=O); I”
DL L0;{if no strength is given}

ELSEIF TRUE;
AO:=A/100;{in meters}
CMS BO A0 LO;

ENDIF;

ENDPROCEDURE;

This procedure takes three parameters: the length of the solenoid in [cm], the

 

strength of the field in [kG’] and the aperture in [cm].

130

APPENDIX B

GRAPHICAL USER INTERFACE FOR
THE TEVATRON SIMULATIONS IN

COSY INFINITY

COSY INFINITY does not have the tools to give the user graphical interface. There—
fore, to make the user interaction with the program for the Tevatron simulations more
efficient, a graphical user interface (GUI) was developed and implemented in Matlab
(see the screenshot of the main working window Fig. B.1). This GUI can be obtained
from the author of the dissertation upon request.

This GUI allows the user to run the Tevatron simulations without knowing the

COSY’s internal language. There are four groups of parameters the user can control:

0 Lattice options allow to select which families of elements are turned on/ off

while calculating the transfer map:

— Skew quadrupoles in dipoles on —— turns on / off the skew quadrupole errors

in dipoles;

— Skew quadrupole correctors on — turns on/ off the skew quadrupole cor-

rectors according to the currently implemented correction scheme;

131

 

hell

5

LU
.m
I

a
.1.
H
F
H
v1
S
O
C

n
d
d

0.
II.

n

n.

r
L
.md
‘1'
e
T

. to B.Erdelyl

acc

IS

 

ion GUI

imulat

Tevatron s

1.

igure B

F

132

— Weak sextupoles on — turns on/off the weak sextupole errors, mainly in

dipoles;
— Feeddown sextupoles on — turns on / off the feeddown sextupole families;

— Chromaticity sextupoles on — turns on / off the sextupole family responsi-

ble for the chromaticity correction;

— Solenoids (in detectors) on — turns on / off the solenoidal fields in the CDF

and DO detector areas;
- Separators on — turns on/off the sets of the separators which bring the

protons / antiprotons to their respective helices.

a Beta function options allow to choose whether the beta functions will be
calculated along with the transfer map:
— Calculate beta functions — turns on/ off the beta function calculation;

— Recalculate initial values —— turns on/ off the recalculation of the initial

values of the Twiss parameters using the periodicity conditions;

-— Current initial values ——- current values of the Twiss parameters, 73; and 73,
are calculated using the formulas 718,31 = (1 + 0:21;,y)/ 1627,21}, can be changed

to the desired values.

0 Miscellaneous options allow to choose the order of calculations and various

logging options:

— Calculation order — usually set to 7 for intermediate tracking results and

11 or 13 for the final results;

133

 

— Filename prefix — tells the program, which filename prefix is used for
.map and .log files, containing the transfer map and the current run log,

respectively;
— Save map — turns on/off the transfer map saving to file;

— Save beta initials — the initial values of the the Twiss parameters will be

saved as defaults for the next run;
— Fixed point on — turns on/off the fixed point calculation; F“

— Tunes on — turns on/ off the tunes calculation.

0 Tracking options allow to change the particle tracking settings:

 

— Sigma — the value of the normal distribution parameter a, the beam is L.

assumed round, so ox = 0y = a;

— Range — characterizes the minimal and maximal radii of the particles

being tracked in 0 units;
— Number of particles — how many particles are being tracked;

— Tracking calculation order — might only be lower than the transfer map

calculation order;
— Turns — the number of turns the particles are being tracked;

— Step — the number of turns after which the particle trace is saved, in other
words, the number of turns after which the dot on the phase space plot is

set;

— Particles distribution — the particles can be launched along the :3 axis, the

 

y axis or along the (0:, y) diagonal, other patterns can easily be arranged;

134

The main working buttons are the “Calculate map and other data”, in which case
only the transfer map and related information such as the beta function values are
calculated; and the “Track” button which triggers the transfer map recalculation if
the transfer map has not been calculated previously, followed by the particle tracking
procedure according to the parameters set. If any of these buttons is pressed Matlab
initiates the COSY code list assembly procedure, substituting all the updated pa-
rameters as necessary, after that COSY is called as an external command to run the
assembled code. As one can see the Matlab GUI gives the user a lot of freedom to
choose the transfer map calculation and tracking parameters.

This Matlab utility was used to obtain all the phase portraits in section 4.1.

135

BIBLIOGRAPHY

[1] PHP: Hypertext Preprocessor. Available at http://www.php.net.

[2] The MAD to COSY Converter. Available at http://cosy.pa.msu.edu/
’ converters/machosy.

[3] Methodical accelerator design. http://mad.home.cern. ch/mad/.

[4] M. Abramowitz and I. A. Stegun (Eds). Handbook of Mathematical Functions
with Formulas, Graphs, and Mathematical Tables, pages 300—302. New York:
Dover, 1972.

[5] S. Ahn and et al. Muon colliders: A scenario for the evolution of the Fermilab
accelerator complex. Technical Report Fermilab-FN-677, Fermilab, 1999.

[6] Y. Alexahin. Private communication.

 

[7] M. M. Alsharoa and et al. Recent progress in Neutrino Factory and Muon
Collider research within the Muon Collaboration. Phys. Rev. ST Accel. Beams,
62081001, 2003.

[8] M. Berz. Introduction to beam physics. Lecture Notes. Virtual Univer-
sity Beam Physics Course, available at http://bt.pa.msu.edu/pub/papers/
bphy861/bphy861.pdf.

[9] M. Berz. Differential algebraic formulation of normal form theory. In M. Berz,
S. Martin, and K. Ziegler, editors, Proc. Nonlinear Effects in Accelerators, pages
77—86. IOP Publishing, 1992.

[10] M. Berz. From Taylor Series to Taylor Models, pages 1—27. American Institute
of Physics CP405, 1997. Chapter in: Nonlinear Problems in Accelerator Physics.

[11] M. Berz. Differential Algebraic Techniques. World Scientific, 1999. Chapter in:
Handbook of Accelerator Physics and Engineering.

[12] M. Berz. Modern Map Methods in Particle Beam Physics. Academic Press, San
Diego, 1999. Also available at http: //bt .pa.msu.edu/pub.

[13] M. Berz. COSY INFINITY version 8.1 user’s guide and reference manual. De-
partment of Physics and Astronomy MSUHEP-20704, Michigan State University,
2002. http : / / cosy . pa . msu . edu/ cosymanu/ index . html.

[14] M. Berz and et al. The COSY INFINITY web page. http://cosy.pa.msu.edu.

136

 

[15] M. Berz, K. Makino, Y.-K. Kim, and P. Snopok. Long-term stability of large
particle accelerators. ECMI Newsletter, 39, 2006.

[16] M. Berz, K. Makino, and W. Wan. An introduction to the physics of beams.
IOP Publishing, forthcoming.

[17] G. I. Budker. Accelerators and colliding beams (in russian). In Proceedings of the
7th International Conference on High Energy Accelerators, volume 1, page 33,
Yerevan, 1969.

[18] D. C. Carey, K. L. Brown, and F. Rothacker. Third-order TRANSPORT with
MAD input. A computer program for designing charged particle beam transport
systems. Technical Report FERMILAB-Pub—98/ 310, Fermi National Acceler-
ator Laboratory, 1998. Also available at http://rampex.ihep.su/manuals/
transport.pdf.

[19] P. Chen. Disruption, beamstrahlung, and beamstrahlung pair creation. Technical
Report SLAC—PUB—4822, SLAC, 1988. Contributed to the DPF Summer Study
High Energy Physics in the 1990’s, Snowmass, Colorado, June 27-July 15, 1988.

[20] P. Chen and K. Yokoya. Disruption effects from the interaction of round e+e—
beams. Phys Rev. D 38, 38(3):987—1000, 1988.

[21] E. D. Courant and H. S. Snyder. Theory of the alternating-gradient synchrotron.
Annals of Physics, 3:1—48, 1958.

[22] B. Erdélyi. Private communication.

[23] B. Erdelyi. Symplectic Approximation of Hamiltonian Flows and Accurate Sim-
ulation of Fringe Field Eflects. PhD thesis, Michigan State University, East
Lansing, Michigan, USA, 2001.

[24] B. Erdélyi and M. Berz. Optimal symplectic approximation of Hamiltonian flows.
Physical Review Letters, 87,11:114302, 2001.

[25] B. Erdélyi and M. Berz. Local theory and applications of extended generating
functions. International Journal of Pure and Applied Mathematics, 11,3:241—282,
2004. Available at httpzl/bt .pa.msu. edu/pub.

[26] N. Gelfand. Reshimming the Tevatron dipoles. Technical Report Beams-doc-
1290—v1, FNAL, August 2004.

[27] C. Johnstone and A. Garren. An IR and chromatic correction design for a 2
TeV Muon Collider. In Proceedings of the Workshop on New Directions for
High-Energy Physics, page 222, 1996.

137

 

[28] C. Johnstone and A. Garren. 50-on-50 GeV Muon Collider storage ring. In Higgs
Factory 2001 Snowmass Report, 2001.

[29] C. Johnstone, W. Wan, and A. Garren. Lattice design for a 50—on-50 GeV Muon
Collider. In Proceedings of the 1999 Particle Accelerator Conference, pages 3065—-
3067, New York, 1999.

[30] C. J. Johnstone and N. V. Mokhov. Optimization of a Muon Collider interaction
region with respect to detector backgrounds and the heat load to the cryogenic
systems. Technical Report Fermilab—Conf—96/ 366, Fermilab, 1996.

[31] S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi. Optimization by simulated
annealing. Science, Number 4598, 13 May 1983, 220, 4598:671—680, 1983. http:
//citeseer.ist.psu.edu/kirkpatrick83optimization.html.

[32] V. Lebedev. OptiM — computer code for linear and non—linear optics cal-
culations, 2004. http : //www-bdnew . f nal . gov/pbar/organizat ionalchart/
lebedev/OptiM/optim.htm.

[33] K. Levenberg. A method for the solution of certain problems in least squares.
Quart. Appl. Math, 2:164—168, 1944.

[34] K. Makino and M. Berz. Remainder differential algebras and their applications.
In M. Berz, C. Bischof, G. Corliss, and A. Griewank, editors, Computational
Difierentiation: Techniques, Applications, and Tools, pages 63—-75. SIAM, 1996.

[35] K. Makino and M. Berz. Solenoid elements in COSY Infinity. Institute of Physics
CS, 175:219—228, 2004.

[36] K. Makino and M. Berz. COSY INFINITY version 9. Nuclear Instruments and
Methods, A558z346—350, 2005.

[37] D. Marquardt. An algorithm for least-squares estimation of nonlinear parame-
ters. SIAM J. Appl. Math, 11:431—441, 1963.

[38] M. A. Martens. Tevatron lattice page. http://www-ap.fnal.gov/~martens/
tev_lattice/tev_lattice.html.

[39] W. M. McKeeman. Algorithm 145: Adaptive numerical integration by Simpson’s
rule. Communications of the ACM, 5(12):604, 1962.

[40] R. E. Meller, A. W. Chao, J. M. Peterson, S. G. Peggs, and M. Furman. Deco-
herence of kicked beams. Technical Report SSC-N-360, SSCL, 1987.

[41] J. A. Nelder and R. Mead. A simplex method for function minimization. The
Computer Journal, 7:308—313, 1964.

138

 

[42] D. Neuffer. Colliding muon beams at 90 GeV. Technical Report FN-319, Fermi
National Accelerator Laboratory, 1979.

[43] R. B. Palmer, A. Sessler, A. Skrinsky, and et al. Muon colliders. In Proceedings
of 9th Advanced I CFA Beam Dynamics Workshop, 1995.

[44] M. L. Shashikant, M. Berz, and B. Erdélyi. COSY INFINITY’S EXPO symplectic
tracking for LHC. IOP CP, 175:299—306, 2004.

[45] P. Snopok. A converter program for Tevatron lattices from OptiM to COSY
INFINITY. Technical Report MSUHEP-40909, Department of Physics and As—
tronomy, Michigan State University, East Lansing, MI 48824, 2004.

[46] P. Snopok, M. Berz, and C. Johnstone. Calculation of nonlinear tune shift
using beam position measurement results. Technical Report MSUHEP-O70421,
Department of Physics and Astronomy, Michigan State University, East Lansing,
MI 48824, 2007.

[47] P. Snopok, C. Johnstone, and M. Berz. Simulation and optimization of the
Tevatron accelerator. In H. M. Biicker, G. Corliss, P. Hovland, U. Naumann,
and B. Norris, editors, Automatic Difi’erentiation: Applications, Theory, and
Tools, Lecture Notes in Computational Science and Engineering. Springer, 2005.

[48] P. Snopok, C. Johnstone, M. Berz, D. A. Ovsyannikov, and A. D. Ovsyannikov.
Study and optimal correction of a systematic skew quadrupole field in the Teva-
tron. In Nuclear Instruments and Methods in Physics Research A, volume 558,
pages 142—146, 2006.

[49] P. Snopok, C. Johnstone, and A. Garren. Recent progress in Muon Collider
lattice design. Technical Report MSUHEP-070420, Department of Physics and
Astronomy, Michigan State University, East Lansing, MI 48824, 2007.

[50] D. J. Summers. Muon acceleration to 750 GeV in the Tevatron, 2007. Neutrino
Factory and Muon Collider Collaboration, UCLA, available at http : //www. cap.
bnl.gov/mumu/conf/MC-O70129/DSummersl-O70129.pdf.

[51] M. Syphers. Private communication.

[52] M. Syphers. Strong transverse coupling in the Tevatron. Technical Report
Beams-doc-1159-v1, FNAL, May 2004.

[53] M. Syphers. Tevatron accelerator physics. Technical Report Beams-doc-1046-v1,
FN AL, February 2004.

[54] F. F. Tikhonin. On the effects with muon colliding beams. Technical Report
P2-4120, JINR, Dubna, 1968.

139

 

[55] D. 'Ifbojevic, K. Y. Ng, and W. Weishi. A lattice for the 50 GeV Muon Collider
ring. In S. Myers and et al., editors, EPA C ’98 Proceedings, pages 362—364, 1998.

[56] D. Tlrbojevic, RB. Palmer, E.D. Courant, J. Gallardo, S. Peggs, S. Tépikian,
and K. Y. Ng. A lattice for the Muon Collider demonstration ring in the RHIC
tunnel. In Proceedings of the 1997 Particle Accelerator Conference, 1997.

[57] S. Turner, editor. CERN Accelerator school, Fifth General Accelerator Physics
Course, volume I. CERN, Geneva, January 1994.

[58] W. Wan. Private communication.

140

  
   

IIIIIIIIIIIIIIIIIIIIIIIIIIIIII

ll[l][l[\[[[l[l][[lll[[\[ll

l
31 5 96