# Kinetic modeling of neutral transport for a continuum gyrokinetic code

T. N. Bernard<EMAIL_ADDRESS>General Atomics, San Diego, CA 92186,
USA  F. D. Halpern  General Atomics, San Diego, CA 92186, USA  M. Francisquez
Princeton Plasma Physics Laboratory, Princeton, NJ 08543, USA  N. R. Mandell
MIT Plasma Science and Fusion Center, Cambridge, MA 02139, USA  J. Juno
Princeton Plasma Physics Laboratory, Princeton, NJ 08543, USA  Department of
Physics and Astronomy, University of Iowa, Iowa City, IA 52242, USA  G. W.
Hammett  Princeton Plasma Physics Laboratory, Princeton, NJ 08543, USA  A.
Hakim  Princeton Plasma Physics Laboratory, Princeton, NJ 08543, USA  G.
Wilkie  Princeton Plasma Physics Laboratory, Princeton, NJ 08543, USA  J.
Guterl  General Atomics, San Diego, CA 92186, USA

###### Abstract

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.

##  I Introduction

Neutral particles are present in the edge plasma of magnetic fusion devices
due to plasma interactions with wall and divertor materials, as well as beam
injection and neutral gas puffs. Atomic and molecular neutrals interact with
plasma particles via electron-impact ionization, charge exchange, radiative
recombination, etc. Neutrals can penetrate into the core and affect plasma
dynamics through these collisional processes. They play a role in plasma
fueling via recycling, whereby ions impinging on the wall are re-emitted as
neutrals, which then become ionized. Furthermore, ion–neutral friction is
necessary to reduce the heat flux to the divertor in detached
scenarios.Stangeby (2000) The effect of neutrals on confinement has been
investigated experimentallyBoivin  _et al._ (2000); Groebner  _et al._ (2002);
Moser  _et al._ (2020); Mordijck (2020) and accurate theoretical and numerical
models are necessary to thoroughly interpret those results.

These experimental studies have generated interest in coupling models of
neutral transport to first-principles plasma turbulence codes. Couplings with
fluid neutral models have effectively modeled some scrape-off layer (SOL)
conditions,Thrysøe  _et al._ (2018); Zholobenko  _et al._ (2021) but are
rigorous only where the neutral mean free paths are shorter than the
characteristic length scale. Comparisons between a fluid neutral code and a
kinetic model were improved recently with a more accurate numerical treatment
of geometry.Dekeyser  _et al._ (2019) Based on typical SOL conditions, neutral
mean free paths can vary from several centimeters to a meter or more.Stangeby
(2000) One characteristic length scale is the size of coherent turbulent
structures in the plasma, or blobs, which are on the order of
centimeters.D’Ippolito, Myra, and Zweben (2011) Kinetic models for neutrals
are necessary to accurately capture both long and short mean free path
neutrals. Because of this, much neutral modeling has been carried out by
kinetic Monte Carlo (MC) codes,Reiter (2002); Stotler and Karney (1994);
Cupini, De Matteis, and Simonini (1983) and these models have been coupled to
fluid turbulence codes like TOKAM2D Marandet  _et al._ (2013) and
TOKAM3X.Galassi  _et al._ (2017); Fan  _et al._ (2019) SOLPS has been coupled
to the six-field two-fluid turbulence code BOUT++ Zhang  _et al._ (2019), and
the XGC framework incorporates a model of MC neutrals.Stotler  _et al._ (2013,
2017); Ku  _et al._ (2018); Hager  _et al._ (2019) MC codes are subject to
statistical noise, which can interfere with the accuracy and convergence of
Eulerian codesJoseph  _et al._ (2017) when coupled to them. Thus, various
continuum kinetic models for neutrals have also been developed and coupled to
first-principles turbulence codes such as GBSWersal and Ricci (2015, 2017);
Wersal, Ricci, and Loizu (2017); Mandell  _et al._ (2020); Mancini  _et al._
(2021) and nSOLT.Russell, Myra, and Stotler (2019); Russell  _et al._ (2021)

Gyrokinetic models, which have long been used to model core plasma turbulence,
have been more recently adapted to the edge and SOL to capture important
kinetic effects, such as parallel transport, trapped particles, nonlinear
wave–particle interactions, etc., that drift-reduced fluid and gyrofluid codes
may model less accurately. Pan, Told, and Jenko (2016); Chang  _et al._
(2017); Shi  _et al._ (2019); Boesl  _et al._ (2019); Dorf and Dorr (2020);
Mandell  _et al._ (2020); Michels  _et al._ (2021) Including neutral
interactions further improves the predictive capabilities of these models.
Thus, in this work we present a continuum kinetic model of neutral transport
coupled to a continuum gyrokinetic solver within the computational plasma
physics framework Gkeyll,Gkeyll (2022) which has previously been used to model
plasma turbulence on open field lines.Shi  _et al._ (2017, 2019); Bernard  _et
al._ (2019); Mandell  _et al._ (2020); Hakim  _et al._ (2020a); Bernard  _et
al._ (2020) The collisionless Vlasov solver in GkeyllJuno  _et al._ (2018);
Hakim and Juno (2020) facilitated the implementation of this model. A
continuum kinetic neutral model avoids the statistical noise issues that are
associated with MC codes and the shortcomings of fluid neutral transport
models. The model only includes atomic neutrals and the electron-impact
ionization and charge exchange processes. We model wall recycling as a
boundary condition in the direction parallel to the magnetic field.
Recombination will eventually be included, when high-recycling divertor
scenarios are considered. For plasmas considered in this work, where
$T_{e}\geq 10$ eV, this interaction is negligible. Other important reactions
to include in future work are radiation, molecular processes, ion–neutral
elastic collisions, and neutral–neutral collisions.

The results presented here comprise the first-of-its-kind coupling of a
continuum gyrokinetic solver to a continuum kinetic model for neutral
transport. We have verified the neutral model against analytic theory in low
dimensional tests with static plasma species. Benchmark tests, also in low
dimensions with static plasma species, show good agreement with the DEGAS2 MC
neutral code.

We have demonstrated proof-of-concept in a high-dimensional simulation in
three physical space dimensions, and two and three velocity space dimensions
for the plasma species and neutral species, respectively (3X2V+3X3V). We model
the National Spherical Torus Experiment (NSTX) SOL with simplified helical
geometry and dynamic plasma species, similar to previous NSTX simulations with
Gkeyll.Shi  _et al._ (2019); Mandell  _et al._ (2020); Mandell (2021) We
compare the results to a simulation without neutrals and with the same
midplane plasma source. As expected, the inclusion of neutrals results in an
increase in the steady state density due to ionization sourcing. A flattening
of the density profile is also observed, and the ratio of the radial to
parallel flux increases. These results support the conclusion in Ref. 26 that
a flatter density profile can form in parameter regimes where the radial flux
becomes more efficient at transporting density relative to the parallel flux.
The electron and ion temperatures decrease due to energy loss occurring
through the ionization process, transfer of energy to neutrals via charge
exchange, and energy exchange via electron–ion collisions. In this simulation,
we find that neutrals do not have a large impact on the parallel heat flux
width. Power balance calculations demonstrate that energy is conserved to
within 5% of the input power.

The article is organized as follows. Section II describes the model equations
used in Gkeyll and the simplified model for neutral transport that has been
implemented. Section III contains results from SOL tests in one spatial
dimension with static plasma species, including comparisons with analytic
predictions and benchmarks with DEGAS2. Section IV contains results from high-
dimensional simulations of a SOL with NSTX parameters and an analysis of the
effect of neutrals. Section V summarizes our findings. The simulations that
produce the results contained in this paper can be reproduced with the input
files made available online (see Appendix A).

##  II Model equations for coupled gyrokinetic plasma and kinetic neutral
species

In this section we describe the equations that Gkeyll solves, including
conduction sheath and wall recycling boundary conditions. The simplified
models for electron-impact ionization and charge exchange that have been
implemented are also explained in detail here.

###  II.1 Plasma species

To model plasma transport on open-field lines in fusion devices, Gkeyll
evolves the full-$f$ gyrokinetic distribution function
$f_{s}(\bm{x},v_{\parallel},\mu,t)$ in the long-wavelength (or drift-kinetic)
limit using the following equation:

$\displaystyle\frac{\partial\mathcal{J}_{s}f_{s}}{\partial
t}+\nabla\cdot(\mathcal{J}_{s}\\{\bm{R},H\\}f_{s})+\frac{\partial}{\partial
v_{\parallel}}(\mathcal{J}_{s}\\{v_{\parallel},H\\}f_{s})$
$\displaystyle=\mathcal{J}_{s}C[f_{s}]+\mathcal{J}_{s}S_{s}.$ (1)

While Gkeyll can handle electromagnetic fluctuations,Mandell  _et al._ (2020)
for simplicity here we will consider only electrostatic fluctuations. The
collision operator $C[f_{s}]$ can include Coulomb collisions, elastic
collisions, and inelastic neutral interactions, and $S_{s}$ is a source term.
Coulomb collisions are modeled using a Dougherty collision operator
Francisquez  _et al._ (2020) and collisional interactions with neutrals are
described in the following subsection. The Jacobian is
$\mathcal{J}=B^{*}_{\parallel}=\bm{b}\cdot\bm{B}_{\parallel}^{*}$, where
$\bm{B}_{\parallel}^{*}=\bm{B}+(Bv_{\parallel}/\Omega_{s})\nabla\times\bm{b}$.
The Poisson brackets are defined by

$\displaystyle\\{F,G\\}=\frac{\bm{B}^{*}}{m_{s}B_{\parallel}^{*}}\cdot\left(\nabla
F\frac{\partial G}{\partial v_{\parallel}}-\frac{\partial F}{\partial
v_{\parallel}}\nabla G\right)$
$\displaystyle-\frac{1}{q_{s}B_{\parallel}^{*}}\bm{b}\cdot\nabla F\times\nabla
G.$ (2)

In the long-wavelength limit, the Hamiltonian is
$H_{s}=\frac{1}{2}mv_{\parallel}^{2}+\mu B+q_{s}\phi$, and the system is
closed by the gyrokinetic Poisson equation

$-\nabla\cdot\left(\frac{n_{i0}^{g}q_{i}^{2}\rho_{\mathrm{s}0}^{2}}{T_{e0}}\nabla_{\perp}\phi\right)=\sigma_{g}=q_{i}n_{i}^{g}(\bm{R},t)-en_{e}(\bm{R},t),$
(3)

which uses a linearized form for the polarization term, similar to the
Boussinesq approximation made in Braginksii fluid codes.

A non-orthogonal field-line-following coordinate system is used for the
configuration space grid, where $z$ is the distance along the field, $x$ is
the radial coordinate, and $y$ is the binormal coordinate. Plasma and neutral
species share the same configuration space grid. For the plasma species,
Dirichlet boundary conditions with $\phi=0$ are used in the radial direction
and periodic boundary conditions are used in the binormal direction.
Conducting-sheath boundary conditions are applied to $f$ in the $z$ direction,
by which one plasma species is partially reflected and the other is completely
absorbed.Shi  _et al._ (2017) The gyrokinetic Poisson equation 3 is solved
across the entire spatial domain, and the values at $z_{\min}$ and $z_{\max}$
(or $x_{\min}$ and $x_{\max}$ in Sec. III) give the electrostatic potential at
the sheath entrance. For simulations presented in subsequent sections, the
sheath potential is positive, and it thus determines a cutoff velocity for the
electron species, below which they are reflected. Ions free stream out of the
domain at the parallel boundaries and are recycled as neutrals when coupled to
the Vlasov neutral solver.

###  II.2 Neutral species

The neutral particle distribution function $f_{n}(\bm{x},\bm{v},t)$ is evolved
using the Vlasov equation

$\frac{\partial f_{n}}{\partial t}+\bm{v}\cdot\nabla f_{n}=C[f_{n}]+S_{n},$
(4)

where $C[f_{n}]$ can represent elastic collisions (not included in present
results), which are modeled by a Bhanagar-Gross-Krook (BGK) operator, and also
inelastic neutral interactions such as ionization and charge exchange. $S_{n}$
is a volumetric source term. As previously mentioned, the neutral distribution
function is evolved on the same configuration space grid as the plasma
species. Though the neutrals are evolved in the field-line-following
configuration space, they are not affected by the magnetic field. The velocity
space grid includes three velocity space coordinates which are orthogonal
($v_{R},v_{Z},v_{\varphi}$), where subscripts ($R,Z,\varphi$) correspond to
the coordinates in a cylindrical system. Simulations in this work assume a
simplified helical magnetic geometry, and the details of the neutral dynamical
equation in this geometry are contained in Appendix B.

A model for wall recycling has been implemented at the boundaries where field
lines terminate at the endplate. These boundary conditions provide a source of
neutrals from the targets that depends on the flux of incident ions. Consider
a simulation with one configuration space dimension, parallel to the magnetic
field, and three velocity space dimensions (1X3V). Coordinate $z$ is parallel
to the magnetic field, and $v_{z}$ is the velocity for neutrals corresponding
to that direction. We define the neutral distribution function in the ghost
cell at the lower ($z_{\min}$) boundary as

$\displaystyle
f_{n}(v_{x},v_{y},v_{z},z=z_{ghost})=\alpha_{b}f_{n}(v_{x},v_{y},-v_{z},z_{\min})$
$\displaystyle+C_{rec}f_{M,rec}(T=T_{n,rec}),$ (5)

where $\alpha_{b}$ is the reflection coefficient. The unit-density Maxwellian
function for recycled neutrals $f_{M,rec}$ is defined by a zero mean flow and
a user-specified temperature of $T_{n,rec}=10$ eV for simulations presented
here. (Typical Franck-Condon temperatures of 2 eV require a finer velocity
space grid and are computationally expensive. Non-uniform velocity space grids
would alleviate this issue and are being explored for future simulations.) The
Maxwellian is scaled such that the magnitude of the flux of incoming neutrals
is equal to the magnitude of the flux of incident ions.

$\displaystyle
C_{rec}\int_{0}^{v_{z,\max}}\int_{v_{y,\min}}^{v_{y,\max}}\int_{v_{x,\min}}^{v_{x,\max}}dv_{x}dv_{y}dv_{z}\;v_{z}\,f_{M,r}$
$\displaystyle\doteq\alpha_{rec}\frac{2\pi}{m}\int_{0}^{\mu_{\max}}\int_{v_{\parallel,i,\min}}^{0}dv_{\parallel}d\mu\;v_{\parallel,i}\mathcal{J}f_{i}(z=z_{\min})$
$\displaystyle C_{rec}=\frac{\alpha_{rec}M_{1z-,i}}{M_{1z+,n}},$ (6)

where $\alpha_{rec}$ is the user-specified recycling fraction, $\mathcal{J}$
is the Jacobian, $\doteq$ denotes ‘weak equality’ according to the modal
discontinuous Galerkin (DG) scheme used in GkeyllHakim  _et al._ (2020b), and
+ or - subscripts refer to partial moments calculated over $v_{z}>0$ or
$v_{z}<0$, respectively. The scaling coefficient $C_{rec}$ is calculated
dynamically in the simulation. An angular dependence in the recycling boundary
conditions will be included in future work. The recycling boundary conditions
in Gkeyll are not currently used to model recycling at the outer wall due to
the zero-flux, conducting wall boundary condition used for the plasma species
at the radial boundaries.

###  II.3 Plasma-neutral interactions

Here we present the details of the simplified neutral interaction models that
have been implemented in Gkeyll. Electron-impact ionization, described by the
process $e^{-}+n\rightarrow i^{+}+2e^{-}-E_{iz}$, can be rigorously modeled by
a collision term containing an integral over velocity space:

$\displaystyle C^{iz}_{e}=C^{iz}_{i}=-C^{iz}_{n}$
$\displaystyle=f_{n}(\bm{v})\int
f_{e}(\bm{v}^{\prime})\sigma_{iz}(|\bm{v}-\bm{v}^{\prime}|)|\bm{v}-\bm{v}^{\prime}|\,d\bm{v}^{\prime}.$
(7)

Assuming that the electron thermal speed is much greater than the neutral
thermal speed, $|\bm{v}^{\prime}|\gg|\bm{v}|$, this simplifies to

$\displaystyle\int
f_{e}(\bm{v}^{\prime})\sigma_{iz}(|\bm{v}-\bm{v}^{\prime}|)|\bm{v}-\bm{v}^{\prime}|\,d\bm{v}^{\prime}$
$\displaystyle\approx\int
f_{e}(\bm{v}^{\prime})\sigma_{iz}(v^{\prime})v^{\prime}\,d\bm{v}^{\prime}=n_{e}\langle\sigma_{iz}v_{e}\rangle,$
(8)

where $\langle\sigma_{iz}v_{e}\rangle$ is the ionization rate parameter, whose
value can be obtained from atomic databases for various elements. In Gkeyll,
the ionization reaction rate is approximated with a fitting formula Voronov
(1997)

$\langle\sigma_{iz}v_{e}\rangle=A\frac{1+P\,(E_{iz}/T_{e})^{1/2}}{X+E_{iz}/T_{e}}\left(\frac{E_{iz}}{T_{e}}\right)^{K}e^{-E_{iz}/T_{e}}\,\times
10^{-6}{\rm m}^{3}/{\rm s},$ (9)

where constants $A$, $P$, $X$ and $K$ are tabulated for elements up to $Z=28$.
$E_{iz}$ is the ionization energy and is 13.6 eV for simulations presented in
subsequent sections. This model has been used previously in Gkeyll as part of
Vlasov-Maxwell simulations of a plasma sheath.Cagas  _et al._ (2017); Cagas
(2018) It was included as a source term in the plasma species equations to
help the sheath relax to a steady state, and lower-energy electrons that
resulted from the process were not accounted for.

We model the ionization collision term for the electron species based on the
neutral model in Ref. 22

$C^{iz}_{e}=n_{n}\langle\sigma_{iz}v_{e}\rangle\left[2f_{M,iz}(n_{e},\bm{u}_{n},v^{2}_{t,iz})-f_{e}\right],$
(10)

where $\bm{u}_{n}$ is the neutral fluid velocity and $v^{2}_{t,iz}$ is the
thermal velocity of the resulting lower-energy electrons, given by
$v^{2}_{t,iz}=v^{2}_{t,e}/2-E_{iz}/(3m_{e})$. Finally, the ion and neutral
collisional ionization terms are given by

$\displaystyle C^{iz}_{i}$ $\displaystyle=n_{e}f_{n}\langle\sigma
v_{e}\rangle,$ (11) $\displaystyle C^{iz}_{n}$
$\displaystyle=-n_{e}f_{n}\langle\sigma v_{e}\rangle.$ (12)

These terms depend on moments of the electron and neutral distribution
function such as $n_{e}$, $v_{t,e}^{2}$, and $\bm{u}_{n}$, which are
calculated at every timestep. The electron temperature appearing in Eq. 9 is
obtained from the thermal speed, $v_{t,e}^{2}=T_{e}/m_{e}$.

The charge exchange process can be described by the collisional term

$\displaystyle C^{cx}_{i}=-C^{cx}_{n}=f_{n}(\bm{v})\int
f_{i}(\bm{v}^{\prime})\sigma_{cx}(|\bm{v}-\bm{v}^{\prime}|)|\bm{v}-\bm{v}^{\prime}|\,d\bm{v}^{\prime}$
$\displaystyle-f_{i}(\bm{v})\int
f_{n}(\bm{v}^{\prime})\sigma_{cx}(|\bm{v}-\bm{v}^{\prime}|)|\bm{v}-\bm{v}^{\prime}|\,d\bm{v}^{\prime},$
(13)

where we have assumed that ions and neutrals have the same mass. Instead of
computing this computationally expensive integral over velocity space, we have
implemented a simplified model based on Refs. 48; 49. The resulting
collisional charge exchange terms are

$C^{cx}_{i}=-C^{cx}_{n}=\sigma_{cx}V_{cx}(n_{i}f_{n}-n_{n}f_{i}).$ (14)

The effective relative velocity $V_{cx}$ is defined by

$V_{cx}^{2}\equiv\frac{4}{\pi}(v_{t,i}^{2}+v_{t,n}^{2})+(\bm{u}_{i}-\bm{u}_{n})^{2},$
(15)

where $\bm{u}_{s}$ and $v_{t,s}^{2}$ are the fluid velocity and squared
thermal speed for species $s$. Cross section data is taken from Ref. 50, and
we currently only consider hydrogen (H-H+) and deuterium (D-D+) reactions
using

$\displaystyle\sigma_{cx,H}=1.12\times 10^{-18}-7.15\times
10^{-20}\ln(V_{cx}),$ (16) $\displaystyle\sigma_{cx,D}=1.09\times
10^{-18}-7.15\times 10^{-20}\ln(V_{cx}).$ (17)

This model currently neglects the perpendicular momentum transfer from
neutrals to the $E\times B$ flow. When a neutral becomes an ion via charge
exchange, the location of the neutral and the guiding center of the resulting
ion differ by a distance of $\rho=\bm{u}\times\hat{\bm{b}}/\Omega_{ci}$. This
process is the same as the pickup ion effect that occurs in the solar
wind.Isenberg (1987) This important feature will be implemented in future
work.

As noted earlier, the plasma species and the neutral species distribution
functions are evolved on the same physical space grid but different velocity
space grids. Since $f_{n}$ appears in the ion collisional terms $C^{iz}_{i}$
and $C^{cx}_{i}$ and $f_{n}$ appears in $C^{cx}_{n}$, an interpolation scheme
in velocity space is necessary to project the distribution functions onto
another phase-space grid. In a discontinuous Galerkin scheme, interpolation
between velocity space grids whose cell and domain boundaries do not overlap
is not trivial, especially when one wishes to avoid conservation and aliasing
errors as done previously.Hakim and Juno (2020) While complexity of this
development effort is outside the scope of this work, it will be implemented
in the future in order to retain all kinetic effects.

For now we avoid such a scheme by assuming that distribution functions are
approximately Maxwellian. Thus, fluid moments $n_{i}$,
$\bm{u}_{i}=u_{\parallel,i}\hat{\bm{b}}$, and $v^{2}_{t,i}$ are taken from the
ion distribution function $f_{i}$ in order to project a Maxwellian function
$f_{M,i}(n_{i},\bm{u}_{i},v^{2}_{t,i})$ onto the neutral Vlasov grid. (Here
$\hat{\bm{b}}=\bm{B}/|\bm{B}|$ is the magnetic field unit vector.) Similarly
fluid moments $n_{n}$, $u_{\parallel,n}=\bm{u}_{n}\cdot\hat{\bm{b}}$, and
$v^{2}_{t,n}$ are used to project a Maxwellian function for neutrals
$f_{M,n}(n_{n},u_{\parallel,n},v^{2}_{t,n})$ onto the ion gyrokinetic grid.
Then the collisional neutral interaction terms become

$\displaystyle C^{iz}_{i}$ $\displaystyle=n_{e}f_{M,n}\langle\sigma
v_{e}\rangle$ (18) $\displaystyle C^{cx}_{i}$
$\displaystyle=\sigma_{cx}V_{cx}(n_{i}f_{M,n}-n_{n}f_{i})$ (19) $\displaystyle
C^{cx}_{n}$ $\displaystyle=-\sigma_{cx}V_{cx}(n_{i}f_{n}-n_{n}f_{M,i}).$ (20)

The models for ionization and charge exchange conserve particles, momentum and
energy in the system, a property that was confirmed in one-dimensional tests
with periodic boundary conditions.

##  III Neutral model verification and benchmarks

The neutral model, as implemented in Gkeyll, has been verified with analytic
theory and benchmarked against the DEGAS2 Monte Carlo neutral code, which is
detailed in this section. In these simulations, the gyrokinetic plasma species
are static (not evolved) and only the neutrals are evolved. Simulation end
times were chosen such that the neutral species achieved a steady state,
verified by the fact that neutral profiles were no longer evolving.

###  III.1 Verification with analytic solution

To verify the neutral interaction models, we derive theoretical predictions of
steady-state profiles for the neutral distribution function in one
configuration space and one velocity space (1X1V) with a static background
plasma. In 1X1V the steady-state neutral Vlasov equation with either charge
exchange or ionization can be rewritten as

$v_{z}\frac{\partial f}{\partial z}=-\nu\left(f(z,v_{z})-g(z,v_{z})\right),$
(21)

where $f$ is the neutral distribution function. When the RHS of Eq. 21
represents charge exchange, $\nu=n_{i}\sigma_{cx}V_{cx}$ and $g$ is
proportional to the ion distribution function

$g(z,v_{z})\equiv\frac{f_{i}(z,v_{z})}{n_{i}}\times\int
dv_{z}^{\prime}f(z,v_{z}^{\prime}).$ (22)

In the case of ionization, $\nu=n_{e}\langle\sigma v_{e}\rangle$ and $g=0$. If
we consider a finite symmetric domain with $z=[-L/2,L/2]$, the general
solution to Eq. 21 is given by the piecewise exponential function

$f(z,v_{z})=\left\\{\begin{array}[]{ll}C_{+}(z,v_{z})e^{-\frac{\nu}{v_{z}}(z+L/2)}&v_{z}>0\\\
C_{-}(z,v_{z})e^{-\frac{\nu}{v_{z}}(z-L/2)}&v_{z}<0,\end{array}\right.$ (23)

where $C_{\pm}$ can be determined from the imposed boundary conditions. A
detailed derivation is given in Appendix C. The zeroth moment of Eq. 23 gives
a theoretical prediction for density, which can be directly compared to
numerical results.

To verify models of neutral interactions in Gkeyll, tests were conducted in
1X1V with either charge exchange or ionization. The recycling coefficient was
$\alpha_{rec}=1.0$ and the reflection coefficient was $\alpha_{b}=0$. A
hydrogen plasma is assumed, and input parameters are the following:
$B_{0}=0.5$ T, $n_{0}=5\times 10^{18}$ m-3, and $T_{e0}=T_{i0}=20$ eV,
$T_{n0}=2$ eV. The connection length is $L_{z}=40$ m. The resolution is
$(N_{z},N_{v_{z}})=(224,32)$. The simulation was run to
$t_{end}=L_{z}/c_{s0}$, where $c_{s0}=\sqrt{T_{e0}/m_{i}}$. Static plasma
conditions are $n_{e}=n_{i}=n_{0}$, $T_{e}=T_{i}=T_{e0}$, and $u_{\parallel
e}=u_{\parallel i}=-c_{s0}$ for $z\leq 0$ or $c_{s0}$ for $z>0$. The neutral
initial conditions are

$n_{n}(t=0)=\left\\{\begin{array}[]{ll}n_{0}({\rm
sech}^{2}(-(L_{z}/2+z-2)/0.2)+10^{-6})&\quad z\leq 0\\\ n_{0}({\rm
sech}^{2}((L_{z}/2+z-2)/0.2)+10^{-6})&\quad z>0,\end{array}\right.$ (24)

with $T_{n}(t=0)=T_{n0}$ and $u_{n}(t=0)=0$. Density initial conditions were
chosen to approximate steady-state conditions, which include an exponential
decay away from the $z$-boundaries, so that the actual steady state could be
achieved more quickly in the code. We considered the case with only charge
exchange and a constant value for the reaction rate
$\sigma_{cx}V_{cx}=2.2\times 10^{-14}$ m3/s. A narrow, shifted Maxwellian
sourced the recycled neutrals at the boundary

$f_{M,rec}=\frac{1}{\sqrt{2\pi v_{t,n0}^{2}}}\;e^{-(v_{z}\mp
c_{s})^{2}/(2v_{t,n}^{2})},$ (25)

where $v_{t,n0}^{2}=T_{n0}/m_{i}$, and the minus and plus signs correspond to
the minimum and maximum boundaries in $z$, respectively. A scaling factor
$C_{rec}=4.98\times 10^{18}$ m-3 was used such that the incoming neutral flux
matched the incident ion flux of the static plasma. The symmetric solution for
the density, Eq. 42, was solved numerically and compared to the steady-state
solution from Gkeyll. This comparison is shown in Fig. 1a. The $L^{2}$-norm of
the error, normalized to the $L^{2}$-norm of the solution predicted by theory,
is shown in Fig. 1b for Gkeyll resolutions $N_{z}=112,168,224$. The error
decreases with increasing resolution.

The same test case was run with ionization only and the resulting density was
compared to the analytic solution, Eq. 42 with $\beta_{i}=0$. This comparison
is shown in Fig. 2a with the associated $L^{2}$-norm of the error in 2b,
demonstrating convergence. Both cases show excellent agreement with
theoretical predictions of density profiles.

Figure 1: Verification of the charge exchange model in Gkeyll. Assuming a
simplified model for the charge exchange cross section (Eq. 20), the
theoretically predicted steady-state density profile is compared with the
resulting profile from Gkeyll in (a). The $L^{2}$-norm of the error relative
to the $L^{2}$-norm of the theoretical prediction is shown in (b). Figure 2:
Verification of the ionization model in Gkeyll. Assuming a simplified
ionization model (Eq. 12), the theoretically predicted steady-state density
profile is compared with the resulting profile from Gkeyll in (a). The
$L^{2}$-norm of the error relative to the $L^{2}$-norm of the theoretical
prediction is shown in (b).

###  III.2 Benchmarks with DEGAS2

Additional tests with static plasma species were carried out and compared with
results from the DEGAS2 Monte Carlo neutral code. These tests were run with
low and high densities with both ionization and charge exchange and also with
each neutral interaction separately. Hydrogen ions are assumed, and input
parameters are $B_{0}=0.5$ T, $n_{0}=5\times 10^{18}\,(10^{19})$ m-3,
$T_{e0}=30$ eV $T_{i0}=60$ eV, $T_{n0}=10$ eV, and $L_{z}=40\,(10)$ m, with
changes to the high density parameters denoted by parentheses. The resolution
is $(N_{z},N_{v_{x}},N_{v_{y}},N_{v_{z}})=(224,16,16,16)$ for the low density
case and $(N_{z},N_{v_{x}},N_{v_{y}},N_{v_{z}})=(448,16,16,16)$ for the high
density case. Both were run to an end time of $t_{end}=3L_{z}/c_{s0}$, where
$c_{s}=\sqrt{T_{e0}/m_{i}}$.

Static plasma parameters are $n_{e}=n_{i}=n_{0}$, $T_{e}=T_{e0}$,
$T_{i}=T_{i0}$, and $u_{\parallel e}=u_{\parallel i}=\frac{c_{s}z}{L_{x}/2}$.
The neutral initial conditions are

$n_{n}(t=0)=\left\\{\begin{array}[]{ll}n_{0}({\rm
sech}^{2}(-(L_{z}/2+z-2)/0.2)+10^{-6})&\quad z\leq 0\\\ n_{0}({\rm
sech}^{2}((L_{z}/2+z-2)/0.2)+10^{-6})&\quad z>0,\end{array}\right.\\\ $ (26)

with $T_{n}(t=0)=T_{n0}$ and $\bm{u}_{n}(t=0)=(0,0,0)$. A Maxwellian neutral
source function $S_{n}$ with density source rate $8\times 10^{21}$ (m3/s)
$n_{0}^{2}$ and temperature $T_{n0}$ is used to approximate a source floor due
to recombination. Recycling wall boundary conditions were used for neutrals
with a recycling coefficient of $\alpha_{rec}=1.0$ and no reflection of
neutrals ($\alpha_{b}=0$). Both high and low density tests were run on 112
cores and took 44 and 183 minutes, respectively.

Figure 3: A comparison of benchmark tests from Gkeyll and DEGAS with static
plasma species that include charge exchange, ionization and a recycling source
at the walls. Steady-state profiles of neutral density (a) and temperature (b)
are shown from the low plasma density case. Neutral density (c) and
temperature (d) are also compared from the high plasma density case.

DEGAS2 was run at identical conditions for the high density and low density
benchmark cases. The main differences are the lack of higher-order information
in the cells and the use of internal tables for the charge exchange cross
section and ionization rates. Also, DEGAS2 is run in 3X3V instead of 1X3V,
albeit with symmetry in physical space. The domain is divided into 896 and
1792 zones, wherein properties are constant, in the low and high density
cases, respectively. This is higher than the Gkeyll resolution because only
cell-averages are available in DEGAS2, whereas the Gkeyll results contain
higher-order spatial information. To minimize noise in the region of low
neutral density, 16M trajectories were used for each of two source types,
volumetric and surface. With 128 cores, each DEGAS2 calculation took about
10-20 minutes. These tests were solving the steady-state equations, while the
Gkeyll tests dynamically evolved to a steady state. This comparison focuses on
benchmarking the model in Gkeyll and not on performance; thus, no efforts were
made to optimize either DEGAS2 or Gkeyll for these tests.

Figure 4: A comparison of benchmark tests from Gkeyll and DEGAS with static
plasma species that include only ionization and a recycling source at the
walls. Steady-state profiles of neutral density (a) and temperature (b) are
shown from the low plasma density case. Neutral density (c) and temperature
(d) are also compared from the high plasma density case.

Since the ionization reaction rate is constant in time and space, the values
calculated from the DEGAS2 model were used in the Gkeyll cases. In DEGAS2,
only the ground state is transported, and the effect of excitation is taken
into account in the effective rate that was used for both codes. This was
$\langle v_{e}\sigma\rangle_{iz}=2.147\times 10^{-14}$ m3/s. Figure 3 compares
the Gkeyll and DEGAS2 results for the cases including both ionization and
charge exchange. Steady-state neutral density profiles are shown in figures 3a
and 3c, demonstrating a slightly larger penetration length for neutrals in
Gkeyll. Neutral temperature profiles are compared in figures 3b and 3d.
Relative to DEGAS2, Gkeyllprofiles differ by up to 40% near the temperature
peaks and less than 10% elsewhere. Overall agreement is better for the higher
density case, seen in Figs. 3c and 3d. In both DEGAS2 and Gkeyll, the neutral
temperature peaks around 80 eV at the radial location where the density
approaches a minimum value. This demonstrates that at high densities neutrals
cannot penetrate far before becoming ionized, except for the more energetic
particles. Since both codes are kinetic, the neutrals can attain maximum
temperature values above the mean temperatures of the static ions and
electrons with which they interact (60 eV and 30 eV, respectively). In the
center of the $z$-domain, the neutrals are born from the source floor with a
temperature of 10 eV. Most are ionized before they can equilibrate to the 60
eV ion temperature via the charge exchange process; hence the temperature is
close to 40 eV.

For charge exchange, DEGAS2 uses the integral collision operator (Eq. 13) with
total cross sections from Refs. 52; 53, while Gkeyll uses a fitting function
to approximate the cross section that depends only on moments of the
distribution function. When only ionization is included and both codes use the
same ionization rate, near perfect agreement results as seen in Fig. 4,
indicating that differences in Fig. 3 arise from the charge exchange model.
Without the charge exchange process, whereby neutrals gain energy from the
ions, the neutral temperature is much lower, the peaks due to the penetration
of high energy neutrals is again visible at the point where the density
profiles approach the minimum value.

The comparison for tests including only charge exchange is shown in Fig. 5. In
5a the Gkeyll penetration length is noticeably greater than that in DEGAS2 and
the temperature peak is broader. Relative to DEGAS2, Gkeyllprofiles differ by
approximately 10% near the temperature peaks. Agreement is better in the high-
density comparison, shown in Fig. 5c-d, since the DEGAS2 charge exchange
operator converges towards a statistical average in the limit of high
collisionality. The Gkeyll model of charge exchange is most appropriate when
neutral densities are greater than $10^{17}$ m-3, since this is where profiles
approach DEGAS2 values. This comparison highlights the importance of an
integral collision operator for the charge exchange model. Including this in
Gkeyll is computationally feasible, particularly if a lower order form of the
distribution function is assumed, i.e. Maxwellian. This will be addressed in
future development work.

Figure 5: A comparison of benchmark tests from Gkeyll and DEGAS with static
plasma species that include only charge exchange and a recycling source at the
walls. Steady-state profiles of neutral density (a) and temperature (b) are
shown from the low plasma density case. Neutral density (c) and temperature
(d) are also compared from the high plasma density case.

##  IV Scrape-off layer simulations using the coupled plasma-neutral model

The neutral model, which evolves neutrals in three physical space plus three
velocity space dimensions (3X3V) has been tested in a coupling with
gyrokinetic simulations of the NSTX SOL, in three physical space dimensions
and two velocity space dimensions (3X2V). These have been extended from
previous Gkeyll NSTX simulationsShi  _et al._ (2019); E. L. Shi (2017);
Mandell  _et al._ (2020); Mandell (2021) and compared to a baseline case
without neutrals. In this section we describe the setup and results of the
simulations.

###  IV.1 Simulation setup

In the baseline case, which is labeled as “no neutrals” in the plots that
follow, electron and ion species were evolved with the electrostatic
gyrokinetic model, Eqs. 1-3, in three configuration space dimensions and two
velocity space dimensions (3X2V). We assume a simplified helical magnetic
geometry, with only open field lines and constant field line length and
curvature throughout the domain, which is described in more detail in Appendix
B. We approximate NSTX SOL conditions for an L-mode discharge, assuming a
singly-ionized deuterium plasma and power into the SOL $P_{SOL}=1.35$ MW.
Relevant parameters for the magnetic geometry include the magnetic field on
axis $B_{axis}=0.5$ T, the major radius at the magnetic axis $R_{0}=0.85$ m,
and the minor axis $a_{0}=0.5$ m. Plasma temperatures $T_{e0}=T_{i0}=40$ eV
are used to set the velocity grid extents and initial conditions.

Figure 6: Plot of the density particle source (m-3/s) for NSTX SOL
simulations, shown in the $(x,z)$ plane. The dashed white line indicates the
“quasi-separatrix” in our open-field line simulations. The dashed orange and
blue line indicates the conducting sheath and recycling boundary conditions
(BCs) applied in $z$ for plasma and neutral species, respectively. In the
radial direction, Dirichlet BCs with $\phi=0$ are applied for plasma species.
For neutral species, reflecting BCs are applied at the outer radial boundary
($x_{\max}$) and absorbing BCs are applied at the inner radial boundary
($x_{\min})$.

Perpendicular dimensions of the simulation box are $L_{x}=50\rho_{s0}\approx
14.6$ cm and $L_{y}=100\rho_{s0}\approx 29.1$ cm, where
$\rho_{s0}=c_{s0}/\Omega_{ci}$, and the box is centered at
$R_{c}=R_{0}+a_{0}$. The parallel length of the box is
$L_{z}=L_{pol}/\sin\chi=8$ m, where $L_{pol}=2.4$ m and
$\chi=\sin^{-1}(B_{v}/B)$ is the pitch angle determined by the ratio of the
vertical field magnitude $B_{v}$ to the total field magnitude $B$. Since the
flux tube simulation box covers a fraction of the SOL, the source is scaled as
$P_{src}=P_{SOL}L_{y}L_{z}/(2\pi RL_{pol})=0.31$ MW. A uniform grid is used
for all species. Resolution for the gyrokinetic species is ($N_{x}$, $N_{y}$,
$N_{z}$, $N_{v_{\parallel}}$, $N_{\mu}$) = (16, 32, 32, 12, 6). Convergence
was checked in $k_{x}$ space, and the resolution is appropriate. The
configuration space extents are $x\in[R_{c}-L_{x}/2,R_{c}+L_{x}/2]$,
$y\in[-L_{y}/2,L_{y}/2]$, and $z\in[-L_{z}/2,L_{z}/2]$. Velocity space extents
are $v_{\parallel}\in[-4v_{t,s0},4v_{t,s0}]$ and $\mu\in[0,6T_{s0}/B_{0}]$,
where $v_{t,s0}=\sqrt{T_{s0}/m_{s}}$ and $B_{0}=B_{axis}R_{0}/R_{c}$. In the
discontinuous Galerkin scheme, solutions are represented using piecewise-
linear $(p=1)$ basis functions on each cell, which effectively doubles the
resolution of the data in each direction. Dirichlet boundary conditions with
$\phi=0$ are applied in the radial direction and periodic boundary conditions
are used in the binormal direction. Conducting-sheath boundary conditionsShi,
Hakim, and Hammett (2015); Shi  _et al._ (2017) are applied in the parallel
direction, which reflect low-energy electrons back into the domain. Ions free
stream out of the domain along $z$. The location of these boundary conditions
is shown in Fig. 6.

Plasma species are sourced with non-drifting Maxwellian distribution functions
to model a particle and heat source at the midplane. The particle density
source is defined as

$S_{n}(x,z)=\left\\{\begin{array}[]{ll}S_{0}\max\left[\exp\left(\frac{-(x-x_{S})^{2}}{2\lambda_{S}}\right),0.01\right],&|z|<L_{z}/4\\\
0&{\rm else},\end{array}\right.$ (27)

where $x_{S}=R_{c}-0.05$ m and $\lambda_{S}=0.005$ m and is shown in Fig. 6 in
the $(x,z)$ plane. The particle source rate $S_{0}$ is calculated for a given
source temperature to achieve the desired power $P_{src}/2$, since power is
divided equally between ions and electrons. The temperature source was
increased relative to previous NSTX Gkeyll simulations to account for heat
loss due to ionization and charge exchange interactions. Due to discretization
the actual source temperature ended up being $T_{e,i}\approx 150$ eV for
$x<x_{S}+4\lambda_{S}$ and $T_{e,i}\approx 130$ eV for $x\geq
x_{S}+4\lambda_{S}$. The plasma species were initialized to a Maxwellian, with
density defined similarly to the source density function, Eq. 27 and a
temperature of $T_{e,i}=50$ eV for $x<x_{S}+3\lambda_{S}$ and $T_{e,i}=20$ eV
for $x\geq x_{S}+3\lambda_{S}$. In Fig. 6 the “quasi-separatrix” is denoted by
the white dashed line, which separates the source region from the SOL region.
Note that this is fixed at a radial value of $x\approx 1.32$ m along the field
line. This is lower than the typical location of the separatrix at the outer
midplane in the experiment ($x\approx$ 1.45 to 1.5 m).Stotler  _et al._
(2015); Scotti  _et al._ (2021)

The Lenard-Bernstein, or Dougherty, collision operator was used to model
plasma collisionsFrancisquez  _et al._ (2020) that include same and cross-
species collisions. A Spitzer collision frequency is calculated from the user-
defined background density $n_{0}=7\times 10^{18}$ m-3 and temperature
$T_{e0}=T_{i0}=40$ eV and is constant in space and time. The collision
frequency was reduced by a factor of 0.01 to reduce simulation run time, since
the timestep is largely dependent on the plasma collisionality. Reduced pitch-
angle scattering reduces the amount of high-energy electrons that enter the
sheath and could thus increase the sheath potential. Reduced resistivity can
also influence blob dynamics. It pushes the simulation towards the sheath-
limited regime, in which currents in the blobs are closed through parallel
currents to the sheath.Krasheninnikov, D’ippolito, and Myra (2008) The
collision frequency can be calculated dynamically based on local density and
temperatures, but this feature was not employed for these simulations. The
background density $n_{0}$ is also used for the ion guiding-center density
$n_{i0}^{g}$ that appears in the ion polarization term of the gyrokinetic
Poisson equation 3.

The case with neutrals is identical to the baseline case described above but
with the addition of neutrals species evolved with the Vlasov model, Eq. 4.
Ionization and charge exchange interactions are included. The neutral grid
resolution is ($N_{x}$, $N_{y}$, $N_{z}$, $N_{v_{x}}$, $N_{v_{y}}$,
$N_{v_{z}}$) = (16, 32, 32, 6, 6, 6), with the same configuration space
extents as the gyrokinetic species and velocity space extents
$v_{x,y,z}\in[-4v_{t,i0},4v_{t,i0}]$ defined in terms of the ion thermal
speed. Piecewise-linear basis functions are also used for the numerical
solution of the neutral distribution function.

Boundary conditions for neutrals are displayed in Fig. 6. In the radial
direction, absorbing boundary conditions are applied at the inner boundary to
model neutrals that penetrate into the core, and perfectly reflecting boundary
conditions are applied at the outer boundary. Since there is no ion flux at
the outer radial boundary in our model, the recycling boundary conditions are
not applied at the outer boundary. Periodic boundary conditions are applied in
the binormal direction, and recycling boundary conditions are applied in the
parallel direction. A recycling rate of $\alpha_{rec}=0.95$ and temperature
$T_{rec}=10$ eV are assumed.

A source floor is applied via a non-drifting Maxwellian distribution function.
A spatially constant particle source rate ${S_{n}=8\times 10^{-21}\,{\rm
m}^{3}{\rm s}^{-1}\,n_{i0}^{2}}$ approximates recombination, since we do not
yet include a model for this process in the code. Neutrals are also
initialized with a Maxwellian, with density profile that decays exponentially
away from the $z$-boundaries to approximate steady-state profiles and a
spatially constant temperature of 10 eV.

A discontinuous Galerkin algorithm that guarantees positivity of the
distribution function and maintains conservation properties has been developed
for the gyrokinetic solverMandell (2021) in Gkeyll but is not yet available
for the Vlasov solver. When recycling boundary conditions are used and no
positivity correction is applied, steep gradients develop in the neutral
distribution function in $z$ as neutrals enter the domain, which results in
unphysical, negative values for the neutral density and temperature.
Therefore, a non-conservative positivity correction is currently employed for
neutrals, which sets the distribution function to zero in a cell where the
cell average of the numerical solution has become negative. This can result in
an increase of approximately 10–20% of neutral particles and energy relative
to conservative simulations without the positivity correction. However, the
unphysical negative values of density and temperature present in the latter
make them unreliable baselines. A conservative positivity-preserving algorithm
for the Vlasov solver is under development and will be included in future
work.

In summary, the nonphysical assumptions in these proof-of-concept simulations
include a collision frequency based on constant values of density and
temperature, which is scaled by 0.01; a temperature of 10 eV for recycled
neutrals; and the non-conservative positivity correction for the neutral
distribution function.

###  IV.2 Results

The simulations were run to 0.4 ms, which is about 4 ion transit times,
$\tau_{i}=L_{z}/(2c_{s0})$. Density profiles were still evolving since the
sources and sinks were not yet balanced. Changes to the profiles are less than
$\lesssim 10\%$ on relevant timescales ($\sim 10$ $\mu$s) based on neutral
collisions and turbulence correlation times. Coulomb collisions were not
considered due to the artificially reduced collisionality, which was described
in the previous section. On the Cori Haswell nodes at NERSC, the case without
neutrals required approximately 18.5k CPU-hours and the case with neutrals
required approximately 93.4k CPU-hours, or about 5 times as many resources.

Snapshots of midplane densities are shown in Fig. 7, with the electron density
from the case without neutrals shown on the left, electron density from the
case with neutrals shown center, and the neutral density shown on the right.
The electron density without neutrals displays a sharper density gradient,
while the case with neutrals displays a density that is more constant
radially. The magntidues are close to NSTX experimental values for L-mode
discharges (see Fig. 7 in Ref. 57). The neutral density is five orders of
magnitude less than the plasma density. Since we do not include recycling
boundary conditions at the outer radial wall (see Fig. 6), neutrals are only
sourced at the endplates and most are ionized prior to reaching the midplane.
Because of this simulation neutral values at the midplane are significantly
less than in experiment.Stotler  _et al._ (2015); Scotti  _et al._ (2021) The
neutral profile is approximately constant radially in the SOL region, whereas
the neutral profile decays exponentially from the outer wall to lower radii.
In the simulation, neutral density is artificially high on the inner radial
boundary due the conducting wall ($\phi=0$) radial boundary conditions used
for the plasma species. Because of this, the plasma density drops at the inner
radial boundary, which results in lower ionization rates. When closed flux
surface geometry is eventually included, plasma density will be higher inside
the separatrix, and that feature should vanish.

Including recycling at the outer radial boundary will likely improve
comparison with experiment. Recycling is included at the $z$-boundaries and
neutral densities at these boundaries, shown by the profile in $z$ in Fig. 8,
approach experimental values.Stotler  _et al._ (2015); Scotti  _et al._ (2021)

Figure 7: Snapshots of density taken perpendicular to the magnetic field at
the midplane. In all plots, a dashed white line represents a “quasi-
separatrix” that separates the source region from the SOL region. The left
plot depicts electron density from the simulation without neutrals. The center
plot depicts electron density from the simulation with neutrals, and the right
plot shows the neutral density. Figure 8: Neutral density profile along $z$,
averaged in $x$, $y$ and time. Figure 9: Steady-state density profiles are
compared. In both plots the black dashed line represents the “quasi-
separatrix” and the gray region denotes the location of the source. Plot (a)
compares the electron densities from the simulations without and with neutrals
on a linear scale. Plot (b) compares electron densities normalized to their
maximum density on a log scale to highlight the flatter density profile in the
case with neutrals. Figure 10: Radial turbulent fluxes, which are calculated
as the product of density fluctuations and radial $E\times B$ velocity
fluctuations and normalized to the mean density, are compared in (a). Electron
parallel particle fluxes, which are calculated from the electron distribution
functions at the $z$ boundary and normalized to mean density, are compared in
(b). The ratio of the fluxes are compared in (c).

Steady-state density profiles were obtained by averaging fields from 0.3 to
0.4 ms. Electron density profiles at the midplane are compared in Fig. 9. The
source region is shaded gray and the “quasi-separatrix” is denoted by the
vertical dashed line. In the simulation with neutrals, ionization increases
density by about a factor of 3, as seen in Fig. 9a. We also observe a
flattening of the density profile to the right of the peak for the case with
neutrals, which is seen most clearly in Fig. 9b, where we have normalized each
density profile to its maximum value and then plotted them on a log scale.
This appears similar to the density shoulder that has been observed
experimentally Boedo  _et al._ (2003); Rudakov  _et al._ (2005); Carralero
_et al._ (2017); Vianello  _et al._ (2017, 2019) and also studied in recent
GBS simulations that include neutral interactions.Mancini  _et al._ (2021)
However, additional Gkeyll simulations with neutrals and scans in the density
source are necessary to make direct comparisons to those results.

In Figure 10, we compare the perpendicular and parallel particle fluxes,
normalized to the density profiles (from Fig. 9a). The radial perpendicular
flux is due mostly to turbulence and is thus calculated from the density
fluctuations and the radial $E\times B$ velocity fluctuations,
$\Gamma_{\perp}=\langle\tilde{n}_{e}\tilde{v}_{r}\rangle_{y,t}$, where
fluctuating quantities are given by $\tilde{A}=A-\langle A\rangle_{y,t}$. The
radial velocity is given by $v_{r}=-(1/B)\partial{\phi}/\partial{y}$. In Fig.
10a, the normalized perpendicular flux is reduced in the simulation with
neutrals. The peak is shifted to higher radii, and the gradient of the flux
profile is shallower. The normalized parallel flux is compared in Fig. 10b and
is calculated as

$\Gamma_{\parallel,e}=\left\langle\left.\int\mathcal{J}f_{e\,h}\dot{\bm{R}}_{h}\cdot\hat{\bm{b}}\,{\rm
d}x\,{\rm d}y\,{\rm d}^{3}\,\bm{v}\right|_{z=L_{z}/2}\right\rangle.$ (28)

The subscript $h$ on denotes that only electrons with kinetic energy large
enough to overcome the potential barrier at the pre-sheath are included in
this calculation. This is consistent with the conducting sheath boundary
conditions described in the previous section. The normalized parallel flux to
the endplates decreases for the case with neutrals; the parallel outflow is
slower since the electron temperature has decreased. In Fig. 10c, the ratios
of the radial to the parallel fluxes are compared. This ratio for the neutrals
case is generally larger than the baseline simulation, suggesting that the
relative decrease in parallel flow leads to density profile flattening,
similar to GBS simulations with neutrals.Mancini  _et al._ (2021)

Figure 11: Steady-state temperature profiles of electron and ion temperature
are compared in (a). Pressure profiles are compared in (b).

Midplane steady-state temperature profiles are shown in Fig. 11a. In the
simulation with neutrals, electron and ion temperatures are both reduced. The
ionization energy requirement reduces the electron temperature. Ion
temperature decreases due to charge exchange and interaction with electrons by
cross-species collisions. The electron temperature gradient is reduced in the
case with neutrals, and the electron and ion temperature profiles have more
similar gradients than in the case without neutrals, indicating a stronger
coupling between the two. The pressure profiles, where $p=n_{e}(T_{e}+T_{i})$,
are shown in Fig. 11b. The case with neutrals is larger due to the increase in
density.

Figure 12: Parallel heat fluxes to the endplates are compared. These profiles
are sums of electron and ion contributions, which were calculated from the
distribution functions using Eq. 29. Heat flux widths were estimated from
exponential fits, $\exp(-x/\lambda_{q})$, at the “quasi-separatrix”, shown by
the green and red dashed curves.

The profiles of parallel heat flux to the endplates are compared in Fig. 12.
The heat flux for each species is calculated as

$Q_{\parallel,s}^{\rm end}=\left\langle\left.\int
H_{s\,h}\mathcal{J}f_{s\,h}\dot{\bm{R}}_{h}\cdot\hat{\bm{b}}\,{\rm d}x\,{\rm
d}y\,{\rm d}^{3}\,\bm{v}\right|_{z=L_{z}/2}\right\rangle,$ (29)

and the total $Q_{\parallel,\rm tot}^{\rm end}=Q_{\parallel,e}^{\rm
end}+Q_{\parallel,i}^{\rm end}$ is plotted. We estimate the width of the
profiles with an exponential fit just outside the quasi-separatrix,
$Q_{\parallel}^{\rm end}(x=1.32)e^{-x/\lambda_{q}}$, and the value for the
case with neutrals is slightly less than the case without neutrals.

The power balance in each simulation is also calculated and shown in Fig. 13.
Both have the same input power at the midplane, indicated by the red line. In
the case without neutrals, the only loss is at the endplates, indicated by the
blue line. The sum of of the input power and loss to endplates is given by the
orange line, which approaches zero. This conforms to expectations for Gkeyll’s
conservative algorithms. In the case with neutrals, some power is lost to
ionization and charge exchange interactions, indicated by the brown line.
Since the input power is the same, the magnitude of losses to the endplates is
slightly reduced. The sum of the source and losses for the neutrals case is
given by the purple line, which also approaches zero, albeit with slightly
worse errors. Presumably this is due to slightly larger fluctuations in the
stored energy, arising from neutral interactions.

Figure 13: The power balance is calculated for the simulations without and
with neutrals. Both have the same input power $P_{SOL}$ indicated by the red
curve. The only power loss in the simulation without neutrals is due to heat
flux out through the endplates (blue curve). The sum of these gives the orange
curve. The simulation with neutrals also includes some power loss due to
ionization and charge exchange (brown curve). Power loss to the endplates is
given by the green curve and with the power input and other losses this sums
to the purple curve. Figure 14: The linear interchange growth rates are
calculated from the pressure profiles and compared in (a), where the
simulation with neutrals shows reduced levels. The normalized density
fluctuations are compared in (b) and are reduced in the simulation with
neutrals.

Since the turbulence is primarily driven by bad curvature, we have calculated
the linear interchange growth rate $\gamma_{I}=c_{s}\sqrt{2/(L_{p}x)}$, where
$c_{s}=\sqrt{(T_{e}+T_{i})/m_{i}}$ is the local sound speed and
$L_{p}=\max\left(0,-\left(\frac{1}{p}\frac{dp}{dx}\right)^{-1}\right)$ is the
pressure gradient scale length. The growth rates are compared in Fig. 14a, and
the curvature drive is reduced in the case with neutrals due to the flattening
of the pressure profile. There is a corresponding reduction in normalized rms
density fluctuation values, shown in Fig. 14b.

Turbulence statistics at the midplane, such as the probability distribution
functions (PDFs) of turbulence fluctuations, have also been investigated. For
the PDF, density fluctuations are measured over the time interval 200–400
$\mu$s at a particular radial value and normalized to the rms value of
fluctuations ($\tilde{n}/\tilde{n}_{rms}$. Examples of PDFs calculated at
$x=1.37$ m are shown in Fig. 15 and compared to a Gaussian distribution
(dashed curve). The case without neutrals is more skewed in the positive
direction and has a longer positive tail. Positive skewness and excess
kurtosis of the density fluctuation PDFs are characteristics of blobby
turbulence.D’Ippolito, Myra, and Zweben (2011) Skewness and excess kurtosis
are calculated as $E[\tilde{n}^{3}]/\sigma^{3}$ and
$E[\tilde{n}^{4}]/\sigma^{4}-3$, respectively, where $E[...]$ is the
expectation value, given by the PDF. The skewness and excess kurtosis of
density fluctuation PDFs are shown in Fig. 16a and Fig. 16b, respectively. The
skewness and kurtosis for the case with neutrals are nearly constant in $x$.
This is consistent with the weaker density gradients and lower fluctuations
observed previously, as there are fewer turbulence events that deviate greatly
from the mean. For example, in the plot of density fluctuation profiles (Fig.
14b), the fluctuation levels from the case with neutrals are nearly constant
in $x$ while they increase more noticeably at large radii for the case without
neutrals case. We note that the increase in skewness and kurtosis at large
radii in the case without neutrals is more consistent with experimental
observations.D’Ippolito, Myra, and Zweben (2011)

Figure 15: Probability distribution functions (PDFs) of electron density
turbulence fluctuations at radial location $x=1.37$ m are compared to Gaussian
distribution (dashed curve). Figure 16: Skewness (a) and excess kurtosis (b)
of electron density fluctuation PDFs are compared.

Radial correlation lengths are compared in Fig. 17a. We find no change in the
radial correlation lengths, except at large radii, where they decrease for the
case with neutrals. This can be used as a proxy to determine blob size.
Autocorrelation times, $\tau_{ac}$, are compared in Fig. 17b, which are
calculated as the $e$-folding time of the autocorrelation functions. The
autocorrelation function is calculated as
$C(\tau)=\langle\tilde{n}(t)\tilde{n}(t+\tau)\rangle/\langle\tilde{n}(t)^{2}\rangle$,
where brackets denote a time average. Autocorrelation times increase in the
simulation with neutrals, suggesting longer temporal coherency of turbulent
structures. Both radial correlation lengths and autocorrelation times were
calculated with data from the interval 200–400 $\mu$s.

Figure 17: Comparisons of radial correlation lengths (a) and autocorrelation
times (b).

The proof-of-concept simulation with neutrals incorporated a self-consistent
source via ionization. Neutral interactions impacted the steady-state profiles
and turbulence statistics in noticeable ways, demonstrating the importance of
incorporating them in predictive modeling. Additional studies have being
carried out in which an additional volumetric source to approximate the
ionization source is added in simulations without the neutral species, in
order to isolate sourcing effects from charge exchange collisions. These
demonstrated that most of the differences in the comparison presented herein
were due to sourcing. However, differences in properties of the blob dynamics
remained. The full details of these results are outside the scope of this
paper and will be presented in a future publication.

##  V Conclusions

We have presented the coupling of a continuum gyrokinetic plasma model to a
continuum kinetic model for neutrals within Gkeyll. A continuum code has the
benefit of reduced noise and improved accuracy at a fixed resolution.
Simplified models of electron impact ionization, charge exchange and wall
recycling have been included. These models have been verified against analytic
theory and show good agreement with DEGAS2 in benchmark tests.

A proof-of-concept has been demonstrated in a high dimensional simulation of
the NSTX SOL with simplified helical magnetic geometry. This simulation
included gyrokinetic electron and ion species in three spatial dimensions and
two velocity space dimensions coupled to Vlasov neutrals in three spatial
dimensions and three velocity space dimensions. In a comparison to a baseline
simulation without neutrals and a midplane source identical to the simulation
with neutrals, the latter exhibited a flatter density profile, reduced
turbulence fluctuations, and longer auto-correlation times.

The flatter density profile contributed to a flattening of the pressure
profile and a corresponding reduction in the linear interchange growth rate.
The normalized density fluctuations were also reduced. In Fig. 16 skewness and
kurtosis were reduced in the simulation with neutrals, also likely due to the
increased density and the flatter density profile on the low-field side of the
domain. The sharper density gradient in the case without neutrals means that
high-density coherent turbulent structures that originate from the source
region deviate more from the background density as they move radially outward,
leading to a stronger positive skewness than in the simulation with neutrals.
Large values of skewness and excess kurtosis are associated with greater
intermittency, i.e. blobby turbulence.D’Ippolito, Myra, and Zweben (2011) The
increased correlation times for the case with neutrals, seen in Fig. 17b,
suggest turbulent structures are more coherent, which may be affecting the
nature of the turbulence and the resulting radial transport of particles..

Collisionality and magnetic geometry also have important effects on SOL
turbulence,Myra, Russell, and D’Ippolito (2006) and subsequent simulations
will include these features more realistically. A reduced Spitzer collision
frequency that was constant in time and space was used in these simulations to
reduce computational costs. The collision frequency can be calculated
dynamically in Gkeyll based on local densities and temperatures, and this
feature will be included in simulations with neutrals. Work is ongoing to
extend the neutral model to be compatible with the general geometry features
available in Gkeyll, which are outlined in Chapter 5 of Ref. 42. Chapter 4 of
the same reference also demonstrates the impact of electromagnetic effects on
blob dynamics, and this could also be studied in simulations that include
neutral interactions.

Furthermore, the effect of neutral interactions in high-recycling or detached
SOL scenarios is important. Modeling these scenarios with Gkeyll will require
the ability to represent low neutral temperatures on the velocity space grid.
A non-uniform grid is currently under development and will help with this
issue. Models for recombination and radiation will eventually be included to
model detached conditions. This will enable us to incorporate new and existing
features such as increased collisionality, realistic geometry and
electromagnetic plasma species in simulations with neutrals.

Some of the benefits of Monte Carlo neutral codes are challenging to reproduce
with a continuum code. MC codes can simulate many species simultaneously,
including molecules and excited metastable states. They can also easily model
short-lived species or those that only occupy part of the computational
domain. On the other hand, continuum codes have the advantage of being able to
model scenarios with a large density variation whose dynamical evolution
significantly impacts physical outcomes. The Gkeyll model can self-
consistently explore the effect of neutral interactions on plasma profiles and
turbulence statistics, while most MC codes couple to reduced models of plasma
transport and do not consider effects on turbulence. With forthcoming GPU
capabilities, Gkeyll will be able to model plasma species, a neutral species,
and an impurity species at reasonable computational cost, and results could be
used to improve reduced models of plasma transport in MC codes.

##  Appendix A Getting Gkeyll and reproducing results

Readers may reproduce our results and also use Gkeyll for their applications.
The code and input files used here are available online. Full installation
instructions for Gkeyll are provided on the Gkeyll website.Gkeyll (2022) The
code can be installed on Unix-like operating systems (including Mac OS and
Windows using the Windows Subsystem for Linux) either by installing the pre-
built binaries using the conda package manager (https://www.anaconda.com) or
building the code via sources. The input files used here are under version
control and can be obtained from the repository at
https://github.com/ammarhakim/gkyl-paper-inp/tree/master/2022_PoP_neutrals.

##  Appendix B Neutrals in simplified helical geometry

We now describe some subtleties of the coordinate system employed in Gkeyll
and show how the Vlasov neutral phase-space grid used in the NSTX simulations
is consistent with various geometric assumptions. In this work and previous
gyrokinetic simulations with the Gkeyll code,E. L. Shi (2017); Shi  _et al._
(2019); Bernard  _et al._ (2019, 2020); Mandell  _et al._ (2020) we have used
a simplified helical model of the magnetic field. The mapping from the
cylindrical coordinates $(R,\varphi,Z)$ to those used in Gkeyll $(x,y,z)$ is
described in Refs. 31; 37 and Chapter 4 of Ref. 42 as:

$x=R,\hskip 14.22636ptz=\frac{Z}{\sin\vartheta}=\frac{L_{c}}{H}Z,\hskip
14.22636pty={R_{c}}{\sin\vartheta}\left(\varphi-\frac{Z}{R_{c}}\cot\vartheta\right),$
(30)

where $R_{c}$ is the center of the $x$-domain, $H$ is the vertical height,
$L_{c}$ is the connection length, $\vartheta$ is the field line pitch angle,
and $\sin\vartheta=B_{v}/B=H/L_{c}$. Nearly all geometrical factors arising
from this non-orthogonal coordinate system are neglected in this simplified
helical model, which effectively takes the $B_{v}\ll B$ limit. Only the
assumption that perpendicular gradients of the distribution function are much
stronger than parallel gradients is retained. By neglecting such geometric
factors, simulations in the non-orthogonal helical field coordinate system
defined by Eq. 30 can also be interpreted as occurring in a simpler orthogonal
coordinate system with only a toroidal magnetic field, as given by the
mapping:

$x=R,\hskip 28.45274pty=Z,\hskip 28.45274ptz=R_{c}\varphi.$ (31)

Figure 18: Diagram of the magnetic geometry corresponding to the geometry
assumptions made in the simulations presented in Sec. IV. (Image from N.
Mandell, PhD Thesis, Princeton University, 2021; licensed under a Creative
Commons Attribution (CC BY) license.)

This is depicted in Fig. 18 which originally appeared in Ref. 42. The
relationship between these two geometry interpretations was pointed out in
Appendices 5.A–5.B of Ref. 42. Note that both coordinate systems are field-
aligned, since gradients parallel to the magnetic field are expressed in terms
of derivatives involving only one coordinate,
$\nabla_{||}=(1/|\vec{B}|)\vec{B}\cdot\nabla\propto\partial/\partial z$. The
difference is that the direction of the magnetic field in physical space is
slightly different in the two coordinate systems. In the geometry of Eq. 31,
the field-line coordinate $z$ is purely toroidal and the binormal coordinate
$y$ is in the vertical direction. The sheath conducting boundary conditions
along the $z$ direction are then effectively being applied along a limiter
that extends in the vertical direction. The simulated field lines start at a
toroidally-localized limiter, but may wrap more (or less) than $2\pi R$ around
the torus before hitting a limiter again. The gradient basis vectors are

$\nabla x=\hat{\bm{R}},\hskip 28.45274pt\nabla y=\hat{\bm{Z}},\hskip
28.45274pt\nabla z=\frac{R_{c}}{x}\hat{\bm{\varphi}},$ (32)

The present method of coupling the Vlasov neutral solver to the gyrokinetic
solver is consistent with this configuration. The neutral species share the
same field-line-following coordinate system for configuration space $(x,y,z)$,
and we assume an orthogonal cylindrical coordinate system for the velocity
space, $(v_{R},v_{Z},v_{\varphi})$. The streaming term in the Vlasov equation
(4) becomes

$\displaystyle\bm{v}\cdot\nabla f$ $\displaystyle=$
$\displaystyle\bm{v}\cdot\left(\nabla x\frac{\partial f}{\partial x}+\nabla
y\frac{\partial f}{\partial y}+\nabla z\frac{\partial f}{\partial z}\right),$
(33) $\displaystyle=$ $\displaystyle v_{R}\frac{\partial f}{\partial
x}+v_{Z}\frac{\partial f}{\partial y}+v_{\varphi}\frac{R_{c}}{x}\frac{\partial
f}{\partial z},$ $\displaystyle\approx$ $\displaystyle v_{R}\frac{\partial
f}{\partial x}+v_{Z}\frac{\partial f}{\partial y}+v_{\varphi}\frac{\partial
f}{\partial z},$

where we have assumed $R_{c}/x\approx 1$ for a flux-tube-like domain localized
around $x=R_{c}$. This configuration is sufficient for this paper in showing
the overall feasibility of this approach and doing a first study of neutral-
turbulence interaction with this code. A more realistic helical geometry with
magnetic shear is available in Gkeyll and has been described in Chapter 5 of
Ref. 42. Eventually, the neutral model will be incorporated in simulations
with more realistic geometry and the streaming term in the Vlasov equation
will be updated to include additional geometric terms that arise from the
associated gradient basis vectors.

##  Appendix C Derivation of 1X1V analytic solution for verification tests

To verify models of neutral interactions, we derive theoretical predictions of
steady-state profiles. First consider the 1X1V neutral equation in the form

$\frac{\partial f}{\partial t}+v_{z}\frac{\partial f}{\partial
z}=-\nu\left(f(z,v_{z})-g(z,v_{z})\right),$ (34)

where

$g(z,v_{z})=\frac{f_{i}}{n_{i}}\int dv_{z}^{\prime}f(z,v_{z}^{\prime}).$ (35)

In the case with only ionization, $f_{i}=0$. We seek a steady-state solution:

$v_{z}\frac{\partial f}{\partial z}=-\nu\left(f(z,v_{z})-g(z,v_{z})\right).$
(36)

Consider a finite, symmetric domain, $z=[-L/2,L/2]$, when solving Eq. 36. We
now have solutions for the left-hand and right-hand sides of the domain:

$f(z,v_{z})=\left\\{\begin{array}[]{ll}C_{+}(z,v_{z})e^{-\frac{\nu}{v_{z}}(z+L/2)}&v_{z}>0\\\
C_{-}(z,v_{z})e^{-\frac{\nu}{v_{z}}(z-L/2)}&v_{z}<0.\end{array}\right.$ (37)

This becomes

$C_{+}(z,v_{z})=C_{+}(-L/2,v_{z})+\int_{-L/2}^{z}dz^{\prime}\frac{\nu}{v_{z}}g(z^{\prime},v_{z})e^{\frac{\nu}{v_{z}}(z^{\prime}+L/2)}$
(38)
$C_{-}(z,v_{z})=C_{-}(L/2,v_{z})-\int_{z}^{L/2}dz^{\prime}\frac{\nu}{v_{z}}g(z^{\prime},v_{z})e^{\frac{\nu}{v_{z}}(z^{\prime}-L/2)},$
(39)

which gives

$f(z,v_{z})=\left\\{\begin{array}[]{ll}\left(C_{+}(-L/2,v_{z})+\int_{-L/2}^{z}dz^{\prime}\frac{\nu}{v_{z}}g(z^{\prime},v_{z})e^{\frac{\nu}{v_{z}}(z^{\prime}+L/2)}\right)e^{-\frac{\nu}{v_{z}}z}&v_{z}>0\\\
\left(C_{-}(L/2,v_{z})-\int_{z}^{L/2}dz^{\prime}\frac{\nu}{v_{z}}g(z^{\prime},v_{z})e^{\frac{\nu}{v_{z}}(z^{\prime}-L/2)}\right)e^{-\frac{\nu}{v_{z}}(z-L)}&v_{z}<0.\end{array}\right.$
(40)

To simplify notation, we introduce $f(-L/2,v_{z}>0)=f^{BC}_{-L/2}(v_{z})$, and
$f(L/2,v_{z}<0)=f^{BC}_{L/2}(v_{z})$ which are the distributions numerically
imposed at the boundary conditions, such that

$f(z,v_{z})=\left\\{\begin{array}[]{ll}\left(f^{BC}_{-L/2}(v_{z})+\int_{-L/2}^{z}dz^{\prime}\frac{\nu}{v_{z}}g(z^{\prime},v_{z})e^{\frac{\nu}{v_{z}}(z^{\prime}+L/2)}\right)e^{-\frac{\nu}{v_{z}}z}&v_{z}>0\\\
\left(f^{BC}_{L/2}(v_{z})-\int_{z}^{L/2}dz^{\prime}\frac{\nu}{v_{z}}g(z^{\prime},v_{z})e^{\frac{\nu}{v_{z}}(z^{\prime}-L/2)}\right)e^{-\frac{\nu}{v_{z}}(z-L)}&v_{z}<0\end{array}\right.$
(41)

The density is thus given by

$n(z)=n^{BC}_{-L/2}(z)+\int_{-L/2}^{z}dz^{\prime}n(z^{\prime})\beta^{+}_{i}(z,z^{\prime})+n^{BC}_{L/2}(z)-\int_{z}^{L/2}dz^{\prime}n(z^{\prime})\beta^{-}_{i}(z,z^{\prime}),$
(42)

where

$\displaystyle n^{BC}_{-L/2}(z)$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}dv_{z}^{\prime}\,f^{BC}_{-L/2}(v_{z}^{\prime})e^{-\frac{\nu}{v_{z}^{\prime}}(z+L/2)}$
(43) $\displaystyle n^{BC}_{L/2}(z)$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{0}dv_{z}^{\prime}\,f^{BC}_{L/2}(v_{z}^{\prime})e^{-\frac{\nu}{v_{z}^{\prime}}(z-L/2)}$
(44) $\displaystyle\beta^{+}_{i}(z^{\prime},z)$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}dv_{z}^{\prime}\,\frac{\nu}{v_{z}^{\prime}}G_{i}(z^{\prime},v_{z}^{\prime})e^{\frac{\nu}{v_{z}^{\prime}}(z^{\prime}-z)}$
(45) $\displaystyle\beta^{-}_{i}(z^{\prime},z)$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{0}dv_{z}^{\prime}\,\frac{\nu}{v_{z}^{\prime}}G_{i}(z^{\prime},v_{z}^{\prime})e^{\frac{\nu}{v_{z}^{\prime}}(z^{\prime}-z)},$
(46)

and $G_{i}(z,v_{z})=f_{i}(z,v_{z})/n_{i}(z)$.

###### Acknowledgements.

T. Bernard acknowledges helpful discussions with J. Boedo and F. Scotti. This
material is based upon work supported by the U.S. Department of Energy, Office
of Science, Office of Fusion Energy Sciences, Theory Program, under Award No.
DE-FG02-95ER54309. J. Juno acknowledges support from a NSF Atmospheric and
Geospace Science Postdoctoral Fellowship (Grant No. AGS-2019828). J. Guterl
acknowledges DOE grant DE-SC0018423. Computing resources were provided by the
Extreme Science and Engineering Discovery Environment (XSEDE), which is
supported by National Science Foundation grant number ACI-1548562 and also by
the National Energy Research Scientific Computing Center (NERSC), which is an
Office of Science User Facility supported under Contract DE-AC02-05CH11231. –
Disclaimer – This report was prepared as an account of work sponsored by an
agency of the United States Government. Neither the United States Government
nor any agency thereof, nor any of their employees, makes any warranty,
express or implied, or assumes any legal liability or responsibility for the
accuracy, completeness, or usefulness of any information, apparatus, product,
or process disclosed, or represents that its use would not infringe privately
owned rights. Reference herein to any specific commercial product, process, or
service by trade name, trademark, manufacturer, or otherwise, does not
necessarily constitute or imply its endorsement, recommendation, or favoring
by the United States Government or any agency thereof. The views and opinions
of authors expressed herein do not necessarily state or reflect those of the
United States Government or any agency thereof.

## Data Availability Statement

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * Stangeby (2000) P. C. Stangeby,  _The Plasma Boundary of Magnetic Fusion Devices_ (Institute of Physics Publishing, 2000).
  * Boivin  _et al._ (2000) R. L. Boivin, J. A. Goetz, A. E. Hubbard, J. W. Hughes, I. H. Hutchinson, J. H. Irby, B. LaBombard, E. S. Marmar, D. Mossessian, C. S. Pitcher, J. L. Terry, B. A. Carreras, and L. W. Owen, “Effects of neutral particles on edge dynamics in Alcator C-Mod plasmas,” Physics of Plasmas 7, 1919–1926 (2000).
  * Groebner  _et al._ (2002) R. J. Groebner, M. A. Mahdavi, A. W. Leonard, T. H. Osborne, G. D. Porter, R. J. Colchin, and L. W. Owen, “The role of neutrals in high-mode (H-mode) pedestal formation,” Physics of Plasmas 9, 2134 (2002).
  * Moser  _et al._ (2020) A. L. Moser, L. Casali, B. M. Covele, A. W. Leonard, A. G. McLean, M. W. Shafer, H. Q. Wang, and J. G. Watkins, “Separating divertor closure effects on divertor detachment and pedestal shape in DIII-D,” Physics of Plasmas 27, 032506 (2020).
  * Mordijck (2020) S. Mordijck, “Overview of density pedestal structure: Role of fueling versus transport,” Nuclear Fusion 60, 082006 (2020).
  * Thrysøe  _et al._ (2018) A. S. Thrysøe, M. Løiten, J. Madsen, V. Naulin, A. H. Nielsen, and J. J. Rasmussen, “Plasma particle sources due to interactions with neutrals in a turbulent scrape-off layer of a toroidally confined plasma,” Physics of Plasmas 25 (2018), 10.1063/1.5019662.
  * Zholobenko  _et al._ (2021) Z. Zholobenko, A. Stegmeir, G. D. Griener, G. Conway, T. Body, D. Coster, F. Jenko, and ASDEX Upgrade Team, “The role of neutral gas in validated global edge turbulence simulations,” Nucl. Fusion 59, 112014 (2021).
  * Dekeyser  _et al._ (2019) W. Dekeyser, X. Bonnin, S. W. Lisgo, R. A. Pitts, and B. LaBombard, “Implementation of a 9-point stencil in solps-iter and implications for alcator c-mod divertor plasma simulations,” Nuclear Materials and Energy 18, 125–130 (2019).
  * D’Ippolito, Myra, and Zweben (2011) D. A. D’Ippolito, J. R. Myra, and S. J. Zweben, “Convective transport by intermittent blob-filaments: Comparison of theory and experiment,” Phys. Plasmas 18, 60501 (2011).
  * Reiter (2002) D. Reiter, “EIRENE—a Monte Carlo linear transport solver,” Available at www.eirene.de (2002).
  * Stotler and Karney (1994) D. Stotler and C. Karney, “Neutral gas transport modeling with DEGAS 2,” Contributions to Plasma Physics 34, 392–397 (1994), https://onlinelibrary.wiley.com/doi/pdf/10.1002/ctpp.2150340246 .
  * Cupini, De Matteis, and Simonini (1983) E. Cupini, A. De Matteis, and R. Simonini, “NIMBUS-Monte Carlo simulation of neutral particle transport in fusion devices,” Brussels: Commission of the European Communities EUR 12, 9 (1983).
  * Marandet  _et al._ (2013) Y. Marandet, P. Tamain, R. Futtersack, P. Ghendrih, H. Bufferand, P. Genesio, and A. Mekkaoui, “Influence of neutral particles on scrape-off layer turbulence with application to the interpretation of fast camera data,” Journal of Nuclear Materials 438 (2013), 10.1016/j.jnucmat.2013.01.107.
  * Galassi  _et al._ (2017) D. Galassi, P. Tamain, H. Bufferand, G. Ciraolo, P. Ghendrih, C. Baudoin, C. Colin, N. Fedorczak, N. Nace, and E. Serre, “Drive of parallel flows by turbulence and large-scale e × B transverse transport in divertor geometry,” Nuclear Fusion 57 (2017), 10.1088/1741-4326/aa5332.
  * Fan  _et al._ (2019) D. M. Fan, Y. Marandet, P. Tamain, H. Bufferand, G. Ciraolo, P. Ghendrih, and E. Serre, “Effect of turbulent fluctuations on neutral particles transport with the TOKAM3X-EIRENE turbulence code,” Nuclear Materials and Energy 18, 105–110 (2019).
  * Zhang  _et al._ (2019) D. R. Zhang, Y. P. Chen, X. Q. Xu, and T. Y. Xia, “Self-consistent simulation of transport and turbulence in tokamak edge plasma by coupling solps-iter and bout++,” Physics of Plasmas 26, 012508 (2019).
  * Stotler  _et al._ (2013) D. P. Stotler, C. S. Chang, G. Park, and S. H. Ku, “Energy conservation tests of a coupled kinetic plasma-kinetic neutral transport code,” Computational Science and Discovery 6 (2013), 10.1088/1749-4699/6/1/015006.
  * Stotler  _et al._ (2017) D. Stotler, J. Lang, C. Chang, R. Churchill, and S. Ku, “Neutral recycling effects on ITG turbulence,” Nuclear Fusion 57, 086028 (2017).
  * Ku  _et al._ (2018) S. Ku, C. S. Chang, R. Hager, R. M. Churchill, G. R. Tynan, I. Cziegler, M. Greenwald, J. Hughes, S. E. Parker, M. F. Adams, E. D’Azevedo, and P. Worley, “A fast low-to-high confinement mode bifurcation dynamics in the boundary-plasma gyrokinetic code XGC1,” Physics of Plasmas 25 (2018), 10.1063/1.5020792.
  * Hager  _et al._ (2019) R. Hager, C. S. Chang, N. M. Ferraro, and R. Nazikian, “Gyrokinetic study of collisional resonant magnetic perturbation (RMP)-driven plasma density and heat transport in tokamak edge plasma using a magnetohydrodynamic screened RMP field,” Nuclear Fusion 59 (2019), 10.1088/1741-4326/ab3be2.
  * Joseph  _et al._ (2017) I. Joseph, M. Rensink, D. Stotler, A. Dimits, L. LoDestro, G. Porter, T. Rognlien, B. Sjogreen, and M. Umansky, “On coupling fluid plasma and kinetic neutral physics models,” Nuclear Materials and Energy 12, 813 – 818 (2017), proceedings of the 22nd International Conference on Plasma Surface Interactions 2016, 22nd PSI.
  * Wersal and Ricci (2015) C. Wersal and P. Ricci, “A first-principles self-consistent model of plasma turbulence and kinetic neutral dynamics in the tokamak scrape-off layer,” Nuclear Fusion 55, 123014 (2015).
  * Wersal and Ricci (2017) C. Wersal and P. Ricci, “Impact of neutral density fluctuations on gas puff imaging diagnostics,” Nuclear Fusion 57 (2017), 10.1088/1741-4326/aa7db0.
  * Wersal, Ricci, and Loizu (2017) C. Wersal, P. Ricci, and J. Loizu, “A comparison between a refined two-point model for the limited tokamak SOL and self-consistent plasma turbulence simulations,” Plasma Physics and Controlled Fusion 59, 1–20 (2017).
  * Mandell  _et al._ (2020) N. Mandell, A. Hakim, G. Hammett, and M. Francisquez, “Electromagnetic full-gyrokinetics in the tokamak edge with discontinuous galerkin methods,” Journal of Plasma Physics 86 (2020).
  * Mancini  _et al._ (2021) D. Mancini, P. Ricci, N. Vianello, M. Giacomin, and A. C. Coroado, “Investigation of the density shoulder formation by using self-consistent simulations of plasma turbulence and neutral kinetic dynamics,” Nuclear Fusion (2021), 10.1088/1741-4326/AC30C9.
  * Russell, Myra, and Stotler (2019) D. A. Russell, J. R. Myra, and D. P. Stotler, “A reduced model of neutral-plasma interactions in the edge and scrape-off-layer: Verification comparisons with kinetic Monte Carlo simulations,” Physics of Plasmas (2019), 10.1063/1.5081670.
  * Russell  _et al._ (2021) D. A. Russell, J. R. Myra, F. Militello, and D. Moulton, “Reduced-model scrape-off layer turbulence (nSOLT) simulations comparing three fueling scenarios,” Physics of Plasmas 28, 092305 (2021).
  * Pan, Told, and Jenko (2016) Q. Pan, D. Told, and F. Jenko, “Fully nonlinear $\delta$f gyrokinetics for scrape-off layer parallel transport,” Physics of Plasmas 23 (2016), 10.1063/1.4964666.
  * Chang  _et al._ (2017) C. S. Chang, S. Ku, A. Loarte, V. Parail, F. Köchl, M. Romanelli, R. Maingi, J. W. Ahn, T. Gray, J. Hughes, B. LaBombard, T. Leonard, M. Makowski, and J. Terry, “Gyrokinetic projection of the divertor heat-flux width from present tokamaks to ITER,” Nuclear Fusion 57 (2017), 10.1088/1741-4326/aa7efb.
  * Shi  _et al._ (2019) E. L. Shi, G. W. Hammett, T. Stoltzfus-Dueck, and A. Hakim, “Full-f gyrokinetic simulation of turbulence in a helical open-field-line plasma,” Physics of Plasmas 26, 012307 (2019).
  * Boesl  _et al._ (2019) M. Boesl, A. Bergmann, A. Bottino, D. Coster, E. Lanti, N. Ohana, and F. Jenko, “Gyrokinetic full-f particle-in-cell simulations on open field lines with PICLS,” Physics of Plasmas 26, 122302 (2019).
  * Dorf and Dorr (2020) M. Dorf and M. Dorr, “Progress with the 5D full-F continuum gyrokinetic code COGENT,” Contributions to Plasma Physics 60, e201900113 (2020).
  * Michels  _et al._ (2021) D. Michels, A. Stegmeir, P. Ulbl, D. Jarema, and F. Jenko, “GENE-X: A full-f gyrokinetic turbulence code based on the flux-coordinate independent approach,” Computer Physics Communications 264, 107986 (2021).
  * Gkeyll (2022) Gkeyll, “The Gkeyll 2.0 Code: Documentation Home,” http://gkeyll.readthedocs.io (2022).
  * Shi  _et al._ (2017) E. Shi, G. Hammett, T. Stoltzfus-Dueck, and A. Hakim, “Gyrokinetic continuum simulation of turbulence in a straight open-field-line plasma,” J. Plasma Phys. 83, 905830304 (2017).
  * Bernard  _et al._ (2019) T. N. Bernard, E. L. Shi, K. W. Gentle, A. Hakim, G. W. Hammett, T. Stoltzfus-Dueck, and E. I. Taylor, “Gyrokinetic continuum simulations of plasma turbulence in the texas helimak,” Physics of Plasmas 26, 042301 (2019), https://doi.org/10.1063/1.5085457 .
  * Hakim  _et al._ (2020a) A. H. Hakim, N. R. Mandell, T. N. Bernard, M. Francisquez, G. W. Hammett, and E. L. Shi, “Continuum electromagnetic gyrokinetic simulations of turbulence in the tokamak scrape-off layer and laboratory devices,” Physics of Plasmas 27 (2020a), 10.1063/1.5141157.
  * Bernard  _et al._ (2020) T. N. Bernard, T. Stoltzfus-Dueck, K. W. Gentle, A. Hakim, G. W. Hammett, and E. L. Shi, “Investigating shear flow through continuum gyrokinetic simulations of limiter biasing in the Texas Helimak,” Physics of Plasmas 27, 062304 (2020).
  * Juno  _et al._ (2018) J. Juno, A. Hakim, J. TenBarge, E. Shi, and W. Dorland, “Discontinuous Galerkin algorithms for fully kinetic plasmas,” Journal of Computational Physics 353, 110–147 (2018).
  * Hakim and Juno (2020) A. Hakim and J. Juno, “Alias-free, matrix-free, and quadrature-free discontinuous galerkin algorithms for (plasma) kinetic equations,” in  _Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis_ (IEEE Press, 2020).
  * Mandell (2021) N. Mandell,  _Magnetic Fluctuations in Gyrokinetic Simulations of Tokamak Scrape-Off Layer Turbulence_ , Ph.D. thesis, Princeton University (2021).
  * Francisquez  _et al._ (2020) M. Francisquez, T. Bernard, N. Mandell, G. Hammett, and A. Hakim, “Conservative discontinuous Galerkin scheme of a gyro-averaged Dougherty collision operator,” Nuclear Fusion 60, 096021 (2020).
  * Hakim  _et al._ (2020b) A. Hakim, M. Francisquez, J. Juno, and G. W. Hammett, “Conservative discontinuous galerkin schemes for nonlinear dougherty–fokker–planck collision operators,” Journal of Plasma Physics 86, 905860403 (2020b).
  * Voronov (1997) G. S. Voronov, “A practical fit formula for ionization rate coefficients of atoms and ions by electron impact: Z = 1-28,” Atomic Data and Nuclear Data Tables 65, 1–35 (1997).
  * Cagas  _et al._ (2017) P. Cagas, A. Hakim, J. Juno, and B. Srinivasan, “Continuum kinetic and multi-fluid simulations of classical sheaths,” Phys. Plasmas 24, 22118 (2017).
  * Cagas (2018) P. Cagas,  _Continuum Kinetic Simulations of Plasma Sheaths and Instabilities_ , Ph.D. thesis, Virginia Polytechnic Institute and State University (2018).
  * Meier (2011) E. T. Meier,  _Modeling Plasmas with Strong Anisotropy, Neutral Fluid Effects, and Open Boundaries_ , Ph.D. thesis, University of Washington (2011).
  * Pauls, Zank, and Williams (1995) H. L. Pauls, G. P. Zank, and L. L. Williams, “Interaction of the solar wind with the local interstellar medium,” Journal of Geophysical Research: Space Physics 100, 21595–21604 (1995).
  * Barnett (1990) C. F. Barnett,  _Atomic Data for Fusion, Vol. 1_ (ORNL, 1990).
  * Isenberg (1987) P. A. Isenberg, “Evolution of interstellar pickup ions in the solar wind,” Journal of Geophysical Research: Space Physics 92, 1067–1073 (1987).
  * Janev and Smith (1993) R. K. Janev and J. J. Smith, “Cross sections for collision processes of hydrogen atoms with electrons, protons and multiply charged ions,” in  _Atomic and Plasma-Material Interaction Data for Fusion_, Vol. 4 (International Atomic Energy Agency, Vienna, 1993).
  * Krstic and Schultz (1998) P. S. Krstic and D. R. Schultz, “Elastic and Related Transport Cross Sections for Collisions among Isotopomers of H+ + H, H+ + H2, H+ + He, H + H, and H + H2,” in  _Atomic and Plasma-Material Interaction Data for Fusion_, Vol. 8 (International Atomic Energy Agency, Vienna, 1998).
  * E. L. Shi (2017) E. L. Shi,  _Gyrokinetic Continuum Simulation of Turbulence in Open-Field-Line Plasmas_ , Ph.D. thesis, Princeton University (2017).
  * Shi, Hakim, and Hammett (2015) E. L. Shi, A. H. Hakim, and G. W. Hammett, “A gyrokinetic one-dimensional scrape-off layer model of an edge-localized mode heat pulse,” Phys. Plasmas 22, 22504 (2015).
  * Stotler  _et al._ (2015) D. Stotler, F. Scotti, R. Bell, A. Diallo, B. LeBlanc, M. Podestà, A. Roquemore, and P. Ross, “Midplane neutral density profiles in the national spherical torus experiment,” Physics of Plasmas 22, 082506 (2015).
  * Scotti  _et al._ (2021) F. Scotti, D. Stotler, R. Bell, B. LeBlanc, S. Sabbagh, V. Soukhanovskii, M. Umansky, and S. Zweben, “Outer midplane neutral density measurements and h-mode fueling studies in nstx-u,” Nuclear Fusion 61, 036002 (2021).
  * Krasheninnikov, D’ippolito, and Myra (2008) S. Krasheninnikov, D. D’ippolito, and J. Myra, “Recent theoretical progress in understanding coherent structures in edge and sol turbulence,” Journal of Plasma Physics 74, 679–717 (2008).
  * Boedo  _et al._ (2003) J. A. Boedo, D. L. Rudakov, R. J. Colchin, R. A. Moyer, S. Krasheninnikov, D. G. Whyte, G. R. McKee, G. Porter, M. J. Schaffer, P. C. Stangeby, W. P. West, S. L. Allen, and A. W. Leonard, “Intermittent convection in the boundary of DIII-D,” Journal of Nuclear Materials 313-316, 813–819 (2003).
  * Rudakov  _et al._ (2005) D. Rudakov, J. Boedo, R. Moyer, P. Stangeby, J. Watkins, D. Whyte, L. Zeng, N. Brooks, R. Doerner, T. Evans, M. Fenstermacher, M. Groth, E. Hollmann, S. Krasheninnikov, C. Lasnier, A. Leonard, M. Mahdavi, G. McKee, A. McLean, A. Pigarov, W. Wampler, G. Wang, W. West, and C. Wong, “Far SOL transport and main wall plasma interaction in DIII-D,” Nuclear Fusion 45, 1589 (2005).
  * Carralero  _et al._ (2017) D. Carralero, M. Siccinio, M. Komm, S. A. Artene, F. A. D’Isa, J. Adamek, L. Aho-Mantila, G. Birkenmeier, M. Brix, G. Fuchert, M. Groth, T. Lunt, P. Manz, J. Madsen, S. Marsen, H. W. Müller, U. Stroth, H. J. Sun, N. Vianello, M. Wischmeier, and E. Wolfrum, “Recent progress towards a quantitative description of filamentary SOL transport,” Nuclear Fusion 57, 056044 (2017).
  * Vianello  _et al._ (2017) N. Vianello, C. Tsui, C. Theiler, S. Allan, J. Boedo, B. Labit, H. Reimerdes, K. Verhaegh, W. Vijvers, N. W. S., Costea, J. Kovacic, C. Ionita, V. Naulin, A. Nielsen, J. Rasmussen, B. Schneider, R. Schrittwieser, M. Spolaore, D. Carralero, J. Madsen, B. Lipschultz, F. Militello, the TCV Team, and the EUROfusion MST1 Team, “Modification of SOL profiles and fluctuations with line-average density and divertor flux expansion in TCV,” Nuclear Fusion 57, 116014 (2017).
  * Vianello  _et al._ (2019) N. Vianello, D. Carralero, C. Tsui, V. Naulin, M. Agostini, I. Cziegler, B. Labit, C. Theiler, E. Wolfrum, D. Aguiam, S. Allan, M. Bernert, J. Boedo, S. Costea, H. D. Oliveira, O. Fevrier, J. Galdon-Quiroga, G. Grenfell, A. Hakola, C. Ionita, H. Isliker, A. Karpushov, J. Kovacic, B. Lipschultz, R. Maurizio, K. McClements, F. Militello, A. Nielsen, J. Olsen, J. Rasmussen, T. Ravensbergen, H. Reimerdes, B. Schneider, R. Schrittwieser, E. Seliunin, M. Spolaore, K. Verhaegh, J. Vicente, N. Walkden, W. Zhang, the ASDEX Upgrade Team, the TCV Team, and the EUROfusion MST1 Team, “Scrape-off layer transport and filament characteristics in high-density tokamak regimes,” Nuclear Fusion 60, 016001 (2019).
  * Myra, Russell, and D’Ippolito (2006) J. R. Myra, D. A. Russell, and D. A. D’Ippolito, “Collisionality and magnetic geometry effects on tokamak edge turbulent transport. I. A two-region model with application to blobs,” Physics of Plasmas 13, 112502 (2006).