Simulation of the Beam-Beam Effects in $e^+e^-$ Storage Rings with a Method of Reducing the Region of Mesh

A highly accurate self-consistent particle code to simulate the beam-beam collision in $e^+e^-$ storage rings has been developed. It adopts a method of solving the Poisson equation with an open boundary. The method consists of two steps: assigning the potential on a finite boundary using the Green's function, and then solving the potential inside the boundary with a fast Poisson solver. Since the solution of the Poisson's equation is unique, our solution is exactly the same as the one obtained by simply using the Green's function. The method allows us to select much smaller region of mesh and therefore increase the resolution of the solver. The better resolution makes more accurate the calculation of the dynamics in the core of the beams. The luminosity simulated with this method agrees quantitatively with the measurement for the PEP-II B-factory ring in the linear and nonlinear beam current regimes, demonstrating its predictive capability in detail.


Introduction
The beam-beam interaction is one of the most important limiting factors determining the luminosity of storage colliders. It has been studied extensively by theoretical analysis [1], experimental measurements [2], and computer simulations [3]. Historically, due to the complexity of the interaction, many approximations, such as strong-weak [4] or soft-Gaussian [5], have been introduced in order to simulate the interaction in a reasonable computing time. The self-consistent simulation of the beam-beam interaction by solving the Poisson equation with a boundary condition has been proposed first to investigate the round beams [6] and then the flat beams [7]. To enhance the accuracy and to reduce the computational overhead, an algorithm (and a code) of the so-called δf method that can handle strong-strong interactions has been introduced [8]. Another self-consistent approach to the beam-beam interaction is to use the Green's function directly [9] or indirectly [10].
In the present paper we will develop a method that takes advantage from both selfconsistent approaches: a smaller region of mesh from the method of using the Green's function and a faster solver for the interior. In order to develop a highly accurate predictive code at the luminosity saturation region, it is necessary to have a fully self-consistent treatment of field-particle interaction at collision. Since we are interested in simulating the Asymmetric e + e − Storage Collider PEP-II [11], which needs to maximize the luminosity and thus the beam current, it is even more crucial that the beam-beam interaction in the large current regime be treated accurately.
In a self-consistent simulation of the beam-beam interaction in storage rings, the beam distributions have to be evolved dynamically during collision with the opposing beam together with the propagation in the rings. During collision, the beam distributions are used at each time sequence to compute the force that acts on the opposing beam.
Since positrons and electrons are ultra-relativistic particles in high energy storage rings, the beam-beam force is transverse and acts only on the opposing beam. Hence, given a beam distribution, we can divide the distribution longitudinally into several slices and then solve for the two-dimensional force for each slice. Self-consistency is achieved by introducing many-body particles in the field that in turn constitute charge-current, the strategy of the particle-in-cell (PIC) procedure (for example, Ref. [12]). In this paper, for simplicity, we use only a single longitudinal slice for a bunch, ignoring any beam-beam effects encompassing over the length of the bunch.

Method
In modern colliders, beams are focused strongly at the interaction point to achieve high luminosity. As a result the transverse dimension of the beam is much smaller than the dimension of the beam pipe at the collision point. Therefore, the open boundary condition is a good approximation for calculating the transverse beam-beam force.

Green's Function
Given a charge density ρ c (x, y), which is normalized to the total charge dxdyρ c (x, y) = Ne, (2.1) where N is the total number of particles, the electric potential φ(x, y) satisfies the Poisson equation with x and y being the transverse coordinates. The solution of the Poisson equation can be expressed as where G is the Green's function which satisfies the equation In the case of open boundary condition, namely the boundary is far away so that its contribution to the potential can be ignored, one has the well-known explicit solution for the Green's function: This explicit solution can be used directly to compute the potential. The main problem of this approach is that it is slow to calculate the logarithm and the number of computations is proportional to the square of the number of macro particles N 2 p . One can reduce N p by introducing a two-dimensional mesh to smooth out the charge distribution [9]. Or to further improve the computing speed, one can map the solution onto the space of spectrum by the Fast Fourier Transformation (FFT) and then calculate the potential [10].

Reduce the Region of Mesh
Another alternative approach is to solve the Poisson equation with a boundary condition [7], because the region (20 µm × 450 µm for PEP-II) occupied by the beam is much smaller than the boundary defined by the beam pipe (2 cm radius) at the collision point. In order to achieve required resolution, a few mesh points per σ of the beam are needed, otherwise the size of mesh is too large for numerical computation.
However, it is unnecessary to cover the entire area with mesh inside the beam pipe since the area is mostly empty. We choose a smaller and finite area of the mesh, which is large enough to cover the whole beam, and by carefully selecting the potential on the boundary, we can obtain the accurate solution inside the boundary.
We denote by φ 1 the solution (2.3) of the Poisson equation. Let φ 2 be the solution obtained by solving the Poisson equation in a two-dimensional area S with the potential prescribed on a closed one-dimensional L bounding the area S (2.6) where (x, y) ∈ L. By definition, we have φ 1 = φ 2 on the boundary L. Let U = φ 1 − φ 2 and use the first identity of Green's theorem [13] in two dimensions where dl is a line element of L with a unit outward normal n. Since U = 0 on L and ∇ 2 U = 0 inside L, we have S (∇U) 2 dxdy = 0, (2.8) implying that U is a constant inside L. We can set U = 0, which is consistent with the value on the boundary. Hence φ 1 = φ 2 . The two solutions are identical.

Field Solver
We adopt the PIC technique to calculate the fields induced by the charge (and current) of the beams self-consistently. The charge distribution of a beam is represented by macro particles. These macro particles are treated as single electron or positron dynamically. In order to compute the field acting on the particles of the opposing beam, we first deposit their charges onto the gird points of a two-dimensional rectangular mesh. We denote by H x the horizontal distance between two adjacent grid points and by H y the distance in vertical direction.

Charge Assignment
We choose the method of the triangular-shaped cloud [15] as our scheme for the charge assignment onto the grid. On a two-dimensional grid, associated with each macro particle, nine nearest points are assigned with non-vanishing weights as illustrated in Fig. 1. We use "0" to denote the first, "+" as the second, and "-" as the third nearest lines. The weights are quadratic polynomials of the fractional distance, r x = δx/H x , to the nearest line w − The coefficients are chosen such that the transition at the middle of the grid is continuous and smooth, and w 0 x + w + x + w − x = 1 which is required by the conservation of charge. In order to retain these properties, the weights of the two-dimensional grid are simply a product of two one-dimensional weights. For example, w 00 = w 0

Poisson Solver
It is crucial to solve the Poisson equation fast enough (within a second on a computer workstation) for the beam-beam simulation, because the radiation damping time is about 5000 turns and several damping times are needed to reach an equilibrium distribution. For the reason of the computing speed, we follow Krishnagopal [7] and choose the method of cyclic reduction and FFT [14]. A five-point difference scheme is used to approximate the two-dimensional Laplacian operator where i and j are the horizontal and vertical indices that label the grid points on the mesh. Truncation errors are of the order of H 2 x and H 2 y . It is worthwhile to mention that, if we use the same number of mesh points per σ in both transverse directions in the case of beam aspect ratio 30:1, the truncation errors in the horizontal plane are dominant. To minimize the errors in our simulation, we select three times more mesh points per σ in horizontal direction compared to the vertical one.

Field
The field E = −∇φ is computed on the two dimensional grid, using a six-point difference scheme The field off the grid is computed with the same smoothing scheme used in the charge assignment to ensure the conservation of the momentum. The fields E x and E y are interpolated between the grid points. They are calculated by using the weighted summation of the fields at the nine nearest points with exactly the same weights used as the charge is assigned.

Track Particles
The motion of a particle is described by its canonical coordinates where P x and P y are particle momenta normalized by the design momentum p 0 .

One-Turn Map
When synchrotron radiation is turned off, a matrix is used to describe the linear motion in the lattice where M is a 4 × 4 symplectic matrix which can be partitioned into blocks of 2 × 2 matrices when the linear coupling is ignored Here M x , and M y are 2 × 2 symplectic matrices. The matrix M x is expressed with the Courant-Snyder parameters β x , α x , and γ x at the collision point where ν x is the horizontal tune. A similar expression is applied in the vertical plane.

Damping and Synchrotron Radiation
Following Hirata [16], we apply the radiation damping and quantum excitation in the normalized coordinates, since it is easily generalized to include the linear coupling. The motion of a particle in the normalized coordinate is described by a rotation matrix which is obtained by performing the similarity transformation where When synchrotron radiation is switched on, we simply replace the rotation matrix R x with following map in the normalized coordinatesx andP x where ηx and ηp x are Gaussian random variables normalized to unity, τ x is the damping time in unit of number of turns and ǫ x is the equilibrium emittance. In the vertical plane, a similar map is applied.

Beam-Beam Kick
Assuming particles are ultra-relativistic and the collision is head-on, the kick on a particle by the opposing beam is given by the Lorenz force where E x and E y are the horizontal and vertical components of the electric field evaluated at the position of the particle. They are computed with the Poisson solver as outlined in the previous section each time two slices of the beam pass each other. And the half of the transverse force is the magnetic force due the beam moving at the speed of light. The energy of the particle, E 0 = cp 0 , appearing in the denominator of the above expressions comes from the normalization of the canonical momenta P x and P y and the use of the s-coordinate, s = ct, as the "time" variable.  Figure 2: The beam-beam kick by a flat Gaussian beam with aspect ratio 30:1 near X axis and Y axis. The dash-dotted curve is the case when φ = 0 is assigned as the boundary condition. The long-dashed curve is the kick when inhomogeneous boundary condition is used. The short-dashed curve is the kick produced by the Erskine-Bassetti formula [17].
A typical beam-beam kick experienced by a particle near the axis is shown in Fig. 2 with the PEP-II parameters, which are tabulated in the next section. As expected based on the derivation in section 2.2, the kick resulted from solving the Poisson equation with the inhomogeneous boundary condition agrees well with the analytic solution. In addition, the agreement demonstrates that the scheme of the charge deposition works well, the mesh is dense enough and the number of macro particles is large enough.
The number of macro particles used to represent the distribution of the beam is 10240. The area of the mesh is 8σ x ×24σ y and there are 15 grid points per σ x and 5 per σ y . There are about 15 macro particles per cell within 3σ of the beam. These parameters are chosen to minimize truncation errors and maximize resolution. The 256×256 mesh is also the maximum allowed by a computer workstation to complete a typical job within a reasonable time.
The discrepancy between the solution with the homogeneous boundary condition, φ = 0, and the analytic one worsen as the beam aspect ratio becomes larger because the actual change of the potential on the horizontal boundaries becomes larger.

Simulation of PEP-II: Validation
An object-oriented C++ class library has been written to simulate the beam-beam interaction using the method outlined in the previous sections. In the library, the beam and the Poisson solver are all independent objects that can be constructed by the user. For example, there is no limitation on how many beam objects are allowed in the simulation and the beams can have different parameters as an instance of the beam class. These features provide us with great flexibility to study various phenomena of the beam-beam interaction.
We will carry out the simulation of beam-beam interaction with the current operating parameters of the PEP-II so that the results of the simulation can be compared with the known experimental observations. As a goal of this study, after a proper benchmarking of the code against the experiment, we shall be able to make predictions on parameter dependence and show how to improve the luminosity performance of the collider.  The parameters used in the simulation are tabulated in Tab. 5.1. The vertical β * y is lowered to 1.25cm [18] from the design value 1.5cm [11]. The horizontal emittance 24nm-rad in the Low Energy Ring (LER) is half of the design value 48nm-rad because the wiggler is turned off to increase the luminosity. The damping time, 9740 turns, in the LER is a factor of two larger than the one in the High Energy Ring (HER) because of the change of the wigglers made during the construction of the machine. The degradation of luminosity from the increase of the damping time was found then to be about 10% based on the beam-beam simulation. The tunes are split and are determined experimentally to optimize the peak luminosity.

Procedure of simulation
The distribution of the beam is represented as a collection of macro particles that are dynamically tracked. The procedure to obtain equilibrium distributions of the two colliding beams is as follows • initialize the four-dimensional Gaussian distribution according to the parameters of the lattice at the collision point and the emittance of the beam. Distributions of two beams are independent and different.
• iterate a loop with three damping times • propagate each beam through corresponding lattice using one-turn map with synchrotron radiation.
• cast the particle distributions onto the grid as the charge distribution with weighting and smoothing.
• solve for the potential on the grid with the Poisson solver.
• compute the field on the grid.
• calculate the beam-beam kick to the particles of the other beam with the field at the position of the particles. The field off the grid is interpolated with the same weighting and smoothing used in the charge deposition.
• save data such as beam size, beam centroid and luminosity.
• end of the loop.
• save the final distributions.
We vary the beam intensity with a fixed beam current ratio: I + :I − = 2:1, which is close to the ratio for the PEP-II operation. At each beam current, we compute the equilibrium distributions.

Beam-Beam Limit
Given equilibrium distributions that are close enough to the Gaussian, we can introduce the beam-beam parameters where r e is the classical electron radius, γ is the energy of the beam in unit of the rest energy, and N is total number of the charge in the bunch. Here the superscript "+" denotes quantities corresponding to the positron and "-" quantities corresponding to the electron. The results of the simulation are shown in Fig. 3. The beam-beam tune shifts for the electron beam are low because of the large beam-beam blowup of the positron beam. At this operating point, the positron is the weaker beam. When I + = 1200mA and I − = 600mA, which is the near the maximum allowed currents when the beams are in collision, the positron beam sizes are σ + x = 260µm and σ + y = 7µm.

Luminosity
Given the two beam distributions, ρ + and ρ − , the luminosity can be written as where n b is the number of the colliding bunches, f 0 is the revolution frequency, and N + , N − are the number of charges in each position and electron bunch, respectively. Since the distribution ρ is normalized to unity dxdyρ(x, y) = 1 (5.3) and proportional to the charge density ρ c , we evaluate the overlapping integral by a summation over ρ + c ρ − c on the mesh. If we assume the distributions are Gaussian, the overlapping integral can be carried out Two methods agree within a few percents. The mesh method gives a higher luminosity than the Gaussian one. We always use the mesh method, since it can be applied to broad classes of distribution.  Figure 4 shows the luminosity of the beams with 415 colliding bunches, which are spaced with every 8 RF buckets and 10% of the gap. The luminosity is beam-beam limited. It also shows that the optimum number of bunches is between 544 and 665 and the luminosity is about 2.3×10 33 cm −2 s −1 given I + = 1200mA. These results quantitatively agree with the experimental observations in the routine operation of the PEP-II. For example, the peak luminosity of the PEP-II is 1.95×10 33 cm −2 s −1 with I + = 1170mA, I − = 700mA, and n b = 665 during the period of June, 2000. The fact that the luminosity value in the simulation is higher than the observation could be explained by the hour-glass effect which is ignored in the simulation.
For a fixed number of bunches, say 554, the simulation shows a maximum luminosity, which is also seen daily in the control room of the PEP-II. From the simulation, we see that the reason for the peaked luminosity is the rapid growth of σ + y once the peak current is passed.
In addition, the simulation predicts that we can reach the design luminosity 3×10 33 cm −2 s −1 by running 829 bunches at the beam current of I + = 1600mA and I − = 800mA. This prediction has not been realized yet at this time. Currently, the total positron current is probably limited below 1200mA by the electron-cloud instability [19]. Once this limitation is removed, we expect to reach the design luminosity with 829 bunches.
There is no particle loss outside the area (8σ x ×24σ y ) covered by the mesh in the first 15 data points. Beyond the 15th points, particle loss is almost about 1%.

Damping Time
Historically, the damping time is typically not considered to be an important parameter for the beam-beam effects. So we make an attempt to reduce the damping time artificially for the LER to speed up the computation. The result is shown in Fig. 5 The only difference of the parameters used in two simulations is the damping time in the LER, which is indicated as the labels in the figure. Indeed, at the low current, the difference of the luminosity is rather small, which is consistent with the simulation performed when the change of the wiggler was made. But the difference grows larger, as the current increases. At the peak luminosity for the PEP-II operation, I + = 1200mA, the difference is about 40%, which is significant.
This result shows for the first time that the damping time is a rather important parameter for the computation of the peak luminosity at high beam currents. Secondly, it points a way to improve the peak luminosity of the PEP-II without the increase of the beam currents, namely to install another wiggler in the LER to reduce the damping time to the original design value.

Discussion
We have developed a hybrid method of solving the potential with an open boundary by using Green's function to fix the potential on a finite boundary and then to solve the Poisson equation for the potential inside the boundary. The method is applied to the simulation of strong-strong interaction of beam-beam effects in PEP-II. The preliminary results of this simulation show a very good quantitative agreement with the experimental observations. Given the simplicity of the two-dimensional model used, the achievement is surprising and remarkable. We have demonstrated that the present code has a highly reliable predictive capability of realistic beam-beam interaction. To further benchmark the code, we need to extend the simulation to include the finite length of the bunch and compare the simulation results directly to the controlled experiments.
This method is quite general. It can be applied to the problem of space charge in three dimensions. It can also be used in the beam-beam interaction of a linear collider. Finally,  it can be applied to any boundary condition to reduce the region of the mesh if Green's function is known.