Bunched-Beam RMS Envelope
Simulation with Space Charge

Christopher K. AllenLos Alamos National Laboratory

1C.K. Allen

Jan 20, 2005

Abstract
We discuss an algorithm for treating the first-order effects ofspace charge in bunched beam RMS envelope simulation. The
envelope simulation is assumed to be of the type employed by
TRACE3D or TRANSPORT. Specifically, transfer matrices
are used to propagate the Courant-Snyder parameters along thebeamline. Thus, to include space charge effects we must
determine a transfer matrix representing the dynamics causedby the self electric fields. Because the space charge matrices
depend upon the Courant-Snyder parameters, and the CourantSnyder parameters in turn depend upon the space charge
matrices, we are also faced with some self-consistency issues.
We present an adaptive technique for applying the spacecharge matrices that maintains a specified level of accuracy.

2C.K. Allen

Jan 20, 2005

Outline1.

Overview
1.
2.

Motivation
Basic Approach2.

Beam Dynamics Model

3.

Self Field Representation

4.

Space Charge Algorithm

5.

Summary3

C.K. Allen

Jan 20, 2005

1. Overview
Motivation
To develop a fast simulation engine for intense charged-particle
beams systems that is first-order accurate.
Such an engine is useful for Initial machine design A (fast) online model for Model Reference Control (MRC)
applications during machine operation

4

C.K. Allen

Jan 20, 2005

1. Overview (cont.)
Approach – Extension of Linear Beam Optics


In a straightforward manner, the linear beam optics model for single particle
dynamics can be extended to the dynamics for the second moments of thebeam.

For intense beams, space charge effects are significant and must be included.
For a beam optics model, this means a matrix Φsc that accounts for spacecharge (linear force!).

For ellipsoidally symmetric beams, we can produce such a Φsc that is almost
independent of the actual beam profile.

Since the second moments depend upon Φsc and Φsc depends upon the secondmoments, we have self-consistency issues. We employ an adaptive
propagation algorithm that maintains certain level of consistency.
5

C.K. Allen

Jan 20, 2005

2. Beam Dynamics ModelOverview
Introduce homogeneous phase space coordinates zLinear beam optics - transfer matrices z(s1) = Φ(s1,s0) z(s0)Moment operator 〈⋅〉Moment matrix τ = 〈zzT〉

λPropagation of moment matrix τ(s)6

C.K. Allen

Jan 20, 2005

2. Beam Dynamics ModelPhase SpaceLet s be the path length parameter
along the design trajectory

Let z(s) denote a particle’s phase spacecoordinate in beam frame at axial
position s.

We parameterize phase space byhomogeneous coordinates in Ρ6×{1}z ≡ (x,x’,y,y’,z,z’,1)T ∈ Ρ6 × {1}

Homogenous coordinates are used because• Translation, rotation, scaling - all be performed by matrix multiplication• Allows coupling between 1st and 2nd order moments in RMS simulation
7

C.K. Allen

Jan 20, 2005

2. Beam Dynamic Model (cont.)
Linear Beam Optics Review
Assume that the particle phase z(s) can be propagated according toz(s) = ΦA(s;s0)z0

where
λ
λz0 is the initial particle phase at s = s0
ΦA(s;s0) is the transfer matrix for the beamline

Notes:λΦA’(s;s0) = A(s)ΦA(s;s0) where A(s) is matrix of applied forcesλIf A is constant, then ΦA (s;s0) = e(s-s0)AλSemigroup property ΦA(s2;s0) = ΦA(s2;s1) ΦA(s1;s0)λΦA(s;s0) is symplectic, i.e., ΦA(s;s0) ∈ Sp(6) ⊂ Ρ7_7
8

C.K. Allen

Jan 20, 2005

2. Beam Dynamics Model (cont.)
TRACE3D, TRANSPORT Type Beam Simulation
Divide beamline into N stages, one element per stage

Entrance position of stage n defined as s ≡ snDefine zn = z(sn)The action of stage n is represented by a transfer matrix Φn(un) = ΦA(sn+1,sn)
where un is the control vector for the element (magnet strengths, etc.)Then

zn+1 = Φn(un)zn

h0z0

Φ0h1
z1

Φ1

h2z2

Φ2

hN-1
z3

zN-

ΦN-1

zN

1u0

u1

u2

uN-1
9

C.K. Allen

Jan 20, 2005

2. Beam Dynamics Model (cont.)
Example – Ideal Steering Dipole
An ideal steering dipole magnet can be model with a transfer
matrix ΦSM of the form
1

0
0

Ö SM (u) ≡  0

0
0

00 0 0 0 0 0 
1 0 0 0 0 Δx′ 0 1 0 0 0 0 
0 0 1 0 0 Δy ′ 
0 0 0 1 0 0 
0 0 0 0 1 0 

0 0 0 0 0 1 

where u ≡ (Δx’, Δy’) are the dipole strengths

10

C.K. Allen

Jan 20, 2005

2. Beam Dynamics Model (cont.)
The RMS Moments
Let the beam be described by the distribution function f : Ρ6×{1} → ΡDefine the moment operator 〈⋅〉 as 〈g〉 ≡ ∫g(z)f(z)d6z

λNow define the following:ϒµ ≡ 〈z〉 ∈ Ρ6×{1} is the mean value vector (beam centroid)ϒτ ≡ 〈z zT〉 ∈ Ρ7×7 is the moment matrixô ≡ zzT x2
 xx′
 xy
=  xy ′ xz
 xz ′
 xxx′

xyxy ′xzxz ′x′ 2x′yx′y ′x′zx′z ′x′yy2yy ′yzyz ′x′y ′yy ′y′2y ′zy ′z ′x′z

yzy ′zz2zz ′x′z ′x′yz ′y ′z ′y′zz ′z ′z ′zz ′yzx x′ y 
y′ 
z 
z′ 

1 

11

C.K. Allen

Jan 20, 2005

2. Beam Dynamics Model (cont.)
RMS Envelope Transport Equations
The dynamics of τ are given byτ(s) = 〈z(s) zT(s) 〉 = 〈(ΦA(s)z0) (ΦA(s)z0)T〉 = ΦA(s)〈z0 z0T〉ΦAT(s)= ΦA(s;s0)τ0ΦAT(s;s0)
NOTES:λλλ

The moment matrices τ(s) carry the 2nd order statistics, or “RMS envelopes”
Transfer matrices ΦA(s;s0) propagate the moment matrices τ (as well as z)
Transfer equations for τ allow coupling between 1st and 2nd order moments

Space Charge:
λTo include space charge we must represent it with a transfer matrix ΦB

Require a linear representation of the self electric fields
12

C.K. Allen

Jan 20, 2005

3. Self Field RepresentationApproach

Use weighted, least-squares expansion the electric self fields
 The weight is the distribution function f itself(Expansion is then more accurate in dense regions) Self force is then represented as Fsc[z;s] = B(s)z B requires determination of field moments 〈xEx〉, 〈yEy〉, 〈zEz〉

λAssume an ellipsoidally symmetric beam in phase space
Can compute these moments in terms of τ
ϒ 〈xEx〉, 〈yEy〉, 〈zEz〉 values only weakly coupled to distribution fϒ

13

C.K. Allen

Jan 20, 2005

3. Self Field Representation (cont.)
Field Expansions
To accommodate a space charge transfer matrix ΦB, each electric self-fieldcomponent Ex, Ey, Ez must have a representation of the form
Ex = a0 + a1x + a2y + a3z
λ(e.g., for x plane)

Multiplying the above equation by the functions {1,x,y,z} then taking moments
 1

 x

 y
 z


xyx2

xyxyy2xz

yzz  a 0   E x xz  a   xE x
 1  = yz  a 2yEx
  2  a   zE x z  3  which we can solve for the ai in terms of the 〈xEx〉, 〈yEy〉, 〈zEz〉
This is the weighted, least-squares, linear approximation for the self fields
14

C.K. Allen

Jan 20, 2005

3. Self Field Representation (cont.)
The self force component Fsc[z] of the dynamics is then given by
 0
 x
 ax
 0

Fsc (z ) = B[ ô( s )]z ( s ) =  a xy
 0
 a xz 00 00 a yx

0 00 a zx

0 00 a yy
0 00 a yz0 00 00 a zy
0 00 a zz
0 00 0  x  x0 a0  x′ 0 0  y  y0 a0  y ′  0 0  z 0 a0z  z ′  0 0  1 

where the aij = aij(τ) are the expansion coefficients

The space charge transfer matrix ΦB(s) is now defined by the equation
ΦB’(s) = B(s) ΦB(s)15

C.K. Allen

Jan 20, 2005

3. Self Field Representation (cont.)
HL

Assume that the distribution f has (hyper) ellipsoidal symmetry.
-1-0.500.51x-1-0.500.51x'00.250.50.751f
Then f has the formx,x'f(z) = f[(z-µ)Tσ-1(z-µ)]
whereλλϒµ ≡ 〈z〉 is the mean-value vector (centroid location)σ ≡ 〈τ〉 - µ µT is the covariance matrix

There exists rigid motion M = RTT that aligns the ellipsoid to coordinate axes
λλT ∈ Ρ3 ⊂ Ρ7 is the translation by (〈x〉,〈y〉,〈z〉)
R ∈ SO(3) ⊂ SO(7) is the rotation which diagonalizes the matrix
 x2 − x 2 xy − x y xz − x zxy − x yy2 − y2yz − y zxz − x z 
yz − y z 2 z2 − z 
16-1-3. Self Field Representation (cont.)
In the ellipsoid coordinates (a,b,c) = φ(x,y,z), all the
cross moments are zero and the only non-zerofield moments are the followingaE aqπ a 2 b 2 c 2 2= Γ( f )a R D (c 2 , b 2 , a 2 ) ,6ε 0
N
2 2 2bEb

qπ a b c 2= Γ( f )b RD ( a 2 , c 2 , b 2 ) ,6ε 0
NcEc

qπ a 2 b 2 c 2 2= Γ( f )c RD ( a 2 , b 2 , c 2 ) ,6ε 0
NG (r ) ≡∞∫ f (s) ds ,Γ( f ) ≡ ∫ G (r ) dr ,r∞220where Γ(f) is almost constant (F. Sacherer) and RD is the elliptic integral
3 ∞dt
R D ( x, y , z ) ≡2 0 (t + x )1/ 2 (t + y )1/ 2 (t + z )3 / 2∫

In these coordinates field approximations have simply formEa ≈aE aa2a

C.K. Allen

Jan 20, 2005

3. Self Field Representation (cont.)
Summary

Aligned ellipsoid to coordinate axes with rigid motion M(τ)

Expanded the self fields as Ei = [〈xiEi〉/〈xi2〉] xiλSolved for the self moments 〈xiEi〉 in terms of matrix τλFor small changes in s self force matrix B is constantϒThen space charge matrix is ΦB(s) = M-1esBM18

C.K. Allen

Jan 20, 2005

4. Space Charge Algorithm
ApproachλλForm a transfer matrix Φ(s;s0) that includes space effects tosecond order (2nd order accurate)
Choose error tolerance ε in the solution (~ 10-5 to 10-7)
Use Φ(s;s0) to propagate τ in steps h whose length isdetermined adaptively to maintain ε19

C.K. Allen

Jan 20, 2005

4. Space Charge Algorithm (cont.)
Assume that A and B are constantϒThen full transfer matrix Φ(s;s0) = e(s-s0)(A+B)

For practical reasons, we are usually given ΦA and ΦB
λλBeam optics provides ΦA(s)In ellipsoid coordinates ΦB(s) has simply form because B2 = 0
ΦB(s) = esB = I + sB

Define Φave(s) = _[ ΦA(s)ΦB(s) + ΦB(s)ΦA(s) ]λ By Taylor expanding Φ(s) = es(A+B) we findΦ(s) = Φave(s) + O(s3)That is, Φave(s) is a second-order accurate approximation of Φ(s)20

C.K. Allen

Jan 20, 2005

4. Space Charge Algorithm (cont.)ϒ

We now have a stepping procedure which is second-orderaccurate in step length h.
Now consider the effects of “step doubling”
λλLet τ1(s+2h) denote the result of taking one step of length 2hLet τ2(s+2h) denote the result of taking twos steps of length hΔ(h) ≡ τ1(s+2h) - τ2(s+2h) = 6ch3
where the constant c = dτ(s’)/ds for some s’∈[s,s+2h]

ϒConsider ratio of |Δ| for steps of differing lengths h0 and h1
h1= h0 [|Δ1|/ |Δ0|]1/3

21

C.K. Allen

Jan 20, 2005

4. Space Charge Algorithm (cont.)
We use the formula h1= h0 [|Δ1|/ |Δ0|]1/3 as the basis for adaptive step sizing
Given |Δ1| = ε, a prescribed solution residual error we can tolerate
For each iteration kλλλ

Let |Δ0| = | τ1(sk+2h) - τ2(sk+2h) |, the residual error
Let h0 = hk be the step size at iteration kLet h1 = hk+1 be the step size at iteration k+1

hk+1 = hk [ε/| τ1(sk+2h) - τ2(sk+2h) |]1/3
where if hk+1 < hk, we must re-compute the kth step using the new steps size hk+1
to maintain the same solution accuracy

22

C.K. Allen

Jan 20, 2005

4. Space Charge Algorithm (cont.)
Example – Simulation of SNS MEBTϒ ε = 10-5
ϒ h0 = 0.03 m
SNS MEBT Simulation
7Horizontal
β (s) (m/rad)65

4

XAL hor
XAL ver32

10
0

0.5

1

1.5

2

2.5

3

3.5

4s (m)

23

C.K. Allen

Jan 20, 2005

Trace3D Verification of XAL Envelope Simulation
SNS MEBT Simulation
8

Implementing
Trace3D RDapproximation

Horizontalβ (s ) (m/rad)765

Trace3DXAL

43

2
1

00

0.5

1

1.5

2

2.5

3

3.5

4s (m)SNS MEBT Simulation
8



Numericallyevaluating RD(a la Carlson)

Horizontal β ( s ) (m/rad)765

Trace3D4XAL
32

1

0
0

0.5

1

1.5

2s (m)

2.5

3

3.5

424

C.K. Allen

Jan 20, 2005

5. Summary
 We can include the effects of space charge in a linear

beam optics model for RMS envelopes
The dynamics are accurate to first order The simulation is fast Key points (from experience)

To maintain self-consistent solution we should employ anadaptive stepping algorithm A second-order accurate stepping methods seems best
 Use accurate numerical evaluation of the elliptic integrals
 Use accurate method for determining rigid motion M25