Macroparticle Collisionality in PIC solver

This paper presents a new numerical technique for the calculation of the Coulomb interaction between nearby charged particles. The efficiency of this method is demonstrated by distributing charged particles in a three dimensional sphere for an amount of particles, Npar , ranging from 1 × 105 to 1.6 × 106. Using the presented method, the computation of the mean field and Coulomb near field interaction can be performed separately.


Introduction
The trajectory of charged particles in an accelerator is determined by the Coulomb law, with i, the test particle, j the index of summation running from 1 to N par , q j the charge of the j-th particle and ϵ 0 the permittivity in vacuum.Equation (1) can be decomposed in two components: 1) The mean field deriving by smoothing of the particle distribution and, 2) the intense field created by the closeness of the particles, which we call the "near field".
The main focus of this paper is on the development of a numerical method for the efficient determination of near field interactions.Methods to reduce the complexity of determining the near Coulomb interaction were introduced in 1987 by L. Greengard and V. Rohklin [1].In that paper, the authors present a method called the Fast Multipole Method (FMM) applied to an N-body system, which handles all pairwise interactions.In this paper, a numerical method is presented that drastically reduces the computation time and load of the near field interactions.

The Mean Field Computation
The distribution of N par charged particles is embedded in a box of the sizes X W , Y W and Z W .The electric field is determined by solving Poisson's equation, where ρ(⃗ x) is the charge distribution and ϕ(⃗ x) the electric potential.The solution of equation ( 2) is The integrand is the product of the charge distribution ρ(⃗ x ′ ) by the distance dependent function G(⃗ x, ⃗ x ′ ) = 1/|⃗ x − ⃗ x ′ |, which is known as the Green's function.Various numerical methods for solving Poisson's equation to determine the electric field are presented in books and publications [2,3].A method that has proven to be the best for computational speed is the spectral method.The numerical solution of the convolution of G and ρ, namely equation (3), is obtained using the Fast Fourier Transformation.Schematically, ρ(⃗ x) and G(⃗ x) are each transformed into their reciprocal space by FFT.By multiplying and performing an inverse Fast Fourier Transformation of the obtained function, the electric potential is computed.However, before the electric field is calculated, the boundary conditions such as the Dirichlet condition must be imposed.Setting the boundary to a fixed value is achieved by introducing relaxation methods (Jacobi, Gauss-Seidel, ADI, SOR, Multigrid Method), which are iterative methods that update the solution until it converges to a steady state, see Ref. [2,4,5].The FFT itself, requires a discretisation of the box into equal sub-volumes.To use the FFT efficiently, the box is split into 2 m − 1 boxes, to have 2 m gridpoints, where m is an integer and can be different in each direction.The grid size is computed with Once the box volume is discretised, the next step is to deposit the charges of the particles to the nearest grid points using special functions [2].A common charge deposition scheme is the Next GridPoint Method (NGP) where the particle charges are fully reassigned to the nearest gridpoint.Another deposition method is the Cloud In Cell (CIC).In this case, the charges are reassigned to the nearest gridpoints using a linear function.The charge deposition scheme illustrated in figure 1 shows how particle charges in the NGP are treated.Note that the spectral method intrinsically smooths the particle distribution as the charges are deposited on a grid.Hence, the field from near particles (near field) is removed.

Determination of the Coulomb Near Field
As mentioned in the introduction, to determine the near field, techniques such as the FMM are used, from which our method was inspired.Just as for the calculation of the mean field, the box containing the beam is divided into a new grid, the Close Encounter Grid (CEG).The new grid has a size of ∆s C,min , equal in all directions, computed according to where E mf,max is the maximum mean field reached.The number of cells for the three directions is N CEG,x , N CEG,y , N CEG,z and computed with The symbol ⌈•⌉ is a function that rounds floating numbers to the next larger integer.

Algorithm
The new algorithm is based on the use of a four-dimensional array C the "near particle"-array.The cells in the CEG are identified by (j, k, l) of C with Each particle is identified by an index p that runs from 1 to N par and is located at (x p ,y p ,z p ).All the CEG are identified in an array of integers while the particle positions are floating numbers.The CEG cell where the particle p is located is found from In the next step, the macro-particles are stored in the array C j,k,l,m according to the following procedure.The particle with index p identifies the CEG cell of indexes (j p , k p , l p ).The indexes (j p ,k p ,l p ) are assigned to the corresponding indexes (j,k,l) of C j,k,l,m .The number of macroparticles N c located in the CEG cell identified by (j,k,l) is N c = C j,k,l,0 .When a new particle is identified to this cell, then we update C j,k,l,0 to N c + 1, and assign the particle identification number p to C j,k,l,(Nc+1) .The fourth index m range from zero to the number of macro-particles the CEG contains.The zeroth cell, for m = 0, in the array C j,k,l,m contains the number of particles at grid position (j, k, l).For m > 0 the particle p index is stored.
For example, consider 3 particles with the indexes p = 2, 4, 7 located in the CEG cell of indexes (56, 12, 43).At the beginning of the loading procedure C j,k,l,0 = 0.
Next, the steps for loading the 3 particles are listed in table 1.The step 1 is the initial situation.The step two is the loading of the particle p = 2.As we load the first particle, we set C j,k,l,0 = 1 and C j,k,l,1 = 2.The third step repeat the same procedure, and the total number of macro-particles in the cell is updated to C j,k,l,0 = 2 with C j,k,l,2 = 4, and so on for the fourth step.This procedure is carried out until all particles have been assigned to a C j,k,l,m .This method reduces the number of particle-particle Coulomb field evaluations in order to assess the near field.Without our methods, the computational load is ∝ N 2 par , applying this procedure reduces the computational load to ∝ N α par , with α ≤ 1.In addition to the advantages previously described, this algorithm can be easily parallelisable as it treats each macro-particle independent from all others.

Numerical Tests
To demonstrate the efficiency of the presented algorithm the following simulation tests have been conducted.A rectangular box of size 0.1 m × 0.1 m × 0.1 m is defined and contain the test particles.The particles are spherically distributed with a radius of 0.01 m.The test particle distribution is taken uniform (FORTRAN random number generator) and Gaussian (algorithm see Ref. [6]).The number of particles ranged from 10 5 to 1.6 × 10 6 .The number of grids in the mean field computation is 128 3 for the uniform and Gaussian distribution as well as for the different N par .In the near field computation the grid number varies according to equation (5).The grid numbers generated for the near field is listed in table 2.   2 shows the computational load to derive the Coulomb field for all particles using Coulomb's law, equation (1), the computational complexity in this case is O(N 2 par ).Instead figure 3 shows the performance of the algorithm we have developed.Figure 4 shows the mean field computational cost of deriving the electric field with a 3D spectral Poisson solver.By comparing figure 3 with figure 4 it is evident that the near field computation of the near field with our method accounts for less than 5% of the computation load from the 3D spectral solver.The orange curve shows the simulation results for the case of a test beam with a normal distribution, while the blue curve is for a test with a particle uniform macro-particle distribution.

Outlook
This work has presented a new method to include in a standard 3D Poisson solver, that usually compute only the mean field, the near field component.Our tests confirm that the performance of the near field method does not affect the overall performance of the 3D solver.Hence, we foresee a future application of "3D solver" + "near field algorithm" into multi-particle simulations for a more comprehensive modeling the effect of intra-beam particle scattering in storage rings.A comparison with the prediction of an analytic models like the Bjorken-Mtingwa or Conde-Martini will allow to validate the method.The challenges in the tracking is to conserve the energy-momentum.By applying the regularization method where the denominator of equation ( 1) is extended by a distance term.This prevents a singularity so that the electric field does not go to infinity.In addition the integration steps will have to be adaptive to avoid macro-particles creating artificial effects.

Figure 1 .
Figure 1.2D Nearest GridPoint (NGP) charge deposition scheme.The deposition is taken place to the closest gridpoint according to the arrows.

Figure 3 .
Figure 3. Reduction of CPU time computing the near field of a uniform (blue) & Gaussian (orange) distribution of all test macroparticles.

Figure 4 .
Figure 4. CPU time computing the mean field with the spectral method for a uniform (blue) & Gaussian (orange) distribution of all test macro-particles.

Table 1 .
Example for j=56, k=12, l=43 and the Values Stored in the C j,k,l,m .

Table 2 .
Number of Grids for the Near Field Computation for the Uniform and Gaussian Distribution.