LOS ALAMOS NATIONAL LABORATORY
LA-UR-02-4979

8/1/02T H E O RY A N D
T E C H N I QU E O F B E A M
ENVELOPES I M U L AT I O NS IMULATION OF BUNCHED PARTICLE
B EAMS WITH E LLIPSOIDAL S YMMETRY AND
L INEAR S PACE C HARGE F ORCES
Christopher K. Allen and Nicholas D. Pattengale
SNS-4, Los Alamos National Laboratory
Los Alamos, New Mexico1

LA-UR-02-4979
Table of Contents
1.

2.

3.

4.

5.

6.

7.

Introduction _____________________________________________________________________ 4
1.1.

Overview____________________________________________________________________ 4

1.2.

Outline _____________________________________________________________________ 5

Physics Background _______________________________________________________________ 5
2.1.

The Design Trajectory_________________________________________________________ 5

2.2.

Phase Space _________________________________________________________________ 6

2.3.

The Longitudinal Phase Plane __________________________________________________ 6

2.4.

Equations of Motion __________________________________________________________ 7

Statistics Background ______________________________________________________________ 9
3.1.

The Density Function _________________________________________________________ 9

3.2.

Liouville's Theorem___________________________________________________________ 9

3.3.

Moments of the Distribution___________________________________________________ 10

3.4.

Moment Equations __________________________________________________________ 10

Beam Physics ___________________________________________________________________ 11
4.1.

The RMS Envelopes _________________________________________________________ 11

4.2.

RMS Emittance _____________________________________________________________ 12

4.3.

Courant-Snyder Parameters for the RMS Ellipse _________________________________ 14

4.4.

The Mean Value Vector and Covariance Matrix __________________________________ 15

4.5.

Linearization of the Self-Fields_________________________________________________ 16

4.6.

Tune Depression ____________________________________________________________ 17

Ellipsoidal Beams and the Equivalent Uniform Beam ___________________________________ 18
5.1.

Ellipsoidally Symmetric Charge Distributions ____________________________________ 18

5.2.

Computation of the Field Moments _____________________________________________ 19

5.3.

The Equivalent Uniform Beam_________________________________________________ 20

Bunched Beam Envelope Equations _________________________________________________ 22
6.1.

Equations for Centroid Motion ________________________________________________ 22

6.2.

Bunched Beam RMS Envelope Equations________________________________________ 22

6.3.

Ellipsoidally Symmetric Beams and the Equivalent Uniform Beam___________________ 23

6.4.

Approximate Form of the Envelope Equations____________________________________ 24

Transfer Matrix Approach _________________________________________________________ 26
7.1.

Transfer Matrices and External Elements _______________________________________ 26

7.2.

Space Charge Impulses _______________________________________________________ 28

7.3.

Equations for Centroid Motion ________________________________________________ 29

7.4.

Equations for the Second-Order Moments _______________________________________ 29

7.5.

Extensions to the Inhomogeneous Case __________________________________________ 30

7.6.

Homogeneous Coordinates ____________________________________________________ 31

2

LA-UR-02-4979
8.

Summary and Conclusion _________________________________________________________ 33

Acknowledgements ___________________________________________________________________ 33
Appendix A: Transfer Matrices for Common Beamline Elements______________________________ 34
A.1

Drift Space _________________________________________________________________ 34

A.2

Quadrupole ________________________________________________________________ 34

A.3

Dipole Steering Magnet (As a Thin Lens) ________________________________________ 35

A.4

Dipole Steering Magnet (As a Thick Lens) _______________________________________ 35

A.5

Bending Magnet_____________________________________________________________ 37

A.6

RF Gap (As a Thin Lens) _____________________________________________________ 39

Appendix B: Self-Field Moment Calculation ______________________________________________ 41
References__________________________________________________________________________ 43

3

LA-UR-02-4979

1. INTRODUCTION
A convenient way to describe particle beam behavior is through the use of statistics. For instance,
instead of a full multiple-particle simulation we can propagate the statistical properties of the beam down
the beamline, such as its various statistical moments. Of course such an approach does not provide as
much information as a multiple-particle simulation, however, it is computationally much less expensive.
The approach is that taken by the simulation codes TRANSPORT [2] and TRACE3D [5]. The current note
describes this type of simulation technique in some detail, in particular, the beam physics background,computation of the self-force effects on the dynamics (i.e., space charge effects), and numericalapproaches.
It is our intention to describe in some detail the theory and technique for simulation of particle beamenvelopes. We hope that this unified presentation will serve as a collective documentation for suchtechniques that are already in common practice. By fully documenting both the theory and numerical
techniques one can more readily develop extensions to the theory and techniques and draw conclusions
toward the validity of the resulting simulations. For example, we could extend the theory here to include
noise from machine errors, or use the theory included here to determine which regimes the simulations aremost accurate. Because our interests are primarily in the control and automation of accelerator systems, it
is important to us that we have a clear understanding of the simulation techniques being employed todevelop robust control applications. The techniques presented here are fast. Thus, not only are they usefulfor off-line accelerator system design, they are also valuable for presenting on-line model references fordeveloping control applications.
Although we cover the theory and technique of beam envelope simulation in detail, we mention
nothing about how one might actually implement such a simulation engine. In a companion article we
present a modern software architecture for performing these simulations [16]. This architecture has alreadybeen implemented and is in current use, embedded into the high-level control system for the Spallation
Neutron Source (SNS) accelerator.
1.1. Overview
We begin with the equations of motion for single particle trajectory under Lorentz forces. Theseequations are tailored to beam physics by using the path length along the design trajectory as theindependent variable, rather than time. By employing Liouville's theorem we find it possible to derive
evolution equations for the beam ensemble's moments from the equations for single particle motion. There
are two unresolved terms in these equations, both determining the effects of the beam's self fields on the
moment dynamics. Specifically, they are the cross-moments of the self electric field and the phase
coordinates. If we assume that the beam has ellipsoidal symmetry in configuration space, it is possible tocompute one of these moments analytically in terms of elliptic integrals (the moment involving the spatial
coordinate). Upon computing this moment we find that it depends only weakly on the actual profile of theellipsoidal distribution. Acknowledging this fact we typically model all laboratory beams with anequivalent uniform beam having the same second moments; the uniform beam is preferred since it has
well-defined boundaries. From there we take two different approaches in the description of the statisticalevolution of the beam.
From the moment equations we can derive a convenient set of equations that describe the behavior of
the rms beam envelopes (assuming a centered beam). We avoid the remaining unknown field moment byintroducing the definition of rms emittance in the equation set. By doing so, however, we implicitly makethe assumption that the rms emittance is either constant or it is a function whose values are known a priori.
This assumption is usually not too restrictive, since it is known that rms emittances are invariants of the
motion if all forces on the beam are linear. The resulting set of equations for rms envelopes can also be
rewritten for the equivalent uniform beam, and is usually seen that way in the literature. They form a
convenient closed set of equations for studying beam dynamics. However, as pointed out they cannot be
used to study nonlinear effects acting on the beam.
The alternate approach to developing a description of the beam's statistical evolution is through the useof transfer matrices. While the previous approach results in a set of coupled ordinary equations that are4

LA-UR-02-4979
continuous in nature, this technique is more of a discrete approach. Generally, the continuous
representation is most convenient for analytic study while the discrete representation is well suited forcomputer simulation. The transfer matrix technique is based upon the observation that whenever all forceson a particle are linear, multiplying the particle's coordinate vector by a matrix can represent the action of
any beamline element. This matrix is determined by the properties of the element and is called the transfer
matrix for the element. If we represent all the second-order moments of the beam by a symmetric matrix,
then the effects of the element on these statistical quantities can be determined by a transpose conjugation
with the transfer matrix. This technique is easy to implement on a computer. A transfer matrix can also be
used to model the action the beam's space charge if we use a linear fit to the self electric field of the beam.
We do this by employing a weighted linear regression, where the weighting factor is the beam distribution.
This gives us a fit that is more accurate in regions of higher beam density. It turns out that this fieldapproximation also conforms to the assumption of constant rms emittance, as done in the previous
description of the beam dynamics. So we essentially have an equivalent representation as before, only thebeam is propagated discretely with transfer matrices.
1.2. Outline
In Section 2 we review the basic physics background necessary for this analysis. Mostly this entails
placing the common physical concepts in the context of accelerator physics. In particular we cover thedesign trajectory, phase space, and the equations of motion for an accelerator system. Section 3 introduces
the basic statistical concepts on which we rely for our simulation technique. There we show how totransform the equations of motion for individual particles into evolution equations for statistical propertiesof the beam. Section 4 describes quantities and parameters particular to particle beam physics. The mostimportant of these concepts is probably the rms envelopes and the Courant-Snyder parameters fordescribing the phase space distribution of the beam. We also cover the mean value vector and thecovariance matrix, convenient methods for representing the statistical properties of the beam. Ellipsoidal
beams are treated in Section 5, that is, beam distributions having ellipsoidal symmetry. Specifically, we
compute the space charge effects for ellipsoidal beams. We see that it is possible to determine the effectsof the self-forces analytically in terms of elliptic integrals for distributions having ellipsoidal symmetry in
configuration space. Sections 6 and 7 cover simulation techniques for bunched beams with ellipsoidal
symmetry. Section 6 treats the envelope equations for bunched beams where we develop the coupled set ofordinary differential equations describing the evolution of the equivalent uniform beam. Section 7 presents
the transfer matrix approach to bunched beam envelope simulation. This technique, although somewhat
less obvious, is more suited to numerical simulation and is also more general than the envelope equation
approach as it can treat coupling between the phase planes. Finally we conclude with Section 8 and present
some transfer matrices for common beamline elements in Appendix A.

2. PHYSICS BACKGROUNDy

Here we present the basic physics
backgroundfor
ouranalysis.
Ultimately, we develop the equationsof motion for individual beam particles
and we do so in a form suitable forbeam physics. We also touch upon thevariousunit
conversions
andconventions used in the beam physics
literature.
2.1. The Design Trajectory
We shall assume a beam transport
or accelerator system that has a
specified design trajectory.
This
configuration is shown in Figure 1.
The distance along this design

synchronous particle

beam
b

x

desig
n
unch

trajectoryz

Figure 1: design trajectory and coordinate system5s

LA-UR-02-4979
trajectory is given by the independent variable s. There also exists the synchronous particle, which has a
specified design velocity v(s) at each point s on the design trajectory. The relative velocity β=v/c andrelativistic parameter γ=1/(1-β2)1/2 are always given with respect to this design velocity, unless otherwise
noted. With respect to the synchronous particle, we construct a system of coordinates (x,y,z), that is, thesynchronous particle is at the origin. The coordinates x,y,z represent displacements from the synchronous
particle in the x, y, z directions, respectively. Locally, the z-coordinate is always aligned with the design
trajectory. Specifically, the tangent vector of the design trajectory always points in the z-direction. Thus,
the xy-plane represents the transverse plane while the z-direction is the longitudinal direction of beam
propagation (in a local sense). Note that the coordinates (x,y,z) are not the inertial frame of the beam, they
are laboratory coordinates that follow the beam.
2.2. Phase SpaceWe form the phase space (or state space) of the particle by considering the normalized momenta(x',y',z') with respect to the synchronous particle. Let p(s)=γmv(s) represent the momentum magnitude of
the synchronous particle. Then the x and y plane relative momentum x' and y' are given by
dx p x γmx&
=
=
,
ds
p
γms&
dy p y γmy&
y' ≡
=
=
,
ds
p
γms&

x' ≡

(1)

where the over dot denotes differentiation with respect to time. For the z plane the situation is different
since the coordinate z is defined to be the difference in longitudinal position from the synchronous particle,
which is traveling at velocity v. Therefore,dz∆v ( ∆s / v )∆ v∆ t
∆z
= lim
= lim
= lim
∆s →0
ds ∆s→0 ∆s ∆s→0 ∆s
∆s
∆v
1 ∆p
=
= 2
,
v
γ p

z′ ≡

(2)

where ∆v is the difference in velocity v of the synchronous particle, ∆p is the difference in longitudinal
momentum p of the synchronous particle, and the last equality comes from relativistic mechanics (seebelow). Thus, the complete phase space coordinates for a particle at location s is given by (x,x',y,y',z,z';s).
For conciseness, it is convenient to denote points in phase space by the vector z. That is,(3)z ≡ ( x, x ′, y, y ′, z, z ′) .

This coordinate space, specifically with the normalized momenta x', y', and z', is also commonly called
trace space in the literature.
2.3. The Longitudinal Phase Plane
Since the z phase plane has several different descriptions in the literature, we shall consider it in moredetail. There are several different coordinate systems commonly used to describe this phase plane.
Because we are considering position and momentum in the direction of propagation it is necessary toconsider relativistic effects when converting between these descriptions. The following equations relatedifferential changes in common relativistic parameters:dγ = βγ 3 dβ ,(4)

d ( βγ ) = γ 3 dβ ,dW = mc 2 dγ = βγ 3mc 2 dβ ,where W is the kinetic energy of the particle given by (γ-1)mc2.

6

LA-UR-02-4979
Commonly used coordinates used in lieu of z are ∆φ and ∆t. The coordinate ∆φ is the RF phasedifference from the synchronous particle, while ∆t is the difference in arrival time from the synchronous
particle. To convert from ∆φ in degrees, or ∆t in seconds, to z in meters use the following formulae:
(5)

∆φ
∆φ
,
= − βλ360
360= − ∆t ( v + ∆v ) = − βc∆t (1 + z ′) ,z = − vT

where T is the period of the RF and λ is the free-space wavelength of the RF. Note that negative values of∆φ indicate a phase difference ahead of the synchronous particle while negative values for ∆t indicate
arrival times before the synchronous particle.
Coordinates that alternately describe z' are ∆p/p and ∆W. The quantity ∆p is the difference inmomentum from the synchronous particle so that ∆p/p is the normalized difference in momentum.
Likewise, ∆W is the difference in kinetic energy from the synchronous particle. To convert from ∆p/p inradians or ∆W in electron volts (Joules) to z' in radians use the relation
z' =

(6)
z′ =

1 ∆p

γ2 p

,

1
∆W
,
γ (γ + 1) W

where W is the beam's kinetic energy.
These relations∆p=mc∆(βγ)=γ3mc∆β=pγ2(∆β/β) and ∆W=mc2∆γ=β2γmc2(∆p/p).

follow

from

the

facts

that

2.4. Equations of Motion
Newton's second law gives the equations of motion for individual particles of the beam. We shallmake all our calculations in the laboratory frame. Recall that the coordinates (x,y,z) are actually thedisplacements from the synchronous particle on the design trajectory. In our case the familiar F=dp/dt
appears as
F=(7)d (γmv )d= βc (γmv )dtds= β 2γmc 2 r ′′ + β 2 mc 2 ( x ′, y ′, z ′ + 1)dγds

where F is the force on a particle, γ is the relativistic factor, m is the particle mass, c is the speed of light,
β≡v/c where v is velocity of the synchronous particle along s, r≡(x,y,s+z) is the position vector, and theprimes indicate differentiation with respect to s. To derive the second line of the above note that thevelocity vector can be expressed as v=βc(x',y',z'+1). The second term of that equation is nonzero only
when particles are accelerated. Many times we can assume that the quantity γ' is negligible, for example intransport systems or between RF gaps.
The Lorentz force law containing all the electromagnetic fields gives the forces F on a beam particle.
We decompose the force vector F into the superposition of the applied external forces Fa and the self-forces
Fs caused by the electromagnetic fields of the beam itself. The applied forces are assumed to be linear bydesign, since it is known that nonlinear forces can degrade beam quality. Thus,(8)

Fa = −Κ (s ) ⋅ r ,where K≡(κx(s),κy(s),κz(s)) represents the action of external beamline components. For example, thefocusing force in the x-direction applied by a magnetic quadrupole lens is given by(9)

Fa , x = − qβc∂B y (0)∂xx,where By is the y-directed magnetic field. Thus, κx=qβc[∂By/∂x].

7

LA-UR-02-4979
The self-forces Fs are the electromagnetic fields produced by the beam itself. We have(10)

Fs = qΕ s (r, s ) + qβc[r ′ + a z ] × B s (r, s ) ,where Es is the self electric field, Bs is the self magnetic field, and az is the unit vector in the direction of
the design trajectory (i.e., it points in the z-direction). The individual velocities around the synchronous
particle, given by r', is typically negligible compared to the collective motion of the bunch. Consequently
we drop the r'. This action leads to magnetic self-fields that are directly proportional to the electric selffields in the perpendicular direction. The reason for this situation is that the collective motion along thedesign trajectory produces the magnetic fields. (Forces from longitudinal magnetic fields will be zero by
the cross product.) We have
(11)

B⊥ =

β
c

a z × Es ,where B⊥ is the magnetic field perpendicular to az. Consequently, the perpendicular force has the form
F⊥=(1-β2)E⊥, where E⊥ is the perpendicular self electric field. The parallel self-force depends only on theparallel self electric field. Collecting these results, the forces in all three directions can be written
Fx = −κ x ( s ) x +(12)

F y = −κ y ( s ) y +

q

γ2
q

γ2

E x ( x , y , z; s ) ,E y ( x , y , z; s ) ,Fz = −κ z ( s ) z + qE z ( x, y , z; s ) .

The resulting equations of motion for a beam particle are

κx
q
γ′
x′ +
x = 3 2 2 Ex ,
2
2
γ
γmc β
γ mc βκyq
γ′y ′′ + y ′ +y = 3 2 2 Ey ,
2
2
γ
γmc β
γ mc βκzqγ′
γ′z ′′ + z ′ +
z=− +
Ez .
2
2
γ γmc 2 β 2γγmc βx ′′ +(13)

Now we simplify the equations with some standard definitions. First, let
( 14)

kα2 ≡

κα
γmc 2 β 2where α ∈ {x, y, z} ,which are the linear proportionality constants for the external focusing. We further reduce the above
equation set by introducing a standard parameter in beam physics, the beam perveance K. This is ameasure of the space-charge effect and we define it as
(15)

K≡

qN
1
q
,
2πε 0 γ 3 β 2 mc 2

where ε0 is the permittivity of free space and N is the total number of particles in the beam bunch. Thisdefinition is somewhat different than that of other authors [21][22]; it is more useful for the present
discussion. Unfortunately, there is no standard definition for the bunched beam perveance, as there is forcontinuous beams. The resulting equations of motion are

8

LA-UR-02-4979
2πε 0
γ′x ′ + k x2 x = K
Ex ,
qN
γ
2πε 0
γ′
y ′′ + y ′ + k y2 y = K
Ey ,
qN
γ
2πε 0
γ′
γ′
z ′′ + z ′ + k z2 z = − + γ 2 K
Ez .
qN
γ
γ
x ′′ +(16)

Note that qN is the total charge of the beam bunch, usually denoted by Q. It can also be expressed in termsof the bunched beam current I as
(17)

Q=I/ f ,where f is the drive frequency of the RF cavity.

3. STATISTICS BACKGROUND
3.1. The Density Function
Assume that the particle beam can be described by a density distribution function f of the particles'
phase space coordinates (x, x', y, y', z, z') and position along the design trajectory s. That is, the entire
ensemble of beam particles is represented by the function f. The function f describes the distribution of the
beam's mass (or charge) in phase space at each location s along the beamline. It can also be interpreted as aprobably density function characterizing the probability that any particle occupies a particular region ofphase space at location s. Thus, the probably that a particle lies in the infinitesimal phase space volume
dxdx'dydy'dzdz' centered at (x, x', y, y', z, z') and located at position s along the beamline is(18)

f ( x, x ' , y , y ' , z , z ' ; s ) dxdx ' dydy ' dzdz ' .

Typically we assume that the particle ensemble is populous enough to be represented accurately by acontinuous f. However, it is always possible to formulate a discrete ensemble by assuming that f is asummation of displaced Dirac delta functions. Finally, note that the function f(x,x',y,y',z,z';s) contains all
the information of the beam but requires an enormous amount of storage space to represent on a computer,
it is a function of seven independent variables.
The evolution of the function f is governed by the Vlasov equation, which incorporates the equationsof motion and the Lorentz force equation. It is a partial differential equation containing partial derivativesof all seven independent variables and is typically intractable to solve in general. We shall avoid the
Vlasov equation by propagating only a small subset of the distribution's moments. For this we require only
the equations of motion and Liouville's theorem.
3.2. Liouville's Theorem
Liouville's theorem states that the total derivate of f with respect to s along particle trajectories is zero,or formally(19)df [ x( s ), x' ( s ), y ( s ), y ' ( s ), z( s ), z' ( s ); s ] = 0 .
dsNote that this is a convective derivative that follows the particle trajectory. The fully expanded derivative
is
(20)

∂f
∂f
∂f
∂f
∂f
∂f
∂f
=0,
z ′′ +z′ +y ′′ +y′ +x ′′ +
x′ +
∂y ′
∂z
∂z ′
∂s
∂y
∂x ′
∂x

where x''=dx'/ds, etc.
This is essentially a statement of the conservation of mass. The practical considerations where thetheorem holds true are for collisionless systems where there are ample enough particles such that the self-

9

LA-UR-02-4979
fields are represented accurately by smooth functions. That is, there is no "graininess" of individual
particles; the distribution behaves much like a fluid. Otherwise, Liouville's theorem is true only for
distribution functions on 6N dimensional phase space, where N is the number of particles.
We mention in passing that by substituting the Lorentz force laws for x'', y'', and z'' in the above weobtain the Vlasov equation. Thus, the trajectories {z(s)} where the above holds actually represent thecharacteristic lines of the Vlasov equation.
3.3. Moments of the Distribution
Moments of the distribution are function averages on phase space, weighted with respect to the
distribution function. For example, let g(x, x', y, y', z, z') be some arbitrary function on phase space, thenthe moment of g, denoted 〈g〉, is given by
1
N
1
=
N

g ≡

(21)

∫ g ( x, x' , y, y' , z, z' ) f ( x, x' , y, y' , z, z' ; s) dxdx' dydy' dzdz' ,∫ g (z) f (z) d z ,
6

where
(22)

N=

∫ f (z; s) d z ,6is the total number of particles in the ensemble and where the integrations are taken over all of phase space.
(In mathematical parlance, we treat f as a measure on phase space.) Note that with this definition 〈g〉 is stilla function of s because f is a function of s. We shall be concerned foremost with the moments of phasespace monomials, that is, moments of the form 〈x〉, 〈x'〉, 〈x2〉, 〈xx'〉, 〈x'2〉 and their counterparts for the otherphase space coordinates. Theses moments represent the evolution of the beam's statistics, in particular thecenter of mass (average position and velocity) and the rms envelopes of the beam.
Notice also that in the above definition 〈⋅〉 is normalized. It is normalized by the total number ofparticles, a factor 1/N, so that the moment of 1 is 1, that is 〈1〉=1.
3.4. Moment Equations
Liouville's theorem enables us to formulate evolution equations for the moments of the beam.
Specifically it allows the moment operator 〈⋅〉 and the differentiation operator d⋅/ds to commute. For
example, if g(z;s) is a function of the phase space coordinates and s, then the derivative of 〈g〉 with respect
to s is given as

(23)

g

′

d
d 1
g =
ds
ds N
= g′
≡

∫ g (z; s) f (z; s) d z = N ∫ [g ′(z; s) f (z; s) + g (z; s) f ′(z; s)]d z ,61

6

Thus, we are able to move differentiation with respect to s within the moment operator and vice-versa.
Now we may average the equations of motion with respect to the beam ensemble then employ
Liouville's theorem to commute the differentiations with respect to s. In the x-plane, we have for the first
momentsx(24)′= x′ ,

2πε 0
γ′′x ′ = x ′′ = −k x2 x −x′ + KEx .
qNγThe quantity 〈Ex〉 is zero for symmetric charge distributions by pair-wise cancellation. Since we assume
that the beam has an ellipsoidally symmetric charge distribution in the sequel, we neglect this term. We
then have the following equations for the first-order moments:

10

LA-UR-02-4979

x

(25)

′

= x′ ,yγ′′
x ′ = −k x2 x −x′γ′= y′ ,zγ′′y ′ = −k y2 y −y′γ′= z′ ,γ′′z ′ = −k z2 z − [ z ′ + 1]
γ

Notice that the above equations are independent of the self-fields and represent the equations of motion
for the beam's center of mass; they form a complete set. Thus, these equations may be propagatedindependently. To simplify the analysis in what follows, we often assume a centered beam to makecomputations easier. We can then displace this centered beam according to the above equations describing
the centroid motion.
Now consider the second-order moments of the form 〈x2〉, 〈xx'〉, 〈x'2〉. Proceeding as before, theequations for the x-plane moments are

x2
(26)

′

= ( x 2 ) ′ = 2 xx '

2πε 0
γ′
′
xx ' = x ′ 2 + xx ′′ = x ′ 2 −xx ′ − k x2 x 2 + K
xE x
γ
qN
′
2πε 0
γ′ 2
x ′ 2 = 2 x ′x ′′ = −2k x2 xx ′ − 2
x′ + 2K
x ′E xγ
qN

There exist similar equations for the y moments. The z second moment equations are
z2
(27)

′

= 2 zz ′

2πε 0
γ′
γ′
′
zz ′ = z ′ 2 −zz ′ − k z2 z 2 −z + γ 2K
zE z
qN
γ
γ
′
2πε 0
γ′ 2
γ′
z ′ 2 = −2k z2 zz ′ − 2z′ − 2z ′ + 2γ 2 K
z ′E zqNγγ

In general we may also have cross moments of the form 〈xy〉, 〈x'y〉, 〈xy'〉, and 〈x'y'〉. The cross moments
would typically result from bending magnets or misaligned beamline elements.
From this point two differing approaches to simulation are typically employed: 1) we use the secondorder moment equations to develop a set of ordinary differential equations describing the rms envelopes ofthe beam, this technique is covered in Section 6. 2) We propagate the entire set of second-order momentsin the form of a symmetric correlation matrix, this technique is covered in Section 7. The later approach ismore general than the former, since it can treat rotated ellipsoids and external coupling between phaseplanes.

4. BEAM PHYSICSIn this section we outline quantities and concepts that are particular to beam physics. Specifically, we
cover rms envelopes, rms emittance, and the rms phase space ellipse.
We also consider a linear
approximation to the self-fields that is particularly important for beam physics calculations.
4.1. The RMS Envelopes
Here we introduce the notion of rms beam envelopes, which is fundamental in much of the beamphysics literature. The rms envelopes of a beam represent the boundary of the beam in a collective, or
statistical sense. As we see in the sequel, for ellipsoidally symmetric beams this statistical behavior isalmost independent of the actual profile of the beam distribution. That is, the rms envelopes of many
beams behave the same regardless of the actual distribution.

First, let us consider the mean values of the particle distribution. We use x , y , z to denote theses firstorder spatial moments, that is11

LA-UR-02-4979(28)

x= x ,

y= y ,

z= z .

These values are the positions of the center of mass in the x, y, and z directions, respectively. That is, they
are the coordinates of the beam centroid.
Now we consider the rms envelopes of the beam.
definitions:

To do so we first introduce the following

1/ 2~
x ≡ x2
,

(29)

1/ 2
~
y ≡ y2,
1/ 2~z ≡ z2.

In the literature these quantities are commonly referred to as the rms envelopes of the beam. It is importantto point out that this is only true for a centered beam. The true rms envelopes of a beam are actually thestandard deviations of the density distribution f. Denoting the standard deviations in the x, y, and z
directions as σx, σy, σz, respectively, they are defined(30)

σ x ≡ (x − x)2

1/ 2

σ y ≡ ( y − y)2

1/ 2

σ z ≡ (z − z)2

1/ 2

[
] ,
=[ y −y ] ,=[z −z ] ,1/ 2= x2 − x 2
2

2

2

2

1/ 2

1/ 2

or
(31)

[

σx = ~
x2 − x2

] ,
1/ 2

[

σ y = ~y 2 − y 2

] ,

[

1/ 2

σ z = ~z 2 − z 2

]

1/ 2

.

Note that when the beam is centered, that is when x = y = z = 0 , the rms envelopes σx, σy, σz are equal to
x, ~
y, ~
z . This is the situation most often encountered in the literature. Typically a centered
the quantities ~beam is assumed to simplify computations, particularly for the moments 〈xEx〉, 〈yEy〉, and 〈zEz〉. A common
approximation in simulation is the assumption that the beam is centered when computing the evolution of
x, ~
y, ~
z ; the offsets x , y , z are computed independently to form the complete beam
the second moments ~state.
4.2. RMS Emittance

One particularly important quantity in beam physics is the rms emittance, usually denoted ε~ . This is afigure of merit indicating the area in each phase plane that the rms beam envelope occupies. The rms
emittance is a function of the second moments, for each phase plane it is defined[] ,ε~ ≡ [ y y ′ − yy ′ ] ,ε~ ≡ [ z z ′ − zz ′ ] .

ε~x ≡ x 2 x ′ 2 − xx ′
(32)

2

2

2

2

2 1/ 2

2 1/ 2y

2 1/ 2z

It is known that whenever all forces acting on the beam are linear and there is no acceleration, the rms
emittance is an invariant of the motion [7][10][12]. An increasing emittance is usually indicative of a loss
in beam quality. Thus, the rms emittance of a beam typically increases whenever unwanted nonlinear
forces are encountered.
When the beam is accelerated the transverse rms emittances always decrease. This condition is simply
an artifact of the definitions of x' and y' and does not imply any decrease in random transverse kinetic

12

LA-UR-02-4979
energy (temperature). Since x'=px/p, it decreases with increasing longitudinal momentum p, which in turndecreases the rms emittance. To see this consider the x phase plane and assume that the self-fields can be
represented as the linear function Ex=ax, where a is a real number. Now simply differentiate the definitionof rms emittance′d ~2
ε x = x 2 x ′2 + x 2
dsγ′= −2 ε~ 2 ,

(33)

γ

x ′2

′

′
− 2 xx ′ xx ′ ,xwhere the second line is obtained by substituting the derivatives from Eqs. (26). The solution to thisequation isε~x2 ( s ) =(34)

ε~x2,0γ 2 ( s),where ε~x ,0 is the initial value. Thus, we can see outright that when the beam energy increases thetransverse rms emittances decrease. Applying the same procedures to the longitudinal plane we get adifferent equation

γ′
γ′
d ~2
ε z = −2 ε~z2 + 2z zz ′ .
γ
γ
ds

(35)

Thus, if the RF phase is misaligned during acceleration, causing a nonzero 〈z〉, it is actually possible toexperience emittance growth in the longitudinal plane.
To alleviate the situation of decreasing transverse emittances with acceleration we may alternately use
the normalized rms emittances, typically denoted ε~n . These emittances are defined as [21]ε~n , x ≡ βγε~x ,(36)

ε~n , y ≡ βγε~y ,ε~n , z ≡ βγε~z .

Performing a similar analysis as with the rms emittances using similar assumptions, we find that thenormalized rms transverse emittances behave as
ε~ ( s ) = β ( s ) ε~ ,
(37)
n, x

n , x ,0where ε~n , x ,0 is an arbitrary constant.

Thus, the transverse normalized rms emittances actually increase withincreasing beam energy. However, once the beam velocity approaches c these emittances level off withincreasing beam energy.
In our phase space coordinates the units of rms emittance are meter-radians. However, recalling that
there are several alternate coordinate systems for the z phase plane, it is common to find emittance values
for this plane given in several different units. When the z phase plane coordinates are (z,∆p/p) the units are
meter-radians. When the coordinates are (∆φ,∆W) the units are degrees-electron volts. Finally, when thecoordinates are (∆t, ∆W) the units are seconds-electron volts. To convert between theses units to meterradians we substitute the coordinate conversion formulas (5) and (6) into the above definition for ε~z . The
results are
ε~ ( m − rad ) = C ε~ ( m - rad ) ,
z

(38)

∆p

z

= C deg ε~z (deg − eV ) ,= C ε~ (sec − eV ) ,seczwhere the conversion factors C∆p, Cdeg and Csec are defined13

LA-UR-02-4979
C ∆p ≡(39)

C deg ≡

1γ2

( no units) ,

βλ
λ
1
1
1
1
=
360 γ (γ + 1) W 360 γ 2 γ 2 − 1 E R

Csec ≡ βc⎛ m − rad ⎞⎜⎜⎟⎟ ,⎝ deg − eV ⎠γc ⎛ m − rad ⎞
1
1
=
⎜
⎟,
γ +1 Wγ 2 − 1 E R ⎝ sec − eV ⎠where ER=mc2/q is the rest energy of the beam particles in electron volts. Note that the beam energy W
must also be given in electron-volts.
4.3. Courant-Snyder Parameters for the RMS Ellipse
When modeling beam distributions, the projections of the distributions onto each phase plane arerepresented ideally by an ellipse. Here we introduce the Courant-Snyder, or Twiss, parameters often foundin the literature. These parameters describe the projections of the beam distribution onto the phase spaceplanes. In particular, we shall describe the ellipse relating to the second-order moments of the beam (therms moments). This ellipse is congruent to the phase space ellipse of the equivalent uniform beam. To
make explicit calculations we consider the x phase plane with the coordinates (x,x'). There are analogous
results for the y and z phase planes.

Consider an ellipse centered in x phase space as shown in
Figure 2. The general equation for such an ellipse can be
expressed(40)x’

γ xε~x−αx

γ x x 2 + 2α x xx ′ + β x x ′ 2 = ε~x ,where αx, βx, γx are known as the Courant-Snyder parameters forthe x phase plane. These parameters are not independent; they
are related through the equation(41)

β xε~xε~xβxxβ xγ x − α x 2 = 1 .

This relation enforces the area enclosed by the ellipse to be πε~x .
It is our desire to relate this ellipse to the rms quantities of thebeam distribution.

Figure 2: rms phase space ellipse

We can solve Eq. (40) for particular values of x and x' in the phase space as shown in Figure 2 (forexample, see [21] and [22]). Consider the projections of this ellipse onto the x and x' axes. We want the
maximum extent of these projections to be equal to the rms quantities for x and x', respectively. Thiscondition leads to the definitions

β x ε~x ≡ x 2
(42)

1/ 2

,

α x ε~x ≡ − xx ′ ,γ x ε~x ≡ x ′ 2

1/ 2

,

where the second equation follows from Eq. (41) and the definition of rms emittance (32).
If the coordinates for the longitudinal phase plane are given as (∆φ,∆W), the conversion of the CourantSnyder parameters to the (z,z') units is given by the following:

14

LA-UR-02-4979
α z ( no units) = α z ( no units),
βλ γ (γ + 1) W
β z ( m/rad) =β z (deg/ eV ) ,

(43)

360
360

1
1
1
1
γ z ( rad/m) =γ z (eV/ deg),
βλ γ (γ + 1) W

where, due to an unfortunate choice of notation, we have β's and γ's representing both relativistic
parameters and Courant-Snyder parameters. The grouping of terms was chosen to make the context clear.
When the longitudinal coordinates are given as (z,∆p/p), the conversion of the Courant-Snyder parameters
to (z,z') units is

α z ( no units) = α z ( no units),
β z ( m/rad) = γ 2 β z ( m/rad) ,(44)

γ z ( rad/m ) =1γ2γ z ( rad/m),

When the longitudinal coordinates are given as (∆t,∆W), the conversion of the Courant-Snyder parameters
to the (z,z') units is given by

α z ( no units) = α z ( no units),
γ +1
β z ( m/rad ) = βcWβ z (sec/ eV ) ,γ
1 γ 1γ z ( rad/m ) =γ z (sec deg).
βc γ + 1 W(45)

4.4. The Mean Value Vector and Correlation Matrix

Here we introduce the mean-value vector z and symmetric 6×6 correlation matrix σ. These quantities
contain the state of the beam statistics up to second order. Maintaining consistency with previous notation,
let the moments 〈x’〉, 〈y’〉, and 〈z’〉 be denoted(46)

x′ ≡ x′y′ ≡ y′andz′ ≡ z′ ,respectively. We accordingly denote the vector of first-order moments as(47)z≡ z ,
= (x

x′

y

y′ z

z ′)T ,

where the superscript T indicates transposition. Thus, z is the column vector of mean values for the phasespace coordinates. Note that z can be a function of the path length parameter s.
When propagating the second-order moment down the beam line it is convenient to work with a matrix
whose elements are the second-order moments, denoted σ. Taking the moment of zzT, which is asymmetric matrix, forms σ. For example, considering only the x phase plane whose coordinates are
represented by the vector( 48)x( s ) ≡ (x ( s ) x ′( s ) )T .

The matrix we are describing is given as 〈xxT〉, a sub-matrix of σ, which we denote by σxx. That is(49)

σ xx ≡ xx

T

⎛ x2
⎜
=⎜
⎜ xx '⎝

xx ' ⎞⎟,
x ′ 2 ⎟⎟⎠

15

LA-UR-02-4979
The matrix σxx is known as the correlation matrix of the x phase plane. We have one such 2×2 matrix foreach phase plane, plus we have matrices containing the cross-correlations between phase planes such as〈xyT〉, 〈xzT〉, and 〈yzT〉. The full second-order moment matrix σ has the block matrix form(50)

σ ≡ zzT⎛ σ xx⎜= ⎜ σ yx
⎜σ
⎝ zx

σ xy
σ yy
σ zy

σ xz ⎞⎟σ yz ⎟ ,σ zz ⎟⎠

where each σαβ is a 2×2 symmetric matrix. By symmetry we must have σαβ=σβα. When the motion is notcoupled between phase planes planes, σ has the simplified block-diagonal form(51)

⎛ σ xx
⎜
σ=⎜ 0
⎜ 0
⎝

0
σ yy
0

0 ⎞
⎟
0 ⎟,
σ zz ⎟⎠

where 0 is the 2×2 zero matrix. Thus, coupling between phase planes is manifested by nonzero values forthe sub-matrices σxy, σxz, or σyz.
The vector z and matrix σ represent all the moments up to order two at each position s along the beamline. Thus, the pair ( z, σ) contains the complete state of the beam's statistical behavior up to order two.
This is a convenient representation for the beam statistics when doing simulation.
The correlation matrix can be expressed in terms of the Courant-Snyder parameters. Referring to Eqs.
(42), the σ matrix for the x phase plane can be written(52)⎛ β x ε~x~
⎝ − α xε x

σ xx = ⎜⎜− α x ε~x ⎞ 1 ⎛ γ x⎟=⎜
γ x ε~x ⎟⎠ ε~x ⎜⎝ α x

−1

αx ⎞
⎟ .
β x ⎟⎠

There are, of course, analogous representations for the y and z phase planes. Using the above fact we canwrite Eq. (40) for the rms phase space ellipse compactly in matrix-vector notation as(53)

x T σ −xx1x = 1 .

4.5. Linearization of the Self-Fields
It has been found that the evolution of the beam's second moments is determined primarily by thelinear part of the forces [19]. The nonlinear parts are usually associated with emittance growth of the beam(see Section 4.1). In this analysis we consider only linear forces. We have already assumed that theexternal forces are linear according to machine design. The self-fields will, in general, have nonlinear
components and, thus, the self-forces also will. Here we describe a linear approximation to the self-fields
Ex, Ey, Ez which is appropriate for our analysis.

Consider the self-field in the x-direction, Ex. We begin by assuming an expansion of the form
(54)

E x ( x ) ≈ a 0 + a1 x .

where a0 and a1 are real numbers independent of x. Terms involving the other phase space coordinates (ofthe form b1y, c1z) are zero with the assumption of symmetry (see Section 5.1). One method of determining
the coefficients a0 and a1 is by a minimum variance, or least squares, estimation technique [6][15]. Here
we seek to minimize the norm of the error in the approximation, which is ||Ex - a0 - a1x||, Ex being the actual
self-field. In order that the expansion be most accurate in areas of higher beam concentration, we choose
the norm ||⋅|| to be(55)∫

g ≡ ⎡ g 2 ( z) f ( z; s ) d 6 z⎤
⎢⎣
⎥⎦

1/ 2

= g2

1/ 2

.

16

LA-UR-02-4979
It is easily checked that this is a valid norm for functions of the phase space coordinates. Using the leastsquares fitting technique, we are left with the following Gram system for the coefficients:
(56)

⎛ 1
⎜
⎜ x
⎝

x ⎞⎛ a 0 ⎞ ⎛ E x ⎞⎟⎜ ⎟ = ⎜⎟.
x 2 ⎟⎜⎝ a1 ⎟⎠ ⎜⎝ xE x ⎟⎠⎠The solution to this system is(57)

a0 = −x

σ x2

xE x ,a1 =1xE x ,

σ x2

however, this can be reduced with some practical considerations.
Since we only use this expansion for space charge effects, we may as well assume that the beam iscentered. This action is equivalent to a simple coordinate translation, which does not affect the shape of theelectric field and, therefore, the moment 〈xEx〉. When the beam is centered (i.e., when x = y = z = 0 andσ =~x ) the values of the expansion coefficients arex(58)

a0 = 0

and

a1 =

xE x
x2

.

Thus, in this case(59)

E x ≈ a1 x .

We have analogous expressions for Ey and Ez so that all the field components have the following linear
approximations:
(60)

Ex ≈

xE x
x

2

x,

Ey ≈

yE y
y

2

y,

Ez ≈

zE z
z2

z.

Note once again that these linear approximations to the true fields are most accurate in the regions ofhighest beam density. They are most useful when considering the collective behavior, or rms behavior, ofthe beam.

4.6. Tune Depression
To demonstrate the usefulness of the material in the previous section we digress somewhat to show
how it is related to a common parameter in beam physics known as the tune depression η. Tune depression
is a parameter indicating the relative effect of space charge compared to the effects of transverse focusing
in a beam channel. To illustrate, return to the original equations of motion (16) for an individual particle.
Employing the linear approximation to the self-field given by Eq. (59), the approximate equation of motionin the x-plane is(61)

x ′′ +

2πε 0 ⎞⎛
γ′x ′ + ⎜⎜ k x2 − K
a1 ⎟⎟ x = 0 .
qN
γ
⎝
⎠

The square root of the quantity in parenthesis is known as the tune-depressed phase advance, or phaseadvance with space charge and is typically denoted k in the literature. The quantity kx in the above is called
the phase advance without space charge and is usually denoted k0. We see that an immediate consequenceof space charge is the reduction of external focusing. Consequently, the frequency of the particles' betatronoscillations, or betatron tune, decreases with increasing space charge forces.
Tune depression η is defined(62)η ≡ k / k0 .

17

LA-UR-02-4979
It is the factor by which the betatron tune frequency is decreased, having a maximum value of 1 (no tune
depression) and a minimum value approaching 0 (complete tune depression). Using the above expansion
we find the tune depression (in the x plane) to be

(63)

⎡
⎤
K 2πε 0 xE x ⎥η = ⎢1 − 2⎢ k x qN
x2 ⎥
⎣
⎦

1/ 2

⎡
xE x ⎤⎥q= ⎢1 − 2⎢ γ κ x x2 ⎥
⎣
⎦

1/ 2

.

From the above relation we see that the tune depression is a function of the self-fields, the beam energy,
and the external fields.

5. ELLIPSOIDAL BEAMS AND THE EQUIVALENT UNIFORM BEAM
Here we define and describe ellipsoidal beams. For these beams it is possible to compute the field
moments 〈xEx〉, 〈yEy〉, 〈zEz〉 analytically in terms of elliptical integrals. The results of these calculations
lead to the concept of the equivalent uniform beam. This notion says that we may model any (ellipsoidal)
beam, at least approximately, by a uniform density beam with the same second moments.

5.1. Ellipsoidally Symmetric Charge
Distributions
We restrict our attention to ensemble distributions
that are ellipsoidally symmetric in configuration space,that is, in (x,y,z) space. For simplicity we assume an
upright, centered ellipsoid at the origin. For an arbitrarilyoriented ellipsoid we can always apply a change of
coordinates to achieve this condition. In this situation the
charge density ρ(x,y,z) has the form
(64)

⎛ x2

y2

⎝a

b

ρ ( x, y , z ) = q F ⎜⎜

+
2

+
2

z2 ⎞
⎟,
c 2 ⎟⎠

z
c

a

b

y

x

Figure 3: reference ellipsoid for ellipsoidally symmetriccharge density

where q is the unit charge, and a, b, c represent the semiaxes of the reference ellipsoid in the x, y, z directions respectively. The function F(⋅) represents the profileof the distribution and is related to the density function f by
(65)

F ( x , y , z; s ) ≡∫∫∫ f ( x, x ′, y, y ′, z, z ′; s)dx ′dy ′dz ′ ,that is, we have integrated out the momentum dependence. This situation is depicted in Figure 3.
We compute the spatial moments of this distribution in order to find a relation between the moments,
the total number of particles, and the semi-axes a, b, c. For example, the total number of particles N is
(66)

N =

∫∫∫

⎛ x2 y2 z2 ⎞F ⎜⎜ 2 + 2 + 2 ⎟⎟dxdydz .
b
c ⎠
⎝a

This integral may be computed by using the change of coordinates(67)

x = ar sin θ cos φ ,y = br sin θ sin φ ,z = cr cosθ ,dxdydz = abc r 2 sin θdrdθdφ ,with the result

18

LA-UR-02-4979(68)

N = 4πabc∞∫ r F (r ) dr ,220

= 2πabc F1 / 2 .

where Fn is the nth moment of F, defined as
(69)

Fn ≡

∞

∫ s F ( s) ds .
n0Note that n is not restricted to integer values. Using a similar approach we can compute the following
moments of the ellipsoidal distribution:
1 = 1,
x = y = z = 0,

(70)

x2 =

a 2 F3 / 2,
3 F1 / 2y2 =b 2 F3 / 2,
3 F1 / 2z2 =c 2 F3 / 2.
3 F1 / 2

These equations provide us expressions of the semi-axes a, b, c in terms of the second spatial moments 〈x2〉,〈y2〉, 〈z2〉, respectively.

5.2. Computation of the Field Moments
It is possible to compute the field moments explicitly (in terms of elliptic integrals) for beams having
ellipsoidal symmetry in configuration space. We start from the following formula for the self-electric
potential of such a bunch [9][13]:(71)

φ ( x, y , z ) =qabc ∞ ∞F (s)ds dt ,
4ε 0 0 T ( x , y , z;t ) t + a 2 1 / 2 t + b 2 1 / 2 t + c 2 1 / 2

∫∫

(

) (

) (

)

where the limit of integration T(x,y,z;t) is definedx2

T ( x , y , z; t ) ≡

+
2

y2

+
2

z2

.
t+a
t+b
t + c2
The formula assumes that the bunch is centered at the (x,y,z) origin. It is found by inverting Laplace's
equation for φ in ellipsoidal coordinates (the variable t is actually an ellipsoidal coordinate). Using Eq. (71)the expressions for the field moments are computed to be (see Appendix A)

(72)

(73)

xE x = Γ( F )qπ a 2 b 2 c 2 2
a RD (c 2 , b 2 , a 2 ) ,6ε 0
NyE y = Γ( F )qπ a b c 2
b R D (a 2 , c 2 , b 2 ) ,6ε 0
NzE z = Γ( F )2 2 2

qπ a b c 2
c RD (a 2 , b 2 , c 2 ) ,
6ε 0
N
2 2 2

∞

∫
Γ( F ) ≡ ∫ G ( r ) dr,

G ( r ) ≡ F ( s ) ds ,
r

∞

2

2

0

where Γ carries the distribution information (it is a function of the profile F), and RD is the Carlson elliptic
integral of the second kind [3][4].
Carlson’s definitions for elliptic integrals are much more convenient in this situation. Moreover, theylend themselves nicely to numeric computation [17]. Definitions for these integrals are as follows:

19

LA-UR-02-4979(74)

R F ( x, y, z ) ≡

1 ∞
dt
,
1
/
2
0
2
(t + x ) (t + y )1 / 2 (t + z )1/ 2

R D ( x, y, z ) ≡

3 ∞
dt
,
2 0 (t + x )1 / 2 (t + y )1 / 2 (t + z )3 / 2

R J ( x, y , z, p ) ≡

∫
∫

3 ∞
dt
.
1
/
2
1
2 0 (t + x ) (t + y ) / 2 (t + z )1 / 2 (t + p )∫Note that only RF is symmetric in its arguments. It is possible to express all the conventional elliptic
integrals (E, F, K, etc.) in terms of RF, RD and RJ.

5.3. The Equivalent Uniform Beam
We would prefer to express Eqs. (73) in a form independent of the parameters a, b, c. Rather we could
use the variables 〈x2〉, 〈y2〉, 〈z2〉 by employing the relations (70). First note that(75)

R D ( x , y , z ) = r 3 / 2 R D ( rx , ry , rz )for any positive real number r. Thus, choosing r=1/c2 we have(76)

RD (a 2 , b 2 , c 2 ) =a 2 b2[, ,1] .
R
D
c2 c2
c3
1

Using this property of RD we can express the self-moments asxE x =(77)yE y =zE z =⎡< y2 > < z2 > ⎤
Q
1
R
,
,1⎥ ,
⎢
D
2
2
3 24πε 0 < x 2 >1 / 2⎣< x > < x > ⎦Λ( F )⎡ < x2 > < z2 > ⎤
Q
1
R
,
,1⎥ ,
⎢
D
2
2
3 24πε 0 < y 2 >1 / 2⎣< y > < y > ⎦Λ( F )⎡< x2 > < y2 > ⎤
Q
1
R
,
,1⎥ ,
D⎢
2
2
3 24πε 0 < z 2 >1 / 2⎣< z > < z > ⎦Λ( F )

where the functional Λ(F) is defined as

(78)

1/ 2

⎡ ∞ 3/ 2⎤r F ( r )dr ⎥
1/ 2
⎢
∞
0
(F )
⎦ .
Λ ( F ) ≡ Γ( F ) 3 / 2 5 / 2 = G 2 ( r 2 )dr ⎣5/ 2
0∞( F1 / 2 )
⎡
⎤
1/ 2
r
F
(
r
)
dr
⎢0
⎥
⎣
⎦

∫

∫
∫

The functional Λ accounts for the effects of the particular beam distribution profile on the beam dynamics.

20

LA-UR-02-4979
Table 1: Ellipsoidally symmetric profiles F and the corresponding Λ(F). The quantity C represents an arbitrary constant.

Distribution

Λ(F)

F(s)

Uniform⎧C⎨⎩0

Parabolic (Waterbag)⎧ C (1 − s )
⎨
⎩0

s ≤1
s >1

Hollow−Cs e
Gaussian

6 3≈ 0.9295
5 5

s ≤1
s >1

10 3
≈ 0.9352
7 7

s

3 5
≈ 0.9462
4 π

2σ 2

−

s

C e 2σ32π≈ 0.9772

Sacherer discovered that the functional Λ is nearly constant with respect to the distribution F [19].
Table 1 lists the values of Λ for several common distributions. We see there that its value varies only by
about 5% for the distributions listed. The implication is that the second-order spatial moments (i.e., the rms
envelopes) of any bunched beam behave approximately the same, irrespective of the particular beamprofile. In other words, the beam dynamics are only loosely coupled to the actual form of the distribution,
so long as it is ellipsoidal. Consequently, we are justified in modeling any laboratory beam with any otherbeam having the same second moments 〈x2〉, 〈y2〉 and 〈z2〉. This notion leads to the concept of theequivalent uniform beam. This bean is the uniform density beam having the same second moment as theactual laboratory beam. Since the uniform beam has well-defined boundaries, this distribution is typically
preferred. Finally, we note in passing that in the case of continuous elliptic beams, the dynamics arecompletely independent of the beam profile.
Recognize that the semi-axes of the reference ellipsoid a, b, c, are actually the envelopes of theuniform beam. This being so, we can use Eqs. (70) to relate the envelopes of the uniform beam and thesecond spatial moments. It is common in the literature to denote the envelopes of the equivalent uniform
beam by X, Y, Z corresponding to the semi-axes a, b, c, respectively. Thus, for the equivalent uniform
ellipsoid we have(79)[

X = 5 x2

] = 5~x[

1/ 2

Y = 5 y2

] = 5 ~y1/ 2[

Z = 5 z2

] = 5~z .
1/ 2

Likewise, it is common to define the effective emittance ε of the equivalent uniform beam in terms of therms emittance according to the following:(80)
ε ≡ 5ε~ε ≡ 5ε~ε ≡ 5ε~ .
x

x

y

y

z

z

The effective emittances represent the areas in the phase planes occupied by the equivalent uniform beam.
We can also determine the second-order moments of the equivalent uniform ellipsoid using theirdefinitions (79) and the definition of rms emittance (32). For the x-plane we find
1 2
X ,
5
1
xx ′ = XX ′ ,
5

x2 =

(81)

x ′2 =

1 2 5ε~x2
X′ + 2 ,5X

where X'=dX/ds is the rate of change of the equivalent beam envelope with respect to s. There areanalogous relations for the other phase planes.

21

LA-UR-02-4979
The envelope X of the equivalent uniform ellipsoid can be expressed in terms of the CourantSnyder parameters through inspection of the extremal points of the trace-space ellipse (for example, see
Reiser [21]). The result is(82)

X = 5β x ε~x ,5ε~X ′ = −α x.

βxThere are similar expressions for the y and z phase planes.

6. BUNCHED BEAM ENVELOPE EQUATIONS
The final result of this section is a set of coupled ordinary differential equations describing theevolution of the equivalent uniform beam, given by Eqs. (86). Although nonlinear, these equations are a
relatively simple set of ordinary differential equations that may be integrated numerically using standard
techniques. The shortcoming of this description is that, aside from space charge, it does not account forexternal coupling between phase planes. That is, coupling such as that from misaligned quadrupoles, skew
elements, etc., cannot be modeled. This drawback is circumvented in the follow section on transfer matrix
methods by considering all the second-order moments.

6.1. Equations for Centroid Motion
The equations for the first-order moments are simply the equations of motion for the centroid of thebeam, which behaves as a single particle. The evolution equations for the average values are found withreference to Eqs. (25). By differentiating the first equation of each set in Eqs. (25) then substituting into
the second we findγ′
x ′ + k x2 x = 0 ,γγ′y ′′ + y ′ + k y2 y = 0,γγ′γ′z ′′ + z ′ + k z2 z = − ,γγx ′′ +(83)

which are a closed set involving only x , y , and z . They are damped harmonic equations and may besolved independently by standard techniques. Recall that the kα2 are functions of s and that it is possible forkα2 to be negative, that is, κα is negative.

6.2. Bunched Beam RMS Envelope Equations
Here we develop a set of coupled ordinary differential equations that describe the evolution of the rms
x, ~
y, ~
z . These equations describe the behavior of theenvelopes of the beam, specifically the quantities ~beam extent in a statistical sense. They can then be used to derive the equations for the equivalent uniform
beam, as is done in the following subsection.
From the second-order moment equations (26) we can derive equations involving only the quantities
~
x, ~
y, ~
z . Recall that ~x ≡〈x2〉½. If we differentiate ~x twice with respect to s then use the relations (26) andx . Proceeding in an
the definition of rms emittance (32) we find a second-order differential equation for ~analogous manner for the y and z planes we are left with the following set of ordinary differential
equations:22

LA-UR-02-4979(84)

2πε 0 xE x
ε~x2
γ′
~
x ′′ + ~x ′ + k x2 ~
x−K
−
= 0,
~
~
qN
x
γ
x3
ε~y2
2πε 0 yE yγ′~y ′′ + ~y ′ + k y2 ~
y−K
−
= 0,
~
~
qN
y
γ
y3
2πε 0 zE z
ε~z2
γ′
γ′
~
z ′′ + ~z ′ + k z2 ~z − γ 2K
−
=− z.
~
3
~
qN
z
γ
γ
z

Note that these equations involve the moments 〈xEx〉, 〈yEy〉, 〈zEz〉. As we have seen, these moments are
x, ~
y, ~
z ; so the above set form a coupled set of ordinary differential equations.
generally functions of the ~Thus, once we determine the field moments we can solve the system using standard techniques (e.g.,
numerically). As they stand, these equations are valid for all beam distributions. In the next subsection we
consider ellipsoidally symmetric beam distributions.

6.3. Ellipsoidally Symmetric Beams and the Equivalent Uniform Beam
Here we substitute the values 〈xEx〉, 〈yEy〉, 〈zEz〉 that we calculated for beams with ellipsoidal symmetry
into the above rms envelope equations (84). The result is

(85)

⎡~
Λ( F ) K 1γ′
y2 ~z 2 ⎤ ε~x2~x ′′ + ~x ′ + k x2 ~
x−
R
,
⎢
D
~ 2 ~ 2 ,1⎥ − ~ 3 = 0 ,γ
x2
3 12 ~⎣x x ⎦ x~2⎡~
Λ( F ) K 1γ′
x2 ~z2 ⎤ εy~−y ′′ + ~y ′ + k y2 ~
y−
R
,
,
1
D ⎢ ~2 ~2 ⎥
~3 = 0,γ
y2
3 12 ~⎣y y ⎦ y⎡~Λ( F ) K 1
γ′
γ′
x2 ~
y 2 ⎤ ε~z2~z ′′ + ~z ′ + k z2 ~z −γ 2R,z.
⎢
D ~ 2 ~ 2 ,1⎥ − ~ 3 = −2~γγ
3 12 z⎣z z⎦ zThese are the rms envelope equations for bunched beams with ellipsoidal symmetry.
Now consider the equivalent uniform beam. Substituting the value of Λ for the equivalent uniform
ellipsoid and using relations (79) and (80) we have
X ′′ +

(86)

ε2
K
γ′
X ′ + k x2 ( s ) X − X R D ( Z 2 , Y 2 , X 2 ) − x3 = 0 ,γ′
2XY ′′ +
Z ′′ +ε 2yγ′KY ′ + k y2 ( s )Y − Y R D ( X 2 , Z 2 , Y 2 ) − 3 = 0 ,
γ′
2
Y

ε2
γ′
γ′
K
Z ′ + k z2 ( s ) Z − γ 2 Z R D ( X 2 , Y 2 , Z 2 ) − z3 = − 5 z ,
γ′
γ′
2
Z

where we have used relation (75) to eliminate the fractional arguments in the elliptical integrals. When
doing numerical computation, however, using the rational arguments for the elliptical integrals is probably
more stable. For example, from Eq. (75) it is probably best to compute RD(X2,Y2,Z2) as(87)() Z1 R ⎡⎢ XZ , YZ ,1⎤⎥ .

RD X 2 , Y 2 , Z 2 =

3

D

⎣

2

2

2

2

⎦

23

LA-UR-02-4979
6.4. Approximate Form of the Envelope Equations
In the literature, the bunched beam envelope equations are sometimes expressed without explicit
reference to elliptic integrals, typically by introduction of a “form factor” and/or approximations in lieu of
the integrals [11]. There, we see algebraic expressions or expressions involving elementary functions
instead of the special functions. One way to achieve this form is with the procedure given below. The
approximation presented here is accurate when the transverse plane envelopes X and Y are approximately
equal.

HL

To simplify the discussion below, we begin byimmediately introducing the form factor ξ. It isdefined byξ ( s) ≡

(88)

x s

1

dt
s ∞
2 0 (t + 1) t + s 2 3 / 2

∫

(

0.8

)

0.6

s
⎧
cos −1 s for s < 1
⎪1 −
2
1 ⎪
s
−
1
=
⎨
s
1 − s 2 ⎪1 −cosh −1 s for s > 1
2
⎪⎩
s −1

0.4
0.2
s
1

2

3

4

5

A plot of this function is shown in Figure 4. From
Figure 4: form factor ξ(s)the figure we see that ξ has a value of 1 at s=0, 1/3
at s=1, then asymptotes toward zero as s approaches
infinity. Note that ξ has the definition given by the integral expression, however, this integral may be
expressed analytically in terms of elementary functions according to the above.
Using the form factor we can approximate the elliptic integrals by the following:(89)

RD ( Z 2 , Y 2 , X 2 ) ≈⎛ Z ⎞⎤31 ⎡⎟⎟⎥ ,⎢1 − ξ ⎜⎜XZ X + Y ⎣⎝ XY ⎠⎦RD ( Z 2 , X 2 , Y 2 ) ≈⎛ Z ⎞⎤31 ⎡⎟⎟⎥ ,⎢1 − ξ ⎜⎜YZ X + Y ⎣⎝ XY ⎠⎦RD ( X 2 ,Y 2 , Z 2 ) ≈3 ⎛ Z ⎞⎟,ξ⎜
XYZ ⎜⎝ XY ⎟⎠where the approximations are most accurate when X≈Y, that is, the axisymmetric case. Substituting theapproximations for the elliptic integral into the envelope equations for the equivalent uniform ellipsoid (86)yields the approximate envelope equations for the equivalent beam. They appear as follows:
X ′′ +(90)

3K 1 − ξ ( Z / XY ) ε x2γ′
− 3 =0,X ′ + k x2 ( s ) X −2γZ(X + Y )X

Y ′′ +

2
3K 1 − ξ ( Z / XY ) ε yγ′− 3 =0,Y ′ + k y2 ( s )Y −2γZ(X +Y )Y

Z ′′ +

3K ξ ( Z / XY ) ε z2γ′γ′
− 3 =− 5 z,
Z ′ + k z2 ( s ) Z − γ 2
2
γ
γ′
XY
Z

where we note the rational argument of the form factor ξ.
For the reminder of this subsection we confirm the approximations of Eqs. (89). We begin byperforming an expansion of the x and y envelopes X, Y about their axisymmetric value R. Employing asmall perturbation parameter ε <<1 (not to be confused with beam emittance) we have24

LA-UR-02-4979(91)

X = R + εX 1 + ε 2 X 2 + K ,Y = R + εY1 + ε 2Y2 + K .

With this expansion, note that
(92)

R=

X +Y
+ O (ε )
2

and

lim

X −Y

ε →0

ε

= X 1 − Y1 .

That is, R can be interpreted as the average value of X and Y to first order. The perturbation parameterallows us to keep track of the order of accuracy in our approximations. As ε increases the beam may
become increasingly more eccentric in the transverse plane. In fact, our final results for this section areexact for axisymmetric beams.
Consider the envelope equation for the z plane. Writing out the elliptic integral explicitly in this casegives the following expression:(93)

RD ( X 2 , Y 2 , Z 2 ) =31
2 ( XY )3 / 2

∫

∞

dt

(t + X / Y )(t + Y / X ) (t + Z 2 / XY )

,

3/ 2

0

where we have pulled a factor XY from each product in the denominator, then applied a change ofvariables. The expression under the radical can be written t 2 + ( X / Y + Y / X )t + 1 with the term inparenthesis expanded as(X − Y ) 2⎛X Y ⎞
4⎜ + ⎟ = 2 + 1 2 1 ε + O(ε ) .
R⎝Y X ⎠
2

(94)
Thus,
(95)

t 2 + ( X / Y + Y / X )t + 1 = (t + 1) + O (ε 2 )and we can approximate the elliptic integral as
RD ( X 2 , Y 2 , Z 2 ) ≈

(96)
=

3 1
2 XYZ

∞

dt

0

2

XY ∫ (t + 1) (t + Z / XY )
Z

3/ 2

,

3 ⎛ Z ⎞⎟ξ⎜
XYZ ⎜⎝ XY ⎟⎠where we have identified the form factor ξ(⋅).
Now consider the transverse plane, in particular, the x plane. Results for the y plane follow in an
analogous manner. The elliptic integral can be written(97)

RD ( Z 2 , Y 2 , X 2 ) =31
2 ( XY )3 / 2∞∫ (t + Z / XY ) (t + Y / X ) (t + X / Y )
dt

1/ 2

0

1/ 2

3/ 2

,

where, again, we have pulled a factor XY from each product in the denominator, then applied a change ofvariables. Now we expand the last two terms in the denominator as(t + Y / X )1 / 2 (t + X / Y )1/ 2 (t + X / Y ) = (t + 1)(t + X / Y ) + O (ε 2 ) ,= (t + 1)(1 / Y )(tY + X ) + O (ε 2 ) ,(98)

R(t + 1)2 + O (ε ) ,Y1( X + Y )(t + 1)2 + O (ε ) .
=
2Y

=

25

LA-UR-02-4979
Thus, we can approximate the elliptic integral as
31
XZ ( X + Y )RD ( Z 2 , Y 2 , X 2 ) ≈

(99)

∞

,
XY ∫ (t + Z / XY ) (t + 1)
Z

dt

2

0

1/ 2

2

⎛ Z ⎞
3
1
⎟
η ⎜⎜
XZ ( X + Y ) ⎝ XY ⎟⎠=where the auxiliary function η(⋅) is definedη ( s) ≡
(100)

∞

s dt

0

2 1/ 2∫ (t + s ) (t + 1) ,2

= 1 − ξ ( s) .
The equality on the second line is found through integration by parts. From this second relation we can
express the elliptic integral for the transverse plane in terms of the original form factor(101)

RD ( Z 2 ,Y 2 , X 2 ) ≈⎛ Z ⎞⎤31 ⎡⎟⎟⎥ .
⎢1 − ξ ⎜⎜XZ X + Y ⎣⎝ XY ⎠⎦

Likewise, for the y plane we have(102)

RD ( Z 2 , X 2 , Y 2 ) ≈⎛ Z ⎞⎤31 ⎡⎟⎟⎥ .
⎢1 − ξ ⎜⎜YZ X + Y ⎣⎝ XY ⎠⎦

Notice that our approximations for the transverse plane elliptic integrals are only accurate to order εwhereas the longitudinal plane approximation is accurate to order ε2. This condition makes sense in thatthe transverse plane approximations are more sensitive to the eccentricity in X and Y than the longitudinalapproximation.

7. TRANSFER MATRIX APPROACH
Here we assume that, to first order, the dynamics of each beamline element n can be represented by amatrix Φn, known as the transfer matrix for the element. Thus, the action of the element on a particle withphase space coordinates z would be given by the matrix-vector product Φnz. In many situations it ispossible to explicitly calculate the transfer matrix for a particular beamline element a priori, either
analytically or numerically using the principles of linear systems. Moreover, it is found that, with aconjugation operation, the rms moments of the beam can be propagated using the same transfer matrix Φn.
For practical numerical simulation we usually separate the action of the external elements and thespace charge effects. That is, the beam is propagated through the element according to its external forces,
and corrections are applied to account for the space charge along the way. We find that it is possible torepresent the action of the linearized space charge forces as a transfer matrix. Consequently, simulating theevolution of the moments up to second order requires determination of the transfer matrices for eachbeamline element and the transfer matrix that accounts for space charge.

7.1. Transfer Matrices and External Elements
Ignore space charge for the moment. Referring to Eqs. (16), the equations of motion for the x plane
can then be put into the matrix-vector form
(103)x ′( s ) = A( s )x( s ) ,

where

26

LA-UR-02-4979( 104)⎛ x( s ) ⎞⎟⎟x ( s ) ≡ ⎜⎜⎝ x ′( s ) ⎠

and

1 ⎞
⎛ 0
A ( s ) ≡ ⎜ − k 2 ( s ) − γ ′( s ) ⎟ .
⎜ x⎟γ ( s) ⎠⎝In principle, if we know kx2(s) and γ(s) there always exists a solution to Eq. (103) of the form [8](105)x ( s ) = Φ( s )x (0) ,where x(0) is the initial value of x(s), and Φ(s) is the fundamental matrix of the system; it is a matrix
function of the independent variable s having the following properties:(106)

Φ′( s ) = A ( s )Φ( s ) ,Φ(0) = I ,where I is the identity matrix. We use this fact to develop a numerical algorithm for simulation. If wedetermine the matrix Φ for each beamline element (say by explicit calculation or by numerical integrationof the above equation for Φ), then the particle coordinates can be propagated by simple matrix

Φ1

Φ0
x0

x1

Φ2

ΦN
x3

x2

xN-1

xN

Figure 5: transfer matrix representation of the beamline

multiplication. (This fact comes from the semi-group property of the fundamental matrix.)
We can discretize the continuous matrix-vector system (105) by considering only points at the entrance
and exit locations of each beamline element. Let Φn denote the fundamental matrix solution for beamline
element n having length ln evaluated at the final location, that is(107)

Φ n ≡ Φ(l n )The constant matrix Φn is known as the transfer matrix for the element n. Now let xn denote the state of
the particle at the entrance of the nth beamline element(108)

x n ≡ x( s n ) ,where sn is the location of the nth element's entrance. This situation is shown in Figure 5. Propagation of
the particle state from one element to the next is given by the set of discrete equations(109)

x n +1 = Φ n x n

n = 0,1,K

Thus, instead of a set of continuous ordinary differential equations we now have a set of discrete transfer
equations to describe the particle motion. By the linearity of electromagnet forces we can describe theaction of each element by its transfer matrix Φn then use a separate process to determine the space-charge
effects down the beamline.
Many times we can compute the transfer matrix for a beamline element analytically. For example, inthe case of an ideal quadrupole lens where the forces are linear throughout and the fringe fields are
negligible, then the transfer matrix is found to be27

LA-UR-02-4979

(110)

⎛
⎜ cos kl⎜
⎜ − k sin kl⎜
0Φ quad = ⎜
⎜
0
⎜
0
⎜
⎜
0
⎝

1
sin klk
cos kl00cosh kl0

k sinh kl

0

0

0

0

0

⎞
0 0⎟
⎟
0
0 0⎟
1
sinh kl 0 0 ⎟ ,⎟kcosh kl 0 0 ⎟
⎟
0
1 l⎟
0
0 1 ⎟⎠
0where k is the quadrupole focusing strength, l is the length of the quadrupole, and where the quadrupole isfocusing in x and defocusing in y. Note that in this ideal case there is no coupling between phase planes. A
description including misalignments would contain nonzero values in the off-diagonal blocks.
7.2. Space Charge Impulses
In the previous subsection we saw that, barring space charge, the transfer matrix for each element canbe used to propagate a particle's state through the element. In the following subsections we find that thesame transfer matrix can also be used to propagate the moments of the beam distribution through theelement. Consequently, we would like to formulate a transfer matrix representation for space charge
effects. Here we develop such a matrix.

A straightforward approximation that simulates the action of space charge is to "kick" the beam at
regular intervals. Specifically, we apply the space charge effect through a section of length ∆s all at once,
as a momentum impulse. This technique is equivalent to modeling space charge as a defocusing thin lens.
Of course to maintain accuracy is it necessary that ∆s be sufficiently small. The magnitude of the impulse
is determined by returning to the Newton's inertia equation
(111)

dpα
= Fs ,αdt

α ∈ {x, y , z} ,where pα is the momentum component in the α direction and Fs,α is the self force in the α direction. To
illustrate the computation we consider specifically the x phase plane. We start with the expression for theforce on a beam particle given by Eqs. (12), considering only the term specific to the self force.
Approximating the derivative dp/dt by finite the differences ∆p/∆t and using the linearized electric field of
Eq. (59) transforms the above equation into the following approximation:(112)∆p x ≈qγ

2

E x ∆t ≈

q

xE x

γ

x2

2

∆t x .

Recognizing that ∆t=∆s/v and converting to our normalized momentum x'≡px/p we have∆x ′ ≈

(113)
≡

1

xE x

q

γ β γmc
2

2

2

x2

∆s x ,∆sx,
f sc, xwhere we have defined the thin lens focal length due to space charge as fsc,x. Thus, the action of spacecharge in the x plane can be written in transfer matrix form as(114)

1⎛ x n +1 ⎞ ⎛⎜ ∆s⎟⎟ = ⎜⎜⎜⎝ x n′ +1 ⎠ ⎜⎝ f sc , x

0⎞ x⎟⎛ n ⎞
1 ⎟⎜⎜ x ′ ⎟⎟ ,⎟⎝ n ⎠⎠

where ∆s is the length over which the space charge effects are being applied.

28

LA-UR-02-4979
In the z phase plane the calculations are somewhat different although the result is the same. There theself-force on a particle is Fs,z=qEz, however ∆z'=(1/γ2)∆(∆p/p) so we end up with the same factor of γ in theforce terms in all planes. Using the value of Λ for the equivalent uniform beam and the beam perveance K
defined by Eq. (15), we compute the thin lens focal lengths as
1
f sc, x

(115)

1
f sc, y1
f sc, z3/ 2⎡< y2 > < z2 > ⎤
,
,1⎥ ,
RD ⎢
2
2
⎣< x > < x > ⎦3/ 2⎡ < z2 > < x2 > ⎤
,
,1⎥ ,
RD ⎢
2
2
⎣< y > < y > ⎦=K ⎡1 1 ⎤2 ⎢⎣ 5 < x 2 > ⎥⎦=K ⎡1 1 ⎤⎢⎥2 ⎣5 < y2 > ⎦=K ⎡1 1 ⎤2 ⎢⎣ 5 < z 2 > ⎥⎦3/ 2⎡< x2 > < y2 > ⎤,1⎥ ,RD ⎢ 2 ,2⎣< z > < z > ⎦The full transfer matrix Φsc that applies the space charge kick over a distance ∆s is then

(116)

1
⎛
⎜
⎜ ∆s / f sc, x⎜0Φ sc = ⎜
0
⎜
⎜
0
⎜
⎜
0
⎝

0
1
0

0
0
1

0 ∆s / f sc, y
0
0
0
0

0
0
0

0
0
0

1
0
0
1
0 ∆s / f sc, z

0⎞
⎟
0⎟
0⎟
⎟.
0⎟
0 ⎟⎟
1 ⎟⎠

Since fsc,x depends upon 〈x2〉 and 〈xEx〉, it must be recalculated whenever the dimensions of the beam bunchchange appreciably. Likewise for the other defocusing strengths fsc,y and fsc,z. Specifically we are sayingthat the space charge transfer matrix is a function of the covariance matrix, that is(117)Φsc = Φsc (σ) .

7.3. Equations for Centroid Motion
We find the equations that propagate the mean value vector simply by taking the moment of Eq. (109)for the full phase coordinate vector zn. Since the transfer matrix Φn does not depend upon the phase spacecoordinates, we get(118)

z n +1 = Φ n z n ,

or
(119)

z n +1 = Φ n z n ,

where
(120)

z n ≡ z( sn ) ,the value sn being the location of the entrance to the nth beamline element. Recall that there are no spacecharge effects for centroid motion. Therefore, we do not apply any space charge kicks when propagatingthese moments. We simply multiply z n by the transfer matrix Φn for each element n in the beam line.
7.4. Equations for the Second-Order Moments

For simplicity, once again consider only the x phase plane. The continuous evolution equation for σxx
can be found by direct differentiation of σxx with respect to s, then applying Liouville's theorem and Eq.
(103). From the definition of σxx we find29

LA-UR-02-4979
σ xx ′ = x ′x T + xx ′ T(121)= Α xx T + xx T A T= Aσ xx + σ xx A T

The solution to this matrix differential equation is given (formally) by(122)

σ xx ( s ) = Φ( s )σ xx (0)Φ T ( s ) ,where Φ(s) is the fundamental matrix solution of Eq. (106), and σxx(0) is the initial covariance matrix forthe x plane. (This fact may be checked by direct differentiation.) The above equation also holds for theentire set of second moment σ when using the fully augmented fundamental matrix Φ(s). Recalling that thetransfer matrix for an element is found by evaluating the fundamental matrix at the end of the element (i.e.,Φn=Φ(ln), where ln is the length of the nth element) leads to the following transfer equation for the secondorder moments:(123)

σ n +1 = Φ n σ n ΦTn ,

where
(124)

σ n ≡ σ( s n ) .

We see that advancing the state of the second-order moments through the nth beamline element isaccomplished by transpose conjugation of σ with the transfer matrix Φn. However, we must still include
the effects of space charge.
Space charge effects are imposed using a transfer equation analogous to (123). However, here weconjugate σ with the space charge transfer matrix Φsc. Recall that the space charge transfer matrix Φsc is afunction of the covariance matrix σ. Therefore, for accurate simulation of long beamline elements it may benecessary to split the beamline element into several subsections, applying the space charge kick after eachsubsection. Specifically, we advance σ by the transfer matrix for a subsection, conjugate by Φsc, and repeat
until we are through the element.
7.5. Extensions to the Inhomogeneous CaseIn some particle beam situations the modeling equations are inhomogeneous, that is, they contain aforcing term. For beamline elements, this is the case for dipole magnets, for example when used as beamsteering magnets (for example, see Appendix A.3 ). Many control systems also have modeling equations of
this type. Let the transfer equations for an inhomogeneous element be given as(125)

z n+1 = Φ n z n + Γ n u nwhere zn is the 6×1 column vector of phase space coordinates, Φn is the 6×6 transfer matrix, Γn is a 6xM
matrix, and un is a column vector of length M. Note that this is still a linear system (in zn and un) only now
we have a forcing term driven by un. The vector un represents some external parameters for the element,
perhaps a control parameter or perhaps an unknown noise source (in this case un would be a randomvariable).
We can compute the propagation equations for σn simply by unwinding the definition and using theabove transfer equation. We get
σ n + 1 = z n +1 z n +1 T ,(126)= (Φ n z n + Γ n u n )(Φ n z n + Γ n u n )T ,= Φnσ n ΦnT + Φn Χ n Γ nT + Γ n Χ nT ΦnT + Γ n Υ n Γ nT ,

where

30

LA-UR-02-4979
Χ n ≡ z nunT ,(127)Υ n ≡ un unT .

The 6×M matrix Χn is the cross-correlation matrix between the phase-space coordinates and the input
vector. The M×M symmetric matrix Υn is the covariance matrix of the drive vector. Note that if un is ascalar value, say un, then Χn= u n z and Υn=un2.
The matrix Γn is given by the structure of the system, representing coupling between the particle state
and the external drive source un. Thus, to use the above transfer equations for σn it is necessary to
determine Xn and Yn. Typically Xn is zero, it represents the correlation between the drive vector un and thecoordinates in phase space z. Most systems are designed to avoid any such correlations. When un is arandom variable (representing noise or other unknown), the matrix Yn is the autocorrelation. For most
systems it can be determined, at least approximately, through measurement. When un is deterministic, Yn
is essentially just un squared, which is easily determined.
7.6. Homogeneous Coordinates
Here we again consider inhomogeneous systems, however, now we assume the special form

(128)

z n +1 = Φ n z n + u n ,where un∈ℜ6 is the external drive vector. That is, the external driving forces can be described as a vector
in phase space. This situation occurs when the effects of an inhomogeneous element behave as translations
in phase space, such as those produced by an ideal dipole magnet. In this situation we can form a
convenient augmented state variable system that has the same form as the homogeneous transfer equations.
Mathematicians typically use homogeneous coordinates to parameterize the projective spaces Pn. They
are also widely used in computer graphics for three-dimensional rendering, since translation, rotation, andscaling can all be performed by matrix multiplication [18]. The n-dimensional real projective space ℜPn
can be described as a set equivalence relations [x0, …, xn] on ℜn+1 where [x0, …, xn]~[wx0, …, wxn] for all
real w≠0, and such that not all the xi are zero. Thus, the points of the project space ℜPn are seen to be thelines in ℜn+1 that pass through the origin. These equivalence classes are known as the homogeneous
coordinates of the projective spaces. (Another equivalent description of the projective space ℜPn is foundby identifying all the antipodal points of the sphere Sn.) The projective space ℜPn can be considered adifferentiable manifold with the atlas consisting of n+1 charts {Ui,φi}i=0n where Ui is the set of equivalence
relations [x0, …, xn] such that xi≠0, and φi:Ui → ℜn is the bijective coordinate map(129)

φi : [ x0 ,K, xi ,K x n ] a ( x0 / xi ,K, xˆ i ,K x n / xi ) .

The caret indicates omission of the coordinate. Note that the union ∪i=0nUi covers all of ℜPn. More aptly,
note that the coordinates of the n+1 chart consist of the following equivalence relations:(130)

φ i−1 ( x 0 ,K, x n −1 ) = [ x 0 ,K x n −1 ,1]∀( x 0 ,K x n −1 ) ∈ ℜ n .

Thus, Ui is seen to be the set of all lines in ℜn+1 passing through the plane {(x0,…,xn)∈ℜn+1⏐xn=1}. We
shall use the homogeneous coordinates of this chart.
Let the augmented phase space coordinate ζ be
(131)

⎡z⎤
ζ ≡ ⎢ ⎥ = (x
⎣1⎦

x′

y

y′ z

z ′ 1)T .

Then system (128) can be written in the form
(132)

ζ n+1 = Θ n ζ n ,where the 7×7 matrix θn is defined as the augmented square matrix

31

LA-UR-02-4979

(133)

⎡Φ
Θn = ⎢ n
⎣ 0

un ⎤
.
1 ⎥⎦

This is exactly the form of the homogeneous transfer equations. Therefore, if we only encounter
inhomogeneous transfer systems of the special form given above, by employing homogeneous coordinates
we can convert back to homogeneous transfer equations.
For envelope calculations consider the modified covariance matrix τ formed from the phase space
vector ζ according to(134)τ ≡ ζζ T .

To simplify the discussion, consider only the x-plane phase space. In the homogeneous coordinates the xplane phase coordinates are given by
(135)

⎡x ⎤
ξ ≡ ⎢ ⎥ = (x⎣1⎦

x ′ 1)T .

Thus, the modified x-plane correlation matrix τxx is formed by(136)

τ xx ≡ ξξ T

⎛ x2⎜
= ⎜ xx ′
⎜
⎜ x
⎝

xx ′
x ′2
x′

x ⎞⎟
⎡σ
x ′ ⎟ = ⎢ xx
T
⎟
⎣x
1 ⎟
⎠

x⎤
.
1 ⎥⎦

We see that in the homogeneous phase space coordinates the covariance matrix τ contains the original
covariance matrix σ plus the mean value vector z . Moreover, transpose conjugation of the covariancematrix τxx by the transfer matrix Θn yields the following augmented system:(137)⎡Φ σ Φ T + Φ n zu T + uz T Φ n T + uu T
τ n +1 = Θ n τ n Θ n T = ⎢ n n n
z nT ΦnT + u nT⎣⎢

Φn z n + u n ⎤⎥.
1⎦⎥

The upper left block is exactly the covariance transfer equation for the inhomogeneous system (128). The
diagonal blocks are the transfer equation, and its transpose, for the mean value vector. Consequently, boththe transfer equation for the mean value evolution and the covariance evolution are included in themodified transfer system.
To make the use of homogeneous coordinates more explicit, consider the case of an ideal dipolecorrecting magnet in the x phase plane. The effect of such a device is to add an impulsive kick in theparticle momentum of intensity ∆x'. Thus, the modified transfer matrix for the x phase plane would appear
as

(138)

⎛1 0 0 ⎞⎜⎟Θ dipole = ⎜ 0 1 ∆x ′ ⎟ .
⎜0 0 1 ⎟⎝⎠The effects of the ideal dipole are found by transpose conjugation of the modified covariance matrix τxx by
the above transfer matrix; we have32

LA-UR-02-4979

(139)

⎛ 2
⎛ 1 0 0 ⎞⎜ x⎜⎟
Θ dipole τ xx Θ dipole T = ⎜ 0 1 ∆x ′ ⎟⎜ xx ′⎜⎜ 0 0 1 ⎟⎜ x
⎝
⎠
⎝
⎛
x2
⎜
= ⎜ x x ′ + x ∆x ′
⎜
⎜
x
⎝

xx ′
x ′2
x′

x ⎞⎟⎛ 1 0 0 ⎞⎜⎟x ′ ⎟⎜ 0 1 0 ⎟ ,⎟1 ⎟⎜⎝ 0 ∆x ′ 1 ⎟⎠⎠xx ′ + x ∆x ′x ′2 + 2 x ′∆x ′ + ∆x ′2
x ′ + ∆x ′⎞⎟′′x + ∆x ⎟ .
⎟
1 ⎟
⎠
x

We can clearly see the effects of an impulsive kick on the entire beam, both in the centroid behavior andthe rms envelope behavior.

8. SUMMARY AND CONCLUSION
Beam envelope simulation as described here is a computationally inexpensive way to obtain particle
beam behavior to first order. Rather than propagating an ensemble of particles then computing the statistics
to determine the rms beam properties, we propagate the rms beam properties directly. We have seen thattheoretically that this is a valid approach. However, there are clearly limitations to simulation results
generated in such a manner. The techniques described here are intended to give a quick insight into generalbeam behavior and not to be used as a tool for detailed dynamics analysis. Consequently the simulation
methods here are appropriate for initial off-line design studies and for on-line model reference control
applications, where fast real-time response is essential.
The validity off our analysis typically comes into concern because, by practical considerations, we areforced to make some simplifying assumptions to implement the simulation procedures. The first
assumption we make is the applicability of Liouville's theorem to form the evolution equations for the
distribution's moments. Technically, Liouville's theorem is valid only in 6N dimensional phase space forcollisionless system, where N is the number of beam particles. However, it is approximately true wheneverthe collective fields of the beam can be accurately described by smooth functions. Thus, we have assumed
that our beam is populous enough, and compact enough for this condition to be true. Another circumstance
where this assumption holds is for relatively cold beams. Both situations assume that the Debye length of
the beam is small enough so that any particle see mostly collective fields and there are few collision-likeencounters. The next critical assumption comes in the form of ellipsoidal symmetry. That is, we assume
that the beam exhibits ellipsoidal symmetry in configuration space. Under this assumption we are able toanalytically compute the effects of space charge on the moment dynamics. Unfortunately it is known that
ellipsoidally symmetric beams are not, in general, stationary beams [21], the only one of this type is theKV, or micro-canonical distribution. Thus, we must assume that our true beam is very close to an
ellipsoidal one. Fortunately ellipsoidal beams can accurately represent many laboratory beams.
The last major assumption is that of linearity. Specifically, we assume not only that all external forces
are linear, but also that internal forces (i.e., self forces) are linear. Consequently any nonlinear effects froma beamline element cannot be modeled, for example, fringe fields, higher order field components, etc. Ourlinear model for the internal fields essentially constitutes the assumption of constant rms emittances. Thisis true for both the rms envelope equations and the transfer matrix method. There is no known practicalmethod for simulating emittance growth using only beam statistics. Currently, only full multiple particle
simulations have this capability. However, it is possible to assign emittance growths based on analytic
approximations or a priori knowledge of emittance values. In the envelope equations we just assign theemittance values directly, for the transfer matrix approach we multiply the second-order covariant by thegrowth factor.

ACKNOWLEDGEMENTS
The authors wish to thank Tom Wangler for many valuable discussions on beam physics and for histime reviewing the manuscript.

33

LA-UR-02-4979

APPENDIX A: TRANSFER MATRICES FOR COMMON BEAMLINE ELEMENTS
We describe here the transfer matrices for some beamline elements in common use. As was done withthe σ matrix, we decompose a general transfer matrix Φ into 2×2 composite matrices according to thefollowing:(A.1)⎛ Φ xx⎜Φ = ⎜ Φ yx
⎜Φ
⎝ zx

Φ xy
Φ yy
Φ zy

Φ xz ⎞⎟
Φ yz ⎟ .
Φ zz ⎟⎠

In most cases only the block diagonal sub-matrices Φxx, Φyy, and Φzz are nonzero. Nonzero values of theoff-diagonal sub-matrices indicate coupling between phase planes.
Many transfer matrices exhibit a semi-group property. Specifically, if we denote the transfer matrix ofan element of length ∆s as Φ(∆s), then the transfer matrix for two such elements is Φ(∆s)Φ(∆s). In general,
for arbitrary lengths ∆s1 and ∆s2 we have(A.2)

Φ( ∆s1 + ∆s 2 ) = Φ( ∆s1 )Φ( ∆s 2 ) .
A.1 Drift Space

The transfer matrix Φ for a drift space of length ∆s is given in terms of its nonzero diagonal blocks(A.3)⎛ 1 ∆s ⎞⎟⎟ ,Φ xx = ⎜⎜⎝0 1 ⎠⎛ 1 ∆s ⎞⎟⎟ ,Φ yy = ⎜⎜⎝0 1 ⎠⎛ 1 ∆s ⎞⎟⎟ ,Φ zz = ⎜⎜⎝0 1 ⎠

where the sub-matrix for the z plane assumes that the coordinates are (z,z'). If the coordinates are (z,∆p/p)
then the sub-matrix is(A.4)⎛ 1 ∆s / γ 2 ⎞⎟.
Φ zz = ⎜⎜1 ⎟⎠⎝0This transfer matrix has the semi-group property.
A.2 Quadrupole
In the case of an ideal quadrupole lens where the forces are linear throughout and the fringe fields are
negligible, the transfer matrix block diagonals are(A.5)⎛
cos klΦ foc = ⎜⎜⎝ − k sin kl

1⎞
sin kl ⎟,k⎟
cos kl ⎠⎛cosh kl
Φ def = ⎜⎜⎝ k sinh kl

1⎞sinh kl ⎟,k⎟cosh kl ⎠⎛1 l ⎞⎟⎟ ,Φ zz = ⎜⎜⎝ 0 1⎠

where Φfoc represents the block diagonal of the focusing plane and Φdef represents the block diagonal of thedefocusing plane of the quadrupole. The quantity l is the length of the quadrupole and k is the quadrupole
focusing strength, it has the values
k mag =qGβγmc

k esq =

1` qV0βa γmc 2(A.6)

where kmag is the value for a magnetic quadrupole and kesq is the value for an electrostatic quadrupole. Here
G is the magnetic field gradient at the beam axis, V0 is the electrostatic electrode potential and a is theaperture of the quadrupole.

34

LA-UR-02-4979

A.3 Dipole Steering Magnet (As a Thin Lens)
We may treat a dipole magnet as a thin lens by again approximating the time derivative of momentum
as a finite difference in the equations of motion. The equation of motion for a particle in a dipole field ofstrength B is
(A.7)

dp
= qv × B ,dt= qβc( x ′, y ′,1 + z ′) × ( B x , B y ,0) ,where we have assumed the magnetic field is constant throughout the magnet.
derivative by a finite difference gives us the equations

Approximating the∆p x = − q(1 + z ′) B y ∆s ,(A.8)

∆p y = q(1 + z ′) B x ∆s ,∆p z = q( x ′B y − y ′B x ) ∆s ,Using the relations ∆x'=∆px/p, ∆y'=∆py/p, ∆z'=(1/γ2)∆pz/p, we can form the transfer system

(A.9)

0
⎛ x n +1 ⎞ ⎛ 1⎜⎟ ⎜1⎜ x n′ +1 ⎟ ⎜ 0⎜ y ⎟ ⎜00
⎜ n +1 ⎟ = ⎜0⎜ y n′ +1 ⎟ ⎜ 0⎜⎟ ⎜00⎜ z n +1 ⎟ ⎜⎜ z ′ ⎟ ⎜ 0 ω ∆s / γ 2x⎝ n +1 ⎠ ⎝

0
0

0
0

1
0
0

0
1
0

0 ω y ∆s / γ 200 ⎞⎛ x n ⎞ ⎛ 0 ⎞⎟⎜ ⎟ ⎜⎟
0 ω x ∆s ⎟⎜ x ′n ⎟ ⎜ ω x ∆s ⎟00 ⎟⎜ y n ⎟ ⎜ 0 ⎟⎟⎜ ⎟ + ⎜⎟,0 ω y ∆s ⎟⎜ y ′n ⎟ ⎜ ω y ∆s ⎟⎟1
0 ⎟⎜⎜ z n ⎟⎟ ⎜⎜ 0 ⎟⎟0
1 ⎟⎠⎜⎝ z ′n ⎟⎠ ⎜⎝ 0 ⎟⎠where
(A.10)

ωx ≡ −

qB y

and

βγmc

ωy ≡

qB x

βγmc

are the cyclotron phase advances in the x and y directions respectively.
Since the contributions from the phase space coordinates are usually small we can often approximate thedipole effect using only the additive term in the above. Referring to Section 7.5, in this case we have Γ=I
and(A.11)

⎛0
⎜
⎜0
⎜0
X=⎜
⎜0
⎜0
⎜
⎜0
⎝

xω x
x ′ω x00xω yx ′ω yyω xy ′ω xzω xz ′ω x000yω yy ′ω y
zω y

0

zω y

0 0⎞
0
⎛0
⎟
⎜
2
0 0⎟
⎜0 ωx
⎜0
0 0⎟
0
⎟, Y=⎜0 0⎟⎜ 0 ω xω y⎜0
0 0 ⎟⎟
0
⎜
⎜
⎟
0 0⎠
0
⎝0

0
0
0 ω xω y
0
0
0

0

ω 2y

0

0

0

0 0⎞
⎟
0 0⎟
0 0⎟
⎟.
0 0⎟
0 0⎟
⎟
0 0 ⎟⎠

A.4 Dipole Steering Magnet (As a Thick Lens)
We present the full solutions for particle trajectories through constant magnetic fields in the transverse
directions. These solutions can be computed analytically and put into transfer matrix form. From the full
solutions one can derive transfer matrices for dipole magnets having fields directed in only one plane.
Starting from the Lorentz forces the differential equations describing the motion are

35

LA-UR-02-4979(A.12)x ′′ = −ω x (1 + z ′),y ′′ = −ω y (1 + z ′),z ′′ = ω x x ′ + ω y y ′ ,where(A.13)

ωx ≡ −

qB y

βγmc

and

ωy ≡

qB x
βγmc

are the cyclotron frequencies in the x and y directions. The above equations can be solved to give us thetransfer matrix(A.14)2⎛⎛ 2⎞⎜ 1 ⎜ ω x sin ω ∆s + ω y ∆s ⎟
c
3
2
⎜
⎜ ωc
ω c ⎟⎠⎝⎜⎛ω2⎜
ω 2y ⎞⎟
x⎜cos ω c ∆s + 2⎜0⎜ ω c2ω c ⎟⎠⎜⎝⎜⎜ 0 ω x ω y (sin ω ∆s − ω ∆s )
c
c
⎜
ω c3
Φ=⎜
⎜
ω xω y⎜0(cosω c ∆s − 1)
ω c2
⎜
⎜
ωx
⎜0
(cosω c ∆s − 1)
⎜
ω c2
⎜
ωx
⎜0
sin ω c ∆s
⎜
ωc
⎝

0

0

1

0
0
0

ω xω y(sin ω c ∆s − ω c ∆s ) 03ωcω xω yω c2(cosω c ∆s − 1)

⎛ ω 2yω x2 ⎞⎟
⎜
s
ω
sin
∆
+
∆s
c
⎜ ω c3
ω c2 ⎟⎠⎝⎛ ω 2yω x2 ⎞⎟
⎜
s
ω
cos
∆
+
c
⎜ ω c2
ω c2 ⎟⎠⎝

ωy(cosω c ∆s − 1)ω c2
ωysin ω c ∆s
ωc

0

0

0
1
0

⎞
ωx
⎟
(
)
s
ω
cos
1
∆
−
c
⎟
ω c2

⎟
⎟
⎟
⎟
⎟
ωy
⎟
(
)
s
ω
cos
1
∆
−
c
⎟
ω c2
⎟,
⎟
ωy
sin ω c ∆s ⎟
ωc⎟⎟sinω c ∆s
⎟
⎟
ωc
⎟
⎟
cos ω c ∆s⎟⎠ωxsin ω c ∆s
ωcand drive vector(A.15)⎞⎛ωx⎜ 2 (cos ω c ∆s − 1)⎟
⎟
⎜ ωc
⎟
⎜ ωx
sin ω c ∆s ⎟
⎜
⎟
⎜ ωc
⎜ωy
(cos ω c ∆s − 1)⎟⎟ ,u = ⎜ω2
⎟
⎜ c
⎟
⎜ ωy
s
ω
∆
sin
c
⎟
⎜
ω
c
⎟
⎜
⎜ sinω c ∆s − ∆s ⎟⎟⎜ω⎜ cos ωc ∆s − 1 ⎟
c
⎠
⎝

where
(A.16)

ω c = ω x2 + ω 2y .

Note that again Γ=I. Note also that in the limit as ∆s→0 the transfer matrix and drive vector approach thethin lens approximation but with the drift component.

36

LA-UR-02-4979(A.17)

∆s
⎛1
⎜
1
⎜0
⎜0
0
lim Φ = ⎜00∆s → ∞⎜⎜00⎜⎜ 0 ω x ∆s
⎝

0
0
1
0
0
0

0
0 ⎞
⎟
0 ω x ∆s ⎟
0
0 ⎟
∆s
⎟,
1
0 ω y ∆s ⎟
0
1
∆s ⎟⎟
ω y ∆s 0
1 ⎟⎠00

A.5 Bending Magnet
In a bending magnet thedesign trajectory is where the
force
from
centripetal
acceleration exactly balances
the magnetic force, it is an arc
with radius R0. To analyze abending magnet we usuallyemployaperturbational
analysis around the design
trajectoryincylindricalcoordinates (r,θ,z). Figure 6
demonstrates how we construct
the beam coordinates (x,y,z) as
perturbations
aroundthesynchronous trajectory of r=R0.

⎛ 0 ⎞
⎟
⎜
⎜ ω x ∆s ⎟⎜ 0 ⎟⎟.
lim u = ⎜∆s →∞⎜ ω y ∆s ⎟
⎜ 0 ⎟
⎟
⎜
⎜ 0 ⎟
⎠
⎝

Bz
y

z

x

s

R0

z
r

θ

α
Clearly the bend angle α isa parameter of the bending
magnet.
The radius ofcurvature R0 is also a parameter
Figure 6: sector bending magnetof the bend; along with thedesign energy it determines thestrength of the magnetic field B0=Bz(R0) on the design trajectory. However, we also require a thirdparameter, the field index n. The field index is defined(A.18)

n≡−R0 ∂B z.
B0 ∂r r = R0

The field index is simply the normalized derivative of the bending field evaluated at the design radius. In
order to provide focusing in both transverse planes it is necessary that the bending field decrease in theradial direction. Kerst and Serber first studied these effects for the betatron [14], they found for focusing inboth transverse planes it is necessary that 0<n<1.
To find the transfer matrix for a bending magnet we start with the equations of motion in cylindricalcoordinates. Since magnetic fields cannot accelerate, γ is constant. We also assume that there is nomagnetic field in the θ direction. With these considerations the equations of motion are

γm&r& − γmrθ& 2 = qrθ&B z ,(A.19)

γmrθ&& + 2γmr&θ& = qz&Br − qr&B z ,γm&z& = − qrθ&Br .

The first-order variation in both magnetic field components Br and Bz can be written in terms of the fieldindex n using the fact that ∇×B=0. We have37

LA-UR-02-4979
B z = B0 − nB0(A.20)
Br = −nB0

r − R0
,
R0

z
,
R0

where we have assumed that Br(z=0)=0 by design. Now we assume a perturbation around the design
trajectory according tor = R0 + εx + O (ε 2 ),(A.21)

θ = ωct + εz+ O (ε 2 ) ,R0z cyl = εy + O (ε 2 ) ,where ε<<1 is a dimensionless parameter used to keep track of the order of approximation. Note that wehave identified the cylindrical coordinate z with subscript cyl whereas the beam coordinate z is unadorned.
Note also that the path length parameter s is related to θ by
(A.22)

s = R0 (α − ω c t )⇒dd= −ω c R0.
dt
ds

Also we have the convenient physical relations
qB0
,
γm
ds &
= θR 0 = ω c R 0 ,v=
dtγmv βγmc
=
,
B0 =
qR0
qR0

ωc ≡
(A.23)

where ωc is the cyclotron frequency.
The transfer matrix solution to the above linear perturbation is given as follows:(A.24)⎛⎜ cos k x ∆s
Φ xx = ⎜⎜ − k sin k ∆s
x
⎝ x

⎛
1
⎞
⎜ cos k y ∆ssin k x ∆s ⎟
=
Φ
kx
yy
⎜
⎟
⎜ − k y sin k y ∆s
cos k x ∆s ⎟⎠⎝1⎞
sin k y ∆s ⎟ky⎟
cos k y ∆s ⎟⎠

2
⎛
⎞ ⎞
⎛
⎜ 1 − k x β ∆s − sin k x ∆s + ⎜1 − 1 ⎟∆s ⎟⎜ k 2R2 ⎟ ⎟Φ zz = ⎜k x3 R02
x 0 ⎠
⎝
⎜0
⎟
1
⎝
⎠

In addition, we also have the off-diagonal blocks(A.25)⎛⎜0Φ xz = ⎜⎜⎜0⎝

R0 (1 − cos k x ∆s ) ⎞⎟⎛ R0 sin k x ∆sk x2
⎟ Φ = ⎜−kxzx⎜
R0 sin k x ∆s ⎟
⎜
0
⎟
⎝
kx
⎠

where

38

R0 (1 − cos k x ∆s ) ⎞
⎟
k x2
⎟
⎟
0
⎠

LA-UR-02-4979(A.26)

kx =

1− n
,
R0

ky =

n
,
R0

∆s = αR0 .
The above description is for a bend in the horizontal plane to the right (positive x direction). For a bend to
the left replace Φxz and Φzx by -Φxz and -Φzx, respectively.
A.6 RF Gap (As a Thin Lens)

Here we treat the effects of a general RF gap as a thin lens. The energy gain in the gap ∆W is given interms of the Panofsky equation [22](A.27)∆W = qE0TL cosφ ,where q is the unit charge, E0 is the longitudinal electric field in the gap, T is the transit time factor (theratio of the energy gained in the RF gap to that of a DC gap of the same field amplitude), and φ is the RFphase of the synchronous particle at the center of the gap. Along with h, the (integer) number of the fieldharmonic, we assume that these are the parameters of the RF gap.
The electric fields in the gap cause longitudinal focusing and radial defocusing, as well as the change
in energy. Here we present the thin lens focusing constants to account for these effects. The values of x', y'and z' tend to decrease with increasing longitudinal momentum simply because they are normalized by thisquantity. However, the unnormalized values px, py, and ∆p are constant with an increase in longitudinal
momentum. Thus, the action of the RF gap is produced by first unnormalizing the momentum components,
transforming by a thin lens, then normalizing the momentum components with respect to the new, larger
longitudinal momentum. For example, for the x phase plane we havep x ,n = pi x n′ = ( βγ ) i mc x n′ ,(A.28)

p x ,n +1 = p x ,n + k x x n ,x n′ +1 =p x ,n +1
pf=

p x ,n +1( βγ ) f mc,

where the subscript i refers to quantities with initial energy and the subscript f refers to quantities with thefinal energy. The following transfer matrix captures this sequence of operations:(A.29)⎛1⎜
Φ xx = ⎜ 0
⎜
⎝

0
⎞⎛ 1
⎟⎜
1
⎟⎜
( βγ ) f mc ⎟⎠⎝ k x

1⎛⎜= ⎜ 1 kx
⎜ ( βγ ) f mc⎝

0 ⎞⎛ 10⎞⎟⎟⎜⎜⎟,
1 ⎠⎝ 0 ( βγ ) i mc ⎟⎠0 ⎞( βγ ) i ⎟ .
⎟( βγ ) f ⎟⎠Likewise, for the other phase planes we have(A.30)10 ⎞⎛⎜ 1 ky()i ⎟βγ
Φ yy = ⎜⎟⎜ ( βγ ) f mc ( βγ ) f ⎟
⎝
⎠
1
0
⎛
⎞
2
⎜ γ2
⎟
(
)
k
γ
βγ
i
i ⎟,
Φ zz = ⎜⎜ γ 2f ( βγ ) f mc γ 2f ( βγ ) f ⎟⎝⎠Note that39

LA-UR-02-4979( βγ ) i = γ i2 − 1,( βγ ) f = γ 2f − 1,

(A.31)

γ i, f = 1 +Wi , fER,

W f = Wi + ∆W ,where ER is the rest energy of the beam particle (mc2). The values of the focusing coefficients are found byintegrating the electromagnetic fields in the gap. We have (approximately, for example see [22])(A.32)

kx
k
=
,
mc mc

ky
mc

=

k
,
mc

kz
k
= −2γ 2
,
mc
mc

where
(A.33)kπ≡ h 2 2 2 qE 0 L sin( −φ ) ,mc
β γ mc λβ , γ are the relativistic parameters evaluated at the average energy(A.34)

W = Wi + ∆W / 2 ,and h is the field harmonic at which the gap is operating.
If the longitudinal coordinates are given as (z,∆p/p) the corresponding transfer block is then
(A.35)

1⎛⎜
Φ zz = ⎜ 1 k z⎜ ( βγ ) f mc⎝

0 ⎞( βγ ) i ⎟ .
⎟( βγ ) f ⎟⎠This value may be used to approximate the previous matrix since γ 2 / γ 2 is typically close to unity.

40

LA-UR-02-4979

APPENDIX B: SELF-FIELD MOMENT CALCULATION
We shall calculate the moment 〈zEz〉. The other field moments follow a similar procedure.
(B.1)

1
N

zE z =

2y2 z2 ⎞⎡ ∂φ ( x, y , z ) ⎤ ⎛⎜ x
F ⎜ 2 + 2 + 2 ⎟⎟dxdydz .
⎥
∂z
b
c ⎠
⎦ ⎝a

∫∫∫ z ⎢⎣−

Substituting the value for potential φ from Eq. (71) we have
⎛ x2
y2
z2 ⎞
⎟
F ⎜⎜
+
+
⎛ x2 y2 z2 ⎞t + a 2 t + b 2 t + c 2 ⎟⎠1 qabc + ∞ + ∞ + ∞ ∞ 2⎝⎜⎟dtdxdydz .
z
F
=
+
+
N 2ε 0 −∞ −∞ −∞ 0
(t + a 2 )1 / 2 (t + b 2 )1 / 2 (t + c 2 ) 3 / 2 ⎜⎝ a 2 b 2 c 2 ⎟⎠∫ ∫ ∫ ∫

(B.2) zE z

Since the integrations are independent we may switch the order so the integration with respect to t is last.
Now employ the change of coordinatesx = r t + a 2 sin θ cos φ(B.3)

y = r t + b 2 sin θ sin φdxdydz = t + a 2 t + b 2 t + c 2 r 2 sin θz = r t + c 2 cos θalong with the appropriate limits of integration to yield(B.4)

zE z =[]

1 qabc ∞ 2π π ∞ 2
r cos 2 θF ( r 2 ) F r 2 + tr 2 u (θ ,φ ) r 2 sin θdrdθdφdt ,N 2ε 0 0 0 0 0∫∫ ∫∫where the function u(θ,φ) is defined(B.5)⎛ cos 2 φ sin 2 φ ⎞ cos 2 θ⎟+.
u(θ , φ ) ≡ sin 2 θ ⎜⎜+2b 2 ⎟⎠c2⎝ a

We now switch the order of integration again and apply the substitution(B.6)

τ = r 2 + tr 2 u(θ , φ )

dt =

dτ
r 2 u(θ , φ )to produce(B.7)

zE z =∞ ∞1 qabc 2π π cos 2 θ sin θdθdφ
r 2 F ( r 2 ) F (τ )dτdr .
0 r2
N 2ε 0 0 0 u(θ , φ )

∫ ∫

∫∫

Using the definitions(B.8)

G(r) ≡∫∞F ( s )dsandΓ( F ) ≡r∞∫ G (r )dr ,220

and the fact that
∞

(B.9)

∫ r F (r )G(r )dr = 4 Γ( F ) ,
2

2

2

1

0

the original integration is reduced to(B.10)zE z =

Γ( F ) qabc 2π π cos 2 θ sin θdθdφ
N 8ε 0 0 0 u(θ ,φ )∫ ∫where we have integrated out the distribution dependence.

41

LA-UR-02-4979
Note that the function u(θ,φ) can be written(B.11)

22⎤⎡ a 2 + b2
12
2 ⎤ ⎡b − a
sin
θ
cos
θ
sin 2 θ ⎥ cos 2φ .
++u(θ , φ ) = ⎢
⎥ ⎢
2 2
2
2 2
c
⎦
⎦ ⎣ 2a b⎣ 2a bUsing this decomposition and the fact thatdφ

2π(B.12)

2π∫ A + B cos 2φ = A − B ,202further reduces Eq. (B.10) tozE z =(B.13)=πcos 2 θ sin θΓ( F ) qπ
abc
dθ
1/ 2
0
N 4ε 0⎡ sin 4 θ ⎛ 1
1 ⎞ sin 2 θ cos 2 θ cos 4 θ ⎤+⎥⎢ 2 2 +⎜ 2 + 2 ⎟
b ⎠
c2
c4 ⎦
⎝a
⎣a b

∫

π
cos 2 θ sin θΓ( F ) qπ
abc
dθ
1/ 2
1/ 2
0
N 4ε 0⎛ sin 2 θ cos 2 θ ⎞ ⎛ sin 2 θ cos 2 θ ⎞⎜⎟ ⎜⎟⎜ a 2 + c2 ⎟ ⎜ b2 + c2 ⎟
⎝
⎠ ⎝
⎠

∫

If we employ the substitution(B.14)cos θ =ct + c2tsin θ =

t+c2

sin θdθ =1cdt
2 t + c 2 3/ 2()

we find that
(B.15)zE z =

Γ ( F ) qπ a 2 b 2 c 2 2 ∞
dt
dt .
c
2
1
/
2
0
4ε 0
N
( t + a ) ( t + b 2 )1 / 2 ( t + c 2 ) 3 / 2zE z =

Γ ( F ) qπ a 2 b 2 c 2 2
c R D [a 2 , b 2 , c 2 ] ,
N
6ε 0

∫

or
(B.16)

where RD is the Carlson elliptic integral of the second kind.

42

LA-UR-02-4979

REFERENCES[1] C.K. Allen, “Bunched Beam Envelope Equations Including Image Effects from a Cylindrical Pipe”,
Phys. Rev. E, Vol. 55, No. 6-B (1997), pp. 7591-7605.
[2] D. C. Carey, K. L. Brown and F. Rothacker, Third-Order TRANSPORT - A Computer Program for
Designing Charged Particle Beam Transport Systems, Stanford Linear Accelerator Center Report No.
SLAC-R-95-462 (1995).
[3] B.C. Carlson, SIAM Journal on Mathematical Analysis, Vol. 8 (1977), pp. 231-242.
[4] B.C. Carlson, Math. Comp., Vol. 49 (1987), pp. 595-606.
[5] K.R. Crandall and D.P. Rusthoi, Trace 3-D Documentation, Third Edition, Los Alamos Code Group
LA-UR-97-886 (May, 1997).
[6] C.N. Dorny, A Vector Space Approach to Models and Optimization (Krieger Publishing, Malabar, FL,1980 - original publication by Wiley, 1975), Chapt. 5.
[7] A.J. Dragt, R.L. Gluckstern, F. Neri, and G. Rangarajan, “Theory of Emittance Invariants”, Frontiersof Particle Beams; Observation, Diagnosis and Correction, Lecture Notes in Physics 343, edited by
M. Month and S. Turner (Springer-Verlag, Berlin, 1988), pp. 94-121.
[8] See any introductory text on ordinary differential equations such as E.A. Coddington and N. Levinson,
Theory of Ordinary Differential Equations (McGraw-Hill, New York, 1955).
[9] R.L. Gluckstern, “Scalar Potential for Charge Distributions with Ellipsoidal Symmetry”, Fermilab
Internal Report TM-1402 (1986).
[10] D.D. Holm, W.P. Lysenko, and J.C. Scovel, “Moment Invariants of the Vlasov Equation”, J. Math.
Phys., Vol. 31, No. 7 (July, 1990), pp. 1610-1615.
[11] P.M. Lapostolle, CERN report AR/Int. SG/65-15, Geneva, Switzerland (July, 1965).
[12] W.P. Lysenko and M. Overley, “Moment Invariants for Particle Beams”, Linear Accelerators andBeam Optics Codes, AIP Conference Proceedings 177, edited by C.R. Eminhizer (AIP, New York,
1988), pp. 323-335.
[13] O.D. Kellogg, Foundations of Potential Theory (Dover, 1953) pp. 192-194.
[14] D.W. Kerst and R. Serber, Phys Rev. 60, 53 (1941).
[15] D.G. Luenberger, Optimization by Vector Space Methods (Wiley, New York, 1969), Chapt. 4.
[16] N.D. Pattengale and C.K. Allen, "The Software Anatomy of a Flexible Accelerator Simulation
Engine", Los Alamos National Laboratory Internal Report LA-UR-02-3761.
[17] W.H. Press, S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery, Numerical Recipes in C, Second
Edition (Cambridge University Press, Cambridge, 1992), pp. 261-269.
[18] W. Schroeder, K. Martin, and B. Lorensen, The Visualization Toolkit, Second Edition (Prentice Hall,
Upper Saddle River, NJ, 1998), Section 3.7.
[19] F.R. Sacherer, “RMS Envelope Equations with Space Charge”, IEEE Trans. Nucl. Sci., Vol. NS-18(1971), pp. 1105-1107.
[20] I. Stakgold, Green’s Functions and Boundary Value Problems (Wiley, New York, 1979), pp. 210-214.
[21] M. Reiser, Theory and Design of Charged Particle Beams (Wiley, New York, 1994), Chapt. 4, p. 151,
p. 196, p. 359.
[22] T.P. Wangler, RF Linear Accelerators (Wiley, New York, 1998), Chapt. 9, p. 278.

43