Issue 76

Beam Dynamics Newsletter

3.18 The Ray-tracing Code Zgoubi - Latest Developments

François Méot, Brookhaven National Laboratory C-AD, Upton, NY, USA.

Introduction

The reader is referred to the ICFA-BD-NL No.43  [1], and to a more recent status article in NIM A  [2], regarding the technical details of the Zgoubi tracking engine and its accuracy, which make it “near-symplectic” - sufficiently so to model millions of turns in large hadron colliders, - and its potential for pushing particles in complex accelerator structures and accelerator optical elements or their field maps.

Extended documentation can be found in the regular FFA workshops  [3] as well as in the SourceForge code development site where in particular many examples of beam lines, FFAs and other ring simulations are maintained  [4]. There is also a comprehensive user’s guide  [5]. Further documentation, including an attempt to demonstrate 6D simulation in all possible types of particle accelerators, can be found in the “Computer Laboratory” Stony Brook University “PHY689” semester series  [6].

A recent activity and source of documentation is the project led by D. Abell and the Radiasoft Company, “New Algorithms in Zgoubi”  [7], including web interfacing  [8], as part of a US DOE computing tools development program in relating to the electron-ion collider R&D  [9].

Finally, the download tool on Zgoubi SourceForge site seems to indicate a steady increase in popularity, with over 3500 downloads from the site, covering about 60 countries, since the code was installed on the web in 2007  [10]. This may be a good indication of its efficiency in solving, at least some classes of, charged particle optics and accelerator design problems.

The latest developments in Zgoubi have found their motivation, as usual, in the requirements for new applications. This will be the focus of the present very brief state-of-the-art review, which starts in § 3.18.1 with a simulation of the 4-pass 150 MeV CBETA ERL at Cornell, of which the return loop is based on a single FFA channel that accepts four different energies  [11]. This simulation of CBETA has benefited from earlier developments in the code when simulating a 7 GeV energy-recovery experiment at the CEBAF recirculating linac  [12]; thus a few words will be said on that, § 3.18.2. As mentioned above, an important topic these days is the EIC R&D. This provides an opportunity to benchmark simulation of spin diffusion in a 21 GeV FFA ERL  [13] and to develop an ergodic approach to the polarization lifetime in an electron storage ring  [14]. This is addressed in § 3.18.3. Finally, a space charge model has been installed in the code for the purpose of studying the effect of inter-particles forces on beam optics in a scaling FFA. The interested reader is referred to related publications  [15]. These diverse accelerator design activities have fostered the development of a series of new functions in Zgoubi (“keywords”) and these are described in § 3.18.4.

Note, finally, that the Zgoubi development program now includes contributions from several others, including the aforementioned  [7], whose work can be found for instance by using the OSTI  [16] or the JACoW accelerator conference  [17] search tools. Developments include interfacing, such as the earlier web-based SIREPO, and other python based tools  [18, 19]. However this article will essentially focus on developments by the code’s author and originator.

3.18.1 CBETA

The design of the Halbach technology FFA cell  [20] for the Cornell-BNL CBETA ERL used, and was validated using Zgoubi  [21]. In the framework of CBETA commissioning, which commences in 2019, Zgoubi is one of the commissioning tools, simulating the ERL components using their field maps  [22] for accurate ray-tracing regarding an accelerator system which amounts to a time-of-flight spectrometer  [23] and as such requires accuracy on path length computation.

(image)

Figure 1: CBETA 150 MeV ERL. The linac (LA) is 36 MeV, and four different energies circulate concurrently in the single-channel return loopc(FA-TA-ZA-ZB-TB-FB): 42, 78, 114 and 150 MeV. Four spreader (SX) and recombiner lines (RX), at the linac downstream and up- stream ends, respectively, ensure proper orbits and optical functions into and out of the return loop.

(image)

Figure 2: OPERA simulation of three CBETA Halbach magnet FoDO cells. The cell includes horizontal and vertical orbit corrector dipoles (iron yoke electromagnets) on top of, respectively, the focusing and defocusing Halbach magnet (the sole QF corrector is highlighted in red here).

CBETA is a four-pass, 150 MeV, 40 mA energy recovery linac, Fig. 1. It uses a single channel loop to recirculate four energies, 42, 78, 114 and 150 MeV, four-passes up, four-passes down. The loop comprises 107 quadrupole-doublet cells, built using Halbach permanent magnet technology, Fig. 2. The arcs (FA-TA and TB-FB sections) use fixed-field alternating gradient (FFA) optics; the straight section (ZA-ZB) is a FoDO channel with all 4 energies on the common quadrupole axis. Spreader (SX) and combiner (RX) sections (4 independent beam lines, S1 to S4 and R1 to R4, respectively) connect the 36 MeV linac to the FFA arcs; they use conventional electro-magnets.

(image)

Figure 3: Field along the reference orbit across one typical dipole of the SX or RX beam lines. Three different field models are stepwise ray-traced: (i) hard edge; (ii) pure dipole field with field fall-offs; (iii) OPERA field map of the dipole. The field integral (\( \int B(s)\,\mathrm {d}s=0.044 \) T m in this particular dipole) is the same in all three models.

   

\( [R_\mathrm {i=1-4,5;j=1-4,6}] \) transport matrix from (i) hard-edge model for reference; (ii) pure dipole field with fringe fields; (iii) OPERA model:

(i) Hard-edge.   Arc length: 22.0606 cm
 1          0.216997      0          0          -3.436E-02
 0           1           0           0          -0.316659
 0           0           0.950276    0.220605    0
 0           0           -0.439585   0.950276    0
-0.31666    -0.034357         0      0           3.6085E-03
(ii) Fringe field.       Arc length: 22.0511 cm
 0.98493     0.216780       0           0          -3.433E-02
-0.13788     0.984953        0          0          -0.314304
 0           0              0.97108     0.220804   -0
 0           0             -0.258899    0.970985    0
-0.31430    -0.034326        0          0           3.7942E-03
(iii) OPERA map.        Arc length: 21.9227 cm
 0.983444       0.215709        0               0            -3.4153E-02
-0.151999       0.983583        0               0            -0.314091
 0              0               0.969607        0.218802      0
 0              0             -0.273559         0.969615      0
-0.314079     -0.0341572        0               0             3.60475E-03

Figure 3 gives an insight in the subtleties of the predictions of ray-tracing, depending on the method used to model the magnetic fields: three different models of the same magnet give three different sets of transport coefficients (only one of the models has proper field simulation, using field maps) which, cumulated over meters of beam line, result in noticeable differences in zero-th and first order optics. On top of this would be high order transport for halo simulations for instance, where accuracy on fields far off-axis is needed. Using field maps allows accurate computation of transverse motion across the Halbach cells, Fig. 4.

Figure 4 is a glimpse at the 4-pass up/4-pass down start-to-end orbit and optics simulations undertaken using Zgoubi, as part of the commissioning of CBETA.

(image)    (image)

Figure 4: Left: 42, 78, 114 and 150 MeV orbits along the re- turn loop, freely propagating through the field map string, without any correction, from their initial theoretical value at the start of the FA arc. Right: 42 MeV betatron functions along R1 after propagating freely from S1 to R1, with merely their initial values imposed at the exit of the linac; differences appear compared to when theoretical initial values (those at FB end) are imposed at start of R1: minor re-tuning however is sufficient to take care of that and get properly matched values at the linac entrance. Markers: records from field map simulations. Thin lines: from Cornell’s BMAD  [24] model for comparison.

3.18.2 ER@CEBAF

In 2015, as part of the electron-ion collider accelerator R&D at BNL  [13], a proposal was made to test high-energy, 5-pass (10 linac passes) energy recovery (ER) at CEBAF  [12], Fig. 5. This new experiment extends the 2003, 1-pass, 1 GeV CEBAF-ER demonstration  [25], and will study ER and recirculating beam dynamics in the presence of synchrotron radiation, provide an opportunity to develop and test multiple-beam diagnostic instrumentation, and can also probe BBU limitations.

(image)

Figure 5: 12 GeV CEBAF recirculator and systems to be added for a 5-pass up/5-pass down 7 GeV ER experiment (in addition to multiple-beam diagnostics).

(image)

Figure 6: The 3-bend phase chicane to be installed along Arc A. It will use main bends recuperated from Cornell’s CESR brightness upgrade.

Dipole magnets recuperated from the CESR upgrade at Cornell in 2018 have been moved to JLab. They will be used for the phase-shift chicane, Fig. 6.

As part of the preparations for the ER@CEBAF experiment, simulations including the effect of synchrotron radiation at high energy (which limits the maximum energy feasible in ER configuration to about 7 GeV) have been undertaken using Zgoubi  [26]. A sample of the findings is given in Fig. 7, a 1-pass simulation (2 linac passes, from 78.5 MeV to 1.4 GeV), which would be the starting step. Complete ER@CEBAF modeling, 5-pass up/5-pass down, is on-going.

(image)

(image)

Figure 7: Top: Optical functions in the 1.4 GeV 1-pass up, 1-pass down ER configuration, from both Elegant  [27] and Zgoubi, superimposed. Bottom: longitudinal phase-space portraits of the ER’ed bunch, observed at exit of south linac, in the presence of SR (78.43 MeV region), and without SR (78.48 MeV region).

3.18.3 Spin diffusion

Zgoubi has been able to treat spin polarization since the late 1980s  [29]. Synchrotron radiation (SR) was included inthe early 1990s  [30], and radiation in beam lines and damping in rings was thoroughly benchmarked in the late 1990s-2010  [31, 32]. However it is only in the context of the EIC R&D that the opportunity presented itself to assess the combined effect of spin diffusion.

The first combined use of spin and SR came in the eRHIC 21 GeV FFA ERL design  [13], where dedicated benchmarking showed that the numerical computation is in very good accord with theory. This is illustrated in Fig. 8  [33, Fig. 24].

(image)

Figure 8: Polarization loss by spin diffusion in eRHIC FFA ERL under the effect of SR induced energy spread  [13], including theoretical expectation  [33, § 5]. \( \sigma _\phi \) is the spin angle spreading with distance in the ERL, from theory (“theor.”) and from tracking (“track”). Horizontal axis: \( a\gamma \) is the number of spin precessions per turn, with \( a=0.00116 \), the electron anomalous magnetic moment, and \( \gamma \) the Lorentz relativistic factor.

A second combined use of spin and SR concerns a more recent study of spin diffusion in the 18 GeV electron storage ring of the eRHIC electron-ion collider. This prompted development of an ergodic approach to the computation of polarization lifetime  [28], which can be summarized as follows.

The evolution of spin polarization with time, from an initial \( \mathrm P_0 \) to an equilibrium \( \mathrm P_{eq} \) (an asymptotic quantity with time constant \( \tau _\mathrm {eq} \) to be determined from the simulations), satisfies

(1) \begin{equation} \label {EqPt} \mathrm {P(t)=P_{eq}(1-e^{-t/\tau _{eq}})+P_0\,e^{-t/\tau
_{eq}}} \end{equation}

This results from (i)  synchrotron-radiation self-polarization (index SP) with time constant \( \tau _\mathrm {SP[sec.]}\approx 99\rho ^2_\mathrm {[m]} \mathrm {R_{[m]}/E_{[GeV]}^5} \) (equating to about 30 min at eRHIC at 18 GeV), asymptotic value \( \mathrm P_{SP} \) (92.4% in a defect-free flat ring), and (ii) polarization loss by diffusion (index D), with time constant \( \rm \tau _D \), such that

(2) \begin{equation} \label {EqTauEq}\rm 1/\tau _{eq}=1/\tau _{SP}+1/\tau _{D} \end{equation}

One has in particular

(3) \begin{equation} \label {EqTauD} \rm P_{eq}=P_{SP}\times \tau _{eq}/\tau _{SP} \end{equation}

The goal in tracking spin motion is (i) to derive polarization lifetime in order to (ii) validate a ring design, including preservation of polarization under the effect of defects, corrections, spin rotators, etc., and (iii) to determine an optimal ring operating point \( \rm a\gamma _{ref} \) that maximizes \( \rm P_{eq} \). It follows that this entails maximizing the diffusion time constant \( \rm \tau _D \).

Usually the polarization lifetime is obtained by tracking a bunch of a few hundred particles over a few thousand turns at store, for each energy bin (one energy bin represents one ring set for a particular beam rigidity - or equivalently, \( \rm a\gamma _{ref} \)). The ergodic method instead tracks a single electron per energy bin, so resulting in a substantial saving in memory and CPU time. A typical outcome is displayed in Fig. 9  [28, Fig. 9]. Note the clear region of minimal spin diffusion in the graph, indicating where the storage ring rigidity setting should be chosen.

(image)   (image)

Figure 9: An ergodic approach to polarization lifetime in the 18 GeV eRHIC electron storage ring. Left, blue curve: depolarization landscape \( \rm \left .S_{y,min}(a\gamma _{ref})\right |_{t\in [0,T]} \) (the smallest value reached by the vertical spin component over the tracking interval \( \rm [0,T] \)) over \( \rm 39.8<a\gamma _{ref}<40.9 \) (200 bins of \( \Delta E=0.44 \) MeV, a single particle per bin), after \( \rm T\approx \)several thousand turns tracking. The right vertical scale (green curve) is the rms width of the energy interval explored by a particle during the tracking. Right: a zoom in on a reduced \( \rm 40.5<a\Delta \gamma _{ref}<40.8 \) interval, showing \( \rm \left .S_{y,min}(a\gamma _{ref})\right |_{t\in [0,T]} \) at (red) 9500 turns and (blue) 450,000 turns (6 s storage).

3.18.4 New keywords

The diverse accelerator design activities using Zgoubi fostered the development of a series of new functions (“keywords”): for instance GOTO, a tool for beam switching, INCLUDE (a similar effect to the Fortran or “include” statement) to concatenate s eparate optical sequence segments, SVDOC to doctor perturbed orbits, SPACECHARGE for simulation of the eponymous effect. All keywords are documented in Zgoubi User’s Guide  [5].

Three are introduced here, as an illustration in relation with the previous sections.

INCLUDE and GOTO

INCLUDE behaves in a similar way to the Fortran or include: it uploads (into the data sequence from which it is called, say a parent file “A”) a part of an external file (say file “B”). This may be convenient for repetitive structures. However the main goal is, rather, to allow the code access to B, or a part of it, properly set, for instance with optimal fields or element positioning as resulting from a prior fitting procedure.

The interest is that B can be worked on independently for optimal setting: when A is then Zgoubi’ed, it will always grab B as optimized.

GOTO can be used to signify a switch in an optical structure. It can for instance be a switch from an injection line into the main ring where multiturn can then occur, followed by a second switch into an extraction beam line, all that in a single Zgoubi run. An example of such a structure: injection+multiturn+extraction, was developed for a start-to-end simulation of the EMMA FFA  [34] (however GOTO was not available at that time). Another obvious use is in recirculating linacs, where GOTO is used to switch the beam between linac, splitters and return arc(s).

The example of CBETA is given below, featuring 4 passes up and 4 passes down, all in a single input data file, using both INCLUDE and GOTO.

In this example, the linac, return FFA loop, spreader and combiner sequences are given in separate files, uploaded by the parent file using INCLUDEs. Doing so has a series of merits :

  • - it shortens the parent input data file, so making the optical structure more apparent,

  • - it allows the linac, optical spreader, combiner and return loop sequences to be picked from files that may actually be Zgoubi’ed independently, for instance, to compute their optical parameters proper.

Labels (such as “S1_#S”, “S1_#E” below) are used for this targeted picking : they define the ends of the sequence to be INCLUDE’d, ignoring the remaining upstream and downstream parts in the file subjected to that INCLUDE (these ends may be specific, for instance, to MATRIX  [5] computations).

GOTO for its part allows throwing the beam to (and getting it back from) the linac, spreader lines, combiner lines, and recirculating loop.

!! OBJECT DEFINITION
'MCOBJET'                   ! MONTE CARLO GENERAOTR
0.14008655052543029e3           ! 42 MeV REFERENCE RIGIDITY
3
10000                       ! 10,000 PARTICLES IN THIS RUN.
2 2 2 2 1 1                                ! 6-D GAUSSIAN DISTRIBUTION
1.69E-04     -2.36E-2 0. 0. 0. 2.7144593 ! CENTER OF DISTRIBUTION
-1.01 0.58    1e-3    3         ! BEAM ALPHA_x, bETA_x
 1.80 1.09    1e-3   3          ! BEAM ALPHA_y, bETA_y
 0. 1. 0. 1                     ! BEAM ALPHA_l, bETA_l
    123456 234567 456789
 'PARTICUL'
ELECTRON


!! ALLOW RANDOMFIELD ERRORS
 'ERRORS'
0 1 123466                        dB(kG)
MULTIPOL{}    1   BP A U 0.d0     0.0   9999


!! RACETRACK SEQUENCE STARTS HERE, AT LINAC ENTRANCE
! Linac
'INCLUDE'
1
./LA/LA_obj5.inc[LA.MAR.BEG\1:LA.MAR.END\1]


!! GOTO SPREADER LINE, DEPENDING ON PASS NUMBER
 'GOTO'
PASS#
S1 S2 S3 S4 S3 S2 S1 DUMP
 'MARKER' RT_SX


!! GOTO FFAG RETURN LOOP.
 'GOTO'
PASS#
FFA FFA FFA FFA FFA FFA FFA FFA
 'MARKER' RT_FFA


!! GOTO COMBINER LINE, DEPENDING ON PASS NUMBER
 'GOTO'
PASS#
R1 R2 R3 R4 R3 R2 R1
 'MARKER' RT_RX


 'MARKER'
 'REBELOTE'     ! DO NEXT PASS: SEND ZGOUBI POINTER BACK TO
7   0.1    99   ! THE TOP (i.e., LINAC ENTRANCE), 7 TIMES.


! DUMP LINE FOLLOWS
'MARKER'    DUMP
....


 'STOP'
!/////////////////////////////
!// job done, job ends here //
!/////////////////////////////
! 4 SPREADER LINES, S1-S4
'MARKER'    S1
'INCLUDE'
1
./SX/S1.inc[S1_#S:S1_#E]
 'GOTO'
GOBACK
......... ETC
'MARKER'    S4
'INCLUDE'
1
./SX/S4.inc[S4_#S:S4_#E]
 'GOTO'
GOBACK


'MARKER'    FFA - FFAG RETURN LOOP, TAKES THE 4 ENERGIES
'INCLUDE'
1
./FFA/FFA.inc[FFA#S:FFA#E]
 'GOTO'
GOBACK


! 4 COMBINER LINES, R1-R4
'MARKER'    R1
'INCLUDE'
1
./RX/R1.inc[R1_#S:R1_#E]
 'GOTO'
GOBACK


........ ETC.


'END'

\( \bullet \) What the S1.inc spreader line INCLUDE file may look like:

'OBJET'      PARTICLE DEFINITION FOR BEAM MATRIX COMPUTATION
1000.
5.01
.001 .01     .001 .01     .001 .0001
0.     0.   0.   0.   0. 1.
-0.975 13.64 -0.975 13.64 0. 1. 0. 0. 0. 0. ! OPTICAL FUNCTIONS AT LINAC


'MARKER' S1_#S            ! INCLUDE GRABS FROM HERE
'INCLUDE'
1
S1_sEQUENCE[LA.MATCH1:S1.MAR.END]
'MARKER' S1_#E            ! INCLUDE GRAB ENDS HER


'FIT'
 ... A set of variables and constraints
 ... which match S1 to FFA loop entrance beta functions.
'MATRIX'     ! DELIVERS TRANSPORT AND BEAM MATRICES, FOR CHECKS.
1 0
'END'
SVDOC

Orbit defects and correction, for the vertical orbit in particular, are paramount in the study of polarization transport and lifetime. Accounting for that requires a random generator (Zgoubi’s “ERRORS” available for that generates all sorts of random field and positioning errors  [5]) and a correction algorithm. This is the role of the new “SVDOC”.

Figure 10 shows an application in the case of the 3.8 km long RHIC polarized proton collider  [35]. SVDOC offers also the possibility of generating, in a single run, an arbitrary number of corrected random orbits, in as many Zgoubi input files ready to run, this is for the purpose of possible statistics.

(image)

Figure 10: Orbit correction in the RHIC lattice, using SVDOC. Empty markers represent the initial random H (squares) and V (circles) orbits, measured at RHIC pickups; the solid markers represent the corrected H and V orbits.

   

\( \bullet \) What SVDOC looks like in a Zgoubi input data sequence  [5]. Note that it is preceded by a “FIT” of which the task, as part of the algorithm, is to find the closed orbit coordinates.

'ERRORS'                      ! GENERATE RANDOM ERRORS. "MULTOPOL"s ARE CONCERNED.
0    2 123456                 ! RANDOM SEED.
MULTIPOL{QF} 1 BP A U 0. 2e-2 3         ! FIELD IN ALL "QF" (PARAMETER "12" IN SVDOC).
MULTIPOL{QD} 1 BP A U 0. 2e-2 3         ! FIELD IN ALL "QD" (PARAMETER "22" IN SVDOC).


 'FAISCEAU'
 'FIT'                  ! CLOSED ORBIT FINDING.
 4                      ! VARIABLES: VARY INITIAL H AND V COORDINATES
Vary Y0; Vary T0; Vary Z0; Vary P0
 4                      ! CONSTARINTS: INITIAL AND FINAL COORDINATES IDENTICAL
Yfinal = Y0; Tfinal = T0; Zfinal = Z0; Pfinal = P0


'SVDOC'                 ! ALGORITHM THAT COMPUTES THE SVD MATRIX AND APPLIES IT
1    0                  ! ON 50 (SEE BELOW) RANDOM DEFECT ORBITS.
PUH{HMON} PUV{VMON}         ! NAMES OF H-, V-, HV-PICKUP FAMILIES, AND
CRH{HKIC} CRV{VKIC}         ! NAMES OF H- & V-CORR FAMILIES, TO GENERATE SVD MATRIX.
1.e-3           1.e-3       ! KICK VALUES FOR COMPUTATION OF SVD MATRIX.
1    2   50     ! "1" TRIGER "ERRORS" KEYWORD; 2 CHANGES PRIOR TO THAT; 50 ORBITS.
ERRORS     12 ! "12": VALUE OF HKIC IN "ERRORS" (it will deteemine the H orbit),
ERRORS     22 ! "22": VALUE OF VKICK IN "ERRORS" (will determine the V orbit).


'SYSTEM'                ! SYSTEM CALL WHICH DELIVERS THE FIGURE ON THE LEFT,
1                       ! PRIOR TO TERMINATING ZGOUBI EXECUTION.
gnuplot <./gnuplot_SVDOrbits.cmd &
'END'
References