A symplectic particle-in-cell model for space-charge beam dynamics simulation

Space-charge effects play an important role in high intensity particle accelerators and were studied using a variety of macroparticle tracking models. In this paper, we propose a symplectic particle-in-cell (PIC) model and compare this model with a recently published symplectic gridless particle model and a conventional nonsymplectic PIC model for long-term space-charge simulation of a coasting beam. Using the same step size and the same number of modes for the space-charge Poisson solver, all three models show qualitatively similar emittance growth evolutions and final phase space shapes in the benchmark examples. Quantitatively, the symplectic PIC and the symplectic gridless particle models agree with each other very well, while the nonsymplectic PIC model yields different emittance growth value. Using finer step size, the emittance growth from the nonsymplectic PIC converges towards that from the symplectic PIC model.


I. INTRODUCTION
The nonlinear space-charge effects from particle interactions inside a charged particle beam play an important role in high-intensity accelerators. A variety of models have been developed to simulate the space-charge effects [1][2][3][4][5][6][7][8][9][10][11][12][13]. Among those models, the particle-in-cell (PIC) model is probably the most widely used model in the accelerator community. The particle-in-cell model is an efficient model to handle the space-charge effects 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 multiparticle dynamics. Violating the symplectic condition in multiparticle tracking might not be an issue in a single pass system such as a linear accelerator. But in a circular accelerator, violating the symplectic condition may result in undesired numerical errors in the long-term tracking simulation. Recently, a symplectic gridless particle model was proposed for selfconsistent space-charge simulation [14]. This model is easy to implement and has a good scalability on massive parallel computers. However, the computational cost of this model is proportional to the number of macroparticles times the number of modes used to solve the Poisson equation. When the number of modes is large, this model becomes quite computationally expensive and one has to resort to parallel computers.
In this paper, we propose a symplectic PIC model for self-consistent space-charge long-term simulation. This model has the advantages of the computational efficiency of the PIC method and the symplectic property needed for long-term dynamics tracking. The multisymplectic particle-in-cell model was proposed to study the Vlasov-Maxwell system in plasmas using a variational method [15][16][17][18]. A symplectic gridless spectral electrostatic model in a periodic plasma was also proposed following the variational method [19]. To the best of our knowledge, at present, there is no symplectic PIC model that was proposed and shown applications for the quasistatic Vlasov-Poisson system in the space-charge beam dynamics study. In this paper, we derived a symplectic PIC model directly from the multiparticle Hamiltonian. We then carried out long-term space-charge simulations and compared this model with the symplectic gridless particle model and the nonsymplectic PIC model in two beam dynamics applications.
The organization of this paper is as follows: After the Introduction, we present the symplectic multiparticle tracking model including the space-charge effect in Sec. II. We present the nonsymplectic PIC model in Sec. III. In Sec. IV, we compare the symplectic PIC model, the symplectic gridless particle model, and the nonsymplectic PIC model with two space-charge simulation examples. Finally, in Sec. V we draw conclusions.

II. SYMPLECTIC MULTIPARTICLE TRACKING MODELS
In the accelerator beam dynamics simulation, for a multiparticle 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 [20][21][22] H ¼ where Hðr 1 ; r 2 ; …; r N p ; p 1 ; p 2 ; …; p N p ; sÞ denotes the Hamiltonian of the system using distance s as an independent variable,φ is the space-charge interaction potential (including both the direct electric potential and the longitudinal vector 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 follow 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 the above formal solution can be written as [23] ζðτÞ ¼ expð−τð∶H 1 ∶ þ ∶H 2 ∶ÞÞζð0Þ Letting 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 The above numerical integrator can be extended to fourth-order accuracy and arbitrary even-order accuracy following Yoshida's approach [24]. This numerical integrator Eq. (7) will be symplectic if both the transfer map M 1 and the transfer map M 2 are symplectic. A transfer map M i is symplectic if and only if the Jacobian matrix M i of the transfer map M i satisfies the following condition: where J denotes the 6N × 6N matrix given by and I is the 3N × 3N identity matrix. For the Hamiltonian in Eq. (1), we can choose H 1 as This corresponds to the Hamiltonian of a group of charged particles inside an external field without mutual interaction among themselves. A single-charged particle magnetic optics method can be used to find the symplectic transfer map M 1 for this Hamiltonian with the external fields from most accelerator beam line elements [21,22,25]. We can choose H 2 as which includes the space-charge effect and is only a function of position. 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. (8), the matrix L needs to be a symmetric matrix, i.e.
Given the fact that L ij ¼ ∂p i ðτÞ=∂r j ¼ − ∂ 2 H 2 ðrÞ ∂r i ∂r j τ, 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 the nonlinear single particle beam dynamics study [26]. If both the transfer map M 1 and the transfer map M 2 are symplectic, the numerical integrator Eq. (7) for multiparticle tracking will be symplectic.
For a coasting beam, the Hamiltonian H 2 can be written as [21] Þ is the generalized perveance, I is the beam current, ϵ 0 is the permittivity of 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 Coulomb 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 perfectly conducting pipe. In this case, the twodimensional 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 perfectly 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. (18)-(21), the electric potential ϕ and the source term ρ can be approximated using two sine functions as [27][28][29][30][31] ρðx; yÞ ¼ X N l l¼1 X N m m¼1 ρ lm sinðα l xÞ sinðβ m yÞ; ð22Þ where ρðx; yÞ sinðα l xÞ sinðβ m yÞdxdy; ð24Þ 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 [27][28][29]. 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 the above expansions into the Poisson Eq. (17) and making use of the orthonormal condition of the sine functions, we obtain where γ 2 lm ¼ α 2 l þ β 2 m . In the simulation, the particle density distribution function ρðx; yÞ can be represented as where SðxÞ is the unitless shape function (also called deposition function in the PIC model) and Δx and Δy are mesh size in each dimension. The use of the shape function can help smooth the density function when the number of macroparticles in the simulation is much less than the real number of particles in the beam. Using the above equation and Eqs. (24) and (26), we obtain and the electric potential as The electric potential at a particle i location can be obtained from the potential as where the interpolation function to the particle location is assumed to be the same as the deposition function. From the above electric potential, the interaction potential φ between particles i and j can be written as Now, the space-charge Hamiltonian H 2 can be written as A. Symplectic gridless particle model In the symplectic gridless particle space-charge model, the shape function is assumed to be a Dirac delta function and the particle distribution function ρðx; yÞ can be represented as Now, the space-charge Hamiltonian H 2 can be written as The one-step symplectic transfer map M 2 of the particle i with this Hamiltonian is given as sinðα l x j Þ sinðβ m y j Þ cosðα l x i Þ sinðβ m y i Þ; sinðα l x j Þ sinðβ m y j Þ sinðα l x i Þ cosðβ m y i Þ: Here, both p xi and p yi are normalized by the reference particle momentum p 0 .

B. Symplectic particle-in-cell model
If the deposition/interpolation shape function is not a delta function, the one-step symplectic transfer map M 2 of the particle i with this space-charge Hamiltonian H 2 is given as Sðy − y i Þ sinðα l xÞ sinðβ m yÞdxdy; where both p xi and p yi are normalized by the reference particle momentum p 0 , and S 0 ðxÞ is the first derivative of the shape function. Assuming that the shape function is a compact local function with respect to the computational grid, the integral in the above map can be approximated as where the integers I, J, I 0 , and J 0 denote the two-dimensional computational grid index, and the summations with respect to those indices are limited to the range of a few local grid points depending on the specific deposition function. If one defines the density related functionρðx I 0 ; y J 0 Þ on the grid as the space-charge map can be rewritten as ðx I 0 ; y J 0 Þ sinðα l x I 0 Þ sinðβ m y J 0 Þ sinðα l x I Þ sinðβ m y J Þ : It turns out that the expression inside the bracket represents the solution of potential on grid using the charge density on grid, which can be written as ðx I 0 ; y J 0 Þ sinðα l x I 0 Þ sinðβ m y J 0 Þ sinðα l x I Þ sinðβ m y J Þ: ð40Þ Then the above space-charge map can be rewritten in a more concise format as In the PIC literature, some compact function such as linear function and quadratic function is used in the simulation. For example, a quadratic shape function can be written as [32,33] The same shape function and its derivative can be applied to the y dimension. Using the symplectic transfer map M 1 for the external field Hamiltonian H 1 from a magnetic optics code and the transfer map M 2 for space-charge Hamiltonian H 2 , one obtains a symplectic PIC model including the selfconsistent space-charge effects.

III. NONSYMPLECTIC PARTICLE-IN-CELL MODEL
In the nonsymplectic PIC model, for a coasting beam, the equations governing the motion of individual particle i follow the Lorentz equations as where r ¼ ðx; yÞ, p ¼ ðv x =v 0 ; v y =v 0 Þ, and a z is a unit vector in the longitudinal z direction. The above equations can be solved using a secondorder leap-frog algorithm, which is widely used in many dynamics simulations. Here, for each step, the positions of a particle i are advanced a half step by Then the momenta of the particle are advanced a full step including the forces from the given external fields and the self-consistent space-charge fields. In the PIC model, the space-charge fields are obtained by solving the Poisson equation on a computational grid and then interpolating back to the particle location. In this study, we employ the above spectral method to attain the space-charge fields as where the charge density ρðx I 0 ; y J 0 Þ on the grid is obtained following the charge deposition Eq. (38). The particle momenta can then be advanced one step using p xi ðτÞ ¼p xi ð0Þþτ qE ext Here, both the direct electric Coulomb field and the magnetic field from the longitudinal vector potential are included in the space-charge force term of the above equation. Then, the positions of the particle are advanced another half step using the new momenta as This procedure is repeated many steps during the simulation.

IV. BENCHMARK EXAMPLES
In this section, we tested and compared the symplectic PIC model with the symplectic gridless particle model and the nonsymplectic PIC model using two application examples. In both examples, the beam experienced strong space-charge driven resonance. It was tracked including self-consistent space-charge effects for several hundred thousands of lattice periods using the above three models.
In the first example, we simulated a 1 GeV coasting proton beam transport through a linear periodic quadrupole focusing and defocusing (FODO) channel inside a rectangular perfectly conducting pipe. A schematic plot of the lattice is shown in Fig. 1 lattice period is 85°. The current of the beam is 450 A and the depressed phase advance is 42°. Such a high current drives the beam across the fourth-order collective instability [34]. The initial transverse normalized emittance of the proton beam is 1 μm with a 4D Gaussian distribution. Figure 2 shows the four-dimensional emittance growth ð ϵ x ϵ x0 ϵ y ϵ y0 − 1Þ% evolution through 200 000 lattice periods from the symplectic gridless particle model, from the symplectic PIC model, and from the nonsymplectic spectral PIC model. These simulations used about 50 000 macroparticles and 15 × 15 modes in the spectral Poisson solver. In both PIC models, 257 × 257 grid points are used to obtain the density distribution function on the grid. The perfectly conducting pipe has a square shape with an aperture size of 10 mm. It is seen that the symplectic PIC model and the symplectic gridless particle model agree with each other very well. The nonsymplectic spectral PIC model yields significantly smaller emittance growth than those from the two symplectic methods, which might result from the numerical damping effects in the nonsymplectic integrator. The fast emittance growth within the first 20 000 periods is caused by the space-charge driven fourth-order collective instability. The slow emittance growth after 20 000 periods might be due to numerical collisional effects. Figure 3 shows the transverse phase space distributions at the end of the simulation (200 000 periods) from these three models. It is seen that the distributions from all three models show similar phase space shapes after the strong collective resonance. However, the core density distributions show somewhat different structures. The two symplectic models show similar structures while the nonsymplectic spectral PIC model shows a denser core. This is also seen in the final horizontal and vertical density profiles shown in Fig. 4.
The accuracy of the nonsymplectic PIC model can be improved with finer step size. Figure 5 shows the 4D emittance growth evolution from the symplectic PIC model and those from the nonsymplectic PIC model with the same nominal step size, from the nonsymplectic PIC model with one-half of the nominal step size, and from the nonsymplectic PIC model with one-quarter of the nominal step size. It is seen that as the step size decreases, the emittance growth from the nonsymplectic PIC model converges towards that from the symplectic PIC model.
In the second example, we simulated a 1 GeV coasting proton beam transport through a nonlinear lattice. A schematic plot of the lattice is shown in Fig. 6. Each  Figure 7 shows the 4D emittance growth evolution with 10 A, 20 A, and 30 A currents from the symplectic PIC model, the symplectic gridless particle model, and the nonsymplectic PIC model. The space-charge tune shift with a 10 A beam current is 0.038, with a 20 A is 0.075, and with a 30 A is 0.113. With a 10 A current, particles stay out of the major third-order resonance; there is little emittance growth. With a 20 A current, the proton beam moves near the third resonance and some particles fall into the resonance, resulting in a significant increase of the emittance growth. With a 30 A current, a lot of particles fall into the resonance resulting in large emittance growth. It is seen from the above plot that all three models predict this current dependent behavior qualitatively. The two symplectic models show very similar emittance growth evolution for all three currents while the nonsymplectic PIC model shows much less emittance growth compared with the other two models. This is probably due to the numerical damping effects in the nonsymplectic PIC model. Figure 8 shows the final transverse phase space distributions after 40 000 turns. All three models show similar phase space shapes. The two symplectic models show closer density distribution than the nonsymplectic PIC model. Figure 9 shows the final horizontal and vertical density profiles from those three models. The two symplectic models yield very close density profiles with non-Gaussian halo tails. The nonsymplectic PIC model yields a denser core but less halo near the tail density distribution.
The accuracy of the nonsymplectic PIC model can also be improved with a finer step size in the above test case with nonlinear resonance. Figure 10 shows the 4D emittance growth evolution from the symplectic PIC model and those from the nonsymplectic PIC model with the same nominal step size, from the nonsymplectic PIC model with one-half of the nominal step size, and from the nonsymplectic PIC model with one-quarter of the nominal step size. It is seen that as the step size decreases, the emittance growth from the nonsymplectic PIC model converges towards that from the symplectic PIC model even in this nonlinear resonance case.

V. CONCLUSIONS
In this paper, we proposed a symplectic PIC model and compared this model with the symplectic gridless particle model and the conventional nonsymplectic PIC model for long-term space-charge simulations of a coasting beam in a linear transport lattice and a nonlinear lattice. Using the same step size and the same number of modes for spacecharge solver, all three models qualitatively predict the similar emittance growth evolution due to the space-charge driven resonance and yield similar final phase space shapes in the benchmark examples. Quantitatively, the symplectic PIC and the symplectic gridless particle model agree with FIG. 9. Final horizontal (left) and vertical (right) density profiles from the symplectic gridless particle model (red), the symplectic PIC model (green), and the nonsymplectic spectral PIC (blue).
FIG. 10. Four-dimensional emittance growth evolution from the symplectic PIC model (red), the nonsymplectic spectral PIC with the nominal step size (green), the spectral PIC with half of the nominal step size (blue), and the spectral PIC with one quarter of the nominal step size (pink) in the nonlinear lattice with a 30 A beam current. each other very well, while the nonsymplectic PIC model yields much less emittance growth value. Using finer step size, the emittance growth from the nonsymplectic PIC converges towards that from the symplectic PIC model.
The computational complexity of the symplectic PIC model is proportional to the O½N grid logðN grid Þ þ N p , where N grid and N p are total number of computational grid points and macroparticles used in the simulation. The computational complexity of the symplectic gridless particle model is proportional to the OðN mode N p Þ, where N mode is the total number of modes for the space-charge solver. On a single processor computer, with a large number of macroparticles used, the symplectic PIC model is computationally more efficient than the symplectic gridless particle model. However, on a massive parallel computer, the scalability of the symplectic PIC model may be limited by the challenges of work load balance among multiple processors and communication associated with the gridbased space-charge solver and particle manager [35]. The symplectic gridless particle model has a regular data structure for a perfect work load balance and a small amount of communication associated with the space-charge solver. It scales well on both multiple processor central processing unit computers and multiple graphic processing unit computers [36].
In this paper, we assumed rectangular perfect conducting pipe boundary conditions and used the spectral method to obtain the space-charge potential. For general shape boundary conditions or open boundary conditions, the space-charge potential on the computational grid can be computed using other numerical methods such as the finite element or the Green's function method. The symplectic PIC model proposed here would still be valid if one uses the resultant potential on the grid.