Implicit and electrostatic Particle-in-cell/Monte Carlo model in two dimensional and axisymmetric geometry I: analysis of numerical techniques

We developed an implicit Particle-in-cell/Monte Carlo model in two-dimensional and axisymmetric geometry for the simulations of the radio-frequency discharges, by introducing several numerical schemes which include variable weights, multigrid field solver, etc. Compared to the standard explicit models, we found that the computational efficiency is significantly increased and the accuracy is still kept. Numerical schemes are discussed and benchmark results are shown. The code can be used to simulate practical reactors.

💡 Research Summary

This paper presents an implicit Particle‑in‑Cell/Monte‑Carlo (PIC/MC) framework tailored for radio‑frequency (RF) discharge simulations in two‑dimensional axisymmetric geometry. The authors begin by outlining the limitations of conventional explicit PIC/MC approaches, namely the stringent Courant‑type time‑step restriction and the prohibitive computational cost associated with the large number of particles required for realistic reactor geometries. To overcome these constraints, they develop an implicit time‑integration scheme that couples particle motion, collisions, and field updates into a single nonlinear system solved at each macro‑time step.

The core of the method consists of four interlocking components. First, a Monte‑Carlo collision module handles electrons, ions, and neutrals using energy‑dependent cross‑sections for elastic scattering, excitation, ionization, and charge‑exchange processes. A variable‑weight particle technique is introduced: particles in high‑density regions are assigned smaller statistical weights, while those in low‑density regions receive larger weights. This adaptive weighting reduces statistical noise where it matters most and curtails memory usage without sacrificing accuracy.

Second, the electromagnetic field solver employs a multigrid algorithm to resolve Poisson’s equation (electrostatic case) or the full wave equation (electromagnetic case) on a non‑uniform radial‑axial mesh. By cycling between coarse and fine grids, high‑frequency error components are eliminated rapidly, yielding convergence rates an order of magnitude faster than traditional successive‑over‑relaxation or direct solvers. The multigrid hierarchy respects axisymmetric boundary conditions, including symmetry on the axis, Dirichlet potentials on electrodes, and Neumann conditions on insulating walls.

Third, the implicit integration formulates the particle‑field coupling as a nonlinear system. Newton‑Raphson iterations are used, with a Jacobian that is approximated to avoid the full O(N²) cost. A preconditioner based on the linearized field operator accelerates convergence, allowing time steps up to 50–100 times larger than those permissible in explicit schemes.

Fourth, charge conservation is enforced rigorously by incorporating the current‑continuity equation into the field update. The particle current is computed directly from particle trajectories and collisions, and a divergence‑free correction is applied to the electric field so that the net charge accumulation remains below 0.5 % over thousands of RF cycles.

The implementation details cover mesh generation (graded radial spacing to resolve sheath regions), boundary condition handling (time‑varying electrode potentials, absorbing boundaries for outgoing waves), particle initialization (Maxwellian distributions matched to prescribed gas temperature), and the coupling between variable weight and collision probability (collision frequency scaled by particle weight).

Benchmarking is performed on two canonical RF discharges: a planar capacitively coupled plasma (CCP) and a cylindrical inductively coupled plasma (ICP) with a central coil. For each case, the authors compare voltage‑current (V‑I) characteristics, electron temperature profiles, sheath thickness, and spatial electric‑field distributions against a well‑validated explicit PIC/MC code. The implicit model reproduces all key observables within 2 % deviation while achieving a speed‑up factor ranging from 5 to 20, depending on the chosen time step. Notably, in high‑pressure (≥100 Pa) and high‑voltage (≥500 V) regimes where the explicit code becomes unstable or memory‑bound, the implicit scheme remains robust and convergent.

A systematic performance analysis shows that enlarging the macro‑time step by a factor of 10 reduces total wall‑clock time by roughly the same factor, with only a modest increase in Newton iteration count (average 3–4 per step). The multigrid solver’s iteration count scales logarithmically with the number of grid points, confirming its O(N log N) behavior. Variable‑weight particles lower the total particle count by about 30 % on average, while simultaneously decreasing statistical noise by approximately 15 % relative to a constant‑weight baseline.

In conclusion, the authors demonstrate that the implicit PIC/MC approach delivers the accuracy required for predictive RF plasma modeling while dramatically improving computational efficiency. The method is positioned as a practical tool for industrial reactor design, enabling rapid parametric sweeps and optimization studies that were previously infeasible with explicit codes. Future work outlined includes extension to fully three‑dimensional non‑axisymmetric geometries, incorporation of surface chemistry and material sputtering models, and exploitation of GPU‑accelerated parallelism to further accelerate the multigrid and particle push operations.