Fast 3D Poisson solvers in elliptical conducting pipe for space-charge simulation

Space-charge effects play an important role in high intensity accelerators. These effects can be studied self-consistently by solving the Poisson equation with the dynamically evolved charge density distribution subject to appropriate boundary conditions. In this paper, two computationally efficient methods are proposed to solve the Poisson equation inside an elliptical perfectly conducting pipe. One method uses a spectral method and the other uses a spectral finite difference method. The former method has a high accuracy and the latter one has a computational complexity of O ð N log ð N ÞÞ , where N is the total number of unknowns. These methods implemented in a beam dynamics tracking code enable the fast simulation of space-charge effects in an accelerator with an elliptical conducting pipe.


I. INTRODUCTION
The nonlinear space-charge effects play an important role in high intensity accelerators by driving beam instability, causing beam emittance growth, halo formation, and even particle loss. These effects can be studied selfconsistently using a particle-in-cell (PIC) method [1][2][3][4][5][6][7][8][9][10][11]. In the PIC method, at each time step, macroparticles that represent the phase space distribution of charged particle beam in the accelerator are transformed from the laboratory frame into the moving beam frame following the relativistic Lorentz transformation and are deposited onto a twodimensional (2D) or three-dimensional (3D) computational grid to attain the charge density distribution in spatial domain. Next, the Poisson equation is solved in the beam frame, the space-charge electric fields are calculated and transformed back to the laboratory frame following the field Lorentz transformation. Then, the space-charge fields are interpolated to the individual macroparticle location from the computational grid. These space-charge fields together with the external accelerating and focusing fields are used to advance the macroparticle momenta through the time step. The updated momenta are used to advance macroparticle positions in the spatial domain. This process is repeated for many time steps until the stopping condition is reached.
To calculate the space-charge fields, one needs to solve the Poisson equation for a given charge density distribution.
A key issue in the PIC simulation is to solve the Poisson equation subject to appropriate boundary conditions efficiently at each time step. In some accelerators such as the Proton Synchrotron at CERN, the conducting pipe that contains a train of charged particle bunches has an elliptical transverse shape. To study the space-charge effects inside those accelerators, an efficient three-dimensional (3D) Poisson solver subject to the transverse elliptical perfectly conducting wall and longitudinal periodic or open boundary conditions is needed.
In previous studies, a number of efficient Poisson solvers have been developed in the community [12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29]. Those solvers will not handle transverse elliptical conducting pipe and longitudinally open or periodic boundary conditions. Finite difference methods were proposed to solve the Poisson equation subject to irregular boundary geometry [30][31][32]. These methods use finite-difference approximation to the differential operator in the Poisson equation. Near the boundary, the curved geometry is approximated by stepwise grid and the boundary conditions are applied on the grid points. However, for an elliptical shape boundary, those approximations can be avoided. By using an elliptical coordinate system as shown in Fig. 1 [33], the curved elliptical boundary can be transformed into a regular rectangular geometry.
In this paper, two new methods are proposed to solve the 3D Poisson equation inside an elliptical conducting pipe with longitudinal periodic or open boundary conditions. One method is based on the transverse pseudospectral and longitudinal spectral Galerkin method. The other method is based on the transverse Galerkin spectral finite-difference method and longitudinal spectral Galerkin method. The former method has the advantage of exponential convergence of a spectral method [34] and can be used in simulations where high precision is needed. The latter method has a second-order accuracy and a computational complexity scaling of OðN logðNÞÞ, and can be used in simulations where quick results are needed.
In the beam frame, the three-dimensional normalized Poisson equation in Cartesian coordinates can be written as: where, ϕ denotes the normalized electric potential inside the beam, ρ denotes the normalized charge density distribution of the beam, and x, y, and z denote the horizontal, vertical, and longitudinal normalized dimensionless coordinates, respectively. The boundary conditions for the electric potential on the elliptical perfectly conducting pipe are ϕðx; y; zÞj ∂Ω ¼ 0 ð2Þ where ðx; yÞ ∈ ∂Ω and satisfies the following elliptical equation: where 2a coshðbÞ denotes the length of the major axis and 2a sinhðbÞ the length of the minor axis. A direct solution to the above Poisson equation subject to the transverse elliptical boundary condition has to handle the irregular curved geometric shape of the boundary. However, this irregular shape of the boundary can be transformed into a regular shape if one uses an elliptical coordinate system as shown in Fig. 1. In the scaled dimensionless elliptical coordinate system ðr; θÞ, the Cartesian coordinates are given by: x ¼ a coshðbrÞ cosðθÞ ð 4Þ y ¼ a sinhðbrÞ sinðθÞ: The Poisson Eq. (1) in the scaled elliptical coordinate system can be rewritten as: and the transverse boundary conditions become: ϕðr; θ þ 2π; zÞ ¼ ϕðr; θ; zÞ: The longitudinal boundary condition can be either periodic if one assumes a train of identical bunches or open if one neglects the interactions among multiple bunches. The organization of this paper is as follows: After the Introduction, the pseudospectral method for a 2D coasting beam inside an elliptical perfectly conducting pipe is discussed in Sec. II; the 3D spectral solver for a bunched beam subject to transverse elliptical boundary and longitudinal periodic or open boundary conditions is presented in Sec. III; the 3D spectral finite-difference solver for the bunched beam is presented in Sec. IV, and conclusions are drawn in Sec. V.

II. THE PSEUDOSPECTRAL SOLVER FOR A 2D COASTING BEAM
A 2D coasting beam with no longitudinal dependency inside a transverse elliptical conducting pipe is considered first. In this case, the Poisson equation in elliptical coordinates is reduced into: and the transverse boundary conditions become: Equation (9) can be solved efficiently using a pseudospectral method.
Using the Chebyshev collocation points along the radial r dimension and the uniform points along the θ dimension, that is, where L þ 1 and M are number of grid points in each dimension respectively, the electric potential can be approximated by: where ϕ lm is the solution of the electric potential at the grid point ðl; mÞ, and p l ðrÞ and f m ðθÞ are corresponding cardinal functions in the r and the θ dimension respectively, and given by [35,36]: Given the above cardinal functions, the differentiation matrix D r in radial dimension is given by: where and the differentiation matrix D θ in the θ dimension is given by: Using these differentiation matrices, the derivative of the electric potential with respect to the radial dimension at grid point r l can be written as: Similarly the derivative with respect to the θ dimension at grid point θ m can be written as: Assuming that Eq. (14) is correct at each grid point and substituting the above representation for the second-order derivative into the original Poisson Eq. (9), the electric potential ϕ lm at each grid point ðl; mÞ can be obtained from the solution of the following algebraic equations: In the above equation, the Dirichlet boundary condition at the pipe wall, i.e., ϕ 0m ¼ ϕ Lm ¼ 0, is used. The solutions of this group of algebraic equations yield electric potentials inside the domain ð−1; 1Þ ⊗ ½0; 2πÞ due to the fact that the Chebyshev collocation points map the original radial domain from [0, 1] to ½−1; 1. These solutions are not all independent of each other. The point ð−r; 2π − θÞ and the point ðr; θÞ in the elliptical coordinate system represent the same point in the Cartesian coordinate. Moreover, there is no need to find the solutions inside the domain ½−1; 0Þ ⊗ ½0; 2πÞ. Using the same location condition, i.e., ϕ jm ¼ ϕ ðL−jÞðM−mÞ for j > L=2, Eq. (23) can be rewritten as: The number of unknowns in the above equation is L=2 × M instead of the original ðL − 1Þ × M in the Eq. (23). This group of linear algebraic equations can be solved using a lower-upper (LU) decomposition method [30]. The above solver was tested using a nearly round pipe with a ¼ 0.036619 and b ¼ 4, which results in a pipe radius close to 1. The normalized charge density distribution is assumed as: The analytical solution of the normalized electric potential in this case is: Figure 2 shows the electric potential distribution inside the elliptical pipe using the above pseudospectral method and the second-order finite difference method of Ref. [37] with 32 × 32 grid points. Both methods give qualitatively similar azimuthal symmetric electric potential distribution. Figure 3 shows the relative errors of the electric potential solution from the above pseudospectral method and from the second-order finite difference method. It is seen that the pseudospectral method yields an order of magnitude smaller error than the second-order finite difference method.
As an application, using the above pseudospectral solver, the electric potential distribution inside a perfectly conducting elliptical pipe was computed for a normalized Gaussian density distribution: Figure 4 shows the electric potential distribution inside the pipe with the horizontal to vertical aspect ratio 1, 2, 4, and 10. Here, 64 × 64 grid points were used and the minor axis length was fixed as 12. As the aspect ratio increases, the electric field distribution becomes more asymmetric and extends further in the vertical direction.

III. THE 3D PSEUDOSPECTRAL SOLVER FOR A BUNCHED BEAM
In most accelerators, the charged particles are longitudinally bunched using rf focusing. In this case, one needs to solve the 3D Poisson equation to study the space-charge effects. First, a case with periodic boundary condition in the longitudinal direction is considered for a train of identical bunches. With this boundary condition, the electric potential satisfies: where P is the longitudinal period length. Given the longitudinal periodic condition, the normalized electric potential and the charge density function can be approximated using Fourier mode expansion: ρðr; θ; zÞ ¼ X For each Fourier mode n, one can use the pseudospectral method along r and θ dimensions as discussed in the preceding section, the solution of the above equation can be written as: . Both x and y coordinates are dimensionless and the normalized electric potential is also dimensionless.
The above equation can be solved using a LU decomposition method for each mode n.
As an application, the electric potential function was calculated for a 3D normalized Gaussian density function: ρðx; y; zÞ ¼ exp − 1 2 Figure 5 shows the normalized electric potential distribution in the middle x − y, y − z, and x − z planes inside an elliptic pipe with horizontal major axis length of 24, vertical minor axis length of 6, and longitudinal periodic length P ¼ 30 using 64 × 64 transverse grid points and 64 longitudinal modes. The electric potential has a longitudinally periodic distribution that decays quickly toward the pipe transverse boundary. If the separation of the bunches inside a bunch train is large and the interactions among bunches are negligible, the longitudinal boundary condition of the electric potential for a single bunch in this dimension can be regarded as an open boundary condition, that is, This open boundary condition can be approximated using a closed Dirichlet boundary condition as long as the longitudinal domain is large enough. The longitudinal boundary condition for the electric potential can be written as: and α n ¼ nπ=ð2TÞ. Substituting the above expansions into the Poisson Eq. (6) and making use of the orthonormal condition of the sine function yield: 1 a 2 ½sinh 2 ðbrÞ þ sin 2 ðθÞ For each mode ϕ n , using the pseudospectral method along r and θ dimensions in the preceding section, the solution of the above equation can be written as: ϕ ljn D 2 θ;jm − a 2 ½sinh 2 ðbr l Þ þ sin 2 ðθ m Þα 2 n ϕ lmn ¼ −a 2 ½sinh 2 ðbr l Þ þ sin 2 ðθ m Þρ n ðbr l ; θ m Þ: ð43Þ The above equation can be solved using the same LU decomposition method for each mode n as the preceding case.
As an application, the electric potential was calculated for the above 3D Gaussian density distribution inside the elliptical pipe with a longitudinal domain size of 100, i.e., T ¼ 50. Figure 6 shows the electric potential distribution in the middle x − y, y − z, and x − z planes inside the elliptic pipe using 64 × 64 transverse grid points and 64 longitudinal modes. The electric potential decays quickly in both longitudinal and transverse dimensions.

IV. THE 3D SPECTRAL FINITE-DIFFERENCE SOLVER FOR A BUNCHED BEAM
The pseudospectral method in the transverse plane derived in the previous section results in a group of algebraic equations that need to be solved using the LU decomposition method. In this section, another method is presented to approximate the transverse Laplacian operator and results in a sparse matrix that can be solved efficiently using an iterative method.
The Poisson Eq. (6) in the elliptical coordinate can be rewritten as Given the periodic condition in the θ dimension, the normalized electric potential and the charge density can be approximated using the Fourier mode expansion in that dimension as: Substituting the above expansions into the Poisson Eq. (6) and making use of the orthonormal condition of the Fourier mode function yields: For the longitudinal periodic boundary condition with a periodic length P, using the above Fourier mode expansion in the longitudinal dimension for each mode m, the above equation is reduced into: ϕ mn ðrÞ þ a 2 4 λ 2 n ϕ ðmþ2Þn ðrÞ þ a 2 4 λ 2 n ϕ ðm−2Þn ðrÞ ¼ρ mn ðrÞ: ð51Þ The above ordinary differential equations can be solved using a finite difference method for each mode m and n.
Here, a second-order finite difference approximation to the above differential operator at radial grid point r l is adopted as This results in the following linear algebraic equations: This group of linear algebraic equations can be solved using a block Gaussian-Siedel iteration method. For the given ϕ lðmþ2Þn , ϕ lðm−2Þn and source termρ lmn , the tridiagonal matrix can be quickly solved for ϕ lmn with a computational cost scaling as OðLÞ. The obtained ϕ lmn is used to approximate the ϕ lðmþ2Þn , ϕ lðm−2Þn for the next iteration. The spectral finite difference method was also applied to the same bunched beam with the Gaussian distribution in Eq. 35. The electric potential solutions in three middle planes are given in Fig. 7. It is seen that these solutions agree with the solutions from the pseudospectral method very well.
For the longitudinal open boundary condition, the same sine function expansion in the longitudinal dimension as in the preceding section and the second-order finite difference approximation in the radial dimension are used. The resultant linear algebraic equation is given as: This group of equations is solved using the same iteration method. The same bunched beam density distribution inside the elliptical pipe as the one in the preceding section was also used to test this method. The electric potentials in three middle plane are given in Fig. 8. These solutions also agree with those from the pseudospectral method very well.

V. CONCLUSIONS
In this paper, two efficient computational methods were presented to solve the Poisson equation inside an elliptical perfectly conducting pipe subject to longitudinal periodic or open boundary. One method uses a transverse pseudospectral and longitudinal spectral Galerkin method. The other uses a transverse Galerkin spectral finite-difference method and the same longitudinal spectral method. The former method has the advantage of exponential convergence of the spectral method with smooth density distribution and can be used in the simulation where high precision of the solution is needed. The latter method has a second-order accuracy and a computational cost scaling as O(Nlog(N)). This method can be used in simulations where the fast return results are needed. The above methods assume a smooth charge density function in the Poisson equation. The numerical errors (noise) in the charge density function obtained from the macroparticles in a PIC simulation will spoil the exponential convergence of the pseudospectral method and has less impact on the finite difference method. However, the numerical noise in the charge density function can be substantially mitigated by using a global spectral basis function approximation (e.g., sine function) to the density distribution function as proposed in some recent studies [38][39][40]. In a future study, we plan to implement those efficient methods into the parallel beam dynamics code suite [4,8], IMPACT, to study the space-charge effects in high intensity accelerators.