Beam-Breakup Instability Theory for Energy Recovery Linacs

Here we will derive the general theory of the beam-breakup instability in recirculating linear accelerators, in which the bunches do not have to be at the same RF phase during each recirculation turn. This is important for the description of energy recovery linacs (ERLs) where bunches are recirculated at a decelerating phase of the RF wave and for other recirculator arrangements where different RF phases are of an advantage. Furthermore it can be used for the analysis of phase errors of recirculated bunches. It is shown how the threshold current for a given linac can be computed and a remarkable agreement with tracking data is demonstrated. The general formulas are then analyzed for several analytically solvable cases, which show: (a) Why different higher order modes (HOM) in one cavity do not couple so that the most dangerous modes can be considered individually. (b) How different HOM frequencies have to be in order to consider them separately. (c) That no optics can cause the HOMs of two cavities to cancel. (d) How an optics can avoid the addition of the instabilities of two cavities. (e) How a HOM in a multiple-turn recirculator interferes with itself. Furthermore, a simple method to compute the orbit deviations produced by cavity misalignments has also been introduced. It is shown that the BBU instability always occurs before the orbit excursion becomes very large.


I. INTRODUCTION
Synchrotron light sources based on Energy Recovery Linacs (ERLs) show promise to deliver X-ray beams with both brilliance and X-ray pulse duration far superior to the values that can be achieved with storage ring technology. This is due to the fact that the emittances in an ERL are largely determined by a laser-driven source, technology which has been improving steadily over the years, and will undoubtedly improve further. To generate high brilliance high flux X-rays it is necessary to accelerate beams to the energies (several GeV) and with the currents (several 100 mA) that are typical in these storage rings. This would require that the linac delivers a power of order 1 GW to the beam. Without somehow recovering this energy after the beam has been used, such a device would be practically unfeasible. Energy recovery [1] can be achieved by decelerating high energy electrons to generate cavity fields which in turn accelerate new electrons to high energy. With this, large beam powers that are not accessible in a conventional linac, can be produced.
Several laboratories have proposed high power ERLs for different purposes. Designs for light production with different parameter sets and various applications are being worked on by Cornell University [2,3], BNL [4], Daresbury [5], TJNAF [6], JAERI [7], the University of Erlangen [8], Novosibirsk [9], and KEK [10]. TJNAF has incorporated an ERL in its design of an electronion collider (EIC) [11] for medium energy physics, while BNL is working on an ERL-based electron cooler [12] for the ions in the relativistic ion collider (RHIC). The work at TJNAF, JAERI and Novosibirsk is based on existing ERLs of relatively small scale.
One important limitation to the current that can be accelerated in such an ERL is given by the beam-breakup (BBU) instability. The size and cost of all these new accelerators certainly requires a very detailed understanding of this limitation.
A theory of BBU instability in recirculating linacs, where the energy is not recovered in the linac, but where energy is added to the beam when it returns after each recirculation turn, was presented in [13]. Such a linac can consists of many cavities and several recirculation turns can be used. This original theory was additionally restricted to scenarios where the bunches of the different turns are in the linac at about the same accelerating RF phase, such as in the so-called continuous wave (CW) operation where every bucket is filled. Tracking simulations [14] compared well with this theory. In the following we therefore refer to it as the CW recirculator BBU theory.
It determines above what threshold current I th the transverse bunch position x displays undamped oscillations in the presence of a higher order mode (HOM) with frequency ω λ . If there is only one higher order mode and one recirculation turn with a recirculation time t r in the linac, the following formula is obtained for T 12 sin ω λ t r < 0: where c is the speed of light, T 12 is the element of the transport matrix that relates initial transverse momentum p x before and x after the recirculation loop, e is the elementary charge, (R/Q) λ Q λ is the impedance (in units of Ω) of the higher order mode driving the instability, Q λ is its quality factor. A corresponding formula had already been presented in [21]. Occasionally, additional factors are found when this equation is stated [15,16,19]. We give a concise derivation which shows that no such additional factors are required.
In the past efforts have been made to derive threshold currents from the analysis of experimental data obtained at the TJNAF FEL-ERL [16,19] using the CW recirculator BBU theory. The data has not been interpreted satisfactorily, and at least part of the reason could be that the underlying theory had not been derived for ERL operation. The theory presented in this paper should therefore be used to extend and improve these previous analysis results, as it is directly applicable to ERL operation.
Here we will derive the general theory of the beambreakup instability in recirculating linear accelerators, in which the bunches do not have to be at the same RF phase during each recirculation turn, similar to what has first been presented in [20]. First we treat the simplest case of one dipole HOM and one recirculation turn. Then we allow many HOMs and many recirculations, and finally we analyze several analytically solvable cases.

II. ONE DIPOLE HOM AND ONE RECIRCULATION
For recirculating linacs, the simplest case of one HOM and one recirculation loop has long been described [17]. Previous theories which assume that the recirculation time is an integral number of RF periods should not be used when investigating ERLs. Next we derive more general formulas that may be applied for arbitrary recirculation times.

A. The Dispersion Relation
In the simplest model of multi-pass beam breakup, bunches are injected into a cavity, which is assumed to have one dipole HOM (e.g. T M 11 -like mode), accelerated in the cavity and then recirculated to pass the cavity a second time before they are ejected. In the case of a twoturn recirculating linac, each bunch would be accelerated on both passes through the linac and ejected to a user area. The RF phase of the bunch would therefore be approximately the same on both passes. In an ERL, the RF phase on the second pass through the cavity is shifted by π with respect to the first pass so that the energy that the bunch gains in the first pass is returned to the cavity during the second pass and the bunch is ejected with reduced energy into a beam dump.
If a dipole HOM is excited in the cavity, then a bunch that enters the cavity on axis experiences a transverse kick and starts to oscillate around the design orbit of the recirculation loop and returns to the cavity with a transverse offset. This offset leads to a change in the energy of the HOM. If it increases the HOM energy, transverse kicks experienced by subsequent bunches will be larger, which will in turn lead to a further growth of the HOM energy once the kicked bunch returns to the linac: an instability develops.
To describe this effect, we use the ideas and nomenclature from [18]. When a current I(t ′ ) passes the cavity during its recovery loop at a time t ′ , the charge I(t ′ )dt ′ with the transverse offset x(t ′ ) excites the dipole HOM, creating a transverse momentum for particles traveling through the cavity subsequently at time t, where the wake function W (τ ) describes the transverse force at time τ after the HOM was excited. The momentum transfer is described by an effective transverse voltage of the HOM, V (t) = c e ∆p x (t). Assuming that all bunches are injected on the cavity's central axis, they do not excite dipole HOMs on their first path through the cavity. However, the effective transverse voltage of the HOM determines what kick ∆p x (t) the bunch sees and what position it will have when it returns to the cavity after the recirculation time t r . The transfer matrix element T 12 maps the transverse momentum p x (t) to x(t + t r ) = T 12 p x (t). Inserting this into Eq. (2) leads to an integral equation for the HOM's effective voltage, To solve this integral equation, one assumes that the current is a continuous stream of short pulses being injected at multiples of an interval between bunches t b , so that the current on the second turn is given by δ D being the Dirac-delta function. Note that t b is an integer multiple of the RF circulation time t 0 = 2π/ω 0 for the RF frequency ω 0 . We write the recirculation time in terms of the time t b between bunches as with an integer n r and δ ∈ [0, 1). For a recirculating linac one has and for an ERL one has δt b ≈ (n + 1 2 )t 0 for some integer n. A "+" sign in Eq. 5 that defines δ may seem more natural but our choice leads to simplified equations.
The HOM voltage at a time t ∈ [nt b + t r , nt b + t r + t b ) is given by Evaluating this at the time t = nt b + t r when the recirculated bunches pass through the cavity leads to In the CW recirculator BBU theory this difference equation was dealt with by assuming that the voltage can be written as V (t) = V 0 e −iωt for t = nt b where a positive imaginary part of ω indicates instability. Note that this does not require V (t) to be a harmonic function, but that it can be a linear combination of harmonic functions with frequencies ω + m 2π t b for integers m. This distinction has not been always made clear and is a potential source of confusion [19]. One obtains the equation The smallest value of the current I 0 for which there is a real ω is the threshold current I th of the instability. We proceed by writing V (t) in terms of its Laplace transform, retaining all possible frequencies in HOM voltage, which automatically enables proper description of arbitrary recirculating configuration: Note that this is not the conventional way of writing a Laplace transform. We have chosen this notation in order to make it appear more similar to a Fourier Transform. It also makes the subsequent notation more similar to the CW recirculator BBU theory. The Laplace transform is used rather than the Fourier transform since we want to analyze the onset of instability where the frequencies ω become complex. With the following definitioñ we obtainṼ SinceṼ Σ (ω) is periodic with 2π/t b , it has a Fourier series, and its Fourier coefficients are V (nt b ). This shows that V Σ (ω) does not vanish. We can therefore choose t = (n + n r )t b in Eq. (7) and sum over n, This finally yields the dispersion relation between I 0 and ω which can be used for all δ. A corresponding derivation, which has treated the beam recirculation in a way that can be applied to ERLs, has been presented in [20]. We believe that this paper should be referenced more often, since it is hardly referenced, whereas the earlier papers with CW recirculator BBU theory, which was not derived for ERLs, are often referenced in the context of ERLs, where they are not strictly applicable.

B. The Far-Field Wake
The sum in the dispersion relation can be computed when the far-field approximation for the wake function is used, With The dispersion relation thus becomes For δ = 0 this becomes the dispersion relation of the CW recirculator BBU theory: This describes the case when the recirculated bunches are in the same buckets as the accelerated bunches. When the recirculated bunches are just between accelerated bunches, then δ = 1 2 , For the case that every bucket is filled, this would be an ERL with t b = t 0 . The dispersion relation for ERLs for other δ is less simple than Eq. (20) and has δt b = (n + 1 2 )t 0 in Eq. (18). This occurs when the decelerating and accelerating bunches are not perfectly centered between each other.
The scale is arbitrary.

C. The Threshold Current
For a given positive current I 0 , the values of ω that satisfy the dispersion relation Eq. (18) will in general be complex. If they all have negative imaginary parts, the beam motion is stable. If one of them has positive imaginary part it will be unstable.
For small currents the beam motion is stable. When the current is increased, at some point, one of these ω will become real. At this point the threshold current is reached. The threshold current I th is therefore the smallest current I 0 for which there is a real ω that satisfies the dispersion relation. To find this current, we note that It is therefore sufficient to investigate ω ∈ [0, π/t b ]. Figure 1 shows I 0 (ω) in the complex plain for ω ∈ [0, π/t b ]. The intersection with the real axis that has the smallest positive value yields the threshold current.

III. APPROXIMATE THRESHOLD CURRENT
It can often be justified to linearize in ǫ = ω λ 2Q λ t b , ǫ ≪ 1, which describes a situation when HOM decay is negligible on the time scale of the bunch spacing t b . This applies to linacs when nearly every RF bucket is filled. The smallest |I 0 | in Eq. (18) is obtained when cos ωt b is close to cos ω λ t b , which occurs whenever ω is close to ω λ,n± = ±ω λ + n 2π t b for any integer n. Due to Eq. (21), all these frequencies lead to the same threshold current. We therefore additionally linearize in ∆ω = ω − ω λ,n± ,  The curve is dark blue in the region where T12 sin ωtr < 0 and light green otherwise. The ω of the instability is given by the intersection of the dotted line with a dark blue solid curve for which ∆ω is smallest. Left: A case of T12 sin ω λ tr < 0. Right: A case of T12 sin ω λ tr > 0.
assuming ∆ωt b ≪ 1. This leads to Within the linearization, the phase factors could be combined to e iωtr . This is not done in order to retain the symmetries of Eq. (21) for ω → ω + 2π/t b and for ω → −ω * , i.e. ∆ω → −∆ω * and ω λ,n± → ω λ,−n∓ . Due to these symmetries, the real current close to ω λ,n± is the same for each of these frequencies. Without loss of generality we therefore use ω λ,n± = ω λ and no longer require the symmetries, Since I 0 must be real, ∆ωt b sin ωt r = ǫ cos ωt r leading to the following two equivalent equations, For this formula to describe the threshold current, it is required that I 0 > 0 and therefore T 12 sin ωt r < 0. Figure 2 shows ∆ωt b and ǫ cot ωt r versus ωt b ∈ [0, π] for t r = 6.88t b . The dotted line and the curve have to meet at a region where T 12 sin ωt t < 0, which is indicated by a dark blue curve.
When n r ǫ ≪ 1, i.e. HOM decay is negligible on the time scale of recirculation time, then ω ≈ ω λ in the region where T 12 sin ω λ t r < 0 and one obtains which is the traditional and commonly used approximation in Eq. (1) which had been derived for δ = 0. In the region with T 12 sin ω λ t r > 0 one obtains ωn r ≈ nπ, which can be used in Eq. (25), For n r ǫ ≫ 1, i.e. when HOM damping is substantial on the recirculation time scale, one again uses Eq. (25) with ωt r ≈ (2n ∓ 1 2 )π, where the + or − sign is determined by the sign of T 12 sin ω λ t r , Note that in this case the threshold current weakly depends on t r and can be estimated simply by I 0 = 2ǫ/K|T 12 |.
One can perform these approximations more accurately, for example, by approximating cot ωt r by a line or second order curve. However, in regions where n r ǫ is not much larger or much smaller than 1, simple formulas cannot be found. Figure 3 (top) shows the threshold current obtained with the approximate analytic solution compared with the threshold current that is found by tracking particles for the simple case of one cavity with one HOM and one recirculation loop. Figure 3 (bottom) compares the same tracking results with a numerical solution of the dispersion relation Eq. (18). The data agrees remarkably well with the approximate formula in the region where ∆ω is small, i.e. where I 0 is relatively small. In the region where the threshold current is relatively large, the agreement with the approximate expression is not satisfactory, however.
To find the threshold current with Eq. (18), the smallest positive real value of I 0 for ω ∈ [0, π/t b ] was found by linearly interpolating 100 points in the region ω ∈

A. Instability Growth Rate
We denote the threshold current byÎ 0 and the real frequencyω ∈ [0, π/t b ] satisfies Eq. (18) for this current. When the current I 0 is slightly larger thanÎ 0 , there is one frequency ω(I 0 ) that satisfies Eq. (18) and is close tô ω. It has a positive imaginary part. All other frequencies ω at which Eq. (18) holds and for which thereforẽ V Σ (ω) might not vanish have an imaginary part that is not positive. In Eq. (10) there are therefore exponentially growing terms. For currents that are only slightly larger than the threshold current, the complex frequency ω(I 0 ) can be expanded with respect to ∆I = I 0 −Î 0 , The rise time per current of the instability is thus given by α = ℑ{ dω dI0 }|Î 0 . The dispersion relation Eq. (18) leads to a long formula. However, using the simplified Eq. (22) leads to Provided parameters are not in the region where the curves diverge in Fig. 3 (top), i.e. sin ω λ t r is not close to zero, the following approximate formulas hold for the growth rate: For n r ǫ ≪ 1 one obtains α = 1

IV. MULTIPLE DIPOLE HOMS AND MULTIPLE RECIRCULATIONS
Recirculating linacs with many cavities and several recirculation loops have been considered early on [13,18]. Here we use the same nomenclature as much as possible. The N higher order modes, which can be associated with different cavities, are numbered by an index i. The N p passes through the linac are numbered by an index I.
The horizontal position and momentum that the beam has at time t in the HOM i during turn I is denoted z I i (t) = (x I i (t), p I i (t)). The transport matrix that transports the phase space vector z J j at HOM j during turn J to z I i is denoted T IJ ij and the time it takes to transport a particle from the beginning of the first turn to HOM i during turn I is denoted t I i . The beam is propagated from after HOM i − 1 to after HOM i by .
This equation can be iterated to obtain the phase space coordinates as a function of the HOM strength that creates the orbit oscillations. With the matrix element T IJ ij = (T IJ ij ) 12 one obtains The strength V i (t) of the HOM i is created by all particles that have traveled through that HOM via the integral where I I i (t) is the current at time t that the fraction of the beam has which passes the HOM i on turn I. Combining this with Eq. (32) leads to the following integraldifference equation: Now the approximation of short bunches is used. The current is given at time t by pulses that are equally spaced with the distance t b , This reduces the integral to a sum, where n(t, t I i ) = Max m (t ≥ mt b + t I i ). Computing leads to The second summation starts at m * = n − n(nt b + t L i , t I i ) which can be simplified by writing i > δ I i and 0 otherwise. Shifting the summation index m now leads to and with δ i (I, L) = γ i (I, L) + δ I i − δ L i , which is between 0 and 1, shifting the index n finally leads to the relation The following sum is equivalent to that in Eq. (17), Equation (42) reduces to If a vector V is introduced that has the coefficients V I i , this equation can be written in matrix form, where the matrix M = W(ω)U is determined by Eq. (44). When all electrons are considered to have the speed of light, t I i − t L i does not depend on the HOM number i and we therefore drop this index and obtain the following matrix coefficients: where Top(x) is the smallest integer that is equal to or larger than x. With Kroneckerδ ik this determines the matrices W and U to be For each frequency ω, I −1 0 is an eigenvalue of M(ω). Since the eigenvalues are in general complex, but I 0 has to be real, the threshold current is determined by the largest real eigenvalue of M(ω). The matrix has the properties and it is therefore again sufficient to investigate ω ∈ [0, π/t b ] to find the threshold current. Note that V Np N never appears on the right hand side of Eq. (42) so that the dimension of M can be reduced by one to N × N p − 1. Furthermore the dimension can be reduced when two fractional parts δ I i and δ J j are equal since then V I i and V J j are identical. Note also that for N = 1 and N p = 2 Eq. (42) reduces to the dispersion relation for one HOM in Eq. (18).

A. Instability Growth Rate
The growth rate of the instability is again computed by first obtaining the threshold currentÎ 0 and the real frequencyω for whichÎ −1 0 is an eigenvalue of M(ω). If this is the kth eigenvalue λ k (ω) of the matrix M(ω), then the growth rate of the instability is given by

V. MULTIPLE HOMS IN ONE CAVITY
The presented theory for multiple HOMs and multiple recirculation turns can in general only be evaluated with computers. However, for some simple situation an analytical understanding is possible.
One such case is an accelerator with one recirculation loop and one cavity in which many HOMs can be excited. The (1, 2) matrix elements that refer to transport between HOMs for the same pass are zero, T JJ ij = 0. All matrix elements that describe the recirculation loop are identical, T IJ ij = T 12 if I = J. The matrix elements in Eq. (46) are then given by (51) Equation (45) becomes with δ = δ(2, 1), A summation over the index i shows that j V 1 j can only be nonzero when Comparing this with Eq. (14) shows that one only has to replace w(δ) of the single HOM by i w i (δ) to arrive at the threshold formula for the multi-HOM case. To find the smallest real I 0 that Eq. (53) can produce for real ω, i w i (δ) has to be maximized. The sum is especially large when the denominator of one of the terms is very small, i.e. when cos ωt b ≈ cos ω λ t b . When the HOM frequencies modulo 2π/t b are sufficiently different for the different HOMs, the maximal absolute value of i w i (δ) will be close to the largest absolute value that any of the w i (δ) could have individually. This is due to the fact that all the w j (δ) for j = i are relatively small for frequencies ω for which the denominator of w i (δ) is small.
For HOM frequencies that are sufficiently different in the above sense, the threshold current for several HOMs therefore does not differ significantly from the threshold current of the worst individual HOM. Figure 4 shows how the threshold current changes when the frequency of one HOM is fixed at a small threshold current with | sin ω 1 t r | = 1 and a second HOM frequency is varied. Superimposed is I th if only the second HOM is present. It is apparent that the HOM that would produce the larger threshold-current if it was solely present only influences the threshold current of the pair when the frequencies ±ω 1 mod 2π/t b and ±ω 2 mod 2π/t b are closer together than about ∆ω λ = ǫ/t b = ω λ /2Q λ .
One can draw the conclusion that in the case of many HOMs they do not interact destructively when the frequencies ±ω λ mod 2π/t b are not very close together. Tracking simulations also demonstrate this effect.
One strategy to increase the BBU threshold current is the introduction of HOM frequency spreads between cavities. As an example, Fig. 5 shows the threshold current found by tracking as a function of uniform frequency spread of 20 HOMs in a single cavity. For all three curves (R/Q) λ Q λ is the same. It is seen that for lower Q λ , the  curve begins to saturate for larger frequency spread than for the high Q λ case. It is also seen that the threshold for frequency spread ω λ /2Q λ is similar in all three cases. Both observations are consistent with the above assertion that HOMs do not interfere when they are further apart than ω λ /2Q λ . The fact that oscillations in the Fig. 5 are smaller for low Q λ is consistent with the conclusion that a broader HOM resonance peak should lead to more overlaps and as a result to less pronounced differences in the threshold for different HOM frequencies.

VI. ONE DIPOLE HOM IN TWO CAVITIES
In order to see how two dipole HOMs interact, we will now analyze a set of two HOMs with one recirculation loop, i.e. N p = 2. We abbreviate δ = δ(2, 1) and α r = e iωnrt b . The dispersion relation is then given by where the matrix N is given by (55) The third row is similar to the first row and eliminating it by similarity transformations leads to the 2 × 2 matrix Even though ω 1 (δ) appears in the denominator, the formula for the eigenvalues of this matrix does not contain such a denominator. Therefore the largest real eigenvalue will again occur at a frequency ω for which one of the HOM frequencies satisfies cos(ωt b ) ≈ cos(ω λ t b ). The term w i (δ) and w i (0) of the other HOM can then again be neglected, so that for sufficiently different HOM frequencies mod 2π/t b the threshold current is again approximately determined by the HOM which would have the smallest I th if there were no other HOMs present.
We therefore now assume that the two HOMs are equal, w 1 = w 2 . Now one can perform an approximation analogous to Eq. (23) leading to .
It is interesting to analyze whether the effect of the two HOMs can cancel. A cancellation could occur most naturally when the linac and the recovery loop are mirror symmetric. A mirror symmetry of the linac means that the beta functions of the first pass, going from low to high energy, are the mirror image of those of the second pass, going from high to low energy. An example of such an optics is shown in Fig. 6.
Since the (1, 2) element of the transport matrix between a region with momentum p 0 and a region with momentum p can be written with Twiss parameters as this symmetry leads to T 22 21 = T 11 21 . An additional mirror symmetry of the return arc leads to T 21 22 = T 21 11 . The eigenvalues of the matrix in Eq. (57) then become, T 21 11 ± T 21 12 (T 21 21 + 2e −iωtr T 11 21 ) .
Since there are two solutions to the quadratic eigenvalue equation, to every eigenvalue that is smaller than 1/I th of a single cavity, there exists in general the one that is larger. Therefore two cavities do not compensate their instabilities, but it is possible to decouple the cavities to the extent that the combined threshold current is just as large as that for a single cavity. For this, one has to choose T 21 12 = 0, i.e. the phase advance of the return arc has to be a multiple of π.
Since the kick of a HOM disturbs the beam most at low energy, the first and the last cavity of an ERL are the strongest contributors to BBU. It seems therefore advisable to adjust the phase advance of the arc to a multiple of π also when the linac has more than two cavities.

A. Multiple Recirculation Turns
One could envision a multi-turn recirculating linac as ERL. The beam would pass the same linac N r times to reach its top energy, and subsequently it would be decelerated in just as many turns through the linac. The current in the cavities would be 2N r times higher than the current that is available at high energy. One could therefore conjecture that the BBU threshold current is N r times smaller than for a one-turn ERL. Here we will show that the threshold current can be significantly smaller than that conjecture and in general can be expected to decrease quadratically with N r . Each bunch passes the linac N p = 2N r times and the matrix in Eq. (46) becomes where the lower indexes have been suppressed since only one HOM is considered. The eigenvalue equation thus becomes An approximation equivalent to that leading to Eqs. (23) and (57) results in Comparing the exponent on the right hand side to those in Eqs. (44) and (46) shows that this can be written as Since the right hand side does not depend on L, all the terms V J e iωt J are equivalent and the condition that they do not vanish is where the double sum goes over all pairs of I and J for which I > J. Similar to the condition obtained from Eq. (23) which is analyzed with Fig. 2, the fact that I 0 has to be real entails the condition In regions where I>J sin(ω[t I − t J ])T IJ > 0, we again obtain the approximation that ω ≈ ω λ . The equation for the threshold current of an N r times recirculating ERL I Nr th corresponds therefore to that of the case without recirculation in Eq. (1), A comparison with Eq. (1) shows that this current is smaller than the one-turn ERL by a factor of up to I>J |T IJ |/|T 12 |. This is in agreement with earlier result presented in [21]. Assuming that all matrix elements are of about equal magnitude, the threshold current in an N r times recirculating ERL is therefore in general smaller by about a factor of N r (2N r − 1). This conclusion is consistent with tracking results for microtrons [22] and for two-turn ERL [23]. The scaling in a particular case, however, can be quite different depending on details of the lattice design, e.g. approximate scaling with N r was reported in [24].

VII. CAVITY MISALIGNMENTS
In the derivation above, it was assumed that the bunches travel along the cavities' symmetry axes when the current is below the threshold for BBU instability. When the cavities are misaligned, the beam will excite dipole higher order modes even below the threshold and the trajectory will be disturbed by these modes.
Let us assume that the ith cavity is misaligned with respect to the path adjustment of the Ith turn by x I 0i , leading to the misalignment vector x 0 . The dipole HOMs that are excited by the beam are now not only due to the beam position fluctuation that is produced by the HOMs themselves, but additionally due to the cavity misalignments. The HOM voltages in Eq. (33) are therefore given by With the manipulations that led to Eq. (37) this leads to Following the derivation to Eq. (44) leads to For all oscillation frequencies that are not multiples of the bunch repetition frequency, ω = 2πl/t b , l is an integer, the term x I 0i ∞ n=−∞ e iωnt b vanishes so that the condition for V L i (ω) to be non-zero is the same as for the BBU instability without misalignments x I 0i . Below threshold, the HOM voltage therefore has the following form: The voltage seen by bunch i on turn L is therefore lt L i so that for the time t L i , Eq. (68) can be written as After performing the summation over m one obtains where the superscript ω = 0 means that w(δ) is computed with Eq. (17) for ω = 0. Comparing with Eqs. (44) and (45) shows that this can be written as with W = W(0) and U from Eqs. (47) and (48). The vector of voltages is therefore given by V = [I 0 WU − I] −1 I 0 W x 0 . Equation (32) can now be used to compute the beams distance from the cavity center at each turn and each cavity. One obtains Evaluating this for a single cavity with a single HOM leads to For δ = 1 2 and with ǫ = ω λ 2Q λ t b one obtains There is a current I 0 at which the denominator becomes 0 and the orbit deviation would become very large. The question arises whether this current is larger than the BBU threshold I th or smaller, so that large orbit excursions would present a new kind of instability.
This problem does not only arise for the single HOM case of Eq. (75) but also for the general case of Eq. (77).
Very large orbit excursions x occur for currents for which the matrix inverse does not exist.
The inverse matrix to be inverted is ( where A is the matrix that diagonalizes W(0)U. We therefore see that 1/I 0 for which the orbit gets very large is given by the eigenvalues of W(0)U. These values are naturally smaller than 1/I th , which is the largest eigenvalue of W(ω)U that is produced for any frequency ω.
This proves that the BBU instability always occurs before the orbit excursion becomes very large.

VIII. TRACKING RESULTS
The tracking code BI (stands for beam instability) was developed to perform studies of beam breakup in recirculating linacs [25]. The algorithm models point charge bunch interactions with HOMs in linacs, taking into account proper time delays between the cavities, transfer maps, etc., allowing BBU simulations due to longitudinal, transverse and other higher order modes in a general linac configuration.
The basic algorithm can be summarized as following. The string of HOMs that a bunch sees in its lifetime between injection and ejection points is represented by a list of pointers to the actual cavities. The proper time delays between cavities is also stored for each pointer. E.g. for N HOMs and N p passes, the list of pointers would be N N p long pointing to N HOMs. This approach allows one to represent any recirculation configuration without limitations. As the train of bunches is injected into the structure, the next instance when any bunch sees any pointer is determined, and the HOM voltage in the corresponding cavity is updated. Then, this bunch is pushed to the next pointer where its coordinates are stored, waiting for its turn in time to be the next bunch going through a pointer. This way no bunches end up ahead of time precluding a situation when a bunch sees a cavity with incorrectly updated HOM fields, i.e. causality is properly realized. Furthermore, the algorithm is general enough to allow modeling of the longitudinal instability where timing between different bunches is no longer kept fixed. The practical realization of this algorithm is relatively fast, allowing the tracking of a complete 5 GeV ERL with 300 HOMs for 0.1 ms in less than a minute on an average per-sonal computer. This duration is sufficient to determine the onset of transverse BBU instability in most practical cases.
The output of the code contains amplitudes of HOM voltages as a function of time, which is used to determine the growth rate of the instability by fitting an exponential. Several successive calls are made to the tracking unit to determine the threshold. The length of the tracking time is estimated from Eq. (30) based on the desired accuracy in threshold determination.

IX. CONCLUSION
For dipole HOMs we have derived the BBU theory for arbitrary recirculation times, so that the theory can be applied to ERLs. The resulting equations have been used to find analytical results for (1) multiple HOMs in one cavity, (2) two equal HOMs for a one-turn ERL, and (3) one HOM for a multiple times recirculating ERL. For (1) the numerical observation that it often suffices to include only the strongest of several different HOM was explained and the distance in frequency was derived for which one can consider two HOMs as different, for (2) it was shown that two cavities do not cancel each others instability but that they can be arranged so that they do not add dangerously, and for (3) was shown that the BBU threshold current for an N r times recirculating ERL is roughly up to N r (2N r − 1) times smaller than that in a corresponding one-turn ERL.
Furthermore, a simple method to compute the orbit deviations produced by cavity misalignments has also been introduced. And it is shown that the BBU instability always occurs before the orbit excursion becomes very large.
Several comparisons with tracking data verify the applicability of the theory and of the tracking program. The conclusions should be useful in determining and optimizing the maximum current for the many ERLs that are currently in design and pre-proposal stages worldwide.

X. ACKNOWLEDGMENTS
We thank Joseph Bisognano and Geoffrey Krafft for careful reading of the manuscript and for their useful comments.