Null collision Monte Carlo simulation model for particle-in-cell method

The article discusses the original null collision approach for modeling collisions in plasma for particle-in-cell method. A distinctive feature of the method is the simulation of particle collisions without the use of trigonometric functions and calculation of pairwise collisions (as is done in classical papers), which significantly increases the speed of the developed algorithms. The developed algorithms are used in other software systems for modeling particle motion in a magnetic field. Comparison of the results of system’s operation in the presence of collisions and in their absence is made.


Introduction
Simulating plasma collisions using the Monte Carlo method is a well-known problem. Many articles have been published on this topic with differing approaches to its solution.
For example, in the classical paper [1], pairwise collisions of particles are simulated based on principles of classical mechanics, which allows conserving the total energy and momentum of the entire system. Unfortunately, for the particle-in-cell method, when the number of tracked particles and the number of iterations are large enough, this approach is excessively laborious, since it requires searching and sorting over the entire array of particles [2]. In addition, the calculation of pairwise collisions is rather slow because of the need to calculate the scattering angles and new velocities of the particles.
Various authors in their papers [3][4] have proposed a number of ideas that allow avoiding the calculation of pairwise or group collisions using the so-called "Null collision" approach. In this approach, there is no calculation of the direct interaction (collision) of two or more particles, but the speed and trajectory of individual particle change. The parameters of these changes (collision frequency, deviation magnitude, new velocity) are determined from other parameters, for example, density, average thermal ion velocity, etc.
Despite the fact that this approach can significantly reduce the program execution time, it nevertheless has its drawbacks. One of the main disadvantages in the case of a three-dimensional problem is the fact that it is required to determine the new coordinates of the particle's motion vector immediately after changing its trajectory, which cannot be done without the use of trigonometric functions (whatever method is used).

 
The calculation of such a number of trinonometric functions on a computer is a relatively slow procedure, which ultimately affects the final speed of program execution [5].
We have proposed and programmatically implemented an approach based on "Null collision", which does not use trigonometric functions at all.

Trigonometry-free null collision Monte Carlo simulation model
The proposed method is based on the fact that the interaction of any two particles in a cell can be represented as the interaction of this particle and the "average" particle in this cell. This eliminates the need to search and sort by an array of particles.
The interaction of particles is represented as the scattering of every particle on a stationary "average" particle (for this, a transition to the corresponding coordinate system is done).
To determine the mean free path of a particle l i , a formula is used from [8] where T i is thermal energy, N k is density in particle's cell, L i is Coulomb logarithm with ionic components.
The essence of the method is easiest to demonstrate in the figure. The place of "impact", and, as a consequence, the angle of deflection of the particle from its initial trajectory is chosen randomly.
The maximum angle of deflection of a particle is taken as 30% (in this case, the principles of classical mechanics are used, specifically, the solution of the problem when one ball hits a stationary other with the same mass. It is easy to show that the angle of deflection of a flying ball cannot exceed π/6).
The algorithm consists of the following steps (see fig. 2): 1. Around the end of the particle velocity vector (green vector), a ball is constructed, the radius of which is chosen so that the deflection angle (between green and red vectors) does not exceed π/6. 2. A cube with sides 2R is built around the ball.  Figure 2. New motion vector of a particle 3. Then we choose 3 random uniformly distributed variables a, b and c (from -1 to 1), which are responsible for a random point inside this cube (normalized by length). We need them to be inside the ball, not inside the cube, so we check the condition a 2 + b 2 + c 2 ≤ 1. If it fails (a point in the cube, but outside the ball), we re-generate a, b and c until the condition is met.
The probability of choosing 3 suitable variables on the first try is approximately 52%, and the average number of steps is approximately 1.9. Given that choosing a random number is a very fast procedure, regeneration practically does not slow down the program.
4. We build a new vector (in the figure in yellow) with added variables a, b and c.
At the end of the procedure for calculating the new parameters of motion for all particles, the energy and momentum of the entire system is normalized. As you can see, trigonometric functions are not used in the algorithm. However, this approach also has drawbacks, since it requires normalizing the parameters of the system at each step of the collision calculation.
Next, we briefly describe another collisional algorithm in order to compare the results of the proposed null collision algorithm.

Algorithm for collisional diffusion calculating
In this section, we will briefly talk about one of the collisional algorithms (then we will compare the results of this algorithm with the algorithm described above).
Let a monoenergetic beam of particles of type a move through a stationary plasma. Due to the interaction with plasma particles, the average beam velocity < v α > decreases in accordance with the equation ( In a plasma with a Maxwellian particle velocity distribution, the frequencies v Here e α and m α are charge and mass of α-type particle, Λ αβ is Coulomb logarithm, u β = 2T β /m β is thermal velocity of β-type particles. Due to angular scattering, the scatter in the transverse component (across the average beam velocity) and the longitudinal components of the velocities of the beam particles increases  Figure 3. Velocity distribution of particles as a result of the null collision algorithm monotonically (4)- (5): We will not dive into the specifics of the realization of this algorithm on a computer, we just note that it was implemented in software and then compared with the above-described null collision algorithm.

Results
First, let's compare the operation of the null collision and collisional algorithms, and in this case we do not pay attention to the speed of operation, but to the obtained distributions.

Comparison of the distributions of null collision and collisional algorithms
Both algorithms were implemented in computer software. The simulation showed that both approaches lead to the Maxwellian distribution of particles (see figures 3 and 4).
It can be seen that the distributions are similar, so the null collision algorithm can be used as a faster alternative.

Comparison of speed with other null collision algorithms
In the classical paper [4], the scattering of a particle is realized in the form of a deviation of its trajectory. In this case, the length of the vector is conserved, which does not require normalization.
This approach was implemented as a test program in order to compare the performance of these two algorithms.   Simulation showed that for a system containing 5000 particles, the time gain reaches 10% due to the absence of trigonometric functions (even taking into account the normalization procedure) compared to the approach described in [4]. However, it should be noted that the result may be different when calculating on specialized computers, orientated on calculation of trinonometric functions.

Application of the procedure in modeling plasma confinement systems
The developed null collision method was used in modelling plasma confinement systems to estimate the motion of particles in a magnetic field taking ion-ion collisions into account [9]. The most significant difference is the distribution of background ions during the operation of such systems, which can be seen in the figures 5 and 6.  Figure 6. Velocity distribution of particles as a result of collisional diffusion algorithm As can be seen in the figures, ion-ion collisions affect the results obtained; therefore, the development of high-speed null collision methods are very important for this kind of research.

Conclusion
In this work, a high-speed algorithm was proposed for simulating ion-ion collisions in plasma for the particle-in-cell method (for one type of ions). The main advantage of the developed algorithms is a very high speed of work due to the absence of the need to calculate trigonometric functions. The developed algorithms were used in plasma confinement modeling systems. The results showed that taking into account ion-ion collisions significantly affects many parameters in the system, such as the distribution of background ions and ion velocities at the cavity boundaries.