Coupling Fluid Neutrals to Gyrokinetic Plasma Dynamics for Edge and <scp>SOL</scp> Turbulence Simulations

ABSTRACT Accurate modeling of turbulent transport in magnetic confinement fusion devices requires extending first‐principles gyrokinetic simulations from the core to the edge and scrape‐off layer (SOL), where additional physics—particularly plasma–neutrals interactions—must be included. Neutrals in these regions reshape radial profiles, influencing gradient‐driven instabilities, and modifying transport across the separatrix and in the SOL. We present a tightly coupled scheme between a continuum gyrokinetic plasma model and a simplified fluid model for neutrals. Specifically, we extend the full‐, electromagnetic, collisional code GENE‐X with a neutrals density evolution equation, where charge exchange is treated as a diffusive process. This hybrid approach represents a first step toward a consistent fluid treatment of neutrals within gyrokinetic plasma simulations. We focus here on the derivation and implementation of Krook‐type source terms for ionization, recombination, and associated radiative losses. These operators are constructed to ensure strict conservation of mass, total momentum, and total energy, including kinetic, binding, and radiative components. The coupling is verified in terms of its numerical convergence and in three‐species relaxation tests (with electrons, deuterium ions, and neutral deuterium), with a particular focus on its impact on plasma moments, as well as, the distribution functions. Results indicate that even with a simplified neutrals model, the coupling introduces kinetic modifications to the plasma, whose influence on macroscopic quantities will need to be assessed in forthcoming turbulence simulations.