Coupling Fluid Plasma and Kinetic Neutral Models with Deterministic Monte Carlo Simulation*

Tokamak boundary plasma simulations often combine fluid plasma models with kinetic neutrals represented by Monte-Carlo methods.  The statistical nature of Monte-Carlo simulations, however, prevents coupled solvers from converging to an exact steady-state equilibrium; instead, solutions typically display stochastic oscillations  that are not only noisy when compared to true equilibria, but are also biased. Naively, this results in an error floor many orders of magnitude above the floating point error. Mathematically, this problem arises from the fact that the result of Monte Carlo simulations does not depend continuously on the plasma parameters. Here, we introduce ``deterministic'' Monte-Carlo simulations that eliminate this discontinuity. The key ingredients are: 1) coupled iterative schemes which fix random seeds every iteration, and 2) finite-size particle methods. We investigate convergence to floating-point accuracy for both one and two-dimensional coupled plasma-neutral simulations using simplified models that can capture charge-exchange, ionization, re-combination, and a variety of boundary behaviors.