Bunched-Beam Envelope Simulation
with Space Charge within the SAD
Environment

Christopher K. AllenLos Alamos National Laboratory

1C.K. Allen

Dec, 2005

Abstract
The capability for simulating the envelopes of three-dimensional(bunched) beams has been implemented in the SAD accelerator
modeling environment. The simulation technique itself is similar tothat of other common envelope codes, such as Trace3D,
TRANSPORT, and the XAL online model. Specifically, we follow
the second-order statistics of the beam distribution rather thantracking individual particles. If we assume that the beam maintainsellipsoidal symmetry in phase space, we can include the first ordereffects due to space charge using a semi-analytic model. This is theattractive characteristic of envelope codes, since it greatly reducescomputational time. This new feature of SAD is implemented
primarily in the SADScript interpreted language, with only a small
portion appearing as compiled code. As such, the simulation does run
slower than other compiled envelope codes such as TRACE3D orXAL, however, as interpreted code it does have the benefit of being
easily modified. We demonstrate use of the new feature and present
example simulations of the J-PARC linear accelerator section.

2

C.K. Allen

Dec, 2005

Outline1.

Overview
1.
2.

Motivation
Basic Approach2.

Envelope Dynamics Review

3.

SADScript Implementation

4.

(Field Calculations)5.

Simulation Results

6.

Issues and Conclusions
3

C.K. Allen

Dec, 2005

1. Overview
Motivation
To have envelope simulation capability for three-dimensional(bunched) beams, including space-charge, within the SAD
environment.

Such an engine is useful for\z Model reference (Fast)z Low energy electron simulationz Proton simulationz Longitudinal effects
4

C.K. Allen

Dec, 2005

1. Overview (cont.)
RMS Envelope – Approach Used Within SAD Environmentzzzz

The simulation principle is that same as that used by Trace3D andTRANSPORT. Specifically, it is an extension of linear beam optics to thesecond-order moment dynamics.
For a beam optics model we require a matrix Φsc to account for the linear partof the space-charge force, it is accurate only over short distances ∆s.
In the SAD environment we are given the full transfer matrix Φn for eachelement n. We must take the Nth root of each Φn where N = Ln/∆s is thenumber of space charge “kicks” to be applied within the element.
The space charge matrix Φsc depends upon the second moments, however, by
the dynamics equations the second moments σ depend upon Φsc. Thus, wehave self-consistency issues and must employ a propagation algorithm thatmaintains a certain level of consistency.
5

C.K. Allen

Dec, 2005

2. Envelope Dynamics Review
RMS envelope simulation is based on the following:zPhase space coordinates z = (x x’ y y’ z dp)Tz

Linear beam optics - transfer matrices zn+1 = Φn znz

Moment operator 〈⋅〉, 〈g〉 ≡ ∫g(z)f(z)d6z

z

Moment matrix σ = 〈zzT〉zPropagation of moment matrix σn+1 = ΦnσnΦnT

6

C.K. Allen

Dec, 2005

2. Envelope Dynamics Review (cont.){The moment matrices {σn} propagate down the beamline according toσn+1 = ΦnσnΦnTwhere the {Φn} are transfer matrices for each lattice elements{{To include space charge effects we must determine the self forces(from σ) then augment the dynamics σn+1 = ΦnσnΦnT accordingly.
The quantity σ ∈ ℜ6×6 is thematrix of second-order momentsof the beam distribution
given byTσ = zz

⎛ x2
⎜
⎜ xx′
⎜
⎜ xy
=⎜
⎜ xy ′⎜⎜ xz⎜⎜ xdp⎝xx′

xyxy ′xz

x′ 2x′yx′y ′x′z

x′yy2

yy ′yz

x′y ′

yy ′

y′2y ′z

x′z

yz

y ′z

z2x′dpydpy ′dp

zdp

xdp ⎞⎟x′dp ⎟⎟ydp ⎟⎟y ′dp ⎟⎟zdp ⎟⎟
2 ⎟dp
⎠7

C.K. Allen

Dec, 2005

3. SADScript Implementation
Overview
zRegister the beamline with the SAD environment using
GetMAIN[latticeFile] where latticeFile is the input deck describing
the machinez

Acquire the set of transfer matrices {Φn} and lengths {Ln} for all
beamline elements from the SAD environmentz

Take the Nth root of each transfer matrix where N is the number of
space charge “kicks” to be applied within the element. This is doneusing the matrix logarithm function.

z

Propagate moment matrix σ through each element using abovetransfer matrix and the space charge matrix Φsc computed for eachstep ∆s8C.K. Allen

Dec, 2005

3. SADScript Implementation (cont.)
Initialization{{We can obtain {Φn} and {Ln}, the lengths of the elements, from calls tothe SAD environment{Φn} = TransferMatrices/.Emittance[Matrix->True];{Ln} = LINE[“LENGTH”];
The initial moment matrix σ0 is built from the initial Twiss parameters
σ0 = CorrelationMatrix6D[{α,β,γ}x, {α,β,γ}y, {α,β,γ}z]⎛ β x ε~x⎜~
⎜ − α xε x
⎜ 0σ0 = ⎜⎜ 0
⎜⎜ 0
⎜ 0
⎝

− α x ε~xγ ε~000
β ε~0− α y ε~y− α y ε~yγ ε~000000x xy y00000y y0β z ε~z− α ε~z z⎞⎟⎟0 ⎟⎟0 ⎟⎟− α z ε~z ⎟γ z ε~z ⎟⎠009

C.K. AllenDec, 2005

3. SADScript Implementation (cont.)
Sub-Dividing Beamline Elements (the Nth root of Φn){ The transfer matrix Φn for an element n has the form
Φn = exp(LnFn)
where Ln is the length of the element and Fn is the generator matrix which
represents the external forces of element n.
{To sub-divide element n, we require the matrix Fn, given byFn = log(Φn)/Ln{

The “sub-transfer matrix” Φn(∆s) for element n can then be computed asΦn(∆s)= exp(∆sFn)

10C.K. AllenDec, 2005

3. SADScript Implementation (cont.)
Transfer Matrices with Space Charge{ Whether using the equations of motion or Hamiltonian formalism, withina section ∆s of a element n we can write the first-order continuous
dynamics asz’(s) = Fnz(s) + Fsc(σ)z(s)where the matrix Fn represents the external force of element n and Fsc(σ)
is the matrix of space charge forces.
{For Fsc(σ) constant, the solution is z(s) = exp[s(Fn + Fsc)]z0.

{Thus, the full transfer matrix including space charge should beΦn = exp[∆s(Fn+Fsc)]11

C.K. AllenDec, 2005

3. SADScript Implementation (cont.)
Numerical Efficiency{For a step size of ∆s, rather than computing the “exact” transfer matrix
Φn(∆s) = exp[∆s(Fn+Fsc)],we compute a transfer matrix which is second-order accurate in ∆s.

{Note thatΦn(∆s) = Φsc(∆s/2) Φn(∆s) Φsc(∆s/2) + O(∆s3)
whereΦn(∆s) = exp(∆sFn)(computed once per element)Φsc(∆s/2) = exp(∆s/2Fsc) = I + ∆s/2Fsc (by idempotency, Fsc2 = 0)

We have reduced a matrix exponentiation at each step ∆s to onematrix addition and two matrix multiplications
12C.K. AllenDec, 2005

3. SADScript Implementation
SADScript Modules and Some Notable Functions{

oldsad/Packages/Scheff.n{ {{sn},{γn},{σn}} = ScheffSimulate[K0, σ0, ∆s:0.01]{ {{sn},{γn},{Φn}} = GetBeamLineElementData[]{ SaveBeamMatrixData[file, {sn}, {γn}, {σn}]{oldsad/Packages/Trace3dToSad.n{ K = ComputePerveance[f, Er, W, Q]{ {α,β,γ}SAD = TraceToSadTransTwiss[{α,β,γ}T3D]{ {α,β,γ}SAD = TraceToSadLongTwiss[f, Er, W, {α,β,γ}T3D]{

oldsad/Packages/TwissUtility.n{ σ = CorrelationMatrix6D[{α,β,γ}x, {α,β,γ}y, {α,β,γ}z]{

oldsad/Packages/MatrixFunctions.n{ F = MatrixLog[Φ]{ Φ = MatrixExp[F]13C.K. Allen

Dec, 2005

4. Field Calculations{Space charge effects are included by assuming thebeam has ellipsoidal symmetry with dimensions
corresponding to the statistics in σ.

zcaf(z) = f(zTσ−1z)

byx{Analytic field expressions for such a bunch
distributions are availableqabc
φ ( x, y , z ) =
4ε 0∞

∫∞

0 x2

∫y2f ( s)z2(t + a 2 )1/ 2 (t + b 2 )1/ 2 (t + c 2 )1/ 2

dsdt++t + a 2 t +b 2 t +c 2

where a, b, c, are the semi-axes of the ellipsoid(depends upon σ) and (x,y,z) are the coordinatesalong the semi-axes
14

C.K. Allen

Dec, 2005

4. Field Calculations (cont.)
Coordinate Transformations{ To apply the previous formula for φ we must rotate to the coordinates ofthe beam ellipsoid semi-axes using a transformation
R ∈ SO(3) ⊂ SO(6).
{Moreover, we require a transformation G = diag(1,1,1,1,γ,1/γ) to convert
longitudinal coordinates from (z,dp) to (z,z’) (momentum to primed){The complete transformation isΛ = RTGTσGR
where the 〈xy〉, 〈xz〉, 〈yz〉 elements of Λ are zero

It is important that this transform is numerically accurate!
15

C.K. Allen

Dec, 2005

4. Field Calculations (cont.)
Take the Linear Part of the Electric Fields
To each electric self-field component Ex, Ey, Ez is expanded in the form
Ex = a1x + a2y + a3zz(e.g., for x plane)

Multiplying the above equation by the functions {x,y,z} thentaking moments
⎛ x2
⎜
⎜ xy
⎜
⎜⎜ xz
⎝xyy2yzxz ⎞⎟⎛ a ⎞ ⎛ xE ⎞1x ⎟⎜⎜⎟yz ⎟⎟⎜ a 2 ⎟ = ⎜ yE x ⎟⎟⎜ ⎟ ⎜z 2 ⎟⎟⎝ a3 ⎠ ⎝ zE x ⎠⎠if 〈xy〉 = 〈xz〉 = 〈yz〉 (transform GR) we have E x ( x, y, z ) =xE xx2xThis is the weighted, least-squares, linear approximation for theself fields
16

4. Field Calculations (cont.)
Field Moments
For ellipsoidal beams having a density distribution f,
the self-field moments can be calculated from φand are given by the following:yE yΛ( f )

[]

[]

[]

Q
< x 2 > R D < y 2 >, < z 2 >, < x 2 > ,3 24πε 0Λ( f ) Q=< y 2 > R D < z 2 >, < x 2 >, < y 2 > ,3 24πε 0xE x =zE z =Λ( f )zcax

Q< z 2 > R D < x 2 >, < y 2 >, < z 2 > ,3 24πε 0

where Γ(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∫

byC.K. Allen

Dec, 2005

4. Field Calculations (cont.)
Space Charge Generator Matrix
Thus, the full space charge generator matrix Fsc(σ) is given by⎛ 0⎜⎜1 / f x⎜−1 ⎜ 0Fsc (σ ) = G R
⎜ 0
⎜
⎜ 0
⎜ 0⎝0000000000000 1/ f y000000000 1/ f z0⎞⎟0⎟0 ⎟ T −1
⎟R G0⎟⎟0⎟0 ⎟⎠

[]

[]

[]

1
K222=
R
<y>
,<z>
,<x>
Df x 2 ⋅ 53 / 2

where1
K222=
R
<z>
,
<x>
,<y>
Df y 2 ⋅ 53 / 21
γ 2K
222=
R
<x>
,<y>
,<z>
D
3/ 2fz 2⋅518

5. Simulation Results
J-PARC BT Line Simulation{ Show SADScript envelope
simulation of the J-PARC
BT line for several cases
Zero current{ 30 mA{ 130 mA{

{Compare the SADScript envelope
simulation to simulations provided
by Trace3Dz

Notable differences
Trace3D does not impose
symplectic condition{ Trace3D can simulate emittance
growth thru RF Gaps (removed){

Longitudinal βz(s)12000 mArad/m

5. J-PARC Simulation

900

T3DSAD

600

300

004080120160200s (m)Horizontal βx(s)Vertical βy(s)
30

2020
rad/m

rad/m

30T3DSAD

10T3DSAD

10004080120s (m)160200004080120s (m)160200Longitudinal βz(s)
2500

2000

rad/m5. J-PARC Simulation
30 mA

1500
SAD
T3D
1000

500

004080120160200s (m)Horizontal βx(s)30Vertical βy(s)3020rad/m

20

SAD
T3D

T3D
SAD
10100004080120s (m)

160

200

04080120s (m)1602005. J-PARC Simulation

Longitudinal βz(s)

8000

6000
rad/m

130 mA

SAD
T3D

4000

2000

004080120160200s (m)rad/mHorizontal βx(s)Vertical βy(s)4040303020SAD
20
T3D

1010SAD
T3D

0004080120s (m)16020004080120s (m)160200C.K. Allen

Dec, 2005

6. Issues{ Computation of the matrix logarithm log(Φn) isexpensive.
Current procedure uses an iterative technique which
computes a matrix exponential exp(Fn) at each stepz The procedure works for Φn close to the identity matrix I{ It is not robust, but it should suffice for symplectic
matricesz{ Computation of the matrix exponential exp(Fn) is

also non-trivialz

Use a Taylor expansion with scaling and squaring.
23

C.K. Allen

Dec, 2005

6. Issues (cont.){

Currently only a simple stepping procedure is employed.
The step size ∆s remains constant throughout simulationz By implementing an adaptive stepping algorithm we can obtainsignificant speedup and maintain a specified level of accuracy inthe solution (see “Bunched Beam Envelope Simulation withSpace Charge”, KEK, Jan 20, 2005.)z{

The longitudinal space charge force seems to be slightly
stronger in all SAD simulations as compared to Trace3D.
I have debugged the SAD code extensively and have not foundany errors in the theory or implementation.
z I believe this condition is simply a result of the difference insimulation architecture, but I may be wrong.
z24

C.K. AllenDec, 2005

6. Conclusions{A major issue is computational speed
For the J-PARC simulation
Trace3D runs on the order of 0.5 secondsz SAD ScheffSimulate runs on the order of 0.5 minutesz

By implementing the adaptive stepping and/or implementing
the computationally expensive functions as compiled code weshould see significant speedup.
Only the elliptic integral function RD(x,y,z) in implemented incompiled code.
z Implementing MatrixExp() as compile code would be the mostcost-effectivez{The small difference in longitudinal dynamics between SAD
and Trace3D may be an artifact of the different approaches.
However, I am not sure.
25

C.K. AllenDec, 2005

3. SADScript Implementation{Initialization
Load the beamline information from the target “input deck”
z We acquire all the transfer matrices {Φn} and all the element
lengths {Ln}z{Space Charge
Compute the partial transfer matrix Φn(∆s) where ∆s = Ln/Nz Compute space charge transfer matrix Φsc for a distance ∆s.
z Combine Φn (∆s) and Φsc for full transfer matrix Φn,sc(∆s)
z{PropagationPropagate σn through element using σi+1 = Φn,sc(∆s)σiΦn,sc(∆s)Trecomputing Φn,sc(∆s) as necessaryz Propagate {σn} through beamline using above procedure for eachelement n.
z26

C.K. Allen

Dec, 2005

3. SADScript Implementation (cont.)
Beam Dynamics with Space Charge
To include space charge effects we must determine the self forces then augment
the dynamics σn+1 = ΦnσnΦnT to include them.
To propagate the moment matrix σn through an element we must compute atransfer matrix Φn,sc that includes space charge.
zzz

Such a transfer matrix can only be computed accurately for a short distance ∆s < Ln
We must divide each beamline element into sub-elements of length ∆s having theappropriate transfer matrix N√Φn
We apply the dynamics σn+1 = N√Φnσn N√ΦnT many times, recomputing

We use a “kick-like” approach - correcting the beam state at regular intervals
through each elementz

We divide each beamline element into sub-elements of length ∆sz

We represent the self force through ∆s as a transfer matrix Φscz

Because the self forces depend upon the beam shape, the matrix self force transfermatrix Φsc depends upon the moment matrix σ27