  We present the first-of-its-kind coupling of a continuum full-f gyrokinetic
turbulence model with a 6D continuum model for kinetic neutrals, carried out
using the Gkeyll code. Our objective is to improve the first-principles
understanding of the role of neutrals in plasma fueling, detachment, and their
interaction with edge plasma profiles and turbulence statistics. Our model
includes only atomic hydrogen and incorporates electron-impact ionization,
charge exchange, and wall recycling. These features have been successfully
verified with analytical predictions and benchmarked with the DEGAS2 Monte
Carlo neutral code. We carry out simulations for a scrape-off layer (SOL) with
simplified geometry and NSTX parameters. We compare these results to a baseline
simulation without neutrals and find that neutral interactions reduce the
normalized density fluctuation levels and associated skewness and kurtosis,
while increasing auto-correlation times. A flatter density profile is also
observed, similar to the SOL density shoulder formation in experimental
scenarios with high fueling.
