Frequency analysis and dynamic aperture studies in ELENA with realistic 3D magnetic fields
aa r X i v : . [ phy s i c s . acc - ph ] A ug Frequency analysis and dynamic aperture studies in ELENA with realistic 3Dmagnetic fields
Lajos Bojt´ar ∗ CERN, CH-1211 Gen`eve 23 (Dated: September 1, 2020)We briefly present recent progress with our algorithm and its implementation called SIMPAdescribed in a previous paper [1]. The algorithm has a new and unique approach to long-term4D tracking of charged particles in arbitrary static electromagnetic fields. Using the improvementsdescribed in this paper, we made frequency analysis and dynamic aperture studies in ELENA. Theeffect of the end fields and the perturbation introduced by the magnetic system of the electroncooler on dynamic aperture is shown. A special feature of this study is that we have not introducedany multipole errors into the model. The dynamic aperture calculated in this paper is the directconsequence of the geometry of the magnetic elements. Based on the results, we make a fewsuggestions to reduce the losses during the deceleration of the beam.
I. INTRODUCTION
The importance of fringe fields in small rings is wellknown and it has been taken into account for multipolemagnets at various degrees for decades [2–5]. The ap-proach to particle tracking we described in [1] naturallyincludes the end fields altogether and for all kinds of ele-ments with the same treatment. The aim of this paper istwo-fold. One is to apply our algorithm described in [1]on the Extra Low ENergy Antiproton (ELENA) ring [6],the 30.4 meters circumference machine built at CERNto decelerate antiprotons. Another aim is to gain insightinto the beam dynamics of ELENA by frequency analysisand dynamics aperture studies.Long-term tracking of charged particles is a fundamen-tal problem of accelerator physics, plasma physics and itis important in astrophysics. In a recent paper [1] we de-scribed a new algorithm allowing long-term symplecticintegration of charged particle trajectories in arbitrarystatic magnetic and electric fields. We introduced sev-eral new ideas and pointed out the importance of tworequirements for a physically valid description of the elec-tromagnetic fields.One important point is that the potentials describingstatic magnetic and electric fields must obey the Laplaceequation everywhere in the domain of interest. This isa necessary condition to get physically valid fields fromthe potentials. Failing to satisfy this condition leads tospurious sources [7] which in turn leads to an energy driftduring the tracking. All 3D grid methods suffer from thisproblem because there is no known interpolation satisfy-ing the Laplace equation while being continuous betweenthe 3D grid cells.In accelerator physics, transfer map methods [8–10]are also applicable to handle realistic 3D electromagneticfields. These are conceptually rather different and theyassume the existence of a closed orbit, which puts re- ∗ [email protected] strictions to the range of problems it can handle. FFAs,stellators, tokamaks for example are difficult to handlewith transfer maps, but our algorithm can treat thosecases without complications.Another approach to deal with the problem of the rep-resentation of 3D fields is to take the Fourier transformof the field in three variables [11–14]. One disadvantageof this method is the high number of harmonics neededwhen the field has detailed local features.The second requirement for a physically valid descrip-tion of the field is the continuity of the potentials. Cut-ting the potentials at some distance, even far from a mag-netic element, leads to a systematic error which accumu-lates during the tracking leading to an energy drift andnon-preservation of the phase space volume. This hasbeen recognized also by others [3, 5, 15] and mitigationshas been used to handle this effect. Our method avoidscutting the potentials by treating the entire volume ofinterest as a whole.Before this study, little was known about the effect ofthe magnetic system of the electron cooler on the ELENAbeam dynamics. More specifically, the effect of the toroidcoils in the cooler was mostly unknown, not due to lack ofeffort, but due to the lack of appropriate tracking tools.A considerable amount of work went into the develop-ment of a tracking tool by colleges at CERN to addressthis issue [16]. The tool was based on a field map givenon a 3D grid. It produced manifestly non-symplectic be-haviour already in a few hundred turns.This triggered the author to devise the algorithm [1]and implement it. The name of the software is SIMPA, anabreviation of symplectic integration through monopolearrangements. At present it is a Java library. A lot ofwork remains to be done before it can be released for gen-eral use. It is intended to become a tool to produce fieldmaps suitable for long-term, charged particle tracking ingeneral, not only in the domain of accelerator physics. Itcontains code to use the generated field maps in actualtrackings. Recent progress [17] made possible to imple-ment an explicit second order symplectic integrator fornon-separable Hamiltonians, which arise when magneticfields are involved.The simulation tools needed for frequency analysis anddynamics aperture studies described in this paper werewritten in Java using the SIMPA library. Both the simu-lation source code and the library source code is availablefrom the author upon request. The SIMPA source codewill be released under an open source license in the fu-ture. People interested in the code are encouraged tocontact the author.The algorithm described in [1] was applied to a sim-ple test machine with a beam region having a constantcircular profile. This was covered with a single line ofoverlapping spheres. See Fig. 4. in [1] for an illustration.Actual machines like ELENA have beam regions withmore complex geometries. Much progress has been madeto extend the algorithm to general volumes. Full descrip-tion of the enhanced algorithm is outside the scope of thiswork and will be described separately, only a summaryfollows.The paper is organized in the following way. Wepresent the progress made to cover general volumes andexplain also how this leads to gain in speed and accu-racy. Then a number of basic checks are presented, com-paring the model of ELENA built with SIMPA to theMAD-X [18] model. In the frequency analysis sectionthe effect of fringe fields and the magnetic system of theelectron cooler is investigated. The dynamic aperturesection quantifies the effect of resonance lines identifiedas relevant by the frequency analysis. In the conclusionswe suggest a few possible improvements for the perfor-mance of ELENA.
II. THE TRACKING ALGORITHMA. A short summary
We recommend reading the previous paper [1] to un-derstand the algorithm in detail, as only a summary isprovided here.Symplectic integrators keep the error resulting fromthe integration itself bounded, but can not cure the er-rors coming from the representation of the fields. Theseare two separate sources of errors. It is crucial to havea physically valid representation of the fields. When anelectromagnetic field not satisfying the Maxwell equa-tions is given to a symplectic integrator the integrationwill be still symplectic, but the result will describe some-thing other than the true physical behavior of the system.A common example is the energy drift due to spurioussources appearing in 3D grid methods [7].The potentials describing the fields must satisfy theLaplace equation and should be continuous everywherein the beam region to obtain accurate results in long-termtrackings. To comply with these conditions, the poten-tials are expressed analytically in terms of their sources.We use hypothetical magnetic and electric monopolesand current loops as sources. These sources are placed outside of the volume of interest, at some distance fromthe boundary, and their strength is set such that theyreproduce the normal component of the magnetic or elec-tric field at the boundary by solving a system of linearequations.The normal component of a magnetic or electric fielddetermines completely the field inside a closed simply-connected volume. For a multiply-connected volume, thisis not always true. If the integral of the magnetic fieldalong a closed curve is not zero, then magnetic monopolesalone can not reproduce the magnetic field. One exampleof this situation is the solenoid magnet. In these cases,the non-conservative part of the vector potential mustbe provided by one or more current loops. The wordnon-conservative should be understood here in a mathe-matical sense, the vector field having a non-zero integralalong a closed loop , that is H v · d r = 0.After the potentials are reproduced at the boundary bythe sources, they can be evaluated anywhere inside thevolume analytically. The potentials satisfy the Laplaceequation and they are continuous. However, this methodis too slow to be practical. We will refer to this part ofthe algorithm as slow evaluation. We need a much fastercalculation of the potentials to be of practical use.Several orders of magnitude improvement can beachieved by using a local description of the potentials.Spherical harmonics scaled appropriately are called solidharmonics. Regular solid harmonics are the canonicalrepresentation for harmonic functions inside a sphere[19]. A key characteristic of the algorithm is the descrip-tion of vector and scalar potentials by solid harmonicsinside a set of overlapping spheres covering the volumeof interest.The potentials satisfy exactly the Laplace equation in-side the spheres. The discontinuity between the spheresdecrease exponentially with the degree of solid harmon-ics expansion and can be easily kept close to machineprecision. The representation of the potentials in termsof solid harmonics is optimal in terms of memory andallows fast evaluation. We will refer to this part of thealgorithm as fast evaluation. The speed of the evaluationis strongly dependent on the degree of the solid harmon-ics expansion. This property is the base for the speedimprovements we will describe later.One might ask what is the point of having a slow eval-uation when we have a fast one. To calculate the solidharmonics coefficients for each ball in the fast evaluationwe need the values of the scalar and vector potentials.These potentials must be continuous and must satisfythe Laplace equation. This is the role of the slow eval-uation. CAD software or measurements can not providethe potentials satisfying these conditions. B. Extension to general domains
In [1] we show cased the algorithm on a simple accel-erator. The beam region had a circular shape everywherewith a constant 3 cm radius all along the reference orbit.This was covered by a single row of overlapping sphereswith a 4 cm radius. This solution to cover the beamregion is sufficient to demonstrate the method, but realaccelerators have a more complex aperture in general.Apart from that, our algorithm aims to apply to a widerange of problems, not only in the domain of particleaccelerators.The usual way to describe a general volume is to giveits boundary surface by triangulation. A common for-mat to describe a triangulated surface is the StandardTessellation Language abbreviated as STL. The bound-ary of the volume of interest is described by an STL filein SIMPA. There are many CAD software able to pro-duce STL files, however, to construct the aperture foran accelerator using CAD software is not a trivial task.There is a utility implemented in SIMPA facilitating theconstruction of the STL file. It takes the design orbit anda list of transfer profiles as the input and extrudes thoseprofiles along the design orbit. For each transfer profile,the user must specify the longitudinal position where theextrusion starts and ends along the orbit. The transverseprofiles are described as polygons. The number of seg-ments for each polygon is restricted to the same value,making the implementation of the triangulation straight-forward.Once the boundary of the beam region is specified byan STL file the next step is to fill it with overlappingspheres, such that the entire beam region is covered with-out gaps. There are constraints concerning the size of thespheres. The slow evaluation calculates the potentialsfrom sources, which are placed at some distance fromthe boundary. The spheres should be small enough to notoverlap with the sources. The solid harmonics expansionworks only if there are no sources inside the spheres.It is practical to place the center of the spheres ontoa regular lattice. This simplifies the construction of thecover as well as finding the spheres when the potentialsare evaluated. There are many possible lattices, but onlythe hexagonal close-packed ( HCP ) and the face-centeredcubic (FCC) lattices have the optimal packing density,about 0.74. We chose the HCP lattice. The coordinatesof the sphere centers on an infinite HCP lattice can beobtained by three simple expressions. We keep only thosespheres from the infinite lattice which are necessary tocover the beam region.The first step is to find the indices of one initial sphereon the lattice, which has its center inside the beam re-gion. Then a 3D version of an algorithm from computergraphics, called flood-fill, is applied to collect all spherecenters inside the boundary. In essence, the flood-fill al-gorithm checks neighbouring lattice locations recursively,if they are inside or outside the boundary. After this step,we have a list of spheres which covers most of the beamregion, but not all. Near the boundary, there can begaps.A second step is necessary to cover the beam regionfully. In the second phase, all spheres kept in the first
FIG. 1. Triangulated boundary of the ELENA beam region.The interior plot is an expanded 3D view of an arc. The reddots are the centers of the covering spheres located on an HCPlattice. step are checked in the following way. If the sphere inquestion has missing neighbours in the lattice, then wetake the half-line from the sphere center to the centerof the missing neighbour and calculate its intersectionpoint with the sphere in question. If this point is in-side the boundary of the beam region, then we add themissing sphere to the cover. We do this for all missingneighbours, hence eliminating gaps. We should note thissecond step can be applied only if the distance of sourcesproviding the potentials are further from the boundarythan the diameter of the covering spheres. In the presentstudy that is the case. Otherwise, some spheres wouldhave sources inside, which is not allowed. A more so-phisticated second step of the covering process could beapplied by removing the restriction to put sphere centersonly to the lattice points and allowing variable spheresizes in the second step. This is a non-trivial problem incomputational geometry, but certainly has a satisfactorysolution. FIG. 1 shows the ELENA beam region and thecenters of the covering spheres. The radius of the sphereswas 1 cm in this study and 75951 of them were neededto cover the aperture.
C. Improvements in speed and accuracy
Covering the volume of interest as described abovemakes the method applicable to a much wider range ofproblems, but this is not the only advantage. In general,the error of the spherical harmonic approximation on thesurface of a sphere is less than the first term omitted fromthe infinite sum of spherical harmonics times a constant[19]. Let us inspect Eq. (9) in [1], which is restated herefor convenience. f ( r, θ, φ ) = 1 R ℓ max X ℓ =0 ℓ X m = − ℓ r ℓ c ℓm ˜ P mℓ (cos θ ) Φ( φ ; m ) . (1)In this equation the value of r is between zero and R ,the radius of the sphere. We can set r = R and calculatethe maximum error of the solid harmonics approximationas ǫR ℓ max , where ǫ is the error due to the first omittedterm in the sum with degree ℓ max + 1. Because of thepresence of the term r ℓ , the biggest possible error in theapproximation strongly depends both on R and ℓ max .For instance, in [1], we used a single row cover with R = 4 cm and ℓ max = 38. In this study R = 1 cm. Ifwe had used the same ℓ max = 38 with R = 1 cm spheres,the maximum error of the approximation ǫR ℓ max wouldhave been a factor 4 smaller, assuming infinite precisionarithmetic and the two spheres having the same center.This is a big number, so we can reduce ℓ max and the errorof the solid harmonics approximation still stays close tomachine precision.In this study we used ℓ max = 12 uniformly for allspheres. This speeds up the evaluation of the potentialsand their derivatives significantly, because the numberof operations to calculate the sum in Eq. (1) is pro-portional to ℓ . On Intel i7-6700 CPU at 3.40 GHz,using a single core with hyper-threading enabled, we got2 . × evaluations per second of the vector potentialand its derivatives with ℓ max = 12. This should be com-pared against 5 . × evaluations with the single rowcover and ℓ max = 38, a factor 4.66 improvement.Further speed improvement can be achieved by chang-ing ℓ max between the spheres to keep the error of theapproximation below a prescribed value. There is alsoample space for optimization by improving memory lo-cality, better use of vector instructions of the CPU, andparallelization. This will be implemented in the future.There is one downside of the reduced sphere diameter,that is the increased memory use. The number of spheresneeded to cover a given volume is proportional to 1/ R .The number of coefficients in a solid harmonics approx-imation is proportional to ℓ , which in turn dependson R . Together, the number of coefficients needed to de-scribe a potential in a given volume with some prescribedaccuracy is n = c ℓ /R , where c is a constant. Theconstant c depends on the volume to be covered. In thepresent study the size of the field map is 448 MB. Us-ing a single row cover with ℓ max = 38 it is only 29 MB.Memory is traded for speed. III. PREPARATIONSA. Preparing the field maps
As a first step, the field of each type of magnet inELENA has to be expressed as a collection of magneticmonopoles or a combination of magnetic monopoles and current loops. This provides continuous and analytic po-tential everywhere in the beam region. To do so, the mag-netic field values have been obtained from the CAD soft-ware OPERA for each magnet at specific points on a sur-face surrounding the beam region. This surface is closeto the poles of the magnet. The magnetic monopoles areplaced outside of the surface at 2.5 cm distance. Then asystem of linear equation was solved to find the strengthsof the monopoles for each magnet type. The relative pre-cision of the reproduction of the magnetic field from thevector potential is typically between 10 − and 10 − .Once each type of magnets in the machine has beenexpressed as a collection of monopoles and current loops,the next step is to assemble the ELENA ring from thesecollections. This has been done by rotating and translat-ing the collections of sources to the correct place. Theentry and exit coordinates of the magnets have been ob-tained from the MAD-X model with the survey com-mand.After the machine is assembled, a particle can betracked to check if the procedure used so far has beencorrect. This tracking is slow because the vector poten-tial of all sources has to be evaluated at each time step.Nevertheless, it is fast enough to track a few 10’s of turnsto calculate the closed orbit and set the values of the orbitcorrectors, which is necessary when the electron cooler isincluded.The magnets of the ELENA ring have been organizedinto magnet groups. Usually, a magnet group consists ofmagnets connected to the same power supply. For eachmagnet group, a field map was produced. These fieldmaps consist of a large number of overlapping spheres,each with some number of coefficients determined bythe maximum degree of solid harmonic expansion as de-scribed earlier. The location and size of the spheres arecalculated once and it is common for all field maps. Thismakes it easy to scale field maps individually and sumthem into a single map before tracking. B. Basic checks
It is reassuring to implement some checks and comparethe results to an established tool when this is possible.Since there is already a MADX-PTC model available forELENA, we use this model as a reference for compari-son. There are several global parameters of the machinewhich can serve as cross-checks, such as the tunes andthe optical functions. We compared the optical β anddispersion functions to those obtained with MAD-X atthe same currents for the quadrupoles. In this compar-ison, the electron cooler magnets and the solenoid com-pensators were turned off in both models as well as allorbit corrector magnets.In our model, the optical functions are obtained fromtracking. This might look more complicated than theusual linear optics model, but it is necessary because wedeal with arbitrary 3D fields. To obtain the Twiss func- b h b v D h b h , b v , D h [ m ] s [m] MAD−X SIMPA FIG. 2. Comparison of the Twiss functions β h , β v and D h calculated by SIMPA and MAD-X at the working point Q h = 2 .
454 and Q v = 1 . tions, we tracked a p with nominal energy while savingthe phase space variables at each time step. The sym-plectic integration method we implemented according to[17] uses canonical variables, but these were convertedto positions and angles. The procedure results in manyPoincar´e maps along the orbit. We then fitted an el-lipse to the phase space variables at each longitudinalposition. From the parameters of the ellipse, the Twissfunctions were calculated by well-known formulas. To getthe dispersion, another particle with a different dp/p wastracked and the center of the fitted ellipses was comparedto the dp/p = 0 case. The β and D functions are plottedin FIG. 2. The agreement with MAD-X is remarkable ifone considers the completely different calculation of theTwiss functions.The tunes were also compared. We chose the workingpoint Q h = 2 .
454 and Q v = 1 . Q h = 2 . Q v = 1 . p with dp/p = 0 .
001 was tracked and compared tothe p with nominal energy to calculate the linear chro-maticity. The tunes for both p were calculated from thetracking results with a very precise frequency analysistool called PyNAFF [21]. The chromaticities we gotare ξ h = − .
72 and ξ v = − .
15, which is significantlydifferent from the results obtained with MADX-PTC, ξ h = − . ξ v = − .
81. One possible explanationfor this discrepancy is that the end fields in MADX-PTCare taken into account by the FINT and HGAP parame-ters [18, 22] for the bending magnets, while our algorithmhandles the end fields nearly exactly for any element. Wehave not made measurements during commissioning withsimilar conditions in ELENA.
IV. FREQUENCY ANALYSISA. Motivation
Small machines have special challenges, which can beneglected in large rings. One of these is the increased in-fluence of fringe fields on the beam dynamics. ELENA isa ring with 30.4 meters circumference. It is mandatory totake into account the effect of fringe fields already duringthe design of machines of this size. The ELENA MADX-PTC model includes the fringe fields for the main bendingmagnets, which is one of the main sources of nonlinear-ities. Another important source of nonlinearities is themagnetic system of the electron cooler, but it is not in-cluded in the MADX-PTC model. Our model includesall fringe fields and also the magnetic field of the electroncooler with high fidelity.
B. The method
According to the KAM theory, invariant tori of a con-servative dynamical system are stable under small per-turbations. These tori describe quasi-periodic orbits inphase space. The motion of the system can be describedby action-angle variables. On invariant tori the actionvariables are constant and the angle variables have fixedirrational and Diophantine frequencies. When the fre-quencies are rational numbers, they are related to eachother by integers and resonance occurs. The conditionfor resonance in a particle accelerator in 4D phase spaceis given as nν x + mν y = p (2),where ν x , ν y are the horizontal and vertical tunes and n, m, p are integers, not all zero. A more detailed ac-cessible explanation can be found in [23] for instance.Although the frequencies where a resonance occurs areknown a priory, their strength or significance is notknown. One can track particles to investigate thestrength of the resonances. Doing so for a large num-ber of turns is time-consuming. Fortunately, there aresome early indicators to detect chaotic behavior. Oneearly indicator is the frequency shift.An early use of frequency analysis for conservativeHamiltonian systems by Laskar [24] was to investigatethe stability of the solar system. Later it was applied toother domains among them accelerator physics [25, 26].For a comprehensive description of the method see [27].The precision of Laskar’s NAFF [28] algorithm is the fea-ture that makes it interesting. While the frequency errorof the standard FFT is proportional to 1 /N , the error ofNAFF is proportional to 1 /N , where N is the number ofsamples. This allows us to determine the characteristicfrequencies of the motion with a small number of turns.The first characteristic frequencies for each plane ofthe motion are the tunes. The drift in the tunes canserve as an early indicator of the long-term stability ofthe motion [29]. For initial conditions corresponding tochaotic trajectories, the frequency can only be defined fora given time interval. Calculating this local frequency fortwo consecutive time intervals and taking their differencewe can calculate tune shifts ∆ ν x,y . For a coasting beam,the figure of merit can be defined as: D = log ( q ∆ ν x + ∆ ν y ) (3) D calculated with this definition can be used to identifystable areas in the tune diagram.To identify promising working points, we scanned alarge range of tunes in the horizontal and vertical planesand for each initial condition the figure of merit D wascalculated. The ELENA machine has three quadrupolefamilies, which give three parameters for the tune scan.There is also a constraint for the dispersion in the elec-tron cooler region, which is fixed to 1 meter. This con-straint selects a 2D plane in the 3D space of possiblequadrupole currents. The two-dimensional tune spacewas scanned between 2 . − .
52 and 1 . − . π mm mrad ] in both planes. Thisparticular value for the emittance was chosen to be largeenough to see clearly the non-linear effects due to thefringe fields and the electron cooler, but without hav-ing too much losses due to the physical aperture of themachine. In this paper emittances and acceptances aredefined as the area of the ellipse corresponding to thesingle particle Courant-Snyder invariant, unless statedotherwise.The tunes were scanned in 160 steps in both directionsgiving 25600 initial conditions. Each particle was trackedfor 300 turns and the phase space variables were saved ata single longitudinal position for each turn. The trackingresult was post-processed with PyNAFF, a Python im-plementation of the NAFF algorithm [28]. To calculate∆ ν x,y , we split the 300 turns into two sets. Then the fig-ure of merit D in Eq. (3) was calculated for each pointin the diagram and plotted as colors. C. Results
FIG. 3 shows the figure of merit D for the bare ELENAmachine with bendings and quadrupoles only, plottedagainst the estimated tunes as colors. The term ‘esti-mated tune’ means the input for the polynomial, map-ping the desired tunes to quadrupole currents, which isslightly different from the actual tunes. The red dotsindicate the initial conditions which resulted in particleloss in less than 300 turns.Instead of the estimated tunes, we can plot D againstthe tunes calculated by pyNAFF as shown in FIG. 4.This gives a precise value for the tunes, but this tune cannot be calculated for the lost particles, that is why FIG. 4has no red dots inside. FIG. 3 and FIG. 4 are somewhatcomplementary.On both plots, many higher-order resonance lines canbe seen. It should be emphasized, these resonance linesare present also in a perfectly manufactured, built, andaligned machine. We have not put any imperfections intothe model. All the resonance lines are the results of thefringe fields of the bending magnets and the quadrupoles.They are the direct consequence of the geometry of themagnets.The next step was to repeat the exercise with the mag-netic elements of the electron cooler included. Those arethe main solenoid coil, two toroid coils at the ends, twosolenoid compensator magnets, and several orbit correc-tors including the two ‘Kyoto type’. We did all trackingsat the lowest energy of 100 KeV, because the effect of theelectron cooler is the biggest there. The cooler has fixedcurrents for all coils during the deceleration cycle.The step size was set to 1 cm while tracking the baremachine. With the electron cooler included, it was re-duced to 0.5 cm. This was needed because the coolerintroduced some spatially fast-changing magnetic fields.FIG. 5 shows D against the estimated tunes and FIG.6 against the numerically measured tunes. It is imme-diately apparent that the magnets of the electron coolerhave a significant effect on the beam dynamics. Manyresonance lines became stronger and wider and particlelosses are more frequent at the strongest resonance lines.It is natural to ask, how much space is available for thebeam in these diagrams. To answer this question we sam-pled 1000 particles from a Gaussian distribution with 2 σ RMS emittances 75 [ π mm mrad ] in both planes, whichis the acceptance of the machine and dp/p = ± − which is about the momentum acceptance given as 2 σ RMS. These particles were tracked with the machine setto the working point q h = 2 .
455 and q v = 1 . ξ h = − . ξ v = − .
55 . Another contribution comes from the be-
FIG. 3. D plotted as colors against the estimated tunes for thebare ELENA machine, consisting of only the main bendingsand the three quadrupole families. The red dots correspondto lost particles.FIG. 4. D plotted as colors against the numerically measuredtunes for the bare ELENA machine, consisting of only themain bendings and the three quadrupole families. tatron amplitude dependence of the tunes. The thirdand probaly the hardest contribution to deal with is thespace charge. The effect of space charge is not inves-tigated in this work, but is clearly an important factorto consider. According to [6], the incoherent tune shiftcan be -0.1 for a short time before extraction when thebeam is bunched. For a coasting beam it is much less,about -0.01, but for a much longer time in the order ofseconds. Even in the case of coasting beam the contri- FIG. 5. D plotted as colors against the estimated tunes forthe ELENA machine with electron cooler. The red dots cor-respond to lost particles.FIG. 6. D plotted as colors against the numerically measuredtunes for the ELENA machine with electron cooler. The ma-genta dots in this plot correspond to the tunes of particleswhich survived 10 turns in a bunch filling the acceptance,without space charge included in the model. bution of the space charge to the tune footprint is sig-nificant. Therefore it is important to choose a workingpoint which have a sufficiently large space around it inthe tune diagram without strong resonance lines, whichis the working point q h = 2 .
455 and q v = 1 . V. DYNAMIC APERTURE
Frequency analysis allowed to compare and select thepromising working points by giving estimates of the rel-ative strength of resonances in the tune diagram. Thenext step is to analyze the effect of these resonances onthe beam. The concept of dynamic aperture quantifiesthe available stable phase space volume for the beam.
A. The method
The stability domain is a connected region in phasespace of initial conditions that remains bounded after N turns [30].The stability region is the area inside the last simply-connected invariant curve in phase space. In 4D, theinvariant curves do not always separate a connected re-gion, but in most practical cases it exists. We will seelater in the results that even if the working point is on aresonance line a connected region can be identified.In 4D, the volume of the stability domain can be cal-culated by the integral Z Z Z Z χ ( x, p x , y, p y ) dx dp x dy dp y , (4)where χ ( x, p x , y, p y ) is a function of initial conditions,with value one if the particle has survived N turns, orzero if it was lost.The dynamic aperture is defined as the radius of thehyper sphere having the same volume in phase space asthe stability domain. In practical calculations the phasespace coordinates ( x, p x , y, p y ) of Eq. (4) are expressed interms of polar coordinates α, r in physical space and thephase space angles υ , υ , then the integral approximatedas a sum over these coordinates.The direct calculation of Eq. (4) by summation over allfour variables is very CPU time demanding. Instead weapplied the method described in [30, 31] as ‘integrationover the dynamics’. The idea is to replace the integrationover the phase space angles υ , υ by averaging over thedynamics. While particles take turns, they sweep the fullrange of phase space angles υ , υ , and this informationis available.The procedure is the following. With fixed initial phasespace angles ¯ υ , ¯ υ , the polar angle α k is scanned in K steps. For each α k the radius of the last stable orbit r ( α k , ¯ υ , ¯ υ ) is sought by increasing the polar radius vari-able in J steps. Then the N turns of the last stable orbitfor each α k is distributed into an angular grid of size M × M , such that each grid element should contain atleast one iterate from the orbit of r ( α k , ¯ υ , ¯ υ ). Then thedynamic aperture is calculated with the formula r d = " π KM M X m ,m K X k =1 r m ,m ( α k , ¯ υ , ¯ υ ) sin (2 α k ) , (5) as given in [30].We found that M could not be always increased abovea rather low value without having at least one empty binin the M × M grid. This is due to the very inhomoge-neous distribution of phase space angles υ , υ , which isa known phenomenon [30, 32]. A low value for M leadsto low accuracy.There is another way to calculate the dynamic aperturebased on normal forms given to the first order in [30] as r NF = " π K K X k =1 { ρ ( α k , ¯ υ , ¯ υ ) + ρ ( α k , ¯ υ , ¯ υ ) } sin (2 α k ) , (6)where ρ ( α k , ¯ υ , ¯ υ ) and ρ ( α k , ¯ υ , ¯ υ ) are the nonlinearinvariants, which are the same as the Courant-Snyderinvariants. Since there is already an implementation ofellipse fit to the phase space variables in SIMPA, thevalues of ρ ( α k , ¯ υ , ¯ υ ) and ρ ( α k , ¯ υ , ¯ υ ) are easily ob-tained. It should be noted however, this sum can leadalso to non-negligible errors if the working point is on astrong resonance line and the form of the phase space plotis significantly different from an ellipse. Nevertheless, itis still the best first-order approximation of the dynamicaperture. For this reason we calculated with both Eq.(5) and Eq. (6) as a cross check.In our numerical experiments, the closed orbit was dif-ferent from the design orbit. The electron cooler has in-troduced an orbit distortion, which was corrected. How-ever, the correction was not perfect. The maximum orbitdeviation was less than 3 millimeters horizontally and lessthan 1 millimeter vertically. It could have been reducedfurther, but doing so for each point we investigated israther tedious. Also, in reality, the orbit of ELENA dur-ing the commissioning was bigger than that. B. Results
The dynamic aperture of six working points has beenevaluated to compare resonance conditions with a non-resonant case. The selected set of tunes is indicated onFIG. 7. These resonance lines were selected because theyare close to the region where the working point of ELENAwas set during the commissioning. The first workingpoint was chosen to be far from strong resonance linesat Q h = 2 . , Q v = 1 . N = 10 turns, K = 20 polar steps, with α k ∈ [0 , π/ J = 20 and the biggest possi-ble M , which is also indicated in TABLE I. The energyof the particles was set to 100 keV. We also calculatedthe errors of the dynamic aperture values in in TABLE Iwith Eq. (10) given in [33], it is ± . m ].Most of the dynamics aperture values in TABLE I areclose to 0.011. At first sight one is tempted to conclude FIG. 7. The numbered white points in the tune diagramindicate the tunes where the dynamic aperture was calculated.TABLE I. Dynamic apertures for the six points in FIG. 7calculated with averaging over the dynamics ( r d ) and withnormal forms ( r NF ). M is the angular grid size for r d calcula-tion. For the tunes falling onto a resonance line, the equationof the resonance is also given.Point r d [ m ] M r NF [ m ] Resonance condition1 0.0115 14 0.0115 NA2 0.012 8 0.0112 Q h − Q v = 13 0.01 6 0.01 3 Q h + 2 Q v = 104 0.0117 3 0.0114 5 Q h = 125 0.0113 9 0.011 − Q h + 3 Q v = 26 0.0113 8 0.0108 Q h + 4 Q v = 8 that the resonance lines of order 4 or 5 are already weakenough to not reduce further the physical aperture of themachine. The dynamic aperture however is defined asthe radius of a four dimensional hyper sphere and smalldifferences in the radius can hide significant phase spacevolume reductions at specific angles of α k .For this reason we found it informative to give also amore detailed description by plotting those initial emit-tances in FIG. 8 for each point in FIG. 7, which survived N = 10 turns. We repeated the same exercise with par-ticles having a momentum deviation dp/p = 2 × − ,which is the maximum momentum spread measured dur-ing the commissioning. The quadrupole currents wereadjusted for each point to have the same tunes as for thecase dp/p = 0. The results are plotted in FIG. 9 .In FIG. 8, 9 it is apparent that some of the resonancelines make the dynamics aperture significantly smaller e v e h e v e h e v e h e v e h e v e h e v e h FIG. 8. Initial emittances plotted in units of [ π mm mrad ]for particle with dp/p = 0. The horizontal axis correspondsto the horizontal initial emittance and the vertical axis to theinitial vertical emittance. Only those initial conditions areplotted which survived N = 10 turns. The numbers in theupper right corners correspond to the numbering in FIG. 7.The lines indicate the last connected initial emittances foreach angle α k . than the physical aperture. The physical acceptance ofELENA is determined by the size of vacuum chambersand the optics.There are a number of interesting features displayed inthese figures, some of them surpising. The coupling reso-nance Q h − Q v = 1 indicated by point two is rather wideand looks strong, but it does not significantly reduce thedynamic aperture. The smallest dynamic aperture wasobtained at point three, due to a 5th order resonance andnot a lower order resonance line. We did not include first,second and third order lines, those should obviusly beavoided during operation. Higher-order resonance lineshowever can not always be evaded. The incoherent spacecharge tune shift in ELENA at the lowest energy is sig-nificant. Its estimated value for a bunched beam at ar-0 e v e h e v e h e v e h e v e h e v e h e v e h FIG. 9. Initial emittances plotted in units of [ π mm mrad ]for particle with dp/p = 2 × − . rival to 100 KeV is ∆ Q = − .
03 [6]. The incoherentspace charge tune shift contributes to the tune footprintof the beam and some of the resonances shown in FIG.7 are unavoidable. Therefore it is important to choose aworking point which gives the most space free of harmfulresonance lines. Point one in FIG. 7 is a good choice, be-cause the surrounding resonance lines Q h − Q v = 1 (onpoint two) and Q h + 4 Q v = 8 (on point six) are ratherforgiving according to FIG. 8, 9.The initial emittances for the test particle to be trackedwere calculated by a linear model, that is the Courant-Snyder theory. It assumes that the values of the phasespace variables lie on an ellipse. This is only approxi-mately true under resonance conditions. Although theangle α k was changed linearly in equal steps, in FIG.8 case 5 for instance, shows distortions. This is becausethe phase space plot deviates significantly from an ellipseand the requested emittance is different from the actualinitial emittance, which is plotted.The dynamic aperture is defined as a function of turns. -0.008-0.006-0.004-0.002 0 0.002 0.004 0.006 0.008-0.015 -0.01 -0.005 0 0.005 0.01 x ’ [ r ad ] x [m] H -0.008-0.006-0.004-0.002 0 0.002 0.004 0.006 0.008-0.02 -0.015 -0.01 -0.005 0 0.005 0.01 0.015 x ’ [ r ad ] x [m] V FIG. 10. Phase space plots of a particle near a 4th orderresonance, lost after 14211 turns.
It is known to decrease with the inverse of the logarithmof the number of turns [33–35]. Indeed it is easy to findinitial conditions which leads to particle loss after morethan 10 turns, an example is given in FIG. 10. VI. CONCLUSIONS
We have applied the SIMPA code for long-term sym-plectic charged particle tracking in arbitrary static elec-tromagnetic fields on the ELENA machine. The algo-rithm described in a previous paper has been extended tohandle complex beam region geometries. A short descrip-tion of the enhanced algorithm was given. The improve-ments not only raised its generality, but also increasedthe speed by a factor 4.66.The frequency and dynamic aperture analysis identi-fied a number of 4th and 5th order resonance lines inthe tune diagram strong enough to reduce the dynamic1aperture below the physical one. What made the fre-quency and dynamic aperture analyses different, is thefact that we have not introduced any multipole error intoour model apart from the multipole components due tothe geometry of the magnets, which are inevitable. Allthe resonances seen in the frequency analysis are directconsequences of the geometry of the fields, even if themagnetic elements are manufactured perfectly.We showed in the frequency analysis section the effectof the magnetic system of the electron cooler on the beamdynamics by comparing the two cases, with and with-out electron cooler. The electron cooler introduced non-negligible magnetic perturbations, strengthening manyresonance lines. To our best knowledge there was nosimilar study done before with electron coolers.
A. Suggestions to improve the performance
The commissioning of the ELENA ring was finishedin November 2018, at the start of the 2.5 years longCERN wide shutdown. The deceleration efficiency wasestimated between 40 - 50 %, somewhat below the 60 %design value. Most of the losses occurred on the rampjust before the lowest energy plateau. This is not sur-prising, since both the incoherent tune shift due to spacecharge and the intrabeam scattering has the biggest ef-fects at low energy. Also the perturbation introduced bythe electron cooler is the strongest at the lowest energy.During the commissioning two working points were ex-plored extensively, Q h = 2 . , Q v = 1 .
45, and Q h =2 . , Q v = 1 .
44. These working points are not far fromthe 4th and 5th order resonances number 3, 4, 5 in FIG.7 and these resonance lines are not harmless according to FIG. 8, 9. The tune footprint of the beam certainlyoverlaps with some of the resonance lines. The incoher-ent tune shift of a bunched beam due to space chargeat arrival to 100 KeV energy is in the order of -0.03, de-pending on the beam intensity and beam size. This is nota negligible value, but we believe it is not big enough todrastically change the results of our analysis. One sup-porting evidence for this assumption, is that we haven’tobserved significant dependence of the deceleration effi-ciency on the beam intensity during the commissioning.After cooling at 100 keV energy, the incoherent tune shiftis in the order of -0.1, but only for the duration of the finalpart of the bunching before extraction, about 10 millisec-onds. No significant losses were observed related to thebunching process during the commissioning. Based onthe frequency and dynamic aperture analysis, we suggestto explore working point number one in FIG. 7, as it givesthe most available space without harmful resonances.The tune footprint of the beam with maximum emit-tances in all planes shown on FIG. 6 is mostly comingfrom the horizontal chromaticity which is ξ h = − . VII. ACKNOWLEDGMENTS
I would like to express my gratitude towards GianluigiArduini, Christian Carli and Massimo Giovannozzi fromCERN BE-ABP for their valuable comments and sugges-tions. [1] L. Bojt´ar, “Efficient evaluation of arbitrary static electro-magnetic fields with applications for symplectic particletracking,”
Nuclear Instruments and Methods A , vol. 948,p. 162841, 2019.[2] G. Lee-Whiting, “End effects in first-order theory ofquadrupole lenses,”
Nuclear Instruments and Methods ,vol. 76, no. 2, pp. 305 – 316, 1969.[3] M. Berz, B. Erd´elyi, and K. Makino, “Fringe field effectsin small rings of large acceptance,”
Phys. Rev. ST Accel.Beams , vol. 3, p. 124001, Dec 2000.[4] E. Forest and J. Milutinovic, “Leading order hard edgefringe fields effects exact in (1 + δ ) and consistent withmaxwell’s equations for rectilinear magnets,” Nuclear In-struments and Methods A , vol. 269, no. 3, pp. 474 – 482,1988.[5] B. D. Muratori, J. K. Jones, and A. Wolski, “Analyticalexpressions for fringe fields in multipole magnets,”
Phys.Rev. ST Accel. Beams , vol. 18, p. 064001, Jun 2015.[6] V. Chohan, C. Alanzeau, M. E. Angoletta, J. Bail-lie, D. Barna, W. Bartmann, P. Belochitskii, J. Bor-burgh, H. Breuker, F. Butin, M. Buzio, O. Capatina,C. Carli, E. Carlier, M. Cattin, T. Dobers, P. Chig- giato, L. Ducimetiere, T. Eriksson, and T. Zickler,
ExtraLow ENergy Antiproton (ELENA) ring and its TransferLines: Design Report . CERN, 04 2014.[7] J. Brackbill and D. Barnes, “The effect of nonzero ∇ . B on the numerical solution of the magnetohydrodynamicequations,” Journal of Computational Physics , vol. 35,no. 3, pp. 426 – 430, 1980.[8] C. E. Mitchell and A. J. Dragt, “Accurate transfer mapsfor realistic beam-line elements: Straight elements,”
Phys. Rev. ST , vol. 13, p. 064001, Jun 2010.[9] D. T. Abell, “Numerical computation of high-order trans-fer maps for rf cavities,”
Phys. Rev. ST Accel. Beams ,vol. 9, p. 052001, May 2006.[10] M. Venturini and A. J. Dragt, “Accurate computationof transfer maps from magnetic field data,”
Nuclear In-struments and Methods A , vol. 427, no. 1, pp. 387 – 392,1999.[11] P. Chang, “A general method for symplectic particletracking in a three-dimensional magnetic field,” tech.rep., SLAC-PUB-8514, 2000.[12] J. Bahrdt, M. Scheer, and G. Wustefeld, “Tracking sim-ulations and dynamic multipole shimming for helical un- dulators,” in Proceedings of EPAC , 2006.[13] M. Titze, J. Bahrdt, and G. W¨ustefeld, “Symplectictracking through straight three dimensional fields bya method of generating functions,”
Phys. Rev. Accel.Beams , vol. 19, p. 014001, Jan 2016.[14] F. Mackay, R. Marchand, and K. Kabin, “Divergence-free magnetic field interpolation and charged particle tra-jectory integration,”
Journal of Geophysical Research:Space Physics , vol. 111, no. A6, 2006.[15] A. J. Dragt,
Lie methods for nonlinear dynamics with ap-plications to accelerator physics . University of MarylandPress, 2018.[16] P. Belochitskii and V. L. Joergensen. Private Communi-cation, 2018.[17] M. Tao, “Explicit symplectic approximation of nonsep-arable hamiltonians: Algorithm and long time perfor-mance,”
Phys. Rev. E , vol. 94, p. 043303, Oct 2016.[18] H. Grote and F. Schmidt, “Mad-x users guide,” 2019. http://madx.web.cern.ch/madx/ .[19] K. Atkinson and W. Han,
Spherical Harmonics andApproximations on the Unit Sphere: An Introduction .Springer, Berlin, Heidelberg, 2012.[20] E. Forest, F. Schmidt, and E. McIntosh, “Introductionto the polymorphic tracking code,”
KEK report , 2002-3.[21] F. Asvesta, N. Karastathis, and P. Zisopoulos, “Py-NAFF: A Python module that implements NAFF.” https://github.com/nkarast/PyNAFF .[22] K. L. Brown,
A first- and second-order matrix theory forthe design of beam transport systems and charged particlespectrometers . Stanford, CA: SLAC, 1972.[23] W. Scandale, “Dynamic aperture,” in
CERN AcceleratorSchool: Course on Advanced Accelerator Physics (CAS) ,pp. 0109–146, 7 1994.[24] J. Laskar, “Secular evolution of the solar system over10 million years,”
Astronomy and Astrophysics , vol. 198, pp. 341–362, June 1988.[25] R. Bartolini, A. Bazzani, M. Giovannozzi, W. Scandale,and E. Todesco, “Tune evaluation in simulations and ex-periment,”
Particle Accelerators , vol. 52, pp. 147 – 177,1996.[26] R. Bartolini, M. Giovannozzi, W. Scandale, A. Bazzani,and E. Todesco, “Precise measurement of the betatrontune,”
Particle Accelerators , vol. 55, pp. 1 – 10, 1996.[27] J. Laskar, “Frequency analysis for multi-dimensional sys-tems: Global dynamics and diffusion,”
Phys. D , vol. 67,no. 13, p. 257281, 1993.[28] J. Laskar, “The chaotic motion of the solar system: A nu-merical estimate of the size of the chaotic zones,”
Icarus ,vol. 88, no. 2, pp. 266 – 291, 1990.[29] E. Todesco, M. Giovannozzi, and W. Scandale, “Fastindicators of long term stability,”
Particle Accelerators ,vol. 55, pp. 27–36, 1996.[30] E. Todesco and M. Giovannozzi, “Dynamic aperture esti-mates and phase-space distortions in nonlinear betatronmotion,”
Phys. Rev. E , vol. 53, pp. 4067–4076, Apr 1996.[31] M. Giovannozzi and E. Todesco, “Numerical methods toestimate the dynamic aperture,”
Particle Accelerators ,vol. 54, pp. 203–212, 1996.[32] J. D. Meiss, “Symplectic maps, variational principles,and transport,”
Rev. Mod. Phys. , vol. 64, pp. 795–848,Jul 1992.[33] M. Giovannozzi, W. Scandale, and E. Todesco, “Dynamicaperture extrapolation in presence of tune modulation,”
Phys. Rev. E , vol. 57, pp. 3432–3443, 1998.[34] M. Giovannozzi, W. Scandale, and E. Todesco, “Inverselogarithmic extrapolation of survival plots in hadron col-liders,” in
Beam Dynamics Newsletters , vol. 12, 1996.[35] M. Giovannozzi, W. Scandale, and E. Todesco, “Predic-tion of long term stability in large hadron colliders,”