Moment-preserving Monte-Carlo Coulomb collision method for particle codes

An expansion of the Boltzmann collision integral (Eqs. 2-3) for small angle Coulomb collisions in an ideal plasma results in the LFP collision operator (see Ch. 11 of [13] and Ch. 13 of Ref. [1]). This operator is of the drag-diffusion type with nonlinear drag and diffusion coefficients. For each particle pair, this process is a diffusive one for the CM scattering angle, $\Theta_{ij}$. The mean square CM scattering angle for a single collision of particle $i\in\alpha$ with particle $j\in\beta$ is computed by averaging $\Theta$ over all possible impact parameters (see Sec. 3.3 of Ref. [12]) and can be expressed as

Here, $b_{\max}=\lambda_{De}$ is the maximum impact parameter for screened Coulomb collisions with $\lambda_{De}$ the plasma Debye length, $\ln\Lambda=\ln\left(2\lambda_{De}/b_{0,ij}\right)$ is the Coulomb logarithm, $b_{0,ij}/2=q_{\alpha}q_{\beta}/\left(4\pi\epsilon_{0}\mu_{\alpha\beta}u_{ij}^{2}\right)$ is the minimum impact parameter with $\mu_{\alpha\beta}=m_{\alpha}m_{\beta}/\left(m_{\alpha}+m_{\beta}\right)$ the reduced mass, and $u_{ij}\equiv|\textbf{v}_{i}-\textbf{v}_{j}|$. The total number of collisions of particle $i$ with particle $j$ in time interval $\Delta t$ is $\pi b^{2}_{\max}n_{j}u_{ij}\Delta t$, where $n_{j}=w_{j}/\Delta V$ is the density associated with particle $j$. The total mean square CM angle for particle $i\in\alpha$ interacting with particle $j\in\beta$ during time interval $\Delta t$ is

Being a diffusive process, the distribution of scattering angles $\Theta_{ij}$ after time $\Delta t$ is Gaussian with variance equal to that given in Eq. 6. The polar scattering angle $\Theta_{ij}$ is thus chosen by sampling from such a distribution. The total change of velocity of particle $i\in\alpha$ is obtained by sequentially performing this process for all $j\in\beta$. To see how this method connects with the LFP equation, the expected value for the change in velocity of particle $i\in\alpha$ is derived below.

The change in $\textbf{u}_{ij}$ parallel to itself after scattering through a small polar angle $\Theta_{ij}$ is $\Delta u_{||}=-2u_{ij}\sin^{2}\left(\Theta_{ij}/2\right)\approx-u_{ij}\Theta_{ij}^{2}/2$. After averaging $\Delta\textbf{u}_{ij}$ over the uniformly distributed azimuthal angle $\phi$, the parallel component is the only non-zero component [13] and we can write

The change in velocity of particle $i\in\alpha$ in the lab frame is related to $\Delta\textbf{u}_{ij}$ by $m_{\alpha}\Delta\textbf{v}_{i}=\mu_{\alpha\beta}\Delta\textbf{u}_{ij}$. The final step in computing the expected value of $\Delta\textbf{v}_{i}$ in time interval $\Delta t$ is to sum Eq. 7 for all $j\in\beta$:

After dividing by $\Delta t$, Eq. 8 is the well-known drag coefficient for the LFP equation (see Eq. 13.28 of Ref. [1]). The tensor diffusion coefficient ($\frac{1}{2}\left<\Delta v_{i}\Delta v_{k}\right>/\Delta t$), can additionally be obtained following a similar process [13]. Note that each term in the sum in Eq. 8 scales linearly with the normalized scattering length $s_{ij}\equiv\frac{1}{2}\left<\Theta^{2}_{ij}\right>_{\text{all }b}$. The same is true for the diffusion coefficient.

The discussion in Sec. 2.1 focused on the transport of a single particle $i\in\alpha$ interacting with species $\beta$. For numerical implementation, we seek to have the correct transport for each particle of each species. The normalized scattering length for particle $i\in\alpha$ interacting with particle $j\in\beta$ via cumulative small angle Coulomb collisions in an ideal plasma, and vice versa, are define as

If the particles are equally weighted ($w_{i}=w_{j}$ for all $i$ and $j$), then $s_{ij}=s_{ji}$ for all $i$ and $j$. One can then loop over all binary pairings (see Fig. 1) and scatter each particle in the pair using the same common scattering angles. This approach obtains the correct scattering physics while also ensuring exact conservation of momentum and energy. However, when the interacting particles have different weights, then $s_{ij}\neq s_{ji}$. The reason is that macro-particles $i$ and $j$ are interacting with a different number of physical particles. One could use independent scattering angles for each particle in the pair. The scattering physics would be correct for each particle, but this approach requires sampling from two different Gaussian distribution functions for each binary pair.

We seek a weighted-particle implementation of the $N{\times}N$ method that only requires sampling from one Gaussian for each binary pair. This is achieved using the following normalized scattering length for the scattering pairs:

where $w_{\max}=\max(w_{i},w_{j})$ is used for the target density. Eq. 10 is the correct value as given in Eq. 9 for the lower weight particle, but it is too large by $w_{\max}/\text{min}(w_{i},w_{j})$ for the larger weight particle. To compensate for this, the velocity vector for the larger weight particle is only updated with a probability $w_{\text{min}}/w_{\max}$. For simplicity, assume that $w_{j}>w_{i}$. Then, after many pairings between particle $i\in\alpha$ and $j\in\beta$, the average value of the normalized scattering length for the larger weight particle is

which is the correct value as given in Eq. 9. This approach is referred to as the rejection method and it is straight-forward to show that, in addition to achieving the correct scattering physics, momentum and energy are also preserved on average [6].

$N{\times}N$ methods are accurate, but the problem is that their computational cost scales with $N^{2}$. The first particle method for Coulomb collisions with computational cost of order $N$ is the seminal work by Takizuka and Abe, which uses $N_{\max}=\max\left(N_{\alpha},N_{\beta}\right)$ pairs per time step. The pairs are formed randomly and are such that each particle from the larger number population scatter exactly once per time step. The particles from the smaller number population scatter $N_{\max}/N_{\text{min}}$ times per time step on average. The normalized scattering length used in this method to the determine the scattering angle for each binary pair is computed using the value for the pair as used for the $N{\times}N$ method given by Eq. 10 multiplied by the number of particles for the population with the least number of particles: $N_{\text{min}}=\text{min}\left(N_{\alpha},N_{\beta}\right)$. This is expressed as

To prove that this method gives the correct normalized scattering length on average for each particle we assume, without loss of generality, that $N_{\beta}>N_{\alpha}$. First, we consider equal weight particles as is considered in TA77. On average after many scatterings, the expected value of the change in velocity of particle $j\in\beta$ (the larger density population) is

which is the correct value with $s_{ji}$ defined in the second expression in Eq. 9. The factor of $1/N_{\alpha}$ in the second expression is the probability that particle $j\in\beta$ is paired with particle $i\in\alpha$. The expected value of the change in velocity of particle $i\in\alpha$ (the smaller density population) is

which is the correct value with $s_{ij}$ defined in the first expression of Eq. 9. The $N_{\beta}/N_{\alpha}$ factor in the second expression is present because each particle from species $\alpha$ is paired $N_{\beta}/N_{\alpha}$ times on average for each scattering event (see Fig. 1). Thus, the scattering physics is correct for all particles on average when using $wN_{\min}/\Delta V$ for the target density. Since the particles are equally weighted here, momentum and energy are conserved exactly for each binary pair.

For weighted particles, Eq. 12 is used, but the velocity of the larger weight particle is only updated with probability $w_{\text{min}}/w_{\max}$, same as it is for the $N{\times}N$ method. The expected value of the change in velocity of particle $j\in\beta$ from the larger number population is

For the last step, we made use of the fact that $P_{j}\left(\alpha_{i},\beta_{j}\right)w_{\max}=w_{i}$ for all cases, where $P_{j}\left(\alpha_{i},\beta_{j}\right)$ is the probability that particle $j\in\beta$ scatters when paired with particle $i\in\alpha$. The expected value of the change in velocity of particle $i\in\alpha$ from the smaller number population is

where again for the last step use is made of the fact that $P_{i}\left(\alpha_{i},\beta_{j}\right)w_{\max}=w_{j}$ for all cases, where $P_{i}\left(\alpha_{i},\beta_{j}\right)$ is the probability that particle $i\in\alpha$ scatters when paired with particle $j\in\beta$.

In summary, the weighted-particle version of TA77 presented here is the same as the equal weight particle version with two changes. First, $w_{\max}$ is used for the density in the normalized scattering length for each binary pair. Second, the velocity of the larger weight particle is updated with probability $w_{\min}/w_{\max}$. The pairing method used here (see Fig. 1), which is the same as in TA77, can also be used for the N97 method. The only difference between these two methods is the formula that takes the normalized scattering length and returns the scattering angle. These methods asymptote to the same variance as $\Delta t$ goes to zero [7].

Everything here is described for inter-species scattering. For intra-species scattering of a species with $N_{\alpha}$ particles in a cell, the pairing is done as described by TA77 and for weighted-particles one uses Eq. 12 with $N_{\text{min}}=N_{\alpha}-1$. When $N_{\alpha}$ is odd, the first three particles are paired in an $N{\times}N$ way and the normalized scattering length is reduced by a factor of two for these pairings to maintain the correct scattering length.

A slightly different method for extending the methods of TA77 and N97 to include weighted particles while preserving the scattering physics on average is described in Ref. [7]. In that work, the normalized scattering length is computed using Eq. 12, but with $N_{\min}$ replaced by $N_{\max}$, and $w_{\max}$ replaced by $w_{\max}/\mathcal{D}_{ij}$, where $\mathcal{D}_{ij}$ is duplicity factor for the particle pair defined by the number of times the particle from the lower number population is paired during a scattering event. Since the average duplicity factor per binary pair is $\left<\mathcal{D}_{ij}\right>=N_{\max}/N_{\min}$, the method described in this work is like that from Ref. [7] when using the average duplicity factor for each particle pair.

Consider two populations, call them $A$ and $B$, of fully ionized carbon 12 atoms ($A=12$, $Z=6$) with $m_{i}=Am_{u}-Zm_{e}$, where $m_{u}$ is the atomic mass unit and $m_{e}$ is the electron mass. The initial density, temperature, and velocity of population $A$ are $n_{A}=10^{19}/\text{cm}^{3}$, $T_{A}=500$ eV, and $U_{A}=+655$ km/s. Those for population $B$ are $n_{B}=10^{20}/\text{cm}^{3}$, $T_{B}=500$ eV, and $U_{A}=0$ km/s. Momentum conservation gives $U_{f}=U_{A}n_{A}/\left(n_{A}+n_{B}\right)=59.55$ for the equilibrium value for the velocity of each species. Energy conservation gives $T_{f}=1.97$ keV for each species. The simulations are performed in quasi 0D using a 1D grid with 180 grid cells. The number of particles in Table 1 below summarizing the numerical setup for this test refer to the number of particles per cell. A fixed value of $\Lambda=10$ is used for the Coulomb logarithm. The particle positions are not advances and there are no fields involved. The simulations are velocity-space only and the results presented below are averaged values on the grid.

The numerical parameters for tests T1a-T1d are summarized in Table. 1. One thing to note is that the number of particles per cell used here is 10 times less than that considered in Ref. [7]. Results for tests T1a (equal particle weights) and tests T1b (equal number of particles) are shown in Fig. 2. The results from test T1a are considered the baseline solution. Results from test T1b shown in this figure, where the ratio of particles weights for the low-density population to the high-density populations is 1:10, are indistinguishable from test T1a. The relaxation of the mean velocity and temperature profiles shown in Fig. 2 can additionally be compared with the same results presented in Fig. 3 of Ref. [7].

The violation in total x-momentum and energy for test T1b with and without the moment-correction algorithm applied are shown in Fig. 3. Although not shown, the velocity and temperature relaxation results shown in Fig. 2 for test T1b are similar independent of using the moment-correction algorithm. The relative error in momentum and energy conservation is between $0.5-1.5$% when not correcting the moments and using unequal weight particles. Momentum and energy conservation are preserved to machine precision with weighted particles when applying the moment-correction method, as seen in the right panel of Fig. 3.

Results from all four test cases described in Table. 1 are shown in Fig. 4. The relaxation of the velocity and temperature profiles for both species match well for all cases. These results show 1) that the weighted scattering algorithm achieves the correct scattering physics for a wide range of weighted-particle scenarios, and 2) the moment-correction method can be applied without significantly altering the scattering physics. It is worth noting here that the original method for doing weighted-particle Coulomb scattering by Nanbu in Ref. [6] does not achieve the correct relaxation for any of the unequal weight particles tests for this problem (see Fig. 6 of Ref. [7]).

Time profiles for the number of binary pairs required to restore energy conservation for test T1d are shown in Fig. 5 for various values of the energy-correction factor $f_{E}$. Results are shown in the middle panel of this figure when using $10\times$ more particles. The number of binary pairs needed to restore energy conservation scales approximately linear with $f_{E}$ and approximately as $\sqrt{N}$. The effects of $f_{E}$ and $N$ on $T_{A}$ from test T1d are shown in Fig. 6. The zoomed in inserts in this figure show the solution converging to the equally weighted particle solution as $f_{E}\rightarrow 0$ and with increasing $N$.

The second test considered is similar to Test 1, but populations A and B are now separated into different species, species $\alpha$ and species $\beta$, and a third population C is added to species $\beta$. Population C is a low-density population with the following initial conditions: $n_{C}=10^{18}/\text{cm}^{3}$, $T_{C}=500$ eV, and $U_{A}=-655$ km/s. The simulation parameters for Test 2 are summarized in Table 2. A variety of different weight combinations for the different populations are considered to rigorously verify the weighted-particle scattering algorithm. Tests T2a-T2e are the same as those considered in Ref. [7], but with $10\times$ less particles, same as for Test 1. Sorting is used for $\alpha-\beta$ and $\beta-\beta$ scattering but is not used for $\alpha-\alpha$ scattering since each of the particles in species $\alpha$ have equal weights.

Time-profiles of the mean velocity and temperature of populations A, B, and C for tests T2a-T2e are shown in Fig. 7. Each of the quantities for each of the populations shown in Fig. 7 are observed to relax to the steady state solutions. Furthermore, the relaxation to the steady state values is observed to agree well for all particle weighting combinations considered. There is a small, but noticeable deviation of the $U_{x}$ profile for population $B$ for test T2e, which is an extreme case where $w_{A}{:}w_{B}{:}w_{C}=1{:}400{:}4$. Furthermore, there are only 100 particles per cell for population $B$ in this test. This small discrepancy goes away with increasing the number of particles per cell.

The ability to maintain rigorous conservation of momentum and energy for Tests T2a-T2e is shown in Fig. 8. The relative change in total x-momentum and energy is shown in this figure to be on the order of machine-precision for each of these tests.

For the final test, electron-ion thermalization of a fully ionized carbon plasma is considered ($A=12$, $Z=6$). The electron and ion distributions are initialized as non-drifting Maxwellians with the following conditions: $n_{i}=1{\times}10^{23}/$cm³, $T_{i}=50$ eV, $n_{e}=6{\times}10^{23}/$cm³, and $T_{e}=150$ eV. The equilibrium solution for the temperature is $T_{f}=136$ eV. The simulations are performed using a 2D periodic domain on a $40{\times}40$ grid with uniform spacing equal to $1.3285$ nm. For all simulations, $\Delta t=3.54{\times}10^{-4}$ fs is used for the time step and $\ln\Lambda=3$ is used as the Coulomb logarithm. This time step is roughly $100$ times smaller than the characteristic electron collision time: $\tau_{e}[\text{s}]=3.44{\times}10^{5}T_{e}^{3/2}[\text{eV}]/\left(Zn_{e}[1/\text{cm}^{3}]\Lambda\right)$.

Simulation results for this test problem are shown in Fig. 9. Here, time profiles for the average electron and ion temperatures on the 2D grid are shown for three cases. The number of ion particles per cell is set to $N_{i}=512$ for each of these simulations. The results in the left panel are obtained using equally weighted particles (T3a), in which case the number of electron particles per cell is $N_{e}=6N_{i}=3072$. The simulations corresponding to the middle and right panels use $N_{e}=N_{i}=512$, in which chase $w_{e}=6w_{i}$. The results in the middle and right panels are obtained without (T3b) and with (T3c) the moment-correction algorithm, respectively. The electron and ion temperature profiles from each of these simulations are observed to match well. Furthermore, the relaxation of $T_{e}$ and $T_{i}$ shown in Fig. 9 are also observed to agree well with results obtained from a 0D Spitzer model:

where $\nu_{T}=2m_{e}/M_{i}/\tau_{e}$ is the thermalization rate.

The ability to use weighted particles while conserving momentum and energy for test T3c is shown in Fig. 10. The left two panels in this figure show the change in total momentum in each direction. The amplitude of the change in momentum for test T3c shown in the middle panel is about 2-3 times larger than that for the simulation with equal weights shown in the left panel. In both cases, the change in momentum is close to machine precision and behaving like a random walk. Similarly, the relative change in energy for test T3c also behaves like a random walk, but with a step size that is about $10{\times}$ larger than for test T1a where equal weight particles are used. The reason for this is unclear. It may be due to the large electron-ion mass ratio. In either case, energy is conserved nearly to machine-precision for weighted-particles with the moment correction algorithm applied.

It is worthwhile to analyze how the energy-correction scheme described above perturbs the species distribution functions. For simplicity, assume $\Delta E$ after momentum is restored is negative. Using Eq. 23 in Eq. 22, we can write

It is additionally assumed that the particles are sorted by weight such that $m_{1}\approx m_{2}$. Then, using Eq. 26, the post-energy-correction particle velocities in Eqs. 19-20 become

We analyze Eqs. 27-28 in several different scenarios. First, assume the velocity of the two particles have similar magnitudes, but are in opposite directions. In this case, the velocity of both particles is scaled by $1+f_{E}/2$. If the velocity of the two particles is similar in magnitude and are in the same direction, then there is essential no modification to the particle velocities because the CM energy is approximately zero. Another case to consider is where the particle speeds are much different, e.g., $|\textbf{v}_{1}|\ll|\textbf{v}_{2}|$. In this scenario, the velocity of the fast particle is modified by $1+f_{E}/4$, while the speed of the slower particle gets enhanced to a value as large as $|\textbf{v}_{2}|f_{E}/4$. In other words, for this scenario the fast particle is only perturbed slightly while the slow particle can be perturbed significantly. We comment that more strongly perturbing slower particles is preferred as perturbations to the lower velocity part of the distribution function typically has less of an effect on the physics.

The energy-correction method works by adjusting the weighted-CM energy of random binary pairs after the momentum correction has been applied. If the violation of energy conservation, $\Delta E$ is large in magnitude compared to the combined weighted-CM energy of all binary pairs, then the method can fail. We find that this can occur when there are very few particles in a cell (e.g., $N$ order 10), and at least one of the particles has a weight that is much larger than the others. When the large weight particle is determined to scatter, it can cause a relatively large violation in energy conservation in the cell that cannot be re-absorbed by adjusting the weighted-CM energy of the very few binary pairs in the cell. We have only encountered this problem with the number of particles in cell being order $10$. For such a low number of particles in a cell, one cannot expect physically accurate results for scattering with weighted particles. The simplest solution is to do nothing and live with a momentum-conserving but not energy-conserving solution. The next simplest solution is to save the momentum vectors of the particles prior to scattering and, if the method fails, reset the particles velocities and skip doing scattering in this cell during this time step.

The moment-correction algorithm increases the computational cost of the overall scattering algorithm. For test T1d, but with $10{\times}$ more particles ($N=82,000$), the relative computational cost for no moment-correction, moment-correction without sorting, and moment-correction with sorting, are is 1.00, 1.23, and 1.72, respectively. The computational cost of applying the moment-correction algorithm itself increases the cost by $23$%. The cost increases by an additional $40$% when sorting the particles. We are currently using the standard std::sort() function, which may not always have the ideal computational cost of order $N\log\left(N\right)$. In the future, we will seek out more efficient sorting algorithms to reduce the computational cost. Furthermore, it is not always necessary to sort all of the particles (see Fig. 5). The computational burden of sorting can also be greatly reduced by only doing partial sorting.

A method for doing weighted-particle Coulomb collisions in a plasma that preserves scattering physics while exactly conserving momentum and energy is described in this manuscript. We first show how to extend $N{\times}N$ scattering methods for equally weighted particles to arbitrarily-weighted particles such that the normalized scattering length is correct for each particle on average. Then, it is straightforward to deduce an expression for the scattering length of binary pairs for order $N$ methods that is a direct generalization of that given by Takizuka and Abe [2] to include weighted particles.

The weighted-particle scattering method alone only preserves momentum and energy on average. A new method is introduced to rigorously restore exact momentum and energy conservation post scatter. The method, which is based on inelastic scattering dynamics, works by forming random binary pairs and making small adjustments to the weighted CM energy of each pair. The method works for both non-relativistic and relativistic particles. The efficacy of the entire algorithm is illustrated with various test problems, including those in which the original Nanbu method for weighted-particle Coulomb scattering [6] is known to produce incorrect results [7].

While this work focused on small angle Coulomb collisions, the extension to weighted particles and the method to restore momentum and energy conservation are not specific to Coulomb collisions in any way. These methods work equally well for other collision processes too, such as those governed directly by the Boltzmann collision integral discussed at the beginning of Sec. 2 where the normalized path length is used as a probability of a collision occurring rather than as a variance for a Gaussian distribution function from which a scattering angle is sampled from.