ADAPTIVE THREE-DIMENSIONAL RMS ENVELOPE SIMULATION IN
THE SAD ACCELERATOR MODELING ENVIRONMENT*C.K. Allen#, LANL, Los Alamos, NM 87545, U.S.A.
K. Furukawa, M. Ikegami, K. Oide, KEK, Tsukuba, Ibaraki 305-0801, Japan
Abstract
The capability for three-dimensional RMS envelope
simulation, including space charge, has been implemented
in the SAD (for Strategic Accelerator Design) [5]accelerator modeling environment used at KEK. The
dynamics within the model are similar to that used byTrace3D [3] and TRANSPORT [1]. Specifically, thematrix of all second-order beam moments is propagatedusing a linear beam optics model for the beamline.
However, the current simulation employs an adaptive
space-charge algorithm. It maintains the integration step
size as large as possible while enforcing a given error
tolerance. We concentrate on the adaptive nature of theRMS simulation, since this is the novel feature.

The complete set of phase space coordinates for a beam
particle, including both position and normalized
momentum, at axis location s is given by (x,x',y,y',z,δ) ∈ℜ6. For conciseness, it is convenient to denote thesepoints in phase space by a vector quantity z; we have(1)

Consider now a distribution of particles forming a threedimensional beam bunch and assume further that thedistribution can be described by an s-dependent functional
f : ℜ6 × ℜ+ → ℜ on phase space. Denote by 〈⋅〉 : ℜ6 → ℜ
the phase space moment operator with respect to f. That is,for any (possibly s-dependent) function g on phase spaceℜ6 we have 〈g〉 ≡ (1/Q)∫g(z;s)f(z;s)d6z where Q ≡∫f(z;s)d6z is the total charge of the bunch.
The phase space coordinates are themselves functions
on phase space and we may take their moments. Thevector z ≡ 〈 z〉 contains all the first-order moments of thedistribution f and its value is the centroid of the beam at
location s. The symmetric, positive-definite matrixτ ≡ zzT(2)

contains all the second moments of f and is known as thecorrelation matrix of the distribution. The matrix τ(s)describes the second-order evolution of f. It is ourindependent variable. Another quantity often seen in theliterature is the covariance matrix σ ≡ 〈(z− z )(z− z )T〉 = τ− zz T . It is the matrix of central, second moments of f______________________________________________* Work supported by KEK under a foreign visiting scientist grant.
#ckallen@lanl.gov

Beam Dynamics
We use a linear beam optics model for the beamline
where each beamline element n may be represented as atransfer matrix Φn. We let τn denote the correlation matrix
at the entrance to element n and τn+1 denote the correlation
matrix at its exit. Then the dynamics equations for {τn}are given by [3]τ n +1 = Φ n, sc τ n ΦTn, sc ,BACKGROUNDz ≡ ( x, x ′, y , y ′, z, δ ) ∈ ℜ6 .

and analogous to the standard deviation, its univariate
counterpart. The matrix σ describes the width of f and is
the origin of the term “RMS envelopes”. For the case of acentered beam where z = 0 , we have σ = τ, the case mostoften treated in the literature.

(3)

where Φn,sc is the transfer matrix for element n including
space charge. It has the formΦ n, sc = exp(LnG n + Ln G sc ) ,(4)where Ln is the length of element n, Gn is the generator
matrix representing external forces exerted by element n,and Gsc is the generator matrix representing the internal
forces of the beam. Equation (4) requires the assumptionthat both Gn and Gsc are constant. This is often the casefor the matrix Gn, but seldom true for the space charge
matrix Gsc; it depends strongly upon the correlation
matrix τ. Thus, if we are to employ Eq. (3) in an accurate
numerical algorithm, we are forced to step a smaller
distance h < Ln, then recompute Gsc = Gsc(τ) as necessary.

Model Quantities
Our algorithm propagates through beamline element n
in steps of varying lengths h ≤ Ln. For each beamline
element n, assume that we are given separate transfer
matrices Φn(h) for the element without space charge and atransfer matrix Φsc(h) for space charge alone, rather than
the full transfer matrix with space charge Φn,sc(h). Thus,
Φn(h) and Φsc(h) are the partial transfer matrices throughn for a distance h.
From the two separate transfer matrices Φn(h) andΦsc(h), we need to obtain the full transfer matrix Φ n,sc(h).
We must first accept that we will not get an exact valuefor Φn,sc(h). The transfer matrix Φsc(h) is computed under
the assumption that Gsc is constant over h, which it is not.
The larger the value of h, the less accurate the value ofΦsc(h). There is no point in computing Φn,sc(h) to highaccuracy, since we are always limited by the length of h.
We employ one particularly convenient fact,Φ n, sc ( h) = Φ sc ( h / 2)Φ n ( h)Φ sc ( h / 2) + O ( h3 ) .

(5)This formula gives us a transfer matrix Φn,sc(h) fromΦn(h) and Φsc(h) which is second-order accurate in h.
Equation (5) can be proven with a direct application of theZassenhaus formula for the product of matrix
exponentials [7].

Computation of the Transfer Matrices
SAD provides the transfer matrix Φn for each element n,
but not the partial transfer matrix Φn(h) that we require forspace charge effects. We must first determine thegenerator matrix Gn = (1/Ln)log(Φn) then compute Φ n(h) =exp(hGn). The matrix logarithm is computed by an
iterative technique based upon the fact that exp(−Ai)Φn →I as Ai → log(Φ n) [4]. This technique is efficient because
repeated computation of the matrix exponential is still
faster than direct computation of the matrix logarithm. It
works well for when Φ n is symplectic, but convergenceoften fails for the general case.
The process for computing the Φsc(h) is covered
elsewhere in the literature (e.g., see [3]). However, wenote some points of interest. Inferring from Eq. (4) wefind Φsc(h) = exp(hGsc). Letting Gsc0 denote the spacecharge generator matrix in the beam frame, then Gsc =L0−1R0−1Gsc0R0L0 where L0 is the Lorentz transform to therest frame and R0 is a rotation that aligns the beam
ellipsoid to the coordinate axes. In this context Gsc0 is= 0. Since
idempotent, that is, (Gsc0)2exp(L0−1R0−1Gsc0R0L0) = L0−1R0−1exp(Gsc0)R0L0, we have()Φ sc ( h) = L−01R T0 I + hG0sc R 0 L 0 ,(6)from the Taylor expansion of the matrix exponential andthe fact that R0T = R0−1. Application of (6) substantiallyreduces the computation cost of determining Φsc(h).

PROPAGATION ALGORITHM
In many existing RMS envelope simulation codes thestep size h is chosen a priori, typically as an integer
divisor of Ln. There is no way of insuring solution
accuracy in that case. In contrast we choose our step
length h dynamically, depending upon an error criterion.
The technique below picks the largest h while maintaining
solution precision. It is based upon the proceduredescribed in Press et. al [6].
Within a given beamline element n and with a given
step size h we define our stepping operator Sh : ℜ6×6 →ℜ6×6 for the correlation matrix τ. Let Φ(h) ≡Φsc(h/2)Φn(h)Φsc(h/2) and define define Sh according toS h (τ ) ≡ Φ( h) τΦT (h) .

(7)

From Eq. (3) we see that Sh steps the correlation matrix
τ(s) at a position s within element n to position s + h. Letτ(s+h) denote the exact solution of τ for an advance froms to s + h. By Eqs. (3) and (5), Sh[τ(s)] produces a secondorder approximation to τ(s+h), or more formally

S h (τ( s) ) = τ( s + h) + h3C ,(8)where C ∈ ℜ6×6 is a constant matrix given by (an abuseof) Taylor’s theorem C = τ ′′′(sˆ) / 3! for some ŝ ∈ [s,s+h].
The term h3C is the remainder for the stepping process
and our objective is to control the magnitude of this value.
To do so, let ||⋅|| be any matrix norm on ℜ6×6 then define

ε ( h) ≡ h3C = h3 C .

(9)We recognize ε(h) as the residual error in ourapproximation of τ(s+h). Assume that we are given an apriori constraint on this error, say ε . Then at each step i,our objective is to find the largest hi such that ε(hi) ≤ ε .
We accomplished this objective through step doubling.
Let τ1(s+2h) ≡ S2h(τ(s)) be the result of taking one stepof length 2h and τ2(s+2h) ≡ Sh[Sh(τ(s))] be the result oftaking two steps of length h. Then we have

τ1 ( s + 2h) = τ( s + 2h) + (2 h) 3 C ,τ 2 ( s + 2h) = τ( s + 2h) + 2( h 3C) ,(10)Let ∆(h) ≡ τ1(s+2h) − τ2(s+2h) so∆( h) = 6 h3 C = 6ε ( h) .

(11)

Then the ratio of ||∆(h)|| for two potentially differing steps
sizes hi and hi+1 is given by∆( hi ) / ∆( hi +1 ) = hi3 / hi3+1 .

(12)This relation is the foundation for generating step lengths.
Assume we are given a step hi and we wish the next step
hi+1 to maintain the error ε , that is ε ( hi +1 ) = ε . By Eq.
(11) we have ||∆∆(hi+1)|| = 6ε(hi+1) = 6ε . Substituting this
into the above and rearranging yields the desired result,[hi +1 = hi 6ε / ∆( hi )]1/ 3 .

(13)

Interpretation of formula (13) goes as follows: although
we used a length hi to step τ from s to s + 2hi, we could
have used a length of size hi+1 and still remained withinerror constraint ε . Since there is no point in recomputing
τ for the different hi+1 (we are already at location s + 2hi),
the implication is we should try a step hi+1 for the next
iteration.
There are a couple of additional points to note here
before concluding the section. The most important is that
if we find hi+1 to be less than hi, we must roll back thecomputation and re-step τ using the smaller step size hi+1.
This procedure is necessary because the condition hi+1 < hi
implies we have violated our error constraint ε(hi) ≤ ε .
Another point of practical concern is that we should notchange the step length if formula (13) suggests a verysmall change in hi+1. If hi+1 is less than hi by only a few
percent, it may not be worth the trouble of recomputing τfor a small gain in accuracy. Thus, we provide the caveat
that, given the small “slackness parameter” δh, we actually
update the next step value hi+1 only if |hi+1 − hi|/hi ≥ δh.
Finally, we note an additional procedure to thealgorithm which provides a potentially modest gain in
solution accuracy. Referring back to Eqs. (10), we cansubtract four times the second equation from the first toyield τ = (4/3)ττ2 − (1/3)ττ1 + O(h4). Use of this formula is
known as internal extrapolation. However, we cannot
monitor its accuracy. Although it is higher order, it may
not be higher accurate, we have no way of knowing. Its
use seldom does harm.
30Horizontal beta(s) rad/m

25

20

Trace3DSAD

15

10

50040

80

120

160

200s (m)Longitudinal beta(s) rad/m

2000

1500

Trace3D
SAD

1000

synchrotron. An H− beam enters the transport at 181 MeV
and 30 mA. The two plots in the figure show thehorizontal and longitudinal β functions, respectively. In
both simulations we have chosen ε = 10−5, δh = 0.05,
initial step size h0 = 3 cm, and used the l1 matrix norm.
Trace3D uses a constant step size of h = 1 cm. In Figure 1
we see very good agreement in the horizontal plane and asmall, but noticeable discrepancy in the longitudinal.
Trace3D appears to produce more space-charge effect in
the longitudinal direction. SAD uses a symplectic
technique for generating the transfer matrices whereas
Trace3D computes Φ n from the equations of motion, butthe nature of this discrepancy is still unknown.

SUMMARY
We have described an adaptive stepping algorithm forpropagating the correlation matrix τ through beamline
elements in the presence of space charge. It is based upon
the discrete transfer equations (3) for τ and the adaptive
step sizing formula of (13). At each iteration i the use offormula (13) keeps our step size hi as large as possible,
however, it does require some overhead. We must
compute three applications of Eq. (7) for a single iteration.
Yet each iterate actually propagates the τ a distance 2hiand, consequently, must be compared to two single
applications of Eq. (7). The resulting computational
overhead is 150% of a single-step algorithm, plus thecomputation of the norm ||∆∆(hi)|| and any roll-backincurred. Thus, the worst-case scenario takes at least 1.5
times the non-adaptive approach (hi constant). What we
gain from the non-adaptive approach is the guarantee ofsolution accuracy ε . Moreover, we are also guaranteed aself-consistency space-charge calculation. In most cases,
we expect a computational advantage of adaptive stepping.
Considering the overall advantages contrasted with thesmall amount of additional code development, theadaptive stepping process appears as a clear benefit inRMS envelope simulation.

500REFERENCES
004080120160200s (m)Figure 1: simulation results for SAD versus Trace3D

SIMULATION RESULTS
An RMS simulation module was developed for theSAD accelerator modeling environment based upon theprinciples here. The module was implemented primarily
in SADScript [5], a scripting language similar in syntax toMathematica [8].
Figure 1 compares the results of the SAD RMS
envelope simulation and Trace3D for the case of the Japan
Proton Accelerator Research Complex (J-PARC) at Tokai,
Japan. The beamline being modeled is the transport
section between the linear accelerator and the 3 GeV[1] K.L. Brown and R. Rothacker, SLAC-91, Rev. 3,
UC-28 (1983).
[2] D.C. Carey, The Optics of Charged Particle Beams(Harwood Academic Publishers, London, 1987).
[3] K.R. Crandall and D.P. Rusthoi, LANL Report LAUR-97-886.
[4] This technique was suggested by E. Forest.
[5] K. Oide, http://acc-physics.kek.jp/SAD/sad.html[6] W.H. Press, S.A. Teukolsky, W.T. Vetterling, and B.P.
Flannery, Numerical Recipes in C, 2nd Ed.
(Cambridge University Press, Cambridge, 1992), pp.
463-469.
[7] D.H. Sattinger and O.L. Weaver, Lie Groups andAlgebras with Applications to Physics, Geometry, andMechanics (Springer-Verlag, New York, 1986), Chap.
3, p. 36.
[8] S. Wolfram, http://www.wolfram.com.