Hamiltonian Formalism for Solving the Vlasov-Poisson Equations and Its Applications to Periodic Focusing Systems and the Coherent Beam-Beam Interaction

A Hamiltonian approach to the solution of the Vlasov-Poisson equations has been developed. Based on a nonlinear canonical transformation, the rapidly oscillating terms in the original Hamiltonian are transformed away, yielding a new Hamiltonian that contains slowly varying terms only. The formalism has been applied to the dynamics of an intense beam propagating through a periodic focusing lattice, and to the coherent beam-beam interaction. A stationary solution to the transformed Vlasov equation has been obtained.


INTRODUCTION
The evolution of charged particle beams in accelerators and storage rings can often be described by the Vlasov-Maxwell equations 1,2 . At high energies the discrete-particle collision term 3 comprises a small correction to the dynamics and can be neglected. Radiation effects at sufficiently high energies for leptons can be a significant feature of the dynamics, and should normally be included in the model under consideration.
The Vlasov-Maxwell equations constitute a considerable simplification in the description of charged particle beam propagation. Nonetheless there are only a few cases that are tractable analytically 1,2 . Therefore, it is of the utmost importance to develop a systematic perturbation approach, able to provide satisfactory results in a wide variety of cases of physical interest.
Particle beams are subject to external forces that are often rapidly oscillating, such as quadrupolar focusing forces, RF fields, etc. In addition, the collective self-field excitations can be rapidly oscillating as well. A typical example is a colliding-beam storage ring device, where the evolution of each beam is strongly affected by the electromagnetic force produced by the counter-propagating beam [4][5][6] . The beam-beam kick each beam experiences is localized only in a small region around the interaction point, and is periodic with a period of one turn.
In this and other important applications, one is primarily interested in the long-time behavior of the beam, thus discarding the fast processes on time scales of order the period of the rapid oscillations. To extract the relevant information, an efficient method of averaging is developed in Sec. 2. Unlike the standard canonical perturbation technique 7,8 , the approach used here is carried out in a "mixed" phase space (old coordinates and new canonical momenta) 9 , which is simpler and more efficient in a computational sense. It should be pointed out that the formalism developed here is strictly speaking non-canonical and in

THE HAMILTONIAN FORMALISM
We consider a N-dimensional dynamical system, described by the canonical conjugate pair of vector variables (q, p) with components q = (q 1 , q 2 , . . . , q N ), The Vlasov equation for the distribution function f (q, p; t) can be expressed as where is the Poisson bracket, H(q, p; t) is the Hamiltonian of the system, and summation over repeated indices is implied. Next we define a canonical transformation via the generating function of the second type according to and assume that the Hessian matrix of the generating function S(q, P; t) is non-degenerate, i.e., det H ij = 0.
This implies that the inverse matrix H −1 ij exists. The new canonical variables (Q, P) are defined by the canonical transformation as We also introduce the distribution function defined in terms of the new canonical coordinates (Q, P) and the mixed pair of canonical variables (q, P) according to f 0 (Q, P; t) = f (q(Q, P; t), p(Q, P; t); t), F 0 (q, P; t) = f (q, p(q, P; t); t).
In particular, in Eq. (8) the old canonical variables are expressed in terms of the new ones, which is ensured by the implicit function theorem, provided the relation (6) holds. As far as the function F 0 (q, P; t) is concerned, we simply replace the old momentum p by its counterpart taken from the first of Eqs. (7). Because we can express the Poisson bracket in terms of the mixed variables in the form Differentiation of Eq. (7) with respect to time t, keeping the old variables (q, p) fixed, yields or Our goal is to express the Vlasov equation (2) in terms of the mixed variables (q, P).
Taking into account the identities we obtain Furthermore, using the relation where we express the Vlasov equation in terms of the mixed variables according to where is the new Hamiltonian.
For the distribution function f 0 (Q, P; t), depending on the new canonical variables, we clearly obtain where the new Hamiltonian K is a function of the new canonical pair (Q, P), such that

ING LATTICE
As a first application of the Hamiltonian formalism, we consider the propagation of a continuous beam through a periodic focusing lattice in a circular ring with radius R. Particle motion is accomplished in two degrees of freedom in a plane transverse to the design orbit.
The model equations consist of the nonlinear Vlasov-Poisson equations 1 where is the normalized Hamiltonian, and q = (x, z). The transverse canonical momenta p = (p x , p z ) entering the Hamiltonian (27) are dimensionless variables which represent the actual transverse momenta of the moving particle scaled by the longitudinal momentum of the synchronous particle 10 .
In addition, R is the mean radius of the accelerator and ψ is a normalized potential related to the actual electric potential ϕ according to where N is the total number of particles in the beam, e b is the particle charge, and ε 0 is the electric susceptibility of vacuum. Moreover, the parameter λ is defined by where β s = v s /c is the relative velocity of the synchronous particle, Lorentz factor, and is the classical radius of a beam particle with charge e b and rest mass m b . The coefficients G x,z (θ) determining the focusing strength in both transverse directions are periodic functions with period Θ.
Following the procedure outlined in the preceding section we transform Eqs. (25) - (27) according to where ǫ is formally a small parameter proportional to the focusing field strength, which will be set equal to unity at the end of the calculation. The next step is to expand the quantities S, K and ψ in a power series in ǫ according to We now substitute the expansions (35) -(37) into Eqs. (33) and (34) and obtain perturbation equations that can be solved order by order.
The lowest order solution is evident and has the form First order O(ǫ): Taking into account the already obtained lowest order solutions (38) and (39), the Hamilton-Jacobi equation (33) to first order in ǫ can be expressed as Imposing the condition that the first order Hamiltonian K 1 be equal to we obtain immediately Here we have introduced the notation Note that since the focusing coefficients are periodic functions of θ they can be expanded in a Fourier series where and Ω = 2π/Θ. Therefore for the quantities G x,z and G x,z (θ) expressed in terms of the Fourier amplitudes, we obtain Second order O(ǫ 2 ): To this order, the Hamilton-Jacobi equation (33) takes the form It is straightforward to solve Eq. (48), yielding the obvious result For the second order potential ψ 2 we obtain the equation or, making use of (39), In Eqs. (49) -(51), G x,z (θ) denotes application of the integral operation in Eq. (44) to because G x,z = 0.
Third order O(ǫ 3 ): To third order in ǫ, the Hamilton-Jacobi equation (33) can be written as The third-order Hamiltonian K 3 is given by the expression Equation (53) can be easily solved for the third-order generating function S 3 . The result is For the third-order electric potential ψ 3 we obtain simply Fourth order O(ǫ 4 ): To the fourth order in the expansion parameter ǫ the Hamilton-Jacobi equation (33) can be expressed as The obvious condition to impose is that the fourth-order Hamiltonian K 4 be equal to With Eq. (58) in hand, it is straightforward to solve the fourth-order Hamilton-Jacobi equation (57) for S 4 . We obtain For the fourth-order electric potential ψ 4 , we obtain the Poisson equation Fifth order O(ǫ 5 ): In fifth order, we are interested in the Hamiltonian K 5 . Omitting algebraic details we find In concluding this section, we collect terms up to fifth order in ǫ in the new Hamiltonian K = K 0 + ǫK 1 + ǫ 2 K 2 + . . . and set ǫ = 1. This gives where the coefficients A u , B u and C u are defined by the expressions The Hamiltonian (62), neglecting the contribution from the self-field ψ 0 , describes the unperturbed betatron oscillations in both horizontal and vertical directions.
It is useful to compute the unperturbed betatron tunes ν x,z in terms of averages over the focusing field-strengths. For a Hamiltonian system governed by a quadratic form in the canonical variables of the type in Eq. (62), it is well-known that the characteristic frequencies ν x,z can be expressed as Keeping terms up to sixth order in the perturbation parameter ǫ, we obtain In terms of Fourier amplitudes of the focusing coefficients, the Eq. (67) can be expressed as For purposes of illustration, we consider a simple F ODO lattice with equal focusing and defocusing strengths +G and −G, and period Θ. We also assume that the longitudinal dimensions θ f of the focusing and defocusing lenses are equal; the longitudinal dimensions θ d of the corresponding drift spaces are assumed to be equal as well. Moreover, For simplicity we consider the horizontal degree of freedom only (the vertical one can be treated in analogous manner). The Fourier amplitudes of the focusing coefficients are where n = 0, 1, 2, . . .. To second order in ǫ, we obtain for the horizontal betatron tune In the limit of infinitely thin lenses, θ f → 0, Eq. (71) reduces to the well-known expression where use of the identity has been made.
It is evident from Eqs. (68) and (71), that the Hamiltonian averaging technique developed here represents a powerful formalism for evaluating the betatron tunes in terms of averages over the focusing field strength.

COHERENT BEAM-BEAM INTERACTION
As a second application of the Hamiltonian formalism developed in Sec. 2, we study here the evolution of two counter-propagating beams, nonlinearly coupled by the electromagnetic interaction between the beams at collision. For simplicity, we consider one-dimensional motion in the vertical (q) direction, described by the nonlinear Vlasov-Poisson equations where is the Hamiltonian. Here λ k is the beam-beam coupling parameter, defined according to 11 Moreover, (k = 1, 2) labels the beam, f k (q, p; θ) is the distribution function, θ is the azimuthal angle, and ν k is the betatron frequency in vertical direction. In addition, R is the mean machine radius, r e is the classical electron radius, N 1,2 is the total number of particles in either beam, V k (q; θ) is the normalized beam-beam potential, β * kq is the vertical beta-function at the interaction point, and L kx is the horizontal dimension of the beam ribbon 12 .
Our goal is to determine a canonical transformation such that the new Hamiltonian is time-independent. As a consequence, the stationary solution of the Vlasov equation (21) where ǫ is again a formal small parameter, which will be set equal to unity at the end of the calculation.
The next step is to expand the quantities S k , K k and V k in a power series in ǫ, analogous to Eqs. (35) -(37), according to where First Order: O(ǫ) where Third Order: O(ǫ 3 ) In third order we are interested in the new Hamiltonian, which is of the form where ζ(z) is Riemann's zeta-function 13 where Integrating Eq. (94) over P we obtain a nonlinear integral equation of the Haissinski type 14 for the equilibrium beam density profile ̺ where Here sgn(z) is the well-known sign-function.
Let us further specify the function G k (K k ) and assume that it is given by the thermal equilibrium distribution 1,10,15 where N k is a normalization constant, defined according to and ε k is the unnormalized beam emittance. The second term in the Hamiltonian (97) can be transformed according to where Expanding the beam density ̺ (3−k) 0 (q 1 + q) occurring in the integral in Eq. (102) in a Taylor series and integrating by parts, we obtain where Substituting (100) and (104) into Eq. (96) we obtain Taking into account that where H n (z) is the Hermite polynomial 13 of order n, we obtain where and Here, Φ(z) is the error function 13 .
In order to determine the normalization constant(s), N k , we utilize the method of Laplace to take the integral of the beam density ̺ These are evidently maxima, since Integrating the beam density (109) over q, we obtain 16 Equation (114) represents two transcendental equations for determining the normalization constants N k . For the beam centroid and the beam size, i.e., the first and the second moments of the beam density (109), we obtain In order to proceed further, we assume that the beam-beam coupling parameter λ k is small, and expand the equilibrium beam density ̺ where ̺ (k) and The main goal in what follows is to determine the normalization constant(s) N k0 . To do so we integrate Eq. (117) over q. As a result of simple algebraic manipulations, we obtain Introducing the new unknowns we can write the two equations for determining M 1,2 as From Eq. (122), as a result of simple algebraic manipulations we obtain the quadratic for M 1 , and the equation for determining M 2 once M 1 is known. Equation (124) has one real double root if and only if the discriminant is equal to zero. This gives Since the scaled normalization constants M 1,2 should be positive we choose Thus we obtain To conclude this section we note that in the case of D = 0 we have two solutions for either M, i.e., Note also that the discriminant D is invariant (does not change) under permutation of b 1 and b 2 . In other words, four different physically realizable situations are possible for a wide range of parameters The inequality in Eq. (131) has been obtained under the condition that both solutions in Eq. (130) are positive. This case corresponds to the so-called "flip-flop" state 17 of the two colliding beams, which is a bifurcated state better to be avoided.

CONCLUSIONS
We have developed a systematic canonical perturbation approach that removes rapidly oscillating terms in Hamiltonians of quite general form. The essential feature of this approach is the use of mixed canonical variables. For this purpose the Vlasov-Poisson equations are transformed to mixed canonical variables, and an appropriate perturbation scheme is chosen to obtain the equilibrium phase space density. It is worthwhile to note that the perturbation expansion outlined in the preceding section can be carried out to arbitrary order, although higher-order calculations become very tedious.
In conclusion, it is evident from the present analysis that the Hamiltonian averaging Finally, we reiterate that the mixed variable Hamiltonian formalism developed in the present analysis can be used to derive amplitude equations, describing processes of formation of patterns and coherent structures in a number of plasma and beam systems in which collective processes are important.

ACKNOWLEDGMENTS
We are indebted to S.A. Heifets for many fruitful discussions concerning the subject of the present paper. It is a pleasure to thank H. Qin for illuminating discussions and comments.
This research was supported by the U.S. Department of Energy.