Numerical solution of the Haïssinski equation for the equilibrium state of a stored electron beam

The longitudinal charge density of an electron beam in its equilibrium state is given by the solution of the Haïssinski equation, which provides a stationary solution of the Vlasov-Fokker-Planck equation. The physical input is the longitudinal wake potential. We formulate the Haïssinski equation as a nonlinear integral equation with the normalization integral stated as a functional of the solution. This equation can be solved in a simple way by the matrix version of Newtons’s iteration, beginning with the Gaussian as a first guess. We illustrate for several quasirealistic wake potentials. Convergence is extremely robust, even at currents much higher than nominal for the storage rings considered. The method overcomes limitations of earlier procedures, and provides the convenience of automatic normalization of the solution.


I. INTRODUCTION
The collective longitudinal motion of electrons or positrons in a storage ring seems to be well described by the Vlasov-Fokker-Planck (VFP) equation, in which the collective force is described by a wake potential which accounts for the electromagnetic environment due to the vacuum chamber.The equation has solutions that are stationary in time, which may or may not be stable under perturbations, depending on the value of the beam current.These are solutions of the Haïssinski equation [1], which may be stated as a nonlinear integral equation or integrodifferential equation.To determine the threshold in current for an instability to appear, one can linearize the Vlasov equation about the Haïssinski solution.Alternatively, one can integrate the VFP equation as an initial value problem in time, with the Haïssinski equilibrium as the initial value [2].In either approach, computation of the Haïssinski solution is an essential first step in determining the threshold for an instability.(Admittedly, one can also get an idea of stability by running a VFP integration from an arbitrary initial value, say a Gaussian, but that will lead to a somewhat ambiguous definition of the threshold.) The method of solution presented here was worked out by the first-named author twenty years ago, but was not published except for a description in words in Ref. [2].Although the method was adopted by a few colleagues, it has not become a standard tool.Since it is quite simple and avoids limitations of other methods, a belated publication seems worthwhile.
The idea of the method will seem obvious to anyone acquainted with ideas of functional analysis [3] and their application in numerical methods [4].The integral equation for the charge density λ is viewed as an equation FðλÞ ¼ 0 on an appropriate function space.The equation is discretized by a numerical quadrature rule for the integrals involved, and then solved by the matrix version of Newton's method.An essential step is to define F so that a solution is automatically normalized.We shall not be concerned with a rigorous basis for discretization, but methods to treat that issue are available [5].
What was not so obvious before implementation is the extremely robust convergence of the Newton iterates.For realistic wake potentials we have never seen a failure of convergence to machine precision in a few iterations, even at currents far beyond the threshold of instability.Here the starting point for the Newton iteration was merely a Gaussian, the zero current solution.
In this paper we include results for high current, in order to explore the mathematical properties of the equation, and to demonstrate the excellent convergence of the iteration.In many physical problems governed by nonlinear equations one finds critical values of some strength parameter, for instance in the buckling of a column at a certain value of the load [6], ( [3], Chap.4).As a function of the parameter a solution may branch into two or more solutions, or become complex, or simply cease to exist in one way or another.It is then natural to look for critical points in the current parameter in the Haïssinski equation.For the wake potentials considered here we find no such points up to very high currents.In fact, we can argue that our solutions are locally unique in the function space considered, since bifurcation can occur only at a singularity of the Jacobian of the system [3].Our Jacobian is always far from singularity.The high current solutions represent unstable equilibria, and will not be realized in the laboratory.

II. SOLUTION OF THE VLASOV-FOKKER-PLANCK EQUATION FOR THE EQUILIBRIUM STATE
We are concerned with longitudinal motion within a single bunch of particles in an electron storage ring.The linearized motion without collective effects is described in terms of the slip factor η, the dimensionless constant which relates the first order change in revolution frequency ω r to a change in momentum P: Here ω 0 and P 0 are the nominal values of revolution frequency and momentum (those for a particle synchronizing with the rf), γ 0 is the nominal Lorentz factor, and α is the momentum compaction factor.(Some authors define η with the opposite sign, and some call η the momentum compaction factor.) The dynamical variables of longitudinal motion are often taken to be Δϕ and ΔE, where Δϕ is the deviation of the rf phase from its synchronous value at the time the particle encounters the rf field, and ΔE is the deviation of the energy from the nominal value E 0 [7].We prefer to work with equivalent dimensionless variables q and p, normalized to be of order 1, and a corresponding dimensionless τ, equivalent to the time [8].We define Here z ¼ s − s 0 ¼ s − β 0 ct is the distance (in arc length s on the reference orbit) to the synchronous particle at s ¼ s 0 , thus positive for a leading particle, and sgnðηÞ is the signum function, equal to 1 for η > 0 and −1 for η < 0. The constant ω s ¼ 2πf s is the circular synchrotron frequency.At first we think of σ z and σ E as some positive constants to render q and p dimensionless and of order 1, leaving to later a specific choice of their values.One can show that Δϕ ¼ −hz=R, where h is the harmonic number and R ¼ C=2π, where C is the circumference of the reference orbit followed by the synchronous particle.From this we can write the differential equations [7] (which approximate a discrete map) in terms of the new variables as follows, The corresponding Hamiltonian is In a storage ring with normal equilibration from synchrotron radiation, in which effects of diffusion balance effects of dissipation, the phase space density in the limit of small beam current is where A is a constant for normalization.Hence we can interpret σ z and σ E as the r.m.s.bunch length and energy spread for weak current, provided that a ¼ 1 or Henceforth we choose σ z and σ E to satisfy (6), whatever their interpretation.The probability density in phase space, normalized to 1, is denoted by fðq; p; τÞ, and the spatial probability density by λðq; τÞ, thus The Vlasov-Fokker-Planck (VFP) equation to determine where t d is the longitudinal damping time.
The single-particle equations of motion, modified to include the Vlasov collective force, are where −q is the linear force from rf and −F is the collective force.Here F is a functional of the distribution f which is assumed to have the form The wake potential W is defined to have the dimension of a potential per unit charge, and to be positive where it causes an energy gain.It follows that the normalized current I has the form where N is the number of particles and ν s ¼ ω s =ω r is the synchrotron tune.The notation Fðq; fð•; τÞÞ is intended to indicate that F depends on all values of f over phase space at time τ.In some models the wake potential is zero in front of the bunch (q > 0) but that is not assumed in the following.
Although the formula (10) usually goes unquestioned, it is in fact not the most general form of the collective force for a time dependent charge density.For a bunch on a curved orbit it does not account for the charge density being different at the retarded time from what it is at the current time [9,10].
In view of ( 9) the VFP equation takes the form where β ¼ ðω s t d Þ −1 , with t d being the longitudinal damping time.The Fokker-Planck terms on the right-hand side account for damping and diffusion due to incoherent synchrotron radiation.We are interested in an equilibrium, a time-independent solution of ( 12), denoted by f 0 ðq; pÞ.We seek such a solution in the Maxwell-Boltzmann form The Fokker-Planck terms add up to zero for this Gaussian function of p, owing to compensation of diffusion by damping.Thus f 0 will be an equilibrium solution provided that the spatial density λ satisfies Any solution of ( 14) may be represented as follows: where the constant A is chosen to enforce the normalization of ( 13).This follows from separation of variables [dλ=λ ¼ −ðq þ FÞdq] and integration.Now it is convenient to reverse the order of integrations, after introducing the integrated wake potential S, where The kernels Wðq − q 0 Þ and Sðq − q 0 Þ may be viewed as giving the response to a delta function source and a step function source, respectively [11].
It follows from ( 15) and ( 17) that a normalized solution of ( 14) must satisfy This nonlinear integral equation ( 18) is our main object of study.It is convenient to rewrite it as where φðqÞ ¼ IλðqÞ and Fðφ;IÞ ¼ φðqÞ with all integrations on ð−∞; ∞Þ.

III. PREVIOUS METHODS OF SOLVING THE HAÏSSINSKI EQUATION
To motivate our method we briefly review techniques in common use, and point out limitations that are avoided by our algorithm.

A. Solution by simple iteration
An obvious approach is to generate a sequence of functions fλ ð0Þ ; λ ð1Þ ; Á Á Ág by the rule where λ ð0Þ is the normalized Gaussian and A has some trial value, say 1= ffiffiffiffiffi ffi 2π p .If the sequence converges, try again with different values of A, searching for a value of A such that the final iterate is normalized to adequate precision.This could be made more convenient by normalizing every iterate; in other words, just apply simple iteration to our Eq.( 18) with embedded normalization, so that Unfortunately, in numerical experience this sequence or (21) fails to converge at larger I, including values of practical interest.Rather, the iterates eventually oscillate between one pattern and another.This failure has no physical significance, as is shown by successful continuation of the solution to large I by other methods, for instance the one we advocate.

B. Solution of the equation in integro-differential form
This method aims to solve the Haïssinski equation expressed as the integro-differential equation of ( 14).This can be done in a simple way only if WðqÞ ¼ 0 for q > 0, a condition that is not strictly true for numerically determined wake potentials for real storage rings.In fact such potentials are nonzero in a small region 0 < q < a.A more serious violation of the condition can occur in the case of coherent synchrotron radiation.Depending on circumstances, it may happen that WðqÞ will be non-zero over a large range of positive q We seek a numerical solution which is strictly Gaussian for q ≥ κ, approximating the actual solution which is asymptotic to a Gaussian.We write λðqÞ ¼ A expð−q 2 =2Þ; q ≥ κ.Then with the above mentioned restriction on W the integro-differential equation to solve is The idea is to start at q ¼ κ, where the right-hand side is known, then integrate backwards in 2N steps of −Δq ¼ −κ=N to q ¼ −κ.If we apply Euler's method, the first two integration steps are as follows: The integral in the last term in (25) can be approximated by the trapezoidal rule as Thus λðκ − ΔqÞ and λðκ − 2ΔqÞ are determined by ( 24), ( 25) and (26).Continuing in a similar way we build up the discretized solution λðκ − iΔqÞ; i ¼ 0; …; 2N, which depends on the constant A in the initial condition.The process must be repeated to search for an A such that the solution is normalized.The solution, unnormalized in general, is well defined for any current I, so if normalization can be achieved we have overcome the restriction to small current required in the iterative method.Unfortunately we see no simple way to automate the normalization.The awkwardness in normalization, and the requirement that WðqÞ vanish for q > 0, are two undesirable features of this method that we wish to avoid.

C. Solution by time-domain integration of the Vlasov-Fokker-Planck equation
Another possibility is to integrate the full VFP equation (12) as an initial-value problem, using the method of local characteristics [2].With fðq;p;0Þ¼expðq 2 þp 2 Þ=2π as the initial value, the solution is expected to converge to the Haïssinski solution at large τ, provided that the current is below the threshold for instability.The disadvantage of this approach is that is does not allow the study of currents above threshold, and it takes much more computer time.It does provide, however, a useful check of the VFP solution algorithm, given a Haïssinski solution from another method.

IV. NUMERICAL SOLUTION OF THE NONLINEAR INTEGRAL EQUATION BY NEWTON'S METHOD
We discretize the equation ( 19) on a uniform mesh of n points q i , running from −κ to κ: We write φ i for the numerical approximation to φðq i Þ, and S i−j for Sðq i − q j Þ.We discretize the integrals by some quadrature rule with weights w i .Then the discretized form of ( 19) is Newton's method defines a sequence of approximations by successive linearizations of the equation.If φ ðpÞ is the pth approximate solution, then φ ðpþ1Þ is obtained from the first order Taylor development about φ ðpÞ : An initial guess φ ð0Þ , sufficiently close to the desired solution, is required.The Jacobian matrix element computed from (28) is Given φ ðpÞ , we compute Fðφ ðpÞ ; IÞ and ∂Fðφ ðpÞ ; IÞ=∂φ from (28) and (30) and then solve the system (29) of n linear equations for x ¼ φ ðpþ1Þ − φ ðpÞ to find the update φ ðpþ1Þ ¼ x þ φ ðpÞ .A convenient criterion for convergence may be stated in terms of a vector norm, for instance We judge convergence by the quantity demanding that it reach a small value, say 10 −14 , as p increases.Of course, one must also check convergence under refinement of the mesh (27).We normally do that just by graphical comparisons, but it could be done more quantitatively.
At sufficiently small current the Gaussian should be a suitable first guess, φ ð0Þ ¼ I expð−q 2 =2Þ= ffiffiffiffiffi ffi 2π p .In practice this choice is good for realistic currents with reasonable wake potentials, in fact at currents considerably higher than realistic.On the other hand, to understand the mathematical properties of the equation it may be useful to go to much higher currents.
An obvious approach to high current is to begin with the Gaussian and increase I in steps, taking a solution at I as the first guess for an attempted solution at I þ ΔI.An improvement to this idea can be achieved at little cost by instead using a linear extrapolation in I: The derivative is found by differentiating the I-dependent equation with respect to I: At the end of the Newton iteration for current I we have in hand both the Jacobian ∂FðIÞ=∂φ and the quantity ∂FðIÞ=∂I (from the second term of (28).Thus it takes only one solution of the linear system (34) to produce the required derivative dφðIÞ=dI for (33).
A convenient way to arrange the code is to make this method of advancing I always available, and so that the case of a single I is merely a special case.Thus one specifies the initial and final values of I, and the number of intermediate values, taken to be evenly spaced.This is convenient for plotting I-dependent quantities such as the centroid position or the r.m.s.bunch length of the Haïssinski distribution, and also for exploring the high current regime.

V. TESTS OF THE METHOD FOR QUASIREALISTIC WAKE POTENTIALS
We consider examples of the wake potential, obtained by solving Maxwell's equations with a quasirealistic model of the vacuum chamber providing the boundary conditions on metallic walls.The ideal wake potential W 0 ðqÞ, often called the delta wake, would be the longitudinal field Eðq; sÞ at a fixed normalized distance q from a point charge circulating on the ideal orbit, averaged in the position s over one turn.In practice the point charge is replaced by a short Gaussian charge distribution (a "driving bunch") to provide the approximated wake potential WðqÞ.This smooth function could be called a "pseudo-Green function" to distinguish it from a true Green function which cannot be smooth.If curvature of the orbit is neglected, W 0 ðqÞ displays "causality", in that it vanishes in front of the point charge (q > 0).In contrast WðqÞ will be nonzero for small q > 0, but will fall off quickly with increasing q.
Pioneering simulations of wake potentials were carried out for the damping rings of the Stanford Linear Collider (SLC).The model was axially symmetric, with the fields being computed by Weiland's code TBCI [12].The boundary conditions were for infinite conductivity of the chamber walls.There were two calculations, one for the original vacuum chamber [11], and one for a new vacuum chamber designed to have smoother walls [13].The latter replaced the original in an attempt to gain a higher threshold in current for a bunch instability.
Computer power and codes for electromagnetics have been greatly improved since the work for the SLC, and the physical model for newer storage rings has necessarily been extended to include coherent synchrotron radiation from curved orbits.A shorter driving bunch was needed, owing to shorter bunches in the new rings, and better electromagnetic codes allowed the inclusion of three-dimensional structures.An example of this more modern effort is a calculation for the low energy ring (LER) of KEKB (which is now out of service) [14].The model included CSR and resistive wall contributions as well as geometric wake fields.This ring allowed a configuration with negative momentum compaction [15], so we want to include that case in the Haïssinski solutions.
A rather different example of an ambitious calculation was for the positron ring of DAFNE at Frascati [16,17].There CSR was not important owing to the large bunch length, but the geometric structures were modeled very carefully.An overview of the SLC and DAFNE cases may be found in Ref. [18].
For each of the examples mentioned we make a cubic spline interpolation of the wake potential data from the relevant simulation, then integrate the spline analytically to make a smooth representation of the integrated potential SðqÞ for input to the integral equation (18).
For a qualitative comparison of the various cases it is useful to see how much the collective force −F resembles that from a linear combination of purely inductive and purely resistive components.The corresponding wake potentials are WðqÞ ¼ aδ 0 ðqÞ (inductive) and WðqÞ ¼ −bδðqÞ (resistive), where a and b are dimensionless positive constants.The corresponding force components are proportional to Z δ 0 ðq − q 0 Þλðq 0 Þdq 0 ¼ λ 0 ðqÞ; Signs are determined by the requirement that there be energy loss from particles at the front of the bunch.We make a weighted least-squares fit to the actual F by minimizing the following integral with respect to a and b: Z λðqÞ I Z Wðq − q 0 Þλðq 0 Þdq 0 − aλ 0 ðqÞ þ bλðqÞ Another way to make a qualitative comparison of cases is to plot the bunch centroid and the r.m.s.bunch length as a function of current.Such plots, along with the fit to the inductive plus resistive wake, will be given for each of our examples.

A. SLC damping ring with the original vacuum chamber
For details of this example see Ref. [11].For an rf voltage of 800 KeV the relevant parameters are as follows: For a typical bunch population of N ¼ 5 × 10 10 the normalized current is I ¼ 0.136 pC=V.
The wake potential WðqÞ computed with a Gaussian driving bunch with σ ¼ 0.5 mm is shown in Fig. 1

(left).
FIG. 1. Results for the SLC damping ring with its original vacuum chamber.Left: Wake potential WðqÞ; Right: Equilibrium charge density for N ¼ 5 × 10 10 (blue) and in the limit of zero current (red).
The finite extent of the driving bunch accounts for the potential not being zero at small positive q.
For solution of the integral equation we choose the weights w i for numerical quadrature to be those for Simpson's method; namely fw i g ¼ ðΔq=3Þð1; 4; 2; …; 2; 4; 1Þ..In the present and following examples we define the mesh in ( 27) with κ ¼ 6 for a mesh extending to 6σ z , and m ¼ 400 for 801 mesh points.We take r ¼ 10 −14 in (32) as the criterion for convergence.For a bunch population of N ¼ 5 × 10 10 , roughly the maximum that was stored in the ring, we get the Haïssinski solution shown in Fig. 1 (right).For comparison we show the Gaussian solution for zero current.The iteration to achieve this solution, beginning with the Gaussian, converged in 11 steps.
The good convergence is found to persist at much higher currents.In Fig. 2 (left) we compare solutions for N ¼ ð5; 10; 20; 30Þ × 10 10 .These solutions all started with the Gaussian as first guess, but the extent of the mesh had to be increased (taking κ ¼ 10) because of the increased bunch lengthening.At the highest current, 39 iterations were required.It is interesting to find that the value of r at the first iterate is always rather large, say 0.25, even at very small current.The same sequence of solutions is obtained by applying the method of continuation in I presented in Sec.IV.After a small step in I very few iterations are needed for convergence.
The very pronounced bunch lengthening in this example corresponds to a flat bottom in the distorted potential well.The well as given by ( 17) is shown in Fig. 2 (right), for the same sequence of currents.Since log λðqÞ ¼ −VðqÞ − log A, the wavy modulations in λ at high current must have a counterpart in VðqÞ.Taking the logarithm makes the modulations too small to be apparent on the scale of the graph of VðqÞ.
The fit to a sum of purely inductive and resistive wakes, plotted in Fig. 3 (right), shows that the inductive character is dominant within the bunch distribution: a=b ¼ 3.47.A purely inductive wake lengthens the bunch while keeping it symmetric about q ¼ 0, while a purely resistive wake makes the bunch lean forward with little change in its length.Accordingly, the bunch form overlayed in Fig. 3 shows relatively little leaning.
In Fig. 3 (left) we show the evolution with current of the normalized bunch length and centroid position.For N ¼ 5 × 10 10 the normalized current is I ¼ 0.136 pC=V.Both σ q and hqi show a steady increase, almost linear at high I.
The computations were done by a FORTRAN code using standard software for the linear algebra to solve for Newton iterates.The latter provides an estimate of the condition number of the Jacobian matrix, which turned out to be acceptably small, ranging from 10 at nominal current to 86 at six times nominal.Thus the iterates are numerically well defined.The computation time was negligible.In the following examples, the convergence and condition numbers were no worse than in the present case, and often better.

B. SLC damping ring with the improved vacuum chamber
This case is reviewed in Refs.[13,18].For an rf voltage of 800 KeV the relevant parameters are the following: For a typical bunch population of 5 × 10 10 we have I ¼ 0.14 pC=V.In Figs.4-6 we see a marked change in comparison to the case of the original damping ring.In Fig. 6 (right) the resistive component is shown to dominate the inductive: b=a ¼ 2.66.The bunch leans forward in the rf bucket, to compensate for the energy loss from the resistive wake field.The bunch length tends to saturate with increasing current, as is seen in Fig. 6 (left).The fall-off of charge density at the leading edge becomes sharper and sharper as the current increases, as is seen in Fig. 5.

C. KEKB low energy ring
This case is reviewed in Refs.[14,15].For an rf voltage of 8 MeV the relevant parameters are the following: For a typical bunch population of 6.6 × 10 10 we have I ¼ 0.0275 pC=V.We plot the charge density for this N in Fig. 7 (right) and the wake potential in Fig. 7 (left).The behavior of the charge density and the potential well with increasing current is shown in Fig. 8.This ring was able to run with negative momentum compaction.In Fig. 7 (right) we plot a charge density for that case in the black curve.This was obtained by changing the sign of I, while keeping all other parameters unchanged; see (11).The steep fall-off at the back of the bunch is typical for negative momentum compaction.
In this example the linear combination of inductive and resistive components provides a remarkably accurate representation of the collective force, as is seen in Fig. 9 (right).The inductive part dominates moderately, with a=b ¼ 2.18, so that we see the inductive pattern of strong bunch lengthening with relatively little forward tipping as the current is increased.Compare the behavior of bunch length versus current in Fig. 9 (left) with the corresponding graph Fig. 6 (left) for the previous example, in which the resistive part dominated.

D. DAFNE positron ring
See Refs.[16,17] for information on this case.For an rf voltage of 250 kV we have remarkably similar in our normalized variables.Again we have a very good fit to a sum of inductive and resistive components, with almost a 2∶1 ratio of inductive to resistive parts, as is seen in Fig. 12 (right).The pattern of bunch forms and bunch length versus current is very similar to that of KEKB.

VI. CONCLUSIONS AND OUTLOOK
We have demonstrated a simple and convenient method to solve the Haïssinski equation, with quasirealistic wake potentials for different kinds of electron storage rings.In work not covered in this report, we have also verified that the method works as well for the broad band resonator model of the wake potential, with similar or better experience regarding convergence.We have also applied the method with the wake from coherent synchrotron radiation, accounting for the "shielding" due to the vacuum chamber.The parallel plate model of the vacuum chamber was invoked in [19], and the toroidal model with resistive wall in [20].
There is scope for mathematical analysis of the Haïssinski equation, which we hope to present in a later paper.One can prove existence and uniqueness of solutions at sufficiently small current, under weak conditions on the wake potential.Also, a critique of previous work on ideal models of the wake potential seems to be in order.The models of purely inductive, purely resistive, and purely capacitive wake potentials involve some interesting mathematical questions that should be re-examined.

FIG. 4 .
FIG. 4. Results for the SLC damping ring with improved vacuum chamber.Left: Wake potential WðqÞ; Right: Equilibrium charge density for N ¼ 5 × 10 10 (blue) and in the limit of zero current (red).

9 FIG. 8 .
FIG. 6. Left: Bunch centroid hqi and r.m.s.length σ q as a function of normalized current, for SLC damping ring with improved vacuum chamber.For N ¼ 5 × 10 10 the value of I is 0.140 pC=V.Right: A fit to FðqÞ by a linear combination of resistive and inductive terms, with a ¼ 3.06 and b ¼ 8.14, at N ¼ 5 × 10 10 .