Comparison of Simulation Methods of Ion-Atomic Collisions in PIC-MC

The main ion-atomic collision treatment methods based on Monte-Carlo simulation are considered and discussed. We have proposed an efficient scheme for simulation of time between collisions taking into account cross-section dependence on ion velocity and random generation of ion velocities and scattering angles after collisions. The developed algorithm of simulation of interval between collisions takes into account the change of relative velocity of ion-atom pair as well as the change of cross-section of collision and atomic concentration. At the same time, unlike the widely used “null-collision” method, both the probability of collision and change of particles’ state which determines this probability are taken into consideration for each particle independently in time. The simulation results according to the techniques proposed are found to be close to the theoretical values of ion drift velocities. It is revealed that the “null-collision” method results in exceeding of drift velocity in strong and intermediate fields. At the same time the proposed method of accumulation of probability under the same conditions gives values close to theoretical ones. In weak fields calculated values of drift velocity in both methods exceed theoretical values to some small extent.


Introduction
The simulation of plasma by Particle-in-Cell Monte-Carlo (PIC-MC) method has received widespread application in calculations of ion current on probe and dust particle and plasma simulation of glow and high-frequency discharges [1]- [6].The important stage in PIC-MC method is simulation of ion collisions with neutral atoms.It includes interval simulation between collisions, simulation of velocity and direction of ion motion after collision.
V. Sysun et al. 1234 The drift time of an ion up to the next collision has probabilistic nature: probability of collisions for the time dt doesn't depend on the previous way and proportional to dt : ( ) ( ) ( ) , n -concentration of atoms, ( ) u t -relative velocity ion-atom, ( ) u σ -cross-section of collision.
Let us assume by this ( ) ′ -probability of absence of collision for the time t .Then the probability of collision for dt after non-collision time t , determining decrease of P′ as t increases to dt , is ( ) ( ) ( ) ( ) ( ) Equation ( 2) is probability of collision since the time t is: ( ) ( ) The following methods of simulation of the interval between collisions are applied.

Constant Time between Collisions [7]-[9]
Equation ( 3) is approximation of constant time between collisions ( ) ( ) This approximation is correct at inverse dependence of cross-section on velocity.
In the supposition based on Equation (3) ( ) Simulation of the random value "τ " is done with the help of uniformly distributed at [0,1] random value "rand".Equation (4) shows the expression for random value "rand" ( ) ( ) In Equation (4) ( ) is replaced with rand , as it is the random value as well.

Constant Path Length [5] [6] [9] [10]
Equation ( 5) is the expression for path length in this case Let us replace a variable in Equation ( 1) , then taking into account Equation (5) we will get Equation ( 6): ( ) For each ion path length is simulated after each collision.Equations (7) are the expressions for ion's path length ( ) ( ) ( ) V. Sysun et al. 1235 The achievement of individual path length by each ion is checked by summing up its ways on each time step The approximation of constant path length takes into account the change of ion velocity, however the crosssection of collisions is considered to be constant here that can be accepted with some approximation at resonant charge exchange.In [9] a problem of the method considered has been revealed.If the method of constant path length was used, the decrease of average ion energy was observed because of more frequent collisions of fast ions transferring energy to atoms and absence of direct Maxwellian process among charged particles in PIC.

Simulation of Collision Probability on Each Time Step [1] [3] [11]
For the time interval t ∆ relative velocity and cross-section are considered to be constant.Equation ( 8) is probability of collision for the time t ∆ : While each ion is simulated, if ⋅ ∆ the collision occurs.Values u and σ cor- respond to the current time point t .This method demands significant increase of computing time that makes it hardly applicable to large ensembles of particles.

"Null-Collision" Algorithm
This algorithm is considered in details in [2] and actively used [2] [4] [12].At first, the maximum value of product ( ) max u σ ⋅ on the whole range of all possible relative ion-atom velocities is determined (or appointed).
Upon the value ( ) max u σ ⋅ the constant conditional collision probability for the time t ∆ is calculated.Equation ( 9) is constant conditional collision probability for the time t ∆ ( ) ( ) Then, in case of total number of ions N on each time step, the number of ion is simulated P 0 N times.For this ion rand is simulated again and the type of collision or the absence of it is defined.Equation ( 10) is the condition of event that k-type of collision occurs: collision doesn't occur.It compensates for the randomness of choice ( ) max uσ .
Treatment methods of collision probability are combined in the "null-collision" technique.At first, the number of the ion is determined in the ensemble, and then the type of collision or its absence is defined through the change of ion state in time.The number of arithmetic operations ( ) times is smaller than in the method of determination of collision probability on each time step.
In [12] "null-collision" method is applied to define the collision probability on each time step.At the collision doesn't occur, and it occurs in the opposite case.
There is a variety of revisions of the "null-collision" algorithm used for simulation of electron-atom collisions when cross-sections are strongly dependent on electrons' energy.Thus, in the method of constant time between collisions the time interval is determined by maximum possible value of ( ) ( ) ( ) then we will have ( ) ( ) for density of probability.Equation ( 12) is collision probability for the time τ : Equations ( 13) are the expressions if individual ion is simulated: Summing up by steps in time lasts up to the reaching of equality in Equation ( 13).In such a way, this algorithm takes into account the change of u as well as the change of σ and n , having approximately the same number of arithmetic operations.At the same time, unlike in the "null-collision" method, both the probability of collision and change of particles' state which determines this probability are taken into consideration for each particle independently in time.According to the algorithm features it is natural to call this technique "the method of accumulation of probability".

Simulation of Energy and Direction of Ion Motion after the Collision
Direction of ion motion after the collision is characterized by taking into account deviation from original direction θ and azimuth angle ϕ .Angle ϕ is equiprobable at symmetrical atoms 2π rand ϕ =
At elastic scattering of ion on atom in these works scattering goes only forward.Equation ( 15) is the formula for θ in this case: Equation ( 16) is energy after scattering: In works [11] [14] at elastic collisions the law of ion-atom interaction is stated.Impact parameter and relative velocity are simulated by randomizer and then depending on the minimal radius of approximation the angles and ion energies are calculated after collision.This method makes simulation process significantly more difficult.At the same time Equations ( 15), ( 16) don't take into account the transfer of energy from atoms to ions that is the most significant at weak fields.
Let us consider Monte-Carlo algorithm of simulation of elastic ion-atom collisions on the model of elastic spheres.In this model in Figure 1 normal and tangential components of ions' velocity after collision are: , where m and M are weights of ion and atom, correspondingly.
Further we consider ion motion in their own gas m M = , then 1 . Angles α and β -original angles  17) is simulation formula for the azimuth angle: For the density of probability we have ( ) , where , it corresponds to Equations ( 15), ( 16) taking into account that 1 rand − is the same random value like rand .However, in weak fields an ion can get significant additional velocity from an atom in collisions.Equation ( 19) is ions' velocity after collision in this case ( ) ( ) ( ) Let us consider atom 2 υ simulation.Equation ( 21) is the absolute velocity Maxwell distribution function: ( ) ( ) Equation ( 20) is the expression for simulation after introducing the new parameter The integral can be replaced by an analytic expression, having approximation with accuracy up to 3%.Equation (23) is the approximation for t υ ( ) ( ) The results of calculations presented in Table 1 reveal the approximation (23) to be close to the exact values.
Equation ( 24) is density of probability β angles and azimuth angle 2 ϕ : It determines . Angle θ determines the ion deviation from its direction before collision.To determine the ion current in the given direction (along the field) it is necessary to know the angle to this direction (angle θ ′ in Figure 2).

Let us have angle 0
θ to the axis х up to the collision, θ and ϕ are angles of scattering.Equation ( 25) is the expression for θ ′ : It is here that we need angle ϕ according to Equation (17).At resonant charge exchange angle θ ′ is simulated immediately.Equations ( 26) are the formulae in this case

Simulation Experiment
To compare "null-collision" method and method of accumulation of probability the simulation experiment of ion motion in constant electric field E in approximation to elastic spheres was carried out.The motion with the charge exchange on atoms and elastic interaction was considered separately.The number of ions was taken as equal to 10 5 .Original ion velocities had Maxwell distribution with the temperature of atoms by Equations ( 23) and (24).The interval of time t ∆ at integrating motion equations is chosen depending on the probability of collision for t ∆ according to Equation ( 9) within ( ) Value n σ was taken as constant.Collision process was simulated according to Equations ( 17 Equation ( 30) is the approximation used in intermediate fields: The results of simulation are given in Table 2 and Table 3 ( )   It is necessary to mention that computing time in "null-collision" method (in the same conditions) turned out to be 2 -3 times less than in method of accumulation of probability.However, in case of large ( ) P t ∆ probability of collision increases more slowly than the velocity achieved for the time between collisions.That results in exceeding of drift velocity in the "null-collision" method.So at P = 0.2; 0.5; 0.8 this exceeding corresponded to approximately 4%; 16%; 40%, respectively.
Method of accumulation of probability in strong and intermediate fields gives values close to theoretical ones.Accuracy is limited only by fluctuations of average velocity increasing when the total number of ions decreases and their temperature goes up.In weak fields calculated values of drift velocity in both methods exceed theoretical values to some extent (up to 5%).It can be caused by the fact that when Equations (28) are obtained integrated square velocity is accepted as an average relative velocity:

Conclusion
The given effective models of calculating time between collisions by the method of accumulation of probability and ion angle velocities' simulation after collision on the basis of solid spheres result in ion drift velocities close to theoretical ones and can be applied in the process of simulation of ion motion in heterogeneous plasma with non-constant concentration of atoms and ions and dependence of cross-section on velocity.

Figure 2 .
Figure 2. Angles of ion deviation: θ -to the direction up to the col- lision, θ ′ -to the direction of the electric field.
the range of(2 -20)  values, achieved between collisions for the average time.At the same time change ( ) max u σ in the given range practically didn't influence the average drift velocity.

Improvement of Methods and Simulation Experiment 3.1. Improvement of Methods of Constant Path Lengths and Time between Collisions-The Method of Accumulation of Probability
Individual time interval is simulated for each electron ( ) 0 ln rand τ τ = − and when this value is achieved then the event of collision should be treated by generation a new rand.If ⋅ < then the collision takes place, otherwise the absence of collision is assumed.V. Sysun et al. 1236 3.

Table 1 .
Comparison of approximated t υ  and exact t υ values. 1239

Table 2 .
Value of drift velocity d υ′ at resonant charge exchange.

Table 3 .
Value of drift velocity d υ′ at elastic collisions.