Equilibrium of an Arbitrary Bunch Train in Presence of a Passive Harmonic Cavity: Solution through Coupled Ha\"issinski Equations

We study the effect of a passive harmonic cavity, introduced to cause bunch lengthening, in an electron storage ring. We derive a formula for the induced voltage from such a cavity with high $Q$, excited by a a sequence of bunches, allowing for arbitrary gaps in the sequence and arbitrary currents. Except for a minor term that can be determined iteratively, the voltage is given in terms of a single mode of the Fourier transforms of the bunch forms, namely the mode at the resonant frequency of the cavity. Supposing that the only wake field is from the harmonic cavity, we derive a system of coupled Ha\"issinski equations which determine the bunch positions and profiles in the equilibrium state. The number of unknowns in the system is only twice the number of bunches, and it can be solved quickly by a Newton iteration, starting with a guess determined by path-following from a solution at weak current. We explore the effect of the fill pattern on the bunch lengthening, and also the dependence on the shunt impedance and detuning of the cavity away from the third harmonic of the main accelerating cavity. We consider two measures to reduce the effects of gaps: 1) distribution of the gaps around the ring to the greatest extent allowed, and 2)"guard bunches"with higher charges adjacent to the gaps, compensating for the charge missing in gaps. Results for parameters of the forthcoming ALS-U light source are presented.


I. INTRODUCTION
In electron storage rings the phenomenon of Touschek scattering often limits the lifetime of a stored beam [1]. This is the aspect of intrabeam scattering in which small transverse momenta are transformed through Coulomb scattering and a Lorentz boost into one large and one small longitudinal momentum in the lab frame, sending both particles outside the momentum aperture of the ring. The effect may be counteracted by reducing the charge density in the beam. One way to do that is to increase the bunch size in the longitudinal direction. This can be done by adding a cavity with resonant frequency close to a low harmonic of the main r.f. frequency, say the third or fourth harmonic. This is often called an HHC (higher harmonic cavity).
The quadratic potential well of a usual r.f. system can be turned into a quartic potential well, by arranging the HHC so as to zero the second and third derivatives of the effective well [2,3]. This condition, often referred to as "ideal", results in a flat-top equilibrium bunch profile with substantial bunch length increase, say by a factor of four or more in cases of interest, and an increase of the Touschek lifetime by a comparable factor. However, the flat top is not necessarily the best configuration, since a further lifetime improvement can be achieved by "over-stretching", which causes the appearance of two peaks in the bunch profile. This must not be carried too far, however, since eventually the average lifetime will degrade rather than improve with over-stretching.
The harmonic cavity may be passive or actively excited, but a natural first step is to consider the less expensive passive option. Our discussion is for the passive case, but our methods could be adapted to the active system. In the passive case, the field induced in the cavity by a bunch train depends strongly on the fill pattern. If the beam has a uniform fill pattern, e.g. all r.f. buckets are filled or all the bunches are separated by a fixed number of empty r.f. buckets, there exists a beam equilibrium with all bunches having the same profile (possibly of the flattop form if the ideal HHC settings are met). However, if there are significant gaps between bunch trains (or a long gap following a single-train beam) the quality of the beam equilibrium can be compromised. Instead of uniform charge distributions along the train, one then sees a variation of the bunch form and centroid along the train.
This may cause severe limitations to the effectiveness of the HHC system, either because of the resulting uneven lifetime or/and because of interference with the functioning of the machine feedback systems used for beam stabilization, and may prevent the attainment of the desired bunch lengthening.
There are several reasons for the presence of gaps in the bunch train. Historically, gaps have been needed for ion clearing. Another demand arises from the requirements of synchrotron light users, who may need different fill patterns for different types of experiments.
Experiments needing precise timing of x-ray pulses generally require more gaps than those asking for high brilliance. In the ALS-U, gaps are needed for on-axis injection from the accumulator ring [4].
In this paper we present a robust and efficient method to evaluate the beam equilibrium for arbitrary HHC settings and beam-fill patterns. Our approach, entailing the numerical solution of a system of non-linear algebraic equations, extends the method introduced in [34] for the determination of single-bunch Haïssinski equilibria in the case of short-range wake fields. It is much faster than macro-particle based methods and, we believe, an improvement on the method recently introduced by T. Olsson et al., [31].
Our immediate objective is to study the effect of the fill pattern on the bunch densities in the equilibrium state. While this is a useful first step with rewarding practical implications, e.g. offering guidance on choice of the HHC design parameters, our final goal is to understand the threshold in current for an instability, and the time-dependent behavior beyond the threshold.
The widened potential well has some benefits: the reduced peak bunch current and increased longitudinal tune spread may lead to the damping of certain instabilities. However, other instabilities may be induced, either through the fundamental or higher order modes of the HHC [23][24][25][26], or by possibly aggravating the effect of higher order modes in the main cavity [27]. The method presented here is an essential ingredient toward the application of mode-analysis techniques to the study of beam stability when HHC's are present.
Beside reports on specific projects as cited above, there are several papers which discuss the issues that concern us in a more or less general way, through theory, simulations, and measurements. Byrd and Georgsson [3] and Hofmann and Myers [2] treated the situation without HHC beam loading (i.e., without the cavity wake field), which is the starting point for the present work. Towne [28] studied stability of stretched bunches in the presence of a broad band impedance together with a high-Q resonator, using Vlasov-Fokker-Planck simulations and measurements at NSLS-VUV. Byrd, De Santis, Jacob and Serriere [29] initiated the study of the impact of gaps in the bunch train. They used the term "transient beam loading", which several authors have adopted. (Since a transient effect is usually thought of as short-lived in time, not the case here, "inhomogeneous beam loading" might be a more descriptive term.) A direct antecedent of our work is the paper of Tavares, Andersson, Hansson, and Breunlin [30] who were concerned with self consistency in the equilibrium bunch densities. The study of this topic was continued by Olsson, Cullinan, and Andersson [31] who developed an iterative scheme to find the equilibrium charge densities. Bassi and Tagger [13] investigated the option of a super-conducting HHC, invoking self-consistent simulations and emphasizing the importance of beam loading in the main accelerating cavity for a full picture.
The content of the paper is as follows: Section II describes our choice of coordinates and the description of the bunch train.
Section III and Appendix A review the equations of motion.
Section IV states the primary formula for the voltage induced by the harmonic cavity, then Section V notes that the induced voltage can be expressed in terms of an effective wake potential, which is represented by a compact formula that is the basis for further work.
Section VI goes on to find an explicit formula for the induced voltage from an arbitrary bunch train, which is in terms of the Fourier transforms of the bunch forms at the resonant frequency of the harmonic cavity.
Section VII states the Vlasov-Fokker-Planck equation, and shows how its steady state solution is given by the solution of coupled Haïssinski equations. Section VIII shows that the mean energy transfer in the equilibrium state is exactly equal to the energy loss per turn. Section IX calculates the integral of the induced voltage, to get the potential wells for the Haïssinski system. Section X describes a Newton iteration for solution of the Haïssinski system, while Section XI gives the associated Jacobian matrix, and Section XII shows how to follow the Newton solution as a function of current.
Section XIII presents numerical results for the parameters of ALS-U and a comparison to a macro-particle simulation. Section XIV estimates Touschek lifetimes as a function of the cavity detuning.
Appendix B discusses the perturbation of the synchronous phase due to the harmonic cavity, and reports that there is no necessity to base the coordinate system on the perturbed phase. Appendix C explains how our general formula for the induced voltage reduces to a known formula in the case where all bunches are identical.

II. CHOICE OF VARIABLES AND DESCRIPTION OF BUNCH TRAIN
Synchrotron motion in a storage ring can be described in terms of the longitudinal coordinate z = β 0 ct − s, the distance to the reference particle. Here s measures position in the laboratory as arc length along a reference trajectory, and the reference particle has position s 0 = β 0 ct at time t. Particles leading the reference particle have z < 0. The opposite sign convention is often adopted, indeed in our own papers.
For a single bunch, z is familiar as the "beam frame coordinate", which is suitable as a phase space coordinate for equations of motion and the Vlasov equation. In the case of many bunches, z is a convenient global coordinate for description of the total charge density, and merely by adding constants to z we can construct local beam frame coordinates for all the bunches. Moreover, z has the convenient property of being proportional to s at fixed t and proportional to t at fixed s. Thus if we wish to demonstrate periodicity in s at a fixed time we have only to demonstrate periodicity in z.
We consider a sequence of n b bunches , giving a total charge density of the form where λ 1 is the wavelength of the main r.f. cavity, and C is the circumference of the ring.
The m j are non-negative integers specifying the filled r.f. buckets. Without loss of generality we take m 1 = 0; then m j ≤ h − 1, where h is the harmonic number, equal to the maximum number of bunches, and hλ 1 = C. We take ρ j (z)dz = 1, and define ξ j as the ratio of the charge in bunch j to the average bunch charge. The leading bunch in a train, having the most negative z, has the highest bunch index, j = n b .
The bunch profiles ρ j (z) are time-independent, since we are concerned with the equilibrium state, and are initially unknown functions to be determined by the condition of equilibrium.
The infinite sequence in (1) is intended to mimic the periodicity of the charge density in a circular storage ring. We have ρ tot (z + C) = ρ tot (z), so that at fixed t the density is periodic in s with period C. At fixed s it is also periodic in t with period C/βc. The idealization of supposing that the charge pattern exists for all t ∈ (−∞, ∞) is justified, given the large storage times of typical machines.
The total voltage seen by a particle at arbitrary z (at an arbitrary distance from the reference particle) is taken to be where k 1 = 2π/λ 1 . In the model to be explored, the induced voltage V r comes only from the lowest mode of the passive harmonic cavity, as excited by the bunch train. The relation of φ 0 to the synchronous phase, the phase at which the cavity supplies the mean energy lost per turn, will be discussed presently.
We define z j , the argument of the density ρ j , as Then by (2) the total voltage as a function of z j is since the first term in (2) is periodic in z with period λ 1 .

III. EQUATIONS OF MOTION
The usual equations of motion for a single particle, subject only to applied r.f., describe oscillations in a potential well with minimum at the location of the synchronous particle.
Since the harmonic cavity broadens the well and shifts its minimum, a natural step would be to modify the equation of motion so that it describes oscillations about the shifted minimum.
On the other hand, this might be an unnecessary complication if the shift is sufficiently small.
The coordinate of the unperturbed problem might provide a perfectly accurate description, even if it is not the distance to the minimum.
We first recall the derivation of the standard equations for a single particle with only applied r.f. We first derive difference equations, referring to to changes over a full turn, and later replace them by differential equations, since the changes are very small. The salient variable of interest is the phase φ of the applied r.f. at the time that the particle crosses the accelerating cavity.
At the n-th turn the r.f. kick at phase φ n restores the energy loss U 0 of the previous turn and also changes the energy of a generic particle from E n to E n+1 : The synchronous phase φ 0 is that for which the energy supplied is exactly U 0 : For stable motion this angle should be in the second quadrant, with the square root defined to be positive. Defining where E 0 is the nominal energy of the ring, we write (5) as The change with n of (∆φ) n depends on the revolution frequency, which in turn depends on (∆E) n , these dependencies being linear to a good approximation. Invoking the definition of the momentum compaction factor α we show in Appendix A that We wish to follow the trajectory of z, namely which is related to the trajectory of φ as follows: hence (∆φ) n = k 1 z n .
The sign in (12) is correct: if φ(t) > φ 0 at time t when the particle arrives at the cavity, it has arrived later than the reference particle, which is to say that s(t) < β 0 ct.
Approximating the difference equations (9) and (10) by differential equations, with dt = T 0 and δ = (E − E 0 )/E 0 , and applying (13) we have Here T 0 = C/β 0 c is the nominal revolution time of the ring. In replacing (9) by (14) we have equated T 0 with the time between successive arrivals at the cavity, but this is correct at best in an average sense, because different particles have different revolution times. This approximation is not usually acknowledged in textbook treatments of the problem.
The generalization of (14) and (15) to account for many bunches and the harmonic cavity is obtained by invoking the total voltage (2) and the replacements z For this we note that the required relation (10) is derived in Appendix A with allowance for the presence of V r . The derivation requires V r (z) = V r (z + C), which is assured in the following formalism.
A new feature is that φ 0 is no longer the synchronous phase, since the induced voltage V r causes an additional energy increment that must be taken into account. We are nevertheless free to choose φ 0 according to (6) and (7), and we shall indeed make that choice. Some nearby value could do as well.
In order to clarify the impact of the shifted synchronous phase, we have also carried out a calculation with the coordinate system shifted accordingly. We conclude that there is no need to work in such a system. This issue is reviewed in Appendix B.

IV. PRIMARY FORMULA FOR THE INDUCED VOLTAGE
At an arbitrary z the induced voltage from the harmonic cavity will be Here W is the wake potential of the cavity, which for sufficiently large Q has the form In this formula ω r = k r c is the circular resonant frequency of the lowest mode of the cavity, R s is its shunt impedance, Q its quality factor, and θ(z) is the unit step function, equal to 1 for z ≥ 0 and 0 otherwise. The θ function is an expression of causality.
The expression (18) satisfies the obvious requirement that V r be periodic with period C.
To see that, evaluate V r (z + C) by changing the integration variable to z = z − C and the summation variable to p = p + 1.
We suppose that the support of any ρ j (z), the region in which it is non-zero, is much less in extent than λ 1 , a condition that is satisfied in any ring of interest for this study.
To proceed it is convenient to translate the variable of integration and reverse the order of integration and summation, so that the formula (18) takes the form To give an idea of typical parameters for the following work, we list in Table I a tentative set of parameters for ALS-U, the forthcoming upgrade of the Advanced Light Source at Lawrence Berkeley National Laboratory.

V. EFFECTIVE WAKE POTENTIAL
The sum over p in (20) can be thought of as an effective wake potential W(z) which is to be convolved with an effective charge density ρ(z), defined as follows: Then the induced voltage may be expressed as Applying (19) and expanding the cosine by the double angle formula we have The θ function requires p ≥ −z/C, but since p is an integer that means where x denotes the ceiling of x, which is the smallest integer greater than or equal to x.
Expressing the sine and cosine of pk r C in terms of exponentials, we find that where with It is convenient to define real polar variables (η(k r ), ψ(k r )) such that Then from (26), (27) and (28) we have Substituting in (25) and applying the double angle formula in reverse we have where The function χ(z), plotted in Fig.1, is periodic with period C and has a sawtooth form, with its value C at jumps defined by the limit from the left. It follows that W(z) and the induced voltage V r (z) defined by (22) are periodic with period C.
In (30) we have an appealing, compact formula for the effective wake potential, which will lead to the induced voltage after a straightforward evaluation of the integral (22). One must keep in mind that the integrand in (22) has a jump at z = z owing to the jump in An important feature of W(z) is its behavior as a function of the detuning parameter ∆k/k 1 , where for a 3rd harmonic cavity ∆k = k r − 3k 1 . Figure 2 shows the real and imaginary parts of the function ζ = η exp(−iψ) for a typical case. The function resembles a Lorentzian resonant line form, but in fact differs substantially from an actual Lorentzian.
The half-width of the peak in the real part is roughly 3/(2Q). Later we shall find that a true Lorentzian with that half-width occurs in the case of a complete fill with h identical bunches. Since one can show that ζ approaches a true Lorentzian as ∆k tends to zero, the maximum of η is exactly at ∆k = 0.
We are now in a position to compute the induced voltage from (22). The density ρ, defined in (21), is zero except for n b isolated peaks, the bunch profiles. We define the interval Ω i which is to contain the support of the i-th bunch, much shorter than the main r.f. wavelength: This is just to say that the support in terms of the beam frame coordinate z i is within the region |z i | ≤ Σ. Note that the elements of Ω i are close to z = −m i λ 1 and therefore decrease with increasing i. Because of the stated restriction on Σ, no two of the Ω i can intersect: Some numerical experimentation may be needed to find an appropriate and economical value of Σ.
We note that V r (z) need be evaluated only at z within the various Ω i , since the collective force enters the dynamics only in those regions, through the Haïssinski or Vlasov equations.
Also, ρ(z ) is non-zero only for z in the same sets. It follows that the function χ(z − z ) in (30) takes on only two values: This follows from the fact that |z − z | < C. Regardless of the fill pattern, |z − z | cannot be greater than a number C − λ 1 + O(Σ). (For instance, if we have only two buckets, C = 2λ 1 , then the distance between a particle in one bucket and a particle in the other is λ 1 plus a quantity of order Σ .) Thus which implies (34).
Let us evaluate V r (z) for z ∈ Ω i , and divide the terms into three groups, those for which z < z , those for which z > z , and the one "diagonal" term in which both z < z and z > z can occur. Thus where The sums are regarded as empty when the lower limit exceeds the upper.
In each term of (37) and (38) we change the integration variable to z j = z + m j λ 1 , expand the cosine by the double angle formula, and recognize the resulting integrals as real and imaginary parts of a Fourier transform at k r . Then we find Carrying out a similar calculation for v > , we may write the sum of the two terms as The factor exp(k r z/2Q) is very close to 1 in the present application. By Table I we have while |z| is at most about 4 · 10 −2 m in (41). Nevertheless, we shall not replace this factor by 1, since we expect to apply our formulas to cases with low Q in future work.
The term v d of (39), representing the force on a bunch due to the field that it itself The Fokker-Planck term on the right hand side is where ω s is the circular synchrotron frequency and t d is the longitudinal damping time.
We seek an equilbrium in which ∂f i /∂t = 0 and f i has the factored Maxwell-Boltzmann Under this hypothesis the Fokker-Planck term vanishes and the spatial density ρ i must By separating variables and integrating we see that a solution must have the form where A i is a normalization constant, just the integral of the numerator over [−Σ, Σ].
Noting that cT 0 = C, we introduce the definitions Recalling (47) we then have (48) expressed as where We refer to U i as the "potential" for the i-th bunch, even though it has the dimension of an energy times a length. The system of equations (50) for i = 1, · · · , n b will be called the coupled Haïssinski equations.

VIII. MEAN ENERGY TRANSFER IN THE EQUILIBRIUM STATE
According to (16) the power transferred to a single particle with coordinate z i in the i-th bunch is The mean value of the power over the equilibrium distribution is obtained from (47) as thus for every i, The first term on the right hand side is the mean energy supplied by the external r.f., while the second term, which is negative, represents the mean energy lost to the harmonic cavity per turn. We automatically have energy balance, on the average, in the equilibrium state.

IX. INTEGRAL OF THE INDUCED VOLTAGE
To express the integral in (51) we define S and C as follows: The large-Q approximation stated here on the right is not used in our code, since we wish to be set up for later applications with small Q Applying this in (42) and we find It remains to calculate with v d from (39). After changing the integration variable in (39) We can avoid the double integral in (59) through an integration by parts. After applying the double angle formula to the cosine, one of the terms comprising (59) takes the form Now in a partial integration the factor exp(−k r ζ/2Q) cos(k r ζ + ψ) is integrated by applying (55) , while the u-integral is differentiated. Proceeding similarly with the other terms, we eliminate all double integrals.

X. SOLUTION OF COUPLED HAÏSSINSKI EQUATIONS BY NEWTON'S METHOD
Let us multiply (50) by exp(k r z i /2Q) and then take the Fourier transform, as in (41).
This yieldŝ where A i is defined in (50) and If the diagonal term in v d were known, the real and imaginary parts of (62) constitute 2n b equations in the 2n b unknowns Reρ j , Imρ j . Defining a notation for the diagonal term, we write (62) more briefly as where F andρ are complex column vectors with n b components, in which For given u d we try to solve (65) by the matrix form of Newton's method, namely Here (∂F/∂Reρ, ∂F/∂Imρ) are complex matrices with elements (∂F i /∂Reρ j , ∂F i /∂Imρ j ).
Lacking any better choice, we begin the process withρ (0) obtained from Gaussians, all with the nominal bunch length: Note that we could not use the direction-independent complex derivative ∂/∂ρ i , since U i is not an analytic function ofρ i , being always real.
To account for the diagonal term we adopt the simple device of computing u d in (67) from the previous Newton iterate. That procedure yields a convergent scheme, and shows that the contribution of the diagonal term is negligible, at least in the present case of a high-Q cavity. If our scheme is later applied to a low-Q case, a more sophisticated method might be needed to determine v d .

XI. EXPRESSION OF THE JACOBIAN MATRIX
Since the exponent U i is linear in the unknownsρ j , it is not difficult to write down the Jacobian, the matrix of the partial derivatives that appear in (67). One must not forget the derivatives of A i , which are essential to ensure that the final ρ i (z i ) are automatically normalized to have unit integral. The complete 2n b × 2n b Jacobian in block matrix form is where with

XII. CONTINUATION IN CURRENT
In contrast to experience with the Haïssinski equation for a single bunch, we shall find that the Newton iteration (67) beginning with (68) does not converge at the desired design current. This must be because the solution at full current deviates extremely from the unperturbed Gaussian, whereas the deviation is relatively small in the single bunch case.
At small current the Jacobian is nearly diagonal and positive definite, since the offdiagonal terms have a factor eN . This augurs well for the success of the Newton iteration at sufficiently small current. It then seems reasonable to get a solution with small current, then take that solution as the starting point for a Newton iteration at a somewhat higher current. If the second iteration converges we can perhaps repeat the process several times to reach the required large current.
The calculation could be made more efficient by extrapolating linearly in current after each successful iteration. This should allow a larger increment in current. Let us define a convenient current parameter such as I = I avg , the average bunch current. Expanding the notation to include I-dependence, and suppressing reference to u d , we write (65) as and differentiate with respect to I to obtain The solution of this linear system for dρ/dI affords the linear extrapolation ρ(I + ∆I) =ρ(I) + dρ(I) dI ∆I .
The extrapolation is not a costly step, since the Jacobian matrix in (78) is already known from the previous Newton iteration.
It is helpful to redefine the unknown asρ = Iρ, so that the current appears linearly in the transformed version of (77), namelỹ Now even the right hand side of the equation to be solved, −∂F /∂I, is an integral already computed during the previous Newton iteration.
Our continuation procedure is an example of a general method for solving nonlinear problems, called path following or executing a homotopy [35]. One follows a known or easily computed solution as a function of a parameter, which could be multidimensional. The continuation could stall or display a bifurcation if a singularity of the Jacobian were encountered.

XIII. NUMERICAL RESULTS FOR PARAMETERS OF ALS-U A. Complete train without gaps
We first present results for a complete train without gaps, for which n b = h = 328, at the nominal current of 500 mA. The calculation starts at a current of 150 mA and proceeds to the desired current in three equal increments, by means of the algorithm of the previous section. The convergence criterion for a Newton iteration is in terms of a sum of normalized residuals of the equation (65) that is to be satisfied:   The shunt impedance and detuning in Table I were chosen to maximize the bunch lengthening at the nominal current, while keeping a flat top in the density. At higher impedance or lower detuning one can achieve a larger bunch lengthening, but at the expense of getting a density with two maxima. The bunch in this situation is sometimes described as being "over-stretched". Figure 4 shows the effect of decreasing the detuning, and lists the corresponding r.m.s bunch lengths. It is important to note that all bunch forms turn out to be the same, merely by putting n b = h, even though the equations contain no explicit constraint that they be the same.
This is gratifying and as it should be by physical intuition, but the mathematical mechanism for it to happen is somewhat obscure.

B. Train with a single gap
Next we consider a train with a single gap of 44 empty buckets, thus 284 filled buckets in a row. The average current, i.e., the total charge divided by the revolution time, is taken to be the same as before, corresponding to the individual bunch charge being larger by a factor 328/284. With an initial average current of 150 mA, increased to 500 mA in three steps, the convergence is even better than in the previous example. In Figure 5 we show representative bunch forms near the front, middle, and end of the train. Each bunch is given as a function of its beam frame coordinate z j = z + m j λ 1 . There is much less bunch lengthening than in the complete fill, and a large centroid shift varying linearly along the train. In Figure 6

C. Train with distributed gaps
The sharp reduction in bunch lengthening induced by a single gap leads to the idea of distributing the empty buckets around the ring as much as possible [33]. This has a chance of resembling more closely the complete fill. For ALS-U the minimum acceptable gap consists of 4 empty buckets, since a gap of 10ns is required to accommodate the rise and fall times of the fast kicker that does on-axis injection from the accumulator ring. With such gaps we need 9 trains of 26 bunches and two of 25 to account for 328 buckets total: 9 × 26 + 2 × 25 + 11 × 4 = 328.
We consider the Case C2 of Ref. [33], in which the two trains of 25 are as far apart as possible. This was found to be slightly more helpful than putting those two side-by-side. Figure 7 shows the result for 6 bunches out of a train of 26, including the initial and final bunches. Figure 8 shows the ratio of the r.m.s. bunch length to the natural bunch length vs. bunch number, while Figure 9 displays the centroid vs. bunch number. Fortunately, the average bunch lengthening now has a value near the case of the complete fill. Furthermore, the big centroid displacement of the single gap case is gone. There is a small and linear centroid displacement along each sub-train, but its magnitude is similar to that of the complete fill.
Although maximal distribution of the gaps is a step in the right direction, it leaves us with a strong variation of bunch form along the train and some highly skewed charge distributions.
We should then look for further means to imitate the complete fill as much as possible.
D. Guard bunches to compensate the damage from gaps.
If bunches at or near the ends of a train are given greater charge, enough to equal the missing charge due to the gap, the inner bunches may feel less perturbation from the gap.
This idea was advanced by Byrd et al. [29] in 2002, and re-invented at the Argonne APS in Only the guard bunches differ from this pattern, as is seen in Fig.10.
Such highly intense guard bunches could suffer a microwave instability or have a reduced lifetime, or be undesirable for their impact on the synchrotron light pattern. It therefore becomes interesting to distribute the guard charge over several bunches. As an example we take 4 guard bunches at the beginning of the train and 4 at the end, each with 50% more charge than the inner bunches (ξ = 1.5). As is shown in Figures 11-13, the inner bunches again are flat topped, while there is a gradual transition in the guard sequence from the end bunch form to the inner.
Rather than a uniform distribution of charge in the guard segment, one could try some kind of taper, for instance a power law with arbitrary exponent, ξ(j) = 1+aj −b , j = 1, 2, · · · n g , ξ(j) = 1+a(n t +1−j) −b , j = n t , n t −1, · · · , n t −n g +1 , for a train of n t bunches with n g guard bunches at either end. We might try a sharply peaked distribution with large b in order to imitate the case of a single guard bunch but with less peak charge. For instance, putting b = 1.7 and n g = 13 in a train with n t = 26 we get the result of Figure 14. The bunch population ξ(j) of (82) is plotted in Fig.15. With 30% less charge in the end bunches we get a pattern very similar to that of a single guard bunch, in that most of the interior bunches are close to flat topped, and have markedly less charge than the four guard bunches with uniform population as considered above ( Figures   11 -13). E. Comparison to a macro-particle simulation We have applied the code elegant [37] to make a macro-particle simulation for comparison to results of the present method. This was part of an exploration of parameter space, and the parameters are different from those of Table I in the following choices: α = 2.07 · 10 −4 , σ δ = 1.14 · 10 −3 , σ z0 = 4.43 mm , U 0 = 3.29 · 10 5 eV R s = 10 6 Ω , Q = 1.67 · 10 4 , δf = 2.27 · 10 5 Hz .
Also, in (16) we take eV sin φ 0 = (9/8)U 0 , following Eq.(B10) in [25]. The simulation used the cavity wake field description provided by elegant, and was done with 10000 macroparticles per bunch. The fill pattern is that of Section XIII C, with 284 bunches and distributed gaps of 4 buckets each.
The agreement is good enough both to provide a check on our semi-analytic scheme and to affirm the viability of a macro-particle simulation.

XIV. INCREASE IN THE TOUSCHEK LIFETIME
Following Refs. [3] and [38] we note that in Piwinski's analysis [1] the Touschek lifetime is inversely proportional to the integral of the square of the longitudinal charge density. Under bunch lengthening the lifetime should then increase by a factor where ρ and ρ 0 are the charge densities with and without the action of the harmonic cavity.
The approximation in (84) consists of neglecting the ratio of momentum acceptances for the two cases, which Byrd and Georgsson [3] judge to be of order 1. Taking ρ 0 to be a Gaussian with the nominal (zero current) bunch length, we compute R by (84) for the case of distributed gaps, as in Section XIII C, without guard bunches. With smaller detuning and consequent over-stretching we get a further increase in lifetime, but with increasing variation along the train. At the smallest detuning of Fig.19, δf = 185 kHz, any advantage of over-stretching is gone, since for half of the train the lifetime is smaller than for δf = 235 kHz.

XV. CONCLUSIONS AND OUTLOOK
We have described an effective scheme to determine the equilibrium state of an arbitrary bunch train, subject to the wake field from a passive harmonic cavity in its fundamental mode. The calculation proceeds by an iterative method with extremely robust convergence.
The computation time is negligible, and the results agree with macro-particle simulations, which are much heavier calculations and also much more noisy.
The quick computation allows a convenient exploration of parameter space, and in particular an examination of schemes to counter the bad effects of gaps in the bunch train.
We have seen that by distributing the empty buckets around the ring as much as possible the bunch lengthening and centroid displacement can be made comparable to those of the complete fill. Although there is then considerable deviation from the flat topped distribution achieved in the complete fill, that turns out not to harm the Touschek lifetime. Also, most of the bunches can be given a flat topped form by invoking guard bunches adjacent to the gaps.
We have adopted a minimal physical model, with the only induced voltage (wake field) coming from a single resonant mode of the harmonic cavity. With this we could illustrate the power of a new technique in the simplest way. The next step toward a realistic model should be to include the induced voltage from the main accelerating cavity (beam loading).
Since our formalism allows any number of resonators, this is a straightforward extension. In fact, we have revised the code to include the main cavity, and have found that the iterative solution works as well as before, with only a factor of two increase in CPU time.
Another refinement that could be significant is to include the effect of the usual short range wake fields from vacuum chamber corrugations. The magnitude of the effect on bunch lengthening can probably be judged by invoking a broad band resonator model of the machine impedance, which is normally applied with a Q of order 1. We have shown how to accommodate a low Q in our formalism, by retaining exponential factors that could have been set to 1 in the present high-Q calculations. We have not presented the full equations for low Q, but those follow after replacing (19) by the well known formula for a broad band resonator [39], and proceeding with nearly the same steps as before. Our iterative determination of the diagonal term v d of Section VI might have to be revised.
One could also include higher order modes of cavities, and whispering gallery modes describing coherent synchrotron radiation [40].
Besides improving the physical model of the equilibrium state, an urgent matter is to study the stability of the equilibrium. This can of course be done by macro-particle simulations, but we would like to appeal as much as possible to direct solution of the Vlasov-Fokker-Planck equation by the method of local characteristics, which proceeds with very low numerical noise [41]. This can be done easily for the case of a complete fill, with only one phase space distribution to contend with. The present study also suggests possible reduced models of trains with gaps, in which identity of some of the bunches would be enforced in one way or another. Our technique of exploiting geometric sums can help to simplify expressions for the induced voltage.
A special point of interest is the effect of over-stretching on thresholds of instability.
We have seen, without accounting for stability, that over-stretching can give an additional increase in the Touschek lifetime.

XVI. ACKNOWLEDGMENTS
We thank several colleagues for advice and references to the literature: Karl Bane, To derive the second equation of motion, note that the azimuthal location of a particle with revolution frequency ω 0 + ∆ω is If the cavity is at θ = 0, the n-th passage of the cavity occurs at time t n such that 2πn = (ω 0 + ∆ω)t n .
At that time the cavity phase is The term 2πnh on the right hand side can be dropped, since it does not affect the applied voltage V 1 sin φ n , nor the induced voltage V r ((φ n −φ 0 )/k 1 −m j λ 1 ). Indeed, under substitution of that term, the argument φ n /k 1 = φ n λ 1 /2π takes on the value nhλ 1 = nC. Since V r is periodic with period C, it is not changed by the presence of the term 2πnh in φ n .
With the definition (8) we then have For highly relativistic particles above transition, the momentum compaction factor α can be written as Passing to the corresponding differential equation, we obtain (15).

Appendix B: Perturbed synchronous phase
To determine the perturbed synchronous phase we put k 1 z i + φ 0 = φ 0i in (16). We then see that the synchronous phase φ 0i for the i-th bunch is defined by It is the phase at which the force is zero, at the center of the distorted potential well.
In (B1) we have a nonlinear equation to solve for φ 0i . If the equation (B1) can be solved, one can work out the dynamics for a new variablez i defined by where φ i is the dynamical phase of the applied voltage when the particle arrives. That is, the applied voltage takes the form V 1 sin(k 1zi + φ 0i ), andz i is zero at the minimum of the distorted potential well.
The scheme now involves a two-part iteration. In an iterate of Part 1 the synchronous phases φ 0i are determined by solving (B1) with a given function V r . In a succeeding iterate of Part 2, those φ 0i are used to calculate the charge densities and thus to form a new value of V r , by the algorithm described in Section X.
We programmed this scheme and found that it converges at moderate current, but runs into difficulty near the design current, because at that current we are getting close to the situation in which (B1) does not have a unique solution, owing to the advent of a doubly peaked charge density. At moderate current the results agree quite precisely with those from the simpler scheme based on the current-independent φ 0 and the original variables z i .
Since the simpler scheme works at any current up to the design current and even far beyond, we have applied it for all further work. It is not necessary to base the coordinate system on the synchronous phases, but they can be found a posteriori as the location of the minima of the distorted potential wells computed using the z i as coordinates.
Appendix C: Reduction to the case with all buckets filled Here our task is to reduce our general formula for the induced voltage to its form when all r.f. buckets are filled. We shall find that the resulting expression agrees to a close approximation with a formula well known in the literature. Thus the following slightly complicated calculation serves as a good check on the preceding work.
In the equilibrium state the bunches will all have the same charge density ρ(z). We adapt the methods of Sections IV and V. The total charge density will be ρ tot (z) = ∞ p=−∞ h j=1 ρ(z + (j − 1)λ 1 + pC) .
By translating the integration variable we get the induced voltage as In the notation of (21) the sum over p in this expression is which can be calculated from the formula (30). To apply the formula we first show that to an excellent approximation. Since 0 ≤ j − 1 ≤ h − 1 and |z|, |z | < Σ, we have From this we can evaluate the ceiling function that appears in the definition (31) of χ. At the lower and upper bounds of its argument from (C5) we have −1 + (λ 1 − 2Σ)/C = 0, since 2Σ λ 1 , 2Σ/C = 1 .
The evaluation (C7) occurs only for j = 1 and only for the part of the integration where z − z > 0. Since there are several hundred terms of similar magnitude in the sum over j, this case may safely be ignored. Thus only the evaluation (C6) occurs, which implies that (C4) is correct.
Note that r is the same as in (28).