Numerical thermalization in $n$-D particle-in-cell simulations

Technical Summary: Numerical Thermalization in n-D Particle-in-Cell Simulations

Problem Statement
Particle-in-cell (PIC) simulations are widely used to model kinetic plasma behavior by representing physical particles as macroparticles with finite shape functions. While often assumed to solve the collisionless Vlasov equation due to the smoothing effects of grid interpolation, PIC simulations frequently suffer from artificially high collisionality. This "numerical thermalization" arises when the macroparticle weight (the ratio of physical particles to macroparticles) is large, amplifying inter-particle forces despite the grid smoothing. In higher dimensions (2D and 3D), achieving a numerical thermalization time ($\tau_{num}$) longer than the physical thermalization time ($\tau_{phys}$) is often intractable, even when the grid resolves the Debye length. This artificial collisionality distorts velocity distributions and affects key quantities of interest (e.g., transport, emission, breakdown) over simulation timescales $t \gtrsim \tau_{num}$, potentially leading to false convergence or incorrect physical predictions.

Methodology
The authors validate a kinetic theory for numerical collision times against electrostatic PIC simulations across 1D, 2D, and 3D dimensionalities.

1. Simulations: Single-species electrostatic PIC simulations were performed with uniform initial conditions. The velocity autocorrelation function (VACF), $Z(t)$, was calculated from particle trajectory data. The numerical collision time was defined as the exponential decay time of $Z(t)$.
2. Parameters: Simulations varied the macroparticle weight ($W$), grid spacing ($\Delta$), and macroparticle shape radius ($R$). Both unfiltered simulations (where $R = \Delta/2$) and filtered simulations (where $R > \Delta/2$ via $k$-space filtering) were conducted to isolate shape effects from grid heating.
3. Theory: The study employs the kinetic theory developed by Okuda, Birdsall, and Langdon. This framework incorporates the macroparticle shape function and weight into a Lenard-Balescu-like kinetic equation. The theory derives collision operators for shaped macroparticles, accounting for the modification of the inter-particle force kernel by the shape function $S(k)$ and the weight $W$.
4. Comparison: Theoretical predictions for the collision time were derived using dynamic, static, and slow-speed approximations of the drag integral. These were compared directly against the decay rates observed in the PIC simulations.

Key Contributions and Results

• Validation of Kinetic Theory: The paper demonstrates that the kinetic theory of Okuda, Birdsall, and Langdon accurately predicts the numerical collision time in 1D, 2D, and 3D PIC simulations. The theory correctly captures the dependence of $\tau_{num}$ on macroparticle weight, shape size, and dimensionality.
• Force Kernel Interpretation: The study clarifies that numerical collisions are driven by an effective inter-particle force where the finite shape function weakens and smooths the force, while the macroparticle weight amplifies it. The theory shows that the ratio of numerical to physical collision times is independent of species mass and charge ratios, depending primarily on $W$ and $R$.
• Dimensionality Effects: The results highlight that reaching the physical thermalization time is significantly more challenging in 3D than in 1D or 2D for macroparticle sizes that resolve the Debye length. In 3D, the collision time scales logarithmically with the plasma parameter $\Lambda$ in the limit of small particle radius ($R \ll \lambda_D$).
• Trade-off Analysis: The paper quantifies the trade-off between reducing numerical thermalization and maintaining physical accuracy. Increasing the macroparticle size ($R$) to reduce collisions (by smoothing forces) simultaneously weakens the mean Vlasov field and alters collisionless processes (e.g., wave dispersion and damping). There exists a "distorted field" limit where further increasing $R$ violates physical constraints on spatial resolution.
• Higher-Order Shapes: The analysis extends to higher-order shape functions (e.g., B-splines). It is shown that increasing the spline order is mathematically equivalent to increasing the effective particle width. Thus, using higher-order shapes does not fundamentally alter the scaling of numerical thermalization but effectively increases the particle size.
• Application to Real Systems: The authors apply the derived constraints to several existing PIC studies (e.g., ICF hotspots, magnetic reconnection, Penning discharges, Hall thrusters, streamer discharges, and ITER neutral beam injectors). In many 3D cases, the numerical collision time is found to be shorter than the physical collision time, suggesting that numerical thermalization may dominate or significantly influence the results unless specific mitigation strategies (like Monte Carlo collisions or extreme grid refinement) are employed.

Significance
The primary significance of this work is the provision of a validated, predictive framework for estimating the timescale of self-consistent Coulomb interactions in PIC simulations. By establishing that kinetic theory accurately predicts $\tau_{num}$ as a function of grid spacing and macroparticle weight, the paper provides PIC users with a means to:

1. Characterize Relevance: Determine a priori whether numerical thermalization will impact the specific physics of interest (e.g., whether $\tau_{num}$ is shorter than the simulation time or the physical collision time).
2. Set Constraints: Define permissible operating regimes for macroparticle weight and shape size to ensure that numerical artifacts do not obscure physical phenomena.
3. Avoid False Convergence: Recognize that convergence with respect to particles per cell does not guarantee the elimination of numerical thermalization; sufficient macroparticle counts are required to resolve the physical collision rate, not just noise features.

The paper concludes that while PIC methods are powerful, the "collisionless" assumption is often invalid in 3D simulations with standard parameters. Users must explicitly account for numerical thermalization to ensure the fidelity of long-timescale transport and collisional processes.