Equilibrium of an Arbitrary Bunch Train in the Presence of Multiple Resonator Wake Fields

A higher harmonic cavity (HHC), used to cause bunch lengthening for an increase in the Touschek lifetime, is a feature of several fourth generation synchrotron light sources. The desired bunch lengthening is complicated by the presence of required gaps in the bunch train. In a recent paper the author and Venturini studied the effect of various fill patterns by calculating the charge densities in the equilibrium state, through coupled Ha\"issinski equations. We assumed that the only collective force was from the beam loading (wake field) of the harmonic cavity in its lowest mode. The present paper improves the notation and organization of the equations so as to allow an easy inclusion of multiple resonator wake fields. This allows one to study the effects of beam loading of the main accelerating cavity, higher order modes of the cavities, and short range geometric wakes represented by low-$Q$ resonators. As an example these effects are explored for ALS-U. The compensation of the induced voltage in the main cavity, achieved in practice by a feedback system, is modeled by adjustment of the generator voltage through a new iterative scheme. Except in the case of a complete fill, the compensated main cavity beam loading has a substantial effect on the bunch profiles and the Touschek lifetimes. A $Q=6$ resonator, approximating the effect of a realistic short range wake, is also consequential for the bunch forms.


I. INTRODUCTION
This is a sequel to Ref. [1], in which we explored the action of a higher harmonic cavity (HHC), a standard component of 4th generation synchrotron light sources, employed to lengthen the bunch and reduce the effect of Touschek scattering. In that work we introduced an effective scheme to compute the equilibrium state of charge densities in an arbitrary bunch train. The train is allowed to have arbitrary gaps and bunch charges. We chose the simplest possible physical model, in which the only induced voltage (wake field) is due to the lowest mode of the HHC. We write V r3 for this voltage, the notation designating "resonator, 3rd harmonic". We recognized, however, that excitation of the main accelerating cavity (MC) by the bunch train produces an induced voltage V r1 of comparable magnitude, the effect described as beam loading. Our excuse for omitting V r1 was that in practice it is largely cancelled by adjusting the rf generator voltage V g through a feedback system. The sum of V r1 and V g should closely approximate V rf , the desired accelerating voltage.
In real machines there are always gaps in the bunch train, and that leads to varying bunch profiles and centroid displacements along the train. At first sight this would suggest that V r1 would be different for different bunches, so that compensation could only be partial, perhaps only manifest in some average sense. On the contrary, we shall calculate an equilibrium state in which the compensation is essentially perfect for all bunches. This happens by an adjustment of the charge densities of all bunches, to new forms that sometimes differ substantially from those without the MC. The adjustment is achieved automatically through a new algorithm presented here. This iterative procedure minimizes a mean square deviation of V r1 + V g from V rf , summed over all bunches, as a function of two generator parameters, which are equivalent to amplitude and phase.
It is not clear that this is a faithful model of the feedback mechanism, which could conceivably amount to a weaker constraint on the bunch profiles. Nevertheless, this study clarifies the mathematical structure of the problem, and appears to be a worthwhile preliminary to a full time-dependent model of the system including a realistic description of feedback.
Beside the main cavity, we also assess the role of the short range wake field from geometric aberrations in the vacuum chamber, and higher order modes in the HHC. These effects are added with the help of improvements in notation and organization of the equations.
Extensive numerical results are reported with parameters for ALS-U. For consistency with the previous work the parameters chosen are partially out of date as the machine design stands at present, but that will not greatly affect our pricipal conclusions. Although the qualitative picture of the model with HHC alone is still in place, there are large quantitative changes. Even then, we underestimate the full effects, because we can only get convergence of our iterative method when the current is a few percent less than the design current.
In Section II we briefly recall our previous algorithm for solving the coupled Haïssinski equations. Section III introduces the improved notation and organization which allows an easy inclusion of multiple resonator wake fields. Section IV enlarges the system of equations to provide a calculation of the diagonal terms in the potential, thus overcoming a limitation of the previous formulation. Section V describes the method for determining the generator parameters so as to compensate the induced voltage in the main cavity. Section VI, with several subsections, reports numerical results for the case of ALS-U [2,3], always making comparisons to results with only the HHC in place. Subsection VI A treats the case of a complete fill, illustrating the compensation of the main cavity in the simplest instance.
Subsection VI B considers a partial fill with distributed gaps, as proposed for the machine.
Subsection VI C is concerned with over-stretching by reduction of the HHC detuning. Subsection VI D explores the effect of the short range wake field with a realistic wake potential.
Subsection VI E checks the effect of the principal higher order mode of the HHC. Subsection VI F presents our closest approach to a realistic model, including the harmonic cavity, the compensated main cavity, and the short range wake, altogether. Subsection VI G examines the effect of the main cavity when there is only a single gap in the bunch train. Section VII reviews our conclusions and possibilities for further work. Appendix A derives the expression for the diagonal terms in the potential, for resonators of arbitrary Q.

II. SUMMARY OF METHOD TO COMPUTE THE EQUILIBRIUM CHARGE DENSITIES
In [1] we derived a set of equations to determine the equilibrium charge densities of n b bunches, which may be stated succinctly as follows: Here I is the average current andρ is a vector with 2n b real components, consisting of the real and imaginary parts ofρ i (k r3 ), where k r3 is the wave number of the lowest resonant mode of the 3rd harmonic cavity. These quantities are defined in terms of the beam frame charge densities ρ i (z), normalized to 1 on its region of support [−Σ, Σ], aŝ where Q is the quality factor of the cavity. The vector in (1) is arranged as follows: Accordingly, F in (1) is a real vector with 2n b components, so that we have 2n b nonlinear algebraic equations in 2n b unknowns, depending on the parameter I.
For the high Q of a typical HHC the quantity (2) is very close to the Fourier transform, but we have persistently written all equations for general Q for later applications involving low-Q resonators.
In (1) the diagonal terms of the induced voltage have been dropped, i.e. the effects on a bunch of its own excitation of the cavity. This omission is justified for the typical high Q of an HHC. Our method to handle the diagonal terms in the general case is introduced in Section IV.
A solutionρ of (1) determines the charge densities by the formula of Eq.(50) in [1], where U i is the potential felt by the i-th bunch, defined in Eq.(51) of [1]. Here µ and A i are constants, and z i is the beam frame longitudinal coordinate of the i-th bunch. The potential U i depends on all components ofρ, on the mean energy loss per turn U 0 , and on the parameters of the applied voltage V rf which we write as We solve (1) by the matrix version of Newton's iteration, defined in (67) of [1]. We begin at small current I, taking all components of the first guess forρ to be the transform (2) of a Gaussian with the natural bunch length. We then continue step-wise to the desired current, making a linear extrapolation in current to provide a starting guess for the next Newton iteration at incremented current. The extrapolation is accomplished by solving for ∂ρ/∂I from the I-derivative of (1):

III. FORMALISM FOR MULTIPLE RESONATORS
The scheme allows the inclusion of any number of resonator wake fields, but to do that conveniently requires some care in notation and organization of the equations. With n r resonators there are 2n b n r = n u unknowns, which we assemble in one long vectorρ : Reρ 1 (k r,nr ), · · · , Reρ n b (k r,nr ) , Imρ 1 (k r,nr ), · · · , Imρ n b (k r,nr ) .
Here k r,n is the resonant wave number of the n-th resonator, and the subscript ofρ denotes as usual the bunch number.
To identify the bunch number and the resonator number for the k-th component of the vector, we define two index maps: ι(k) which gives the bunch number and r(k) which gives the resonator number. Namely, Here x , the ceiling of x, is the least integer greater than or equal to x. We also need two projection operators: P re (k) which is equal to 1 if k corresponds to a Reρ and is zero otherwise, and P im (k) which is equal to 1 if k corresponds to a Imρ and is zero otherwise.
These are expressed in terms of the ceiling of k/n b as follows: The potential U j (z) for bunch j, generalizing Eq.(51) of [1] to allow n r resonators, is stated as The first term in (11) is e times the integral of the applied voltage, now called the generator voltage and written as At x 1 = cos φ 1 , x 2 = sin φ 1 this reduces to the desired V rf of (5). In an amplitude-phase representation we have The first term in (12) represents the diagonal contributions, the effect on bunch j of its own excitation of the resonators, as opposed to excitation by the other bunches which is described by the second term. By writing the latter as a simple matrix-vector product we greatly simplify the calculation of the Jacobian of the system, making it formally the same for any number of resonators.
Referring to Eqs.(28), (51), (55), (56), (57), (58) of [1], we can write down the matrix elements M (z) i,k in the second term of (12). For this we introduce a notation appropriate for labeling by the index k of (7). Functions of k, defined via the index maps, are labeled with a tilde:k where The result for the matrix from (51) and (57) of [1] is seen to be (noting that ω r /k r = c) In the present notation the system of coupled Haïssinski equations, generalizing (66) of [1], takes the form The normalization integral appearing in the first term is We require the Jacobian matrix [∂F j /∂ρ k ] for the solution of (18) by Newton's method, assuming that the diagonal terms are fixed. This is found immediately from (12), (18), and The compact expressions in (17) For the work of the following section we also need the induced voltage from the main cavity, which we designate as the first resonator in the list (n = 1). For the j-th bunch this takes the form The diagonal term v d 1j can be evaluated in terms of integrals derived in Appendix A.

IV. THE FULL SYSTEM OF EQUATIONS WITH DIAGONAL TERMS
Through (18) we have a system of n u algebraic equations for determination ofρ, provided that the diagonal terms in U i are given. The latter are functionals of the charge densities ρ i (z i ), from which it follows that (18) can be stated in vector notation as On the other hand, the ρ i (z i ) are determined in turn as solutions of integral equations provided thatρ is given. The integral equations are like normal single-bunch Haïssinski equations, but with a background potential determined byρ, namely In vector notation The potential U i depends on the ρ i through its diagonal terms, in the first sum in (12). Our procedure will be to interleave the solution of (22) at fixedρ, by the usual Newton method, with the solution of (24) at fixedρ . If this algorithm converges we shall have consistency between ρ andρ and a solution of the full system.
It turns out, most fortunately, that the solution of (24) is obtained by plain iteration as would be applied to a contraction mapping, In our application this usually converges to adequate accuracy in just one step, or three at most, and takes negligible time.
This scheme based on (22) and (24) is used in all calculations reported below. It replaces the method used in [1], which was to evaluate the diagonal terms from the value of ρ from the previous Newton iterate. That works only for high-Q resonators, so is not adequate for handling the short range machine wake.
We wish to choose (x 1 , x 2 ) so as to minimize, in some sense, the difference for all i = 1, · · · , n b . A reasonable and convenient choice for an objective function to minimize is the sum of the squared L 2 norms of the quantities (26). With a normalizing factor to make it dimensionless and of convenient magnitude that is The region of integration [−Σ, Σ] is the same as that used in the definition of the potential Note that the minimum of f cannot be strictly zero, since V r1 is sinusoidal with wave number k r1 , whereas the other terms are sinusoidal with a slightly different wave number Let us adopt the vector notation x = (x 1 , x 2 ) with norm |x| = |x 1 | + |x 2 |. The equations to solve now depend on x, having the form To avoid notational clutter we suppress reference to the diagonal terms, leaving it understood that a solution of (28) forρ actually involves the scheme of the previous session. As usual we solve forρ, for an increasing sequence of I-values. The scheme will be to minimize f (x) at each I, thus providing a new x = arg min f to be used at the next value of I. As will now be explained, the minimization will also be done iteratively, so that we have an x-iteration embedded in theρ-iteration.
We wish to zero ∇ x F , which is to find x to solve the equations To solve (29) a first thought might be to apply Newton's method, starting at some low where with By (30) we have an update x 0 → x which establishes the pattern of the general iterate x (k) → x (k+1) . This will be carried to convergence in the sense |x (k+1) − x (k) | < x , with a suitable x to be determined by experiment. Each iterate requires a value for V r1i and for ∇ x V r1i , which we compute numerically by a divided difference, Thus one x-iteration requires threeρ-iterations to provide the necessary values of V r1i (which are constructed fromρ). The firstρ iteration to find V r1i (z, x 1 , x 2 ) produces aρ which is a very good guess to start the remaining two iterations to make the derivatives, which then converge quickly.
The choice of ∆x in (34) requires a compromise between accuracy and avoiding round-off error. We found that ∆x = 10 −4 was widely satisfactory, whereas success with smaller values depended on the circumstances.

VI. NUMERICAL RESULTS WITH AND WITHOUT THE MAIN CAVITY
As in [1] we illustrate with parameters for ALS-U [2,3], the forthcoming Advanced Light Source Upgrade. Although the machine design is not yet final, one provisional set of parameters for our main cavity (actually the effect of two cavities together) is as follows: Here the shunt impedance R s and quality factor Q are loaded values, the unloaded values divided by 1 + β, with coupling parameter β = 7.233. We take these parameters for the main cavity, otherwise keeping the same parameters as in [1], Table I. Thus we take U 0 =217 keV, even though a value of 330 keV may be contemplated for the set (35).

A. Complete Fill
We first take the case of a complete fill, thus n b = h = 328. The average current is to be 500 mA, which we reach in 8 steps starting from 200 mA. The CPU time is 15 minutes, rather than 20 seconds for the calculation without the main cavity. The increase is mostly due to a much slower convergence of theρ-iteration, the x-iteration being a minor factor in CPU time. To save time we gave x the rather large value of 0.05, but then made a refinement to x = 10 −6 at the final current, in an extra 2 minute. The steepness of the objective function f (x1, x2) of (27) is extraordinary, having values around 10 4 in the sequence with x = 0.05 while falling to a value close to 1 after the refinement. An interesting question is how this steepness would be reflected in a feedback system.
The result for the charge density, shown in the blue curve of Fig.1, is quite close to the result without the main cavity, shown in red. It should be emphasized that there is no explicit constraint requiring all bunches to be the same. We have computed 328 bunches separately, and have found that they all come out to be the same. This constitutes a good check on the correctness of the equations and the code. In Fig.2 we show the compensation mechanism. The sum of the generator voltage V g and the induced voltage V r1 from the main cavity is the orange curve. The latter deviates from the desired effective voltage V rf by less than 2%, as is seen Fig.3.
The phasor of the generator voltage moves closer to π/2 and its magnitude

B. Partial fill C2 with distributed gaps
Next we take a partial fill with distributed gaps, labeled as fill C2; see Section XIII-C of [1]. There are 284 bunches in 11 trains, with 4 empty buckets between trains. There are 9 trains of 26 and 2 of 25, with the latter positioned at opposite sides of the ring. All bunches have the same charge. As in the preceding example we start the calculation at low average current and advance in steps trying to reach the desired 500 mA. The convergence of iterations is at first similar to that of the preceding case, but begins to falter around 430 mA average current, at which point the convergence of theρ-iteration becomes problematic.
By taking smaller and smaller steps in current we can reach 496 mA, but beyond that point the Jacobian matrix of the system appears to approach a singularity, as is indicated by its estimated condition number having a precipitous increase, from 700 at the last good solution to 2900 at a slightly higher current. Nevertheless, the x-iterations continue to converge as long as theρ-iterations do. In the following, graphs are plotted for the maximum achievable current, stated in figure captions. Now the plots of V g and V r1 and their sum look exactly the same as in Fig.2, for every bunch. The minimization of f (x 1 , x 2 ) has caused the bunch forms to rearrange themselves so that the compensation is essentially perfect for every bunch. The deviation of V r1 + V g from V rf , scarcely visible on the scale of Fig.2, varies from bunch to bunch, but is still less than 3% for all bunches. The corresponding results for the bunch centroids is seen in Figs.7. Again the deviation from the case without the main cavity is quite substantial.
The main point of practical interest is the increase in Touschek lifetime achieved through the bunch stretching caused by the HHC. Again, the MC has a sizeable effect in reducing the lifetime and in causing a larger variation along a train. This is shown in Fig.8 which gives the ratio of the lifetime τ to the lifetime τ 0 without the MC.
We next consider the same fill pattern with 11 trains, but with a taper in the bunch charges putting more charge at the ends, according to a power law as shown in Fig. 15 of [1]. This is an example of invoking guard bunches to reduce the effect of gaps. As is seen in Figs.9 and 10, the guarded inner bunches, which resemble that of the complete fill, are little affected by the MC. The strong asymmetry between the front and back of the train is perhaps surprising, but it should be noticed that Fig.10 already shows an appreciable front-back asymmetry. The strong amplification of this asymmetry by the MC is in line with its big effects seen generally.  peak. In our case a decrease from df = 250.2 kHz to df = 235 kHz produces a double peak in the model without the MC at full current, as is seen in Fig.4 of [1]. We would like to know how this setup looks with the compensated main cavity in play. Not surprisingly, the convergence of our iterative solution breaks down at a lower current than in the case of the normal detuning; the stronger the bunch distortions the poorer the convergence. With  Even at a current significantly less that the 500 mA design current the distortion due to the main cavity is quite large, which leads to the conclusion that the main cavity must be included in a realistic simulation of over-stretching.

D. Effect of the short range wake field
The short range wake field from various unavoidable corrugations in the vacuum chamber retains importance in the latest storage rings, in spite of the best efforts to reduce it. Since it can cause substantial bunch distortion in the absence of an HHC, we would like to know how much it affects the operation of the HHC. A result for the longitudinal wake potential at ALS-U, from a detailed computation by Dan Wang [4], is shown in Fig.13. The corresponding impedance, is plotted in Fig.14.
In Ref. [1] we suggested that a low-Q resonator wake could be treated on the same footing as the high-Q resonators, and for that reason we wrote all equations for a general value of Q.
We recognized, however, that the diagonal term in the potential would now be dominant, while being nearly negligible in the high-Q case. It could not be treated by the method used in [1], but is easily handled by the presently adopted method of Section IV.
For the equilibrium state, the impedance at f > 20 GHz is irrelevant, even though it could have a role out of equilibrium. This assertion follows from the fact that the frequency spectrum of our calculated charge densities never extends beyond 15 GHZ, no matter which wake fields are included. Consequently, a reasonable step is to concentrate on the first big peak at 11.5 GHz. The wake potential in our equations (defined in (19) of [1]) is based on an impedance as follows, which is of Lorentzian form with half-width Γ/2: Figures 15 and 16 show a fit to (38) with parameters as follows: f r = 11.549 GHz , R s = 5730 Ω , Q = 6 .
The fit is rough in the imaginary part, but probably good enough to estimate the magnitude of the effect of the short range wake.
Discussions of low-Q resonator models in the literature usually invoke the impedance of , often with Q near 1. As is illustrated in Figures 15 and 16, in our case with Q = 6 the LRC model does not give a better fit than the simpler Lorentzian, except for enforcing Z(0) = 0. At the expense of some complication our equations could be modified to accommodate the LRC form, but that appears to be unnecessary, at least in the present example.
Henceforth, the impedance from (38) and (39) will be referred to as SR (short range).
Taking first a complete fill, and including just the HHC and SR, we get the result of Fig.17. Next we consider the partial fill C2 with distributed gaps as treated in the previous section. Figures 18 and 19 show the results for HHC+SR and HHC alone. As expected, the effects of SR are more pronounced in the partial fill than in the complete fill. Correspondingly, the maximum current achieved is 472.6 mA. As in previous cases we expect a substantially larger effect at the design current of 500 mA. At least for the equilibrium state in ALS-U, it appears that the HOM can be neglected.
The role of HOM's in longitudinal coupled-bunch instabilities is discussed in Ref. [6].
F. The full model: HHC+MC+SR.
We are now prepared to include the harmonic cavity, the compensated main cavity, and the short range wake, altogether. The convergence of the Newton sequence suffers even more than in the previous cases, and the continuation in current reaches only I av = 471.9 mA.
Effects seen at this current must severely underestimate what can be expected at 500 mA, because of the strong variation near the design current that we have observed in every case.
For fill C2 we see the charge densities in Figures 21 and 22. bunch lengthenings is shown in Fig.25.

VII. CONCLUSIONS AND OUTLOOK
Continuing the investigation of Ref. [1] we have extended the physical model to include the effect of the main accelerating cavity in its fundamental mode, previously omitted. We introduced a new algorithm to adjust the parameters of the rf generator voltage so as to compensate the voltage induced in the cavity by the beam, thus putting the net accelerating voltage at a desired value. When the cavity is excited by a bunch train with gaps this compensation implies a modification of the bunch profiles, which is produced automatically in our scheme.
We illustrated the outcome for parameters of the forthcoming ALS-U storage ring, revisiting examples treated in [1] without the main cavity. The results are similar in modo grosso, but there are significant quantitative differences, especially in cases of overstretching of bunches. Generally speaking there is more bunch distortion and less symmetrical patterns in the bunch trains, and the rms bunch lengthening is a bit smaller and much more variable along the train. Correspondingly, the Touschek lifetime increase is smaller and more variable over a train.
We have not tried to model the feedback system that compensates the beam loading in practice. Our aim was only to show the theoretical existence of an equilibrium state with precise compensation in place.
It was disappointing, and somewhat surprising, to find that the Newton iteration to solve the coupled Haïsinski equations encounters convergence difficulties at large current (near the design current) when either the main cavity wake or the short range wake is added to the HHC wake.
A colleague suggested that the failure of convergence might hint at an instability. One should be cautious about such an idea. The issue here is just the existence of an equilibrium.
An equilibrium may or may not be stable under time evolution, so stability is a different issue.
Our failure to find an equilibrium in some cases of high current may be due to a failure of technique, not necessarily an indication that no equilibrium exists. At high current we are trying to achieve convergence of the Newton iteration close to a singularity of the Jacobian, but not squarely on the singularity. In this case it is crucial to have a starting guess sufficiently close to a solution, but in practice the required degree of closeness is unknown.
We made some efforts to improve the guess by a seemingly careful continuation in current from the last good solution, but there was no clear success.
A likely remedy for the convergence failure is to return to the conventional formulation of the Haïssinski equations as integral equations for the charge densities, in place of the present formulation as algebraic equations for Fourier amplitudes. For a single Haïssinski integral equation discretized on a mesh in z-space, the Newton iterative solution is ultra-robust, converging at currents far beyond realistic values [7]. It seems likely that similar good behavior will hold for the coupled integral equations. The size of the discretized system does not grow with the number of resonator wakes, in contrast to the present system, and the full z-space description of the short range wake could be invoked in place of the low-Q resonator model.
To make this z-space formulation feasible on modest computer resources we can assume that all bunch sub-trains are identical, and all separated by identical gaps. For the ALS-U this would mean artificially increasing the harmonic number from 328 to 330, and having 11 trains of 26 separated by gaps of 4 buckets. Then we have 26 independent charge densities, which can adequately be described by 100 mesh points each. Thus the Jacobian of the Newton iteration is 2600 × 2600, a modest size that will yield a very quick computation.
Moreover, this identical train model would make it feasible to do a time-domain solution of the coupled Vlasov-Fokker-Planck equations by the method of local characteristics (discretizing the Perron-Frobenius operator) [8]. This could answer the urgent question of stability of the equilibria, and provide a window to the dynamics out of equilibrium. Incidentally, determination of the instability threshold by the linearized Vlasov system could also be attempted.
It was gratifying to find that the sub-iteration to enforce the main cavity compensation converged very quickly whenever the main iteration converged. It can be employed in the same way in the proposed z-space system. Also, the formalism for multiple resonator wakes will still be advantageous in the z-space scheme.