MESA Technical Note: Beam Breakup Instability Threshold Current
MMESA Technical Note: Beam Breakup Instability Threshold Current
S. Glukhov ∗ and O. Boine-Frankenheim TEMF TU Darmstadt, Darmstadt
C. Stoll and F. Hug
Institut f¨ur Kernphysik, JGU Mainz (Dated: February 2, 2021)MESA (Mainz Energy recovery Superconducting Accelerator) is an energy recovery linac (ERL)which is under construction at Johannes Gutenberg University in Mainz. It will be operated inexternal beam (EB) mode with 150 µ A electron beam at 155 MeV and energy recovery (ER) modewith 1 mA (first stage) and later 10 mA (second stage) electron beam at 105 MeV. An importantfactor which may limit performance of the machine is a beam breakup (BBU) instability whichmay occur due to excitation of higher-order modes (HOMs) in superconducting RF cavities. Thiseffect occurs only when the injected beam current exceeds a threshold value. The aim of the presentwork is to develop a software for reliable determination of the threshold current in MESA, findmain factors which may change its value and finally make a decision concerning capability of MESAoperation at 10 mA and need for additional measures for suppressing BBU instability.
INTRODUCTION
The schematic MESA layout in ER mode is presentedin Fig. 1. The machine consists of the injection arcT0, two accelerating modules connected with recircula-tion arcs T1, T2, T3 and a pseudo internal target (PIT)arc where experiments are performed. Each acceleratingmodule consists of two 9-cell TESLA cavities connectedas shown in Fig. 2. 5 MeV electron bunches are injectedthrough T0 and move through T1, T2 and T3 gaining12.5 MeV on each cavity pass. The length of PIT arc isadjusted in a way that bunches arrive to the next cavitywith a 180 ◦ phase shift (in decelerating phase) and thenmove through T3, T2 and T1 losing 12.5 MeV on eachcavity pass [1]. FIG. 1. Schematic view of MESA.FIG. 2. CAD model of MESA accelerating module.
The BBU instability mechanism is as follows. Electronbunches periodically injected into ERL excite parasiticHOMs in superconducting RF cavities. These HOMsmay act back on the bunches on subsequent passes.Dipole modes have the most impact because among other non-axisymmetric modes their field is concentrated closerto the cavity axis. A bunch kicked by a HOM receives anadditional transversal momentum. When the bunch re-turns into the same RF cavity after recirculation this mo-mentum transforms into a transversal offset. A transver-sally shifted bunch changes HOM voltage in the cavity,the change depends on the HOM voltage value of theprevious recirculation. This may lead to unstable mo-tion [2, 3].
MAIN DEFINITIONS
For the purpose of BBU simulations an accelerator canbe represented as a sequence of arcs and thin RF cavitiesat the end of each arc in the order seen by one bunch onits way from injector to beam dump. Each cavity mayappear several times in this sequence depending on num-ber of recirculations. Regions between two cavities of thesame accelerating module are straight but for simplicitythey will be also referred to as arcs in the following text.Each cavity is assigned with a certain number of dipoleHOMs. Adjacent HOMs form a pair with close frequen-cies and nearly orthogonal polarization. In all subsequentsimulations the same (and only one) HOM pair is excitedin all cavities. The beam is represented as a sequence ofpointlike on-energy equidistantly injected bunches. It isnecessary to clearly define the characteristics for each ofthese entities:1. Accelerator: • accelerated particle rest energy E ; • RF voltage period t RF ; • number of cavities N c ; • number of HOMs N m in each cavity( N m = 2 in all simulations); a r X i v : . [ phy s i c s . acc - ph ] F e b • number of arcs N a .2. Arc: • × T ( i ) ; • length L i ; • reference particle energy E i ; • reference particle velocity v i = c (cid:113) − ( E /E i ) ; • reference particle time of flight ∆ t i = L i /v i .3. Dipole HOM: • frequency f λ or ω λ = 2 πf λ ; • Q-factor Q λ ; • damping factor Γ λ = ω λ Q λ ; • shunt impedance ( R/Q ) λ ; • reduced shunt impedance( R/Q ) (cid:48) λ = ( R/Q ) λ ( ω λ /c ) ; • polarization angle θ λ .4. Beam: • injector current I ; • number of injected bunches N b ; • time interval between bunches t b ( t b = t RF for the MESA case); • bunch charge q = I t b .5. Bunch: • X = ( x, p x , y, p y , z, p z ) T ; • longitudinal coordinate Z along reference or-bit in laboratory frame.6. Combined entities: • average particle momentum in the cavity atthe end of the i -th arc p i = ( E i v i /c + E i +1 v i +1 /c ) / • HOM λ voltage in the j -th cavity V ( j ) λ .A small remark concerning HOM properties is nec-essary. An on-axis particle can be kicked by a HOMtransversally; a line perpendicular to the direction of thekick is called polarization axis. The polarization angle isan angle between the horizontal axis (in accelerator co-ordinates) and the normal to the polarization axis. Theprime in ( R/Q ) (cid:48) λ (measured in Ω / m ) is introduced todistinguish it from ( R/Q ) λ (measured in Ω) used in [4]and [5]. The usage of ( R/Q ) (cid:48) λ leads to simpler formulaebecause the impedance of the dipole HOMs is quadrat-ically dependent on the distance from the polarizationaxis and this is the proportionality coefficient. BUNCH TRACKING
When a bunch passes the i -th arc followed by the j -thcavity, its coordinates are transformed as follows: X (cid:55)→ T ( i ) X + 1 p i (cid:88) λ (cid:61) V ( j ) λ (0 , cos θ λ , , sin θ λ , , T . The HOM voltages in the cavity also change: V ( j ) λ (cid:55)→ V ( j ) λ e ( iω λ − Γ λ )∆ t i + q (cid:18) RQ (cid:19) (cid:48) λ ( x cos θ λ + y sin θ λ ) ,λ = 1 . . . N m . (1)The most obvious tracking algorithm is to unroll thelattice onto the positive semiaxis, place initially all thebunches on negative semiaxis equidistantly and thenmove them through the lattice with their velocities stop-ping each time one of them encounters a cavity. Usuallyin simulations N b > , therefore it is important to usea tracking algorithm which has complexity not more thanlinear by this parameter. On each step a bunch shouldbe found which is the first to encounter a cavity, but notall the bunches should be checked, only the first bunchin each arc. When a bunch leaves an arc, the next onebecomes the first. Therefore number of bunches to checkis always not more than N a .The longitudinal position of all bunches also shouldnot be updated on each step. It is enough to have a lab-oratory time counter and an individual counter for eachbunch which shows when its position was updated lasttime. A similar counter is also necessary for each cavityto show when HOM voltages in it were updated last time.The difference between this counter and laboratory timeshould be used in (1) instead of ∆ t i . BEAM BREAKUP INSTABILITY THEORY(SIMPLEST CASE)
The main equations of the BBU theory were derived in[4] and [5]. The approach developed in [4] is much sim-pler but suitable only for the case of recirculating linacwhere initially injected and all the recirculated bunchesare at the same RF phase. Laplace transformation of V λ ( t ) is introduced in [5] to overcome this limitation. Itwill be shown in this section that such a complication isnot necessary and the approach used in [4] can be eas-ily extended to the case of an arbitrary RF phase of therecirculated bunches.Consider a charge q e moving at the speed of light c through an accelerator parallel to its axis with a trans-verse offset d and a test particle moving on-axis with thesame speed and time delay t . Transversal wake function W t ( t ) is a transversal Lorentz force acting on the testparticle integrated over path of interest per unit chargesand transversal offset: W t ( t ) = 1 q e d (cid:90) path ( E x ( z, z/c + t ) − cB y ( z, z/c + t )) dz . For point charges in an RF cavity this becomes: W t ( t ) = (cid:88) λ (cid:18) RQ (cid:19) λ ω λ c e − Γ λ t sin( ω λ t ) , t (cid:62) , t < , (2)where index λ iterates over all HOMs which can be ex-cited in the cavity.Consider an ERL with one recirculation turn andone thin RF cavity containing one dipole HOM. TheHOM voltage V λ ( t ) had given additional deflection tothe bunch; this deflection transformed into additionaltransversal offset d ( t ) after one recirculation: d ( t ) = T q e pc V λ ( t − t r ) , (3)where t r is recirculation time. On the other hand, a beamcurrent I ( t ) passing the cavity at a transversal position d ( t ) creates the following deflecting voltage: V λ ( t ) = (cid:90) t −∞ W t ( t − t (cid:48) ) I ( t (cid:48) ) d ( t (cid:48) ) dt (cid:48) . (4)The beam is a series of short bunches: I ( t ) = ∞ (cid:88) m = −∞ I t b δ D ( t − mt b ) , (5)where δ D is the Dirac delta-function.According to [5], integer ( n r ) and fractional ( δ ) partof the ratio t r /t b can be introduced as follows: t r = ( n r − δ ) t b , (cid:54) δ < . (6)Beam current through the cavity is schematically shownin Fig. 3. FIG. 3. Beam current through the cavity (schematic view).
For the MESA case all arcs should have δ = 0, exceptof the PIT arc which should have δ = . This providessteady acceleration of the bunches from the injector tothe experimental area and steady deceleration of the used bunches on the way to the beam dump. However formu-lae for stability analysis can be obtained for arbitraryvalue of δ .Stability analysis implies searching for a harmonic so-lution for V λ ( t ) but it can not be a pure harmonic func-tion V λ e − iωt because at least it is not smooth at t = nt b where n is integer. As a workaround V λ ( t ) can be con-sidered as a discrete function on a uniform grid with aperiod t b . Samples should be taken at the time points t m = ( m + δ ) t b for integer m (when bunches enter thecavity for the first time) because V λ ( t m ) appear in theright-hand side of (4) after substitution of (3) and (5).Then V λ ( t ) contains harmonics with frequencies ω + πnt b and can be represented as Fourier series: V λ ( t ) = ∞ (cid:88) n = −∞ V n e − i (cid:16) ω + πntb (cid:17) t . (7)Fourier amplitudes V n are not important for the stabilityanalysis because all harmonics are either stable or unsta-ble depending on the sign of (cid:61) ω .Substituting (2), (3) and (5) into (4) and taking intoaccount that summation index should take integer valuesone obtains: e − iωt = K i t/t b − δ (cid:88) m = −∞ e − Γ λ t e Γ λ mt b ×× (cid:0) e iω λ t e − iω λ mt b − e − iω λ t e iω λ mt b (cid:1) e − iωmt b e iωt r , where K = T q e pc (cid:18) RQ (cid:19) λ ω λ c I t b = T q e p (cid:18) RQ (cid:19) (cid:48) λ I t b . (8)Here the only difference to [4] is presence of δ in theupper summation limit. After summation of geometricseries this becomes:1 = K i e iωn r t b − δ Γ λ t b ×× (cid:18) e iδω λ t b − e ( − Γ λ + iω + iω λ ) t b − e − iδω λ t b − e ( − Γ λ + iω − iω λ ) t b (cid:19) (9)or K = 2 e − iωn r t b + δ Γ λ t b ×× cos ( ωt b + i Γ λ ) − cos( ω λ t b ) e Γ λ − iωt b sin( δω λ t b ) − sin(( δ − ω λ t b ) (10)This is a dispersion relation connecting injector currentand effective HOM voltage frequency. The result is ex-actly the same as in [5]. BEAM BREAKUP INSTABILITY THEORY
Several generalizations are necessary:1. multiple recirculations;2. multiple cavities;3. multiple HOMs in each cavity;4. two transversal degrees of freedom.From now on the HOM set of the whole machine shouldbe considered. If a HOM with the same properties is ex-cited in two different cavities then it should be treated astwo different HOMs. Therefore, HOMs should be num-bered continuously throughout the machine. HOM num-bering is not affected by the number of recirculations andthe machine topology in general. Vice versa, transportmatrices depend strongly on the topology but not on theHOM properties. Transport matrix from the I -th pass ofthe l -th cavity to the J -th pass of the m -th cavity maybe denoted as T IJlm . Values t r , n r and δ should be re-placed with t IJlm , n IJlm and δ IJlm correspondingly. Also forconvenience particle momentum at the I -th pass of the l -th cavity may be denoted as p Il .HOM polarization should be taken into account tostudy bunch dynamics with two transversal degrees offreedom. Assume that dipole HOMs λ and µ with po-larization angles θ λ and θ µ are excited in cavities at thebeginning and at the end of a linear lattice region withtransport matrix T IJ ( lm ) . At the entrance a bunch re-ceives a transversal kick proportional to the distance tothe first polarization axis in a direction perpendicular tothis axis. At the exit this kick transforms into additionaldistance to the second polarization axis with the follow-ing proportionality coefficient: (cid:98) T IJλµ == T IJl ( λ ) m ( µ )12 cos θ λ cos θ µ + T IJl ( λ ) m ( µ )34 sin θ λ sin θ µ ++ T IJl ( λ ) m ( µ )32 sin θ λ cos θ µ + T IJl ( λ ) m ( µ )14 cos θ λ sin θ µ . Here the fact is explicitly indicated that the index of aHOM unambigiously defines the index of the cavity con-taining it. This expression should be used in (3) insteadof T (which corresponds to θ λ = θ µ = 0).If summation over HOMs from (2) is introduced into(4), then from (9) follows that inverse injector current I is an eigenvalue of N c N m × N c N m matrix W ( ω ) withelements W µλ ( ω ): W µλ ( ω ) = t b q e i (cid:88) I ; J>I (cid:98) T IJλµ p Il ( λ ) (cid:18) RQ (cid:19) (cid:48) µ e iωn IJl ( µ ) m ( λ ) t b ×× (cid:32) e ( iω λ +Γ λ ) δ IJl ( µ ) m ( λ ) t b − e ( − Γ λ + iω + iω λ ) t b − e ( − iω λ +Γ λ ) δ IJl ( µ ) m ( λ ) t b − e ( − Γ λ + iω − iω λ ) t b (cid:33) . COMPLEX CURRENT PLOT
The eigenvalues of matrix W ( ω ) can not be found ana-lytically, thus indirect way is used. There are no infinitely growing solutions for I = 0, this means that all ω ’s havenegative imaginary parts. An instability emerges when I exceeds a threshold value, this means that one of the ω ’sgets a positive imaginary part. Therefore, ω should bevaried along the real axis within some region which spansthe frequencies ω λ of all excited HOMs. In general, theresulting current values will be complex, then the valuewith zero imaginary part and the smallest positive realpart is the threshold current.This leads to complex current plot technique when I ( ω ) values are plotted as a curve on the complex planeand then the closest to zero intersection of this curvewith the positive real semiaxis is searched for. But incalculations this curve is represented as a set of pointswith varying step between them. Moreover, the numberof curves is equal to the number of HOMs, and eigen-value solver does not guarantee the order of the eigen-values. This means that their order may be different foradjacent ω values. Without special countermeasures thiswould lead to jumps between curves and possible spuri-ous intersections. The following algorithm was used tosolve the problem: for each I on the current step theclosest value from the previous step is found (values withsmaller | I | are considered first), and values from the cur-rent step are sorted accordingly. Then the order of thecurves is preserved and for sufficiently small step in ω intersection I ∗ with real axis between adjacent values I and I can be found as follows: I ∗ = (cid:60) I (cid:61) I − (cid:60) I (cid:61) I (cid:61) I − (cid:61) I , if (cid:61) I (cid:61) I (cid:54) . PARAMETER SPREAD
This section is devoted to the estimation of differ-ent uncertainties which may affect the threshold current.Cavity manufacturing tolerances may lead to a HOM pa-rameters spread. Beamline element misalignments andpower supply errors may lead to a transport matrix ele-ments spread.The frequencies and Q-factors of the first 36 dipoleHOMs were measured twice for 4 MESA cavities (No. 7,8, 9 and 10): during cold tests at DESY Hamburg andhorizontal tests at HIM [2]. The basic values of frequen-cies and Q-factors were taken from [6]. The maximaldeviations of measured frequencies from the basic valuesare presented in Fig. 4. The maximal ratios of measuredQ-factors and the basic values are presented in Fig. 5.In a cylindrically symmetric cavity the polarizationangle of the first dipole HOM in a pair is arbitrary inthe range − π/ . . . π/
2, the second is exactly orthogo-nal. In presence of high frequency power couplers andHOM couplers these angles have certain values, more-over the exact orthogonality is lost. However, produc-tion tolerances may randomize polarization angles backbut deviation from orthogonality may remain. Shuntimpedances and polarization angles were calculated withCST MICROWAVE STUDIO. Deviations from orthogo-nality within HOM pairs are presented in Fig. 6.
FIG. 4. Measured frequency spread.FIG. 5. Maximal Q-factor ratio in measurements and simu-lations.FIG. 6. Deviation from orthogonality in dipole HOM pairs.
Quadrupole gradient and quadrupole roll errors maylead to matrix elements variation. For each configura-tion separate lattice file with random error values is gen- erated, then the matrix elements are calculated using elegant [7]. General tolerance to quadrupole gradienterrors should be taken as a value of ∆
K/K , usually itdoes not exceed 10 − . General tolerance to quadrupoleroll errors usually is not more than 1 mrad.Other lattice related errors do not affect the transportmatrix elements significantly because the optics is linear.Reference orbit offsets and cavity misalignments mayonly change the equilibrium HOM voltage (if any) butnot its frequency and, therefore, can not affect thresholdcurrent. It was shown in [5] that the beam deflectioncaused by HOMs may lead to large reference orbit dis-tortions but this happens at much higher current valuesthan the BBU threshold current. WORST CASE CONFIGURATION
Assume that each HOM interacts only with itself ondifferent passes but not with other HOMs (of the sametype but excited in different cavities and with slightlydifferent parameters due to cavity manufacturing toler-ances). Therefore, for each HOM the worst combinationof cavity manufacturing tolerances which leads to mini-mal possible threshold current can be found.Introduce in (10) the following notation: (cid:15) = Γ λ t b , ω = ω λ + ∆ ω . If ∆ ωt b (cid:28) (cid:15) (cid:28) K ≈ − ωt b cos( ωt r ) + (cid:15) sin( ωt r ) ++ i ( (cid:15) cos( ωt r ) − ∆ ωt b sin( ωt r ))) .K is real, therefore∆ ωt b ≈ (cid:15) ctg( ωt r )and K ≈ − (cid:15) sin( ωt r ) . For one cavity with one HOM, multiple recirculationsand two transversal degrees of freedom this leads to: I λth ≈ − c q e (cid:16) RQ (cid:17) λ Q λ ω λ (cid:88) I ; J>I (cid:98) T IJλλ p Il ( λ ) sin (cid:16) ( ω λ +∆ ω ) t IJl ( λ ) l ( λ ) (cid:17) . (11)The Q-factor and the shunt inpedance appear only in acommon factor, thus their effect on the threshold currentis straightforward. So, the maximal possible Q-factorvalue should be chosen and the denominator of the secondfraction in (11) should be maximized by absolute valueover variables ∆ ω and θ λ .The matrix elements T and T are always small be-cause MESA has no transversal coupling, thus all ex-trema are at θ λ = 0 and θ λ = π/
2. Then only ∆ ω isto be varied. The relative HOM polarization is not im-portant because interaction between HOMs is neglected.Therefore, two configurations for each HOM pair shouldbe checked: the first HOM is horizontally polarized, thesecond is vertically polarized and vice versa. The onewith smaller threshold current is the worst case configu-ration. SIMULATION RESULTS
In Fig. 7 threshold currents for 100 different randommachine configurations are shown, with two worst caseconfigurations added. The frequency, Q-factor and de-viation from orthogonality spread for each HOM pair ischosen in accordance to the maximal spread obtainedin measurements, HOM polarization angle is random,quadrupole gradient relative error is 10 − , quadrupoleroll error is 1 mrad. The technique for predicting theworst case configuration proposed in the present paperseems to be quite reliable. FIG. 7. Threshold current for 100 random machine configu-rations and 2 worst case configurations.
The minimal threshold current I th = 68 .
28 mA was ob-tained for the worst case configuration with HOMs No.39 and 40 at 1.87 GHz. Figure 8 shows the central part ofcomplex current plot for this particular case, the thresh-old current value is marked with a red point. Currentvalues of I = 67 mA and I = 70 mA were chosen todemonstrate stable and unstable motion. For both valuesthe HOM voltages in all cavities over time are shown inFig. 9 and 10. It can be seen that after a relatively shorttransition process HOM voltages begin to constantly fall(below threshold) or rise (above threshold).This is an example of a visual check performed for sev-eral tens of machine configurations to make sure that theresults obtained with both techniques are consistent. Inprinciple, the threshold current may be found using only bunch tracking with varying current (e.g. with bisectionmethod, as in the bi code by I. Bazarov [8]) but careshould be taken because the exponential fit of the volt-ages or final coordinates may give different decrementsigns depending on the considered time interval. FIG. 8. Complex current plot for minimal threshold currentconfiguration. The threshold current value is marked with ared point.FIG. 9. HOM voltages time dependence for a current belowthreshold.FIG. 10. HOM voltages time dependence for a current abovethreshold.
CONCLUSION
A set of scripts in Python were developed for calcula-tion of the BBU instability threshold current using twodifferent techniques: a complex current plot technique(based on the well-known BBU instability theory) anddirect bunch tracking. The results obtained with bothtechniques are consistent for various dipole HOM andlattice configurations in MESA. A method proposed forpredicting the worst case configuration with the small-est threshold current works also well. With the mostdangerous uncertainties taken into account, the thresh-old current is not less than 68 mA which is well abovethe 10 mA design current for the second stage of MESAoperation in ER mode. Therefore no additional measuresneed to be taken to avoid negative consequences of theBBU effect.The authors thank W. M¨uller and W. Ackermann forconsultations in TESLA cavity and CST MICROWAVESTUDIO related topics, A. Khan and D. Simon forMESA lattice files, H. De Gersem for assistance in prepa-ration of the present paper. ∗ [email protected] [1] F. Hug, K. Aulenbacher, R. Heine, B. Ledroit, and D. Si-mon, MESA - an ERL Project for Particle Physics Exper-iments, in (2017) p. MOP106012.[2] C. Stoll, F. Hug, and D. Simon, Beam Break Up Simula-tions for the MESA Accelerator, in (2018) p. MOPSPP009.[3] C. Stoll and F. Hug, Beam Breakup Simulations for theMainz Energy Recovering Superconducting AcceleratorMESA, J. Phys. Conf. Ser. , 012111 (2019).[4] G. Krafft, J. Bisognano, and S. Laubach, Calculat-ing beam breakup in superconducting linear accelerators(1990).[5] G. H. Hoffstaetter and I. V. Bazarov, Beam-breakup in-stability theory for energy recovery linacs, Phys. Rev. STAccel. Beams , 054401 (2004), arXiv:physics/0405106.[6] W. Ackermann, H. De Gersem, C. Liu, and T. Weiland, Eigenmode Calculations for the TESLA Cavity Consid-ering Wave-Propagation Losses through Fundamental andHigher-Order Mode Couplers , Tech. Rep. (2015).[7] M. Borland, elegant: A Flexible SDDS-Compliant Codefor Accelerator Simulation, in6th International Com-putational Accelerator Physics Conference (ICAP 2000)