Theory of coupled-bunch longitudinal instabilities in a storage ring for arbitrary rf potentials

We present a theory of coupled-bunch longitudinal instabilities that is based on the coupled set of Vlasov equations governing the particle distribution function, and can be applied to arbitrary longitudinal potentials. We find that the coupled-bunch growth rate is given by a dispersion relation that is parametrized by the eigenvalues of the linear (harmonic) coupled-bunch matrix problem. Our theory therefore treats the wakefield-driven source of the instability and the effect of Landau damping together and on equal footing, and also indicates that the stabilizing effects of Landau damping can approximately be compared to the instability growth rates in a harmonic potential that has the same bunch length. We then apply the theory to a weakly nonlinear oscillator and to a quartic potential that is relevant for ultralow emittance storage rings that employ bunch-lengthening systems. In the latter case we find that the theory compares quite favorably with particle tracking simulations for the parameters of the planned upgrade of the Advanced Photon Source.


I. INTRODUCTION
Coupled-bunch instabilities, which may arise when many bunches interact over many turns, can lead to multibunch oscillations, an increase in effective emittance, and a resulting reduction in storage ring performance. Classic results regarding longitudinal coupled-bunch stability can be found in many references, most notably the early works of [1,2] and the subsequent books [3,4]. However, these references (and the many cited therein) are typically restricted to simple harmonic motion in all three planes; in particular, most theoretical approaches neglect the effects of Landau damping that arise when the longitudinal motion is nonlinear. While this approximation is valid for short bunches in single rf systems as is often the case, it does not describe multiple rf systems that are used to lengthen the bunch and increase Landau damping [5][6][7][8][9][10][11], and it fails for the extreme stretching and highly nonlinear (possibly quartic) longitudinal potentials planned for many ultralow emittance light sources (see, e.g., [12]).
Early work describing multibunch stability in the presence of two rf systems tuned to a quartic potential was described in Ref. [13], in which the beam response was expanded via a Sacherer-like modal decomposition, and the coupling between azimuthal modes was neglected. Very general theories of collective instabilities in arbitrary potentials were developed in Ref. [14] and later in Ref. [15]. Our approach is similar in some respects to [14], but we have restricted our attention to coupled-bunch instabilities and long-range wakefields which, in our opinion, makes the final results much easier to use. Our theory can also be thought of as a generalization to the results of Thompson and Ruth [16], who showed how to formulate the coupled-bunch problem as a simple matrix equation. Indeed, our final equation involves the eigenvalues of the same coupling matrix, but these eigenvalues enter as parameters of a dispersion relation for the complex normal mode frequencies. Hence, we incorporate the potentially unstable coupled-bunch motion with the nonlinear effects of Landau damping in a unified and self-consistent framework.
In this way our approach also differs from that of Refs. [17,18], who determined stability by comparing the usual coupled-bunch growth rate to a separately calculated Landau damping rate for the longitudinal potential under consideration. Their resulting stability threshold is generally in good agreement with our approach, but our theory predicts some differences due to the fact that it is self-consistent. We point out these differences for the case where the multibunch instability is driven by a higher-order rf cavity mode, and support these predictions with tracking simulations.
We develop the theory in Sec. II, starting from the set of coupled Vlasov equations and concluding with the coupledbunch dispersion relation (21). We then proceed in Sec. III to apply the theory to three potentials. First, we reproduce wellknown results when the force is linear and the motion is simple harmonic. Next, we solve the dispersion relation for a weakly nonlinear oscillator, finding that while the stability limit is reasonably well described by subtracting the corresponding Landau damping rate from the multibunch growth rate, the shape of the growth rate curve is somewhat changed and the effect of Landau damping decreases with the strength of the instability. Our final example applies the theory to a double rf system tuned to flatten the potential into a quartic oscillator. In this case the general phenomenology is similar to that of the weakly nonlinear oscillator of Sec. III B, but the effects are more pronounced. We then show that the theoretical predictions closely match the results of tracking simulations for the APS-U.

II. THEORY
Our theory investigates the longitudinal motion in phase space, which we parametrize using the position in the bunch z ¼ s − v 0 t and its conjugate momentum p z ¼ −δ ¼ −ðγ − γ 0 Þ=γ 0 . Here, the co-moving coordinate z is defined in terms of the independent variable along the ring s, the particle arrival time t, and the reference velocity v 0 which is nearly the speed of light c, while p z is the negative of the scaled energy deviation from the reference Lorentz factor γ 0 ; our definition of p z ¼ −δ ensures that the Hamiltonian is positive for particles above transition.
We assume that the bunches in the storage ring are wellseparated from each other, so that the distribution function Fðz; p z ; sÞ of each bunch satisfies its own Vlasov equation. Hence, we have for each n, where 0 ≤ n ≤ N b − 1 and N b is the number of bunches in the ring. As indicated, the single particle equations of motion can be obtained from the Hamiltonian Hðz; p z ; sÞ. We will take the Hamiltonian to incorporate the longitudinal single bunch potential (including rf focusing and, potentially, single-bunch collective effects) and the long-range wakefields driven by, for example, higher-order modes in the rf cavities. We write the potential due to the former as V 0 ðzÞ and the latter as V wake ðz; sÞ, so that the Hamiltonian is where α c is the momentum compaction factor (for simplicity we assume that the we are above transition and the slip factor α c − 1=γ 2 ≈ α c ) and we will deal with V 0 ðzÞ later. We describe the long-range potential V wake ðz; sÞ via the wakefield W k ðzÞ, which quantifies the one-turn energy loss of a test particle at z due to a drive particle at z ¼ 0. We will use the fact that W k ðz > 0Þ ¼ 0, which is generally true for ultrarelativistic particles due to causality but may also apply for the long-range wakefield in a low energy ring, and note that the total energy loss for any particle is obtained by summing the contributions from all N b bunches in the ring over all previous turns.
To obtain an explicit expression for V wake ðz; sÞ, we normalize the distribution function F n such that R dzdp z F n ¼ 1, express the number of particles in bunch n as N part n , and write the equilibrium centroid spacing between bunch n and j to be L n;j , with L n;j > 0 if j > n and L n;j ¼ −L j;n if j ≤ n. Then, the potential acting at position z in bunch n due to the other particles in bunch j is given by Here, χ j ¼ e 2 N part j =γmc 3 T 0 is the coupling strength, with T 0 the revolution time in the ring, m the particle mass, and e the fundamental unit of charge. The sum over l gives the contribution of the wakefield driven by F j at all previous turns, and the total long-range potential for bunch n is found by summing over all bunches, At this point our description of the multibunch dynamics has been quite general, and we are still rather far from any solution. To determine the stability of the Vlasov equation (1) under the influence of the long-range potential (3)-(4) we will begin by making two important approximations: first, we will linearize the Vlasov equation about its equilibrium; second, we will assume that the long-range potential varies slowly over the length of one bunch. Both of these approximations have been used extensively in the literature, but we have found a novel way to apply them that results in a set of equations that can be solved rather easily.
We linearize the Vlasov equation by writing F n ðz; p z ; sÞ 1 F n ðz; p z Þ þ f n ðz; p z ; sÞ, whereF n is the equilibrium distribution that is independent of s (i.e., static), while f n is the s-dependent (dynamic) perturbation, and in some sense jf n j ≪ jF n j. Separating out the equilibrium and perturbed distribution results in a similar division of the long-range potential, so that we have V wake ðz; sÞ ¼V wake ðzÞ þṼ wake ðz; sÞ; ð5Þ whereV wake andṼ wake involve a sum overF j and f j , respectively. To make further progress we employ our second approximation that the wakefield varies slowly over the length of the bunch by Taylor expanding W k as [16] If we insert the expansion (6) into the static potential V wake , we find that the first and last terms give rise to an additional energy loss of each bunch that we can compensate for by adjusting the accelerating rf phase and/or voltage. Meanwhile, the second term results in an additional longitudinal focusing force that can be absorbed into V 0 with an appropriate redefinition of the synchrotron frequency.
On the other hand, because R dẑdp z f n ðẑ;p z ; sÞ ¼ 0, only the third term from the expansion (6) contributes to the perturbed potentialṼ wake . Hence, the long-range wakefield due to the perturbation can be written as Having absorbed the equilibrium wakefieldsV wake ðzÞ into the definitions of V 0 ðzÞ, the Hamiltonian now reads Hðz; p z ; sÞ ¼ α c 2 p 2 z þ V 0 ðzÞ þṼ wake ðz; sÞ whereṼ wake ðz; sÞ is given by (7). Now, we assume that we have solved the unperturbed problem, meaning that we have found the canonical transformation from ðz; p z Þ to the action-angle variables ðΦ; IÞ of the Hamiltonian H 0 . Under this transformation H 0 ðz; p z Þ → H 0 ðIÞ and F j ðz; p z Þ →F j ðIÞ, while the potential In terms of the action-angle variables, the set of linearized Vlasov equations simplify to where the oscillation frequency ωðIÞ ¼ c∂H 0 =∂I. Now that we have linearized the coupled Vlasov system, we can isolate the time dependence by defining f n ðΦ; I; sÞ ¼f n ðΦ; IÞe −iΩs=c ð11Þ for complex Ω. In addition, we will eventually be interested in a coupled set of equations for the centroid positions, which we introduce with the notation Using the exponential time dependence (11), the centroid definition (12), and the long-range potential (9), we find that the Vlasov system (10) where it should be understood that the wakefield is to be evaluated at ξ ¼ −ðlcT 0 þ L n;j Þ. We solve the linearized Vlasov system (13) using a method similar to that presented in [14,19,20], which we apply by writing the left-hand side as Then, we multiply both sides by ce −iΩΦ=ω =ω and integrate over Φ from Φ 0 to Φ 0 þ 2π. Sincef n ðΦ þ 2π; IÞ 1 f n ðΦ; IÞ we find that To evaluate the right-hand-side we expand z as a Fourier series in the angle Φ 0 , writing zðΦ 0 ; IÞ ¼ P m e imΦ 0 z m ðIÞ with z −m ¼ z Ã m and z 0 ¼ 0. In this case we can integrate over Φ 0 , so that after cancelling common terms we obtain Finally, we obtain a closed-form expression for the beam centroids hzi n by multiplying both sides of (16) by z and integrating over all phase space. In doing this, we simplify the right-hand-side using Z 2π Then, using χ j ¼ e 2 N part j =γmc 2 T 0 we find that To arrive at (18), the only assumption we made beyond linearization was that the long-range wakefield varies slowly over the length of the bunch. The result is a system of equations for the bunch centroids in terms of equilibrium properties, so that any instability-driven shape distortions are not included (in other words, we neglect any intrabunch coupled motion). This is often a good approximation in rings where the bunch length is much shorter than the rf-wavelength, but in other cases one could consider generalizing our approach to include some number of the higherorder moments hz 2 i n , hz 3 i n , etc., but with a commensurate increase in the complexity of the resulting equations. Now, we will make two additional approximations in order to simplify the interbunch coupling into an easy-tosolve matrix problem. To do this, we first assume that the shape of the equilibrium distribution function is independent of the bunch number, so thatF n ¼F for all n. This assumption allows us to cleanly separate the multibunch coupling from the specifics of the longitudinal potential, but also means that the following analysis may not describe large variations in the bunch profile that can arise in nonuniform fill patterns with passive harmonic systems (for those interested, we indicate how to relax this restriction in Appendix A). Second, we will assume that we can consistently approximate e ilΩT 0 by e ilhΩiT 0 , where hΩi is real and with magnitude equal to the average synchrotron frequency. This approximation is generally applicable to instabilities near threshold, in which case the wakefield results in a small perturbation to the longitudinal motion and the predicted growth rates are much smaller than hΩi.
In addition, there are other situations when setting e ilΩT 0 ≈ e ilhΩiT 0 is valid because the matrix to be defined in Eq. (19) turns out to be approximately independent of hΩi. This occurs, for example, if the long-range wakefield is driven by a higher-order rf cavity mode whose linewidth is much larger than hΩi, as is the case for the APS-U example to be treated and discussed further in Sec. III C 2.
Under the two approximations just mentioned, the coupling between bunch centroids is given by a matrix that can be diagonalized in the usual way. To connect our approach to that of Ref. [16], which assumes a quadratic potential for V 0 ðzÞ, we will introduce the matrix M whose components are given by where σ δ is the rms energy spread, σ z is the rms bunch length, and σ t ¼ σ z =c. At present the prefactor above may seem somewhat arbitrary, but it has been chosen so that M may be computed using codes developed for single rf systems. At any rate, using these assumptions and definitions we can write the set of equations for the centroid positions (18) as The second line in (20) is a dimensionless complex function of Ω that contains the dependence on the longitudinal potential. This part will describe the effects of Landau damping when the dynamics is nonlinear and the oscillation frequency depends on amplitude (action). On the other hand, the first line encompasses the coupling between bunches as given by M. This part only depends on V 0 via the fixed characteristic frequency hΩi, so that it can be treated using the usual methods developed for harmonic potentials. In particular, diagonalizing the matrix M will result in an uncoupled set of equations for the multibunch normal modes. It turns out that this diagonalization can be done analytically for a "uniform fill" in which each bunch has the same number of particles and is equally spaced throughout the ring (i.e., when N part j and L j;jþ1 are independent of j). If this is not the case, one can proceed with a usual eigensolver or with a program like CLINCHOR [21]. Either way, the diagonalization entails finding a matrix U such that UMU −1 is diagonal, so that ðUMU −1 Þ n;j ¼ λ n δ n;j with δ n;j the usual Kronecker delta. The coupled-bunch eigenvector τ then has components τ n ¼ P j U n;j hzi j , and for each associated eigenvalue λ n we get the following equation for the complex frequency Ω:

III. APPLICATIONS TO SPECIFIC LONGITUDINAL POTENTIALS
In this section we will apply the theory developed in Sec. II to several specific longitudinal potentials. To do this in a uniform and concrete manner, we will assume that the equilibrium distribution is exponential in the energy H 0 ðIÞ so thatF In electron storage rings the damping and diffusion due to synchrotron emission naturally drivesF to be a Gaussian function of p z , while this assumption must be evaluated more carefully for proton machines. The exponential-inenergy equilibrium Eq. (22) implies that and our dispersion relation (21) simplifies to A. Harmonic potential of a single rf system For the first example we will assume that the longitudinal potential is given by the lowest-order harmonic potential of a single rf system, so that V 0 ðzÞ ¼ ω 2 s z 2 =c 2 , and the transformation to action-angle variables is given in any number of books (e.g., [22,23]) as Making the transformation (25) results in the unperturbed Hamilonian H 0 ðIÞ ¼ ω s I=c, which in turn implies that the oscillation frequency is independent of action and that the unperturbed distribution function is an exponential function of I, with the average action is given by hIi ¼ σ z σ δ . In addition, the only nonzero Fourier coefficients of the longitudinal position are z 1 ðIÞ ¼ z −1 ðIÞ ¼ σ z ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi I=2hIi p , and the integration over action in the dispersion relation (24) can be computed as Then, we find that the dispersion relation (24) simplifies to Ω 2 − ω 2 s ¼ 2λ n ω s . To convert this into the usual expression, we note that for small growth rates we have Ω ≈ hΩi ¼ AEω s . Taking the positive sign implies that Hence, the nth coupled bunch mode is unstable if ℑðλ n Þ > 0. Note that if we had chosen hΩi¼−ω s the matrix M → M Ã and we would have found that Ω ≈ −ðω s þ λ Ã n Þ. In other words, the complex frequency has the same imaginary part (growth rate) as that in Eq. (28), but the real part of Ω is of opposite sign. This reflects the fact that the centroid position is a real number, but otherwise gives no additional physics, so in what follows we will focus on solutions with ℜðΩÞ > 0.

B. Weakly nonlinear potential
In this example we will expand on our previous results by perturbing the harmonic potential of Sec. III A as We will assume that the quadratic-in-action term is small, so that the dimensionless parameter jbj ≪ 1. This nonlinear term could come from a higher-harmonic rf system employed to introduce additional Landau damping, or could simply originate from the next-order term in the expansion of the pendulum-like rf potential. Our approach here will be to add the new physics of Landau damping which arise when b ≠ 0 in the simplest manner possible. In particular, the H 0 of (29) implies that the particle frequency ωðIÞ ¼ ω s ð1 þ bI=hIiÞ, and the dominant new effect comes from the addition of two poles in the dispersion relation at The pole at I AE is relevant, meaning that it can cross the integration contour along I ≥ 0, when ℜðΩÞ ≈ AEω s , respectively. As argued previously we need only investigate the case with ℜðΩÞ > 0, so that in what follows we only consider the pole I þ . When jbj ≪ 1 the longitudinal physics other than the frequency's dependence on I can be approximated by the b → 0 limit, and we takē Inserting (31) and ωðIÞ ¼ ω s ð1 þ bI=hIiÞ into the dispersion relation (24) results in the following integration: In the first line we have dropped all b ≠ 0 corrections except that which gives the pole at I þ and the associated new physics of Landau damping, while the final Eq. (32) comes from the expectation that Ω ≈ ω s so that Equation (33) defines the dispersion relation of the complex coupled bunch frequency Ω as a function of the eigenvalue λ n . However, evaluating the integral is complicated by the singularity at x ¼ ζ; in particular, integrating along the real line as indicated in (33) results in a dispersion relation that is discontinuous as ζ goes from having a slightly positive to negative imaginary part, a fact that can be seen from the Sokhotski-Plemelj theorem for a < 0 < b with P denoting the principal value. The resulting jump by 2πiζe −ζ in Eq. (33) is unphysical, since the dispersion relation should define a smooth complex function of ζ in terms of λ n . This seeming paradox was resolved by Landau [24], who noted that for a well-posed initial-value problem the Laplace transform is naturally convergent provided ℑðΩÞ > 0, so that as written the dispersion relation (33) must also have ℑðζÞ > 0. Landau then showed that to make sense of Eq. (33) for arbitrary complex ζ requires analytically continuing the integral to ℑðζÞ < 0. Operationally, this means that one must deform the integration contour in Eq. (33) to be always below the pole at x ¼ ζ, in which case the dispersion relation is a smooth function of ζ. The preceding discussion has been a brief presentation of the mathematics behind Landau damping as it applies to multibunch instabilities; good accounts of the physics of Landau damping can be found in, e.g., [3,25]. Once the Landau contour for the integration has been specified, the dispersion integral in Eq. (33) can be evaluated in terms of the exponential integral EiðζÞ as We expect that Landau damping will be unable to contain the coupled bunch instability when the growth rate is larger than the frequency spread, ℑðλ n Þ > jbω s j, in which case we will generally also have jΩ − ω s j > jbω s j. When the instability is in some sense strong these inequalities will be well-satisfied, so that expanding (35) for jζj ≫ 1 yields Solving the resulting quadratic equation for Ω − ω s implies that Hence, we find that when the instability is strong the coupled-bunch mode frequency is shifted by both ℜðλ n Þ and the nonlinear term ∼bω s , while the growth rate is given by ℑðλ n Þ. If we had solved Eq. (35) to next order in jbω s =λ n j we would have found that the reduction in the growth rate due to the nonlinearity scales as jbω s =λ n j; in other words, the effect of Landau damping actually decreases as the coupled-bunch matrix growth rate increases. While this particular fact may be somewhat surprising, it is not surprising that Landau damping is irrelevant to Eq. (37), since the general prescription is that Landau damping can only counter instabilities whose growth rate is of the order of the nonlinear frequency spread jbω s j.
To investigate the dynamics when Landau damping plays an important role, we now consider the case when Landau damping just balances the destabilizing long-range wakefields, so that the coupled bunch mode frequency Ω is purely real and ℑðζÞ ¼ 0. Then, we can identify the Landau damping rate by the imaginary part of λ n for which Eq. (35) is satisfied when ζ is real. We start by assuming that the eigenvalue λ n would give pure growth for a simple harmonic potential, λ n ¼ iν with ν real and greater than zero. Then, solving the imaginary part of Eq. (35) implies that 1 ¼ ζ 0 e −ζ 0 Eiðζ 0 Þ or ζ 0 ≈ 1.347. This in turn means that the collective bunch motion oscillates with a frequency Ω ≈ ω s ð1 þ 1.347bÞ; while plugging ζ 0 into the real part of (35) shows that Landau damping counteracts the matrix growth rate Equation (39) agrees with the physical intuition that the Landau damping rate is of order the frequency spread jbω s j. Interestingly, it is also about 12% larger than that of Ref. [26,27]. Further analysis requires a more complete model of λ n ; for this purpose we will assume that the multibunch instability is driven by a single higher order mode (HOM) with resonant frequency ω HOM and quality factor Q. For simplicity we also assume that ω s ≪ ω HOM =2Q ≪ ω 0 ¼ 2π=T 0 ; the first condition implies that the width of the resonance is much larger than the synchrotron frequency, while the second means that the HOM damping time is much longer than the revolution time so that the HOM overlaps with only one revolution harmonic. In this case the unstable matrix eigenvalue can be approximated by where the dimensionless frequency difference from a revolution harmonic is ϖ ¼ ðω HOM − Nω 0 Þð2Q=ω HOM Þ for integer N. The strength of the instability is proportional to the total beam current I tot and the HOM shunt impedance R s , being characterized by the growth rate We can now use Eqs. (35) and (40) to investigate multibunch stability as a function of the HOM parameters ω HOM and R s . In Fig. 1(a) we plot the coupled-bunch growth rate as a function of the HOM frequency difference from a revolution harmonic for four values of the HOM strength as characterized by the maximum matrix growth rate ν. When the matrix growth rate equals the characteristic nonlinear frequency shift, ν ¼ bω s , the magenta line in Fig. 1(a) shows that the curve of ℑðΩÞ just barely crosses the instability threshold near jϖj ≈ 0, which is consistent with Eq. (39). As the matrix growth rate (i.e., HOM strength R s ) increases, the coupled bunch mode is unstable over a wider range of ω HOM and with a larger growth rate, as expected. Interestingly, the instability curve is not symmetric about the revolution harmonic, and is instead skewed toward HOM frequencies just below Nω 0 . This differs from the usual case of a quadratic potential, which has a maximum growth rate when ω HOM ¼ Nω 0 þ ω s ; since we have assumed that ω s ≪ ω HOM =2Q this small upward shift satisfies 0 < ϖ ≪ 1 and is not noticeable on Fig. 1(a).
To further emphasize these points, Fig. 1(b) compares the weakly nonlinear instability curves with ν ¼ 4bω s (red) and ν ¼ 2bω s (dark blue) to predictions for the same matrix growth rate ν but in a purely quadratic potential (magenta and cyan lines, respectively). The latter simple harmonic oscillator (SHO) curves are nearly symmetric about ω HOM ¼ Nω 0 and have a maximum growth rate ℑðΩÞ ¼ ν when ϖ ≈ 0. This is in clear contrast to the case when the oscillation frequency depends on amplitude. For example, when ν ¼ 2bω s the nonlinearity completely stabilizes the coupled bunch mode when ϖ ≳ 1, but provides much less effective damping when ϖ < 0. At the revolution harmonic ϖ ¼ 0, the nonlinear growth rate is approximately 0.6bω s lower than that of the quadratic potential SHO 2, which is a somewhat smaller reduction than the Landau damping rate of 0.9bω s predicted in Eq. (39). The Landau damping at ϖ ¼ 0 is even less effective when the matrix growth rate doubles to 4bω s : Fig. 1(b) indicates that in this case the nonlinearity only reduces the growth rate by about 0.4bω s . These findings are consistent with the discussion following Eq. (37), in which we noted that the effective Landau damping rate becomes less as ν increases.
In summary, Fig. 1 shows that the coupled bunch growth rate in a nonlinear potential cannot be strictly reduced to a single Landau damping rate; the nonlinearity reduces the growth rate much more significantly when ω HOM > Nω 0 , and becomes less effective as the strength of the instability grows. Nevertheless, we find that the usual method to assess stability by subtracting the Landau damping rate from the predicted matrix growth rate ν does appear to provide a useful quantitative estimate when the predicted instability growth rates are small.

C. Quartic potential of a double harmonic rf system designed for bunch lengthening
We now conclude our examples by calculating the multibunch growth rates for a quartic longitudinal potential, V 0 ðzÞ ∝ z 4 , and by comparing the theory to tracking simulations. The motion in a quartic potential is fundamentally nonlinear, in that as the amplitude decreases to zero so too does the oscillation frequency. Nevertheless, we will find that the multibunch stability in a quartic potential shares the main features with the weakly nonlinear one that we discussed in the previous section: the growth rate curves are skewed toward negative detunings ϖ < 0, and the effects of Landau damping decrease as the strength of the instability increases. Here, however, these effects can no longer be considered as a perturbation of the usual harmonic oscillator.

Instability theory in a quartic potential
A quartic potential may arise in a storage ring by employing an rf system with two or more frequencies.
In the simplest example we imagine having two rf systems, with a main cavity at the fundamental frequency ω rf , and a harmonic cavity at frequency hω rf for some integer h. Then, the total potential is We Taylor expand (42) assuming jhω rf z=cj ≪ 1, and obtain a quartic potential by choosing the parameters V h =V 1 , ϕ 1 , and ϕ h for a given V 1 such that the quadratic and cubic terms vanish while the linear term just cancels the equilibrium energy loss. The resulting z 4 potential is of particular importance for next generation ultralow emittance storage rings including MAX-IV [28] and the APS-U [28,29], which rely on harmonic rf systems to increase lifetime and decrease emittance growth due to intrabeam scattering by stretching the bunch length.
Having made an appropriate choice of parameters such that the linear, quadratic, and cubic terms in the Taylor expansion of Eq. (42) have been removed, we now have V 0 ðzÞ ∝ z 4 . We will find it convenient to parametrize the resulting quartic oscillator in terms of the rms bunch length σ z and energy spread σ δ as follows: Here, ΓðxÞ is the usual Gamma function, and the second line expresses the potential in terms of the longitudinal action I (for details see Appendix B). Using the Hamiltonian (44) we find that so that the longitudinal frequency scales with the oscillation amplitude to the one-third power. The dynamics of a quartic oscillator is fundamentally nonlinear, since the oscillation frequency is proportional to a power of the oscillation amplitude. For this reason it is less obvious what multibunch growth rates are reduced by the effects of Landau damping. From a facility perspective, it is natural to compare the coupled bunch growth rates of the quartic potential due to the double rf system with those of the usual quadratic potential obtained when only the main rf cavities operate (i.e., to compare the dynamics of a double rf system to that of (42) with V h ¼ 0). However, in this case the benefits of Landau damping are countered by the smaller mean synchrotron frequency in the quartic potential, so that the growth rates in a double rf system with flattened V 0 ðzÞ are typically larger than those when the harmonic cavity is turned off. On the other hand, we have shown that the matrix growth rate that enters Eq. (24) is derived from a harmonic potential whose equilibrium bunch length is the same as that of the combined main and harmonic systems. Hence, if one wants to understand the coupled dynamics in terms of a growth rate that is reduced by Landau damping, the appropriate theoretical comparison is between the quartic oscillator and this fictitious harmonic potential with identical σ z . Now, we proceed to calculate the dispersion equation (24) associated with multibunch stability in a quartic longitudinal potential. Doing so requires two more quantities beyond the Hamiltonian (44) and the amplitudedependent frequency (45), namely, the unperturbed distribution functionFðIÞ, and the Fourier coefficients of longitudinal position z. The first of these merely requires finding the appropriate normalization for the exponentialin-energyF given in Eq. (22); as we show in Appendix B, while the Appendix also shows that to a good approximation we need only the lowest Fourier coefficient of z given by Now we insert the quartic quantities (45)-(47) into the dispersion relation (24) and simplify. In particular, we improve matters by changing the integration variable to and by introducing the dimensionless frequency In equilibrium the synchrotron frequency for a harmonic potential satisfies ω s ¼ α c σ δ =σ t , so that ζ ≈ Ω=ω s for this ω s . Using (48) and (49) we find that the multibunch dispersion relation for the quartic potential associated with a double rf system is where BðζÞ comes from analytically continuing the dispersion relation along the Landau contour, and is given by A similar dispersion relation was also derived in [14]. The procedure for determining multibunch stability is essentially the same as that detailed previously: we first determine the matrix growth rate λ by solving the coupledbunch eigenvalue problem for a (fictious) harmonic potential whose bunch length σ t is the same as that of the quartic potential under consideration; then, we insert the λ with largest imaginary part into Eq. (50) and numerically solve the dispersion integral to find ζ ¼ Ωσ t =α c σ δ ; the coupled bunch growth rate is given by the imaginary part of Ω, with ℑðΩÞ > 0 indicating instability. In the next subsection we turn to illustrating and validating this theory using an example from the APS-U, in which we compare predictions of Eq. (50) to those of tracking simulations.

Comparing theory to tracking for the APS-U
To fully understand how we compare simulation results to theory, we begin with an introductory description of the APS-U parameters and our tracking simulation techniques, and then proceed to compare the growth rates predicted from theory to those found in simulation. The APS-U plans to retain twelve 352 MHz main rf cavities that are presently in operation at the APS. These cavities have been extensively characterized, and five HOMs per cavity have been identified that may drive coupled-bunch instabilities. We list these HOMs and their parameters in Table I. In general, the precise values of f HOM are not known and may differ between cavities by amounts greater than the revolution frequency f 0 ¼ 1=T 0 , so that a conservative stability analysis of the full system typically randomizes the values of f HOM for each cavity over the range of one revolution harmonic [30,31]. Doing this for the full contingent of 5 × 12 ¼ 60 HOMs leads to a forest of unstable conditions depending upon the overlap of possibly several HOMs that have T 0 f HOM close to an integer [32]; here, for clarity we will assume that only one 921 MHz cavity HOM is close enough to a revolution harmonic to drive an instability.
The APS-U lattice is a seven-bend achromat based upon the hybrid design of Ref. [33], which was first scaled to fit the APS footprint, then modified to incorporate reverse/anti-bends to lower the emittance [34,35], and finally extensively optimized [36]. The longitudinal lattice and rf parameters used in the simulations are listed in Table II; while the APS-U lattice has continued to evolve, the present (and we expect final) parameters are very similar. Note in particular that the list include a bunch-lengthening rf cavity operating at the fourth harmonic (f h ≈ 1.4 GHz), and that the main and harmonic parameters are chosen to flatten the longitudinal potential to z 4 .
We simulate the APS-U ring by tracking 50 000 particles per bunch in 48 equally-spaced bunches through a number of elements using the tracking code ELEGANT [37]. Particle coordinates are advanced around the ring via the ILMATRIX element, which allows for fast, symplectic particle tracking TABLE I. Parameters of the five potentially problematic HOMs for the APS-U. Note that all of these HOMs satisfy f HOM =2Q ≫ hf s i, so that in a single rf system the matrix growth rate is given by (40) and the strongest instability occurs when f HOM is a harmonic of the revolution period.
(kHz) through a periodic cell including chromatic and amplitude-dependent tunes, beta functions, dispersion, and path-length. ILMATRIX does this by computing a linear matrix for each particle that is determined from the initial particle coordinates and user-supplied parameters including the Twiss parameters, tunes, dispersion, etc, and how these quantities depend on the particle energy (giving chromatic effects through third order in p z ) and on the transverse coordinates (modeling the amplitudedependent tunes and path length). For our purposes here the transverse dynamics are largely irrelevant and are retained only for completeness, while the primary quantity of interest is the momentum compaction α c (the second order momentum compaction is also included in the simulations, but does not affect the results).
In addition, the ELEGANT simulations model the rf acceleration and focusing via two prescribed RFCA elements that are tuned as listed in Table II. This is a simplification of the passive (beam-driven) higherharmonic cavity planned for the APS-U, although this approximation should not significantly affect the instability physics. Finally, the 48 equally-spaced bunches communicate with each other to possibly drive multi-bunch instabilities via the long-range RFMODE element, which models a single 921 MHz HOM with parameters listed in Table I using the fundamental theorem of beam loading and phasor rotation. In the simulation we first wait 5000 turns for the beam to approach an equilibrium, and then slowly introduce the HOM by ramping R s from zero to its given value over 5000 turns. We then simulate the dynamics for an additional 15 000 to 50 000 turns to see if an instability develops, and fit the growth in energy centroid oscillations to an exponential. We originally also tried to fit a "Landau damping rate" by exciting the coupled-bunch motion and then letting it damp, but were unable to drive the marginally stable coupled-bunch oscillation without also reducing the synchrotron tune spread and, hence, the effects of Landau damping.
Having set the stage, it's now time to compare results from our ELEGANT simulations to theoretical predictions for the APS-U. We will make these comparisons using simulations that both omit and then include the damping and diffusion due to synchrotron emission. To ensure that all simulations have the same longitudinal distribution and synchronous phase for identical rf cavity settings, those without damping and diffusion replace the physics of synchrotron emission with a uniform energy loss U 0 for all particles once per turn. Furthermore, since our theory was derived from the Vlasov equation it does not strictly apply when synchrotron emission is included; we will find that approximating the effects of synchrotron radiation by a simple damping works well if the predicted growth rate is relatively large, but becomes less good near threshold when the instability growth rate is comparable to the radiation damping rate.
Our first example compares the theoretical coupledbunch instability threshold for a single HOM to that found in ELEGANT simulations with no synchrotron emission. We vary the strength of the instability by changing the shunt impedance R s , but otherwise use the parameters of the HOM from Table I whose frequency f HOM is close to 3389f 0 ≈ 920.6 MHz; as indicated by Table I, this HOM satisfies f HOM =2Q ≫ hf s i, so that the corresponding coupling matrix M is approximately independent of hΩi.
To validate the theory we will compare the simulated and theoretically predicted threshold ℑðλ thresh Þ, defined to be the minimum value of the growth rate associated with the HOM λ given in Eqs. (40)-(41). This ℑðλ thresh Þ can be interpreted as the effective Landau damping rate, since it corresponds to the expected growth rate for a simple harmonic potential with the same bunch length σ t . We determine λ thresh from the simulations by slowly increasing the shunt impedance in steps, and noting the value of R s for which a coupled-bunch instability is first observed.
The results for this "Landau damping rate" are summarized in Fig. 2(a), where we have scaled ℑðλ thresh Þ by the quantity α c σ δ =σ t ≈ 1.24h2πf s i to make this comparison in terms of dimensionless parameters. We see that the theory well predicts the instability threshold over a wide range of HOM detuning when there is no synchrotron damping. Furthermore, Fig. 2(a) shows that the effects of Landau damping decrease at large detunings and are not symmetric about ϖ ¼ 0; on resonance the theory predicts an instability threshold of 0.19 α c σ δ =σ t , which agrees reasonably well with the Landau damping rate of 0.176 α c σ δ =σ t that is given by Bosch et al. [18].
We continue our comparisons of theory to simulation in Fig. 2(b) with more traditional plots of the predicted growth rate as a function of the detuning from resonance ϖ ¼ ð2Q=ω HOM Þðω HOM − Nω 0 Þ, with Nf 0 ¼ N=T 0 ¼ 3389f 0 ≈ 920.6 MHz. For this we again apply no synchrotron radiation damping and use the 921 MHz HOM parameters of Table I (with fixed R s ). Figure 2(b) shows that the theory and simulation agree quite well over the entire range of instability, and that the predicted growth rates are consistently lower than those of the fictitious simple harmonic oscillator (SHO) with matching bunch length σ t ≈ 51 ps which was used to compute the matrix growth rate λ. Furthermore, Fig. 2(b) also shows that the growth rate curves in the quartic potential are skewed toward negative frequency detunings, much like the weakly nonlinear example of the previous section. In particular, the maximum growth rate occurs when f HOM is below the revolution harmonic such that ϖ ≈ −1=5, which is significantly larger than the average synchrotron frequency hf s i ≈ 0.035ðf HOM =2QÞ. Furthermore, the frequency shift with the largest growth rate is in the opposite direction to that in a simple harmonic potential, since ℑðΩÞ is maximized in an SHO when f HOM ¼ Nf 0 þ f s for integer N. Finally, we close by comparing the simulated and theoretical growth rate curves in the case where we include synchrotron emission in the simulations and theory, the latter of which we model by subtracting the synchrotron damping rate of 49 1=s from the computed growth rate. This is not strictly valid since the physics of synchrotron emission involve damping and diffusion that are described by a Fokker-Planck rather than Vlasov equation, and we FIG. 2. Comparison of theory (blue) to ELEGANT simulations (red) of the predicted and observed multi-bunch growth rate for no synchrotron radiation damping and longitudinal potential V 0 ∝ z 4 . Panel (a) illustrates the effective Landau damping by plotting the threshold matrix growth rate as a function of normalized HOM detuning; the Landau damping rate is ∼ðα c σ δ =σ t Þ=10 ∼ hf s i. Panel (b) plots the simulated and theoretical growth rate for the 921 MHz HOM in Table I, and compares them to the theory of a simple harmonic oscillator (SHO) of same bunch duration σ t ≈ 51 ps.  Tables I and II, i.e., is identical to Fig. 2(b) but for the inclusion of synchrotron emission, which we model in the theory by simply subtracting off the synchrotron damping rate of 49 1=s. The agreement is quite good. Panel (b) shows that the agreement becomes less impressive as one approaches the instability threshold. The red line shows that halving the shunt impedance (and, hence λ) agrees reasonably well; on the other hand, while subtracting the damping time from the theory predicts stability when R s is cut by one-third, the simulations show an instability near resonance. will find that the agreement between theory and simulation suffers somewhat as the synchrotron damping rate becomes comparable to theoretical growth rate.
We plot the growth rate curves including synchrotron emission in Fig. 3. Figure 3(a) plots the growth rate curves from theory and simulation for the HOM parameters listed in Table I, showing very good agreement between the two. Hence, we predict that the APS-U will have to contend with longitudinal coupled-bunch instabilities driven by cavity HOMs whose growth rates may approach 450 1=s. Presently, the APS-U plans to deal with this instability using a combination of a longitudinal feedback system and cavity temperature tuning, in which the rf cavity temperature is adjusted to shift f HOM away from a revolution harmonic. In addition, Fig. 3(a) shows that the physics of synchrotron emission can be approximated by the damping term e −s=cτ z when theoretical growth rate is much larger than 1=τ z .
In contrast to this, panel (b) of Fig. 3 indicates that the simple damping approximation becomes less appropriate as the theoretical growth rate becomes comparable to 1=τ z . The red and magenta lines plot the growth rates obtained from ELEGANT when λ is decreased by reducing the shunt impedance R s by a factor of 2 and 3, respectively, while keeping R s =Q constant. In the former case when R s → R s =2 the agreement to the blue theory line is still adequate, but when the R s → R s =3 ELEGANT predicts an instability growth rate near resonance of ℑðΩÞ ≈ 1=τ z ≈ 50 1=s. This implies that the effects of Landau and synchrotron damping are not additive near the instability threshold, and that accurately predicting the growth rate requires a more sophisticated theory.

IV. CONCLUSIONS
We have derived a theoretical description of longitudinal coupled-bunch instabilities for arbitrary longitudinal potentials. We find that the coupled-bunch growth rate Ω is given by the dispersion relation (24), where λ is the matrix growth rate of an associated harmonic potential having the same bunch length. Hence, the theory combines the wakefielddriven multibunch instability with the effects of Landau damping in a uniform and self-consistent manner. In the limit of simple harmonic motion Eq. (24) reduces to the usual, well-known formula, while for more general potentials its solution agrees reasonably well with the usual approach that separately calculates growth and damping, provided the former growth rate is computed using the same bunch length. We applied the theory to HOM-driven coupled-bunch instabilities, and showed that in a nonlinear potential the growth rate is an asymmetric function of the HOM frequency difference from a revolution harmonic, being skewed towards negative frequency detuning in contrast to the usual simple harmonic theory. Finally, we applied the theory to a quartic potential relevant for the bunch-lengthening at the APS-U, and showed that its predictions agree well with simulations using the tracking code ELEGANT.

ACKNOWLEDGMENTS
We would like to thank M. Borland Our theory as presented becomes seemingly intractable if the background distribution varies across bunches: in general finding Ω requires solving an N b × N b matrix dispersion integral. Here, we use methods developed in Ref. [38] to show that a rather simple generalization to the dispersion relation (21) is possible if the ring is filled with two or more essentially identical trains of bunches. Within each bunch train the backgroundF n ðIÞ can vary according to bunch number, but the approximate periodicity permits analytic calculation of the matrix eigenvalues λ. Such bunch patterns may arise due to rf transients in the main and harmonic cavities when gaps are introduced to combat ion instabilities. For this calculation, we let 1=ðMT 0 Þ be the minimum spacing between bunches (equispaced bunches have M ¼ N b ), and define in terms of which Eq. (A2) can be written as in the sum of the dispersion relation (24).