A Propagator-based Multi-level Monte Carlo Method for Kinetic Neutral Species in Edge Plasmas

We propose and investigate a new multi-level Monte Carlo scheme for numerical solutions of the kinetic Boltzmann equation for neutral species in edge plasmas. In particular, this method explicitly exploits a key structural property of neutral particle dynamics: the prevalence of frequent collisions for which the outgoing velocity is determined by local plasma parameters. Using this property, we derive a multi-level algorithm based on collision event propagator and show, both analytically and through numerical experiments, that it reproduces the results of standard Monte Carlo methods. We further demonstrate that, in the context of coupled plasma-neutral edge simulations employing correlated Monte Carlo, the proposed scheme retains trajectory correlation to machine precision as the system evolves, whereas conventional methods exhibit rapid decorrelation. These results indicate that the propagator-based multi-level Monte Carlo scheme is a promising candidate for use in fully implicit Jacobian-free Newton-Krylov (JFNK) solvers for coupled plasma-neutral systems.

💡 Research Summary

This paper introduces a novel multi-level Monte Carlo (MC) scheme designed to efficiently and accurately solve the kinetic Boltzmann equation for neutral species in the edge plasma of magnetic fusion devices. The core innovation lies in explicitly leveraging a fundamental physical property of neutral particle dynamics in this environment: the Velocity-Localizing Property (VLP). Specifically, during frequent resonant charge-exchange collisions, the outgoing velocity of a neutral particle is effectively determined by the local plasma parameters (e.g., taking on the thermal velocity of a local ion). This causes the velocity distribution to “collapse” to a known function of spatial position immediately after such collisions.

The authors construct a mathematical framework around this observation using kernel methods. They define two families of integral operators (kernels): the source-correction operator J_σ, which advances the spatial location distribution of particles after σ collisions, and the standard estimation operator K_σ, which tallies the time spent by particles during that evolution. Under the VLP assumption, J_σ can be separated into a product of a known velocity distribution N_p(x)(v) and a spatial kernel j_σ(x, x’). This critical reduction allows the action of J_σ to be represented by a “propagator” matrix P_σ that operates only on the real-space grid, completely bypassing the curse of dimensionality associated with the full phase space.

The algorithm proceeds in two conceptual levels. First, in the “source correction” level, the method solves for a “virtual source” distribution s⋆. This is the particle distribution which, when used as the initial condition for a single application of the standard MC estimator K_σ, yields the true equilibrium distribution n⋆ directly. Thanks to the VLP, s⋆ is obtained by solving a manageable linear system on the real-space mesh: (Id - P_σ) s⋆ = s_0, where s_0 is the physical source (e.g., from walls). This can be solved via direct linear algebra or as the fixed point of the contraction mapping s ↦ s_0 + P_σ s. Second, in the “reconstruction” level, a single, standard MC simulation (applying K_σ) is run starting from the computed virtual source s⋆ to obtain the final, detailed equilibrium distribution n⋆(x, v).

Numerical experiments in simplified 2D tokamak geometries confirm that the proposed propagator-based method reproduces the same equilibrium results as conventional MC techniques. Its most significant advantage is demonstrated in the context of coupled plasma-neutral simulations using correlated MC. The new scheme maintains trajectory correlation to machine precision as the coupled system evolves. This is because the propagator matrix P_σ provides a deterministic mapping for particle transitions between spatial cells for a fixed random seed. In stark contrast, conventional MC methods exhibit rapid trajectory decorrelation over time due to their inherent randomness at every collision step.

This preservation of correlation is a game-changer for compatibility with advanced solvers. It enables the formation of stable Jacobian-vector products, which are essential for Jacobian-Free Newton-Krylov (JFNK) methods. JFNK solvers can dramatically accelerate convergence to steady-state in large-scale implicit simulations but have been largely incompatible with traditional kinetic MC due to its noise. Therefore, the propagator-based multi-level MC scheme bridges a critical gap, offering a promising path toward high-fidelity, fully implicit, and computationally tractable integrated simulations of edge plasmas, which is a vital step for predictive fusion energy research.