A Symplectic Multi-Particle Tracking Model for Self-Consistent Space-Charge Simulation

Symplectic tracking is important in accelerator beam dynamics simulation. So far, to the best of our knowledge, there is no self-consistent symplectic space-charge tracking model available in the accelerator community. In this paper, we present a two-dimensional and a three-dimensional symplectic multi-particle spectral model for space-charge tracking simulation. This model includes both the effect from external fields and the effect of self-consistent space-charge fields using a split-operator method. Such a model preserves the phase space structure and shows much less numerical emittance growth than the particle-in-cell model in the illustrative examples.


I. INTRODUCTION
In high intensity accelerators, the nonlinear space-charge effect from charged particle interactions inside the beam has significant impact on beam dynamics through the accelerator. It causes beam emittance growth, halo formation, and even particle losses along the accelerator. To study the space-charge effect, multi-particle tracking has been employed to dynamically follow those charged particles through the accelerator. In the accelerator community, most of those multi-particle tracking codes use particle-in-cell (PIC) method to include the space-charge effect self-consistently in the simulation [1][2][3][4][5][6][7][8][9][10][11].
The particle-in-cell method is an efficient method in handling the space-charge effect self-consistently. It uses a computational grid to obtain the charge density distribution from a finite number of macroparticles and solves the Poisson equation on the grid at each time step. The computational cost is linearly proportional to the number of macroparticles, which makes the simulation fast for many applications. However, those grid based, momentum conserved, PIC codes do not satisfy the symplectic condition of classic multi-particle dynamics. Violating the symplectic condition in multi-particle tracking might not be an issue in a single pass system such as a linear accelerator. In a circular accelerator, violating the symplectic condition may result in undesired numerical errors in the long-term tracking simulation. This issue together with the numerical grid heating was brought up during the 2015 space-charge workshop at Oxford [12]. A gridless spectral based macroparticle model was suggested by the author at the workshop to mitigate the numerical grid heating and to satisfy the symplectic condition of particle tracking.
Multi-symplectic particle-in-cell model was proposed to study Vlasov-Maxwell system and electrostatic system in plasmas using a variational method [13][14][15][16][17]. To study the space-charge effect in high intensity beams, a quasi-static model is normally employed. In the quasi-static model, a moving beam frame is used to contain all charged particles through the accelerator. The Poisson equation is solved in the beam frame to obtain electric Coulomb fields from the charged particles. These electric fields are transformed to the laboratory frame through the Lorentz transformation. The space-charge forces acting on each individual particle include both the electric fields and the magnetic fields, which is different from the electrostatic model that includes only electric fields. To the best of our knowledge, at present, there is no symplectic self-consistent space-charge model available in the accelerator community. In this paper, following the idea suggested at the Oxford workshop, we present a two-dimensional and a three-dimensional symplectic quasi-static multi-particle tracking model for space-charge simulations. The model presented here starts from the multi-particle Hamiltonian directly and uses a gridless spectral method to calculate the space-charge forces.
The organization of this paper is as follows: after the introduction, we present the symplectic multi-particle tracking model including the space-charge effect in Section II; We present a symplectic space-charge transfer map for a 2D coasting beam in Section III and a symplectic space-charge map for a 3D bunched beam in Section IV; We discuss computational complexity of the proposed model in Section V and draw conclusions in Section VI.

II. SYMPLECTIC MULTI-PARTICLE TRACKING WITH SPACE-CHARGE EFFECTS
In the accelerator beam dynamics simulation, for a multi-particle system with N p charged particles subject to both a space-charge self field and an external field, an approximate Hamiltonian of the system can be written as [18][19][20]: where H(r 1 , r 2 , · · · , r Np , p 1 , p 2 , · · · , p Np ; s) denotes the Hamiltonian of the system using distance s as an independent variable, ϕ is related to the space-charge interaction potential between the charged particles i and j (subject to appropriate boundary conditions), ψ denotes the potential associated with the external field, r i = (x i , y i , θ i = ω∆t) denotes the normalized canonical spatial coordinates of particle i, p i = (p xi , p yi , p ti = −∆E/mC 2 ) the normalized canonical momentum coordinates of particle i, and ω the reference angular frequency, ∆t the time of flight to location s, ∆E the energy deviation with respect to the reference particle, m the rest mass of the particle, and C the speed of light in vacuum. The equations governing the motion of individual particle i follows the Hamilton's equations as: Let ζ denote a 6N-vector of coordinates, the above Hamilton's equation can be rewritten as: where [ , ] is the Poisson bracket. A formal solution for above equation after a single step τ can be written as: Here, we have defined a differential operator : H : as : H : g = [H, g], for arbitrary function g. For a Hamiltonian that can be written as a sum of two terms H = H 1 + H 2 , an approximate solution to above formal solution can be written as [21] ζ(τ ) = exp(−τ (: Let exp(− 1 2 τ : H 1 :) define a transfer map M 1 and exp(−τ : H 2 :) a transfer map M 2 , for a single step, the above splitting results in a second order numerical integrator for the original Hamilton's equation as: Using the above transfer maps M 1 and M 2 , a fourth order numerical integrator can also be constructed as [21]: where α = 1 − 2 1/3 , and s = τ /(1 + α). An even higher order accuracy integrator can be obtained following Yoshida's approach [22]. Assume that M 2n denotes a transfer map with an accuracy of order 2n, the tranfer map M 2n+2 with (2n + 2)th order of accuracy can be obtained from the recursion equation: where z 0 = 1/(2 − 2 1/(2n+1) ) and z 1 = −2 1/(2n+1) /(2 − 2 1/(2n+1) ). The above numerical integrator Eqs. 7-9 will be symplectic if both the transfer map M 1 and the transfer map where J denotes the 6N × 6N matrix given by: and I is the 3N × 3N identity matrix. For the given Hamiltonian in Eq. 1, we can choose H 1 as: A single charged particle magnetic optics method can be used to find a symplectic transfer map M 1 for this Hamiltonian with the external fields from most accelerator beam line elements [19,20,23].
We can choose H 2 as: which includes the space-charge effect and is only a function of positions. For the space-charge Hamiltonian H 2 (r), the single step transfer map M 2 can be written as: The Jacobi matrix of the above transfer map M 2 is where L is a 3N ×3N matrix. For M 2 to satisfy the symplectic condition Eq. 10, the matrix L needs to be a symmetric matrix, i.e.
Given the fact that L ij = ∂p i (τ )/∂r j = − ∂ 2 H2(r) ∂ri∂rj τ , the matrix L will be symmetric as long as it is analytically calculated from the function H 2 . This is also called jolt-factorization in nonlinear single particle beam dynamics study [24]. If both the transfer map M 1 and the transfer map M 2 are symplectic, the numerical integrator Eqs. 7-9 for multi-particle tracking will be symplectic. In the following sections, we will derive the symplectic space-charge transfer map of H 2 for a two-dimensional coasting beam and for a three-dimensional bunched beam.

III. SYMPLECTIC SPACE-CHARGE MAP FOR A COASTING BEAM
In a coasting beam, the Hamiltonian H 2 can be written as [19]: where K = qI/(2πǫ 0 p 0 v 2 0 γ 2 0 ) is the generalized perveance, I is the beam current, ǫ 0 is the dielectric constant in vacuum, p 0 is the momentum of the reference particle, v 0 is the speed of the reference particle, γ 0 is the relativistic factor of the reference particle, and ϕ is the space charge Coulomb interaction potential. In this Hamiltonian, the effects of the direct electric potential and the longitudinal vector potential are combined together. The electric Coulomb potential in the Hamiltonian H 2 can be obtained from the solution of the Poisson equation. In the following, we assume that the coasting beam is inside a rectangular perfect conducting pipe. In this case, the two-dimensional Poisson's equation can be written as: where φ is the electric potential, and ρ is the particle density distribution of the beam. The boundary conditions for the electric potential inside the rectangular conducting pipe are: where a is the horizontal width of the pipe and b is the vertical width of the pipe. Given the boundary conditions in Eqs. 20-23, the electric potential φ and the source term ρ can be approximated using two sine functions as [25][26][27][28][29]: where where α l = lπ/a and β m = mπ/b. The above approximation follows the numerical spectral Galerkin method since each basis function satisfies the boundary conditions on the wall [25][26][27]. For a smooth function, this spectral approximation has an accuracy whose numerical error scales as O(exp(−cN )) with c > 0, where N is the number of the basis function (i.e. mode number in each dimension) used in the approximation. By substituting above expansions into the Poisson Eq. 19 and making use of the orthonormal condition of the sine functions, we obtain where γ 2 lm = α 2 l + β 2 m . In the multi-particle tracking, the particle distribution function ρ(x, y) can be represented as: where δ is the Dirac function. Using the above equation and Eq. 26 and Eq. 28, we obtain: and the electric potential as: From the above electric potential, the interaction potential ϕ between particles i and j can be written as: The one-step symplectic transfer map M 2 of the particle i with this Hamiltonian is given as: Here, both p xi and p yi are normalized by the reference particle momentum p 0 . Using the symplectic transfer map M 1 for the external field Hamiltonian H 1 from an optics code and the transfer map M 2 , one obtains a symplectic multi-particle tracking model including the self-consistent space-charge effect following Eqs. 7-9.
As an illustration of above symplectic multi-particle tracking model, we simulated a 1 GeV coasting proton beam transporting through a rectangular perfect conducting pipe with a FODO lattice for transverse focusing. The initial transverse density distribution is assumed to be a Gaussian function given in Fig. 1. We computed the electric field along the x axis using the above direct gridless spectral solver with 15 × 15 modes and the electric field from a second order finite difference solver with 129 × 129 grid points. The results are shown in Fig. 2. The solution of the spectral solver agrees with that of the finite difference solver very well with even 15 × 15 modes due to the fast convergence property of the spectral method. The relative maximum field difference (normalized by the maximum field amplitude) between two solutions is below 2%. The use of the mode number in this example is somewhat empirical. It depends on the physical problem to be solved. If one knows about the smallest spatial structure of the problem, one can choose the mode number with the wavelength to resolve this spatial structure. Without knowing the detailed structure in the particle density distribution, one can use a trial-and-error method until the appropriate solution is attained. Figure 3 shows the proton beam root-mean-square (rms) envelope evolution through 20 FODO lattice periods. The FODO lattice used in this example consists of two quadrupoles and three drifts in a single period. The total length of the period is 1 meter. The zero current phase advance is about 87 degrees and the phase advance with 100 A current is about 74 degrees. The relatively low intensity beam used in this example is to avoid space-charge driven resonance and to separate the numerical emittance growth from the physical emittance growth in the simulation. The symplectic integrator is good for long term tracking since it helps preserve phase space structure during the numerical integration. Figure 4 shows the stroboscopic plots (every 10 periods) of x − p x and y − p y phase space evolution of a test particle through the last 20, 000 periods of the total 100, 000 lattice periods including the selfconsistent space-charge forces. As a comparison, we also show in this figure the phase space evolution of the same initial test particle using the standard momentum conserved PIC method and the second order finite difference solver for space-charge calculation [19,30]. Qualitatively, these two models show similar shapes in phase space. However, looking into the details of the phase space, they have quite different structures. The single particle phase space from the PIC model shows a dense core while the phase space from the symplectic multi-particle model shows a nearly hollow core. Figure 5 shows the 4-dimensional emittance growth ( ǫx ǫx0 ǫy ǫy0 − 1)% evolution from the symplectic model and that from the PIC model. It is seen that the symplectic model has a much smaller emittance growth than the PIC model. This emittance growth is a numerical artifact due to the small number of macroparticles (50, 000) used in the simulation, which was studied in references [31][32][33]. The small number of macroparticles introduces numerical errors in the computing of the electric potential and results in the artificial emittance growth. A more detailed study of the numerical emittance growth associated with this new method is under way and will be reported in future publication. The apparent non-zero emittance growth at the beginning is due to the charge redistribution of the initial Gaussian distribution within a much shorter time scale (not visible in the plot) compared with the total plotting time scale of 100, 000 periods.

IV. SYMPLECTIC SPACE-CHARGE MAP FOR A 3D BUNCHED BEAM
In a 3D bunched beam, the Hamiltonian H 2 can be written as [34]: where κ = q/(lmC 2 γ 2 0 β 0 ), l = C/ω is the scaling length, and β 0 = v 0 /C. The above Hamiltonian includes both the electric potential and the longitudinal magnetic vector potential. The electric potential in the beam frame can be obtained from the solution of a three-dimensional Poisson's equation: where ρ is the charge density distribution in the beam frame. The boundary conditions for the electric potential inside the rectangular perfect conducting pipe are: where a is the horizontal width of the pipe, b is the vertical width of the pipe. The solution of the 3D Poisson equation subject to the above boundary conditions was studied with several numerical methods [28,29]. To obtain a fast analytical solution of the electric potential we use an artificial boundary condition in this study. Here, the longitudinal open boundary condition is approximated by a finite domain Dirichlet boundary condition: φ(x, y, z = 0) = 0 (43) φ(x, y, z = c) = 0 where c is the length of the domain that is large enough so that the electric potential goes to zero at both ends of the domain. The choice of the length of the domain depends on how fast the electric potential vanishes outside the beam. From the reference [29], we know that the solution of the electric potential for each transverse mode can be written as: where γ 2 lm = (lπ/a) 2 + (mπ/b) 2 . This solution decreases exponentially as a function of z outside the beam. This suggests that a short distance (in the unit of aperture size) might be sufficient to have the electric potential approach to zero.
Given the boundary conditions in Eqs. 37-44, the electric potential φ and the source term ρ can be approximated using three sine functions as: ρ(x, y, z) = where α l = lπ/a, β m = mπ/b, γ n = nπ/c. Substituting the above expansions into the Poisson Eq. 36 and making use of the orthonormal condition of the sine functions, we obtain where Γ 2 lmn = α 2 l + β 2 m + γ 2 n . In the multi-particle tracking, the charge density ρ(x, y, z) can be represented as: where w is the charge weight of each individual particle and δ is the Dirac function. Using the above equation and Eq. 48 and Eq. 50, we obtain: and the electric potential as: Nn n=1 1 Γ 2 lmn sin(α l x j ) sin(β m y j ) sin(γ n z j ) sin(α l x) sin(β m y) sin(γ n z) From the above electric potential, we obtain the interaction potential between particles i and j as: Nn n=1 1 Γ 2 lmn sin(α l x j ) sin(β m y j ) sin(γ n z j ) sin(α l x i ) sin(β m y i ) sin(γ n z i ) (54) Given particle's spatial coordinates (x i , y i , θ i ) in the laboratory frame, the particle spatial coordinates in the beam frame (x i , y i , z i ) can be obtained under the following approximation: Now, the space-charge Hamiltonian H 2 can be written as: Nn n=1 1 Γ 2 lmn sin(α l x j ) sin(β m y j ) sin(−γ n lγ 0 β 0 θ j ) sin(α l x i ) sin(β m y i ) sin(−γ n lγ 0 β 0 θ i ) (56) The one-step symplectic transfer map M 2 of the particle i with this Hamiltonian is given as: where both p xi and p yi are normalized by mC.
As an illustration of above symplectic model, we simulated a 1 GeV, 3D bunched proton beam transporting through a periodic focusing channel. The initial transverse and longitudinal density profiles of the beam are shown in Fig. 6. The beam has a 3D Gaussian distribution with a longitudinal to transverse aspect ratio of three. Figure 7 shows the relative transverse electric field difference along the x-axis and the longitudinal electric field difference along the z-axis using above gridless spectral method with 15 × 15 × 15 modes and using a spectral-finite difference solver [28] with 129 × 129 × 257 grid points. It is seen that even with only 15 modes in each direction, the gridless spectral solver produces space charge fields in good agreement with the fields from the spectral-finite difference solver with finer resolution. The relative maximum field differences (normalized by the maximum field amplitude along the axis) are below 2% in both directions. Figure 8 shows the transverse and the longitudinal rms envelope evolution through a periodic focusing channel. Each period of the focusing channel consists of two transverse uniform focusing elements, two longitudinal uniform focusing elements, and four drifts. The total length of the period is one meter. The zero current phase advance in a single period is about 86 degrees in the transverse dimension and 40 degrees in the longitudinal direction. The phase advance with 0.1 A average current at 100 MHz RF frequency is about 81 degrees in the transverse direction and 39 degrees in the longitudinal direction. Figure 9 shows the six-dimensional rms emittance growth ( ǫx ǫx0 ǫy ǫy0 ǫz ǫz0 − 1)% evolution through the periodic focusing channel from the above symplectic gridless spectral model and from the standard PIC method with the spectral-finite difference Poisson solver. It is seen that the symplectic spectral model gives much less numerical emittance growth than the standard PIC method in the simulation using 160, 000 macroparticles.

V. COMPUTATIONAL COMPLEXITY
The gridless symplectic multi-particle spectral model can be used for long-term tracking study including spacecharge effects. The computational complexity of this model scales as O(N mode × N p ), where N mode is the total number of modes. The standard PIC model can have a computational cost of O(N p ) + O(N grid logN grid ) when an efficient Poisson solver is used, where N grid is the total number of grid points. This suggests that the PIC model would be faster than the symplectic multi-particle spectral model on a single processor computer. However, the symplectic multi-particle spectral model is very easy to be parallelized on multi-processor computer. One can distribute all macroparticles uniformly across processors to achieve a perfect load balance. By using a spectral method with exponentially decreasing errors, the number of modes N mode can be kept within a relatively small number, which significantly improves the computing speed. Figure 10 shows the parallel speedup of the symplectic multi-particle spectral model as a function of the number of processors for a fixed problem size, i.e. ∼ 50, 000 macroparticles and 15 × 15 modes in the 2D model and ∼ 160, 000 macroparticles and 15 × 15 × 15 modes in the 3D model. It is seen that the speedup increases almost linearly for both models. This is because both models have perfect load balance among all processors. The only communication involved in these models is a global reduction operation to obtain the density distribution in the frequency domain. The decrease of the speedup in the 2D case might be due to the specific computer architecture used in this timing study, which has 24 shared memory computing cores inside a node. Outside the node, the communication among processors (cores) is slowed down due to the across node communication. The above scaling results show that the symplectic multi-particle spectral space-charge tracking model can have a good scalability on multi-processor parallel computers and is especially suitable for tracking simulations on large scale supercomputers or GPU computers.

VI. CONCLUSIONS
In this paper, we proposed a new symplectic multi-particle tracking model for self-consistent space-charge simulation. This model uses a gridless spectral method to calculate the space-charge potential and while avoiding the error associated with numerical grid in the standard PIC model. It also shows much less numerical noise driven emittance growth than the PIC method for long term simulation. Even though the computational cost of the symplectic spectral model is higher than the PIC method on a single processor computer, the proposed model scales well on multi-processor parallel computers. It has a perfect load balance and uniform data structure, which is suitable for GPU parallel implementation. The new symplectic multi-particle spectral model enables researchers to carry out long term tracking studies including space-charge effects.
The symplectic space-charge transfer map presented in this paper assumes a rectangular perfect conducting pipe. A transverse open boundary condition might be approximated using this model by moving the conducting wall away from the beam. For a general boundary condition, it is quite difficult to obtain an analytical expression of the electric potential from an arbitrary density distribution. For a round perfect conducting pipe, a Fourier mode and a Bessel mode might be used to approximate the particle density distribution and the electric potential of the Poisson equation in a cylindric coordinate system [28]. However, it takes more time to compute the Bessel function expansion than the simple sine function expansion.
The symplectic space-charge model presented here also assumes a straight conducting pipe. A study of the solution of the Poisson equation in a bended conducting pipe using the Frenet-Serret coordinate was done in reference [35] and shows that for a large normalized bending radius (bending radius/transverse aperture size), e.g. 100, there is barely any difference between the straight pipe solution and the bended pipe solution. This condition (large normalized bending radius) can be satisfied in most circular accelerators. Thus, the symplectic space-charge model in this paper can still be used for space-charge simulation in circular machines.