. .r:
5.243.... .1.
: ..\ .23.: .r .
. if; I: . 5..
$1. 2.?
.3. r ..

.3

 

 

 

a-s'e-o-o-Q-o

 

 

. 5.5.

4 i I 3., v...
2.2:,

..:)..
ail l: f,

 

0...

.ug. ,

34.. ., .T
s 4 53.x:

 

 

Lamas '

205\

This is to certify that the
dissertation entitled

HIGH ORDER FINITE ELEMENT METHODS TO COMPUTE
TAYLOR TRANSFER MAPS

presented by
SHASHIKANT MANIKONDA

has been accepted towards fulfillment
of the requirements for the

Ph.D. degree in Physics and Astronomy

 

 

H Q‘VKK EM

Major Professor’s Signature
8] 29/ oé:
l \

Date

 

MSU is an Affinnative Action/Equal Opponunity Institution

LIBRARY
Michigan State
University

 

 

 

 

II-.-I-Oun-9-0-0-n-o-o-I-c-<-o-o-9-I-a---.--.— - .I-I-O-0-0-U-O-O-O-D-l-I-l--O—O-t-I-I-C-'-'-O-U-O-I-I-O-I-l-¢-I-I-I-l-

 

PLACE IN RETURN BOX to remove this checkout from your record.
To AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

2/05 p:/ClRC/DateDue.indd-p.1

 

 

 

HIGH ORDER FINITE ELEMENT METHODS
TO COMPUTE TAYLOR TRANSFER MAPS

By

Shashikant Manikonda

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

2006

ABSTRACT

HIGH ORDER FINITE ELEMENT METHODS
TO COMPUTE TAYLOR TRANSFER MAPS

By

Shashikant Manikonda

In beam physics, map methods are important techniques for the design and analysis
of lattice structures. The computation of the transfer map for an electric or magnetic
element requires the multipole decomposition of the field in the region where the beam
passes. In the first part of this dissertation we present new techniques to extract the
multipole decomposition of the electric or magnetic field from the measured field data
or from the knowledge of the current distribution.

The new high precision technique developed to obtain the multipole decomposi-
tion of the field from the measured field data solves the Laplace equation using the
Helmholtz vector decomposition theorem and differential algebraic methods. This
technique requires the field to be specified on a closed surface enclosing the volume
of interest. The method outperforms the conventional finite difference and finite el-
ement methods in both the execution speed and the precision achieved. We extend
this technique to obtain a verified solution to the Laplace equation by using the Tay-
lor model methods. We also parallelize this technique and implement it on a high

performance cluster.

We then present a new high precision technique to find the magnetic field of an
arbitrary current distribution. The technique uses the Biot-Savart law and differential
algebra methods to compute the magnetic field. Using this technique we develop new
computational tools to design accelerator magnets.

Both these techniques can also be combined to solve the Poisson equation when
the source distribution is specified inside a volume and the field is specified on the
surface enclosing the volume. Besides providing a natural multipole decomposition
of the field both these tools have the unique advantage of always producing purely
Maxwellian fields.

We demonstrate the utility of these techniques in solving practical problems by
applying them to real life applications. We present the design and analysis of a
novel combined function multipole magnet with an elliptic cross section that can
simplify the correction of aberrations in the large acceptance fragment separators for
radioactive ion accelerators. We then apply the Laplace field solver to the measured
magnetic field data of the dipole magnet of the MAGNEX spectrometer and extract
the multipole decomposition of the magnetic field. Finally, we present the linear
and high order ion optic simulations for the proposed design of the superconducting
fragment separator (Super-F RS) and also apply the field solver technique to extract
the transfer map for the magnetic field data obtained through the TOSCA simulation

for the Super-F RS quadrupole magnet.

TABLE OF CONTENTS

LIST OF TABLES vii
LIST OF FIGURES ix
1 Introduction 1

2 Background information 3
2.1 Vector field and the Helmholtz theorem ................. 3
2.1.1 The Helmholtz theorem for Euclidean three-space ....... 4

2.1.2 The Helmholtz theorem for a finite volume .......... 7

2.1.3 The Helmholtz theorem for time dependent vector fields . . . 10

2.1.4 Maxwell’s equations from the Helmholtz theorem ....... 14

2.2 Laplace’s and Poisson’s equations .................... 16
2.2.1 Fundamental solution and Green’s functions .......... 17

2.2.2 Mean-value formulas ....................... 19

2.2.3 An analytical solution of the 2D Laplace equation ....... 20

2.2.4 Analytical and numerical solutions of the 3D Laplace equation 21

2.3 The Differential Algebra an ...................... 26
2.3.1 Solution of PDEs using DA ................... 28
2.3.2 Taylor transfer maps ....................... 28

2.4 Taylor Models .............................. 32
2.4.1 The Taylor Model integration scheme .............. 33
2.4.2 An example of Taylor Model expansion ............. 34

iv

3 3D Laplace and Poisson solver using DA 40

3.1 The Laplace solver ............................ 42
3.1.1 Methods using boundary data .................. 43
3.1.2 The two dimensional case .................... 44
3.1.3 The three dimensional case ................... 48

3.2 Verified solution of the 3D Laplace equation .............. 59
3.2.1 An analytical example: the bar magnet ............. 62

3.3 Parallel implementation of the Laplace solver ............. 65

3.4 Magnetic field due to arbitrary current distribution .......... 67
3.4.1 Field computations using the Biot-Savart law and DA . . . . 69
3.4.2 Line current example: circular loop ............... 70
3.4.3 Surface current example: spherical coil ............. 73
3.4.4 Three dimensional current distribution ............. 76

3.5 Extraction of transfer map from surface field data or current distribution 81

3.6 Summary ................................. 83

4 Applications 87

4.1 The conceptual design of a quadrupole magnet ............ 87
4.1.1 2D design of the quadrupole magnet .............. 90
4.1.2 3D design of the quadrupole and the fringe field analysis . . . 97

4.2 The dipole magnet of the Catania MAGNEX spectrometer ...... 109

4.3 Ion-optic simulations for the Super-FRS ................ 112
4.3.1 The pre-separator branch .................... 114
4.3.2 The main-separator ........................ 118
4.3.3 The energy buncher ........................ 118
4.3.4 The Super-FRS quadrupole magnet .............. 120

4.4 Summary ................................. 124

APPENDICES 128

A Comparison of the GICOSY and the COSY INFINITY beam physics

codes 129
B The code converter 134
C New COSY INFINITY tools for beam physics simulations 136

BIBLIOGRAPHY 139

vi

LIST OF TABLES

2.2.1 Analytic solution functions for different coordinate systems where the
method of seperation of variables can be applied ........... 22

2.4.1 The Taylor model expansion of the kernel. The coefficients up to second
order are shown ......................... 37

2.4.2 The Taylor model expansion of the kernel. The third order coefficients,
the reference points and the remainder bound interval are shown. . . . 38

2.4.3 The Taylor model expansion of the kernel. The third order coefficients,
the reference points and the remainder bound interval are shown. . . . 39

3.1.1 A sample eighth order Taylor expansion in two surface variables. . . . 56

3.1.2 A sample eighth order Taylor expansion in two surface variable and three
volume variables. ............ . . . . . ..... . 58

3.3.1 The code for the parallel algorithm. . . . ............. 68

3.4.1 A eighth order Taylor expansion of the analytic formula for the 2 com-
ponent of the magnetic field of a circular coil on the central axis. . . . 70

3.4.2 A sixth order Taylor expansion of the 2 component of the magnetic field
computed for the circular loop example ............... 72

3.4.3 A fourth order Taylor expansion of the 2 component of the magnetic
field computed for the spherical coil example. . . . . ........ 76

3.5.1 The difference between the analytic and the computed Taylor maps. The
difference map is shown to the third order .. . . . . . ....... 84

3.5.2 The difference between the analytic and the computed Taylor maps. The
difference map shows some of the fourth and the fifth order terms. . . 85

4.1.1 The center position of the current carrying coils in the first quadrant . 89

vii

4.1.2

4.1.3

4.2.1
4.3.1
4.3.2
4.3.3
4.3.4
4.3.5

4.3.6

4.3.7

A.O.l
A.O.2

B.0.1

The fifth order Taylor expansion of the magnetic field about the point
(0.0,0.0,0.0) for the current configuration producing a pure quadrupole
field .............................. 92

The fifth order Taylor expansion of the magnetic field about the point
(0.0,0.0,0.0) for the current configuration producing a pure hexapole

field .............................. 93
(a)Areas mapped in the dipole (b)Planes mapped in the dipole. . . . 110
The Super-FRS design parameters ................. 113
The relevant first order matrix elements for the pre—separator stage. . . 115
The Non-linear resolution at the focus PF2. ............ 117
The Non-linear resolution at the focus MF2 ............. 118

The extracted Taylor transfer map for the Super-FRS quadrupole. The
Taylor transfer map is shown to second order ............. 124

The extracted Taylor transfer map for the Super-FRS quadrupole. The
third order coefficients of the Taylor transfer map are shown ...... 125

The extracted Taylor transfer map for the Super-FRS quadrupole. Some
of the fourth and fifth order terms of the Taylor transfer map are shown. 126

System parameters used for the computation of the transfer maps. . . . 131

The first column shows the partical optical element. The column 2 shows
the RMS average of the coefficients of the difference map between GICO

and COSY INFINITY for the 5th order computation. And the third
column provides comments about the observed difference. ...... 133

The table shows the equivalent COSY INFINITY procedure calls for the
important GICO particle optical elements. . . ........... 135

viii

2.3.1

2.4.1

3.1.1

3.1.2

3.1.3

3.1.4

3.1.5

3.1.6

3.1.7

LIST OF FIGURES

Flow chart for extracting the Taylor transfer maps from either the elec-
tric or the magnetic field information or both. ...........

The figure shows a volume element centered at ($0,310, 20) element and
a surface element (x5, y3, 23). ..................

(a) The plot shows the dependency of n (r) on the radius r. (b) The plot
shows the dependency of 77 (7') on the radius as the radius r approaches
the boundary. 10,000 error sets (Ne) around the circle or radius R = 2
were chosen for the analysis. We show the plots for both the best and
the worst case scenario ......................

(a) The geometric layout of the bar magnet, consisting of two bars of
magnetized material. (b) The magnetic field By on the center plane of
the bar magnet. 130 = IT and the interior of this magnet is defined by
-0.5 g :1: _<_ 0.5, lyl s 0.5, and cs 3 z s 0.5 .............

The error for the field calculated for the bar magnet example with rectan-
gular grid for individual points (0, 0, 0), (0.1, 0.1, 0.1) and (0.2, 0.2, 0.2)
(0.3, 0.3, 0.3) ...........................

The error for the field calculated for the bar magnet example with cylin-
derical grid for individual points (0, 0, 0), (0.1, 0.1, 0.1) and (0.2, 0.2, 0.2)
(0.25, 0.25, 0.25) .........................

The error for the field calculated for the bar magnet example
with speherical grid for individual points (0,0,0), (0.1,0.1,0.1) and
(0.15, 0.15, 0.15) (0.2, 0.2, 0.2). ..................

The average error for the field calculated for the bar magnet example for
finite elements of width 0.4 around points (0,0,0) and (0.1, 0.1, 0.1). . .

The plot shows the dependency of the average error on the length of the
volume element. ........................

47

51

52

53

3.1.8

3.2.1

3.2.2

3.4.1

3.4.2
3.6.1

4.1.1
4.1.2
4.1.3

4.1.4

4.1.5
4.1.6
4.1.7

4.1.8

4.1.9

The plot shows the dependency of the average error on the number of
volume element. ........................ 60

The remainder interval width versus the surface element length for inte-
gration over a single surface element and vanishing volume size. . . . . 63

The remainder interval width versus the length of volume element for y

component of the magnetic field .................. 64
The cross section of a flux ball consisting of a sphere with winding on

its surface. .......................... 73
Magnetic field lines for the spherical coil ............... 74
The flow chart for extracting transfer maps from the measured magnetic

field data or the source distribution. ............... 83
The cross section view of the quadrupole. . . ..... . . . . . . 89
The layout of the racetrack coils in the first quadrant ......... 90

The operational plot for the quadrupole and the octupole. The coeffi-
cients are computed at the horizontal half aperture .......... 95

The operational plot for the hexapole and the decapole. The coefficients

are computed at the horizontal half aperture ............. 96
The schematic digram of a current coil. .............. 97
The three dimensional layout of quadrupole coils ........... 98

(a) The plot for B»; for the data on the plane 2 = 0 m. (b) The plot for
By for the data on the plane z = 0 m. . . . . ........... 100

(a) The plot for Bz for the data on the plane 2 = 0 m. (b) The vector
plot of the two dimensional field on the plane 2 = 0 m. ....... 101

(a) The plot for 8;; for the data on the plane 2 = 0.25 m. (b) The plot
for By for the data on the plane 2 = 0.25 m ............. 102

4.1.10 (a) The plot for B; for the data on the plane 2 = 0.25 m. (b) The vector

plot of the two dimensional field on the plane 2 = 0.25 m. ...... 103

4.1.11 (a) The plot for 8;; for the data on the plane 2: = 0.5 m. (b) The plot

for By for the data on the plane 2 = 0.5 m .............. 104

4.1.12 (a) The plot for B; for the data on the plane 2 = 0.5 m. (b) The vector

plot of the two dimensional field on the plane z = 0.5 m. ...... 105

4.1.13 (a) The plot for Ba; for the data on the plane 2 = 1 m. (b) The plot for

By for the data on the plane z = 1 m ................ 106
4.1.14 (a) The plot for Bz for the data on the plane 2 = 1 m. (b) The vector

plot of the two dimensional field on the plane .2 = 1 m. ....... 107
4.1.15 The plot of the magnetic field on the midplane, y = 0 m. Only the

magnetic field in the first quadrant is shown ............. 108

The dipole magnet of the MAGNEX spectrometer. ......... 108

4.2.1
4.2.2

4.2.3
4.3.1

4.3.2

4.3.3
4.3.4
4.3.5
4.3.6

4.3.7

C01

C02

The layout of the measurement grids in different regions of the dipole
magnet ............................. 109

The contour plot of magnetic field errors for the region 1 and plane A. . 111

The layout of the Super-FRS. Spatially separated rare-isotope beams are
delivered to the experimental areas via three different branches. . . . . 112

The figure shows the pre—separator layout and the projections of the

trajectories into the x and y planes ................. 116
Envelopes and the dispersion line for the pre—separator stage ...... 117
Ion-optical layout of the energy buncher. . . . ........... 119
The TOSCA model for the quadrupole magnet ............ 121

The difference between the relative error of the y component of the
magnetic field on the mid plane. The figure only shows the results for
the first quadrant. ....................... 122

The rms average difference between the TOSCA simulation result and
the new Laplace solver technique versus the volume element length for
four volume points (0.0,0.0,0.0), (-0.1,-0.025,-0.2), (-0.2,—0.05,-0.4) and
(-O.3,-0.075,—0.6). ....................... 123

The x—a phase space plots at the dispersive focal plane PF2 for a beam of
40 mm mrad and three different momenta of Ap/ p = $2.5 ‘70. The plot
(a) shows the phase space before second order aberration corrections,
and (b) shows the phase space after second order aberration corrections. 137

Ion-optic transmisson plot for the Super-F RS high energy branch. . . . 138

xi

CHAPTER 1

Introduction

The advent of computers has provided new means to solve the problems in physics.
Traditionally, the choice of the technique to numerically solve a PDE is driven by the
factors faster execution and the minimal use of the computational resources. Both
these factors are purely practical limitations due to the limited time and resources. In
the traditional numerical techniques the requirement of fast execution and minimal
use of resources can only be achieved at the cost of limiting the precision of the
numerical result. However, for many problems in physics that investigate phenomena
and processes at high energy (TeV) or in relatively small (nano/femto) length and
time scales, highly accurate results are an absolute necessity which the traditional
numerical methods can not provide. For instance, the Large hadron collider (LHC)
accelerator being built at CERN is designed for an energy of 14 TeV at the interaction
point and luminosity of 10350m_2 sec‘1 [5], and this requires magnets to be designed
with a magnetic field precision of 935: ~ 10—4 [37]. The modeling and study of
such instruments require high precision numerical tools. This leads to an additional
constraint for modern numerical techniques of obtaining very high precision results.
The traditional numerical techniques can not achieve all of the three conditions at

the same time, and hence fail to solve many modern problems. This leads us to

investigate new numerical techniques that are not only fast and make efficient use of
the computational resources, but also give high precision results.

One of the limitations with the conventional techniques comes from the fact that
the mathematical functions cannot be directly used on a computer. The treatment
of a function is done based on the treatment of numbers, and as a result, virtually
all the classical numerical algorithms are based on the mere evaluation of functions
at specific points. One way to overcome this difficulty is through the use of the local
Taylor expansion of a function about a point. We are then able to extract more
information about the function than just the value at a specific point. Once again
due to the limitation of the computational resources it is necessary to truncate the
Taylor expansion.

Algebraic operations like +, —, - and composition can be defined on truncated Tay-
lor expansions, leading to an algebra called the truncated power series algebra (TPSA)
[7]. The power of TPSA can be further enhanced by the introduction of derivations 8
and their inverses, corresponding to the differentiation and integration on the space
of functions. The resulting structure is called a Differential Algebra (DA) [67, 47].
The Differential Algebra provides a framework to develop techniques and algorithms
to use a truncated Taylor expansion of a function on a computer. The numerical
techniques based on DA have the unique advantage of getting high accuracy at a very
small cost of the execution time and the computational resources compared to the
traditional techniques.

In beam physics, the DA techniques were first introduced by M. Berz 1989 [12, 10,
9, 11] to compute the high—order Taylor expansions of the transfer maps. The beam
physics codes based on the DA techniques have become indispensable tools to the

design and analysis of accelerator lattices. In recent years the DA techniques have

been applied to solve DAEs, ODEs and PDEs [50, 42, 43, 28].

CHAPTER 2

Background information

In this chapter we present the background behind the work presented in this disser-
tation. We start by defining a vector field in the 3-dimensional Euclidean space E3
and discussing the relevant properties. We then discuss the Laplace equation and the
Poisson equation and present some of their properties. We also present a brief survey
of the analytic and numerical techniques to solve the Laplace and the Poisson equa-
tions. We then discuss the background of the numerical techniques that we utilize in

this dissertation, namely the Differential Algebra (DA) and the Taylor Model (TM).

2.1 Vector field and the Helmholtz theorem

Let V be a bounded region in the 3—dimensional (3D) Euclidean space IE3. A vector
field on the IE3 is a function E that assigns to each point (m,y,z) in V a three-
dimensional vector 5’ (17, y, z). A vector field that has zero curl everywhere is called
an z'r'rotatz'onal field. Such a field can be expressed as a gradient of a scalar field.
This scalar field is called a scalar potential. A vector field that has zero divergence
everywhere is called a solenoidal field. Such a field can be expressed as a curl of a

vector field. This vector field is called a vector potential.

2.1.1 The Helmholtz theorem for Euclidean three-space

The Helmholtz theorem [40, 59, 44, 62, 61] expresses any general vector field as a sum

of an irrotational field and a solenoidal field over all of a Euclidean three-space.

Theorem 1 (The Helmholtz theorem for Euclidean three-space) A general
continuous three-vector field defined everywhere in a Euclidean three-space, that along
with its first derivatives vanishes sufficiently rapidly at infinity, may be uniquely rep-

resented as a sum of an irrotational part and a solenoidal part.

The theorem imply that any vector field E (F) can be written as
[3' = —v¢,, + v x A}, (2.1.1)

where A} is a vector potential and (in is a scalar potential. A simple proof of this
statement follows directly from a well-know vector identity for an arbitrary vector
field,

—v217=t7x (in?) —V(V-l7)
Now by choosing B = —V2V, (in = V - l7 and At = V x V, we get the equation
(2.1.1). However, this assumes that there is always a solution to the vector Poisson

equation, B = —V2V. We now propose the following

Proposition 2 For any vector field E (F) that vanishes fast for large r, and satisfies

a Holder condition, a vector potential V given by

.. B’ (r’) ,

V (7") = -— [IE3 mdfl, (2.1.2)
is the solution to a vector Poisson equation

-o

[3' (r) = —v2v (7*). (2.1.3)

We now prove the proposition above. We first note that the volume integral in the
equation (2.1.2) exists only if the vector field B (F) vanishes fast for large r. Such a
vector field has a compact support in IE3, as a result the integral over all of Euclidean

three-space in the equation (2.1.2) can expressed as an integral on finite volume

-' J
_. B (r ) ,
V 47r[1"'- r l
where volume V is in IE3. We now take the Laplacian of V in equation (2.1.4) to get
~ .I
2* 2 B(r) ’
V V (F) = —V / —-—J—df2 , (2.1.5)
V 47r [F— r [

 

 

Br’
= _/V ERVZ 77—174] dn'.

 

Since the Laplacian is with respect to the unprimed coordinates, r, the integral and
the vector E can be brought out of the Laplacian. The final step is proving the

relation

~ J
_/VB_(Tlv2 _}_ dn’—_-B'(7=). (2.1.6)

_. J
47’ r—rl

 

Once we prove the equation above, we can use the above equation and the equation
(2.1.5) to complete the proof that the vector potential V (7") given by equation (2.1.1)
is the solution to the vector Poisson equation (2.1.3). To prove the equation (2.1.6)
we need all the three components of the vector field B,- (i = 1, 2,or 3) to satisfy a

Holder condition.

Definition 3 [61] Let r0 2 [F2 — F1] be the distance between two points F1 and F2.

If three positive constants a, n, c exist such that
[Bi (F2) — Bi (Fill < (1713’,

for all points F1 and F2 for which r0 < c, the quantity B,- is said to satisfy a Holder

condition.

The 8,; is also said to be Holder continuous. We are now ready to prove the relation
expressed in equation (2.1.6).

Proof. We note that

1 I

V2 _——J— =OVF§£F,
F—rl

 

thus the volume integral in equation (2.1.6) is zero except for the contribution from
the singular point 1" = r1. As 7‘! approaches F, the distance [1" — fl I tends towards
zero. We surround this singular point by a small sphere of radius 6, with surface 56
and volume V5. Since E (r!) satisfies a Holder condition, we can choose 6 so small
that for all values of 7‘! inside the sphere, B is essentially equal to its value B (7?!) at
the singular point. Then the integral in the equation (2.1.6) becomes

.. _1 _.
B (7‘) 1 , B (n 1 ,

f V2 ————— d0 = — V2 dfl . (2.1.7)

V 47r ;_ ill 47r V5

 

 

 

 

Using the relations,

 

 

1 I 1

V d J[ V a J] ,
T—T‘ T—T‘

2 I _ ’2 I

V —l”| ‘V “’l’
T—T T—T

 

and the divergence theorem on the volume integral we obtain

BOO—l) V2 1 d9, = Eff; (if!) -V, ——}-— dS, (2.1.8)
47" V6 F— fill 47f S 7"— 7"![
e

 

 

 

 

where n. (1"!) is the unit vector normal to the surface 3.; at f! . We note that

 

 

 

I 12
d5 = 17—77 dw,
. J
, 1 ”10")
V = ,
17—74] 7-:__.,:’[2

 

 

substituting this in the equation (2.1.8), we prove the relation expressed in the equa-

 

 

 

 

 

 

 

 

tion(2.1.6)
~ .1 . .1
B r “ -n r 2
1 I B I I
/()V2 dc: (flfif‘.____( F—f’dw
V 471' 7-:_,,-,’| 47f _. J2
5. _.~[
_gm 4
= dw=—B' .
4, f (n
Se
I

We see that the vector field B (7"!) must satisfy a Holder condition the equation
(2.1.7) to be valid. This requirement is stronger than continuity but less stringent
than differentiability. All infinitely often differentiable functions in Co0 satisfy this

condition.

2.1.2 The Helmholtz theorem for a finite volume

For a bounded region the statement can be modified as follows

Theorem 4 (The Helmholtz theorem for a finite volume) A general continu-
ous three-vector field that is defined everywhere in a finite volume V of a Euclidean
three-space and whose tangential and normal components on the bounding closed sur-
face S are given may be uniquely represented as a sum of an irrotational part and a

solenoidal part.
Once again, the theorem implies that any vector field E (if) can be written as

We now proceed to prove theorem 4, and obtain the explicit expressions for the vector
potential At and the scalar potential qbn in terms of the essential characteristics of the

vector field, namely, divergence, curl, discontinuities, and boundary values. We adapt

the derivation from [62]. By using the equation (2.1.6) a vector function E (2:, y, 2)

can be represented as

-* J
= *lIBiI)V2(|r_lIlld‘/’

~ .I
= —v2/V—1:(—T—Zl—dvl. (2.1.9)

 

 

Using the vector identity
Vx (Vxe') =V(V-B) 4723‘,

we rewrite equation (2.1.9) as

i ) /

 

 

 

 

 

 

~ 3(7") I 3'0") I
B(F)=Vx Vx/ -———,—dV —V V- —————7—dV . (2.1.10)
1 V47rr—rl J V47rlr—rl
\ .31? I \ <2 I

We first consider the divergence term (25,

~ J
g: fl.v( 1 )dv'. (2.1.11)

 

v47r

We note that

 

substituting this in the equation (2.1.11) and using the divergence theorem for the

second part of the volume integral,

gives the desired form of the scalar potential a). The curl term At in the equation

(2.1.10) can be written as

~ J
21’, = Vx/———B(T) dv’,
V r r

 

 

1 .. I I 1 I
= —/ B (F) x V I dV. (2.1.12)
4” V [r—r I

We note that

.. I I 1
VX = —B(F)XV ,

 

 

 

 

 

 

 

 

 

using this in equation (2.1.12) we get

_, v' [3' J B J
At: /V 47:}: _(;[)dvl_4in/1/VIX[:—%dv"

 

For the second integral on the right hand side of the above equation we now show the

following §(J .. g _l
1 VIX r)dV/=—1 was. (2.1.13)

2.. I lI-I’l lr-I’l

 

To prove this result we consider a constant vector C and apply the divergence theorem

to the quantity C - V, x (B (1"!) / [7" - 7"[),

 

 

 

fd-V'x éOVdV' = —/V'-Cx g TJ)dv'
I I—Il I lI-I’l ’
~ J
= —inxB(T,)-r‘id8,,
s [“7"]
= _é.ffl_(_:2d5’,
5 [7‘4]

since the constant vector 6’ is arbitrary we prove the relation (2.1.13). Using this

relation leads to the desired form of the vector potential At,

Iii =/VV, x §(#)dv'+ji-T-”—lai(r—Vds'.

I I
47r|r—r[ 47rlr—rl

 

Hence, the Helmholtz identity for any vector field B (F) in a finite volume is

 

 

I -. _1 _. .. J
3'0?) : _V [V -B(T’)dv,_fn‘B(T’)dS,
V47rr—r S47rr—r

 

 

 

 

 

+Vx /V,XB(7J)dv'+fEfiJ—)—ds'
V

I
47rr—r[ 47r[r—r|

 

2.1.3 The Helmholtz theorem for time dependent vector

fields

Recent works [30, 82, 81] have extended the Helmholtz theorem to time dependent
vector fields. We present a summary of some of those recent works here.

It is possible to obtain the Helmholtz theorem and an explicit form of the scalar
and vector time dependent potentials by replacing the three dimensional analysis on
the Euclidian three space in the sections 2.1.1 and 2.1.2 with a four-vector equivalent
on the Euclidian four space IE4, or the Minkowski four-space R3+1 [81, 80]. However,
the approach we describe below is more suitable for our purposes, as it allows us to
not only obtain the Helmholtz theorem for the time dependent vector field but also to
obtain Maxwell’s equation starting from the Helmholtz theorem and the continuity
equation.

Let f (F, t) be a time dependent scalar function and let g (F, t) be a time dependent
scalar function defined by

_j _‘ J
g(F,t) z/fV ’t_ “T VG) dv', (2.1.14)

_. J
47rr—rl

 

 

 

10

where c is a constant. We introduce the notation, If 2 F — F], R = [If], i/J (R) =
1/ (47rR), and r = t — R/c. Using this notation, the equation (2.1.14) can be written

as
I I
g(F, t) = ff (FJ) w(R)dV. (2.1.15)
The gradient of g and its divergence, i.e., the Laplacian of g can be given as
I
W (F. I) = / (WI + W) W .

V29 (F. t)

j (W2 f + 2v f - Vi/J + fvzw) dv'

The gradient of the scalar function f and w are given by

_ 16f.
Vf — -Z-a-;6R, (2.116)
1 .
Vii) = "EV“?R,

and taking the divergence of V f , we obtain

182f 28f

2 ___
V f _ ?E’ — fig; (2.1.17)

Using the equations (2.1.16) and (2.1.17), the Laplacian of g can be given as
J
2 ~ 16247”) ’ J 2 ’
V g(r,t) = ——————-—z,/I(R)dV + f r ,r V de. (2.1.18)
c2 Brz
We apply the equation (2.2.4) derived in the section 2.1.2 to the three dimensional
vector field with one of its components being f (F, r) and other two being zero to

obtain

/ f (14,7) v2¢dv' = — f (F, t). (2.1.19)
Now, by changing the differentiation with respect to r in the first term on the right
hand side of the equation (2.1.18) to differentiation with respect to t, and bringing
the differentiation with respect to t outside the integral, and using equations (2.1.19)
and (2.1.15) we obtain
_1_529 (7‘? t)
02 6t2

11

V2907»): —I(r.t>-

We can now rewrite the equation above as

2
(V2 — $567?) g(F, t) = —f(F, t). (2.1.20)

We now consider a time dependent vector function E (F, t) and define a vector

potential A (F, t) by

.. .. I I
A(F, t) = [F (F,r) 7,12(R)dV, (2.1.21)
and using the equation (2.1.20) component—wise for the vector A, we get
2
2 I a .. ., _. _,
V ——— A r,t =—F r,t . 2.1.22
( C, 6,2) < > < > < )

We now define a scalar potential (15 associated with 13’ (F, t) such that it is a solution

to
16¢ (I: t)
c2 at

 

= —\7 - 211m). (2.1.23)

The equation above is the well-known Lorenz gauge condition. We now use the
standard vector identity for the double curl of a vector field on the vector A and use

the equations (2.1.23) and (2.1.22) to get

Vx (v x 2177110) = V(V-A(F,t)) —V2A(F,t)

 

 

Rearranging the equation above we arrive at the Helmholtz theorem for time dependent
vector fields,

Ezcl—a-(V¢+%%>+VX(VXA).

Now we derive additional relations which would facilitate the derivation of

Maxwell’s equations from the Helmholtz theorem in the section 2.1.4. Starting from

12

the equation (2.1.21), we evaluate the divergence of A to obtain

V-A(F,t) = [V [¢(R)F(r~’,r)] dv'

/ [—V7’zp-13"+IIIV-13[dVI

= f[_v'.(vfi)+uv'-F+wv-F]dv'
= j—v'-(vfi)dv'+/[v'-F+v-F’[vdv'.

\_v——-’
I

Using the divergence theorem, we can write the integral I as

I=/—v'.(vfi)dv'=—]i¢fi.d§'.

Assuming that the vector function E vanishes sufficiently fast as r —> 00, we see that

the integral I vanishes. Thus, we have

V-A(F,t) =/[v’.fii(r~",7) +v.fi(r’,v)] deI.

We rewrite V, - 1:" and V . E as

 

 

3
I _. I
v -F(r,7) = Z , + —7
i=1 _ 613i 6T 62:2.
- ~ J
3 BF-(r T)
_. _‘l t I 67'
V'F(7’,T) = 2 Tag].
i=1[_ 1

We note that (Br/023;) = — (87/827,). Using this, we can conclude that

'3 3131'. J,
V.,If(f"t) = / 24:2 del
i=1 555.-

 

? I .. I
/ V7. - F] de , (2.1.24)

I I
where the notation VT is the gradient operator with respect to $2. with r kept constant,

for i = 1, 2,3. The scalar potential in the Lorenz gauge equation (2.1.23) can be

13

expressed in terms on another scalar function pa (r, t) as

I I
(Ia: t) = f p. (m) II (adv. (2.1.25)
and also from the equation (2.1.20), the scalar potential 95 (F, t) also satisfies
2
2 1 3 .. _ s
(V — 32-5?) ¢(T,t) — —,0a (r,t). (2.126)

Using the equations (2.1.24) and (2.1.23) and assuming that the time derivative can

be exchanged with the spatial integral, we get

810a (7?, t)

__ 2 ._. ~
8t _ cV F(r,t). (2.1.27)

2.1.4 Maxwell’s equations from the Helmholtz theorem

We will now obtain Maxwell’s equations using the Helmholtz theorem for time de-
pendent vector fields and the continuity equation

39 (F. t)

V-J(r,t)+ 8t

 

= 0, (2.1.28)

where I (F, t) is the current density vector and p (F, t) is the charge density at point F
and at time t.

We start from the Helmholtz theorem for time dependent potentials

- 16 62? -
F = as (V¢+ a?) +v x (V x A). (2.1.29)

We choose E (F, t) to be uof (F, t), where #0 is permeability of the free space. Using
the equation (2.1.21), we see that the vector potential A can be written as

K: [Mm/I

4rrR ’

and using the equation (2.1.27), we get

floflg’i) : —C2V'#0j(fit)

= —IIoc2v.J‘(r.t) =uoc2%,’fl. (2.1.30)

14

In the equation (2.1.30) we use the continuity equation (2.1.28). We now attribute
a meaning to the constant c as the speed of light, and c = 1/,/,uoe , and 60 is the

permittivity of free space. We thus get

 

(2.1.31)

Using the equation (2.1.25), we define d), the electric scalar potential, as
(J l
p r ,r ,
",t = ——dV .
¢ (T ) / 41reOR
Additionally, we now define two new quantities. The electric field intensity E (F, t) is

defined as _.
E (F, t) = — (Vd (F, t) + 8_A5(:,_t)) , (2.1.32)

and the magnetic induction B (F, t) is defined by
B(F,t) = V x A(F,t).
From the definition of B (F, t) it follows that the divergence of B (F, t) vanishes, that
is,
v - B (F, t) = 0. (2.1.33)
Now the equation (2.1.29) becomes the Ampere-Maxwell equation

vXB(I=:t)=IIoJ(I=:I>+i26 0" ).
c at

 

(2.1.34)

Taking the curl of both sides of the equation (2.1.32) and using the definition of

B (F, t), We get
_aB(I'~:t)

VXE(F,t)= at

 

(2.1.35)

Taking the divergence of both side of the equation (2.1.32) and using the equations

15

(2.1.23), (2.1.26) and (2.1.31), we get

 

 

V~E<I~zt> = —v2¢<r. #5559759]
_ 2 _2 _iacfl'fit)
_ —v¢(r, ) mic? at]
_ _ 2__1___3: 7... _p(F.t)
_ (v 62&2)¢(,t)_ .0 . (21.35)

We have derived Maxwell’s equations (2.1.36), (2.1.35), (2.1.33), and (2.1.34) start—

ing from the Helmholtz theorem (2.1.29) and the continuity equation (2.1.30).

2.2 Laplace’s and Poisson’s equations

A vector field that is both solenoidal and irrotational and is continuously differentiable
is called a Laplacian vector field. Since this field is irrotational, there is a potential (25
such that

B: _ V¢7

and since the field is solenoidal and continuously differentiable, this potential satisfies

Laplace’s equation,
9% 6% 6%

+— +— =0, 2.2. 1
6:1: 2 6y2 fizz ( )

the corresponding inhomogenous equation is called Poisson’s equation,

v2¢— _

—V% = f. (2.2.2)

Laplace’s and Poisson’s equations describe a wide variety of physical processes,
e.g. the gravitational potential, the electrostatic potential, the magnetostatic poten-
tial, chemical concentrations, steady-state heat conduction, steady-state diffusion of a
solute, steady-state irrotational flow of the incompressible fluids, and the probabilistic
investigation of Brownian motion. In a typical interpretation, (b denotes the density
of some quantity in equilibrium, as in the case of an electric potential. Laplace’s

equation is equivalent to the requirement that there be no maximum or minimum of

16

qt inside the volume V and the value of (15 at a point be equal to the average of the
values of d at the neighboring points.

A 02 function that satisfies Laplace’s equation is called a harmonic function. In
this section we discuss properties of harmonic functions which are relevant to the work
that we present. All theorems presented are valid over R”, but since the problems
that we address are all in two or three dimensions we restrict our attention to the
special case n = 2 or n = 3. We also discuss the established techniques to solve

Laplace’s and Poisson’s equations. The material presented here is adapted from [36].

2.2.1 Fundamental solution and Green’s functions

In the section 2.1.1 we proved that for any vector field B (F) the vector potential V (F)

 

given by
_, g (I!) ,
V(F) = —/ ——-7dfl , (2.2.3)
1E3 47r F — r [
is a solution to the vector Poisson equation B = -—V2V. In the equation above,

irrespective of what the vector field B (7") is, the kernel of the volume integral,
1/ (4n [F — F I), will always remain the same. This observation leads us to define

the fundamental solution (I) of the Laplace equation as

<1)(I‘”)=1/(47r|7"l)-

We can also arrive at a fundamental solution by finding an explicit solution to the
Laplace equation with certain symmetry properties [36]. For a two dimensional case it
can be easily shown that a fundamental solution is given by CD (F) = — (log lFl) / (2n).
The fundamental solution (I) has a property that it is invariant under rotation and it
is singular at origin. Using the notation for the fundamental solution the equation

(2.2.3) can be written as a convolution integral

17m = —/E3<1> (75— I") ['3’ (Fl) 39'.

17

We see that the recognizing the fundamental solution of the Laplace equation al-
lows us to construct the solution to the more complicated Poisson equation using a
convolution integral.

The Green’s function is a fundamental solutions that also satisfy boundary condi—
tions or initial conditions. To discuss Green’s function it will be useful to first present

Green’s formula.

Theorem 1 (Green’s Formulas) Let d, i/I E C'2 on a closed domain U C IE3. Then

2. fU VI; . wdv = — fav ¢V2de + faU aggds
3. III [3v2¢ — M729] (iv = fan [33% — 53%] d5

All three Green’s formulas above can be easily proved using the divergence theorem
[36]. Greens’s formulas along with the definition of Green’s function can be used to
prove the uniqueness of the solution for certain boundary value problems. We now

formally define Green’s function.

Definition 2 (Green’s function) Green’s function for the region U is of the form

G (F, I") = <2 (77— r") — \II (5', Fl) , (2.2.4)

I
where (I) (F — F ) is a fundamental solution and \II is a harmonic function.

The boundary conditions will depend on the problem. For the special case of the
Dirichlet and the Neumann boundary condition, it can be shown using the definition
of Green’s function and Green’s formulas that the unique solution can be found for

the Poisson problem —V2¢ = p in U. The solution to this problem can be expressed

 

For the Dirichlet boundary conditions where the potential d is specified on the bound—

ary 6U, Green’s function can be chosen to be
I
G' (F,F >— — 0 on the boundary 8U,

and for the Neumann boundary conditions where the normal component of the gra-
dient of the potential (15 is specified on the boundary 8U, Green’s function can be

chosen to be

III , >
—-5,—- = 0 on the boundary BU.
V

Green’s function is uniquely determined by the equation (2.2.4) and the boundary
conditions. Green’s function is symmetric with respect to F and 7"], G (F, Fl) =
G (7"! , F) . Except for few simple geometries, like the half-plane and the sphere, it

is usually very hard to find Green’s functions. For many problems Green’s function

may not exist or may not be uniquely determined.

2.2.2 Mean-value formulas

Consider an open set U C E3 and suppose (25 is a harmonic function within U. Let
8 (F0, R) denote a ball at F0 with radius R in U. The mean-value formulas declare
that 45 (F0) equals both the average of <1) over the ball (95 (F0, R) and the average of ()5

over the entire ball S (F0, R).

Theorem 3 (Mean-value formulas for 3D Laplace ’3 equation). Ifqb E C 2 (U) is har-

monic, then

¢ (F0): 4—7r—LRZ f dds =4——7r:l‘f3 [$713110 (1X19 (2.2.5)
65(F0, R)

for each sphere of radius R with center at point F0, 5' (F0, R) C U.

The converse of this theorem is also true. One of the consequences of the mean

value theorem is that the maximum or a minimum of a Harmonic function can occur

19

only on the boundary. If the function has a maximum or minimum inside U then the
function is just a constant. This result provides another way to establish uniqueness

of solutions in certain boundary value problems for Poisson’s equation.

Regularity

The regularity theorem says that if (I) E C2 is harmonic, then necessarily (15 E C 00.
Thus harmonic functions are automatically infinitely often differentiable. Even
though the Laplace equations itself has only second order partial derivatives it has

the interesting feature that the all partial derivatives of the solution d) exist.

Analyticity

An analytic function 96 is an infinitely differentiable function such that the Taylor

series around any point 2:0

 

To) = Z " if”) (as—x0)”.
n=0 '

is convergent for :1: close enough to 2:0, and its value equals (f2 (2:). The analyticity

theorem states
Theorem 4 (Analyticity) Assume (15 is harmonic in U. Then (i is analytic in U.

This property allows us to express a harmonic function qt as a Taylor series.

The last important property is the Harnack’s inequality, which assert that the
values of a non-negative harmonic function within V are all comparable. The value
of a harmonic function (I) cannot be very small or very large at any point of V unless

(Z) is very small or very large everywhere in V.

2.2.3 An analytical solution of the 2D Laplace equation

Let 2 = :z:+iy, dz = d$+idy and d2 = dx—idy. Then, (6/62) = % ((6/61) + i (6/6y))
and (6/62) 2 % ((6/ 6:17) — i (6/6y)). The Laplacian operator in 2D can be given as

20

A = 4 (6/62) (6/6z). For any C1 complex function f on a complex domain S2 with

boundary 652, the Cauchy-Pompeiu integral representation can be written as

f (20) — 1 Mdz — if Q dedy. (2.2.6)

— 27; 39 z - 2() 1r z — 20
If the function f is analytic then it satisfies the Cauchy-Riemann equation (6 f / 62) =

0, and the equation (2.2.6) reduces to the standard Cauchy’s formula

f(20) —' —1- HZ)

-— , dz.
2m 39 z — 20

 

Cauchy’s formula tells us that the value f at any point 20 in Q is completely deter-
mined by the value of f on the boundary 6D. This can also be used to develop a

numerical scheme to solve the 2D Laplace equation.

2.2.4 Analytical and numerical solutions of the 3D Laplace
equation

A solution of Laplace’s equation is uniquely determined if the value of the function
is specified everywhere on the boundary (Dirichlet boundary conditions) or the nor-
mal derivative of the function is specified everywhere on the boundary (Neumann
boundary conditions) or a linear combination of the solution and its normal deriva-
tive is specified on the boundary (Robin boundary conditions). For many problems
neither of the three boundary conditions above is suitable. We may then use the nat-
ural boundary conditions. Natural boundary conditions usually set the solution to a
distinct value at infinity (asymptotic boundary conditions). The kind of boundary
condition can vary from point to point on the boundary, but at any given point only
one boundary condition can be specified. When the region on which the PDE prob-
lem is posed is unbounded, one or more of the above boundary conditions is usually

replaced by a growth condition that limits the behavior of the solution at infinity.

21

 

ytlc so ution ns systems w the
method of seperation of variables can be applied.

The method of separation of variables

The method of separation of variables is a suitable technique for determining solutions
to linear PDEs, usually with constant coefficients, when the domain is bounded in
at least one of the independent variables. It turns out that in the three—dimensional
Laplace’s equation, there are some coordinate systems in which the solution takes
on the form R(§1,§2,€3) ~X1 (£1) - X2 (52) - X3 (£3), where the additional factor
R is independent of the separation constants. Laplace’s equation can be solved by
separation of variables in 11 coordinate systems [76, 59]. The form these solutions take
is summarized in Table 2.2.1. In addition to these 11 coordinate systems, separation
can be achieved in two additional coordinate systems by introducing a multiplicative
factor.

In addition to the method of separation of variables, the finite Fourier transform
and the power series method can be used to find an analytic solution to the Laplace

and the Poisson boundary value problems.

Numerical methods

For all the electromagnetic problems that cannot be solved analytically, numerical

methods are the only way to proceed. Most of the industrial packages are based on

22

one of the three classes of numerical techniques [6]:

The finite difference method (FDM) discretizes the differential operator at
each point of a rectangular grid covering the entire region of interest. The differential
operator is approximated by an algebraic expression (difference formula) with refer-
ence to the adjacent nodes. This leads to a large system of linear equations and the
solution requires inversion of large and sparse matrices. The FDM usually utilizes
uniformly spaced grids. The results usually are less accurate than other methods.

The finite element method (FEM) is based on division of the volume of space
in which the Laplace or Poisson equation is satisfied into small volumes (the finite
elements). Within each finite element a simple polynomial is used to approximate
the solution. To obtain the polynomial a variational formulation over the volume is
used. The variational quantity to be minimized is the total electrostatic or magne-
tostatic energy stored in a region. Element geometries and unknowns are expressed
by polynomials with nodal values as coefiicients. Relating these approximations to
the operator equation through minimizing the energy functional yields the solution
at the nodes. The FEM utilizes either uniform or nonuniform grids and it is possible
to implement automatic mesh size control. Even though, this method is considered

superior to FDM, but still suffers from some drawbacks. Among them are:

a The finite element techniques requires the mesh be extended to a reasonable
distance with either potential or derivative boundary conditions applied to the
outer surface, so that the truncation has an insignificant effect on the region of

interest. This usually leads to a large number of finite elements.

0 For problems where the magnetic or the electric field in the region of inter-
est differ by a large ratio (”106 ) in the maximum and minimum value, the
FEM elements have to start with extremely small element size and gradually
adapt themselves to much larger element size. This will increase the number of

unknowns substantially.

23

Most of these techniques using F DM or FEM utilize relatively low approximation
order and provide solution as a data set in the region of interest. They also require
a prohibitively large number of mesh points and careful meshing. Both FDM and
FEM are geared to solve electric and magnetic potentials. Since in beam physics
applications the knowledge of the values of electric or magnetic field is required,
these values have to be extracted by numerical differentiation of the potential. But
the numerical differentiation process is very sensitive to numerical, truncation and
round off errors.

The boundary element method (BEM) or source based field models the field
inside of a source free volume due to a real sources outside of it can be exactly
replicated by a distribution of fictitious sources on its surface. The error due to
discretization of the source falls off rapidly as the field point moves away from the
source. Since the discretization is done only on the surface, the dimensionality of the
problem is decreased by one. It makes the modeling of the problem easier and more
user friendly. The trade off here is that the matrices generated by BEM are usually
smaller and denser matrices. One technique that falls in this category is the image
charge method. This requires proper choice of planes / grids to place point charges (or
Gaussian distribution) and solve a large least square fit problem to find the charges.
This involves a lot of guess work and computation time involved in getting the right
solution. Knowing these charges, the potential and field can be directly computed
everywhere in space. In problems where extreme ratios exist between smallest and
largest details of the structure, BEM will be the method of choice. Finally, since the
field and the potential everywhere in space are being computed from the actual charge
distribution on the boundaries, it will result in extremely accurate results when this
method is applied to particle ray tracing.

BEM methods are only applicable to problems for which Green’s functions can be

defined. This places considerable restrictions on the range and generality of problems

24

to which boundary elements can usefully be applied. And once again these methods
usually compute the potential and not the field.

New methods that use the Helmholtz vector decomposition theorem are being used
in recent years to overcome the difficulties encountered by the BEM. The techniques
based on this method have the added benefits of giving the field directly and are
particularly suitable for beam physics applications. The references [63, 74, 73, 72, 75]
describe techniques based on this method to compute multipole decomposition of the
fields. Other applications based on this method include vector tomography [60] and
incompressible flows [33].

In beam physics, the detailed simulation of particle trajectories through magnets
in spectrographs and other large acceptance devices requires the use of detailed field
information obtained from measurements. Likewise, for high energy accelerators like
the LHC, higher order description of the beam dynamic via one-turn maps is required
to study the long term beam stability [69, 35]. The construction of such high order
one-turn truncated Taylor maps [17] requires the precise information of the electro-
magnetic field in the individual electromagnetic components (quadrupoles, dipoles,
sextupoles etc.) of the lattice.

It is commonly known that for a device that satisfies midplane symmetry, the entire
field information can be extracted from the data in the midplane of the device [17].
However, it is well known that this method has limitations in accurately predicting
nonlinear field information outside the immediate vicinity of the midplane because
the extrapolation requires the computation of higher order derivatives of in-midplane
data, which is difi'icult to do with accuracy if the data is based on measurements.
Thus it is particularly useful to employ techniques that rely on field measurements
outside the midplane. In particular, in modern particle spectrographs it is common to
measure the fields on a fine mesh on 2 to 4 planes outside the midplane. These data

have frequently been used to model the overall field as a superposition of point-charge

25

fields of so—called image charges [32, 18]. However, the computational effort required
for this approach is large, as it requires the inversion of a matrix with a dimension
equal to that of the number of image charges.

However, the out-of-plane field measurements in essence provide field data on the
top and bottom surfaces of a box containing the region of interest through which the
beam passes. If the planes extend outward far enough to a region where the fringe
field becomes very small, or can easily be modeled, and inwards far enough that the
field becomes rather homogenous, field data are known on an entire surface enclosing
the region of interest. The method we present can extract the field information as a
multipole expansion in the volume of interest if a discrete set of field measurements

are provided on a closed surface enclosing the volume of interest.

2.3 The Differential Algebra an

For real analytic function f in v variables, we form a vector that contains all Taylor
expansion coefficients at 53' = 0 up to a certain order n. The vector with all the
Taylor coefficients is called the DA vector. Knowing this vector for two real analytic
functions f and 9 allows to compute the respective vector for f + g and f - 9, since the
derivatives of the sum and product function is uniquely defined from those of f and
g. The resulting operations of addition and multiplication lead to an algebra, the
so-called Truncated Power Series Algebra (TPSA) [7]. The power of TPSA can be
enhanced by the introduction of the derivations 6 and their inverses, corresponding
to the differentiation and integration in the space of functions. This resulted in
the recognition of the underlying differential algebraic structure and its exploitation,

based on the commuting diagrams for the addition, multiplication, and differentiation

26

and their inverses:

f.g —T—» F,G f,g —T—» no
+.—[ [ea-J] [6W (2.3.1)
figTFgG ffg—TFgG
f L F
Will [50,551

6f,6‘1f —T—> 3015,0611:

In the equation above the operation T is extraction of the Taylor coefficients of
prespecified order n of the function. Thus, the operation T can be used to create
a DA vector from a function. The operation T is an equivalence relation, and the
application of T corresponds to the transition from the function to the equivalence
class comprising all those functions with identical Taylor expansion to order n. The
symbols ED, 9, G, O, 60 and 661 denote operations on the space of DA objects
which are defined such that the commuting relation expressed in the equation (2.3.1)
holds.

Using the Differential Algebra and analytic formulae to compute the n—th order
derivative of an univariant fundamental function, like sin z, cos 2:, log 2:, tans: etc., we
can compute the derivatives up to order n of functions in v variables. The detailed
description of obtaining the n-th order Taylor expansion for a multivariate function
that is expressible on a computer is described in [17]. Also, composition of two
multivariate functions can be defined using the DA. Many problems involving the use
of DA techniques can be formulated as fixed-point problems, such as the inversion of a
multivariate function. Here fixed point theorems can be applied to show that existence
of the solution, and this also provides a practical means to obtain the solution [39].

As mentioned before the DA techniques have been widely applied in beam physics,
asteroid problems and other problems involving the solution of ODEs or PDEs. The

focus of the work we present is finding a solution to the Laplace and Poisson equation

27

using DA techniques. In this context it is worthwhile to look at the techniques that

already use DA to solve PDEs [17].

2.3.1 Solution of PDEs using DA

The complicated PDEs for the fields and potentials stemming from the representation
of Maxwell’s equations in particle optical coordinates could be iteratively solved in
the DA framework to any order in finitely many steps by rephrasing them in terms

of a fixed-point problem [17, 54]. For example, consider the general PDE

6 6d 6 6d 6 6(2) _
016:1:(@620)+b16y<b26y)+cl6z(cz6z) — 0’

where a1, b1, cl, a2, b2, and 02 are functions of :13, y and z. The equation above can be

written in fixed point form

_ i 92 __ 912 .32 :13 :92
¢(:r,y,z) — ¢Iy:0+/yb2 {6yly:0 [g(bl 6:1: (£12627) + b16z (626z>)}'

If d and 6¢/ 6y are specified on the y = 0 plane then it is possible to iteratively
calculate various high orders in y. Techniques based on the fixed point scheme for
PDEs can be used to solve the Laplace equation when the potential d) and the normal
derivative 66 / 6y are specified e.g. on the midplane of a magnet. Also, by considering
the Laplace equation in cylindrical coordinates, it is possible to devise a technique
to obtain the magnetic field of a magnet with cylindrical symmetry when the field is

specified on the central axis.

2.3.2 Taylor transfer maps

An ensemble of particles with similar coordinates is called a beam. Detailed under-
standing of the beam dynamics requires the study of the motion of the reference
particle as well as the motion of the particle in the relative coordinates. The position

and momenta are usually sufficient to describe the motion. The state vector is given

28

by
2(8) : (wipivryapyrzapZ)7

where the point (:r, y, z) and the momentum vector (pm, py, p2) gives the position and
the momentum of a particle, and the arclength s along the reference orbit is used as
the independent variable. The space spanned by the state vector is called phase space
and the volume of the phase space is called emittance. At each point on the reference
orbit it is possible to define an unique orthogonal coordinate system, denoted by
(em, éy, es), satisfying a certain set of conditions [50, 17]. In this coordinate system

the motion of the particles in the beam can be described using relative coordinates,

I x )

which are given by

a =Px/P0
" y
Z s) =
( b=Py/P0
l= [C(t—to)

 

 

\5= (E-Eol/EOI
where the position (2:, y) describe the position of the particle in the local coordinate
system, p0 is a fixed momentum and E0 and to are the energy and the time of flight
of the reference particle, a and b are the momentum slopes, E is the total energy,
and It has a dimension of velocity which makes l a length like coordinate. The point
Z = 0 corresponds to the reference particle.

The transfer map or transfer function M relates initial conditions at so to the

conditions at 3 via
2 (s) = M (so. 9 (2180)) .
Transfer maps are origin preserving, M (0) = 0, and the transfer maps have the

property that
M (31,32) 0M (30,31) = M (50,32),

which says that the transfers maps of systems can be built up from the transfer maps

of the pieces. For a deterministic system (a unique solution exists) the transfer map

29

is the flow of ODEs
_ = f(f, s), (2.32)

in the independent variable 3. For most dynamical systems, the transfer map can not
be represented in a closed form. For weakly non-linear systems, like an accelerator
system, the map can be expanded as a Taylor series. Usually, the Taylor expansion
converges rapidly. Implementation of such a map on a computer would require the
map to be truncated at a certain order. From the implementation point of view the
Taylor transfer map is an array of DA vectors with each DA vector being an array
storing the coefficients of the truncated Taylor series expansion of the final phase
space coordinate in terms of the initial phase space coordinates. The Taylor transfer
maps can be used to replace and speed up element-by—element tracking, to look at
aberration content, or to monitor a design process [45]. A detailed discussion of the

properties and use of the Taylor transfer maps can be found in [17].

Generating transfer maps

To generate a transfer map we have to solve the equation 2.3.2. For beam physics
applications, the set of ODEs describing the equations of motion when the reference

trajectory is restricted to a plane and the energy is conserved are given by

:c =a(1+hx)—

W=b1+h)
::{( (1+h) p—1—+np—9— }i,
)1+nops v0

B B
—+b—ip-Q——i[ -(1+h:z:)+hm
1 “l" 770 P: X30 XmO Ps XmO Ps

::[__1_+77P_0E_y_+ Ba: _ (132 P0
=0

 

—] (1+ hsr),
1+ no P3 XeO XmO aXmO ps

, (2.3.3)

30

 

Beam Physics

 

 

 

 

 

 

 

 

ODE's
. . ODE Integrators Magnetic
Dx'fereinrgal (Runge-Kutta) AND/OR Electric
3 DA Integrators Field Information
Taylor Map

 

Figure 2.3.1. Flow chart for extracting the Taylor transfer maps from either the
electric or the magnetic field information or both.

where h is the radius of curvature, the ratio pO/ps is given by

1

IE: 77(2+77) .TE—aQ—IP -2
Ps 710(2+770) m3 ’

the ratio 77 of the kinetic to rest mass energy is given by 77 =
(E — eV (:r,y,s))/ (mcz), the XmO = 100/ (ze) is the magnetic rigidity and Xe0 =
(vaO) / (ze) is the electric rigidity. In the equation (2.3.3) Bx, By, B; and EC, Ey, Ez
represent the 2:, y and 2 components of the magnetic field and the electric field.
Using the differential algebra and the knowledge of the magnetic and the electric
field we can solve the equation (2.3.3) by using one of the ODE integration schemes.
Figure 2.3.1 shows the flow chart for the extraction of transfer maps. Traditionally,
it was only possible to extract the transfer maps to low orders (3 3), using the
perturbative analytic approach. However, the introduction of differential algebra to
beam physics has made it possible to extract maps to in principle arbitrary order.

Further, the direct availability of the derivation 6,- and their inverses 6i— 1 allows

31

to devise efficient numerical integrators of any order. The code COSY INFINITY
makes it possible to use DA integrators or a Runge—Kutta integrator of order eight
with automatic step size control based on a seventh-order scheme for this purpose.
The Runge—Kutta integration can be carried out not only with real numbers but also
for DA vectors. At each time step the eighth order Runge—Kutta scheme requires the
evaluation of the function f at thirteen points to compute the Taylor transfer map
[51, 64, 48], which in turn requires the electric and magnetic field information at these

thirteen points.

2.4 Taylor Models

Taylor model methods are newly developed techniques that unify many concepts
of high-order computational differentiation with verification approaches covering the
Taylor remainder term. These were developed as a better alternative to the conven—
tional interval methods that have limited practical applicability. The reasons for the
failure of conventional interval methods for large dimensional and domain size prob-
lems are discussed in [50, 42, 53]. The results obtained with Taylor model methods
include verified optimization, verified quadrature and verified propagation of extended
domains of initial conditions through ODEs, and approaches towards verified solution

of DAEs and PDEs.

Definition 1 (Taylor Model) Let f : D C IR” ——-> R be a function that is (n + 1)
times continuously partially difierentiable on an open set containing the domain D.
Let 170 be a point in D and P the n — th order Taylor polynomial of f around 1:0. Let

I be an interval such that
f(:r:) E P(:r — 230) + I for all :r E D.
Then the pair (P, I) is called an n — th order Taylor model of f around 2:0 on D.

32

For the problems we discuss in this work the domain D is always [—1, 1]” C IR”. A
full theory of Taylor model arithmetic for elementary operations, intrinsic functions,
initial value problems and functional inversion problems has been developed; see [50,
42, 53] and references therein. The arithmetic and the computational implementation
is performed in such a way that the interval enclosure is mathematically rigorous,
taking into account all round-off and threshold cut-off errors due to floating point
arithmetic. Details about the verified implementation of arithmetic operation in the
code COSY INFINITY can be found in [66, 53].

For the purposes of the further discussion, one particular intrinsic function, the
so—called antiderivation, plays an important role. We note that a Taylor model for
the integral with respect to variable i of a function f can be obtained from the Taylor
model (P, I) of the function by merely integrating the part Pn_1 of order up to order

n — 1 of the polynomial, and bounding the n-th order into the new remainder bound.

2.4.1 The Taylor Model integration scheme

We start out by defining the indefinite integration for a Taylor model in one variable
and then proceed to define a special finite integral that we are going to use extensively

in this work.
Definition 2 For an n-th order Taylor model T = (P, I) and k = 1,. . . ,v, let
37k
Qk=A P(n_1) ((131,...,33k_1,§k,33k+1,...,1‘v) dék.
The antiderivative 6;1 of T is defined by
_1 _
6k (P, I) _ (Qk, (B (PM — P(n_1)) + I) -2).
More details about the implementation of the anti-derivation operation can be

found in [22].

33

Let T” = (P, I) , be a n — th order Taylor model in v variables. While solving PDEs

we commonly encounter the finite integration

1
Ja—l) : /_1T”($1,...zrk_1.€k,l‘k+1,~-,17v)dék,

where J(v_1) is the resulting Taylor model of order n in (v — 1) variables, and k =
1,. . . ,v. We now describe the steps to compute the Taylor model J (”-4) using the

antiderivative 66—1 1 operator.

1. Split the Taylor Model T v = (P, I) in to a polynomial P of order n and interval
I.

2. Construct a new Taylor model G” = (P, 16), where 16 = [0,0].

3. Apply the antiderivative operator 65—1 1 to the Taylor model C”.

E 6

subtract them to obtain a new Taylor model R”_1.

4. Evaluate the Taylor models 611G” (at £1 = 1) and 6116'” (at £1 = —1) and

5. Add the interval 2 - I to the Taylor model R”"1 to give the Taylor model after

integration with respect to the £1 variable.

The steps 1 through 5 can be repeated for the integration in more than one variable.

2.4.2 An example of Taylor Model expansion

For a rectangular surface enclosing the volume of interest we now Show the Taylor
model expansion of the kernel 1/ (4n [F — Fl [) over one surface element and one vol-
ume element. A point, F = (:13, y, z), inside the volume element centered at ($0, yo, zO)
is described by

x = $0 + 0.5-A1Axv,

y = yo + 0.5 - AgAyv,

z = 20 + 0.5 - )3sz,

34

where A111), Ayv and sz are the dimensions of the rectangular volume element and
I I I
the parameters A1, A2, A3 6 [—1,1]. A point, F! = (a: ,y ,z ), on the surface element

centered at (1:3, ys, zs) is described by

IL‘ = (133 + 0.5 ' A4AI3,
I

y = ys + 0.5 . A5Ays,
I

Z : 287

where A273 and Ays are the dimensions of the rectangular surface element and the
parameters A4, A5 6 [—1, 1]. A schematic diagram of the rectangular volume element,
rectangular surface element and the rectangular box enclosing the region of interest
is shown in Figure 2.4.1.

For illustration we choose the volume element to be centered at (1, 1, 1) and surface
element at (2,2,3) and A2,; = Ays = 1/16 and Axv = Ayv = sz = 1/8. The
Taylor model expansion of the kernel with respect to the parameters A1, A2, A3, A4
and A5 using Taylor model is given in tables 2.4.1 and 2.4.2. In the representation
of the Taylor model expansion, the entries in the first column provide the number
assigned to each of the coefficients in the Taylor model expansion to easily identify
them. The entries in the second column provide the numerical value of the coefi‘icients.
The entries in the fourth through eighth provide the expansion orders with respect to
the parameters A1, A2, A3, A4 and A5. The total order for each coefficient is the sum
of all the orders in columns four through eight, which is given in the third column.

The Taylor Model integration of the kernel over the surface element, dS, can be

1
dd
hs4wrl_.=lly“’/_11/_141IIZI_

Table 2.4.3 shows the Taylor model integration of the Taylor model expansion of the

expressed as

WdA4dA5

 

 

kernel given in tables 2.4.1 and 2.4.2.

35

 

 

 

 

 

 

 

 

 

 

 

?
(XSYSZS)

/: L
/ /

/ /
/
, /
E /
E /
(XOeYOeZo) *4 ,
2* ______ T ....... _g .......... my ............ , Y
I!" i J /
I". .I
I".

Figure 2.4.1. The figure shows a volume element centered at (11:0, yo, 20) element and
a surface element (x5, y3, 23).

36

I COEFFICIENT ORDER EXPONENTS
1 0.2012846454507073E-01 0 0 0 0 O 0
2 0.1529273937765291E-03 1 1 O 0 O 0
3 0.1529273937765291E-03 1 O 1 O 0 0
4 0.2334154957641759E-03 1 0 0 1 O 0
5 -.3823184844413227E-04 1 O 0 0 1 0
6 -.3823184844413227E-04 1 0 O 0 0 1
7 -.7724385987298169E-06 2 2 0 O 0 O
8 0.3485629176768297E-05 2 1 1 0 0 0
9 -.7724385987298169E-06 2 0 2 O O O
10 0.5320170848751610E-05 2 1 O 1 O O
11 0.53201708487516108-05 2 0 1 1 O 0
12 0.1544877197459632E-05 2 0 O 2 O O
13 0.3862192993649084E-06 2 1 0 O 1 O
14 -.8714072941920742E-06 2 0 1 0 1 0
15 -.13300427121879033-05 2 O 0 1 1 0
16 -.8714072941920742E-06 2 1 0 O O 1
17 0.3862192993649084E-06 2 0 1 0 0 1
18 -.1330042712187903E-05 2 O 0 1 0 1
19 -.4827741242061355E-07 2 O 0 O 2 0
20 0.2178518235480185E-06 2 0 O O 1 1
21 -.4827741242061355E-O7 2 0 O O O 2

Table 2.4.1. The Taylor model expansion of the kernel. The coefficients up to second
order are shown.

37

22 -.3526083774525404E-07
23 0.8876341263195013E-08
24 0.8876341263195013E-08
25 -.3526083774525404E-07
26 0.1354809982277133E-07
27 0.20210181967026673-06
28 0.1354809982277133E-07
29 0.9690617197256710E-07
30 0.9690617197256710E-07
31 -.9032066548514218E-08
32 0.26445628308940538-07
33 -.4438170631597507E-08
34 -.2219085315798754E-08
35 -.6774049911385670E-08
36 -.5052545491756669E-07
37 -.2422654299314178E-07
38 -.2219085315798754E-08
39 -.4438170631597507E-08
40 0.26445628308940533-07
41 -.5052545491756669E-07
42 -.6774049911385670E-08
43 -.2422654299314178E-07

44 -.6611407077235133E-08
45 0.55477132894968832-09
46 0.84675623892320923-09
47 0.11095426578993773-08
48 0.1109542657899377E-08
49 0.1263136372939167E-07
50 0.55477132894968833-09
51 -.6611407077235133£-08
52 0.8467562389232092E-09
53 0.5509505897695944E-09

54 -.1386928322374221E-09
55 .1386928322374221E-09
56 0.55095058976959445-09

wwwwwmono.)wwwwwwwwwwwwwwwwwwwwwwwwwww
COOOOCHCOHOOHOOHOHMOOHOHMOOHOHMOHMW
OOOOOHOOHOOHOOHOMHOOHOMHOOHOMHOOJMHO
OOOOHOOHOOHOOMHHOOOMHHOOOQJMMHHHOOOO
OHNWOOOHHHMMMCOOOOCHHHHHHCOOOCOCOOO
wMHOMMMHHHOOOHHHHHHOOOOOOOOOOOOOOOO

 

 

VAR REFERENCE POINT DOMAIN INTERVAL
1 0.000000000000000 [-1.000000000000000 , 1.000000000000000 ]
2 0.000000000000000 [-1.000000000000000 , 1.000000000000000 ]
3 0.000000000000000 [-1.000000000000000 , 1.000000000000000 ]
4 0.000000000000000 {-1.000000000000000 , 1.000000000000000 ]
5 0.000000000000000 [-1.000000000000000 , 1.000000000000000 ]
REMAINDER BOUND INTERVAL
R [-.6787457395636954E-007,0.17693690345864713-006]

13******#**********1*******3|!******IIIIIUOHI********************

Table 2.4.2. The Taylor model expansion of the kernel. The third order coefficients,
the reference points and the remainder bound interval are shown.

38

 

I COEFFICIENT ORDER EXPONENTS

1 0.1965667222668859E-04 0 0 0 0 0 0

2 0.1493431579848917E-06 1 1 0 0 0 0

3 0.1493431579848917E-06 1 0 1 0 0 0

4 0.2279448200822031E-06 1 0 0 1 0 0

5 -.7543345690720868E-09 2 2 0 0 0 0

6 0.3403934742937790E-08 2 1 1 0 0 O

7 -.7543345690720868E-09 2 0 2 0 0 0

8 0.5195479344483994E-08 2 1 0 1 0 0

9 0.5195479344483994E-08 2 0 1 1 0 0

10 0.1508669138144172E-08 2 0 0 2 0 0

VAR REFERENCE POINT DOMAIN INTERVAL

1 0.000000000000000 [-1.000000000000000 , 1.000000000000000 ]
2 0.000000000000000 [-1.000000000000000 , 1.000000000000000 ]
3 0.000000000000000 [‘1.000000000000000 , 1.000000000000000 ]
4 0.000000000000000 [-1.000000000000000 , 1.000000000000000 J
5 0.000000000000000 [-1.000000000000000 , 1.000000000000000 I

REMAINDER BOUND INTERVAL

R [-.1628051142682157E-008,0.2054209443924109E-008]

t*t*#10!It*********It***1!!!18**********************************

Table 2.4.3. The Taylor model expansion of the kernel. The third order coefficients,
the reference points and the remainder bound interval are shown.

39

CHAPTER 3

3D Laplace and Poisson solver

using DA

In this chapter we describe a new method to solve the Poisson equation

V2¢ (7") p (F) in the bounded volume 0 C 1E3, (3.0.1)

V4607) = fl?) on the surface 60.

In comparison to the Neumann problem where only the component of the gradient
normal to the surface is specified, the boundary condition specified in the equation
(3.0.1) requires the full gradient to be specified on the surface. However, as we will see
later that the numerical scheme we develop requires all components of the gradient to
be specified on the boundary. In practice this poses no problem as we often measure
or have information about the complete field rather than just the normal component.

In the equation (3.0.1) the source density p (“r“) and the normal component of the

gradient 9' (“I") are related by

r 3x: fi-“I‘ 2m ..
fflpw [89 gm . (302)

where 13. is the unit vector perpendicular to the surface (99 at I". The equation (3.0.2)

states that the net field across the surface of the volume 9 is equal to the total

40

amount of field created at the sources inside the volume {2, and this follows directly
from application of the divergence theorem to the equation (3.0.1). The potential

(151 (7") due to the source density p (f) in vacuum is given by

._ 1 9(7’) 3
WWW/Md 1‘-

The potential (pl (7") satisfies the Poisson equation

 

V2¢1 (7") = p (7") in the bounded volume Q C IE3. (3.0.3)
Now, subtracting the equation (3.0.3) from (3.0.1), we get

v2 (a (F) — 31 (a) = 0 in the bounded volume a c 1E3,

V ((25 (F) — (251 (7’)) 9'0"”) — val (f) on the surface 89.

Let a new potential 1p (7") and a vector field ff) be defined by

1N7?) = (Mil—(film.
§(F)-V¢1(F)-

k.”
3
n

The potential 1b (7") is then the solution to the Laplace equation

V2¢ (F) = 0 in the bounded volume Q C IE3, (3.0.4)

th (F) = f0") on the surface 80.

We have succeeded in splitting our original problem of solving the Poisson equation
(3.0.1) into two independent problems, the solutions of which can then be combined
to get the final solution. In the first problem we find the scalar potential 1/2 (77') that
satisfies the Laplace equation inside the volume Q enclosed by surface 652, when
the gradient of the potential, th (F) = f0“) is specified on the surface 60. In the
second problem we compute the potential $1 (7") due to the source distribution p (f)
in vacuum. Both these elements can then be combined to get the solution to the

Poisson equation.

41

In the section 3.1.1 we first discuss the benefits of using the boundary data and
present the analytic closed form solution for the 2D case that can be easily found
by application of Cauchy’s integral formula. We then use a 2D example to highlight
the advantages of the methods that use the boundary data to compute the solution.
In section 3.1.3 we present the theory and the implementation of the new scheme
to find the solution of the 3D Laplace equation when the gradient of the solution
is specified on the surface enclosing the volume of interest . This scheme is based
on the Helmholtz theorem and the tools of the code COSY INFINITY [17, 23, 24].
In the section 3.1.3 we present an application of this new scheme to a theoretical
bar magnet problem. We also extend the theory to find the verified solutions of
the Laplace equation 3.2. The implementation and the results of the application to
the analytic bar magnet problem are discussed. The rest of this chapter describes a
new DA based technique to compute the magnetic field due to an arbitrary source

distribution.

3.1 The Laplace solver
The 3D Laplace equation
V21!) (7") = 0 in the bounded volume 0 C IE3 (3.1.1)

is one of the important PDEs of physics, describing among others the phenomenol-
ogy of electrostatics and magnetostatics. In many typical applications, not only the
normal derivative of ’l/J but indeed the entire gradient 62,0 is known on the surface;
for example, in the magnetostatic case the entire field E = 6gb is measured, and
not merely the component normal to the surface under consideration. Thus the field
computation problem can be viewed as solving a boundary value problem for the

three dimensional (3D) Laplace equation for the field, i.e. to obtain the solution of

42

the PDE

Vzip (f) = O in the volume 0 C IE3

where V2,!) (7") = f0“) is specified on the surface ('30.

As discussed in the chapter 2 the existence and uniqueness of the solution for the 3D
Laplace equation case can easily be shown through the application of Green’s formulae
or by using mean-value formulas. But, the analytic closed form solutions for the 3D
Laplace equation case can usually only be found for special problems with certain
regular geometries where a separation of variables can be performed. However, in
most practical 3D cases, numerical methods are the only way to proceed. Frequently
the finite difference or finite element approaches are used to find the approximations of
the solution on a set of points in the region of interest. But because of their relatively
low approximation order, for the problem of precise solution of PDEs, the methods
have very limited success because of the prohibitively large number of mesh points

required. Furthermore, direct validation of such methods is often very difficult.

3.1.1 Methods using boundary data

Boundary data methods such as those utilized below are based on a description of
the interior field in terms of particular surface integrals involving the surface data.
These approaches have various advantages. Firstly, the solution is analytic in terms
of the interior variables, even if the boundary data fail to be differentiable or are even
piecewise discontinuous; all such non~smoothness is removed after the integration is
executed. Hence a Taylor polynomial approximation in terms of interior variables
can be performed; and we expect that a Taylor approximation of a certain order will
provide an accurate approximation over suitable domains.

Secondly, since for the PDEs under consideration here the solution functions are

known to assume their extrema on the boundary because of analyticity or harmonic-

43

ity, a method that uses boundary data is expected to be robust against errors in
those boundary data with errors in the interior not exceeding the errors on the sur-
face. Thirdly, if the boundary data given have statistical errors, such errors have a
tendency to even average out in the integration process as long as the contributions
of individual pieces of integration are of similar significance. Thus we expect the
error in the computed field in the interior to be generally much smaller than the
error in the boundary data. This ensures that the methods using boundary data are

computationally stable .

3.1.2 The two dimensional case

As an introduction to the general approach, we begin with the discussion of the 2D
case, the theory of which can be fully developed in the framework of elementary
complex analysis, and which also describes the situation of static electric or magnetic
fields as long as no longitudinal field dependence is prasent. It is based on the use
of Cauchy’s integral formula stating that if the function f is analytic in a region

containing the closed path C, and if a is a point within C, then

f (a) = 5:72 0 £42 (3.1.2)

where the integral denotes the path integral over C. Cauchy’s formula is an integral
representation of f which permits us to compute f anywhere in the interior of C,
knowing only the value of f on C. This integral representation of f is also the
solution of the 2D Laplace equation for the primitive of (Re( f), - Im( f )) with the
function f specified on the path C.

Now, suppose a random error of 6 (z) is introduced in the measured data around
the path C. Then by the equation (3.1.2) we can compute the error E (a) introduced

in the computation of f (a) at some point 0 inside C as

E(a)——1— Mdz—f(a)— 1 f2 5(2) dz (3.1.3)

_27rz' C z—a _27rT z—a

 

44

We note that while E (a) is given by a Cauchy integral, E need not be analytic since
6(2) need not assume the function values of an analytic function. In fact, if it would,
then it already would be uniquely specified on any dense subset S of C, which removes
the freedom for all values of E on points on C that are not in S.

While the error E itself may be bounded in magnitude, if the integral is approx-
imated by one of the conventional numerical quadrature methods, the result can
become singular as the point a approaches the boundary C. This case may limit the
practical use of the method and needs to be studied carefully. As an example, we
consider the case of quadrature based on adding the terms of a Riemann sum, i.e.

the approximation

 

 

NZ 2,
533% 6(2) dz z —-1—.Z 5( 3),.» . (zj — zjnl) = E(a) (3,1,4)

z—a 2m j=1(2j‘0‘(

where the N 2 points zj are spaced equidistantly around C; since C is closed, 20 = z N 2'
By studying the approximation 13‘ (a) as the point a approaches the boundary C, we
can analyze the stability of the method with respect to the discretization of the path
C.

As an example, we choose the path C as a circle of radius R enclosing the region of
interest. We assume a random error of 6 (z) (i10‘2) is introduced in the measured
data around C. The point a is given by r - exp (2' - qfi) and the points zj are given
by R - exp (2' - 27rj/Nz) for j = 0,. . . ,Nz. Letting 6m(zj) denote the error assigned
to point zj in error set m, for each of these error sets we express the Riemann sum

Cm (a) for point a by

_ 1 N" izj5m(zj) 27r
Cm(O)—2—W2o§ (zj_a(r)) ON;

We then form the average of the magnitude of the error over Ne error sets to obtain

 

1 Ne
W“) = V Z lCm(a)|-
‘3 1

m:

45

Note that 17(r) still depends on the phase (1). However, in the statistical limit there
is apparently invariance under rotation by exp (i - 27r/Nz); and one quickly sees that
there are two limiting cases for the choice of the phase. These are the case (p = 0,
where the a will eventually collide with the zj for j = N z as r ——> R and thus a
”worst case” divergence will appear, and the case (1‘) = 27r/2Nz, in which case the
a will approach the mid point between zj for j = N z and zj for j = 1 as r -—> R.
Choosing sufficiently fine discretization of the path and sufficiently many error sets
6m, the quantity 17(7') for these two cases will be a good measure for the accuracy
that can be achieved with the surface integral method.

For our specific example, we choose random errors of maximum magnitude 10—2
at N z = 10,000 points on the circle of radius R = 2. For each value of r, we perform
the computation for a total of Ne = 10,000 error sets. The results of this analysis
are shown as plots in Figure 3.1.1a and Figure 3.1.1b for the two cases that represent
the ”worst case” and the ”best case” situation.

We first observe that sufficiently away from the surface, the expected smoothing
effect is happening, and the errors in the function values are indeed well below the
errors assumed on the surface. A rough quantitative analysis shows that this error is
about two orders of magnitude below the surface data error, corresponding well with
the statistically expected decrease of the error by 1 / m. As a approaches the curve
closer than 10-3, in the ”best case” situation, the error rises to about 10-2, which is
because now only nearby grid points contribute to the sum and thus the smoothing
effect disappears. In the ”worst case” scenario, divergence actually happens; but the
average error is still at the level of the original random error of 10—2 for values of I"
that are only about 10“4 away from the radius 2.

So overall we see that the method performs significant smoothing, and even with
the simplest discretization as a Riemann sum, good accuracy is maintained even as

we approach C. We note in passing that with more sophisticated quadrature meth—

46

 

4 " Worst Case —.— -
Best Case w... .....

L09( '0 (r))

 

u
in

 

 

 

Radius (r)
(a)

 

4 Wdrst Case I _._ d
Best Case "—4.--.

L09( 71 (I'))

 

 

 

 

 

Figure 3.1.1. (a) The plot shows the dependency of 77 (r) on the radius r. (b) The plot
shows the dependency of 17 (r) on the radius as the radius r approaches the boundary.
10, 000 error sets (Ne) around the circle or radius R = 2 were chosen for the analysis.
We show the plots for both the best and the worst case scenario.

47

ods, for example those based on Gaussian methods [17], the divergence effect can be

significantly controlled.

3.1.3 The three dimensional case

The scheme we use for the 3D case is based on the Helmholtz vector decomposition
theorem discussed in the sections 2.1.1 and 2.1.2. We begin by representing the
solution of the PDE via the Helmholtz theorem, which states that any vector field T3”

which vanishes at infinity can be written as the sum of two terms

 

 

E(a) = 6x “,(5)+v¢,,(r), where (3.1.5)
1 fi(")-—§(f) 1 63¢)
(anti?) = — S. - 8 ds——— —_,——43-dv, and
47r 5Q '13 — :csl 47r 9 la: — xv]
- 1 fi(f)x§(f) 1 vx'é’w)
Adi") = ““ 3.. _. S d8+-— —_.—_.—v—dV.
47r an |a: —— xsl 47r Q la: — xvl

Here 80 is the surface which bounds the volume 9. 533 denotes a point on the
surface 69, and it}, denotes a point within 0. ii is the unit vector perpendicular to
09 that points away from Q. 6 denotes the gradient with respect to 55v-

The first term is usually referred to as the solenoidal term, and the second term as
the irrotational term. Because of the apparent similarity of these two terms to the
well-known vector- and scalar potentials to R, we note that in the above representa-
tion, it is in general not possible to utilize only one of them; for a given problem, in
general both on and [It will be nonzero.

For the special case that R = 6V, we have P x R = O; furthermore, if V is a
solution of the Laplace equation 62V 2 0, we have P - If = 0. Thus in this case, all
the volume integral terms vanish, and 45,, (53') and 4t (55) are completely determined

from the normal and the tangential components of R on the surface 69 via

 

 

1 " “ If “
cm (a?) = 47 [an ”(72” 53l$8)d8 (3.1.6)
Km?) = —i ”(535) X B (maze. (3.1.7)

477 89 lf— fsl

48

For static electric or magnetic fields without sources in 9, which are characterized
by the Laplace problem that we are studying, the divergence and the curl of the
field vanish and hence these fields can be decomposed into irrotational and solenoidal
parts. For any point within the volume (2, the scalar and vector potentials depend
only on the field on the surface 89. And due to the smoothing properties of the
integral kernel, the interior fields will be analytic even if the field on the surface data

fails to be differentiable.

Surface integration and finite elements via DA

Since the expressions (3.1.6) and (3.1.7) are analytic, they can be expanded at least
locally. The idea is now to expand them to higher orders in BOTH the two components
of the surface variables :33 and the three components of the volume variables :3". The
polynomial dependence on the surface variables will be integrated over surface sub-
cells, which results in a highly accurate integration formula with an error order equal
to that of the expansion. The dependence on the volume variables will be retained,
which leads to a high order finite element method. By using sufficiently high order,
high accuracy can be achieved with a small number of surface elements, and more
importantly, a small number of volume elements. We describe the details of the
implementation in the following.

The volume 9 is subdivided into volume elements. Using the prescription for the
surface field, the Taylor expansion of the field is computed at the center of each volume
element. The final solution inside the overall volume is given as local expansions of
the field in different volume elements.

To find the local expansions for each volume element, we first split the domain of
integration 89 into smaller elements I‘i. Ffom the surface field formula we extract
an approximate Taylor expansion in the surface variables is about the center of the

surface element. Then the integral kernel 1/ IF —— PSI and the field B on the surface

49

are Taylor expanded in the surface variables 7"}; about the center of each surface
element. We also Taylor expand the kernel in the volume variables 7" about the center
of the volume element. The final step is to integrate and sum the resulting Taylor
expansions for all surface elements. Depending on the accuracy of the computation
needed we choose step sizes, order of expansion in 7' (:r, y, z), and order of expansion
in 1‘3 (:13, y, z).

All the mathematical operations to perform the expansion, surface integration,
curl and divergence were implemented using the high-order multivariate differential
algebraic tools available in the code COSY INFINITY [23, 24, 17] which automat-
ically leads to the respective field representation to any order without any manual

computations.

An analytical example: the bar magnet

As a reference problem we consider the magnetic field of an arrangement of the two
rectangular iron bars with inner surfaces (y = iyo) parallel to the midplane (y = 0 m)
as shown in Figure 3.1.2a. The interior of these uniformly magnetized bars, which are
assumed to be infinitely extended in the iy-directions, is defined by: 1131 S a: 3 3:2,
|y| 3 yo, and 2:1 3 2 S 22. From this bar magnet one can obtain an analytic solution

for the magnetic field B (:c, y, z) - see for example [31] - and the result is given by

2 P

B (any, z) z 0' E (_1)2+] arctan l _Z__]_ + arctan _74.
y 47r , + Y -R_
i7j=1 _ + ij - ij

2 ' 2- R7.

BO 2: ' ' .7 + 2]

Ba: (xaya 2) = _ (_1)2+] 1n +

47f i.j=1 _ 23' + Rz'j

2 - X- R."-

BO ' ' J + z]

B. (e.g. z) = — 2 3 (4)”? 1n ————-—+

4“ i,j=1 _ XI + Rij

 

1
where X,- =:c—:c,-, Yi =y0iy, Z7; =z—zz', and RE:- = (Xi2+YJ-2+Zi)2. We

note that because of the symmetry of the field around the midplane, only even order

50

BY

  
   
 

 
  
  

 

0.2 , {Am
4 ‘~'. ‘1
. ey'w 1;.» X W -
M . '
° 0' 3:33)? X ‘\
o T M Y K‘
. i?» "V51 I

 

(a) (b)

Figure 3.1.2. (a) The geometric layout of the bar magnet, consisting of two bars
of magnetized material. (b) The magnetic field By on the center plane of the bar
magnet. 30 = 1T and the interior of this magnet is defined by —0.5 g :r S 0.5,
lyl S 0.5, and -0.5 g z s 0.5.

terms exist in the Taylor expansion of this field about the origin. The mid plane field
of such a magnet is shown in Figure 3.1.2b.

The method presented in the section 3.1.3 is valid for any volume enclosed by a
smooth surface. For this example we choose three cases, (a) a volume enclosed by a
cube, (b) a volume enclosed by a cylinder, and (c) a volume enclosed by a sphere. For
these casas we first study the performance of the surface integration method. To this
end, we consider a cube of edge length 0.8 m, a cylinder of length 0.8 m and radius
0.4 m, and a sphere of radius 0.4 m. The center of all three geometries coincides with
the center of the interior of the uniformly magnetized bars. The six surfaces of the
cube are each subdivided into a 44 x 44 mesh. The surface of the cylinder and the
sphere are also subdivided so that the number of cells are same as in the cube case.
On each of the mesh cells, the contribution from the Helmholtz integral is expanded

using differential algebraic tools [17], and the resulting polynomial is integrated.

51

Figure 3.1.3 shows the accuracy of the predicted field, compared with the exact
solution, as a function of the order of expansion within the surface mesh cells. Results
are shown for the points (0, O, O), (0.1, 0.1,0.1), (0.2, 0.2, 0.2) and (0.3, 03,03) for the
cube case. It can be seen that at order six, an accuracy of approximately 10‘12 is
reached, which is very high compared to conventional numerical field solvers. The

corresponding results for the cylinder and sphere case are shown in figures 3.1.4 and

 

 

 

 

 

3.1.5.
'2 l 1 l l I I
Error at point (0.0,0.0,0.0) _._.
Error at point (0.1,0.1,0.1) m..-
4 iii)? ........... ‘3‘ Error at point (O.2,0.2,0.2) -...-. -
Error at point (0.3,0.3,0.3) —--~a .
-6 b
A -8 -
§
m
V -10 -
8’
_l
-12 -
-14 _
* -;- a:
‘16 1 1 l 1 1 1

Order

Figure 3.1.3. The error for the field calculated for the bar magnet example
with rectangular grid for individual points (0,0,0), (0.1,0.1,0.1) and (O.2,0.2,0.2)
(0.3, 0.3, 0.3).

We note that a change from an even order to the next higher order does not
produce significant change in the error, which is due to the specific symmetry of
the magnetic field and the resulting fact that even orders dominate in the Taylor
expansion. Also, this study clearly demonstrates that the method works for all smooth
surfaces. Henceforth we only present the results for the cube case.

For the next example, we split the volume inside the bar magnet into 4 x 4 x 4 finite

52

 

 

 

 

 

'2 l I 1 fi I I
Error at point (0.0,0.0,0.0) —+—
Error at point (0.1.0.1 ,0.1) —---—-*
-4 [e Error at point (O.2,0.2,0.2) ..... .
-6 e
A -8 -
§
LIJ
v -10 -
8’
_l
-12 _
-14 _
-16
2 3 4 5 6 7 8 9

Order

Figure 3.1.4. The error for the field calculated for the bar magnet example
with cylinderical grid for individual points (0,0,0), (0.1,0.1,0.1) and (O.2,0.2,0.2)
(0.25, 0.25, 0.25).

53

 

 

 

 

 

Error atrpoint (00.00.00) I l—o—
l “-31-:-"-'?-‘-?-rr- Error at point (0.1,0.1,0.1) -------><
Error at point (015,015,015) 1!
'6 ’ , Error at point (0.2, .2, .2) a ‘
-8 —
g -10 -
21,
8’
_J
-12 -
-14 _
-16 l 1 1 l 1 1
2 3 4 5 6 7 8 9

Figure 3.1.5. The error for the field calculated for the bar magnet example with

speherical grid for individual points (0,0,0), (0.1,0.1,0.1) and (015,015,015)
(0.2, 0.2, 0.2).

54

 

Log(RMS error)

 

 

 

 

-1o 1 1

Order

Figure 3.1.6. The average error for the field calculated for the bar magnet example
for finite elements of width 0.4 around points (0,0,0) and (01,01, 0.1).

elements of width :l:01 m. Within each of the elements, a Taylor expansion in the
three volume variables is carried out, resulting in a polynomial representation of the
field within the finite element cell. The polynomial representation is used to evaluate
the field at 1000 randomly chosen points within the cell, and comparing the result
with the analytical answer. Figure 3.1.6 shows the resulting RMS error for finite
elements centered around (0,0,0), (0.1,0.1,0.1), (O.2,0.2,0.2) and (03,03, 0.3). The
plot for the finite element centered at (0.3, 03,03) shows the behavior of the RMS
error as we approach the boundary. It can be seen that the method remains stable as
we approach the boundary. For the finite elements well within the volume of interest,
it can be seen that at order 7, an accuracy of approximately 10"6 is reached.

We see that the method of simultaneous surface and volume expansion, all of which
can be carried out fully automatically using differential algebraic tools [17] imple—

mented in the code COSY [23, 24], leads to accuracies that are significantly higher

55

than those of conventional finite element tools, even when unusually large finite ele-
ments are used.

For purposes of illustration, we now show in the table 3.1.1 the Taylor expansion of
the field given by the equation (3.1.6) and calculated using the DA tools of COSY over
one surface element for a particular point frozen inside the volume of interest. The
center of the surface element is at (—039, —0.39, 0.4) and the point is at (0.1, 01,01).
The surface element is described by (—0.39 + 0.5/\xx3, —0.39 + 0.5Ayy3, 0.4), where

Ax, Ay represent the length and width of the surface element and 2:3, ys E [—1, 1]. In

I COEFFICIENT ORDER EXPDNENTS

1 0.1430015055365947E-01 0 O O O 0 0
2 0.6922600731781813E-03 1 O O O 1 O
3 -.9437452710153340E-03 1 0 O 0 O 1
4 -.1561210105220474E-04 2 0 0 O 2 0
5 -.4471499751575185E-04 2 0 O O 1 1
2O -.3232493054085583E-07 5 O O O 1 4
21 0.6156849473575023E-O7 5 0 O 0 0 5
22 0.8960505971632865E-10 6 O O 0 6 0
23 0.1890553337467643E-08 6 0 O O 5 1
24 -.9792219471281489E-09 6 0 O 0 4 2

41 -.2417698920592542E-10 8 O O 0 4 4
42 0.7717865536738434E-10 8 O 0 O 3 5
43 -.2649803372019223E-11 8 O O O 2 6
44 -.2561415687161454E-10 8 0 O O 1 7
45 0.8506329051477273E-10 8 O 0 0 O 8

 

Table 3.1.1. A sample eighth order Taylor expansion in two surface variables.

the representation of the Taylor expansion in $3 and ys in the table 3.1.1, the entries
in the first column provide the number assigned to each of the coefficients in the

Taylor expansion to easily identify them. The entries in the second column provide

56

the numerical value of the coefficients. The entries in the fourth, fifth and the sixth
columns provide the expansion orders with respect to the volume variables(:r, y, 2).
And the entries in the seventh and eighth column provide the expansion orders with
respect to the surface variables ($3,313). The total order for each coefficient is the
sum of all the orders in columns four through eight, which is given in the third
column. Since we compute the Taylor expansion about a particular point (0.1, 01,01)
frozen in the volume of interest in two surface variables ($3,313), we notice that the
entries in column four, five, six are all zero. It can be seen that in this expansion,
the contributions of higher order terms depending on the surface variables decrease
rapidly, and thus the expansion shown would lead to a result of very high accuracy.

We now present the Taylor expansion of the contribution of the equation (3.1.6)
for one surface element and over one volume element inside the volume of interest
in the table 3.1.2. The center of the surface element is at (—0.39, —0.39,0.4) and
the center of the volume element is at (0.1, 0.1, 0.1). The surface element and the
volume element can be fully described by (—0.39 + 0.5)‘3233, —0.39 + 0.5Ayys,0.4)
and (01 -+- 0.5pxx, 0.1 + 0.5pyy, 0.1 + 0.5pzz), respectively, where Ax, Ay represent
the length and width of the surface element, and pm, py, pz represent the length,
width and height of the volume element, and x3,y3,:c,y,z 6 [—1,1]. In this case
the coefficients of the Taylor expansion depend on both the surface ($3,313) and the
volume variables (:5, y, z). The coefficients depending only on the surface variables
and the coefficient of the zeroth order term are same as in the previous example of the
expansion in just the surface variables. Once again we notice that the contributions of
higher order terms decrease rapidly for higher order, showing that also the expansion
in volume variables leads to a very accurate representation.

We now study the error dependency on the size(length) of the volume element, or
equivalently the number of volume elements chosen for the computation. For the order

of computation 3,5,7 and 9, Figure 3.1.7 and Figure 3.1.8 provide the dependence of

57

I COEFFICIENT ORDER EXPONENTS
1 0.1430015055365947E-01 O 0 0 O O 0
2 -.9590481459719686E-02 1 1 0 0 0 0
3 -.9590481459719686E-O2 1 O 1 O 0 O
4 -.9768082968233012E—O2 1 O 0 1 0 0
5 0.6922600731781813E-03 1 0 O O 1 O
6 -.9437452710153340E-03 1 0 O O 0 1
454 -.4509222359486833E-07 6 0 1 0 O 5
455 -.3067430813781439E-07 6 0 O 1 O 5
456 0.8960505971632865E-10 6 O O 0 6 0
457 0.1890553337467643E-08 6 0 O 0 5 1
458 -.9792219471281489E-09 6 0 0 O 4 2
1283 -.2417698920592547E-1O 8 O 0 0 4 4
1284 0.7717865536738462E-10 8 O 0 O 3 5
1285 -.2649803372019148E-11 8 O 0 O 2 6
1286 -.2561415687161455E-10 8 O 0 O 1 7
1287 0.8506329051477271E-10 8 O O O O 8

 

Table 3.1.2. A sample eighth order Taylor expansion in two surface variable and three
volume variables.

the average error on the length of the volume element and the total number of volume
elements. As an example, for cell lengths of 0.1, which leads to a total number of only
550 finite elements, an accuracy of 10’10 can be reached with a ninth order method.
Similarly, for a seventh order method with a cell length of 0.2, corresponding to 125
boxes, accuracies of about 10'6 can be reached. Compared to the conventional 3D
Laplace solvers which typically utilize in the order of 106 cells to achieve accuracies

in the order of 10—3, these results are quite promising.

58

 

Order 5 2---.»
Order 7 ..... 1...-
o r- Order 9 a 1

 

 

 

C __.x-—-x
O _ x. .)(---X"‘* *
3'45 **X'**.x* *lr’*"*. “
U) at ”kg”- a!“ " l‘ I!
2 x— x .x J!" I” 1‘ D
E Jr .*._ I” ”r - a” a ....... B 8
g -6 »- ‘ ,1”.- a“ B a ....... B _,
" 43 ......
" ”I (30 a
n ___....a
-8 1- pm “'8 a
.rBl‘
1’8
r
-10 E 1 1 1 1 1 1
0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45

Volume Element Length

Figure 3.1.7. The plot shows the dependency of the average error on the length of
the volume element.

3.2 Verified solution of the 3D Laplace equation

For various practical problems, very precise and verified solutions of this PDE are
required; but with conventional finite element or finite difference codes this is difficult
to achieve even without validation because of the need for an exceedingly fine mesh
which leads to often prohibitive CPU time. The method presented in the section
3.1 can used to build an alternative approach based on high-order quadrature and
a high-order finite element method to find a verified solution to Laplace’s equation.
Both of the ingredients become accessible through the use of Taylor model methods
[53, 50] and the corresponding tools in the code COSY INFINITY [24, 23]. The
solution in space is first represented as a Helmholtz integral over the two-dimensional
surface. The latter is executed by evaluating the kernel of the integral as a Taylor

model [53, 50] of both the two surface variables and the three volume variables inside

59

 

Order 3 —+—

OMmS m*m
0mm7 w*m
0 0mm9 Maw.“

 

 

 

 

E
a) '4 l
a) “k l” ''''
2 *— --------- x ---------------- x- -----
5 ---------------------- x
c» -6 .1- "“ q
3 .1 --------- g .............

_8 l— B G d

"“8"“
o 100 200 300 400 500 600

Number of volume elements

Figure 3.1.8. The plot shows the dependency of the average error on the number of
volume element.

the cell of interest. Finally, the integration over the surface variables is executed as a
mere polynomial integration, resulting in a local Taylor model of the solution within
one cell. The final solution is provided as a set of local Taylor models, each of which
represents an enclosure of a solution for a sub-box of the volume of interest. Examples
of the method and the precision that can be achieved will be given.

To obtain a verified solution to Laplace’s equation we start from equations 3.1.6 and
3.1.7. Using the fact that iff 75 is, we have 6(1/ If — it's-l) = — (if — f3) / If — isl3,
and similar relationships, it is possible to explicitly obtain the gradient of the scalar

potential, and with some more work the curl of the vector potential; the results have

60

the explicit form

 

 

van (:13) = _Z; 89 If— M3 3 (3.2.1)
__+
g g 1 (55— 58) x (ma-as) x B (553))
v x A, (a?) = —- f 3 d3 (3.2.2)
47f an |:z:'—:E3|

From the equation (3.1.5) we know that the field inside the volume of interest is just
a sum of the irrotational (equation (3.21)) and the solenoidal (equation (3.22)) part.
This is then the solution for the magnetic field as surface integrals. But to numerically
integrate the kernel and get the verified solution as the local Taylor model we need
a specialized numerical scheme. In the next section we introduce one such scheme
based on the Taylor models of the code COSY INFINITY [24, 23]. The anti-derivation
operation on the Taylor models discussed in section 2.4 will be extensively used in

implementation of the scheme.

Solution of the Helmholtz Problem using Taylor models

In the following, we develop a verified method based on Taylor model methods to
determine sharp enclosures of the field H and the potential 2,!) utilizing the Helmholtz
method.

Utilizing Taylor model arithmetic, the following algorithm now allows to solve the

Laplace equation for the Helmholtz problem.

1. Discretize the surface 69 into individual surface cells S,- with centers 3,; and the

volume 9 into volume cells V]- with centers vj.
2. Pick a volume cell Vj.

3. For each surface cell 82-, evaluate the integrands in the equation (3.2.1) and
(3.2.2), the so—called ”kernels”, in Taylor model arithmetic to obtain a Taylor
model representations in BOTH the surface variables of S,- AND the volume

variables of Vj, i.e. in a total of five variables.

61

4. Use the Taylor model anti-derivation operation twice to perform integration

over the surface variables of each cell Si-

5. Add up all results to obtain a three dimensional Taylor model enclosing the

field I? over the volume cell Vj.

6. If a verified enclosure of the potential «p to I? over the volume cell V]- is desired,

integrate the field B over any path using the anti-derivation operation.

As a result, for each of the volume cells Vj, Taylor model enclosures for the fields
I? and potentials 1,0 are obtained. All the mathematical operations to evaluate these
Taylor Models and surface integration are implemented using the Taylor Model tools
available in the code COSY INFINITY [24, 23].

Apparently the computational expense scales with the product of the number of
volume elements and the number of surface elements; of these, the number of volume
elements is more significant because of their larger number. In practice one observes
that when using high-order Taylor models, a rather small number of volume elements

is required, in particular compared to the situation in conventional field solvers.

3.2.1 An analytical example: the bar magnet

Once again, we consider the bar magnet example described in the section 3.1.3. Now
we compute the validate solution to Laplace’s equation with the boundary conditions

described in this example.

Results and analysis As a first step in the analysis of the influence of the dis-
cretization of the surface and volume on the result, we study the contributions of
the surface elements towards the remainder interval part of the total integral. The
volume expansion point is chosen as 7" = (.1, .1, .1) , and the size of the volume box

around it is chosen zero. Thus after the surface integration, the polynomial part of

62

the dependence on volume vanishes except for the constant term, and the accuracy
is only limited by the width of the surface element, which after integration over the
surface variables influences the width of the remainder bound. We plot the width
of the remainder interval versus surface element length for the scalar potential in
Figure 3.2.1. The center of the surface element is chosen as F3 = (.034, .011, .5). It
is observed that for high orders, the method quickly reaches an accuracy of around
10"16 for about 25 surface subdivisions, which correspond to about 210 z 1000 sur-
face element cells per surface. Under the assumption that each of these surface cells
brings a similar contribution, the accuracy due to the surface discretization will be

in the range of approximately 6 - 1000- 10"16 < 10"”.

 

0 I I I I I I I I I
Order 8 —+—
Order 7 -------><-
Order6 ----- *-- “A
_5 _Order 5 *3 ,,..~--;:.,-1|

Order 4 ---+--- -‘_..7;:'-;:,-,~5g
Order 3 -- "-0- -- -* ,..--'"“ ,..--",‘;"_,’,.-;:"".3'
Order 2 W43 ”Muir” .t ’.,-" .E"::"’I/

-10 ~ “j 'x ‘

LOG10(lnterval Width)
5:

 

 

 

 

-13 -12 -11 -10 -9 -8 -7 -6 -5 -4 -3
LOG2(Surface Element Length)

Figure 3.2.1. The remainder interval width versus the surface element length for
integration over a single surface element and vanishing volume size.

We now study the dependency of the polynomial part and width of the remainder
interval of the magnetic field on the volume element length. In all these plots the

surface element length is kept fixed at 1/128. Figure 3.2.2 shows the remainder

63

interval width for the 3; component of the magnetic field versus volume element lengths
for different orders of computation. The other components of the magnetic field

exhibit a similar behavior.

 

' ' l I I 1 fl r
Order8 -—+——
2 Order? ----»--- —
Order6 ----- .....
Order 5 a .4
0 Order4 "w. 4
order 3 -. ".O' .... __ u. 0. ..A‘ -" fl- -'.',’-.": -"’:.;‘."'
Order2 ----~¢---- A H . .7; 575.5%
§ . ”." .... , 3”” I!
a -4 _ ........... x z
E
2
E
5 -6 - .
2'?
_l
-8 - 4
-10 1

 

 

 

 

-7.5 -7 -6.5 -6 -5.5 -5 -4.5 -4 -3.5 -3
LOGZ(Vqume element length)

Figure 3.2.2. The remainder interval width versus the length of volume element for
y component of the magnetic field.

We see that a verified accuracy in the range of 10‘4 can be achieved for a vol-
ume element width of around 10-1, corresponding to a total of around 1000 volume
elements. This number compares very favorably to the above-mentioned numbers
for the commercial code TOSCA [3, 4]. An accuracy in the range of 10"7 can be
achieved for a width of around 10—1'4, corresponding to a total of around 200,000
volume elements.

Overall, we see that the method of simultaneous surface and volume expansion
of the Helmholtz integrals leads to verified tools for the solutions of PDEs which
when executed in the Taylor model arithmetic can lead to very sharp enclosures. It is

obvious that the method can be generalized to other surface-integral based approaches

64

to the solution of PDEs.

3.3 Parallel implementation of the Laplace solver

Recently the code COSY INFINITY has been ported to support large scale computa-
tions for beam dynamics simulation and design optimization, including verified global
optimization. Test results with some of the COSY tools adapted to parallel execu—
tion have shown that they scale linearly to about 1000 processors [46]. Along this
line of high performance computing efforts with COSY INFINITY, a parallel imple-
mentation of the Laplace solver has been practiced on the NERSC (National Energy
Rwearch Scientific Computing Center) IBM RS6000 Seaborg Cluster consisting of
6080 processors [2].

In our implementation of the Laplace solver, we first discretize the surface enclosing
the volume of interest into surface elements, and then the magnetic field contribu-
tion of each surface element at a given observation point or volume is computed
independently. We then sum up the magnetic field contributions of all the surface
elements to obtain the magnetic field at the observation point or volume. The large
summation over all the surface elements can be trivially parallelized. However, if
trivially parallelized, the large summation over all the surface elements at the end
may take significant time and render the algorithm inefficient. Also, to utilize the
computational resources productively and minimize the cross communication between
processors, which may potentially slow down the computation, we need to contrive
an efficient algorithm.

An efficient parallel algorithm to some extent depends on the architecture of the
cluster that we use. For example, the Seaborg cluster, on which we implement our
parallel algorithm, has 380 computing nodes with each node having 16 processors.

Processors on each node have a shared memory pool of 16 to 64 GBytes. The commu-

65

nication between the processors within a node is much faster than the communication
between the processors of different nodes. Hence, care must be taken in designing an
algorithm to minimize the communication between the processors of different nodes.
Also, the Seaborg cluster supports several optimization modes that further improve
the performance of the algorithm. The detailed discussion is beyond the scope of this
work, but can be found in [2].

We now highlight some of the aspects of the parallel algorithms for the Laplace
solver. Let NPR be the number of processors we choose to use on the cluster for
parallel execution, and let N SP be the number of surface elements over which the
magnetic field is specified. The master processor will divide the computation of N S P
surface points over NPR processors. The master processor itself is a part of the
NPR processors.

If the summation over all surface elements is trivially parallelized, each processor
computes the partial sum of the magnetic field contributions for the assigned surface
elements, and sends the result of the partial sum back to the master processor, where
these partial sums are summed up to get the total magnetic field. The number of
communications between processors is equal to k x NPR, where k is a constant
depending upon the problem.

We now describe a new efficient algorithm for parallelization of the summation over
all the surface elements. The master processor will divide the computation of N S P
surface points over NPR processors. The NPR processors are split into groups of N 1
processors each, leading to N2 = N PR/N 1 groups. The master processor assigns
one of the processors in a group as the master processor for the group, and these
processors are referred to as sub—masters. The master processor itself is a part of the
NPR processors and also a part of the N2 sub-masters. The number of processors

N1 in each group is given by

N1 = INT (2- NPR) ,

66

where the INT function finds nearest integer smaller than are equal to (2 - \/N_PR).
However, if NPR is not exactly divisible by N 1, then we decrease the value of N1
by one, N1 = (N 1 — 1). We repeat the process till NPR is exactly divisible by N 1.
This ensures that both N 1 and N2 are integers. Each processor computes the partial
sum of the magnetic field contributions for the assigned surface elements, and sends
the computed result of the partial sum back to the sub-master processor. Each sub-
master processor computes the partial sum of the magnetic field contributions for each
group and sends the result of the groups partial sum back to the master processor,
where these group partial sums are summed up to get the final field. The number of
communications between the processors is roughly equal to k x (N 1 + N2). Since,
(N 1 + N2) << N1 x N2, the communication time is greatly reduced as compared to
the trivial parallelization case.

The above algorithm is implemented using the parallel loop block (PLOOP) avail-
able in the code COSY INFINITY. The steps are highlighted in Table 3.3.1. We
use two nested parallel loop blocks, one over the N2 groups and the other over N1
processors in each group. At the end of a parallel loop a communication option can
be specified to gather computed results from all the processors that have participated
in the parallel loop block on to only one processor, which once again minimizes the
cross communication between the processors. In the section 4.3.4 we present exam-
ple of the parallel implementation to compute the magnetic field for the Super-FRS

quadrupole magnets.

3.4 Magnetic field due to arbitrary current distri-
bution

In this section we describe a method to compute the magnetic field of an arbitrary

current distribution using DA techniques. The motivation to look at this problem

67

{Loop over N2 groups}
PLOOP JJ 1 N2;

{Loop over N1 processors}
PLOOP II 1 N1;

{Evaluate the processor number PP}
PP:=II+(JJ-l)*N1;

[Code to identify the surface elements JBEG through JEND
for which the processor PP will evaluate the partial sum]

{Loop to compute the partial sum of the scalar and vector potential
contributions over surface elements JBEG through JEND}
LOOP IL JBEG JEND;

[Code to compute the scalar and the vector potential
contribution of a surface element IL.]

ENDLOOP;

{End the parallel loop over the group of N1 processors and
send the results to sub-master processor using communication mode 4}
ENDPLOOP 4 PN1_SCLPOT PN1_VECPOT;

{Loop to evaluate group partial sum of N1 processors}
LOOP II 1 N1;

[Summation to get group partial sum GN2_ SCLPOT and GN2_ VECPOT]

ENDLOOP;

{End the parallel loop over the N2 groups and send the results to master processor}
ENDPLOOP 4 GN2_SCLPOT GN2_VECPOT;

{Loop to evaluate sum over N2 groups}
LOOP JJ 1 N2;

[Summation to get sum SCLPOT and VECPOT]

ENDLOOP;

[Code to evaluate the divergence of SCLPOT and the curl of VECPOT and sum them
to get the magnetic field]
Table 3.3.1. The code for the parallel algorithm.

68

comes from the need to design and optimize an accelerator magnet in a simple and in
an efficient way. Once again the DA based techniques have the advantage of providing
the complete multipole decomposition of the field. It is also straight-forward to verify
the that the field computed in the source free region indeed has identically vanish-
ing divergence and curl. Hence, it is guaranteed that the field computed is always
Maxwellian. Previous attempts to compute the magnetic field using DA techniques

are described in [34, 68].

3.4.1 Field computations using the Biot-Savart law and DA

We first describe a general frame work that uses the Biot-Savart law to compute the
field due to an arbitrary current distribution. The Biot-Savart law for line, surface

and volume currents is given by

 

 

B(m ~_—. —/ 7‘ "‘ dl, (3.4.1)
4” l (r—f’)

B(F) = _/ 7‘27‘ da, (3.4.2)
4" 5 F—f’)

 

“ _ fig F—r
B(F) _ 4% /V dr, (3.4.3)

where f is a current vector, and the surface current density vector 1? is the current per
unit width perpendicular to the current flow, and the volume current density vector
j is the current per unit area perpendicular to the flow of the current. The vector
F— 7"! points from the current element at f, to the observation point r' where we want
to compute the magnetic field. The vector 1272'. 7:; represents the unit vector in this
direction. The magnetic constant, #0: is the permeability of vacuum. In SI units,

the value is exactly expressed by #0 = 47r x 10‘7 NA’2. In most cases numerical

integration is usually required to find the total magnetic field at any point.

69

The general idea is to discretize the domain and express the integral over the domain
as the sum of integrals over smaller intervals. Depending on the type of the problem,
we express the domain in terms of one, two or three parameters and scale them such
that the new domain is a box [—1,1]"’, where n is number of parameters. We then
Taylor expand the kernel in terms of the previously defined parameters and integrate
it. We first present two analytical examples to describe the DA approach to obtain
the multipole decomposition of the magnetic field for the line and surface currents.
And, in the section 3.4.4 we describe the technique to find the magnetic field for the

volume current distribution.

3.4.2 Line current example: circular loop

For a circular loop with its axis oriented in the 2 direction, the analytic formula for

the magnetic field on the z—axis is given by

—. 2 21
B(0,0,z) = g—i—Va
(22 + R2) 2
We choose an example with radius R = 0.4 m and current of I = 1 A. For this

example, the eighth order expansion of the 2 component of the magnetic field at z = 0

is given in Table 3.4.1.

I COEFFICIENT ORDER EXPDNENTS
1 0.1570796326794896E-05 O O 0 0
2 -.1472621556370215E-04 2 0 O 2
3 0.1150485590914231E-03 4 0 O 4
4 -.8388957433749600E—03 6 0 O 6
5 0.5898485695605188E-O2 8 O O 8

 

Table 3.4.1. A eighth order Taylor expansion of the analytic formula for the 2 com-
ponent of the magnetic field of a circular coil on the central axis.

70

In the representation of the Taylor expansion above, the entries in the first column
provide the number assigned to each of the coefficients in the Taylor expansion to
easily identify them. The entries in the second column provide the numerical value
of the coefficients. The entries in the fourth, fifth and the sixth columns provide the
expansion orders with respect to the observation point (:r, y, z). The total order for
each coefficient is the sum of all the orders in columns four through eight, which is
given in the third column.

We now use a DA based approach to solve this problem. Let the point r’ = (:r, y, 2)
be the observation point. We discretize the length of the current loop into N = 4000
current elements of length (27rR) /N. Let (can, yn, 2n) describe a point inside the nth

current element, where n = 1, . . . , N. And let 3 be a parameter such that

xn = R - cos (A3 - (n - 0.5 + 0.55)), (3.4.4)
yn = R - sin (A3 - (n — 0.5 + 0.53)),

211:0,

where A3 = 27r/N and s E [—1, 1]. By varying s we can now get all points inside any
given current element which is centered at fl, = (R - cos (n - A3) , R - sin (n - A3) ,0).
We can now compute the contribution due to the current element at r; using the
Biot-Savart law for a line current, the equation (3.4.1), by first expressing the inte-
gral in terms of the new parameter s using the equation (3.4.4), and then expanding
the kernel in terms of the variables 31,3], 2 and s. We then integrate the resulting
polynomial with respect to s in the interval [—1, 1]. Finally, we sum up the contribu-
tion due to all the current elements to get the resulting field at (:1), y, z) .The Taylor
expansion of the :2 component of the magnetic field, E (:13, y, z) , computed using the

DA framework available in the code COSY INFINITY is given in Table 3.4.2.

Note that the expansion given in Table 3.4.2 reduces to the expansion for the

71

I COEFFICIENT ORDER EXPONENTS
1 0.1570796326794921E-05 0 0 0 O 0
2 0.7363107781851075E-05 2 2 0 O O
3 0.7363107781851056E-05 2 0 2 0 O
4 -.1472621556370156E-04 2 O O 2 O
5 0.4314320965821285E-04 4 4 O O 0
6 0.8628641931856381E-04 4 2 2 O O
7 0.4314320965821278E-04 4 0 4 0 0
8 -.3451456772742676E-03 4 2 O 2 O
9 -.3451456772742692E-03 4 0 2 2 O
10 0.1150485590914342E-03 4 O O 4 0
11 0.2621549198046754E-O3 6 6 0 0 O
12 0.7864647594140242E-03 6 4 2 O 0
13 0.7864647594140274E-03 6 2 4 0 O
14 0.2621549198046748E-03 6 O 6 O 0
15 -.4718788556484144E-02 6 4 O 2 0
16 -.9437577112968266E-02 6 2 2 2 O
17 -.4718788556484138E-02 6 0 4 2 0
18 0.6291718075312196E-02 6 2 0 4 O
19 0.6291718075312182E-02 6 O 2 4 O
20 -.8388957433749353E-03 6 0 0 6 O

 

Table 3.4.2. A sixth order Taylor expansion of the 2 component of the magnetic field
computed for the circular loop example.

analytic case, Table 3.4.1, when .r = 0,y = O. The result that we get is valid as
long as the observation point does not lie on the current loop. This case will require
special treatment. It is also easy to verify that the magnetic field computed here is
divergence less and curl free.

The DA technique that we develop here can easily be generalized to compute the
magnetic field due to the current on a wire of arbitrary shape and orientation in
space. We only need the discretization information and the parametrization over

each current element to compute the magnetic field.

72

 

 

 

 

 

 

 

 

 

 

 

</ f
RdeyJ ; W:
I:
\ : /
\ > J
4- ------------- b

Figure 3.4.1. The cross section of a flux ball consisting of a sphere with winding on
its surface.

3.4.3 Surface current example: spherical coil

For the case of a properly wound spherical coil the magnetic field inside has the
property that it is uniform. The current density on the surface can be thought of as
a sinusoidally distribution between the north and south poles of a sphere. Current
loops of appropriately varying diameter, spaced evenly as projected onto the .2 axis,
automatically simulate such a distribution. The coil, with a radius R and a wire
carrying the current I , is shown in Figure 3.4.1.

To deduce the surface current density representing this winding, note that the
density of turns on the surface is the total number, N, divided by the total length,
2R, and so the number of turns in the incremental length dz is (N /2R) dz. Since
z = rcosH , a differential length dz corresponds to an angular increment d6/ dz =
--R sin 0. Therefore, the number of turns in the differential length RdB as measured
along the periphery of the sphere is (N / 2R) sin 6. With each turn carrying the current

I, the surface current density is

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 3.4.2. Magnetic field lines for the spherical coil.

The magnetic field intensity for such a surface current distribution is shown in
Figure 3.4.2. The field inside the sphere is constant.

The magnetic field inside and outside the sphere is as follows.

R = #03}; (cos (6’) if — sin (0) i9): fit—g évéz, r < R
a N1 3
B = MOISR (g) (2 cos (0) ET + sin (6]) i9), r > R

The exterior lines of the magnetic field intensity are those of a dipole, while the
interior field is uniform. Figure 3.4.2, shows the magnetic field lines.

We consider a specific case where we choose the radius of sphere R = 0.4 m and the
total current, NI, such that 21%- = 1000 A / m. The analytic formula then gives the
2 component of the field inside the sphere to be g7!" x 10—4 = 0.8377580409572781 x
10’03 Tesla. We now use a DA based approach to solve this problem. Let the
point 1” = (:13,y, 2) be the observation point. Let N2; = Ny = 88 be the number of
subdivisions in azimuthal 6 and polar ¢ angles, leading to N7; x Ny surface elements.

Let (X,- j, Y, j, Z,- j) describe a point inside the (i, j)-th surface element and Kid

74

describe the surface current on the (i, j)-th surface element, where i = 1, . . . , N3; and
j = 1,... ,Ny. The point (R, 623%) in spherical coordinates is equivalent to the
point (XM’ YM’ ZN) in the Cartesian coordinates The surface elements and the

current at (i, j) can now be described in terms of two parameters, 1:3 and ys, such

that

<25.- —-— Ax - (2' — 0.5 + 06:22.9) (3.4.5)
Xi,j = R - COS ($2) - sin (Qj)

YiJ = R-sin(¢i)'5in(6j)

Zi,j = R - cos (Gj)
_. NI , . . .
Kid 2 ER ~s1n (Oj) -(—— 8111 (925i) :5 + cos (45,-) y),

where Ax = 27r/Nx, A3, = 27r/Ny and :cs,ys E [—1,1]. By varying 11:3,y3 we can now
get all points inside any given surface element. We can now compute the contribution
due to the surface element at if, using the Biot-Savart law for the surface current, the
equation (3.4.1), by first expressing the integral in terms of the new parameters 1:3
and ys using the equation (3.4.5), and then expanding the kernel in terms of the three
volume variable as, y, z, and the two surface variables ass and y3. We then integrate
the resulting polynomial with respect to (1:3,y3) over the domain [—1, 1] x [—1,1].
Finally, we sum up the contribution due to all the surface elements to get the resulting
field at (:r, y, z) . The 2 component of the magnetic field, R (:13, y, z) , computed using
the DA framework available in the code COSY INFINITY is given in Table 3.4.3
Note that the zeroth order term in the Taylor expansion given in Table 3.4.3 agrees
with the values computed using the analytic formula. The result that we get is valid
as long as the observation point does not lie on the surface. This case will require
special treatment. It is also easy to verify that the magnetic field computed here is

divergence less and curl free.

75

MAGNETIC FIELD (Bx,By,BZ)

0.000000000000 0.000000000000 0.8377580409801E-03 00000
0.000000000000 0.000000000000 -0.7490871758478E-12 20000
0.000000000000 0.000000000000 -0.7490871081318E-12 02000
-0.1496412647114E-11 0.000000000000 0.000000000000 10100
0.000000000000 -0.1496412522919E-11 0.000000000000 01100
0.000000000000 0.000000000000 0.1496363218034E-11 00200
0.000000000000 0.000000000000 0.1406498322650E-08 40000
0.000000000000 0.000000000000 0.2812995931364E-08 22000
0.000000000000 0.000000000000 0.1406498327655E-08 04000
0.5625993567287E-08 0.000000000000 0.000000000000 30100
0.000000000000 0.5625991024637E-08 0.000000000000 21100
0.5625991025586E-08 0.000000000000 0.000000000000 12100
0.000000000000 0.5625993571062E-08 0.000000000000 03100
0.000000000000 0.000000000000 -0.1125198581064E-07 20200
0.000000000000 0.000000000000 -0.1125198581536E-07 02200
-0.7501323883764E-08 0.000000000000 0.000000000000 10300
0.000000000000 -0.7501323883789E-08 0.000000000000 01300
0.000000000000 0.000000000000 0.3750661939904E-08 00400
0.1072194029532E-15 0.000000000000 0.000000000000 31100

 

Table 3.4.3. A fourth order Taylor expansion of the 2 component of the magnetic
field computed for the spherical coil example.

Once again, the DA technique that we develop here can easily be generalized to
compute the magnetic field due to an arbitrary surface current distribution. We only
need the discretization information and the parametrization over each surface element

to compute the magnetic field.

3.4.4 Three dimensional current distribution

We discuss the three-dimensional case by using the example of the current in the
straight wire with a finite cross section. This particular case is of practical importance
and we will use this to develop tools to design new accelerator magnets. These tools

will then be used to design a 3D model of a quadrupole magnet in the section 4.1.

76

Magnetic field computation for a wire with a rectangular cross section

using DA

We first develop a framework to express a point inside a finite length current wire of
a rectangular cross section that is oriented in an arbitrary direction in the 3D space.
This would then facilitate the implementation of the Biot-Savart law and Ampere’s
law in the differential algebraic (DA) framework.

Let a1 and fig be two unit vectors on a plane (7‘21 74 712), then any point on the plane
can be given as A161 + A262, where the parameters A1 and A2 are any real numbers.
The unit vector a} = 73.1 x 1‘22 along with two points on the plane completely defines
the plane formed by the unit vectors hl and 13.2. In the special case where 1‘21 is
perpendicular to 712 (n1 432 = 0) and A1, A2 6 [—1,1], the vector (Alfil + A2732) /2
describes a point inside a unit square lying on the plane and centered at the origin.
We can translate and scale this unit box to describe a rectangle lying on this plane.

The vector 17p defined by
VP = 17¢ + 0.5-(A1b-fi1+ /\2w . 712), (3.4.6)

describes a point inside a rectangle of the length b and width 11) and centered at
a point 176. All points inside the rectangle are completely described by the two
parameters (A1, A2) and are independent of the actual dimensions of the rectangle
and its orientation in the three dimensional space. In other words, we have projected
the area of the rectangle oriented in an arbitrary direction in the three dimensional
space on to a square described by the parameters (A1, A2).

For the case of a rectangular box of height I, let the cross section be defined by the
unit vectors in and ’ft2. The center point of the rectangular cross section lies on a
straight line Whose direction is given by the unit vector 71 1' A point on this line can
be expressed by

VG: VCO+/\3'fii,

77

where the vector 1760 is the center of a face on the rectangular box that is perpendicular
to the unit vector iii, and A3 6 [0,1]. The parameters (A1, A2, A3) then completely

describe a rectangular box. Any point inside the box is given by
9,?” = 1760 + 0.5-(A1b-n1+ A211) 432) + A341,]. (3.4.7)

Once again, we have projected the volume of the rectangular box oriented in an
arbitrary direction on to a box defined by the parameters (A1, A2, A3) . The equations
(3.4.6) and (3.4.7) will be used in the next two sections to express the Biot-Savart

law and Ampere’s law in the DA framework.
Magnetic field of a wire with a rectangular cross section using the Biot-
Savart law and DA

We can choose the direction of the current in the wire to coincide with the unit vector
n 1' Using the Biot—Savart law, we can write the magnetic field at an observation

point F as

-'_ "boa:
[3(5) —"—‘9 I-—0- 111/31 [)(ni:(r VP, ))-di3[d,\2-d/\1. (3.4.8)

 

47rb w

 

 

To perform integration with respect to A3 in the equation (3.4.8), we first split the
domain of integration into smaller intervals. Let the length I be divided into N parts

of the size h = l / N . Then the parameter A3 can be written as
A3 (i) = (i +0.5 - v) - h,

where i = 0.5,. . . , (N — 0.5), and v 6 [—1,1]. The position of a point inside the box

in terms of the new parameters (A1, A2, v) is given by

171903: (i, A1, A2,'U) = VCO + if? - 7A2] + 0.5-(Alb-fl1-l- A211) - fig + Uh 41]). (3.4.9)

78

The equation (3.4.8) can now be written as

44.4. 3°54 / 12“”

b 2
4.=05 :—(V],:0$(i)[

 

 

(3.4.10)
The parameters (A 1, A2, v) become the DA variables with respect to which the kernel
of the integral in the equation (3.4.10) is expanded to a high order. In the DA
framework it is straightforward to perform the volume integral over the resulting
polynomial representation.

In the DA frame work it is also straightforward to treat the wires with rectangular
cross sections that have curved endings rather than the straight endings that are
perpendicular to the direction of the current flow. Let the two surfaces at the start
and the end of the wire carrying the current be expressed as A3 = g(A1,A2) and
A3 = f (A1, A2). We can use the equation (3.4.10) to find the magnetic field due to

this new configuration by noting

 

“(132) : f(/\1,/\2)1;9(/\14/\2)’

where the step h is now a function of the parameters (A1, A2). The equation (3.4.9)

can be modified to express a point inside the rectangular box with curved endings as

+0.5-(A1b°ft1 +A2w-i12+v-h(A1,A2) “[3]).
When the surfaces are just inclined planes, the functions f and g are just linear
combinations of the parameters A1 and A2. This special case is useful in the imple-

mentation of the numerical tool to compute the magnetic field due to a coil carrying

a current with a rectangular cross section .

79

Magnetic field of an infinitely long wire with a rectangular cross section

using Ampere’s law and DA

Once again, we can choose the direction of the current in the infinitely long wire with
a rectangular cross section to coincide with the unit vector iii. The closest distance
r 1 between the observation point r' and the line passing through the point 17p and in

the direction fl] is given by

44(44)-(NF-nun,

where 17p is given by the equation (3.4.6). Ampere’s law can then be written as

B(=4) /11/_ #_0 Io ("ix((T‘Vpl/l(’"‘VP)I))d,l3,2 (3.4.1,)

127rbw r]

The equation (3.4.11) can be used to compute the magnetic field of an infinitely
long wire with a rectangular cross section . The parameters (A1, A2) become the DA
variables with respect to which the kernel of the integral in the equation (3.4.11) is ex-
panded to high order, and DA integration is performed over the resulting polynomial

representation.

COSY INFINITY tools for magnetic field computations

Due to their frequent use in the magnet design, a dedicated set of tools has been
written for the rectangular cross section wire and coil in the code COSY INFINITY
[23, 24]. These tools use the differential algebraic framework available in COSY IN-
F INITY [23, 24] to Taylor expand, integrate and evaluate the kernels appearing in
the equations (3.4.11) and (3.4.10). In addition to providing highly accurate results
in the form of the local Taylor expansion of the magnetic field, the DA based im-
plementation has a unique advantage of easily obtaining the curl and divergence of
the magnetic field at any given point. This offers one way to quickly verify if the

magnetic field satisfies Maxwell’s equations.

80

One of the tools computes the magnetic field of a finite wire of a rectangular cross
section. A finite wire can have the current entrance and exit planes inclined to the
central axis or the direction of the current flow. Also, a tool to compute the field for
an infinitely long wire with a rectangular cross section has been implemented. Using
orders around 10, accuracy of about 14 digits can be achieved using these tools.
Another tool to compute the magnetic field of a current coil of a rectangular cross
section has also been implemented. More detail about the current coil implementation
will be given in the section 4.1.2.

In general, a numerical tool to compute the magnetic field for an infinitely often
differentiable current distribution can also be implemented using concepts above and

the DA framework available in the code COSY INFINITY [23, 24].

3.5 Extraction of transfer map from surface field
data or current distribution

In the sections 3.1 and 3.4 we have successfully demonstrated that we can obtain local
expansion of the magnetic field by either using the surface data or the current distri-
bution. We now use the technique described in section 2.3.2 to demonstrate that we
can successfully extract transfer map from the surface field data. The field informa-
tion at 13 point of each time step in the Runge—Kutta integrator is supplied as local
expansion of the magnetic field computed using the techniques we have developed.
We present the results for a test case of theoretical quadrupole magnetic field for
which the analytic transfer map is know. The quadrupole field varies linearly with the
distance from the magnet center, R (113, y, s) = (qu, kqx, 0), where kq is a constant and
has units of T/ m. The quadrupole magnet focuses the beam along one plane while
defocusing the beam along the other plane. We note that the radius of curvature,

h = 0, for the quadrupole magnet. The equation (2.3.3) can now be linearized to

81

obtain

a: = a,
y' = b.
I k
a = ———q-a:,
XmO
I k
b _ _‘Ly,
XmO
I k 1

 

Let w :2 i/kq/XmOa and D = —k/(vo(1+770)(2+170)) and L be length of the

quadrupole. It is straight forward to solve the above system of ODEs to obtain the

first order transfer map for the quadrupole, which is given as

 

 

cos (wL) sin (wL) /w 0 0 0 0 \
—w sin (wL) cos (wL) 0 0 0 0
M = 0 0 cosh (wL) sinh (wL) /w 0 0
0 0 w sinh (wL) cosh (wL) 0 0
0 0 O 0 I D
K 0 0 0 0 0 1 )

It is also possible to obtain similar analytic formulae for the high order coefficients
[27].

We compare the transfer map computed using the analytical formulas and the trans-
fer map computed using Runge—Kutta integrator with the magnetic field computed
using the Laplace solver we have developed. For the case a quadrupole of length 0.2
m, 0.1 m aperture and field strength at the aperture of 6.725 x 10‘2 T we computed
the difference between the transfer maps. The result is presented in tables 3.5.1 and
3.5.2. We can see that the agreement between the analytic and computed Taylor maps
is very good, and this once again demonstrate that the magnetic field generated by
the techniques we have developed are highly accurate. This examples also serves and

a proof that we can extract transfer maps from the measured surface field data using

82

the techniques we have developed. In the section 4.3.4 we will present the example of

the Super-FRS quadrupole magnet for which we will extract the transfer map using

the magnetic field data on a closed surface.

3.6 Summary

.' Measured Field Data:
0. AND/OR I
' Source Distribution:

1

 

 

 

 

 

 

 

 

 

 

 

 

 

B Ph c Helmholtz Theorem Differential
“OmDE'YSi s AND/OR Algebraic
s Biot Savart Law Framework
- - ODE Integrators Magnetic AND/OR
D’i‘fferinual (Runge- -Kutta) Electric Field
ge ra DA Integrators Information
Taylor Map

Figure 3.6.1. The flow chart for extracting transfer maps from the measured magnetic
field data or the source distribution.

A new technique for finding the multipole expansion solution of the three dimen-
sional Poisson equation using the surface data and the current distribution inside the
volume enclosed by the surface has been developed. Since this new technique uses the
field information on the surface enclosing the volume of interest and is implemented
using the high-order multivariate differential algebraic tools, the accuracy achieved

is much higher than that of conventional field solvers. If the data on the surface

83

.7127632E-13 -.7115593E-12 .0000000E+00 .0000000E+00 .0000000E+00 100000
.4718448E-14 -.7105427E-13 .0000000E+00 .0000000E+00 .0000000E+00 010000
.0000000E+00 .0000000E+00 .7149836E-13 .7143869E-12 .0000000E+00 001000
.0000000E+00 .0000000E+00 .4718448E-14 .7127632E-13 .0000000E+00 000100
.0000000E+00 .0000000E+00 .0000000E+00 .0000000E+00 .7057580E-15 200000
.0000000E+00 .0000000E+00 .0000000E+00 .0000000E+00 .3560542E-13 110000
.0000000E+00 .0000000E+00 .0000000E+00 .0000000E+00 .2366163E-14 020000
.0000000E+00 .0000000E+00 .0000000E+00 .0000000E+00 .7064949E-15 002000
.0000000E+00 .0000000E+00 .0000000E+00 .0000000E+00 .3588926E-13 001100
.3567199E-13 .7027697E-15 .0000000E+00 .0000000E+00 .0000000E+00 100001
.4732326E-14 .3557007E-13 .0000000E+00 .0000000E+00 .0000000E+00 010001
.0000000E+00 .0000000E+00 .3581380E-13 .7045552E-15 .0000000E+00 001001
.0000000E+00 .0000000E+00 .0000000E+00 .0000000E+00 .2400857E-14 000200
.0000000E+00 .00000008+00 .4787837E-14 .3571015E-13 .0000000E+00 000101
.1250188E-12 .1237387E-11 .0000000E+00 .0000000E+00 .0000000E+00 300000
.2769872E-13 .3679666E-12 .0000000E+00 .0000000E+00 .0000000E+00 210000
.1039776E-12 .5195115E-13 .0000000E+00 .0000000E+00 .0000000E+00 120000
.9048318E-14 .3291182E-13 .0000000E+00 .0000000E+00 .0000000E+00 030000
.0000000E+00 .OOOOOOOE+00 .3758715E-12 .3729083E-11 .0000000E+00 201000
.0000000E+00 .0000000E+00 .6323985E-13 .7413065E-12 .00000008+00 111000
.0000000E+00 .0000000E+00 .3319892E-13 .3679442E-13 .0000000E+00 021000
.3321907E-12 .3294149E-11 .0000000E+00 .0000000E+00 .0000000E+00 102000
.2018379E-13 .3273745E-12 .0000000E+00 .0000000E+00 .0000000E+00 012000
.0000000E+00 .0000000E+00 .1109808E-12 .1101801E-11 .0000000E+00 003000
.0000000E+00 .0000000E+00 .2223924E-13 .3691161E-12 .0000000E+00 200100
.0000000E+00 .0000000E+00 .7609278E-13 .8903396E-13 .0000000E+00 110100
.0000000E+00 .0000000E+00 .1665335E-15 .3047020E-13 .0000000E+00 020100
.3292882E-13 .6518573E-12 .0000000E+00 .0000000E+00 .0000000E+00 101100
.6717630E-13 .1018835E-12 .0000000E+00 .0000000E+00 .0000000E+00 011100
.0000000E+00 .0000000E+00 .2599943E-13 .3276678E-12 .0000000E+00 002100
.0000000E+00 .0000000E+00 .0000000E+00 .0000000E+00 .1057490E-14 200001
.0000000E+00 .0000000E+00 .0000000E+00 .0000000E+00 .5331369E-13 110001
.0000000E+00 .0000000E+00 .0000000E+00 .0000000E+00 .4746203E-14 020001
.00000008+00 .0000000E+00 .0000000E+00 .0000000E+00 .1057172E-14 002001
.3793775E-13 .5576818E-13 .0000000E+00 .0000000E+00 .0000000E+00 100200
.2498002E-15 .4481073E-13 .0000000E+00 .0000000E+00 .00000008+00 010200
.0000000E+00 .0000000E+00 .1101168E-12 .4617591E-13 .0000000E+00 001200
.0000000E+00 .0000000E+00 .0000000E+00 .0000000E+00 .5388008E-13 001101
.2672016E-13 .5264869E-15 .0000000E+00 .0000000E+00 .0000000E+00 100002
.4732326E-14 .2664362E-13 .0000000E+00 .0000000E+00 .0000000E+00 010002
.0000000E+00 .0000000E+00 .2686198E-13 .5285215E-15 .0000000E+00 001002
.0000000E+00 .0000000E+00 .1000589E-13 .3795380E-13 .0000000E+00 000300
.0000000E+00 .0000000E+00 .0000000E+00 .0000000E+00 .4787837E-14 000201
.0000000E+00 .0000000E+00 .4815592E-14 .2678565E-13 .0000000E+00 000102

Table 3.5.1. The difference between the analytic and the computed Taylor maps. The
difference map is shown to the third order .

84

.0000000E+00 -.5380633E-15 .0000000E+00 .0000000E+00 .1222132E-14 400000
.0000000E+00 -.2578822E-15 .0000000E+00 .0000000E+00 .6219813E-13 310000
.0000000E+00 .0000000E+00 .0000000E+00 .0000000E+00 .1498060E-13 220000
.0000000E+00 .0000000E+00 .0000000E+00 .0000000E+00 .5182920E-13 130000
.0000000E+00 .0000000E+00 .0000000E+00 .0000000E+00 .4468648E-14 040000
.0000000E+00 .0000000E+00 .1829428E—15 .2166037E-14 .0000000E+00 301000
.0000000E+00 .0000000E+00 .0000000E+00 .8101576E-15 .0000000E+00 211000
.2389887E-15 .3119767E-14 .0000000E+00 .0000000E+00 .3914068E-15 202000
.0000000E+00 .7771467E-15 .0000000E+00 .0000000E+00 .1665009E-12 112000
.0000000E+00 .0000000E+00 .0000000E+00 .0000000E+00 .9770657E-14 022000
.0000000E+00 .0000000E+00 .1748517E-15 .2071154E-14 .0000000E+00 103000
.0000000E+00 .0000000E+00 .0000000E+00 .2469727E-15 .0000000E+00 013000
.0000000E+00 .5266477E-15 .0000000E+00 .00000005+00 .1087287E-14 004000
.2226095E-13 .4383921E-15 .0000000E+00 .0000000E+00 .0000000E+00 100003
.4767020E-14 .2219503E-13 .0000000E+00 .0000000E+00 .0000000E+00 010003
.0000000E+00 .0000000E+00 .2240254E-13 .4406926E-15 .0000000E+00 001003
.0000000E+00 .0000000E+00 .0000000E+00 .0000000E+00 .5058454E-14 000400
.0000000E+00 .0000000E+00 .2059464E-13 .5711100E-13 .0000000E+00 000301
.0000000E+00 .0000000E+00 .0000000E+00 .0000000E+00 .7202572E-14 000202
.0000000E+00 .0000000E+00 .4773959E-14 .2233597E-13 .0000000E+00 000103
.5283285E-10 .5267736E-09 .0000000E+00 .0000000E+00 .0000000E+00 500000
.1758265E-10 .2631307E-09 .0000000E+00 .0000000E+00 .0000000E+00 410000
.3701574E-11 .7016947E-10 .0000000E+00 .0000000E+00 .0000000E+00 320000
.4778458E-12 .1071568E-10 .0000000E+00 .0000000E+00 .0000000E+00 230000
.9821050E-13 .8942905E-12 .0000000E+00 .0000000E+00 .0000000E+00 140000
.1018630E-13 .5400845E-14 .0000000E+00 .0000000E+00 .0000000E+00 050000
.0000000E+00 .0000000E+00 .2646860E-09 .2644290E-08 .0000000E+00 401000
.0000000E+00 .0000000E+00 .7049192E-10 .1056648E-08 .0000000E+00 311000
.0000000E+00 .0000000E+00 .10763315-10 .2114512E-09 .0000000E+00 221000
.0000000E+00 .0000000E+00 .1961116E-13 .3858557E-15 .0000000E+00 001004
.0000000E+00 .0000000E+00 .1680600E-13 .1921857E-14 .0000000E+00 000500
.0000000E+00 .0000000E+00 .0000000E+00 .0000000E+00 .1550149E-13 000401.
.0000000E+00 .0000000E+00 .3178013E-13 .7154130E-13 .0000000E+00 000302
.0000000E+00 .0000000E+00 .0000000E+00 .0000000E+00 .9658940E-14 000203
.0000000E+00 .0000000E+00 .4787837E-14 .1955467E-13 .0000000E+00 000104

 

Table 3.5.2. The difference between the analytic and the computed Taylor maps. The
difference map shows some of the fourth and the fifth order terms.

85

enclosing the volume of interest can be given exactly, then in principle, arbitrarily
high accuracy limited only by the computational resources available can be achieved
by this new technique. In practical situations where the field data on the surface
enclosing the volume of interest are experimentally measured, the discretization of
the surface and the errors in the experimentally measured field data may limit the
accuracy achieved, but because the method is naturally smoothing, the accuracy is
expected to exceed that of the measurements. In addition, we have also developed
framework to compute the magnetic field of an arbitrary current distribution using
the Biot-Savart law. Using these techniques we can now extract the Taylor transfer
maps from the measured field data or the source data. Figure 3.6.1 shows the flow-
chart to extract the Taylor transfer maps from the measured magnetic fields data or

the source distribution.

86

CHAPTER 4

Applications

In this chapter we present examples to demonstrate the utility of the new field compu-
tation techniques we have developed. We first present the design of the new quadru-
pole magnet with elliptic cross section and tunable high order multipoles. Next, we
present the application of the Laplace solver to the realistic case of measure mag-
netic field data obtained for the MAGN EX spectrometer dipole magnet. Finally, we
present the ion-optic simulations for the Super-FRS and present the field computa-
tions for the quadrupole magnet utilizing the high performance parallel computing

environment.

4.1 The conceptual design of a quadrupole magnet

For charged particle beams that are wider in the dispersive plane than the trans-
verse plane it is cost efficient to utilize magnets that accept beams with elliptic cross
sections . We now present the conceptual design of a quadrupole magnet with an
elliptic cross section and with tunable high order multipoles. The design consists of

18 superconducting racetrack coils placed on two hollow concentric rhombic prism

87

support structures.

The analysis of this design will use the new numerical tools that we have described in
the section 3.4.4 to compute the local Taylor expansion of the magnetic field starting
from the Biot-Savart law and Ampere’s law. In this section we first discuss the case
where the magnet is considered to be infinitely long. This leads to purely transverse
magnetic fields (2D fields). In the section 4.1.1 we present results of numerical
computations for this 2D case. We also discuss the feasible range of multipole field
strength that can be achieved with this design. Finally, in the section 4.1.2 we present
the results for a finite length magnet and discuss the fringe field effects.

A combination of superconducting racetrack coils is used to produce the desired
magnetic field inside an elliptic cross section. By the proper choice of dimensions,
current density, and placement of these coils, various combinations of the quadrupole
field and the higher order multipole fields can be achieved. The support structure
that holds these coils in place consists of two concentric hollow rhombic prisms, with
the ratio of the diagonals of the rhombus equal to 2. The cross section view showing
the arrangement of the current coils on the support structure is shown in Figure 4.1.1.
The superconducting racetrack coils on the inner rhombic prism produce quadrupole
and octupole fields. The racetrack coils on the outer rhombic prism produce hexapole
and decapole fields, and also allow for a limited dipole field for correction purposes.
Figure 4.1.1 shows the layout of the various coils. The signs "+" and "-" indicate the
direction of the current to produce a positive multipole term.

Due to symmetry in the design about the central axis, it is sufficient to describe
only one quarter of the magnet. Figure 4.1.2 shows one quarter of the cross section.
The positions and the direction of the coils in this quarter are specified in Table 4.1.1.
All the coils have square cross section with thickness of 0.1 m.

In Table 4.1.1 the quantities Q11, Q12, H11, H12, H13, H13 are the magnitude of

the currents. For any given configuration of the inner current coils, the currents

88

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 4.1.1. The cross section view of the quadrupole.

     

Table 4.1.1. center p031 current quadrant

89

 

Dipole coil

 

 

 

Hex coil

 
  
 
 

Deca ooll

  
 

4..
Quad coll

 

 

 

Figure 4.1.2. The layout of the racetrack coils in the first quadrant.

(Q11, Q12) can be used as parameters to obtain different quadrupole and octu-
pole strength. Similarly, for any given configuration of the outer coils, the currents
(H11, H12) can be used as parameters to get the desired hexapole and decapole field
strength.

From the construction point of view it is cost efficient if the same type of coils can
be used. We use current coils of the same shape and size to generate the quadrupole
and hexapole fields. Also, we use the same type of current coils for octupole and

decapole fields.

4.1.1 2D design of the quadrupole magnet

For a 2D design the current coils can be considered to be of infinite length, thus
avoiding any fringe field effects. In this case a coil can be viewed as two current wires

of infinite length and finite cross section which are separated by certain distance. The

90

currents in these wire are equal in magnitude but are opposite in direction.

The magnetic field of the magnet described in the previous section is a superposition
of the fields produced by individual coils. In a 2D case the field of an individual coil is
again a superposition of fields generated by two infinitely long current wires of finite
cross sections . In the next section we show the results of the computation that has

been performed using the infinite wire tool.

Pure quadrupole configuration

A pure quadrupole configuration can be achieved by the following setting of the

currents.

QI1= 3965729.315790645 QI2= 80626.58244399808
HI1= 0.0000000000 HI2= 0.0000000000 HI3= 0.0000000000

The fifth order Taylor expansion of the magnetic field about the point (0.0, 0.0, 0.0)

that we obtained is shown in Table 4.1.2.

Pure hexapole configuration

A pure hexapole configuration can be achieved by the following setting of the currents.

QI1= 0.0000000000 QI2= 0.0000000000

HI1= 499792.9930948258 HI2= 938956.9619500297 HI3= 448723.7676402168

The fifth order Taylor expansion of the magnetic field about the point (0.0, 0.0, 0.0)

that we obtained is shown in Table 4.1.3.

91

MAGNETIC FIELD (Bx,By,BZ)

Bx By 32
0.44408920985018—15 0.000000000000 0.000000000000 000000
0.000000000000 46.09565826333 0.000000000000 100000
46.09565826333 0.000000000000 0.000000000000 010000
-0.7105427357601E-14 0.000000000000 0.000000000000 200000
0.1421085471520E-13-0.1421085471520E-13 0.000000000000 110000
0.000000000000 0.1136868377216E-12 0.000000000000 210000
0.000000000000 -0.2842170943040E-13 0.000000000000 120000
0.7815970093361E-12 0.000000000000 0.000000000000 030000
-0.142108547152OE-13 0.000000000000 0.000000000000 400000
0.000000000000 0.5684341886081E-13 0.000000000000 310000
-0.2273736754432E-12-0.17053025658245-12 0.000000000000 130000
-0.4263256414561E-13 0.5684341886081E-13 0.000000000000 040000
0.000000000000 -3956.535097021 0.000000000000 500000
-19782.67548511 0.000000000000 0.000000000000 410000
0.000000000000 39565.35097021 0.000000000000 320000
39565.35097021 0.000000000000 0.000000000000 230000
0.000000000000 -19782.67548511 0.000000000000 140000
-3956.535097021 0.000000000000 0.000000000000 050000

 

 

Table 4.1.2. The fifth order Taylor expansion of the magnetic field about the point
(0.0,0.0,0.0) for the current configuration producing a pure quadrupole field.

Operational plots

We have mentioned that the currents (Q11, Q12, H11, H12) can be used as parame-
ters to get the desired quadrupole and higher order multipole strength. However,
there is a maximum limit on the current density that the superconducting coils can
support. This puts a limit on the maximum quadrupole and other multipole field
strength that can be achieved. Because of the fact that each multipole is achieved
by superimposing the fields of several coils, this leads to Operating diagrams showing
achievable multipole settings.

To study this situation in detail, we now look at how the multipole strength depends

on the currents. The matrix given in the equation (4.1.1) relates the multipole field

92

MAGNETIC FIELD (Bx,By,BZ)

Bx By Bz

0.000000000000 -0.8604228440845E-15 0.000000000000 000000
-0.4718447854657E-15 0.000000000000 0.000000000000 100000

0.000000000000 -0.3885780586188E-15 0.000000000000 010000

0.000000000000 -18.24792017395 0.000000000000 200000
-36.49584034790 0.000000000000 0.000000000000 110000

0.000000000000 18.24792017395 0.000000000000 020000
-0.1776356839400E-14 0.000000000000 0.000000000000 210000
-0.1776356839400E-14 0.2664535259100E-14 0.000000000000 120000
0.2220446049250E-14-0.4440892098501E-15 0.000000000000 030000
-0.1776356839400E-14-0.3463895836830E-13 0.000000000000 400000
"0.25934809855243-12 0.000000000000 0.000000000000 310000

0.000000000000 0.4547473508865E-12 0.000000000000 220000
0.1705302565824E-12 0.1953992523340E-13 0.000000000000 130000
0.4440892098501E-14-0.4174438572591E-13 0.000000000000 040000
0.1776356839400E-14 0.8881784197001E-15 0.000000000000 500000
0.1421085471520E-13 0.1421085471520E-13 0.000000000000 410000
0.2842170943040E-13 0.3552713678801E-13 0.000000000000 320000
0.2131628207280E-13 0.000000000000 0.000000000000 230000
-0.5684341886081E-13-0.7105427357601E-14 0.000000000000 140000
-0.4440892098501E-14-0.3552713678801E-14 0.000000000000 050000

 

Table 4.1.3. The fifth order Taylor expansion of the magnetic field about the point
(0.0,0.0,0.0) for the current configuration producing a pure hexapole field.

strength at the horizontal half aperture to the currents in the coils for the specific case

of a horizontal half aperture of 0.5 m and a vertical half aperture of 0.25 m . In the

notation B(J1111)’ the superscript denotes the "y" component of the magnetic field and
the subscript (1111) gives the exponent in transport notation. Thus, By is the

(1111)

coefficient of 1:4 in the Taylor expansion of the "y" component of the magnetic field,
or the decapole term in the expansion. The equations (4.1.2) provide relationships
between the coefficients of other multipole terms in the Taylor expansion of the field

to the principle multipole coefficients Bl/ll)’ Bljlll) and By

(1111)'

93

 

 

 

 

 

 

 

3,9 0 0 —0.25137 —0.04316 +0.37029' 'Q11‘
51) +5.76974 +2.40063 0 0 0 Q12
3511) = 0 0 —3.89914 —2.08907 -—1.45431 . H11
3111) —0.40613 +15.44685 0 0 0 H12
Bi 0 0 +1.66569 —2.32478 +2.99743 H13
_ (1111)_ - - ' ‘
(4.1.1)
y __ '11
13(22)— 801) (4.1.2)
y __ y
3(122)“ 35(111)
y
3(1122) _ By _ By
" 5 — (2222)_ (1111)
a: __ y
3(2)_B(1)
x __ y
302) ‘ 28(11)
1'
B(112) : —B‘” 2 By
3 (222) (111)
a: _ a: _ y
B(1112) — ‘B(1222) — 4B(1111)

Operational plot for the quadrupole and the octupole fields

Horn the equation (4.1.1) it can be seen that the quadrupole field strength and the
octupole field strength are coupled via the currents (Q11, Q12). We vary both of these
current densities in the range [—108, 108] A / m2 and plot the resulting octupole field
strength and the quadrupole field strength; the results are shown in Figure 4.1.3.
This plot gives the possible values of the quadrupole and octupole strength that can

be achieved with the configuration of the coils described in the section 4.1.

94

Quadrupole VS Octupole

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

20 I I I I I I I I I
15» f' H‘ "- --- v- ~
101.. .... -.. ¢ ........- .A. A ..... fifi _,

A “ v v

6') - W“

x 5" ‘ = v - A 1 _

V h

o-t ‘—kA A A- _

é ' - ‘ r A 1.4

O 0'- ‘:A A 3 “ -;A - ‘ -‘

O 9:- v A -

31 ' " ' r = ~—~A

3 A A v'

8 _5,_ 4- - 1 - 1 - a
'10_ H r 5'"; 7 Av 7.:' fA' A d
-151— "v—‘Aé A AAA‘A' ' L_- A 1

1 1 1 1 1 '
-20 l 1 1 l

 

 

-10 -8 -6 -4 -2 O 2 4
Quadrupole Cofticlent (11)

Figure 4.1.3. The operational plot for the quadrupole and the octupole. The coeffi-
cients are computed at the horizontal half aperture.

95

Hexapole VS Decapole

 

Decapole Coefficient (x 4)
o
I

 

 

 

 

Hexapole Coefficient (x 2)

Figure 4.1.4. The operational plot for the hexapole and the decapole. The coefficients
are computed at the horizontal half aperture.

Operational plot for the hexapole and decapole fields

Horn the equation (4.1.1) it can also be seen that the dipole, hexapole and decapole
field strength are coupled via the currents (H I 1, HI 2, H 13). However, under normal
operation, we have a strict requirement of zero dipole field for this magnet. The dipole
field is set zero by the proper choice of the current H13. Once again we vary the
current densities of all currents in the range [—108, 108] A / m2 and plot the decapole
field strength versus the hexapole field strength; the results are shown in Figure 4.1.4.
This plot gives the possible values of the hexapole and decapole strength that can be

achieved with the configuration of the coils described in the section 4.1.

96

Figure 4.1.5. The schematic digram of a current coil.

Optimization of the operational region

The matrix given in the equation (4.1.1) is influenced by the geometrical design para-
meters of the system. In order to optimize the operational region of the currents and
the fields, we need to find the optimal geometric configuration of the coils described
in the section 4.1, where the optimal design is defined as the one that would decouple
the influence of the octupole coil current on the quadrupole component of the field
and vice versa. And, at the same time maximize the coupling strength of the cur-
rent in quadrupole coils on the quadrupole component of the field, and maximize the
coupling strength of the current in octupole coils on the octupole component of the
field. The same type of optimization is also required for the hexapole and decapole

components of the field.

4.1.2 3D design of the quadrupole and the fringe field analy-

815

Figure 4.1.5 shows a current coil of a rectangular cross section . This coil is made by

joining four wires having rectangular cross section with the current entrance and exit

97

planes tilted by 45° in opposite directions. This ensures that there is a continuos flow
of current through the coil. The magnetic field for each of these wires is calculated
using the COSY tools described in the section 3.4.4.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 4.1.6. The three dimensional layout of quadrupole coils.

The 3D layout of the quadrupole with four current coils is shown in Figure 4.1.6.
We note that when the length of the magnet is large compared to the aperture of the
magnet, the magnetic field at the center of the magnet is identical to the magnetic
field obtained by the 2D design in the section 4.1.1. We also verify that the curl and
the divergence vanishes at all points inside the magnet.

To present the results we once again consider the design presented in the section
4.1.1, but with a modification that the length of the magnet is finite. We consider a
magnet of length 1 m, extending from —0.5 m 3 2 S 0.5 m. We now look at the mag-

netic field generated by this coil configuration at four different planes perpendicular

98

to the central axis, 2 = 0, z = 0.25 m, z = 0.5 m and z = 1.0 m.

From figures 4.1.7 and 4.1.8 we can see that at the plane 2 = 0 the magnetic field
agrees with the textbook model where 8;; = Cy, By = Cr. In the z direction the
magnitude of the magnetic field is of the order ~ 10—16, which is zero for all practical
purposes.

Now, as we start going away from the center, we start seeing deviation from the
ideal behavior (2 = 0). Figures 4.1.9 and 4.1.10 show the plots for the magnetic field
at the quarter length of the magnet, z = 0.25 m. There is no significant deviation
from ideal behavior in the :1: and y components of the magnetic field. In the z direction
we notice that the magnetic field is nonzero. However, the magnitude is still very
small compared to the components Bx and By.

Figures 4.1.11 and 4.1.12 show the plots for the magnetic field at the entrance and
the exit of the magnet, for z = 0.5 m. Now we see that the magnetic field is different
compared to the field at the center of the magnet. We see that the strength of 8;;
and By falls by a factor of five. The 82 component is very strong, which is almost
three times as large as Bx or By.

Figures 4.1.13 and 4.1.14 show the magnetic field on a plane 0.5 m away from
the entrance and exit of the magnet. The overall field falls off significantly and its
magnitude is ~ 10—1 tesla. Finally, the plot 4.1.15 shows the magnetic field on the
first quadrant of the magnet on the y = 0 m plane. The region stretches from the
center of the magnet to 0.5m (half length of the magnet) outside the magnet in both

:1: and 2 directions. Here the fringe field that falls off in the region can be clearly seen.

99

Bx for Quadrupole at -0

5
2
O

 

 

(

By for Quadrupole at -0

. w _ _ .
N W _ a n
u 1 a 1
n w _ _ u
u m _ _ u
u m . . "
86420

—-_..-____.

5
1..
0

 

 

Figure 4.1.7. (a) The plot for Bx for the data on the plane 2: = 0 m. (b) The plot for

By for the data on the plane 2 = 0 m.

100

82 for Quadrupole at —0

82
29-016

 

 

(a)

Bx‘i+By*j vector at 2:0

 

 

 

 

 

 

 

 

 

 

 

 

 

 

0.15

liq lllllllllll
flflffillll llll l l l l l l llllllll
"'l"""' ffffff l hhhhhhhh

l lllllllll
.71 lllllllllllllllll l lllll
l lllllllll
fl'lflflfllll l l l l l l l llll‘t“\\\
flo/I/el/ollll l l l l ll‘l“t\\\\
//o/,/elel”clrl lllllllll t\\1\\\\\\\\\\\
////e//”’ llllllllllll ‘t\§\|\\\\\\\\
/////////I llllllllll \\\\\\\\\\
III/III]! llllllllllll \\\\\\\\\
///////III llllllllll \\\\\\\\\\
III/III! rrrrrrrrrrrrrr \\\\\\\\
’IIIII’I’I ooooooooooo \\\\\\\\\
’p””.; ............. uss~a-sx
....... ............... ..4.‘._a
I ..a...... ............. ....... l
~.~....... ............. .....v.
s‘sssssss ............. yoorrarr
\\\\\\\\\\\ oooooooooo viii/tilt
\\\\\\\\\\\ 11111111111 [III/III
\\\\\\\\\ lllllllllll Ill/lll/l/
\\\\\\\\\ llllllllllll l////////
\\\\\\t\\\ lllllllllll II////////
\\\\\\\\\\1\1 lllllllllll I’ll/ll/I
\\\\\\|\§\I\O llllllllllll [I’ll/ll/l
\\\\\\\\‘\\|\\l l l lilJl/l/l/l
\u\1\\\t\t\t\1 llllllllll 11 1111/1/9/
\t\1\t\tlil| ill i i l i l l f {fill/Old
\\ lllllllllllllllllllll
1111111111
\Q\Ilhlhl llllllllll l l l i I l 14 flillll'
P
0
w
.
nw

0.25

(b)

Figure 4.1.8. (a) The plot for 82 for the data on the plane 2 = 0 m. (b) The vector

plot of the two dimensional field on the plane z = 0 m.

101

Bx for Quadrupole at 2:25

 

0.25

 

 

(a)

By for Quadrupole at Z=.25

 

 

-0.15
0.25

(b)

Figure 4.1.9. (a) The plot for Ba; for the data on the plane 2 = 0.25 m. (b) The plot
for By for the data on the plane 2 = 0.25 m.

102

82 for Quadrupole at Z=.25

.

m _
m.
_ .
. .
m .

 

 

(a)

Bx*i+By*j vector at Z=.25

 

 

 

 

 

 

 

 

 

 

 

 

 

41 if lllllll‘
1?"””99“11 illlvfflffifil
111111111119111 .1 lllllll

lllllllllll
ff fffffffffffff 1 ‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘‘
[011 ffffffffffff l lllll‘l\\l\\\\\\
f/ll/Olllfflflf 1 1 ll lll ‘t\\\a\l\s\\
l/f/f/III’ llllllllll ‘\\\\\I\\\\\
///////’11 llllllllll \\\\\\\\\\
////////1/ llllllllll \\\\\\\\\\
l/l/l/I/I,’ llllllllll \\\\\\\\\
I/l/l/l/III llllllllll \\\\\\\\\
I/I/II/III rrrrrrrrrrr \\\\\\\\\
xxx/’1’), 1111111111111 \‘\\\\\\
»””,;a. ............. .\\‘\s\s
L....... ............... .....a.

I .......... ........... ........ l
--....u .............. ....vv.
\ss\\xse\-s ........... cacrrr”
xxxxxxxxsss 1111111111 y’all/1”
\\\\\\\\\\\ llllllllllll III/Ill
\\\\\\\\\\\ 1111111111 ll/I/I/l/
\\\\\\\\\\ lllllllllll (Ill/Ill/
\\\\\\\\t\\\ lllllllllll /////////
\\\\\\\\\\ llllllllll //////////
\1\I\I\\\\|\Q\t\t\|\t llllllllll I’ll/l/l/
\Y‘Y‘N‘N‘l‘lllllllllf’ll’I/ll/
\\\\\t\\\\\t\tttll lllllll 1 1111111111/
\X\ 111111111111111111111111111111
llllllllllllllllllllllll
111111111
\u\\l1|lhltl 11111111111111
l 1111111111
b
0
.6
O.

-0.15

0.25

-0.25

(b)

(b) The

0.25 m.

Figure 4.1.10. (a) The plot for B; for the data on the plane 2

vector plot of the two dimensional field on the plane 2 = 0.25 m.

103

Bx for Quadrupole at 2:05

 

 

Bx ——
1.5 -----------
0.5
o _ .........
-1
-1.5
-2 _—
0.25
(a)
By for Quadrupole at 2:05
By ——
2 -. ....... .-
0 ...........
_1 _ .........

 

 

(b)

Figure 4.1.11. (a) The plot for Bay; for the data on the plane 2 = 0.5 m. (b) The plot
for By for the data on the plane z = 0.5 m.

104

82 for Quadrupole at Z=0.5

 

 

(a)

Bx’i+By*j vector at Z=0.5

 

 

Ill —
Ilauh!|""'ll'lll|||-|"--oeno‘
Insaetl'III'I'IFII'IIIIIOOI
acOIOIOIIIIIIIIIIIIIIIIIIIIe

ad...--""'III"-‘II“““"0

0‘IF'FIFIIII'II'IIII‘IIIQQ

OJIIF"""""I'I‘I““ \\\\\

lif""""""-‘|‘ ““““““““

\\\\\\\\\\\\\\\\\\\\\\\\\\

\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

t l l \\\\\\\\\\\\\\\\\\\\\\\\\\\
0 I t I I llllllllllllllllllllllll
o o t ooooooooooooooooooooooooo s Q
a o ooooooooooooooooooooooooo q g

y o o o ........................ s
v . . ........................ . .
c o c ...........................
Q Q . ooooooooooooooooooooooooooo
\ ~ s s ooooooooooooooooooooooo a

\ ‘ - \\\\\\\\\\\\\\\\\\\\\\\\ a

llllllllllllllllllllllllllllll
IIIIIIIIIIIIIIIIIIIIIIIIIIIII
lllllllllllllllllllllllllllll
llllllllllllllllllllllllll

toIII'-‘|l‘“-‘|-'l“““‘lil

uvIl0||“h|l|n.-'llu.l.-‘l“1410l
InnIl---‘ll|-|0lll0101|l.n.ll“‘0‘00
\Ioea-l-||1||¢l.lllllllu“"“-ll

\‘ouellfi’hfi‘fillllllllllllllllhllIIIe.cl

P

 

 

0.15

0

-0.15

0.25

-0.25

(b)

(b) The

Figure 4.1.12. (a) The plot for B; for the data on the plane z = 0.5 m.

vector plot of the two dimensional field on the plane 2 = 0.5 m.

105

Bx for Quadrupole at 2:1

 

 

 

Bx _—
0.02 ---.-._,_.-

'0.02

0.25
(8)
BY for Quadrupole at 2:1

By _—

O.1
0 -._._ --.-

-0.05 ------ -
-O.1

 

 

(b)

Figure 4.1.13. (a) The plot for Ba; for the data on the plane 2 = 1 m. b The plot
for By for the data on the plane 2 = 1 m.

106

82 for Quadrupole at —1

 

 

(a)

Bx‘i+By*j vector at 2:1

 

 

llllllllllllllllllllllllllll
\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\\\\\\\\\\\\\\\\\\\\\\\\\\\
llllllllllllllllllllllllllllll
IIIIIIIIIIIIIIIIIIIIIIIIIIIIII
llllllllllllllllllllllllllllll
llllllllllllllllllllllllllllll
IIIIIIIIIIIIIIIIIIIIIIIIIIIIII
eeeeeeeeeeeeeeeeeeeeeeeeeeeeee
oooooooooooooooooooooooooooooo
oooooooooooooooooooooooooooooo
..............................
..............................
..............................
.............................
..............................
oooooooooooooooooooooooooooooo
oooooooooooooooooooooooooooo
eeeeeeeeeeeeeeeeeeeeeeeeeeeeee
eeeeeeeeeeeeeeeeeeeeeeeeeeeeee
eeeeeeeeeeeeeeeeeeeeeeeeeeeeee
IIIIIIIIIIIIIIIIIIIIIIIIIIIIII
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
IIIIIIIIIIIIIIIIIIIIIIIIIIIIII
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
111111111111111111111111111111
IIIIIIIIIIIIIIIIIIIIIIIIIIIIII

\“"1‘“-‘I‘Il'l’falfl'n7'1'l”,

 

 

0.15

0.25

(b)

Figure 4.1.14. (a) The plot for B; for the data on the plane 2' = 1 m. (b) The vector

plot of the two dimensional field on the plane 2

=1m.

107

By for Quadrupole at Y=0 plane

  
   

..
‘ O.
Q. ..
......O.
----
v #04:... 0

blbomaaa

 

 

- .. -.. ...“:
\V....~.~.~.~.~.~.~.~

       
  

Figure 4.1.15. The plot of the magnetic field on the midplane, y = 0 m. Only the

magnetic field in the first quadrant is shown.

 

Figure 4.2.1. The dipole magnet of the MAGNEX spectrometer.

108

 

Figure 4.2.2. The layout of the measurement grids in different regions of the dipole
magnet.

4.2 The dipole magnet of the Catania MAGNEX
spectrometer

We now address a practical application of the Laplace solver technique to magnetic
spectrometers. The trajectory reconstruction method [21] is one of the important
tools to study magnetic spectrometers. Good computational modelling of the dipole
magnet is very important for this tool to work, and this is particularly so for modern
large-aperture devices such as MAGNEX at INFN, Catania, Italy [70, 49, 29]. Figure
4.2.1 shows the MAGNEX spectrometer configuration.

For purposes of measurement economy, magnet builders usually provide the mag-
netic field only on few separate horizontal planes within the dipole, while the com—
putational treatment of the device requires the knowledge of the field in all of space.
The MAGNEX dipole was divided into a number of volumes defined by areas and
planes as shown in Figure 4.2.2. Four areas were mapped as indicated in Table 4.2.1;
Areas 1 and 4 comprise the Effective Field Boundary regions of the magnet at the

entrance and at the exit where the field undergoes a sudden variation due to the

109

fringe field effects, whereas region 2 and 3 represent the central region of the magnet.
This subdivision is the result of the need of different grid sizes over the mapped area
in order to limit the measurement time. For each of the regions the measurements

were taken on seven different planes as shown in Table 4.2.1.

cm

area en
area en

area ex1
area ex1

   

. m1
Table 4.2.1. (a)Areas mapped in the dipole b Planes mapped in the dipole.

The magnetic measurement were organized so that the RMS error (AB; / B) i =
:r, y, 2 at any mesh point inside the working volume of the magnet was not greater
than 5 x 10-4. The field measurement error due to the error of measuring the Hall
probe voltage was AB = 21:5 x 10‘5 Gauss. The main source of the B measurement
error were assumed to be the errors of positioning the Hall probe [71, 65, 41].

Utilizing that sufficiently outside the dipole the fields will vanish, it is thus possible
to provide field data over the surface of a finite box enclosing the region of interest, and
thus to apply the Laplace solver technique to obtain a field representation everywhere.
We use this method to compute the fields in the region 1 and plane A of the dipole
magnet. The contour plot of the resulting relative errors is plotted over the region
1 in Figure 4.2.3. The region where the sharp valley is observed coincides with the
physical boundary of the dipole magnet.

Figure 4.2.3 we can see that the relative error is not as small as we would expect
from our approach. Since the data was measured in several planes and in several
overlapping regions on each plane it is very hard to extract the magnetic field data
on a closed surface. To demonstrate our approach we have used only the data for

the entrance of the dipole magnet (region 1) and extracted data on a surface of a

110

Log(Error)

 

X axis

Figure 4.2.3. The contour plot of magnetic field errors for the region 1 and plane A.

111

Low-Energy

S U per-F RS Branch

  
  
    
 
 
 

:- K Energy
Buncher

  

Main-Separator

Pro-Separator High-Energy Branch

- = "ea—W
,— wfi~$ Rlng Branch
; «_Production
Target

Focusing :
3’3“” 50 m

 

Driver
Accelerator

Figure 4.3.1. The layout of the Super-FRS. Spatially separated rare-isotope beams
are delivered to the experimental areas via three different branches.

rectangular box enclosing the region of interest. Since the data was measured on only
seven different planes the data on two surfaces, :1: = i0.5 m, of the rectangular box
enclosing the region of interest are discretized with very large step size. This may
lead to large errors from the surface elements on these planes (3: = $0.5 m). Also,
compared to the case where the whole data is enclosed by one closed surface the
smoothing effect is limited for the rectangular box enclosing only the data on region
1 for seven different planes. The difference we observe may also be attributed to the

errors in the experimental measurements and may need independent verification.

4.3 Ion-optic simulations for the Super-FRS

The superconducting fragment separator (Super-FRS) [38] being built at GSI, Ger-
many will be the most powerful in—flight separator for exotic nuclei up to relativistic
energies. The Super-FRS is a large-acceptance superconducting fragment separator

with three branches serving different experimental areas. The Super—F RS uses the

112

 

Parameters
Max. Magnetic RigidityTp ax 20Tm
Momentum Acceptance AP P :l:2.5%
Angular Acceptance ¢y¢x :t40mrad2
Wmentum Resolution 1500
Table 4.3.1. The Super-FRS design parameters.

 

 

 

 

 

 

B p — Ap — B p method, where a two-fold magnetic rigidity analysis is applied in front
of and behind a specially shaped energy degrader. As a consequence the Super-FRS
consists of a two-stage magnetic system, the pre-separator and the main-separator,
each equipped with a degrader. Figure 4.3.1 shows the layout of the proposed Super-
FRS. Table 4.3.1 summarizes the design parameters for the Super—F RS.

The ion-optic design for the Super-FRS has been primarily studied using the beam
physics code GICOSY [19]. The code GICOSY computes the transfer maps to the
fifth order. However, for some electric and magnetic elements the GICOSY code
computes transfer maps to only the third order, and hence the maps generated are
only correct to the third order. The first, second and third order design studies have
been conducted using the code GICOSY. However, for high acceptance and high
resolution machines like the Super-FRS, the high order transfer map computations
are necessary to study the effect of the high order nonlinearities on the resolution
achieved. The code COSY INFINITY [23, 24] is an arbitrary order beam physics
code and hence can be utilized to perform the high order design studies of the Super-
FRS.

The theory and the implementation details of the GICOSY and the COSY INFIN-
ITY codes can be found in [77, 78, 26, 79, 20, 15, 14, 16, 25, 52, 13, 7, 8, 12] and
references therein. Both GICOSY and COSY INFINITY codes are written in the
Fortran 77 language, and they themselves provide a programming environment to
describe a particle optical system conveniently and in an intuitive manner. However,
due to the differences in their approach to the description of the optics, and their

implementation in the Fortran 77 language, they may provide different results for

113

some of the ion or particle optical elements. Though they may be very small, these
differences may be of relevance, for example, when analyzing instruments of high
precision, like fragment separators, where the higher order aberration corrections are
important to get the desired resolution. A systematic study of the differences in the
description of the important particle optical elements in the GICOSY and the COSY
INFINITY codes has been done and the results are presented in the Appendix A.
Also, to facilitate code reusability, a code converter to translate a GICOSY program
into the COSY INFINITY program has been developed and is discussed in Appendix
B. Using this program, all GICOSY files were converted to COSY INFINITY files
to perform the high order simulations. It has been verified that with same choice of
units and the scaling factor and ignoring the fringe fields, the results from the COSY
INFINITY and GICOSY codes agreed to good accuracy to third order.

The code COSY INFINITY has been optimized to do single particle simulations
for large multiturn accelerators like the LHC and the FNAL. There is a plethora of
tools available in the code COSY INFINITY to study and present results for these
multiturn accelerators. However, to present the results for an ensemble of particles
for a single pass accelerator like the Super-FRS, we need to implement new output
tools. New tools developed for the purpose are discussed in Appendix C.

We now present the results of the simulation for the pre-separator stage, the main-
separator stage and the energy buncher using the code COSY INFINITY. We then
apply the Laplace solver technique to the magnetic field data obtained from the
TOSCA simulation of the Super-FRS quadrupole and extract the transfer map for

this magnet.

4.3.1 The pre—separator branch

The standard ion-optical layout of the pre—separator stage is presented in Figure

4.3.2a. The most relevant first order matrix elements are given in Table 4.3.2 at four

114

     

Table 4.3.2. The evant matrlx or t pre-separator stage.

focal planes PFl, PF 2, PF 3, and PF 4. In the first focal plane a focus is realized only
in the x direction, whereas in the symmetrical mid plane, the position of the first
degrader system, foci are required in both the x and y coordinates. Also a parallel
dispersion line is required at PF2. The system is mirror symmetric with respect to
PF2, firstly to have the necessary three foci to achieve the achromatic condition at
PF 4, secondly to minimize the geometrical image aberrations. Figure 4.3.3 shows the
envelopes and the dispersion line for the primary-beam emittances of 4071 mm mrad
and Ap/p of 2.5 % . The target spot size is assumed to be :l:1 mm and i2 mm in
the x and y directions, respectively.

The benefit of the design is that the higher order aberrations are small except for
the chromatic contributions, especially (:13, ad) which is responsible for an enormous
focal plane tilt. The tilt angle at PF2 would be about 7 degrees without second-order
corrections by means of hexapole magnets, while the required tilt is 90 degrees. Such a
tilt would make a high resolution degrader operation very difficult or even impossible.
Therefore, hexapole magnets are used to correct this deficiency. However, as soon as
this correction is done, induced higher order aberrations become a major challenge.
Applying hexapole and octupole corrections, about 80 ‘70 of the first-order momentum
resolving power at the central plane of the pre-separator was regained.

Usually, the resolution is computed using formulas which only consider the impor-

tant aberration terms. However, we can also compute the true nonlinear resolution

115

7]

 

\
/

 

 

 

 

 

 

 

 

 

 

 

1"1—1 —1 W'-

 

 

 

 

 

 

 

 

 

 

 

 

hirl/1 L m

[K \ a _.l

y W J..OJ

7 r

a /. L P

t\l N h e

.m.

In, 1% x

8

km) in

J/t 1w me
—\\ u

 

 

(a) Layout of the Pre-Separator

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

1 _m

— f ..l

/ M

1 u ....OJ

\ r.

_ \ l— p

n L e

.m

p

/ h

_ J _T

\ )

_ _10\
u .

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

NI
1]./~—

Figure 4.3.2. The figure shows the pre—separator layout and the projections of the
116

trajectories into the x and y planes.

 

0.3 I I I l l
X-Envelop —-—
Y-Envelop ---------

o 2 r DiSpersion Plot

0.1

 

 

 

 

o
-o.1 ~ «
-o.2 » 1
-o.3 . 4 . . 1 1
0 1o 20 30 40 so so

Figure 4.3.3. Envelopes and the dispersion line for the pre-separator stage.

   

Table 4.3. . on- ution at PF2.

by launching N particles, uniformly distributed in the phase space, and compute
the spot size of the resulting beam distribution. The nonlinear resolution computed
using this is a good measure of the overall nonlinear effects. Table 4.3.3 shows the
nonlinear resolution achieved at the focus PF2 by changing the order of computation.
The strength of sextupoles and octupoles were fitted to maximize the resolution while
satisfying other necessary conditions. We see that the resolution does not drop sig-
nificantly with the increase in the order of computation. In the pre-separator stage

the required resolution can thus be achieved with the current lattice settings.

117

   

Table 4.3. . on- utlon at MF2.

4.3.2 The main-separator

The main-separator consists of 4 dipole stages with focusing elements in front of and
behind each dipole magnet system. The main-separator also has 4 focal planes, MF1,
MF2, MF 3 and MF4 to accommodate an achromatic system with a degrader station
in the central focal plane MF2.

An analysis similar to the pre-separator stage is performed for the main-separator
stage. Table 4.3.4 shows the non-linear resolution achieved at the focus MF2 by chang-
ing the order of computation. Once again the strength of sextupoles and octupoles
were fitted to maximize the resolution while satisfying other necessary conditions. For
the main-separator stage the resolution changes significantly with change of order of
computation. The main-separator stage may require further analysis to correct for

high order aberrations.

4.3.3 The energy buncher

After production and in-flight separation in the pre-separator and the main separator
stages of the Super—FRS, the secondary beam can be distributed to either of three
experimental areas: a high-energy cave, a storage ring complex and a dedicated low-
energy branch (LEB), as shown in Figure 4.3.1. In the LEB the exotic species can be
slowed-down to ion energies ranging from a few hundred MeV/ u down to rest. The
LEB is equipped with an energy-bunching stage in order to reduce the energy spread

of the secondary ions originating from the fragment production process and from

118

 

Stopplng cell

1......

ES 83

 

Figure 4.3.4. Ion-optical layout of the energy buncher.

119

straggling in the target and the degraders. For the energy buncher the simulations
were performed using the code COSY INFINITY.

The energy buncher consists of a dispersive ion-optical stage with a large split di-
pole magnet system, quadrupole and hexapole magnets. The magnetic quadrupole
triplet in front of the dipole magnet is needed to properly illuminate the field vol-
ume of the dipole magnet to reach the required resolving power and to focus the
secondary beam onto a monoenergetic degrader. The layout of the system and the
trajectories are shown in Figure 4.3.4. The quadrupole triplet behind the monoener-
getic degrader guides the exotic nuclei into the gas cell or any other detector array. Its
main characteristics are: Bpmax = 10 Tm, 63; = 300 mm mrad, ey = 200 mm mrad,
Ap/ p = 21:2.5%, (transverse and longitudinal acceptance). Under these conditions a

momentum resolving power of R = 600 can be achieved.

4.3.4 The Super-FRS quadrupole magnet

Most of the current ion-optic codes use the analytic models for the magnets and
then apply an assumed fringe field models to obtain the transfer map. But, with
the Laplace solver techniques we have developed it is possible to extract the transfer
maps directly from the magnetic field information. We now use the example of the
Super-FRS superferric quadrupole to extract the magnetic field and the transfer map.
The Super-FRS superferric quadrupole is designed with effective length of 0.8 m with
the usable horizontal aperture of $0.2 m and the vertical aperture of :l:0.1 m. Using
the TOSCA package, the magnetic field data was obtained on a rectangular surface
grid enclosing the region of interest. The surface was discretized with a step size of
5 mm, leading to a discretization of 80x40x320 surface elements. The relative field
computation error was AB / B = 70 x 10-4. Figure 4.3.5 shows the T OSCA model
for the magnet. We now use this surface date to compute the magnetic field inside

the volume and compare the result with the TOSCA simulation for the points inside

120

 

Figure 4.3.5. The TOSCA model for the quadrupole magnet.

the volume.

Figure 4.3.6 shows the difference between the TOSCA and the computed magnetic
field for the relative error of the y component of the magnetic field on the midplane of
the quadrupole magnet. Due to double midplane symmetry it sufficient to look at the
results only for the one quarter of the magnet. It can be seen that the results agree
~ 10‘3, which is within the TOSCA model computation error of AB / B = 70 x 10‘4.

We now study the error dependence on the size (length) of the volume element,
or equivalently the number of volume elements chosen for the computation. For the
ninth order computation Figure 4.3.7 displays the dependence of the rms average
error on the length of the volume element. As an example, for cell lengths of 0.1,
with 410 volume data points, the rms error of 8 x 10“5 is observed for the volume
element at (0,0, 0) with a ninth order computation.

We now extract the transfer map from the Super-FRS quadrupole magnetic field

121

Log(Difference in relative error)

 

Computational error in B_Y component
-4.

 
   
  
  
  
  
   

-2.5 - -3.5 ---------
-3 - ........

'3 ‘ -2.8 --
-2.7

-3.5 -1
-4 .
-4.5 -
-5 ..
-5.5 -

 

 

0.1
X axis (m)

0.2

Figure 4.3.6. The difference between the relative error of the y component of the
magnetic field on the mid plane. The figure only shows the results for the first
quadrant.

 

11". '

122

 

1 l r l I l 1
Sr b po'nt (0.0,0.0,0.0) _1__
6r b poht (-0.1,-0.025,-0.2) ----x .....
0 — Sr b poht (0.2,-005,04) ..... 41-.-. 1
Sr b poht (0.3,-0.05.06) a
,El
-1 — x x _

 

 

 

 

.2 _
1% -3 . -
rm
.4.“
5 4 ~ 1
3’

-5 — _

'6 L 1 1 l 1 1

-8 -7 -6 -5 4 -3 .2 -1

loablta'mp '39)

Figure 4.3.7. The rms average difference between the TOSCA simulation result and
the new Laplace solver technique versus the volume element length for four volume
points (0.0,0.0,0.0), (-0.1,-0.025,-0.2), (-0.2,-0.05,-0.4) and (-0.3,-0.075,—0.6).

123

data. We solve the ODEs of motion as described in the section 2.3.2. We use the
multipole expansion solution of magnetic field computed through the Laplace solver
technique to compute the transfer map for the quadrupole. Tables 4.3.5, 4.3.6 and

4.3.7 present the transfer map computed using the magnetic field data on a closed

surface.
-0.4705674 -1.394826 0.000000 0.000000 0.000000 100000
0.5581815 -0.4705674 0.000000 0.000000 0.000000 010000
0.000000 0.000000 3.837901 4.272580 0.000000 001000
0.000000 0.000000 3.213394 3.837901 0.000000 000100
0.000000 0.000000 0.000000 0.000000 1.000000 000010
0.000000 0.000000 0.000000 0.000000 0.3989286 000001
0.1284348E-14 0.1535115E-14 0.000000 0.000000 -0.4476261 200000
0.1159401E-14 0.9402369E-15 0.000000 0.000000 0.4865291E-01 110000
-0.1197808E-14-0.3977569E-14 0.000000 0.000000 -0.1627172 020000
0.1930759E-13 0.5931886E-13 0.000000 0.000000 ~2.059670 002000
0.3353931E-13 0.9565057E-13 0.000000 0.000000 -3.933253 001100
0.4768188 -0.4891389 0.000000 0.000000 0.000000 100001
0.1259816 0.4768188 0.000000 0.000000 0.000000 010001
0.000000 0.000000 -1.858375 -0.9955222 0.000000 001001
0.1398535E-13 0.3825810E-13 0.000000 0.000000 -1.984025 000200
0.000000 0.000000 -2.589889 -1.858375 0.000000 000101
0.000000 0.000000 0.000000 0.000000 -0.2995974 000002

Table 4.3.5. The extracted Taylor transfer map for the Super-FRS quadrupole. The
Taylor transfer map is shown to second order.

The magnetic field data does not extend to the region far enough for the fringe

 

field to vanish. Hence the transfer map computed may not accurately represent the

real transfer map. Further analysis is required to obtain the true transfer map.

 

4.4 Summary I!

In this chapter we have successfully demonstrated the utility of the magnetic field
computation techniques we have developed. These techniques were applied to ana-
lyze a new design of a quadrupole magnet, and to find the multipole decomposition of

the magnetic field of the MAGNEX spectrometer dipole magnet, and the quadrupole

124

-5.023497
-8.056533
-6.019470
-2.428532
0.000000
0.000000
0.000000
26.87556
41.20143
0.000000
0.000000
0.000000
0.000000
29.08193
56.42227
0.000000

-1.219005
-4.472267
-4.481049
-3.420222
0.000000
0.000000
0.000000
33.22080
55.84502
0.000000
0.000000
0.000000
0.000000
48.08264
81.12632
0.000000

-0.1159274E-15 0.000000

-0.2781943E-14-0.4909952E-14
0.6655356E-15 0.7763230E-15
-0.2643101E-13-0.4295697E-13

7.497146 18.20560
19.69457 29.28498
0.000000 0.000000
-0.5441379E-13-0.1018488E-12
-0.2454539 0.3186799
-0.2452654 -0.2454539
0.000000 0.000000
0.000000 0.000000
*0.2705811E-13-0.5399731E-13
0.000000 0.000000
0.000000 0.000000

0.000000
0.000000
0.000000
0.000000
39.46047
44.92503
24.62299
0.000000
0.000000
-39.97053
21.56656
25.13783
17.11116
0.000000
0.000000
-68.87387
0.000000
0.000000
0.000000
0.000000
0.000000
0.000000
-37.98898
0.000000
0.000000
0.000000
1.639435
-5.482791
0.000000
2.527908
0.000000

0.000000
0.000000
0.000000
0.000000
50.88769
63.97251
41.94641
0.000000
0.000000
~97.22287
27.15591
36.62275
29.54703
0.000000
0.000000
-200.9709
0.000000
0.000000
0.000000
0.000000
0.000000
0.000000
~141.6864
0.000000
0.000000
0.000000
0.8273349
-32.98447
0.000000
1.639435
0.000000

0.2575156E-15
0.2217327E-15
-0.5442071E-15
0.000000
0.000000
0.000000
0.000000
0.2950176E-13
-0.1074050E-13
0.000000
0.000000
0.000000
0.000000
0.4487431E-13
-0.2178973E—13
0.000000

0.4923260
0.1642214
0.2524218

3.829388

0.1510141E-13
-0.1102911E-13

.000000
.142136
.000000
.000000
.000000
.000000
.433709
.000000
0.2497318

Cub-OOOOQO

 

 

Table 4.3.6. The extracted Taylor transfer map for the Super-FRS quadrupole. The

third order coefficients of the Taylor transfer map are shown.

125

300000
210000
120000
030000
201000
111000
021000
102000
012000
003000
200100
110100
020100
101100
011100
002100
200001
110001
020001
002001
100200
010200
001200
001101
100002
010002
001002
000300
000201
000102
000003

-0.2363760E-14-0.7663791E-14 0.000000 0.000000 -2.827376 400000
0.1297069E-13 0.6407553E-13 0.000000 0.000000 -4.523454 310000
0.5458716E-13 0.1335168E-12 0.000000 0.000000 -3.289673 220000
0.7936041E-13 0.1489082E-12 0.000000 0.000000 -1.611439 130000
0.5246487E-13 0.1048281E-12 0.000000 0.000000 -0.2505324 040000
0.2214878E-12 0.7629676E-12 0.000000 0.000000 -38.00841 202000

-0.9787271E-12-0.1752240E-11 0.000000 0.000000 -39.51802 112000

-0.2419524E-11-0.5075344E-11 0.000000 0.000000 -36.59077 022000
0.1620236E-11 0.2544577E-11 0.000000 0.000000 14.11275 000400

0.000000 0.000000 11.19360 66.04237 0.000000 000301
0.4103333E-13 0.6831805E-13 0.000000 0.000000 -7.455121 000202
0.000000 0.000000 -2.549064 -1.544657 0.000000 000103
0.000000 0.000000 0.000000 0.000000 -0.2185405 000004
1.572399 13.97146 0.000000 0.000000 0.6739899E-15 500000
18.97848 46.14729 0.000000 0.000000 0.2043959E-13 410000
46.30717 86.55577 0.000000 0.000000 0.6126932E-13 320000
40.04862 79.90497 0.000000 0.000000 0.7543286E-13 230000
18.17722 39.67904 0.000000 0.000000 0.4902103E-13 140000
2.335429 6.248855 0.000000 0.000000 0.5839144E-14 050000
0.000000 0.000000 52.34289 16.25874 0.000000 401000
0.000000 0.000000 “34.47570 -185.8744 0.000000 311000
0.000000 0.000000 -4.958071 -96.76526 0.000000 221000
-7.871988 7.690329 0.000000 0.000000 0.1030780E-12 100202
53.63167 53.02820 0.000000 0.000000 -0.4518421E-13 010202
0.000000 0.000000 -79.35436 -311.0261 0.000000 001202

-0.9529577E—13-0.1265356E-12 0.000000 0.000000 18.30773 001103

-0.7612110E-01 0.1853087 0.000000 0.000000 0.000000 100004

-0.2909687 -0.7612110E-01 0.000000 0.000000 0.000000 010004

0.000000 0.000000 1.494996 0.6985110 0.000000 001004
0.000000 0.000000 -201.7121 -64.37471 0.000000 000500
-0.5328690E-11-0.6519420E-11 0.000000 0.000000 -52.19016 000401
0.000000 0.000000 -15.52484 -102.2664 0.000000 000302
-0.5614135E-13-0.8194864E-13 0.000000 0.000000 11.06390 000203
0.000000 0.000000 2.597161 1.494996 0.000000 000104
0.000000 0.000000 0.000000 0.000000 0.1966990 000005

 

 

:- -. r ".' l-f'O.f..I-w

 

Table 4.3.7. The extracted Taylor transfer map for the Super-FRS quadrupole. Some
of the fourth and fifth order terms of the Taylor transfer map are shown.

126

magnet of Super-F RS. We have also extracted the transfer map for the Super-F RS
quadrupole from the magnetic field data. In addition, we have also presented the re-
sults of the high order ion-optic simulations for the per-separator, the main-separator

and the energy buncher of the Super—FRS.

127

Appendices

128

 

 

APPENDIX A

Comparison of the GICOSY and
the COSY INFINITY beam

physics codes

Both the GICOSY [19] and the COSY INFINITY [23, 24] codes are widely used
computer codes to analyses and simulate the particle optical systems. The GICOSY
code is a combination of the GIOS [77] and COSY 5.0 [20] codes. The code GICOSY
computes the transfer maps to the fifth order utilizing the magnetic and electric
elements from the code COSY 5.0 [20]. We start out by first giving brief description
of the both codes. We then compare the maps computed by both the codes for
important magnetic and electric elements and summarize our findings.

The GICOSY code provides description of the ion optical system using matrices up
to fifth order, dipoles with inhomogenieties, quadrupoles and all kinds of multipoles
both magnetic and electrostatic, precise treatment of the fringing fields with fringing
field integral method or direct tracing of matrix elements using differential algebra use
of variables, fitting with different methods, plots of the system, trajectories, envelopes,

fields and of phase space distributions, and others. The GICOSY code creates matrix

129

 

 

files of all optical elements that can be used in the Monte-Carlo simulation program
MOCADI [1].

The code COSY INFINITY is an arbitrary order beam dynamics simulation and
analysis code. It allows the study of accelerator lattices, spectrographs, beam lines,
electron microscopes, and many other devices. It can determine high-order maps of
combinations of particle optical elements of arbitrary field configurations. The ele-
ments can either be based on a large library of existing elements with realistic field
configurations including fringe fields, or described in detail by measured data. Analy-
sis options include computation of high-order nonlinearities; analysis of properties of
repetitive motion via chromaticities, normal form analysis, and symplectic tracking;
analysis of single—pass systems resolutions, reconstructive aberration correction, and
consideration of detector errors; and analysis of spin dynamics via computation of
spin maps, spin normal form and spin tracking.

Differences in approach and implementation:

The code COSY INFINITY code can in principle computes the transfer maps for
the particle optical elements to arbitrary order, restricted only by the computational
resources available. The coefficients of the transfer maps are computed to machine
precision using the differential algebraic techniques. In the GICOSY code the transfer
maps for most of the particle optical elements are computed to the maximum of 5th
order. But, for some of the optical elements the transfer maps are only computed to
lower orders. The GICOSY code uses TPSA [7] package to compute the composition
of the maps and other operations. However, the coefficients of a transfer map for an
individual particle optical element is computed using explicit formulas. The constants
appearing in these formulas were entered to only ten digit precision. Hence, the
overall accuracy of the maps computed using these coefficients can only be in the
order ~ 10-10. We expect to see differences of this magnitude in our comparison of

the transfer maps generated using the GICOSY and the COSY INFINITY codes.

130

 

Parameters GICOSY and COSY INFINITY
Order of computation 5

Fringe Fields Off

Output Coordinates Symplectic coordinates

witfiEnergy as coordinates

caling Scaled in time

amu = 1.6605402 x 10-27 kg
e0 = 1.6021773349 x 10-19 0

Natural constants CO = 299792458 X 108 IDS—1
amu = 931.4943307 MeV
71 = 3.141592653589793

Table A.0.1. System parameters used for the computation of the transfer maps.

 

 

 

 

 

 

 

 

 

 

Comparison of maps for important particle optical elements:

To compare the equivalent particle optical elements in the GICOSY and the COSY
INFINITY codes, we look at the difference of the maps computed using both the
codes. These maps are computed to the 5th order, the maximum order supported
by the GICOSY code. It would be impractical and also of little use to mention all
the coefficients of this difference map and study each coefficient separately. For the
purpose of comparison of the maps it would be more useful and practical to consider
the RMS average (63:23) of the coefficients of this difference map. This quantity
would clearly highlight the differences in the maps, if any. In an ideal case where the
equivalent particle optical elements in both codes agree perfectly for all the coefficients
of the difference map, the RMS average (675.2213) will be of magnitude 3 10-10.

Prior to computing the transfer maps all the system parameters are adjusted to
reflect same setting in both the codes. These parameters are listed in Table A.O.l.
The GICOSY code was run with computation mode 1, in which the transfer matrixes
of the optical elements are calculated by routines generated by the program Hamilton
[26]. Where as the COSY INFINITY code was run with default settings.

Table A.O.2 shows the result of the comparison of drift, magnetic and electric
particle Optical elements. The first column shows the particle optical element. The
second column shows the difference between the GICOSY and the COSY INFINITY

codes for the 5th order computation of the transfer map. And, the third column

131

 

 

provides comment about the observed difference.

The transfer maps for many of the elements in the GICOSY and COSY INFINITY
codes are same. But, for some we notice that the maps agree to only lower orders,
suggesting that in the GICOSY code these elements were implemented to lower orders.
Fiom Table A.O.2 we see that electric hexapole is implemented to the fourth order in
the GICOSY code. And, the electric quadrupole and the magnetic multipole elements
are implemented to third order in the GICOSY code. And, the magnetic sector is
implemented to only second order in the GICOSY code. Also, it can be seen clearly
from Table A.O.2 that the description of the electric sector and electric multipole
elements are not same in the GICOSY and the COSY INFINITY codes. Further
analysis is needed to find the possible reasons for this difference. Both the GICOSY
and the COSY INFINITY codes provide several options for specifying the fringe field
for the magnets. But, no options were found that would provide exactly the same

description of fringe fields.

132

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Drift 2.74 x 10'12 Identical to 5th order
Magnetic Elements

Dipole 2.21 x 10‘10 Identical to 5th order
Magnetic Sector 3.10 x 10“01 Identical to 2nd order
Quadrupole 3.08 x 10‘11 Identical to 5th order
Hexapole/Sextupole 2.23 x 10"Tr Identical to 5th order
Octupole 7.43 x 10‘12 Identical to 5th order
DecaPole 9.05 x 10'09 Identical to 5th order
Duodecapole 2.24 x 10'“07 Some differences
Multipole 1.86 x 10-01 Identical to 3rd order
Electric Elements

Electric Sector 3.19 x 10—OI Not Identical
Quadrupole 1.54 x 10"01 Identical to 3rd order
Hexapole/Sextupole 6.93 x 10—02 Identical to 4th order
Octupole 2.75 x 10‘12 Identical to 5th order
DecaPole 7.03 x 10"10 Identical to 5th order
Duodecapole 5.66 x 10—11 Identical to 5th order
Multipole 722.79 Not Identical

 

 

 

Table A.O.2. The first column shows the partical optical element. The column 2 I]
shows the RMS average of the coefficients of the difference map between GICO and 1' ‘
COSY INFINITY for the 5th order computation. And the third column provides
comments about the observed difference.

 

133

APPENDIX B

The code converter

The converter program is based on the current release of the GICOSY code [19] (as of
Sept. 2005) and the COSY INFINITY version 8.1 [23, 24] code. The converter code
is implemented using the Perl programming language. The important beamline ele
ments are translated into the respective elements in the code COSY INFINITY; these
include drifts, multipoles, superimposed multipoles, and bends. Table B.0.1 shows
the equivalent COSY INFINITY procedure calls for the important GICOSY particle
optical elements. Some elements supported by the GICOSY code are translated to
drifts or commented out and may have to be adjusted manually. All comments in
the GICOSY code are converted to comments in the COSY INFINITY code. The
GICOSY code does not require the name and size of the variables to be declared. On
the other hand, in the code COSY INFINITY needs all the variables to be declared.
In the current converter most of the variable name are extracted from the GICOSY
code and are declared in the converted COSY INFINITY code. The loops and fitting
routines available in the GICOSY code are commented out in the converted COSY
INFINITY code. These have to be manually converted to the equivalent COSY IN-
FIN ITY code.

134

 

var
var
var . '

 

Table B.0.1.
the important GICO particle optical elements.

procedure calls for

 

135

APPENDIX C

New COSY INFINITY tools for

beam physics simulations

For the convenience of presenting and analyzing the spectrometer simulation results

the following tools have been implemented in COSY INFINITY:

1. The beam profile generator: This would set the beam profile to either an elliptic
or a rectangular shape, which can then be used for phase space plot and the
transmission plot. An elliptic or a rectangular distribution of the particles is

generated using the pseudo random number generator available in Fortran 77.

2. The phase space plot: This can be used to plot the phase space by choosing
x-a or y-b. In addition this can also be used to plot between any two particle

optical coordinates.

3. The transmission plot: This tool launches a beam of 1000 randomly chosen
particles and looks at transmission through the beam line. The loss comes from
the limited aperture of the elements. Final coordinates after each element are
computed and compared with the aperture of the element. If the aperture is

smaller then the particle coordinates the particle is considered lost. This tool

136

 

.»--J

 

0.3 E-l

0.9 E-l

 

(a) Before correcting aberrations

 

0.3 E-l

M

 

(b) After correcting aberrations

Figure C01. The x-a phase space plots at the dispersive focal plane PF2 for a beam
of 40 mm mrad and three different momenta of Ap/ p = :l:2.5 %. The plot (a) shows
the phase space before second order aberration corrections, and (b) shows the phase
space after second order aberration corrections.

137

100 1 1 1 1

 

 

 

 

 

 

ist order trackir'ig ’ _—
2nd order tracking .......
3nd order tracking ........
4th order tracking ~---—-----
5th order tracking -4------
80 ’- ‘~:..'.:.:..'.t..'.::..'..'.7..'..'.:.7..1... ..1 ------ _
l I .. I - - I "a3423;113:551235.1. 11:17:...
"'""""“--'.‘.‘.::.'.‘.“.‘1::;{---.
*1
60 ~ (1'1“ q
o\°
.5
(I)
.9
1..
20 ~ _
0 l l l P 1 1 J; 1
0 20 40 60 80 100 120 140 160 1 30

Length/m

Figure C.0.2. Ion-optic transmisson plot for the Super-FRS high energy branch.

 

can be modified to include slits and other elements that influence the trans—

mission in a trivial manner. All transmission losses come from purely particle

optical effects and is not suitable for modeling reaction products. Example of

transmission plot for the Super-FRS high energy branch is shown in the Figure

C02.

4. Beam envelope and matrix element plots: This plots the beam envelOpe and a

matrix element versus the length along the optical axis. Figure 4.3.3 shows an

example of the beam envelope and the matrix element plot.

5. MOCADI input. This tool generates the COSY INFINITY transfer map output
in the GICOSY output format. This can then be used to create input for the
Monte Carlo simulation code MOCADI. The input file for MOCADI is generated
using the GICOSY and transfers maps from the COSY INFINITY.

138

 

BIBLIOGRAPHY

[1] MOCADI 3.0. http://www-linux.gsi.de/ weick/mocadi/. www html page.

[2] Seaborg - IBM SP. http://www.nersc.gov/nusers/resources/SP/. www html
page.

[3] The Tosca Reference Manual. Technical Report, Vector Fields Limited, 24 Bank-
side, Kidlington, Oxford, OX5 1JE, England.

[4] The Tosca User Guide. Technical Report, Vector Fields Limited, 24 Bankside,
Kidlington, Oxford, OX5 IJE, England.

[5] The Large Hadron Collider - Conceptual Design. Technical Report AC/95—05
(LHC), CERN (1995).

[6] A. ASI. Boundary element method (bem) for charged particle optics. vol. 4510,
pp. 138—147. SPIE (2001).

[7] M. BERZ. The new method of TPSA algebra for the description of beam dy-
namics to high orders. Technical Report AT-61ATN-86-16, Los Alamos National
Laboratory (1986).

[8] M. BERZ. The method of power series tracking for the mathematical description
of beam dynamics. Nuclear Instruments and Methods A258, 431—437 (1987).

[9] M. BERZ. High-order description of accelerators using differential algebra and
first applications to the SSC. In “Proceedings, Snowmass Summer Meeting”,
Snowmass, Co (1988).

[10] M. BERZ. Differential algebra - a new tool. Technical Report LBL—27033,
Lawrence Berkeley Laboratory, Berkeley, CA (1989).

[11] M. BERZ. Differential algebra - a new tool. In “Proceedings of the 1989 IEEE
Particle Accelerator Conference”, p. 1419. IEEE (1989).

[12] M. BERZ. Differential algebraic description of beam dynamics to very high
orders. Particle Accelerators 24, 109 (1989).

139

 

 

[13] M. BERZ. Computational aspects of design and simulation: COSY INFINITY.
Technical Report 740, National Superconducting Cyclotron Laboratory, Michi-
gan State University, East Lansing, MI 48824 (1990).

[14] M. BERZ. Computational aspects of optics design and simulation: COSY IN-
FINITY. Nuclear Instruments and Methods A298, 473 (1990).

[15] M. BERZ. COSY INFINITY, an arbitrary order general purpose optics code.
Computer Codes and the Linear Accelerator Community Los Alamos LA-
11857-0, 137 (1990).

[16] M. BERZ. COSY INFINITY. In “Proceedings 1991 Particle Accelerator Con-
ference”, San Francisco, CA (1991).

[17] M. BERZ. “Modern Map Methods in Particle Beam Physics”. Academic Press,
San Diego (1999). Also available at http: / / bt.pa.msu.edu/ pub.

[18] M. BERZ, B. FAWLEY, AND K. HAHN. High order calculation of the multipole
content of three dimensional electrostatic geometries. Nuclear Instruments and

Methods A307, 1 (1991).

[19] M. BERZ, B. HARTMANN, K. LINDEMANN, A. MAGEL, AND H. WEICK. GI-
COSY Ion Optical Program. http://www-linux.gsi.de/ weick/gico/. www html

page.

[20] M. Benz, H. C. HOFMANN, AND H. WOLLNIK. COSY 5.0, the fifth order
code for corpuscular optical systems. Nuclear Instruments and Methods A258,
402—406 (1987).

[21] M. BERZ, K. JOH, J. A. NOLEN, B. M. SHERRILL, AND A. F. ZELLER. Re-
constructive correction of aberrations in nuclear particle spectrographs. Physical
Review C 47,2, 537 (1993).

[22] M. BERZ AND K. MAKINO. Verified integration of ODEs and flows using differ-
ential algebraic methods on high-order Taylor models. Reliable Computing 4(4),

361—369 (1998).

[23] M. BERZ AND K. MAKINo. COSY INFINITY beam physics manual. Technical
Report MSUHEP-60804, Department of Physics and Astronomy, Michigan State
University, East Lansing, MI 48824 (2006). see also http://cosy.pa.msu.edu.

[24] M. BERZ AND K. MAKINO. COSY INFINITY programmer’s manual. Technical
Report MSUHEP-60803, Department of Physics and Astronomy, Michigan State
University, East Lansing, MI 48824 (2006). see also http://cosy.pa.msu.edu.

140

[25]

[26]

[27]

[28]

[29]

[30]

[31]

[32]

[33]

[34]

M. BERZ, K. MAKINo, K. SHAMSEDDINE, C. HOFFSTATTER, AND W. WAN.
COSY INFINITY and its applications to nonlinear dynamics. In M. BERZ,
C. BISCHOF, G. CORLISS, AND A. GRIEWANK, editors, “Computational Dif-
ferentiation: Techniques, Applications, and Tools”, pp. 363—365, Philadelphia
(1996). SIAM.

M. BERZ AND H. WOLLNIK. The program HAMILTON for the analytic solution
of the equations of motion in particle optical systems through fifth order. Nuclear
Instruments and Methods A258, 364—373 (1987).

K. L. BROWN. A first- and second-order matrix theory for the design of beam
transport systems and charged particle spectrometers. Technical Report 75,
SLAC (1982).

Y. F. CHANG AND C. F. CORLISS. ATOMFT: Solving ODEs and DAEs using
Taylor series. Computers and Mathematics with Applications 28, 209—233 (1994).
see also http://www.mscs.mu.edu/ georgec/Pubs/journ.html#1994a.

A. CUNSOLO, F. CAPPUZZELLO, A. V. BELOZYOROV, A. ELANIQUE, A. FOTI,
A. LAZZARO, O. MALISHEV, L. MELITA, W. MITTIG, C. Nocrrono,
P. ROUSSEL-CHOMAZ, V. SHCHEPUNOV, D. VINCIGUERRA, A. YEREMIN,
AND J. S. WINFIELD. Magnex: A large acceptance magnetic spectrometer for
EXCY T. Nuclear and Condensed Matter Physics - Sixth Regional Conference
513(1), 270—273 (2000).

A. M. DAVIS. A generalized helmholtz theorem for time-varying vector fields.
American Journal of Physics 74(1), 72—76 (2006).

R. DEGENHARDT AND M. BERZ. High accuracy description of the fringe fields of

particle spectrographs. In “Proceedings 1993 Particle Accelerator Conference”,
Washington, DC (1993).

R. DEGENHARDT AND M. BERZ. High accuracy description of the fringe field
in particle spectrographs. Nuclear Instruments and Methods A427, 151-156
(1999).

F. M. DENARO. On the application of the helmholtz-hodge decomposition in
projection methods for incompressible flows with general boundary conditions.

International Journal for Numerical Methods in Fluids 43(1), 43—69 (2003).

B. ERDELYI. “Symplectic Approximation of Hamiltonian Flows and Accurate
Simulation of Hinge Field Effects”. PhD thesis, Michigan State University, East
Lansing, Michigan, USA (2001).

141

[35] B. ERDELYI, M. BERZ, AND K. MAKINO. Detailed analysis of fringe field
effects in the large hadron collider. Technical Report MSUCL-1129, National
Superconducting Cyclotron Laboratory, Michigan State University, East Lans-
ing, MI 48824 (1999).

[36] L. C. EVANS. “Partial Differential Equations”. American Mathematical Society,
Rhode Island (1998).

[37] S. FARTOUKH AND O. BRU N ING. Field Quality Specification for the LHC Main
Dipole Magnets. Technical Report LHC Project Report 501, CERN (2001).

[38] H. GEISSEL, M. WINKLER, H. WEICK, K.-H. BEHR, G. MUNZENBERG, H. S1-
MON, K. SUMMERER, B. ACHENBACH, L. AHLE, T. AUMANN, J. AYSTO,
M. BERZ, A. BLEILE, D. BOUTIN, G. BREITENBERGER, A. BRUNLE,
P. DENDOOVEN, H. EMLING, G. FEHRENBACHER, M. GLEIM, W. HULLER,
H. IWASE, A. JOKINEN, C. KARAGIANNIS, R. KAISER, A. KELIC,
B. KINDLER, G. KLAPPICH, T. KUBO, N. KURZ, K. KUSAKA, H. LEIBROCK,
J. LETTRY, B. LOMMEL, S. MANIKONDA, I. MOORE, G. MORITZ, D. MOR-
RISSEY, C. MUHLE, J. A. NOLEN, I. PSCHORN, T. RADON, R. RON-
NINGEN, G. SAVARD, C. SCHEIDENBERGER, Y. SHATUNOV, B. SHERRILL,
M. SVEDENTSOV, N. TAHIR, K. WENDT, M. YAVOR, A. YOSHIDA, X. YU,
A. ZELLER, AND O. ZURKAN. Technical Report on the Design, Construction,
Commissioning and Operation of the of FAIR. Technical Report, G81 (2005).

[39] J. GROTE, M. BERZ, AND K. MAKINO. High-order DA methods for the deter-
mination of Poincare sections. Nuclear Instruments and Methods (in print).

[40] H. HELMHOLTZ. Uber Integrale der hydromechanischen Gleichung, welche den
Wirbelbewegungen entsprechen. J. reine angew. Math. (Crelle’s Journal) 55, 25

(1858).

[41] K. N. HENRICHSEN. Overview of magnet measurement methods. In “CERN
Accelerator School 1997” (1997).

[42] J. HOEFKENS. “Verified Methods for Differential Algebraic Equations”. PhD
thesis, Michigan State University, East Lansing, Michigan, USA (2001).

[43] J. HOEFKENS, M. BERZ, AND K. MAKINO. Efficient high-order methods for
ODEs and DAEs. In G. CORLISS, C. FAURE, A. GRIEWANK, L. HASCOET,
AND U. NAUMANN , editors, “Automatic Differentiation of Algorithms from Sim-
ulation to Optimization”, pp. 343—348. Springer (2002).

[44] J. D. JACKSON. “Classical Electrodynamics”. Wiley, New York (1975).

142

[45] J .IRWIN AND A.DRAGT. “Taylor Maps, Entry in ’Handbook of Accelerator
Physics and Engineering’, M. Tigner and A. Chao (Eds.)”. World Scientific,
New York (1999).

[46] Y.-K. KIM AND M. BERZ. Parallel constructs in COSY INFINITY. Technical
Report MSUHEP-060805, Michigan State University (2006).

[47] E. R. KOLCHIN. “Differential Algebra and Algebraic Groups”. Academic Press,
New York (1973).

[48] I. KUBLER. Master’s thesis, Justus Liebig Universitéit GieBen, 6300 GieBen,
West Germany (1987).

[49] A. LAZZARO, A. CUNSOLo, F. CAPPUZZELLO, A. FOTI, C. NOCIFORO,
S. ORRIGO, V. SHCHEPUNOV, J. WINFIELD, AND M. ALLIA. Computational
aspects of the trajectory reconstruction in the MAGNEX large acceptance spec-
trometer. In BERZ AND MAKINO, editors, “Computational Accelerator Physics
2002”. Institute of Physics Conference Series 175 (2004).

[50] K. MAKINO. “Rigorous Analysis of Nonlinear Motion in Particle Accelerators”.
PhD thesis, Michigan State University, East Lansing, Michigan, USA (1998).
Also MSUCL-1093.

[51] K. MAKINO. The COSY 8th order Runge Kutta integrator. Neu-
trino Factory/Muon Collider Notes MUC-NOTE-COOL-THEORY—
0238, Fermi National Accelerator Laboratory (2002). see http://www—
mucool.fnal. gov / notes / notes.html. '

[52] K. MAKING AND M. BERZ. Recent applications of COSY to nonlin-
ear beam dynamics problems. In “Proc. of the APS/DPF/DPB Summer
Study on the Eiture of Particle Physics (Snowmass 2001)”, no. T510 (2002).
http: / / www.slac.stanford.edu / econf / C010630 / papers / T510.PDF.

[53] K. MAKING AND M. BERZ. Taylor models and other validated functional in-
clusion methods. International Journal of Pure and Applied Mathematics 6,3,

239—316 (2003).

[54] K. MAKINO, M. BERZ, D. ERREDE, AND C. J. JOHNSTONE. High order

map treatment of superimposed cavities, absorbers, and magnetic multipole and
solenoid fields. Nuclear Instruments and Methods A 519, 162—174 (2004).

[55] S. MANIKONDA, J. NOLEN, M. BERZ, K. MAKING, AND G. WEIZMAN. Design
of a superconducting quadrupole with elliptical acceptance and tunable higher
order multipoles. Technical Report MSUHEP-060412, Michigan State University
(2006).

143

[56] S. L. MANIKONDA AND M. BERZ. A highly accurate high-order method to

solve the Neumann boundary condition problem for the 3D Laplace equation.
International Journal of Pure and Applied Mathematics 23,3, 365—378 (2005).

[57] S. L. MANIKONDA AND M. BERZ. Multipole expansion solution of the Laplace
equation using surface data. Nuclear Instruments and Methods 258, 175—183

(2006).

[58] S. L. MANIKONDA, M. BERz, AND K. MAKINO. A highly accurate high-order
validated method to solve the Neumann boundary condition problem for the 3D
Laplace equation. Transactions on Computers 4-11, 1604—1610 (2005).

[59] P. M. MORSE AND H. FESHBACH. “Methods of Theoretical Physics, Part I and
Ira (1953).

[60] N. F. OSMAN AND J. L. PRINCE. 3d vector tomography on bounded domains.
Inverse Problems 14(1), 185-196 (1998).

[61] H. B. PHILLIPS. “Vector analysis”. J. Wiley & Sons, inc.; London, Chapman &
Hall, limited, New York (1933).

[62] R. PLONSEY AND R. E. COLLIN. “Principles and Applications of Electromag-
netic Fields”. McGraw Hill, New York (1961).

[63] P.L.WALSTROM. Soft-edged magnet models for higher-order beam-optics map
codes. Nuclear Instruments and Methods A519 1-2, 216—221 (2004).

[64] P. J. PRINCE AND J. R. DORMAND. High order embedded Runge—Kutta for-
mulae. Journal of Computational and Applied Mathematics 7, 67—75 (March
1981)

[65] S. I. REDIN, N. M. RYSKULOV, G. V. FEDOTOVICH, B. I. KHAZIN, G. M.
BUNCE, G. T. DANBY, J. W. JACKSON, W. M. MORSE, R. PRIGL, AND
Y. K. S. ET AL. Radial magnetic field measurements with a Hall probe device

in the muon (g—2) storage ring magnet at BNL. Nuclear Instruments and Methods
A 473, 3, 260—268 (2001).

[66] N. REVOL, K. MAKINO, AND M. BERZ. Taylor models and floating-point
arithmetic: Proof that arithmetic operations are validated in COSY. Journal of
Logic and Algebraic Programming 64/ 1, 135—154 (2004).

[67] J. F. RITT. “Differential Algebra”. American Mathematical Society, Washing-
ton, D.C. (1950).

[68] G. SABBI. Magnetic field analysis of HGQ coil ends. Technical Report TD—97—
040, Fermilab (1997).

144

[69] M. L. SHASHIKANT, M. BERz, AND B. ERDELYI. COSY INFINITY’s EXPO
symplectic tracking for LHC. IOP CP 175, 299—305 (2002).

[70] V. A. SHCHEPUNOV, A. CUNSOLO, F. CAPPUZZELLO, A. FOTI, A. LAZZARO,
A. L. MELITA, C. NOCIFORO, AND J. S. WINFIELD. Numerical computation
of arbitrary order transfer maps and reconstructive correction of aberrations in
the large acceptance spectrometer MAGNEX. Nuclear Instruments and Methods
B 204, 447—453 (May 2003).

[71] I. VASSERMAN AND S. SASAKI. Comparison of different magnetic measure-
ment techniques. In “Proceedings of 13th International Magnetic Measurement

Workshop” (2003).

[72] M. VENTURINI, D. ABELL, AND A. DRAGT. Map computation from mag-
netic field data and application to the LHC high-gradient quadrupoles. eConf
C980914, 184—188 (1998).

[73] M. VENTURINI AND A. DRAGT. Computing transfer maps from magnetic field
data. In “IEEE Particle Accelerator Physics Conference” (1999).

[74] M. VENTURINI AND A. J. DRAGT. Accurate computation of transfer maps from
magnetic field data. Nuclear Instruments and Methods A427, 387—392 (1999).

[75] P. WALSTROM, A. DRAGT, AND T. STASEVICH. Computation of charged-
particle transfer maps for general fields and geometries using electromagnetic
boundary-value data. In “IEEE Particle Accelerator Conference 2001” (2001).

[76] E. W. WEISSTEIN. Laplace’s Equation. From MathWorld—A Wolfram Web
Resource. http://mathworld.wolfram.com/LaplacesEquation.html. www html

page.
[77] H. WOLLNIK, J. BREZINA, AND M. BERz. CIOS—BEAMTRACE, a program

for the design of high resolution mass spectrometers. In “Proceedings AMCO-7”,
p. 679, Darmstadt (1984).

[78] H. WOLLNIK, J. BREZINA, AND M. BERz. GIOS—BEAMTRACE, a computer
code for the design of ion optical systems including linear or nonlinear space
charge. Nuclear Instruments and Methods A258, 408 (1987).

[79] H. WOLLNIK, B. HARTMANN, AND M. BERZ. Principles behind CIOS and
COSY. AIP CP 177, 74 (1988).

[80] D. A. WOODSIDE. Uniqueness theorems for classical four-vector fields in Euclid—
ean and Minkowski spaces. Journal of Mathematical Physics 40( 10), 4911—4943
(1999).

145

[81] D. A. WOODSIDE. Classical four-vector fields in the relativistic longitudinal
gauge. Journal of Mathematical Physics 41(7), 4622—4653 (2000).

[82] K.-H. YANG. The physics of gauge transformations. American Journal of
Physics 73(8), 742—751 (2005).

146

   

u[[1]][1][][]]][[[5]]