Collision routine

Rather than analyzing the individual trajectories of any particle, a collision can generally occur for each pair of particles which are within the same sub-cell.

The collision rate within a cell is calculated as follows: The product of the relative velocity and the corresponding cross section $v_\textrm{r}\cdot\sigma$ multiplied with the time step $\delta t$ spans a cylindrical volume. The ratio of this volume and the cell volume is the collision probability $p$:

$$p = \dfrac{v_\textrm{r}\,\sigma(E)\,\delta t}{V_\textrm{cell}},\qquad E = \frac{1}{2}m_\textrm{r}\,v_\textrm{r}^2.$$

Thus, if we assume that the collision cross section $\sigma$ and all relative velocities $v_\textrm{r}$ are constant, the number $N_\textrm{coll}$ of expected collisions in one cell between two species A and B can be simply calculated by:

$$N_\textrm{coll} = \left\{ \begin{array}{r l} \dfrac{N_\textrm{A}\,N_\textrm{B}\,N_\textrm{Real}\,v_\textrm{r}\,\sigma\,\delta t}{V_\textrm{Cell}} & \textrm{A}\neq\textrm{B}\ \ \dfrac{N_\textrm{A}\,N_\textrm{B}\,N_\textrm{Real}\,v_\textrm{r}\,\sigma\,\delta t}{2\,V_\textrm{Cell}} & \textrm{A}=\textrm{B} \end{array} \right..$$

The additional parameters in this equation are:

• $N_\textrm{A},\,N_\textrm{B}$ = The number of simulation particles of species A and B within the cell
• $N_\textrm{Real}$ = The scaling factor between real and simulated particle number (which is assumed to be the same for species A and B, here)

NTC algorithm

In reality, all relative velocities and corresponding cross sections are different. In this case, the so called "NTC-Method"¹ is applied, which is a variation of the acceptance-reject method of statistics: For each collision type, the maximum product $\textrm{max}(v_\textrm{r}\cdot\sigma)$ is determined from the sampled data of $\sigma(E)$. Then, the number of pair selections $N_\textrm{sel}$ is determined as maximum possible collisions in one cell (i. e. assumed that the product $v_\textrm{r}\cdot\sigma$ is maximal for all possible collision pairs in one cell):

$$N_\textrm{sel} = \left\{ \begin{array}{r l} \dfrac{N_\textrm{A}\,N_\textrm{B}\,N_\textrm{Real}\,\textrm{max}(v_\textrm{r}\cdot\sigma)\,\delta t}{V_\textrm{Cell}} & \textrm{A}\neq\textrm{B}\ \ \dfrac{N_\textrm{A}\,N_\textrm{B}\,N_\textrm{Real}\,\textrm{max}(v_\textrm{r}\cdot\sigma)\,\delta t}{2\,V_\textrm{Cell}} & \textrm{A}=\textrm{B} \end{array} \right..$$

In the follwing, $N_\textrm{sel}$ collision pairs are randomly selected within the sub-cell. The actual relativ velocity and energy is calculated from:

$$v_\textrm{r} = (\vec{v}_\textrm{A}-\vec{v}_\textrm{B})^2,\qquad E = \frac{1}{2}m_\textrm{r}\,v_\textrm{r}^2,\qquad m_\textrm{r} = \dfrac{m_\textrm{A}\,m_\textrm{B}}{m_\textrm{A}+m_\textrm{B}}\,\,\textrm{(reduced mass)}.$$

For a selected collision pair, the collision is actually performed, if the following random experiment holds:

$$R < \dfrac{v_\textrm{r}\,\sigma(E)}{\textrm{max}(v_\textrm{r}\cdot\sigma)},$$

where $R$ is a random number in [0, 1].

General NTC implementation with reaction products

We denote $<\!v_\textrm{r}\,\sigma\!>$ as the mean product of relative velocity and cross section of a collision type within one cell. Further, we have species A with weighting factor $N_\textrm{RA}$, species B with weighting factor $N_\textrm{RB}$ and reaction product C with weighting factor $N_\textrm{RC}$. It shall be assumed that $N_\textrm{RA}>N_\textrm{RB}$. Further, in the current simulation volume, there may be a scaling factor $S_\textrm{R}\neq 1$ for all weighting factors. Under these conditions, the physical number $N_\textrm{col,phys}$ of collisions within a cell becomes:

$$N_\textrm{col,phys} = \left\{ \begin{array}{r l} \dfrac{N_\textrm{A}\,N_\textrm{RA}\,N_\textrm{B}\,N_\textrm{RB}\,S_\textrm{R}^2<\!v_\textrm{r}\,\sigma\!>\delta t}{V_\textrm{Cell}} & \textrm{A}\neq\textrm{B}\ \ \dfrac{N_\textrm{A}\,N_\textrm{RA}\,N_\textrm{B}\,N_\textrm{RB}\,S_\textrm{R}^2<\!v_\textrm{r}\,\sigma\!>\delta t}{2\,V_\textrm{Cell}} & \textrm{A}=\textrm{B}\ \end{array} \right..$$

To determine the number of collisions occurring for the simulation particles, this quantity is to be divided by their respective weighting factor $N_\textrm{RA}$ or $N_\textrm{RB}$ and by $S_\textrm{R}$:

$$N_\textrm{col,A} = \dfrac{N_\textrm{A}\,N_\textrm{B}\,N_\textrm{RB}\,S_\textrm{R}<\!v_\textrm{r}\,\sigma\!>\delta t}{V_\textrm{cell}} = \dfrac{N_\textrm{RB}}{N_\textrm{RA}}\cdot\dfrac{N_\textrm{A}\,N_\textrm{B}\,N_\textrm{RA}\,S_\textrm{R}<\!v_\textrm{r}\,\sigma\!>\delta t}{V_\textrm{cell}}$$ $$N_\textrm{col,B} = \dfrac{N_\textrm{A}\,N_\textrm{B}\,N_\textrm{RA}\,S_\textrm{R}<\!v_\textrm{r}\,\sigma\!>\delta t}{V_\textrm{cell}}.$$

With $N_\textrm{RA}>N_\textrm{RB}$, the quantity $N_\textrm{RA}$ is identical with NR_Preselect in the PICMC code, and $\dfrac{N_\textrm{RB}}{N_\textrm{RA}}$ is denoted as F1_Postselect. F1_Postselect is the probability that the post-collision velocities are applied to the first species with the heigher statistical weigth. The corresponding factor F2_Postselect for the lightweigth species is 1.0. If we now have reaction products of species C the number of reaction products created by the collision becomes:

$$N_\textrm{rp,C} = \dfrac{1}{N_\textrm{RC}\,S_\textrm{R}}\cdot\dfrac{N_\textrm{A}\,N_\textrm{B}\,N_\textrm{RA}\,N_\textrm{RB}\,S_\textrm{R}^2<\!v_\textrm{r}\,\sigma\!>\delta t}{V_\textrm{cell}} = \dfrac{N_\textrm{RB}}{N_\textrm{RC}}\cdot\dfrac{N_\textrm{A}\,N_\textrm{B}\,N_\textrm{RA}\,S_\textrm{R}<\!v_\textrm{r}\,\sigma\!>\delta t}{V_\textrm{cell}}.$$

For $N_\textrm{RB}>N_\textrm{RC}$, multiple reaction products are created for each collision, if $N_\textrm{RB}<N_\textrm{RC}$ only a fraction of $\dfrac{N_\textrm{RB}}{N_\textrm{RC}}$ of the collisions result into a reaction product.

Within the PICMC code this is handled in the following functions:

• void DSMC_Volume::pvm_set_collision(...)
• void DSMC_Volume::collision()

In void DSMC_Volume::pvm_set_collision(...) the above described pre- and postcollision selection factors are defined as follows:

// picmc.cxx line 3376 ff
void DSMC_Volume::pvm_set_collision(int i1, int cflag1,
                                    int i2, int cflag2,
                                    double max_eV,
                                    double alpha, double thresh_eV,
                                    int nrp, int *rp, int ncs, double *cs) {
  /* ... */

  double &NR1 = DSMC_Base::pp[i1].nreal,
         &NR2 = DSMC_Base::pp[i2].nreal;

  plist -> NR_Preselect =  (NR1 > NR2) ? NR1      : NR2;
  plist -> F1_Postselect = (NR1 > NR2) ? NR2/NR1  : 1.0;
  plist -> F2_Postselect = (NR1 > NR2) ? 1.0      : NR1/NR2;

  /* ... */

}

Generally, the post-collision selection factor for species A/B can be written as:

$$F_\textrm{(A,B),postselect} = \dfrac{N_\textrm{RA}\,N_\textrm{RB}}{N_\textrm{R,preselect}}\cdot\dfrac{1}{N_\textrm{R(A,B)}}.$$

The generation of reaction products can be found in DSMC_Volume::collision(...), here the most relevant lines are as follows:

  /* ... */

  minnr = minnr = pp[spec_1].nreal*pp[spec_2].nreal / plist->NR_Preselect;

  /* ... */

 if(pp[plist->rp[i]].nreal < minnr) {
   for(j = 0; j  < (minnr/pp[plist->rp[i]].nreal); j++)
     p3d[plist->rp[i]]->insert(x, y, z, vx, vy, vz, 0.0, icell);

   //statistical handling of reminder if nreal not whole-numbered
   if(random.fetch() < ((minnr/pp[plist->rp[i]].nreal) - j + 1))
      p3d[plist->rp[i]]->insert(x, y, z, vx, vy, vz, 0.0, icell);
 } else {
   if(random.fetch() < (minnr/pp[plist->rp[i]].nreal))
     p3d[plist->rp[i]]->insert(x, y, z, vx, vy, vz, 0.0, icell);
 }

 /* ... */

Handling of frozen species

In some scenarios, species can be frozen, i. e. they are loaded into memory and take part within the collision but do not move themselfes. Also, any post-collision velocity is not applied to frozen species. An example is a transport simulation using the species $\textrm{H}, \textrm{H}_2, \textrm{H}_2^+, \textrm{e}^-$ where $\textrm{H}_2^+, \textrm{e}^-$ are frozen. The purpose of such a simulation is to split the time scales in order to speed up the simulation. In an RF parallel plate plasma with the considered species, the population of $\textrm{H}_2^+$ and $\textrm{e}^-$ becomes stable already after a few microseconds, while the stabilization time of the $\textrm{H}$ population (as reaction product from dissociation of $\textrm{H}_2$) is larger by many orders of magnitude. Given the species A and B, then for a frozen species B the collision numbers as described above become:

$$N_\textrm{col,A} = \dfrac{N_\textrm{A}\,N_\textrm{B}\,N_\textrm{RB}\,S_\textrm{R}<\!v_\textrm{r}\,\sigma\!>\delta t}{V_\textrm{cell}} = \dfrac{N_\textrm{RB}}{N_\textrm{RA}}\cdot\dfrac{N_\textrm{A}\,N_\textrm{B}\,N_\textrm{RA}\,S_\textrm{R}<\!v_\textrm{r}\,\sigma\!>\delta t}{V_\textrm{cell}}$$ $$N_\textrm{col,B}\stackrel{!}{=}0.$$

Now, it turns out that computation time can be saved by using $N_\textrm{RB}$ instead of $N_\textrm{RA}$ as NR_Preselect. In this case, also the post-collision production factor of collision products is modified, the number of collision products becomes:

$$N_\textrm{rp,C} = \dfrac{1}{N_\textrm{RC}\,S_\textrm{R}}\cdot\dfrac{N_\textrm{A}\,N_\textrm{B}\,N_\textrm{RA}\,N_\textrm{RB}\,S_\textrm{R}^2<\!v_\textrm{r}\,\sigma\!>\delta t}{V_\textrm{cell}} = \dfrac{N_\textrm{RA}}{N_\textrm{RC}}\cdot\dfrac{N_\textrm{A}\,N_\textrm{B}\,N_\textrm{RB}\,S_\textrm{R}<\!v_\textrm{r}\,\sigma\!>\delta t}{V_\textrm{cell}},$$

i. e. the collision selection factor is now $\dfrac{N_\textrm{RA}}{N_\textrm{RC}}$ instead of $\dfrac{N_\textrm{RB}}{N_\textrm{RC}}$. In the above mentioned PICMC code this is handled automatically. However this is problematic, if $\dfrac{N_\textrm{RA}}{N_\textrm{RC}}$ is much higher than 1 - in our current example, we have $N_{\textrm{R,H}_2}=5\cdot10^{12}$ and $N_\textrm{R,H}=10^5$. Thus, the creation of reaction products is a very rare incident in the PICMC simulation, but if this incident occurs, $5\cdot10^7$ reaction products are created at once!

As a correction, the pre-selection factor has to be chosen in a way, that one collision product is generated for each collision:

$$N_\textrm{R,Preselect} = \dfrac{N_\textrm{RA}\,N_\textrm{RB}}{N_\textrm{RC}}.$$

The post-collision selection factors are determined as stated above. If multiple reation products occur, the minimum of their weighting factors should be inserted in the denominator. Now, a simple way to speed up the computation is to increase the weighting factor $N_\textrm{RC}$ of the reaction product.

The frozen species are treated in DSMC_Volume::pvm_set_collision(...) as follows.

// picmc.cxx line 3376 ff
void DSMC_Volume::pvm_set_collision(int i1, int cflag1,
                                    int i2, int cflag2,
                                    double max_eV,
                                    double alpha, double thresh_eV,
                                    int nrp, int *rp, int ncs, double *cs) {
  /* ... */

  double &NR1 = DSMC_Base::pp[i1].nreal,
         &NR2 = DSMC_Base::pp[i2].nreal;

  plist -> NR_Preselect =  (NR1 > NR2) ? NR1      : NR2;
  plist -> F1_Postselect = (NR1 > NR2) ? NR2/NR1  : 1.0;
  plist -> F2_Postselect = (NR1 > NR2) ? 1.0      : NR1/NR2;

  double NR_tmp,
         NR_RP = (NR1 > NR2) ? NR1 : NR2;
  int irp;
  for(irp=0; irp<nrp; ++irp) {
    NR_tmp = DSMC_Base::pp[rp[irp]].nreal;
    if (NR_tmp < NR_RP) NR_RP = NR_tmp;
  }

  if(plist->cflag1 & 0x02) {
    plist -> NR_Preselect = NR1 * NR2 / NR_RP;
    plist -> F1_Postselect = 0.0; // frozen
    plist -> F2_Postselect = NR_RP/NR2;
  }

  if(plist->cflag2 & 0x02) {
    plist -> NR_Preselect = NR1 * NR2 / NR_RP;
    plist -> F1_Postselect = NR_RP/NR1;
    plist -> F2_Postselect = 0.0; // frozen
  }
  /* ... */

}

Footnotes

1. N o T ime C ounter Method, see e. g. Bird, G.; Monte Carlo Simulation in an engineering context, In: Progr. Astro. Aero. 74 (1981) 230-55 ↩