A New Model for the Collective Beam-Beam Interaction

The Collective Beam-Beam interaction is studied in the framework of maps with a ``kick-lattice'' model in the 4-D phase space of the transverse motion. A novel approach to the classical method of averaging is used to derive an approximate map which is equivalent to a flow within the averaging approximation. The flow equation is a continuous-time Vlasov equation which we call the averaged Vlasov equation, the new model of this paper. The power of this approach is evidenced by the fact that the averaged Vlasov equation has exact equilibria and the associated linearized equations have uncoupled azimuthal Fourier modes. The equation for the Fourier modes leads to a Fredholm integral equation of the third kind and the setting is ready-made for the development of a weakly nonlinear theory to study the coupling of the pi and sigma modes. The pi and sigma eigenmodes are calculated from the third kind integral equation. These results are compared with the kick-lattice model using our weighted macroparticle tracking code and a newly developed, density tracking, parallel, Perron-Frobenius code.


Introduction
In this paper we introduce a new model for the collective beam-beam interaction for hadron beams in 4D transverse phase space (2 degrees-of-freedom). This both generalizes and simplifies the work of [1,2,3,4] on the collective beam-beam interaction in high energy colliders. In addition, it extends the preliminary 1 degree-of-freedom collective case in [5]. Our model is based on the classical method of averaging generalized to maps and collective forces. We do not distribute the beam-beam force around the ring as is usually done. The technique we introduce should be of general interest for studies of Vlasov systems with a localized perturbative collective force.
In Section 2, we discuss our basic kick-lattice model for the evolution of the 4D phase space densities of the two beams. In Section 3, we briefly review the basic averaging theory which is generalized in this paper. Previously, this averaging theory was applied to the weak-strong beam-beam in one and two degrees-of-freedom [6,7]. The equations of the kick-lattice model will be transformed to a standard form for the method of averaging in Section 4 and the general "averaged Vlasov equation" (AVE) will be derived (See equation (24)). We then introduce the special case we treat in this paper, namely the case where the tunes of the two beams are identical and are non-resonant. In this case the AVE has the property that any function of the action only is an equilibrium solution. In Section 5, we linearize about these equilibria and discuss the linearized equations and the associated third kind integral equation. In addition, we compare our integral equation with the analogous integral equations which arise in the standard plasma problem and in the beam dynamics problems concerning the longitudinal dynamics with wake fields for a coasting beam and a bunched beam. In Section 6 we present numerical results for the π and σ mode eigen-problems for an axially symmetric Gaussian equilibrium and compare these results with simulations on the exact model of Section 2. In Section 7 we give a summary and point to future work. An appendix is included which gives a first principles calculation of the beam-beam force.

The Kick-Lattice Model in 4D Phase Space
To describe the evolution equations, we refer to the bunches as "unstarred" and "starred", and for every quantity X describing the unstarred bunch, the quantity X * describes the starred bunch. The evolution equations are symmetric: the equation for the starred bunch is obtained from the unstarred bunch by interchanging starred and unstarred quantities, so we mostly state only the equation for the unstarred bunch.
We consider two counter-rotating particle bunches, which collide head on at a single interaction point (IP). The electromagnetic interaction at the IP is determined up to a proportionality factor by the dimensionless "potential" Φ, which satisfies the Poisson equation −∆Φ = 2πρ * . Here ρ * (x, y) is the spatial density (normalized to one), and the potential Φ[ρ * ] : R 2 → R is given by where G(x, y) = − ln( x 2 + y 2 /σ) = − ln( x 2 + y 2 ) + ln(σ) is the Green's function. In the following we will omit the scaling factor σ which is in principle needed for dimensional correctness but which can be chosen completely arbitrarily since it does not contribute to the beam-beam kick.
Letting n refer to the state of the system just before the IP, particles in the unstarred bunch change their phase-space position u = (x, y, p x , p y ) T according to the map The associated phase space density ψ n evolves via ψ n+1 (u n+1 ) = ψ n (u n ), or which is easily inverted to give Here M is a stable linear symplectic map representing the linear lattice, ζ is the beam-beam factor, P = " 1 0 0 0 0 1 0 0 « projects phase space on configuration space, and the spatial and phase-space densities are related by ρ * (x, y) = R 2 ψ * (x, y, p x , p y ) dp x dp y .
The beam-beam factor, which is derived in the appendix, is ζ = 1+ββ * β+β * 2N * βγ r p , where the absolute value of r p = qq * 4πǫ 0 mc 2 is the classical particle radius (as long as only elementary particles or ions of the same charge state are involved), N is the number of particles, q is the particle charge, γ is the Lorentz factor associated with β, and m is the particle mass. For all modern colliders, i.e. in the limit β,β * → 1, ζ can be approximated by ζ ≈ 2N * γ r p . The evolution law for the starred beam is obtained by replacing M by M * ,ψ * by ψ and ζ by ζ * where starred and unstarred are interchanged in β, γ, N and m.
Equation (2) can be written more compactly as where J 2k = " 0 k I k −I k 0 k « is the unit symplectic matrix. We note here that a map is said to be symplectic if the Jacobian, M, of the map satisfies M T J M = J . We have written the kick in a "Hamiltonian form" because eventually a transformed Φ will be a Hamiltonian for a flow.

Map Averaging and Error Bounds
Here we give an overview of the averaging formalism, which we generalize in this paper, and briefly discuss error bounds. We consider the autonomous "kick-rotate" map in R 2 with the small parameter ǫ. This is a model for the one degree of freedom weak-strong beambeam interaction and was discussed in [6]. The transformation leads to the non-autonomous map where H(v, θ) = V (v 1 cos(2πθ) + v 2 sin(2πθ)). This is in a standard form for the method of averaging in which the transformed dependent variable, v, is slowly varying. If ν is irrational, then from Weyl's equidistribution theorem [8] the average of H(v, nν) over n exists and is given byH It is therefore natural to ask, for what values of ν are solutions of (10) approximated by solutions of the averaged map Even though the maps in (8) and (10) are symplectic, the averaged map is not. However the averaged flow associated with (11) and defined bẏ is Hamiltonian, and it is easy to show that |w n − w(n)| = O(ǫ) over O(1/ǫ) times. Since (12) is autonomous, (11) can be viewed as the Euler method for numerically integrating (12). Approximating equation (10) with (11) is considered in our previous work [6], where we introduce the concept of a far-from-low-order-resonance zone for ν. This zone is formed by removing a finite number of intervals centered on low-order rationals, therefore ν needs to satisfy only finitely many Diophantine conditions, and does not need to be irrational, which makes the formalism much more useful in the applications. The error bound |v n −w n | = O(ǫ) is obtained without the usual near-identity-transformation and is unchanged asymptotically if an O(ǫ 2 ) term is added to equations (8,10). In [7], we extend the formalism of equations (8)(9)(10)(11)(12) to the weak-strong beam-beam interaction in 2 degrees-of-freedom. The equation of motion corresponding to (8) is just equation (2) with ρ * n replaced by the spatial density of the strong beam. We are working out the details of the averaging theorem in this more complicated case with two frequencies [9]. Some ingredients of our approach can be found in [10]. We generalize this to the collective beam-beam interaction in the next section. Map-Averaging for Vlasov Systems We will begin by transforming (2) using a representation of the solution to the unperturbed, ζ = 0, problem. The new coordinates will be slowly varying if ζ is small. As in the previous section we could proceed by letting u = M n q which gives The averaged equation then becomes whereF [ρ * ](w) denotes the n-average of M −n J 4 ∇ u Φ[ρ * ](PM n w). The transformation to (13) turned out to be a major advance in Sobol's implementation of the Perron-Frobenius (PF) method [11] (See [12] and [13] for a discussion of the PF method). In addition, (13) is well suited for an error analysis which is in progress [14]. However, an action-angle transformation may be better suited to understand approximate equilibria and the associated linear analysis that we do here and that is how we will proceed.
The action-angle transformation from u = (x, y, p x , p y ) T to slowly varying coordinates v = (Θ x , Θ y , J x , J y ) T is given by Note that for fixed J and Θ these are solutions of the equations of motion with ζ = 0, that is without the beam-beam force.
Equation (2) becomes where The integral in (20) is taken over [0, 2π] in the Θ's and over [0, ∞) in the J's. Since the transformation is symplectic it is also volume preserving. Thus the (Θ, J)-density is given by Ψ n (v) = ψ n (u), and its evolution law is Clearly v n and v * n are slowly varying for ζ and ζ * small, and it follows that the transformed densities Ψ and Ψ * are slowly varying. Thus (19) is in a standard form for averaging and we now follow the procedure laid out in the previous section. The averaged map problem is obtained from (19) by replacing H by the appropriate n-averageH and dropping the O(ζ 2 ) term. The associated averaged flow problem is autonomous and has the Hamiltonian forṁ Thus the averaged Vlasov equations for Ψ and Ψ * become  (20) by interchanging the starred and unstarred parameters β and ν. System (24) is the new model referred to in the title. Since ζ and ζ * are small one immediate advantage of (24) over (2) is that the step size in a numerical integration of (24) can be O(1/max(ζ, ζ * )) which is much larger than one turn.
At this stage the problem is general with parameters (ν x , ν y , ν * x , ν * y , β x , β y , β * x , β * y , ζ, ζ * ) and the correct averaged HamiltonianH depends on the relation between the four tunes. Here we discuss the case ν x = ν * x and ν y = ν * y because (i) we wish to compare and contrast our results with [3,4] and (ii) it simplifies the calculation of the average. In this case G(· · ·) in (20) can be rewritten as and the phases ϕ x and ϕ y are easily determined from the trigonometry involved. If in addition we consider the case where ν x and ν y are non-resonant (in the sense that k x ν x + k y ν y = k 0 ⇒ k x = k y = k 0 = 0), then the averaging over nν x and nν y can be done separately and each average can be replaced by the associated integral. Thus the averaged Hamiltonian becomes whereḠ(D x , D y ) := 1/(2π) 2 T 2 G(D x sin t x , D y sin t y ) dt x dt y and as before v = (Θ, J). From our experience with the non-collective case [6,7], we expect this to be valid for ν x , ν y far from low-order resonances; work on the error estimates is in progress [14]. Note that H = H * if β x = β * x and β y = β * y . Because of the convolution structure of the Θ integral in (27) (see (25)) any function Ψ * = Ψ * e (J) results inH being independent of Θ. It follows that any pair of densities Ψ e and Ψ * e that are independent of Θ are an equilibrium pair for (24). Since ψ n (x, y, p x , p y ) = Ψ n (Θ, J) ≈ Ψ e (J) = Ψ e 1 2 , we see that (2) and equivalently (19) have quasiequilibria given the averaging approximation. In [5], we have verified this in the 2D phase space case. Figure 1 shows the evolution of two action densities at a hundred out of more than 10 5 turns under the exact map, the 2D analog of (4), using the PF method for tracking phase space densities. The parameters are ν = ν * = √ 5 − 2 and ξ = 3 × 10 −3 . The red crosses represent the evolution of a quasi-equilibrium, namely the centered Gaussian Ψ 0 (Θ, J) = Ψ e (J) = 1/(2πε)e −J/ε (for ε = 1), while the green ×'s stand for a Gaussian that was initially shifted by 1σ x , giving a Θ-dependent density. Since the red crosses for different discrete J lie on top of each other, Ψ e hardly evolves over the 10 5 turns which is consistent with the averaging result. The Θ-dependent density however, fluctuates as is indicated by the band of green ×'s. Thus we have strong evidence for the existence of quasi-equilibria, however we believe that (4) does not admit exact equilibria. This is in contrast to the lepton case, see [15].

The Linearized Equations
To linearize about an equilibrium, we set Ψ Plugging into (24) and linearizing gives Because of the convolution structure ofH in (27), the Fourier modes of (28) are uncoupled, but we do not pursue (28) in that generality here. Instead we investigate solutions, which exist when Ψ e = Ψ * e , ζ = ζ * , andH =H * . Letting f ± = Ψ 1 ± Ψ * 1 , we obtain where f ± are the so-called σ (sum) and π (difference) modes respectively and we have scaled the ζ into the time by defining τ = ζt. We remind the reader that the densities Ψ 1 and Ψ * 1 in (28) are coupled, while f + and f − in (29) are not. Moreover, the σ and π modes can be intuitively interpreted as in-phase and 180 • out-of-phase perturbations respectively to the starred and unstarred beam. Note that t in (28) is dimensionless (same as n in (22)) whereas τ in (29) has the dimension of a length. Note also that the beam current, which is a factor of ζ, has been scaled out of the problem. Thus there can be no threshold for instability in our problem.
Unfolding the Poisson brackets we obtain where To analyze equation (30), we expand To obtain an equation for the Fourier coefficients, we define with Fourier coefficients given by Since K is real and even in Θ x and Θ y , K k is real, and since K k is real, An easy calculation using (35) and (27) gives Thus the Fourier modes determined by (29) are uncoupled and are given by where f e (J) := (2π) 2 Ψ e (J) is the equilibrium action density ( R 2 + f e (J)dJ = 1). There are two standard approaches to analyzing (37): the Laplace transform approach and the eigenvalue approach.
Taking the Laplace transform of (37) we obtain for ℑω sufficiently large. Here we use −iω instead of the usual Laplace variable s,f k (J, ω) denotes the Laplace transform of f k , and Ω(J) :=H(Ψ e )(J) ′ . In the eigenvalue approach, we look for solutions of the form which gives The Laplace transformed equation (38) has a non-homogeneous term, however, the left hand sides of equations (38) and (40) are identical. These equations are of the form and are called Fredholm integral equations of the third kind (see p.2 [16]). If a is bounded away from zero it can be transformed into a Fredholm integral equation of the second kind. Thus the primary interest in this third kind equation is in the case where a has zeros and this is our case. The case where a has zeros is complicated by the fact that there are generalized solutions which are difficult to represent numerically. Equations (38) and (40) have analogues in both plasma physics and other beam dynamics contexts. For example, Crawford and Hislop [17] discuss the standard plasma problem in the periodic case, the case of this paper, summarizing both the Landau and the van Kampen-Case solutions ( [18], [19], [20]). Jackson [21] gives a nice presentation of the Landau approach in the non-periodic case. The third kind integral equations are given in equations (17) and (23) of [17] and in equation (3.5) of [21]. As is well known, the plasma problem leads to Landau damping and growth for certain equilibria depending on the size of the average density. The stability analysis is facilitated by the dispersion function which uncouples the calculation of the poles of the solution from the calculation of the density itself.
Two standard beam dynamics problems concern the longitudinal dynamics with wake fields for a coasting beam and for a bunched beam. The coasting beam case is completely analogous to the periodic plasma problem including a dispersion function and possible Landau damping for small beam current and an instability threshold at some critical current after which there is Landau growth. A recent discussion of the coasting beam problem in the context of coherent synchrotron radiation is given in [22]. The third kind integral equation is given in equation (27) of that paper. The emphasis in [22] is on the threshold for instability which occurs when a zero of the dispersion function reaches the real axis as the current increases from small values. Landau damping is not discussed as it is not important for the stability discussion and furthermore would require an analytic continuation into the lower half ω plane which would require assumptions on the equilibrium (see p.7 in [22]).
The bunched beam case is more complicated as the Fourier modes do not decouple. Furthermore, it appears at first sight that the calculation of the instability threshold must be done in combination with the calculation of the density. However, Warnock has introduced a regularization transformation which eliminates the continuous spectrum. The resulting equation is then discretized leading to a determinant condition, independent of the density, which is analogous to the dispersion relation. A convergence theorem would then make this rigorous. This is discussed in [23], where equation (11) is very similar to our equation (38). However the kernel of the integral equation is much different and, in fact, one expects a stability threshold. It is in this paper that Warnock introduces his regularization transformation which eliminates the continuous spectrum and thus eliminates the numerical problem of representing generalized functions numerically. More recent progress on the regularization is given in [24].
Our equations (38) and (40) are simpler than the longitudinal bunched beam equations in that the Fourier modes are uncoupled. Also, our case is rather special in that it does not depend on the beam current as mentioned above. In fact the π and σ eigen-modes are neutrally stable if k · ∇f e (J) = 0. In this case, the transformation leads to whereK k (J, J ′ ) := |k · ∇f e (J)| 1/2 K k (J, J ′ )|k · ∇f e (J ′ )| 1/2 and η ± := ± sgn(k · ∇f e ). Since the kernelK k is real and symmetric, an eigenvalue ω must be real. It is in fact a remarkable feature of the linearized AVE (43) for the case of ν = ν * , (ν x , ν y ) far-from-low-order-resonance and for equilibria with k · ∇f e (J) = 0, that despite the presence of an amplitude dependent tune shift and a collective force, the modes show neither damping nor growth but instead are stable. In [26,13] we have given numerical evidence that the modes are indeed remarkably stable even in the fully nonlinear regime of tracking with equation (2). We have tried to show that eigensolutions (ω,f ± k ) of (40) for equilibria, which do not satisfy the above condition, k · ∇f e (J) = 0 (e.g. densities with two humps), must have ω ∈ R, but have been unsuccessful. This leaves open the possibility of complex eigenvalues. Since these must come in complex conjugate pairs, there is the possibility of linearly (and thus nonlinearly) unstable solutions.
We have assumed that the two beams have the same nonresonant tunes and this is probably the reason that the eigensolutions are neutrally stable. Previous work indicates that when the tunes are different (a so-called tune split) or near-to-low-order resonance there can be Landau damping or growth. In [25], the authors develop a perturbation procedure in the near resonance case and argue that there are regions of stability and instability (see Figure 2 of [25] ). Landau damping is also discussed in [4]. In [26], we have seen evidence for Landau damping in the 2D phase space case (See Figures 5-7 of [26] ). A future problem for us is to determine the explicit form of the averaged equations in this case. Equation (43) and the associated Laplace transformed equation can be rewritten as (T − ωI)ĝ ± = h. Since T is symmetric, we are looking for conditions which ensure that it is a bounded selfadjoint operator on an appropriate Hilbert space. Such operators have a well developed spectral theory. For example, the spectrum is a compact subset of the real line contained in the interval [m, M] where both m and M are finite spectral values and all spectral values are either in the point spectrum (eigenvalues) or the continuous spectrum, thus the residual spectrum is empty [27]. Numerical results, in the section to follow, indicate that for k T = (1, 0) or (0, 1) the spectrum is the interval [0, k · Ω(0)] for the σ mode with 0 an eigenvalue and [0, k · Ω(0)] ∪ {ω π } for the π mode, with ω π > k · Ω(0) an eigenvalue. A possible explanation for the stability of the modes using the notion of the "Landau resonance" is that the modes can not resonantly exchange energy with a macroscopic fraction of the particles in the beam. The σ mode tune lies at the edge of the incoherent (continuous) spectrum towards infinite orbital amplitudes. Any sensible phase space density falls off rapidly at large amplitudes (or even has compact support) so that the fraction of particles with tunes in resonance with the σ mode tune vanishes. The π mode tune, in the studied parameter regime, lies considerably outside the incoherent spectrum and is thus even less able to dissipate its energy among the single particle trajectories or to draw energy from them.

Numerical Results for π and σ Modes
The analysis of the spectrum for equation (43) gives important information about the π and σ modes. In this section, we discuss properties of solutions of (43) and give our results on the numerical solution of this eigenvalue problem.
If ω is outside of the range of k · Ω(J), (43) can be reduced to an integral equation of the second kind by a simple algebraic transformation. Conversely, if ω is in the range of k · Ω(J), then (43) it must be treated as a third kind equation. Such equations have not been studied as extensively as Fredholm integral equations of the first and second kind. A review of work up to 1973 and new results are given in [28] and more recent results are contained in [29]. Most recently, we have become aware of [30,31,32]. However, to our knowledge, the case when J is 2-dimensional has not been discussed nor have convergent numerical schemes been developed. As mentioned in the previous section the plasma problem and the coasting beam problems are of this type and have been studied extensively, however these are particularly simple.
We now discuss the commonly used, straightforward discretization for integral equations of the third kind as applied to our special case and give our numerical results. At the end of this section we will discuss progress on work toward a convergent scheme. In the straightforward approach, J is put on a mesh and the integral is approximated by a simple quadrature method. This leads to a finite dimensional matrix eigenproblem and seems to lead to reasonable results. This approach has been used in 1D in the beam-beam interaction in [3,4], in the longitudinal bunched beam case in [33] and by us in the 1D beam-beam interaction, [13]. In [13] we obtained excellent agreement between the FFT spectra of the dipole modes in full blown simulations and the eigenvalues of a one degree-of-freedom version of (43).
We consider the special case of axially symmetric Gaussian beams, where with πε being the rms emittance, and horizontal dipole modes where k = (1, 0) T . With these choices equation (43) takes the explicit form We transform the actions I x = J x /(1 + J x ) and I y = J y /(1 + J y ), thereby mapping R + → [0, 1), and use a 60×60 mesh. Even though the linearized averaged Vlasov equation (45) reduces the number of independent variables from four as in (24) to two, the evaluation of the functions Ω(J) andK k (J, J ′ ) is in fact computationally expensive. The computation of Ω = ∇H[Ψ e ] involves a 6-fold integral at each point of the 2D mesh in J andK k involves a 4-fold integral at each point of a 4D mesh in (J, J ′ ). Although we found a way to slightly simplify the calculation for general Ψ e , and reduce the 6-fold integral to a 5-fold integral, going to larger meshes is quite expensive. However, for the important particular choice of (44) we found a very simple formula involving modified Bessel functions:H This formula has been known in the context of the weak-strong approximation of the beambeam tune shift, [34]. It is straightforward to prove that lim J→0 Ω x (J) = lim J→0 Ω y (J) = 1, and that the ranges of Ω x (J), Ω y (J) are both the interval (0, 1). The latter is also the range of the continuous spectrum of (43) . Ω x (J x , J y ) is shown in Figure 2, and the spectrum of the finite dimensional approximation of (45) is shown in Figure 3. The plot suggests that (45) has a continuous spectrum, common to both the σ and π modes, which coincides with the range of Ω x (J). In addition, the σ mode has a discrete eigenvalue ω = 0, and the π mode has a discrete eigenvalue at ω ≈ 1.21.   Figure 3: The dipole spectrum of (45). The black curve is the continuum, common to both modes, the red circle is the discrete eigenvalue for the π mode and the blue diamond is the discrete eigenvalue for the σ mode.   Figure 6: An eigenfunction corresponding to an eigenvalue from the continuous spectrum of (45). and 5 show that these eigenvalues corresponds to regular eigenfunctions. FFT spectra obtained by tracking, with our newly written parallel Perron-Frobenius code [11], which tracks the phase space densities directly in 4D phase space, and with our 4D weighted macro particle tracking code [13], shows pronounced peaks at tunes that correspond to these discrete eigenvalues. This indicates an excellent agreement between these three completely different approaches. We consider this a major feat. Even though our results are quite satisfactory, we are interested in a convergent procedure. A well established theory (See e.g., [35]) guarantees that the above described finite dimensional approximation of the integral equation converges as long as the operator is compact. However, as we mentioned before, Figure 3 indicates the presence of a continuous spectrum, which a compact operator can not have. In addition, the numerically computed "eigenfunctions" associated with the continuous spectrum show a singular behavior, as illustrated in Figure 6. This is to be expected [28,29] and thus these third kind integral equations are not well suited for numerical analysis in their standard form. Therefore further numerical and analytical analysis requires a special treatment of equation (43). Such difficulties are common for equation (41) where a(x) vanishes at least once inside of its domain. The recent work for the longitudinal bunched beam Vlasov equation in [23,24] mentioned above, suggests that the continuous spectrum can be eliminated by the Warnock regularization transformation g(J, ω) = (ω − k · Ω(J))ĝ(J, ω). This gives the nonlinear eigenvalue problem in ω: where the solutions (g (J, ω), ω) for η = ±1 give the σ and π eigenmodes respectively. As mentioned before, equation (41) appears in other physics applications, and analytical results have been published in [28,29,30,31,32]. We have begun a study of this problem in our particular case. Specifically, the existence and uniqueness of solutions of equation (43) in the one and two-dimensional cases has been addressed in [36]. In [36], we consider functions which are continuous except for pole-like singularities at J such that (ω − k · Ω(J)) = 0, and interpret the integrals over the singularities in the principal value sense. Under certain assumptions, we proved a version of the Fredholm alternative theorem: the equation has a unique solution for any right-hand side iff the homogeneous version of this equation has only the trivial solution.
In addition, this framework may provide a numerically convergent scheme for solving (43). Alternatively, we are looking for conditions so that discretization of the Warnock transformed problem will lead to a convergent method. The theory of [36] is a first step.

Summary and Future Work
An important aspect of collective beam-beam theory has been the study of the so-called π and σ modes. The pioneering works of [1,2,3,4] represent a major advance; however approximations are made, the validity of which we would like to understand better. For example, in [3,4], the starting point is a Vlasov equation with a delta function kick and the beam-beam kick is distributed around the ring by smoothing the delta function. The Vlasov equation is linearized about a function which is only an approximate equilibrium and action-angle variables are introduced. The Fourier modes in angle are not uncoupled at this stage but a horizontal dipole mode proportional to exp(iΘ x ) is assumed. This leads to an inconsistent equation. To obtain a consistent equation for this mode an average over the angle variables is taken which, after a Fourier transform in time, leads to an integral equation, the analog of equation (43). In contrast, we start with the kick-lattice model to properly handle the delta function kick. Then we make only one approximation, the averaging approximation, and in addition, as stated earlier, we believe we can give an upper bound on the error of approximation. Our AVE has exact equilibria and thus our exact problem has quasi-equilibria in good agreement with our simulations. The linearization about these equilibria leads to an equation, which in contrast to the above, has uncoupled Fourier modes. The Fourier modes satisfy an integral equation that is easily transformed to a formally selfadjoint problem. We are looking for conditions such that the associated operator is bounded and selfadjoint, a case which has a well developed theory. The standard computation of the π and σ dipole mode frequencies discussed in Section 6 is in good agreement with density tracking based on equation (4), using both the PF and the weighted macro-particle tracking methods. However the standard numerical approach to numerically solving (43) does not converge as the mesh size decreases beyond some limit and we are searching for convergent algorithms such as that suggested in [23,36].
In summary, we have introduced a new model, the averaged Vlasov equation (24), for the collective beam-beam interaction in two degrees of freedom, which we believe has significant potential for deepening our understanding of this important collective effect. Equation (24) was derived in the spirit of the rigorous analysis in [6]. We believe similar error bounds can be derived, thus we believe (24) gives a good approximation to the basic dynamics of (2). In fact, we have checked the one degree-of-freedom analog of (24) with two aspects of a full-blown density tracking approach, the existence of quasi-equilibria and the calculation of the π and σ mode eigenvalues, with excellent results. More importantly, we have checked the two degree-of-freedom AVE with two full-blown simulation codes and have also found excellent agreement in the calculation of π and σ mode eigenvalues. Thus we have confidence in the model.
We have demonstrated its usefulness as a tool for calculating π and σ mode eigenvalues and for clarifying the existence of quasi-equilibria. In the case of leptons, progress has been made on the question of the existence of an equilibrium for the exact model [15]. However, it seems likely that exact equilibria do not exist in the hadron case as the underlying dynamics is likely to be chaotic. In addition, the AVE may lead to a faster algorithm for calculating the density evolution. This is because the beam-beam parameters, ζ and ζ * , are small and thus the time step in numerical integration of (24) can be O(1/max(ζ, ζ * )), which is significantly larger than one turn. We propose to investigate this potential speed up by developing a numerical procedure to integrate (24). As another example, the AVE will be useful in taking the next step beyond the linear theory to investigate coupling between the π and σ Fourier modes. In the Laplace-picture, we may be able to use the fixed point iteration scheme discussed for the plasma problem in [37] or in the eigen-picture presented here we may be able to use the van Kampen-Case approach [19,20]. In the latter case, the work of [17] may be useful. Finally, we can investigate several other effects such as those discussed by Alexahin [4]. Some topics we are considering for future work are: (i) a study of the near-to-low-order resonance case as we do in [6] (ii) adding another degree of freedom to study the effect of synchrotron motion into the dynamics of equation (2) and (iii) a study of the effect of a tune split by letting ν * x = ν x + ζa and ν * y = ν y + ζb and then applying our averaging formalism. Here (a, b) allow us to vary the tune split in units of the kick parameter. As in [6] we expect bifurcations as a and b vary. Items (i) and (iii) likely includes the possibility of both Landau damping and growth.
Our main point is that we now have a model in which many important collective beam-beam interaction effects can be studied in a more systematic way then was previously available.
we will consider the relative velocities of the particles with respect to the synchronous particle of their bunch as non-relativistic.
We consider the momentum change of an unstarred particle moving through the starred bunch. In the kick approximation, we assume that the transverse spatial coordinates are not changed during the interaction, i.e. r(t) = (x, y, z + ut) T . Thus the change in transverse momentum of an unstarred particle passing through the starred beam measured in F * is where u is the speed of the unstarred particle in the starred frame, and E * ⊥ is the transverse electric field of the starred bunch in F * . Note that β is basically the speed of an unstarred particle in the lab frame and that B * is approximately zero in F * .
Since the 3-momentum P * is part of a 4-vector and the boosts involved are all in the longitudinal direction ∆P ⊥ = ∆P ′ ⊥ = ∆P * ⊥ .
One easily shows by direct integration that which is independent of z. Thus where ρ * 2 = R ρ * 3 dz. The interchange of limits going from (51) to (52) is justified for ρ * 3 decaying sufficiently fast at ∞. (The singularities in the integral are integrable.) Since ∆P ′ ⊥ = ∆P * ⊥ the kick in the lab frame is Since the speed of the unstarred synchronous particle in the lab frame is βc, we have p 0 = mβγc. Thus the kick is ∆p where ζ = 2N * qq * 4πmǫ 0 c 2 1 βγ which is the justification for equation (2).