A symplectic particle-in-cell model for space-charge beam dynamics simulation
AA symplectic particle-in-cell model for space-charge beam dynamics simulation
Ji Qiang ∗ Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA
Space-charge effects play an important role in high intensity particle accelerators and were studiedusing 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 modeland a conventional nonsymplectic PIC model for long-term space-charge simulation of a coastingbeam. 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 shapesin the benchmark examples. Quantitatively, the symplectic PIC and the symplectic gridless particlemodels agree with each other very well, while the nonsymplectic PIC model yields different emittancegrowth value. Using finer step size, the emittance growth from the nonsymplectic PIC convergestowards that from the symplectic PIC model.
I. INTRODUCTION
The nonlinear space-charge effects from particle interactions inside a charged particle beam play an important rolein high-intensity accelerators. A variety of models have been developed to simulate the space-charge effects [1–13].Among those models, the particle-in-cell (PIC) model is probably the most widely used model in the acceleratorcommunity. The particle-in-cell model is an efficient model to handle the space-charge effects self-consistently. It usesa computational grid to obtain the charge density distribution from a finite number of macroparticles and solves thePoisson equation on the grid at each time step. The computational cost is linearly proportional to the number ofmacroparticles, which makes the simulation fast for many applications. However, those grid based, momentum con-served, PIC codes do not satisfy the symplectic condition of classic multi-particle dynamics. Violating the symplecticcondition in multi-particle tracking might not be an issue in a single pass system such as a linear accelerator. Butin a circular accelerator, violating the symplectic condition may result in undesired numerical errors in the long-termtracking simulation. Recently, a symplectic gridless particle model was proposed for self-consistent space-charge sim-ulation [14]. This model is easy to implement and has a good scalability on massive parallel computers. However, thecomputational cost of this model is proportional to the number of macroparticles times the number of modes used tosolve the Poisson equation. When the number of modes is large, this model becomes quite computationally expensiveand has to resort to parallel computers.In this paper, we propose a symplectic PIC model for self-consistent space-charge long-term simulation. This modelhas the advantages of the computational efficiency of the PIC method and the symplectic property needed for long-term dynamics tracking. Multi-symplectic particle-in-cell model was proposed to study Vlasov-Maxwell system inplasmas using a variational method [15–18]. A symplectic gridless spectral electrostatic model in a periodic plasma wasalso proposed following the variational method [19]. To the best of our knowledge, at present, there is no symplecticPIC model that was proposed and shown applications for the quasi-static Vlasov-Poisson system in space-charge beamdynamics study. In this paper, we derived a symplectic PIC model directly from the multi-particle Hamiltonian. Wethen carried out long-term space-charge simulations and compared this model with the symplectic gridless particlemodel 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 multi-particle trackingmodel including the space-charge effect in Section II; We present the nonsymplectic PIC model in Section III; Wecompare the symplectic PIC model, the symplectic gridless particle model, and the nonsymplectic PIC model withtwo space-charge simulation examples in Section IV; and draw conclusions in Section V. ∗ Electronic address: [email protected] a r X i v : . [ phy s i c s . acc - ph ] J a n II. SYMPLECTIC MULTI-PARTICLE TRACKING MODELS
In the accelerator beam dynamics simulation, for a multi-particle system with N p charged particles subject to botha space-charge self field and an external field, an approximate Hamiltonian of the system can be written as [20–22]: H = N p (cid:88) i =1 p i / N p (cid:88) i =1 N p (cid:88) j =1 qϕ ( r i , r j ) + N p (cid:88) i =1 qψ ( r i ) (1)where H ( r , r , · · · , r N p , p , p , · · · , p N p ; s ) denotes the Hamiltonian of the system using distance s as an independentvariable, ϕ is related to the space-charge interaction potential between the charged particles i and j (subject toappropriate 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 ) the normalizedcanonical 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 oflight in vacuum. The equations governing the motion of individual particle i follow the Hamilton’s equations as: d r i ds = ∂H∂ p i (2) d p i ds = − ∂H∂ r i (3)Let ζ denote a 6N-vector of coordinates, the above Hamilton’s equation can be rewritten as: dζds = − [ H, ζ ] (4)where [ , ] is the Poisson bracket. A formal solution for above equation after a single step τ can be written as: ζ ( τ ) = exp( − τ (: H :)) ζ (0) (5)Here, we have defined a differential operator : H : as : H : g = [ H, g ], for arbitrary function g . For a Hamiltonianthat can be written as a sum of two terms H = H + H , an approximate solution to above formal solution can bewritten as [23] ζ ( τ ) = exp( − τ (: H : + : H :)) ζ (0)= exp( − τ : H :) exp( − τ : H :) exp( − τ : H :) ζ (0) + O ( τ ) (6)Let exp( − τ : H :) define a transfer map M and exp( − τ : H :) a transfer map M , for a single step, the abovesplitting results in a second order numerical integrator for the original Hamilton’s equation as: ζ ( τ ) = M ( τ ) ζ (0)= M ( τ / M ( τ ) M ( τ / ζ (0) + O ( τ ) (7)The above numerical integrator can be extended to 4 th order accuracy and arbitrary even-order accuracy followingYoshida’s approach [24]. This numerical integrator Eq. 7 will be symplectic if both the transfer map M and thetransfer map M are symplectic. A transfer map M i is symplectic if and only if the Jacobian matrix M i of thetransfer map M i satisfies the following condition: M Ti JM i = J (8)where J denotes the 6 N × N matrix given by: J = (cid:18) I − I (cid:19) (9)and I is the 3 N × N identity matrix.For the Hamiltonian in Eq. 1, we can choose H as: H = N p (cid:88) i =1 p i / N p (cid:88) i =1 qψ ( r i ) (10)This corresponds to the Hamiltonian of a group of charged particles inside an external field without mutual interactionamong themselves. A single-charged particle magnetic optics method can be used to find the symplectic transfer map M for this Hamiltonian with the external fields from most accelerator beam line elements [21, 22, 25].We can choose H as: H = 12 N p (cid:88) i =1 N p (cid:88) j =1 qϕ ( r i , r j ) (11)which includes the space-charge effect and is only a function of position. For the space-charge Hamiltonian H ( r ),the single step transfer map M can be written as: r i ( τ ) = r i (0) (12) p i ( τ ) = p i (0) − ∂H ( r ) ∂ r i τ (13)The Jacobi matrix of the above transfer map M is M = (cid:18) I L I (cid:19) (14)where L is a 3 N × N matrix. For M to satisfy the symplectic condition Eq. 8, the matrix L needs to be a symmetricmatrix, i.e. L = L T (15)Given the fact that L ij = ∂ p i ( τ ) /∂ r j = − ∂ H ( r ) ∂ r i ∂ r j τ , the matrix L will be symmetric as long as it is analyticallycalculated from the function H . This is also called jolt-factorization in nonlinear single particle beam dynamicsstudy [26]. If both the transfer map M and the transfer map M are symplectic, the numerical integrator Eq. 7 formulti-particle tracking will be symplectic.For a coasting beam, the Hamiltonian H can be written as [21]: H = K N p (cid:88) i =1 N p (cid:88) j =1 ϕ ( r i , r j ) (16)where K = qI/ (2 π(cid:15) p v γ ) is the generalized perveance, I is the beam current, (cid:15) is the permittivity of vacuum, p is the momentum of the reference particle, v is the speed of the reference particle, γ is the relativistic factor of thereference particle, and ϕ is the space charge Coulomb interaction potential. In this Hamiltonian, the effects of thedirect Coulomb electric potential and the longitudinal vector potential are combined together. The electric Coulombpotential in the Hamiltonian H can be obtained from the solution of the Poisson equation. In the following, weassume that the coasting beam is inside a rectangular perfectly conducting pipe. In this case, the two-dimensionalPoisson’s equation can be written as: ∂ φ∂x + ∂ φ∂y = − πρ (17)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: φ ( x = 0 , y ) = 0 (18) φ ( x = a, y ) = 0 (19) φ ( x, y = 0) = 0 (20) φ ( x, y = b ) = 0 (21)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 approximatedusing two sine functions as [27–31]: ρ ( x, y ) = N l (cid:88) l =1 N m (cid:88) m =1 ρ lm sin( α l x ) sin( β m y ) (22) φ ( x, y ) = N l (cid:88) l =1 N m (cid:88) m =1 φ lm sin( α l x ) sin( β m y ) (23)where ρ lm = 4 ab (cid:90) a (cid:90) b ρ ( x, y ) sin( α l x ) sin( β m y ) dxdy (24) φ lm = 4 ab (cid:90) a (cid:90) b φ ( x, y ) sin( α l x ) sin( β m y ) dxdy (25)where α l = lπ/a and β m = mπ/b . The above approximation follows the numerical spectral Galerkin method since eachbasis function satisfies the boundary conditions on the wall [27–29]. For a smooth function, this spectral approximationhas 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 PoissonEq. 17 and making use of the orthonormal condition of the sine functions, we obtain φ lm = 4 πρ lm γ lm (26)where γ lm = α l + β m .In the simulation, the particle distribution function ρ ( x, y ) can be represented as: ρ ( x, y ) = 1 N p N p (cid:88) j =1 S ( x − x j ) S ( y − y j ) (27)where S ( x ) is the shape function (also called deposition function in the PIC model). Using the above equation andEq. 24 and Eq. 26, we obtain: φ lm = 4 πγ lm ab N p N p (cid:88) j =1 (cid:90) a (cid:90) b S ( x − x j ) S ( y − y j ) sin( α l x ) sin( β m y ) dxdy (28)and the electric potential as: φ ( x, y ) = 4 π ab N p N p (cid:88) j =1 N l (cid:88) l =1 N m (cid:88) m =1 γ lm sin( α l x ) sin( β m y ) (cid:90) a (cid:90) b S ( x − x j ) S ( y − y j ) sin( α l x ) sin( β m y ) dxdy (29)The electric potential at a particle i location can be obtained from the potential on grid as: φ ( x i , y i ) = (cid:90) a (cid:90) b φ ( x, y ) S ( x − x i ) S ( y − y i ) dxdy (30)where the interpolation function from the grid to the particle location is assumed to be the same as the depositionfunction.From the above electric potential, the interaction potential ϕ between particles i and j can be written as: ϕ ( x i , y i , x j , y j ) = 4 π ab N p N l (cid:88) l =1 N m (cid:88) m =1 γ lm (cid:90) a (cid:90) b S ( x − x j ) S ( y − y j ) sin( α l x ) sin( β m y ) dxdy (cid:90) a (cid:90) b S ( x − x i ) S ( y − y i ) sin( α l x ) sin( β m y ) dxdy (31)Now, the space-charge Hamiltonian H can be written as: H = 4 π K ab N p N p (cid:88) i =1 N p (cid:88) j =1 N l (cid:88) l =1 N m (cid:88) m =1 γ lm (cid:90) a (cid:90) b S ( x − x j ) S ( y − y j ) sin( α l x ) sin( β m y ) dxdy (cid:90) a (cid:90) b S ( x − x i ) S ( y − y i ) sin( α l x ) sin( β m y ) dxdy (32) A. Symplectic Gridless Particle Model
In the symplectic gridless particle space-charge model, the shape function is assumed to be a Dirac delta functionand the particle distribution function ρ ( x, y ) can be represented as: ρ ( x, y ) = 1 N p N p (cid:88) j =1 δ ( x − x j ) δ ( y − y j ) (33)Now, the space-charge Hamiltonian H can be written as: H = 4 π K ab N p N p (cid:88) i =1 N p (cid:88) j =1 N l (cid:88) l =1 N m (cid:88) m =1 γ lm sin( α l x j ) sin( β m y j ) sin( α l x i ) sin( β m y i ) (34)The one-step symplectic transfer map M of the particle i with this Hamiltonian is given as: p xi ( τ ) = p xi (0) − τ πK ab N p N p (cid:88) j =1 N l (cid:88) l =1 N m (cid:88) m =1 α l γ lm sin( α l x j ) sin( β m y j ) cos( α l x i ) sin( β m y i ) p yi ( τ ) = p yi (0) − τ πK ab N p N p (cid:88) j =1 N l (cid:88) l =1 N m (cid:88) m =1 β m γ lm sin( α l x j ) sin( β m y j ) sin( α l x i ) cos( β m y i ) (35)Here, both p xi and p yi are normalized by the reference particle momentum p . B. Symplectic Particle-In-Cell Model
If the deposition/interpolation shape function is not a delta function, the one-step symplectic transfer map M ofthe particle i with this space-charge Hamiltonian H is given as: p xi ( τ ) = p xi (0) − τ πK ab N p N p (cid:88) j =1 N l (cid:88) l =1 N m (cid:88) m =1 γ lm (cid:90) a (cid:90) b S ( x − x j ) S ( y − y j ) sin( α l x ) sin( β m y ) dxdy (cid:90) a (cid:90) b ∂S ( x − x i ) ∂x i S ( y − y i ) sin( α l x ) sin( β m y ) dxdyp yi ( τ ) = p yi (0) − τ πK ab N p N p (cid:88) j =1 N l (cid:88) l =1 N m (cid:88) m =1 γ lm (cid:90) a (cid:90) b S ( x − x j ) S ( y − y j ) sin( α l x ) sin( β m y ) dxdy (cid:90) a (cid:90) b S ( x − x i ) ∂S ( y − y i ) ∂y i sin( α l x ) sin( β m y ) dxdy (36)where both p xi and p yi are normalized by the reference particle momentum p , and S (cid:48) ( x ) is the first derivative of theshape 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: p xi ( τ ) = p xi (0) − τ πK ab N p N p (cid:88) j =1 N l (cid:88) l =1 N m (cid:88) m =1 γ lm (cid:88) I (cid:48) (cid:88) J (cid:48) S ( x I (cid:48) − x j ) S ( y J (cid:48) − y j ) sin( α l x I (cid:48) ) sin( β m y J (cid:48) ) (cid:88) I (cid:88) J ∂S ( x I − x i ) ∂x i S ( y J − y i ) sin( α l x I ) sin( β m y J ) p yi ( τ ) = p yi (0) − τ πK ab N p N p (cid:88) j =1 N l (cid:88) l =1 N m (cid:88) m =1 γ lm (cid:88) I (cid:48) (cid:88) J (cid:48) S ( x I (cid:48) − x j ) S ( y J (cid:48) − y j ) sin( α l x I (cid:48) ) sin( β m y J (cid:48) ) (cid:88) I (cid:88) J S ( x I − x i ) ∂S ( y I − y i ) ∂y i sin( α l x I ) sin( β m y J ) (37)where the integers I , J , I (cid:48) , and J (cid:48) denote the two dimensional computational grid index, and the summations withrespect to those indices are limited to the range of a few local grid points depending on the specific deposition function.If one defines the charge density ρ ( x (cid:48) I , y (cid:48) J ) on the grid as: ρ ( x I (cid:48) , y J (cid:48) ) = 1 N p N p (cid:88) j =1 S ( x I (cid:48) − x j ) S ( y J (cid:48) − y j ) , (38)the space-charge map can be rewritten as: p xi ( τ ) = p xi (0) − τ πK (cid:88) I (cid:88) J ∂S ( x I − x i ) ∂x i S ( y J − y i )[ 4 ab N l (cid:88) l =1 N m (cid:88) m =1 γ lm (cid:88) I (cid:48) (cid:88) J (cid:48) ρ ( x I (cid:48) , y J (cid:48) ) sin( α l x I (cid:48) ) sin( β m y J (cid:48) ) sin( α l x I ) sin( β m y J )] p yi ( τ ) = p yi (0) − τ πK (cid:88) I (cid:88) J S ( x I − x i ) ∂S ( y I − y i ) ∂y i [ 4 ab N l (cid:88) l =1 N m (cid:88) m =1 γ lm (cid:88) I (cid:48) (cid:88) J (cid:48) ρ ( x I (cid:48) , y J (cid:48) ) sin( α l x I (cid:48) ) sin( β m y J (cid:48) ) sin( α l x I ) sin( β m y J )] (39)It turns out that the expression inside the bracket represents the solution of potential on grid using the charge densityon grid, which can be written as: φ ( x I , y J ) = 4 ab N l (cid:88) l =1 N m (cid:88) m =1 γ lm (cid:88) I (cid:48) (cid:88) J (cid:48) ρ ( x I (cid:48) , y J (cid:48) ) sin( α l x I (cid:48) ) sin( β m y J (cid:48) ) sin( α l x I ) sin( β m y J ) (40)Then the above space-charge map can be rewritten in more concise format as: p xi ( τ ) = p xi (0) − τ πK (cid:88) I (cid:88) J ∂S ( x I − x i ) ∂x i S ( y J − y i ) φ ( x I , y J ) p yi ( τ ) = p yi (0) − τ πK (cid:88) I (cid:88) J S ( x I − x i ) ∂S ( y J − y i ) ∂y i φ ( x I , y J ) (41)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]: S ( x I − x i ) = − ( x i − x I h ) , | x i − x I | ≤ h/ ( − | x i − x I | h ) , h/ < | x i − x I | ≤ / h ∂S ( x I − x i ) ∂x i = − x i − x I h ) /h, | x i − x I | ≤ h/ − + ( x i − x I ) h ) /h, h/ < | x i − x I | ≤ / h, x i > x I ( + ( x i − x I ) h ) /h, h/ < | x i − x I | ≤ / h, x i ≤ x I y dimension.Using the symplectic transfer map M for the external field Hamiltonian H from a magnetic optics code and thetransfer map M for space-charge Hamiltonian H , one obtains a symplecti PIC model including the self-consistentspace-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: d r i ds = p i (44) d p i ds = q ( E i /v − a z × B i ) (45)where r = ( x, y ), p = ( v x /v , v y /v ), and a z is a unit vector in the longitudinal z direction.The above equations can be solved using a 2 nd order leap-frog algorithm, which is widely used in many dynamicssimulations. Here, for each step, the positions of a particle i are advanced a half step by: r ( τ / i = r (0) i + 12 τ p i (0) (46)Then the momenta of the particle are advanced a full step including the forces from the given external fields andthe self-consistent space-charge fields. In the PIC model, the space-charge fields are obtained by solving the Poissonequation on a computational grid and then interpolating back to the particle location. In this study, we employ theabove spectral method to attain the space-charge fields as: E x ( x I , y J ) = − ab N l (cid:88) l =1 N m (cid:88) m =1 α l γ lm (cid:88) I (cid:48) (cid:88) J (cid:48) ρ ( x I (cid:48) , y J (cid:48) ) sin( α l x I (cid:48) ) sin( β m y J (cid:48) ) cos( α l x I ) sin( β m y J ) (47) E y ( x I , y J ) = − ab N l (cid:88) l =1 N m (cid:88) m =1 β m γ lm (cid:88) I (cid:48) (cid:88) J (cid:48) ρ ( x I (cid:48) , y J (cid:48) ) sin( α l x I (cid:48) ) sin( β m y J (cid:48) ) sin( α l x I ) cos( β m y J ) (48)where the charge density ρ ( x I (cid:48) , y J (cid:48) ) on the grid is obtained following the charge deposition Eq. 38. The particlemomenta can then be advanced one step using: p xi ( τ ) = p xi (0) + τ ( qE extx v − qB exty ) + τ πK (cid:88) I (cid:88) J S ( x I − x i ) S ( y J − y i ) E x ( x I , y J ) p yi ( τ ) = p yi (0) + τ ( qE exty v + qB extx ) + τ πK (cid:88) I (cid:88) J S ( x I − x i ) S ( y J − y i ) E y ( x I , y J ) (49)Here, both the direct electric Coulomb field and the magnetic field from the longitudinal vector potential are includedin the space-charge force term of the above equation. Then, the positions of the particle are advanced another halfstep using the new momenta as: r ( τ ) i = r ( τ / i + 12 τ p i ( τ ) (50)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 andthe 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 thousandsof lattice periods using the above three models.In the first example, we simulated a 1 GeV coasting proton beam transport through a linear periodic quadrupolefocusing and defocusing (FODO) channel inside a rectangular perfectly conducting pipe. A schematic plot of thelattice is shown in Fig. 1. It consists of a 0 . . FIG. 1: Schematic plot of a FODO lattice. magnet with a single period. The total length of the period is 1 meter. The zero current phase advance throughone lattice period is 85 degrees. The current of the beam is 450 A and the depressed phase advance is 42 degrees.Such a high current drives the beam across the 4 th order collective instability [34]. The initial transverse normalizedemittance of the proton beam is 1 µm with a 4D Gaussian distribution.Figure 2 shows the four dimensional emittance growth ( (cid:15) x (cid:15) x (cid:15) y (cid:15) y − ,
000 lattice periodsfrom the symplectic gridless particle model, from the symplectic PIC model, and from the nonsymplectic spectralPIC model. These simulations used about 50 ,
000 macroparticles and 15 ×
15 modes and 257 ×
257 grid points in
FIG. 2: Four dimensional emittance growth evolution from the symplectic gridless particle model (red), the symplectic PICmodel (green), and the nonsymplectic spectral PIC (blue). the spectral Poisson solver. The perfectly conducting pipe has a square shape with an aperture size of 10 mm. It isseen that the symplectic PIC model and the symplectic gridless particle model agrees with each other very well. Thenonsymplectic spectral PIC model yields significantly smaller emittance growth than those from the two symplecticmethods, which might result from the numerical damping effects in the nonsymplectic integrator. The fast emittancegrowth within the first 20 ,
000 periods is caused by the space-charge driven 4 th order collective instability. The slowemittance 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 thesethree models. It is seen that the distributions from all three models show similar phase space shapes after the strongcollective resonance. However, the core density distributions show somewhat different structures. The two symplecticmodels show similar structures while the nonsymplectic spectral PIC model shows a denser core. This is also seen inthe 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 emittancegrowth evolution from the symplectic PIC model and those from the nonsymplectic PIC model with the same nominalstep size, from the nonsymplectic PIC model with one-half of the nominal step size, and from the nonsymplectic PICmodel with one-quarter of the nominal step size. It is seen that as the step size decreases, the emittance growth fromthe 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. Aschematic plot of the lattice is shown in Fig. 6. Each period of nonlinear lattice consists of 10 periods of 1 m linearFODO lattice and a sextupole magnet to emulate a simple storage ring. The zero current horizontal and vertical tunesof the ring are 2 .
417 and 2 .
417 respectively. The integrated strength of the sextupole element is 10 T/m/m. Theproton beam initial normalized transverse emittances are still 1 um in both directions with a 4D Gaussian distribution.Figure 7 shows the 4D emittance growth evolution with 10 A, 20 A, and 30 A currents from the symplectic PICmodel, the symplectic gridless particle model, and the nonsymplectic PIC model. The space-charge tune shift with10A beam current is 0 . . . rd order resonance, there is little emittance growth. With 20 A current, the proton beam moves near the 3 rd resonance, some particles fall into the resonance and results in significant increase of the emittance growth. With 30A current, a lot of particles fall into the resonance and results in large emittance growth. It is seen that from theabove plot, all three models predict this current dependent behavior qualitatively. The two symplectic models showvery similar emittance growth evolution for all three currents while the nonsymplectic PIC model shows much lessemittance growth compared with the other two models. This is probably due to the numerical damping effects in thenonsymplectic PIC model.Figure 8 shows the final transverse phase space distributions after 40 ,
000 turns. All three models show similarphase 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 modelsyield very close density profiles with non-Gaussian halo tails. The nonsymplectic PIC model yields denser core butless halo near the tail density distribution.
FIG. 3: Final transverse phase space from the symplectic PIC model (top), from the symplectic gridless particle model (middle),and from the nonsymplectic spectral PIC model (bottom).
V. CONCLUSIONS
In this paper, we proposed a symplectic PIC model and compared this model with the symplectic gridless particlemodel and the conventional nonsymplectic PIC model for long-term space-charge simulations of a coasting beam ina linear transport lattice and a nonlinear lattice. Using the same step size and the same number of modes for space-charge solver, all three models qualitatively predict the similar emittance growth evolution due to the space-chargedriven resonance and yield similar final phase space shapes in the benchmark examples. Quantatively, the symplecticPIC and the symplectic gridless particle model agree with each other very well, while the nonsymplectic PIC modelyields much less emittance growth value. Using finer step size, the emittance growth from the nonsymplectic PIC0
FIG. 4: Final horizontal (left) and vertical (right) density profiles from the symplectic gridless particle model (red), thesymplectic PIC model (green), and the nonsymplectic spectral PIC (blue).FIG. 5: Four dimensional emittance growth evolution from the symplectic symplectic PIC model (red), the nonsymplecticspectral PIC with the nominal step size (green), the spectral PIC with half of the nominal step size (blue), and the spectralPIC with one quarter of the nominal step size (pink).FIG. 6: Schematic plot of a one-turn lattice that consists of 10 FODO lattice and one sextupole element. 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 com-putational 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 ofmacroparticles used, the symplectic PIC model is computationally more efficient than the symplectic gridless particlemodel. However, on a massive parallel computer, the scalability of the symplectic PIC model may be limited bythe challenges of work load balance among multiple processors and communication associated with the grided basedspace-charge solver and particle manager [35]. The symplectic gridless particle model has regular data structure forperfect work load balance and small amount of communication associated with the space-charge solver. It scales wellon both multiple processor Central Processing Unit (CPU) computers and multiple Graphic Processing Unit (GPU)computers [36].1 FIG. 7: 4D emittance growth evolution with 10A, 20A, and 30A beam current from the symplectic PIC model (top), thesymplectic gridless particle model (middle), and the nonsymplectic spectral PIC model (bottom).
ACKNOWLEDGEMENTS
Work supported by the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. We would like tothank Drs. A. Friedman, C. Mitchell, R. D. Ryne, and J. Vay for discussions. This research used computer resourcesat the National Energy Research Scientific Computing Center. [1] A. Friedman, D. P. Grote, and I. Haber, Phys. Fluids B 4 , 2203 (1992).[2] H. Takeda and J. H. Billen, Recent developments of the accelerator design code PARMILA, in Proc. XIX InternationalLinac Conference, Chicago, August 1998, p. 156.[3] S. Machida and M. Ikegami, in AIP Conf. Proc 448, p.73 (1998).[4] F. W. Jones and H. O. Schoenauer, Proc of PAC1999, p. 2933, 1999.[5] J. Qiang, R. D. Ryne, S. Habib, V. Decyk, J. Comput. Phys. , 434, 2000.[6] P. N. Ostroumov and K. W. Shepard. Phys. Rev. ST. Accel. Beams 11, 030101 (2001).[7] R. Duperrier, Phys. Rev. ST Accel. Beams , 124201, 2000.[8] J. D. Galambos, S. Danilov, D. Jeon, J. A. Holmes, and D. K. Olsen, F. Neri and M. Plum, Phys. Rev. ST Accel. Beams , 034201, (2000).[9] H. Qin, R. C. Davidson, W. W. Lee, and R. Kolesnikov, Nucl. Instr. Meth. in Phys. Res. A 464, 477 (2001).[10] G. Franchetti, I. Hofmann, M. Giovannozzi, M. Martini, and E. Metral, Phys. Rev. ST Accel. Beams , 124201, (2003).[11] J. Qiang, S. Lidia, R. D. Ryne, and C. Limborg-Deprey, Phys. Rev. ST Accel. Beams , 044204, 2006.[12] J. Amundson, P. Spentzouris, J. Qiang and R. Ryne, J. Comp. Phys. vol. 211, 229 (2006).[13] http://amas.web.psi.ch/docs/opal/opal user guide.pdf. FIG. 8: Final transverse phase of the 30A beam from the from the symplectic PIC model (top), the symplectic gridless particlemodel (middle), and the nonsymplectic spectral PIC model (bottom).[14] J. Qiang, Phys. Rev. Accel. Beams , 014203, (2017).[15] J. Xiao, J. Liu, H. Qin, and Z. Yu, Phys. Plasmas , 102517 (2013)[16] E. G. Evstatiev and B. A. Shadwick, J. Comput. Phys. , p.376 (2013).[17] B. A. Shadwick, A. B. Stamm, and E. G. Evstatiev, Phys. Plasmas , 055708 (2014).[18] H. Qin et al., Nuclear Fusion , 014001, (2016).[19] S. Webb, Plasma Phys. Control. Fusion , 034007, (2016).[20] E. Forest, “Beam Dynamics: A New Attitude and Framework,” The Physics and Technology of Particle and PhotonBeams, Vol. 8 (Harwood Academic Publishers, Amsterdam, 1998).[21] R. D. Ryne, “Computational Methods in Accelerator Physics,” US Particle Accelerator class note, 2012.[22] A. J. Dragt, “Lie Methods for Nonlinear Dynamics with Applications to Accelerator Physics,” 2016.[23] E. Forest and R. D. Ruth, Physica D , p. 105, 1990. FIG. 9: Final horizontal (left) and vertical (right) density profiles from the symplectic gridless particle model (red), thesymplectic PIC model (green), and the nonsymplectic spectral PIC (blue).[24] H. Yoshida, Phys. Lett. A , p. 262, 1990.[25] http://mad.web.cern.ch/mad/.[26] E. Forest, J. Phys. A: Mathe Gen. p. 5321, 2006.[27] D. Gottlieb and S. A. Orszag, Numerical Analysis of Spectral Methods: Theory and Applications, Society for Industrialand Applied Mathematics, 1977.[28] B. Fornberg, A Practical Guide to Pseudospectral Methods, Cambridge University Press, 1998.[29] J. Boyd, Chebyshev and Fourier Spectral Methods, Dover Publications, Inc. 2000.[30] J. Qiang and R. D. Ryne, Comp. Phys. Comm. , p. 18, 2001.[31] J. Qiang, Comp. Phys. Comm. , p. 122, 2016.[32] R.W. Hockney, J.W. Eastwood, Computer Simulation Using Particles, Adam Hilger, New York, 1988.[33] C. K. Birdsall and A. B. Langdon,Plasma Physics Via Computer Simulation,Taylor and Francis, New York, 2005.[34] I. Hofmann, L. J. Laslett, L. Smith, and I. Haber, Part. Accel. 13, 145 (1983).[35] J. Qiang and X. Li, Comp. Phys. Comm.181