Numerical heating in particle-in-cell simulations with Monte Carlo binary collisions

The binary Monte Carlo (MC) collision algorithm is a standard and robust method to include binary Coulomb collision effects in particle-in-cell (PIC) simulations of plasmas. Here we show that the coupling between PIC and MC algorithms can give rise to (nonphysical) numerical heating of the system that significantly exceeds that observed when these algorithms operate independently. We argue that this deleterious effect results from an inconsistency between the particle motion associated with MC collisions and the work performed by the collective electromagnetic field on the PIC grid. This inconsistency manifests as the (artificial) stochastic production of electromagnetic energy, which ultimately heats the plasma particles. The MC-induced numerical heating can significantly impact the evolution of the simulated system for long simulation times ($\gtrsim {10}^3$ collision periods, for typical numerical parameters). We describe the source of the MC-induced numerical heating analytically and discuss strategies to minimize it.

Figure 1

Comparison between (a) pure PIC, (b) pure MC, and (c) PIC-MC simulations of a 1D uniform electron-proton plasma in thermal equilibrium. The total energy conservation ($\mathrm{\Delta }\epsilon /{\epsilon }_0$) and change in electron ($\mathrm{\Delta }{\epsilon }_e/{\epsilon }_0$) and proton ($\mathrm{\Delta }{\epsilon }_p/{\epsilon }_0$) kinetic energies are shown for each case. The change in electromagnetic energy ($\mathrm{\Delta }{\epsilon }_{\text{EM}}/{\epsilon }_0$) is shown only for the PIC and PIC-MC simulations, since the self-consistent electromagnetic field is turned off in the pure MC case. Note that the range of the vertical scale of the PIC-MC case is 10 times larger than the pure PIC and pure MC cases.

Figure 2

Self-consistent radiation emission by a single electron undergoing MC collisions with plasma ions in 2D PIC-MC simulations. The colormaps show the patterns of the out-of-plane component of the electric field ($E_z$), representing the electromagnetic emission of a single electron. The electron is traveling from left to right with velocity $v=\langle v\rangle =\sqrt{8/\pi }{v}_{\mathrm{the}}$, where $v_{\mathrm{the}}$ is the thermal electron velocity, and is represented by the black circle; its trajectory is traced by the black line. (a) Fields produced by the electron in a pure PIC simulation. Also shown are fields produced in PIC-MC simulations for (b) ${\nu }_{ei}=0.04{\omega }_{pe}$ and (c) ${\nu }_{ei}=0.2{\omega }_{pe}$.

Figure 3

MC-induced numerical heating in PIC-MC simulations of thermal equilibrium plasma. The circles represent measurements from OSIRIS PIC-MC simulations and the blue, orange, and red colors represent 1D, 2D, and 3D simulations, respectively. The measured heating rates are in agreement with the theoretical prediction of Eq. (7) (solid lines). The heating rate coefficients $C_D$ have been determined for 1D, 2D, and 3D simulations. Note that the intrinsic numerical heating rate associated with the PIC algorithm has been subtracted from that of PIC-MC simulations to obtain these measurements.

Figure 4

Deterioration of energy conservation in a PIC-MC simulation due to MC-induced numerical heating (black solid curve) and comparison with theoretical predictions (dashed and dash-dotted curves). The deterioration of energy conservation of a pure PIC simulation with the same numerical parameters is represented by the blue solid curve and remains below $2\%$ for the simulated time. Note that the black solid curve corresponds to the difference in global energy variation between the PIC-MC and pure PIC simulations, in order to isolate the heating contribution associated with the PIC-MC coupling from the innate heating in PIC simulations.

Figure 5

MC-induced numerical heating in a two-stream unstable plasma. (a) Early time evolution (${\omega }_{pe}t<400$) of the system, where the exchange between particle kinetic energy (green solid curve) and electrostatic energy (orange solid curve) due to the development of the instability is observed. During this time, identical global energy conservation curves are obtained for both PIC-MC (blue solid) and pure PIC (blue dashed) curves. This is because the effective MC-induced heating time is much larger than this timescale. When the system relaxes to thermal equilibrium at around ${\omega }_{pe}t\simeq 2\times {10}^4$, the effective electron-ion collision frequency increases and the MC-induced heating rate also increases. After this time we observe in (b) the linear increase in particle energy and associated deterioration of energy conservation at the rate given by ${\mathrm{\Gamma }}_{\mathrm{MC}}\left(T_{eq}\right)$ (black dash-dotted curve), where $T_{eq}$ is the temperature reached at thermal equilibrium.

Figure 6

(a) Fourier spectrum of white noise current density fluctuations that mimic the effect of MC collisions and (b) Fourier spectrum of electric field fluctuations produced by the same noisy current density distribution. These spectra are obtained by integrating Maxwell's equations in one dimension using a finite-difference time domain solver in the presence of a noisy current density distribution. (c) Using a low-pass filter on the electromagnetic fields at every $30{\nu }_e{i}^{-1}$ (every 500 time steps), a large fraction of MC-induced electromagnetic fluctuations is suppressed, leading to a strong reduction in the MC-induced heating rate. This is illustrated for a 1D spatially uniform electron-proton plasma in thermal equilibrium with the same physical and numerical parameters used in Fig. 1. The black solid (dashed) and blue solid curves correspond to the PIC-MC simulation without (with) filtering on the fields and a pure PIC simulation without MC collisions, respectively.

Figure 7

Illustration of the effect of using different pair selection strategies for the MC collisions on the MC-induced numerical heating. In this example, we simulate a 1D spatially uniform electron-proton plasma in thermal equilibrium with the same physical and numerical parameters used in Fig. 1. However, in this case collisions were only performed between electrons. Electron-proton and proton-proton collisions were turned off. Shown are the deterioration of energy conservation for the PIC-MC simulation with the standard random pairing strategy (black solid curve), the PIC-MC simulation with NN pairing (black dashed curve), and the pure PIC simulation (blue solid curve). The MC-induced heating is strongly reduced when using the NN pairing strategy, remaining only the innate PIC heating level. Note that the NN pairing strategy only reduces the MC-induced heating when collisions are performed between particles of the same charge-to-mass ratio.