# Scaling laws for electron kinetic effects in tokamak scrape-off layer
plasmas

D<EMAIL_ADDRESS>S Mijin2 M Wigram3 F Militello2 R Kingham1

(1Blackett Lab., Plasma Physics Group, Imperial College London, London, United
Kingdom
2UKAEA, Culham Science Centre, Oxon, United Kingdom
3MIT Plasma Science and Fusion Center, Cambridge, MA 02139, USA
)

###### Abstract

Tokamak edge (scrape-off layer) plasmas can exhibit non-local transport in the
direction parallel to the magnetic field due to steep temperature gradients.
This effect along with its consequences has been explored at equilibrium for a
range of conditions, from sheath-limited to detached, using the 1D kinetic
electron code SOL-KiT, where the electrons are treated kinetically and
compared to a self-consistent fluid model. Line-averaged suppression of the
kinetic heat flux (compared to Spitzer-Härm) of up to 50% is observed,
contrasting with up to 98% enhancement of the sheath heat transmission
coefficient, $\gamma_{e}$. Simple scaling laws in terms of basic SOL
parameters for both effects are presented. By implementing these scalings as
corrections to the fluid model, we find good agreement with the kinetic model
for target electron temperatures. It is found that the strongest kinetic
effects in $\gamma_{e}$ are observed at low-intermediate collisionalities, and
tend to increase at increasing upstream densities and temperatures. On the
other hand, the heat flux suppression is found to increase monotonically as
upstream collisionality decreases. The conditions simulated encompass
collisionalities relevant to current and future tokamaks.

##  1 Introduction

The region of unconfined plasma at the edge of tokamaks, called the scrape-off
layer (SOL), is the barrier between the hot core plasma and the solid surfaces
which make up the inside of the reactor. It is necessary to understand plasma
transport in this region, which occurs primarily parallel to the magnetic
field lines, so that accurate predictions can be made for future devices and
steps can be taken to mitigate heat fluxes which may exceed material
constraints.

Transport in SOL plasmas is often treated with fluid models, where a
Braginskii-like set of transport equations [1] may be solved. However, the
presence of steep temperature gradients parallel to the magnetic field, as
would be expected in reactor-class devices, means heat transport (particularly
for the electrons) is dominated by particles with collisional mean free paths
which are long relative to the temperature gradient length scale and so
becomes ‘non-local’. This can be quantified with the upstream collisionality
parameter $\nu_{u}^{*}=L/\lambda_{u}$ [2], defined as the ratio of the
parallel SOL length $L$ and the upstream mean free path $\lambda_{u}$.
Conditions where $\nu_{u}^{*}$ is small and temperature gradients are large
may not be described accurately by a fluid model.

This effect has been explored in recent years [3, 4, 5, 6, 7, 8, 9, 10, 11],
where it is now well-documented that kinetic suppression of the heat flux can
result in steeper temperature gradients and lower target temperatures when
compared to a fluid model. Somewhat less understood is the region of operating
parameter space where such effects may become important, and the consequences
for the overall energy balance at equilibrium (i.e. how energy going into the
SOL makes its way out). It is still unclear whether kinetic effects in
parallel transport pose a significant uncertainty in modelling approaches for
future devices.

Here we present kinetic and fluid simulations of a 1D SOL plasma model, across
a wide range of the relevant parameter space (input power, plasma density and
connection length), in order to assess and understand kinetic deviations from
fluid model predictions. The model will briefly be presented in Section 2,
followed by an explanation of the simulations that have been carried out in
Section 3. We will then summarise the results (Section 4), highlighting the
areas in which kinetic effects are (and are not) seen. Following a discussion
of the observed results (Section 5), we present scaling relationships in terms
of basic SOL parameters for the main kinetic effects seen - enhancement to the
sheath heat transmission coefficient and suppression of the parallel
conductive electron heat flux - in Section 6. These will be used to reproduce
electron temperature profiles from kinetic simulations in a (corrected) fluid
model. A prediction is made for the strength of these kinetic effects in ITER.

##  2 Kinetic and fluid modelling with SOL-KiT

SOL-KiT is a fully implicit 1D plasma code which has been used to study
kinetic effects in parallel electron transport in the SOL [12, 3, 13]. Here a
very brief outline of SOL-KiT is presented, and the reader is referred to [14]
for more details of the model.

In kinetic mode, SOL-KiT solves the VFP equation for electrons in a hydrogenic
plasma along the direction parallel to the magnetic field (the $x$-axis),

$\frac{\partial f(x,\vec{v},t)}{\partial t}+v_{x}\frac{\partial
f(x,\vec{v},t)}{\partial x}-\frac{e}{m_{e}}E\frac{\partial
f(x,\vec{v},t)}{\partial v_{x}}=\sum_{\alpha}C_{e-\alpha},$ (1)

where $f(x,\vec{v},t)$ is the electron velocity distribution, which is a
function of space, velocity and time. $E$ is the electric field along $x$,
$m_{e}$ is the electron mass, $-e$ is the electron charge and $v_{x}$ is the
electron velocity along $x$. The right hand side consists of Fokker-Planck
collisions (electron-electron and electron-ion) and Boltzmann collisions
(electron-neutral). A spherical harmonic decomposition in velocity space is
used to solve this equation as outlined in [15]. Azimuthal symmetry is assumed
about the $x$-axis so that the magnetic field may be ignored.

The $x$-axis spans from the midplane at $x=0$ (‘upstream’), to the plasma
sheath boundary at $x=L$, where $L$ is the domain length, which is half the
connection length (the parallel distance between two strike points in a
divertor SOL). Power enters the plasma upstream, and leaves at the sheath or
through plasma-neutral collisions. The included collisional processes are
electron-impact ionisation and excitation; the inverse of these processes
(three-body recombination and collisional de-excitation); and resonant charge
exchange (CX) between the ions and neutrals. In addition, radiative
recombination and de-excitation are modelled, which allows energy to leave the
plasma-neutral system.

In fluid mode, moments of equation (1) are solved instead, allowing for a
direct comparison between a fluid and kinetic treatment. Evolution equations
for the electron temperature $T_{e}$, flow velocity $u_{e}$ and density
$n_{e}$ are solved. The $T_{e}$ equation is closed with Braginskii/Spitzer-
Härm heat flow [1, 16],
$q_{\parallel,e}=-\kappa_{e}\frac{\partial\left(kT_{e}\right)}{\partial
x}+0.71n_{e}kT_{e}\left(u_{e}-u_{i}\right)$.

Quasi-neutrality is enforced by setting the ion density $n_{i}=n_{e}$. The
parallel electric field $E$ is evolved with Ampère-Maxwell’s law,
$\frac{\partial E}{\partial t}=-(j_{e}-j_{i})/\epsilon_{0}$, where the ion and
electron currents are $j_{i,e}=\pm en_{i,e}u_{i,e}$. The implicit time
integration used by SOL-KiT means this results in ambipolarity when using
timesteps large relative to the plasma oscillation period.

In kinetic mode, where fluid electron quantities are required, for example
$u_{e}$ in calculating $E$, the appropriate velocity moments of $f$ are taken.

In order to provide a realistic background on which to solve electron
transport, SOL-KiT also models the hydrogenic ions and neutral atoms, both of
which are treated with fluid models (discussed more in the next section).
Atomic processes including ionisation, recombination, excitation and de-
excitation are handled by solving a collisional-radiative model (CRM) for the
neutral atoms alongside the transport equations, using fundamental cross-
sections and rates from Janev [17] and NIST [18]. The fully time-dependent CRM
coupled with neutral transport gives a non-coronal model. Particle sources and
the effect of electron-neutral collisions on the electrons are evaluated with
the inelastic Boltzmann collision operator. In fluid mode, the kinetic
collision operators are evaluated with a Maxwellian electron distribution. At
present no impurity species are treated by SOL-KiT, although this is planned
for a future update to the code.

The upstream boundary is reflective. At the sheath boundary, the Bohm
criterion is applied so that plasma flow reaches the sound speed $c_{s}$. All
plasma particle flux across this boundary is lost and recycled as atomic
neutrals, which are placed in the last spatial cell. In kinetic mode, the
logical boundary condition [19] is applied to the electron distribution, where
the forward-going part of the distribution is reflected and truncated at some
velocity $v_{c}$, which is calculated iteratively to ensure equal electron and
ion fluxes. In fluid mode, this same boundary condition manifests via the heat
flux at the sheath entrance, $q_{sh,e}=\gamma_{e}kT_{e,t}\Gamma_{t}$, for
target temperature $T_{e,t}$ and particle flux $\Gamma_{t}$, and where

$\gamma_{e}=2-0.5\ln\left(2\pi\left(1+T_{i,t}/T_{e,t}\right)m_{e}/m_{i}\right)$
(2)

is the sheath heat transmission coefficient [2].

###  2.1 Extensions to SOL-KiT

For this study, this SOL-KiT model described in [14] has been extended to
provide a more realistic background plasma on which to solve the electron
transport, as well as to provide some reduction in compute time. As some of
these improvements have not yet been documented elsewhere, they are outlined
here.

####  2.1.1 Ion temperature equation

Firstly, an ion temperature equation has been added to the model, allowing for
ion energy transport independently from the electrons as well as additional
channels for ion-electron and ion-neutral energy transfer. More detail on this
change is provided in [13], omitted here for brevity. The ion heat flow is
$q_{\parallel,i}=-\kappa_{i}\frac{\partial\left(kT_{i}\right)}{\partial x}$,
using the Spitzer-Härm $\kappa_{i}$ [16]. The sheath boundary condition on the
ion temperature equation is equivalent to the one for the fluid electron
model, $q_{sh,i}=\gamma_{i}kT_{i,t}\Gamma_{t}$, with $\gamma_{i}=2.5$ [2]. In
fluid mode, ion-electron energy transfer is treated with
$Q_{ie}=-Q_{ei}=-\frac{3m_{e}}{m_{i}}\frac{n_{e}k}{\tau_{e}}\left(T_{i}-T_{e}\right)$
[1], where $\tau_{e}$ is the electron collision time. In kinetic mode, we
instead take the energy moment of the Fokker-Planck collision operator
assuming Maxwellian ions, $Q_{ei}=\int\frac{1}{2}m_{e}v^{2}C_{ei}d\vec{v}$.
The two approaches are equivalent for Maxwellian electrons [13].

Figure 1: Geometry of the quasi-2D neutral transport model; $x$ and $y$ are
the parallel and perpendicular directions respectively; $\phi$ and $\theta$
are the toroidal and poloidal directions; $\alpha$ is the pitch angle of the
magnetic field. In an $x-y$ coordinate system, assuming toroidal symmetry, $y$
coordinates are transformed into $x$ coordinates by translation in the $\phi$
direction (shown by the $x^{\prime}-y^{\prime}$ axes).

####  2.1.2 Fluid neutral model

Secondly, velocity and temperature equations have been added for the
(hydrogenic) neutral particles, in order to capture the important processes of
energy and momentum transfer between ions and neutrals via particle exchange
from ionisation, recombination and CX. In order to capture neutral transport
in the direction perpendicular to the magnetic field lines, we use the
assumption of toroidal symmetry to project perpendicular velocities back on to
the parallel direction[20, 21]. The geometry is shown in Figure 1. For a
velocity $\vec{u}=u_{\parallel}\hat{x}+u_{\perp}\hat{y}$, where $\hat{y}$ is
the unit vector perpendicular to the magnetic field, we can obtain an
effective velocity in the parallel direction,
$u_{eff}=u_{\parallel}+u_{\perp}/\tan\alpha$, where
$\alpha=\tan^{-1}(B_{\theta}/B_{\phi})$ is the pitch angle of the magnetic
field lines relative to the wall. Similarly, we can relate perpendicular
gradients to parallel gradients via $\frac{\partial}{\partial
y}=\frac{1}{\tan\alpha}\frac{\partial}{\partial x}$. With this, we obtain a 1D
neutral transport model which captures something of the cross-field neutral
transport. The parallel neutral velocity equation is

$\begin{split}\frac{\partial u_{n\|}}{\partial t}=-u_{n,eff}\frac{\partial
u_{n\|}}{\partial x}-\frac{1}{m_{n}n_{n}}\frac{\partial p_{n}}{\partial
x}-\frac{S_{n}}{n_{n}}u_{n\|}+\frac{1}{m_{n}n_{n}}R_{n\|},\end{split}$ (3)

and the perpendicular velocity equation is

$\begin{split}\frac{\partial u_{n\perp}}{\partial t}=-u_{n,eff}\frac{\partial
u_{n\perp}}{\partial x}-\frac{1}{m_{n}n_{n}\tan\alpha}\frac{\partial
p_{n}}{\partial
x}-\frac{S_{n}}{n_{n}}u_{n\perp}+\frac{1}{m_{n}n_{n}}R_{n\perp},\end{split}$
(4)

where $p_{n}=n_{n}kT_{n}$ is the neutral pressure, $S_{n}$ is the neutral
particle source and $R_{n,\parallel}$, $R_{n,\perp}$ are the parallel and
perpendicular momentum sources. These contain contributions from particle
sources and CX. We use a similar form of the CX friction term outlined by
Meier[22], along with velocity-dependent CX cross-sections from Janev[17]. For
calculating $R_{n,\perp}$, the plasma is assumed to be stationary in the
perpendicular direction. The pitch angle $\alpha$ in practice determines the
degree to which neutrals are able to transport upstream in this 1D model. For
normal incidence, $\alpha=90^{\circ}$, $u_{n,eff}=u_{n\parallel}$ and the
model is 1D along $x$. The connection of this model with the pressure-
diffusion neutral model in[23], and a similar quasi-2D transport model
implemented in SD1D[21], can be seen by dropping all but the pressure gradient
and friction terms in (4). The perpendicular velocity is then
$u_{n,\perp}\propto\frac{1}{\tan\alpha}\frac{\partial p_{n}}{\partial x}$,
which can be inserted into $u_{n,eff}$ to yield a convective-diffusive model,
where the diffusion coefficient is enhanced by a geometric factor proportional
to $1/\tan^{2}\alpha=(B_{\phi}/B_{\theta})^{2}$.

SOL-KiT solves a continuity equation for each excited state of a neutral
hydrogenic species, denoted with subscript $b$, each with its own density
$n_{b}$ and particle source $S_{b}$. Total neutral density is therefore
$n_{n}=\sum_{b}n_{b}$, and $S_{n}=\sum_{b}S_{b}$. The additional velocity
equations (3 & 4) yield a modified continuity equation for each species,

$\frac{\partial n_{b}}{\partial
t}=-\frac{\partial\left(n_{b}u_{n,eff}\right)}{\partial x}+S_{b}.$ (5)

The neutral temperature $T_{n}$ is also evolved,

$\displaystyle\frac{\partial kT_{n}}{\partial t}=-u_{n,eff}\frac{\partial
T_{n}}{\partial x}+\frac{2}{3}\left[\frac{Q_{n}}{n_{n}}\right.-$
$\displaystyle\left.kT_{n}\frac{\partial u_{n,eff}}{\partial
x}-\frac{S_{n}}{n_{n}}\left(\frac{3}{2}kT_{n}-\frac{1}{2}m_{n}u_{n\|}^{2}\right)\right.$
(6)
$\displaystyle-\left.\frac{1}{n_{n}}\left(1+\frac{1}{\tan^{2}\alpha}\right)\frac{\partial
q_{n}}{\partial x}-\frac{u_{n\|}R_{n\|}}{n_{n}}\right],$

where $q_{n}$ is the neutral heat flow,
$q_{n}=-2.4\left(\frac{n_{n}T_{n}}{m_{n}\nu_{CX}}\right)\frac{\partial
T_{n}}{\partial x}$ [24] and $\nu_{CX}$ is the CX collision frequency. The
energy source term $Q_{n}$ contains contributions from CX collisions and the
energy transfer from particle exchange during ionisation and recombination. As
in the case of CX friction, for the CX contribution to $Q_{n}$ we use a
similar form to that in[22].

The boundary condition on the fluid neutral model is a simplified version of
the approach outlined by Horsten et al.[25, 23], where the neutral
distribution at the boundary is taken to be made up of a set of one-sided
distributions from the incident, reflected and recycled neutrals. We assume
100% reflection of the neutrals incident on the wall and no ‘fast’ recycling
of ions, i.e. they are thermally released from the wall as molecules which
quickly dissociate. The ion particle flux crossing the sheath, $\Gamma_{i,t}$,
is recycled and re-enters as a neutral flux at the same location. These
recycled particles are assumed to be emitted from the wall mono-energeticaly
at the Frank-Condon dissocation energy $T_{FC}=3$eV. Neutrals incident on the
wall are reflected back with some loss of energy. The net energy flux of
neutrals at the boundary is therefore (for 100% plasma recycling)
$q_{n,t}=\gamma_{n}kT_{n}n_{n}c_{n,s}-kT_{FC}\Gamma_{i,t}$, where
$c_{n,s}=\sqrt{kT_{n}/m_{i}}$ is the neutral sound speed at the boundary and
$\gamma_{n}$ is a neutral wall heat transmission coefficient. Here we take
$\gamma_{n}=0.25$ as a reasonable estimate derived from energy reflection
coefficients for deuterium incident on a tungsten wall from the TRIM
database222http://www.eirene.de/html/surface_data.html.

####  2.1.3 Atomic state bundling

The final change is to introduce bundling on the neutral states, so that a
reduced number of states are evolved without significantly altering the
plasma-neutral physics. This is desirable as it reduces the number of
electron-neutral collision operators in (1) which need to be computed, which
increases with the square of the number of evolved neutral states and
therefore represents a code bottleneck.

We will define a cut-off neutral state index $\tilde{b}$, and evolve all
states from $b=1$ to $\tilde{b}-1$ as normal, and group all higher states
(from $\tilde{b}$ to $b_{max}=30$) into a single bundle $\beta$, with bundle
density $n_{\beta}=\sum_{b\in\beta}n_{b}$. For this bundle, we assume a
Boltzmann state distribution at the electron temperature, giving the density
ratio
$n_{b}/n_{\tilde{b}}=(b/\tilde{b})^{2}e^{-\left(\varepsilon_{\tilde{b}}-\varepsilon_{b}\right)/kT_{e}}$
for $b\geq\tilde{b}$, where $\varepsilon_{b}$ is the ionisation energy of
level $b$, which for hydrogen is $\varepsilon_{b}=13.6\mathrm{eV}/b^{2}$.

All particle, momentum and energy source terms relating to this bundle can
then be computed by evaluating the appropriate collision operator with bundle-
averaged cross-sections, $\langle\sigma\rangle_{\beta}$ and ionisation
energies, $\langle\varepsilon\rangle_{\beta}$. These quantities are functions
of the electron temperature, and are precomputed on a relevant range of
$T_{e}$ and interpolated at runtime. For example, the electron energy sink
from ionisation of neutrals in $\beta$ is

$\begin{split}Q_{\beta}^{ion}&=\langle\varepsilon\rangle_{\beta}S_{\beta}^{ion}\\\
&=\langle\varepsilon\rangle_{\beta}n_{\beta}n_{e}\langle
K\rangle_{\beta}^{ion}\\\
&=\langle\varepsilon\rangle_{\beta}n_{\beta}n_{e}4\pi\int_{\infty}dvv^{3}\frac{f_{0}(v)}{n_{e}}\langle\sigma\rangle_{\beta}^{ion},\end{split}$
(7)

where $S_{\beta}^{ion}$ is the ionisation particle source from $\beta$,
$\langle K\rangle_{\beta}^{ion}$ is the bundle-averaged ionisation rate
coefficient, $n_{e}$ is the electron density and $f_{0}$ is the isotropic part
of the electron distribution. We also include an effective collision operator
to capture excitation collisions within a bundle by assuming the net
excitation rate is equal to the rate of radiative de-excitation.

Using $\tilde{b}=5$, this bundling method is found to predict neutral
densities and radiative losses to within a few percent of that predicted by
directly evolving all states, while providing around an order of magnitude
reduction in compute time. Further details of this bundling approach will be
published separately.

##  3 Parameter scan simulations

Figure 2: SOL-KiT simulations carried out for this study, where each black dot
represents a simulation at a given $q_{in}$ and $\langle n_{0}\rangle$. The
exact values of $n_{u}$ and $T_{e,u}$ are taken from the kinetic simulations.
For context, lines of constant collisionality are shown at $\nu^{*}_{e,u}=$
25, 50 and 100.

Two SOL parameters which we have some degree of control over in tokamaks (and
which determine SOL behaviour) are $T_{u}$ and $n_{u}$, the plasma temperature
and density measured at some upstream location. In these simulations, we vary
these by changing the input power flux to the SOL from the core, $q_{in}$, and
the initial line-averaged plasma density, $\langle
n_{0}\rangle=\frac{1}{L}\int_{0}^{L}n_{e}(t=0)dx$, where the plasma is fully
ionised at initialisation. The total number of particles (plasma plus
neutrals) stays constant in the simulations due to 100% recycling of the
target particle flux. The input heat flux is distributed uniformly across the
heating region, $L_{heat}$, as a volumetric energy source, $Q_{in}$, such that
$q_{in}=\int_{0}^{L_{heat}}Q_{in}dx=L_{heat}Q_{in}$.

Of interest in this study is how conditions upstream determine the electron
transport, and a useful measure of this is the electron upstream
collisionality parameter $\nu_{e,u}^{*}$, defined as the ratio of the
connection length $L$ to the electron Coulomb mean free path upstream
$\lambda_{ee,u}$ [2],

$\nu_{e,u}^{*}=\frac{L}{\lambda_{ee,u}}\simeq
10^{-16}\frac{Ln_{u}}{T_{e,u}^{2}}$ (8)

for $T_{e,u}$ in [eV] and $n_{u}$ in [m-3]. Note that this differs slightly
from some forms of $\nu^{*}=L/\lambda_{u}$ employed in the literature (e.g.
[4]), and $\nu_{e,u}^{*}$ here will typically be smaller than collisionality
defined in terms of total plasma temperature because $T_{i,u}>T_{e,u}$
generally.

For a deuterium plasma, a number of density scans were performed at different
input powers and connection lengths. Connection lengths ranged from $L=11.93$m
to 30.97m, input powers from $q_{in}=4$MWm-2 to 128MWm-2, and densities from
$\langle n_{0}\rangle=1.0\times 10^{19}$m-3 to $1.4\times 10^{20}$m-3. With
these input parameters, the simulations cover $\nu_{e,u}^{*}$ from 6.3 to
203.6. At the lowest collisionalities the plasma is sheath-limited, while
detachment is reached at the highest values of $\nu^{*}_{e,u}$ (measured by
rollover of the target particle flux).

$q_{in}$ is distributed over approximately the first third of the domain and
spread equally between the ions and electrons; 100% of plasma particles lost
to the sheath are recycled as neutrals, and the pitch angle used in the
neutral model was $\alpha=15^{\circ}$. 100 spatial cells were used, which were
spaced logarithmically with higher resolution close to the target. For the
simulations with the longest connection length, the spatial grid widths ranged
from 2.28m upstream to 1.05mm at the target. In velocity space (for kinetic
electron runs), a geometric grid of 80 cells was used up to a velocity of
$\simeq$ 12$v_{th,0}$,where $v_{th,0}$ is the thermal velocity of electrons at
a reference temperature of 10eV. The resolution was higher at low velocities,
such that grid widths ranged from 0.05$v_{th,0}$ to 0.35$v_{th,0}$. In the
kinetic runs, the kinetic equation (1) was solved up to the spherical harmonic
$l_{max}=3$.

To reach equilibrium, determined by when the power and particle balance has
converged, the kinetic simulations with SOL-KiT each take a few weeks running
on 8 CPUs, while the fluid simulations typically take a day or less on 4
cores.

These simulations are situated on the $T_{e,u}$-$Ln_{u}$ plane along with
lines of constant $\nu_{e,u}^{*}$ in Figure 2. For reference, present-day
tokamaks (JET, DIII-D, etc.) operate with $Ln_{u}\simeq 10^{20}-10^{21}$m-2
and $T_{e,u}\simeq 20-60$eV. Future devices like ITER and DEMO will operate
with $Ln_{u}\sim 1\times 10^{21}-4\times 10^{21}$m-2 and $T_{e,u}\sim
150-300$eV[26, 27]. Simulating such regimes kinetically is computationally
demanding, but this study represents an attempt to explore kinetic effects in
regimes beyond those attainable in existing devices. The simulations presented
here fall short of reaching the highest values of $T_{e,u}$ expected in future
tokamaks, but do encompass reactor-relevant values of $Ln_{u}$, detached
conditions, and a broad range of upstream collisionalities.

##  4 Results

(a) $q_{in}=4$MWm-2, $L=11.93$m.

(b) $q_{in}=16$MWm-2, $L=11.93$m.

(c) $q_{in}=80$MWm-2, $L=18.79$m.

(d) $q_{in}=128$MWm-2, $L=30.96$m.

Figure 3: Target electron temperatures, $T_{e,t}$, and ion fluxes,
$\Gamma_{i,t}$, for four of the density scans at different input powers and
connection lengths. Results with fluid and kinetic electrons are shown. The
lowest input power runs exhibit detachment, indicated by rollover of
$\Gamma_{i,t}$.

We start by displaying in Figure 3 the target temperatures and particle fluxes
of four of the density scans carried out. Rollover of the target flux, an
indicator of detachment onset, is expected when particle, momentum and power
losses are sufficient to reduce the target flux despite increasing plasma
density. Only the lowest-power run reaches flux rollover here, while all other
runs are partially or fully attached. The absence of rollover at high input
powers is not unexpected given the lack of impurity radiation in the SOL-KiT
model, which would provide additional power dissipation. It can be seen in
Figure 3 that target temperatures are lower when the electrons are treated
kinetically, with the biggest differences in both absolute and relative terms
occurring at low upstream densities. However, there are only small differences
observed in the rollover behaviour in the one density scan which does reach
detachment, 3(a), where a kinetic treatment of the electrons does not change
the position of detachment onset when varying $n_{u}$, i.e. the qualitative
behaviour of $\Gamma_{i,t}$ is the same. There is a decrease in $\Gamma_{i,t}$
with kinetic electrons at high upstream densities, up to around 15%, while
there is a slight increase at low densities.

The reduction in target temperatures in Figure 3 is a reflection of the
suppression of the parallel conductive electron heat flux which is observed in
kinetic simulations, as has been observed in other kinetic studies of parallel
transport [6, 12, 28, 7]. This can be seen in Figure 4(a), which shows
temperature profiles for two simulations at low collisionality, along with
differences in target electron temperatures across all simulations in Figure
5. Figure 4(b) shows the ratio of the kinetic to Spitzer-Härm heat flux
calculated on the kinetic plasma profiles in 4(a). This suppression of the
heat flux relative to that predicted by a fluid treatment, where for a given
heat conductivity $\kappa$ the heat flux is $q_{\parallel,e}=-\kappa\nabla T$,
arises due to fast electrons not depositing their energy locally due to their
large mean free path relative to the temperature gradient length scale. This
means that a steeper temperature gradient is required to achieve the same heat
flux along the SOL, which is fixed by $q_{in}$. There is an uptick in $T_{e}$
close to the wall in the least collisional simulations, which is seen in the
hotter $T_{e}$ profile in Figure 4(a) and as a spike in the heat flux ratio in
Figure 4(b).

In Figure 6 we can see the accumulation of fast electrons from upstream at the
target, where an electron energy distribution from just in front of the sheath
boundary is shown. There is a clear enhanced high-energy tail, while the
thermal bulk is close to the local Maxwellian.

(a) Fluid and kinetic electron temperature profiles.

(b) Ratio of the kinetic to Spitzer-Härm conductive heat flux for the kinetic
simulations in (a).

Figure 4: Kinetic heat flux suppression resulting in steeper temperature
gradients and lower target temperatures for two low collisionality simulations
($q_{in}=64$MWm-2, $L=11.93$m). Figure 5: Difference in target temperatures
for kinetic vs. fluid electrons across all simulations, $\Delta
T_{e,t}=T_{e,t}^{fluid}-T_{e,t}^{kinetic}$. Figure 6: Electron energy
distribution (isotropic part) close to the wall in a SOL-KiT simulation
($q_{in}=64$MWm-2, $L=11.93m$, $\langle n_{0}\rangle=5\times 10^{19}$m-3).
Dashed line is the local Maxwellian. A prominent high-energy tail and
thermalised bulk can both be seen. $T_{e}=7.3$eV, $n_{e}=3.2\times
10^{20}$m-3.

In Figure 7, we show the kinetic enhancement of the sheath heat transmission
coefficient, $\gamma_{e}$, shown as
$\Delta\gamma_{e}=\gamma_{e}^{kinetic}-\gamma_{e}^{fluid}$. Maximum and
minimum values of $\Delta\gamma_{e}$ seen here are 4.96 and 0.47. Differences
in $\gamma_{e}$ here arise because, in kinetic mode, $\gamma_{e}$ is
calculated self-consistently from the logical boundary condition on the
electron distribution, whereas in fluid mode it is calculated from fluid
quantities in the classical way (2). For reference, in fluid mode typically
$\gamma_{e}\simeq 4.8$, and this varies slowly with SOL conditions. In Figure
7, there is a non-monotonic behaviour in $\Delta\gamma_{e}$, where the
classical value is approached at both high and low collisionalities. Similar
behaviour was seen in a power scan in [12] and in a collisionality scan in
[29]. The largest differences occur at low-intermediate collisionalities, but
there is an additional increase in magnitude of this effect along lines of
constant collisionality, moving towards larger $T_{e,u}$ and $n_{u}$. This can
be seen by tracing along the red $\nu^{*}_{e,u}=25$ line in Figure 7, where
simulations at higher $T_{e,u}$ have larger $\Delta\gamma_{e}$. Additionally,
even at the highest collisionalities, where we would expect good agreement
between fluid and kinetic predictions, there is a residual
$\Delta\gamma_{e}\simeq 0.5$. It would therefore appear that convergence of
$\gamma_{e}$ to the fluid value is slow as a function of $\nu_{e,u}^{*}$.

Figure 7: Enhancement to the electron sheath heat transmission coefficient,
$\Delta\gamma_{e}=\gamma_{e}^{kinetic}-\gamma_{e}^{fluid}$, across all
simulations.

Given the enhancement in $\gamma_{e}$ for kinetic electrons, it is natural to
investigate the heat lost to the sheath boundary,
$q_{sh,e}=\gamma_{e}kT_{e,t}\Gamma_{i,t}$, where $T_{e,t}$ is the electron
temperature at the target. This is shown in Figure 8, where the variation of
$q_{sh,e}$ for kinetic and fluid simulations with $n_{u}$ is plotted for
simulations grouped by connection length and input power.. In contrast to the
kinetic enhancement in $\gamma_{e}$, we see that $q_{sh,e}$ is generally in
good agreement for kinetic and fluid simulations. This is perhaps not
surprising, since $q_{sh,e}$ is to a large extent fixed by $q_{in}$, as well
as the fact that kinetic enhancement in $\gamma_{e}$ may be offset by the
reduction in target temperatures (Figure 4(a)). However, this does show that
the overall power balance in these simulations (for example, how much power is
radiated away by electron-neutral collisions) is broadly unchanged despite
modifications to the conductive transport as well as behaviour at the
boundary.

Figure 8: Sheath heat flux from electrons, $q_{sh,e}$, for kinetic and fluid
simulations, grouped into density scans at different connection lengths and
input powers.

Finally, it is worth commenting that the findings in [12] and [5], that
electron-neutral reaction rates are well-approximated by Maxwellian-averaged
rates, is replicated in this study. Differences in rate coefficients for
deuterium ionisation and line radiation are negligible. Some differences in
the total line-integrated particle source do exist in kinetic simulations
(which takes into account ionisation and recombination as well as multi-step
processes involving ex/de-excitation), but these are all under 10%, and are
driven purely by differences in temperature profiles.

##  5 Discussion

The detachment behaviour observed in these simulations contrasts somewhat with
the observation of flux rollover at $n_{u}\sim 2\times 10^{19}$m-3 in the
study by Dudson et al. in[21], which uses a similar simulation setup for the
SD1D fluid code with $q_{in}=50$ MWm-2 and $L=30$m. With the model presented
here, $q_{in}$ must be reduced to below 25 MWm-2 to reach detachment at this
connection length. A separate investigation has highlighted that an
underestimation of the ionisation rate at high electron densities in SD1D may
be the part of the reason for this difference[30], as well as the fact that
SOL-KiT does not include impurity radiation or flux tube expansion.

The uptick in $T_{e}$ close to the target in the least collisional simulations
is seen in the hotter $T_{e}$ profile in Figure 4(a) and in a spike in the
heat flux ratio in Figure 4(b), where $q_{\parallel,e}^{kinetic}$ can become
negative in some simulations. This is the reason why the least collisional
simulations do not have the largest differences in target temperatures (Figure
5) despite exhibiting higher flux suppression, where the uptick in $T_{e}$
close to the wall to some extent cancels the increased temperature gradient.
This feature comes primarily from the perpendicular electron temperature,
which has been seen elsewhere [31, 32], where it is proposed to be related to
an enhancement in the parallel flux of the perpendicular electron energy.

The results in Figure 3 suggest that a kinetic treatment of electron transport
does not predict easier access to detachment. However, it should be said that
only the most collisional density scan considered here actually reaches
detachment, which is where kinetic effects are weakest, and it is possible
that differences in detachment thresholds may exist in conditions with higher
$T_{e,u}$.

The largely unchanged target flux behaviour in simulations with kinetic
electrons, along with broadly similar heat loads to the walls (Figure 8), is
indicative of the fact that a kinetic electron treatment does not
significantly change the particle, momentum or power balance at equilibrium in
this 1D SOL model. This is despite strong heat flux suppression (Figure 4) and
enhancement of the sheath heat transmission coefficient (Figure 7). This can
be understood as resulting from the fact that heat transport is primarily
determined by the input power $q_{in}$. While a modified temperature profile
is needed in kinetic mode to achieve the same parallel heat flux in these
simulations, this is compensated by an enhanced $\gamma_{e}$ which gives a
similar $q_{sh,e}$, leaving the power balance broadly unchanged. In addition,
differences in the temperature profile are insufficient to significantly
change the particle source from electron-neutral interactions.

This power balance behaviour would not necessarily continue to be the case in
the presence of strong radiation sinks from impurities, where modified
temperature profiles and reaction rates could lead to differences in overall
energy transport. This is the subject of an ongoing study.

The unchanged power balance despite the presence of kinetic effects in
parallel heat flux and $\gamma_{e}$ suggests that any attempt to capture
kinetic effects in a fluid framework would need to consider both phenomena. As
such, approaches which treat only the modified heat flux [4] or the boundary
condition [33] may not provide better predictive power than a purely fluid
model.

The strong enhancements to $\gamma_{e}$ are a result of the modified potential
drop across the sheath when calculated kinetically. This depends on $v_{c}$,
the cut-off velocity at which the electron distribution at the sheath is
truncated. The value of $v_{c}$ is set to ensure ambipolar particle flux, and
is therefore somewhat sensitive to a strongly enhanced tail of the electron
distribution, of the kind observed in Figure 6. This may therefore have
consequences for measurements of the electron heat flux, which requires
knowledge of $\gamma_{e}$, and the electron temperature, where electrons are
assumed to be close to Maxwellian. In [34], Tskhakaya et al. observed large
discrepancies in simulated Langmuir probe measurements of $T_{e}$ due to the
departure of the electrons from Maxwellian. Similarly to the study presented
here, the largest discrepancies were seen at intermediate collisionalities.

It is worth noting that the most collisional kinetic simulations here exhibit
a $\gamma_{e}$ which is $\sim 0.5$ higher than the classical value (2). Since
the electron and ion momenta are treated separately and used to solve for the
electric field in this model, this residual discrepancy is not from the pre-
sheath acceleration of the ions. Instead, it is a result of the departure of
the electrons from Maxwellian at the wall.

As discussed, and shown in [12] and [5], electron-impact ionisation rates of
hydrogen are very well approximated by a Maxwellian distribution. This is
unsurprising when considering the distribution shown in Figure 6, which is
non-Maxwellian in the tail but very close to Maxwellian in the thermal bulk.
Given the energy threshold of inelastic processes involving hydrogen are all
at or lower than 13.6eV, the Maxwellian bulk electrons dominate the rates.
This does suggest however that inelastic processes with threshold energies
$\gtrsim$ 50eV (e.g. ionisation of high-Z impurities) may exhibit strong
kinetic enhancement due to the presence of this enhanced tail. An ongoing
study is currently investigating this.

Figure 9: Electron energy distributions at two locations, upstream
($T_{e}=39.9eV$) and close to the wall ($T_{e}=8.0$eV). Electrons are close to
Maxwellian upstream, and the fast tail survives to some extent further
downstream. The gradient of the tail, related directly to the temperature on
these axes, is the same in both cases.

The enhanced tail of the electron distribution seen close to the wall in these
simulations is an imprint from upstream, as can be seen clearly in Figure 9,
where the upstream distribution is plotted alongside that close to the target.
Therefore, if we assume the tail of the distribution at the target has
temperature $T_{e,t}^{tail}=T_{e,u}$, then two conditions existing
simultaneously produce a ‘strongly enhanced’ tail, which can lead to strong
kinetic effects as discussed. These are $T_{e,u}\gg T_{e,t}$ and
$\nu^{*}_{e,u}$. In this study, we see that the imprint can survive up to
moderate values of $\nu^{*}_{e,u}$ and hence drive kinetic effects, for
example in the peak enhancement to $\gamma_{e}$ occurring at
$\nu^{*}_{e,u}\simeq 20$ (Figure 7). It is the interplay of upstream
collisionality and parallel temperature drops which determines the strength of
this imprint. For tokamak edge plasmas with large $T_{e,u}$ as well as
significant power dissipation via impurities, we might expect both of these
conditions to be satisfied.

Contrary to the heat flux suppression, which appears to be a monotonic
function of $\nu^{*}_{e,u}$, the enhancement to $\gamma_{e}$ is more complex.
It peaks at $\nu^{*}_{e,u}\simeq 20$, but also appears to increase for
increasing $T_{e,u}$ at constant $\nu^{*}_{e,u}$. If this behaviour can be
extrapolated to reactor-class devices then we may expect significant
deviations from classical values of $\gamma_{e}$. This is discussed further in
the next section.

##  6 Scaling relationships for observed kinetic effects

Any attempt at capturing kinetic effects at equilibrium in a fluid model of
the scrape-off layer would appear to need to capture both modifications to the
heat flux and enhancement to the sheath heat transmission coefficient. While
models do exist for the former [35, 36, 37], they do not typically provide a
self-consistent method for calculating modifications to the boundary
behaviour. In [38], Tskhakaya et al. provide fits to the modifications to
$\gamma_{e}$ and parallel heat flux for the time-dependent response to a
simulated edge-localised mode (ELM). Here, we present fits to the kinetic
modifications to $\gamma_{e}$ and $q_{\parallel,e}$ seen at equilibrium across
a range of $T_{e,u}$ and $n_{u}$ (and hence $\nu^{*}_{e,u}$), presented as
functions of basic SOL parameters.

The approach taken in developing the relationships presented here has been to
parameterise each simulation in terms of SOL quantities which are either
control parameters or are easily obtainable from experiment or simulation. We
then quantify the two strongest kinetic effects observed in the kinetic
simulations, which we will call kinetic factors. These are the line-averaged
heat flux suppression

$f_{\kappa_{e}}=\frac{1}{L}\int
q_{\parallel,e}^{kinetic}/q_{\parallel,e}^{SH}dx,$

and the enhancement to $\gamma_{e}$,

$\Delta\gamma_{e}=\gamma_{e}^{kinetic}-\gamma_{e}^{fluid}.$

We have then used a least squares fit, allowing a set of fit parameters to
vary, to find combinations of these SOL parameters which provide the best
predictive power for the kinetic factors. In general, it has been found that
there are several forms of the relationships between these kinetic factors and
SOL parameters which provide a comparably good fit. In this case, we have
opted to present fits which are functions of a small number of SOL parameters
to avoid over-fitting the data and keep the approach simple, and which also
have physically realistic asymptotic behaviour (for example, $f_{\kappa_{e}}$
should be close to 1 at high collisionality).

The fits presented here are in terms of control/upstream SOL parameters only.
If there is knowledge of the target conditions, for example in a SOL fluid
code, relationships with very good predictive power have been found, in
particular for $\Delta\gamma_{e}$. One example is provided below. However, it
has been found that some of these do not extrapolate well to conditions where
the upstream collisionality is low and power dissipation is large, for example
due to impurity radiation. Conditions of this nature have not been simulated
yet with SOL-KiT, although this is planned for future work.

Figure 10: The variation of $f_{\kappa_{e}}$ as a function of $\nu_{e,u}^{*}$
in simulations alongside the fit in equation (9).

The proposed relationship for the heat flux suppression is

$f_{\kappa_{e}}=a\exp{(b(\nu_{e,u}^{*})^{c})}+d,$ (9)

with

$a=0.696,\quad b=-8.059,\quad c=-1.074,\quad d=0.260.$

This is a function of $\nu_{e,u}^{*}$ only, which was found to be sufficient
for an accurate prediction of $f_{\kappa_{e}}$. This is plotted against the
simulation data in Figure 10, where the RMS error on the fit is 0.01. There is
clearly significant heat flux suppression at low upstream collisionalities. If
this behaviour can be extrapolated beyond the region of $\nu_{e,u}^{*}$
explored here, this will lead to increased temperature gradients in sheath-
limited regimes where the plasma would otherwise be expected to be nearly
isothermal. This can be seen to some extent in the temperature profiles of the
low collisionality runs in Figure 4(a).

For the enhancement to the electron sheath heat transmission coefficient, we
have

$\Delta\gamma_{e}=\frac{a(q_{in})^{b}\exp(c\nu_{e,u}^{*})}{1+d\exp(e\nu_{e,u}^{*})}+0.5$
(10)

with

$a=9.93\times 10^{-4},\quad b=0.186,\quad c=0.514,\quad d=2.62\times
10^{-4},e=0.553.$

This is shown in Figure 11 and compared to simulation data for several values
of $q_{in}$. The RMS error on this fit across all simulation data is 0.18. It
can be seen that this fit captures both the peak to $\Delta\gamma_{e}$ at
$\nu_{e,u}^{*}\sim 20$ and that fact that the peak increases slowly with
$q_{in}$. It also drops to 0.5 for $\nu_{e,u}^{*}\lesssim 7$, but there is
some disagreement with the simulation data in this regime so precise values
may be higher. This low collisionality behaviour is somewhat speculative, but
it is reasonable to expect that $\Delta\nu_{e,u}^{*}$ will at any rate be
small at low collisionalities. This is because temperature gradients, which
appear to be a necessary condition for significant enhancement to
$\gamma_{e}$, are likely to be smaller at very low collisionalities than for
intermediate $\nu_{e,u}^{*}$. This is, however, complicated somewhat by the
increased heat flux suppression predicted at low collisionalities.

By including conditions at the target into the fits, and considering
relationships of the form $\Delta\gamma_{e}=a_{0}\prod_{i}X_{i}^{a_{i}}$,
where $a_{i}$ are fit parameters and $X_{i}$ are SOL parameters, it was found
that

$\Delta\gamma_{e}=a\left(\frac{n_{u}}{10^{20}}\right)^{b}(T_{e,u})^{c}(T_{e,t})^{d}(q_{in})^{e},$
(11)

with

$a=3.768,\quad b=-1.123,\quad c=-1.603,\quad d=-0.886,\quad e=1.556,$

is a good fit to the data, with an error of 0.10. However, it was found that
this predicts very large values of $\Delta\gamma_{e}$ when upstream
collisionality is small but temperature gradients are large. As such, it is
recommended to use (10) in regimes not studied here, despite slightly worse
agreement with the SOL-KiT data. It is also worth noting that no relationship
of this kind was found to perform better for $f_{\kappa_{e}}$ than (9).

Figure 11: Fit to $\Delta\gamma_{e}$ from equation (10) as a function of
$\nu_{e,u}^{*}$ for several values of $q_{in}$. Data from simulations with the
same $q_{in}$ is shown alongside.

The scalings presented here are straightforward to implement in a fluid model.
Because of the spatial variation in heat flux suppression observed in these
simulations (see Figure 4(b)), implementing the line-averaged quantity
$f_{\kappa_{e}}$ as a prefactor to the Spitzer-Härm conductivity may not yield
accurate temperature profiles, but should be adequate for predicting the
overall power balance and target temperatures. These scalings may also be used
in simple analytical SOL models such as the modified two-point model [39] or
the Lengyel model for predicting detachment onset with radiating impurities
[40].

To test these relationships, we have implemented (9) and (10) in the fluid
version of SOL-KiT, using self-consistent values of $\nu_{e,u}^{*}$, to
calculate the modifications to
$\gamma_{e}=\gamma_{e}^{fluid}+\Delta\gamma_{e}$ and
$q_{\parallel,e}=f_{\kappa_{e}}q_{\parallel,e}^{SH}$. In Figure 12(a), we
compare the target temperatures in kinetic simulations with those in fluid
simulations, with and without the kinetic corrections. Temperature profiles
for a particular simulation are shown in Figure 12(b). Agreement with
$T_{e,t}$ is good, with the RMS error reduced from 43.3% to 8.2% with the
addition of the kinetic corrections. Temperature profiles show that agreement
is improved at the upstream and target locations, but differences exist in the
rest of the domain and the $T_{e}$ uptick close to the wall is not captured.
This is expected due to the line-averaging of $f_{\kappa_{e}}$.

If we assume these scalings can be applied for larger values of $T_{e,u}$ and
$Ln_{u}$, then for the plasma profiles obtained with the ITER scenario
modelled in [26] (‘standard transport’ case; profiles shown in Figure 13(b)),
Figure 13(a) shows the expected values of kinetic factors as a function of
radial distance at the outer midplane. Here, $q_{in}=800$MWm-2 at the
separatrix has been assumed, with a radial decay width of $\lambda_{q}=3$mm.
The collisionality ranges from $\nu_{e,u}^{*}\sim 6-25$. We see values of
$\Delta\gamma_{e}$ up to 2.59 at the peak, representing an enhancement of
around 50% over the classical value, and over 50% heat flux suppression just
beyond the separatrix. It would be straightforward to implement (9) and (10)
in a code such as SOLPS-ITER to explore the significance of such effects.

(a) Target electron temperatures in kinetic simulations (x-axis) compared with
fluid / fluid with kinetic corrections.

(b) Fluid, kinetic and fluid with kinetic corrections temperature profiles for
a simulation with $q_{in}=64$MWm-2, $Ln_{u}=0.3\times 10^{21}$m-2.

Figure 12: Adding kinetic corrections to SOL-KiT fluid mode using equations
(10) and (9).

As an additional test of how these scaling relationships behave in reactor
relevant conditions, the scaling equations were implemented into a two-point
model (2PM) [2] numerical solver, with $f_{\kappa_{e}}$ and $\Delta\gamma_{e}$
scalings incorporated as modifications to the $\kappa_{e}$ and $\gamma$ terms.
This modified 2PM is solved iteratively to find converged self-consistent
solutions with the kinetic factors, and applied to SPARC-like SOL conditions
[41] covering a $\nu_{e,u}^{*}$ range of $\sim$0.2-100. The low collisionality
limiting behaviour for the scalings described above was observed with this
tool. Significantly steeper temperature gradients, with higher $T_{e,u}$ and
lower $T_{e,t}$, were observed in the 2PM solutions at intermediate and low
collisionality, as a result of the added $f_{\kappa_{e}}$ and
$\Delta\gamma_{e}$ factors. Full results of this 2PM study will be published
elsewhere.

A caveat to the fits presented is that the plasma model in SOL-KiT does not
currently include flux tube expansion (or other SOL geometry effects) or
contributions from molecules or impurity species. The former will alter plasma
behaviour, in particular with respect to detachment, and there may also be a
kinetic effect due to the mirror force on the electrons. The latter will
represent additional particle, momentum and energy sources/sinks. An ongoing
project to redevelop SOL-KiT with a more flexible physics model, as well as
improved computational efficiency and parallelisation, should make it possible
to study kinetic effects in the presence of such additional physics.

(a)

(b)

Figure 13: Radial profiles at the outer midplane of (a) kinetic factors and
(b) plasma profiles for the ITER scenario modelled in [26].

##  7 Conclusion

We have presented kinetic studies of electron transport in tokamak scrape-off
layer plasmas across a range of input powers and densities, under steady-state
conditions.

One of the primary aims of this study has been to validate the local
approximation in fluid models, which are frequently employed to model SOL
plasmas. We see that, for SOL equilibria, a kinetic treatment of the electrons
does not change qualitative behaviour in terms of the particle flux to the
target with this plasma model, as shown in Figure 3, despite changes to the
electron temperature profiles and reductions in the target temperatures
(Figure 4).

Typically, it has been assumed that the classical value of $\gamma_{e}$ is
valid at equilibrium, but here there are differences of up to 98%, as shown in
Figure 7. We provide a qualitative understanding of this enhancement in terms
of an imprint of the fast electrons from upstream on the electron distribution
at the target. The presence of this enhanced tail is predicted to have
significant impacts on collision rates for inelastic processes with threshold
energies $\gtrsim$ 50eV, for example the ionisation of plasma impurities. This
is the subject of an ongoing study.

The enhancement of $\gamma_{e}$ and reduction in $q_{\parallel,e}$ at
equilibrium is shown to follow scalings based on basic SOL parameters, (9) &
(10). The performance of these fits is shown in Figures 10 and 11. To test the
ability of these corrections to capture kinetic effects in SOL simulations, we
have shown that implementing them in the fluid version of SOL-KiT does improve
agreement with the fully kinetic $T_{e}$ profiles. While there are caveats to
the use of these scalings outside of SOL-KiT simulations, particularly in
relation to the aspects of SOL physics not included in the model used here, it
does suggest it is viable to capture kinetic effects at equilibrium in studies
of future devices, either in fluid codes or reduced analytical models.
Extrapolating to the ITER tokamak for example does predict significant kinetic
effects, suggesting at the very least that further study into non-local
parallel transport in reactor-class tokamaks is warranted.

To the authors’ knowledge, this is the first attempt at providing relatively
simple scaling laws for kinetic effects in equilibrium SOL plasmas. The
principle that kinetic effects are directly related to basic descriptors of
SOL conditions is a potentially useful approach to analysing their
significance.

The modifications to $\gamma_{e}$ and $q_{\parallel,e}$ in conjunction with
good agreement in power balance and target particle flux behaviour (discussed
in Section 5 and shown in Figure 3), suggest that both effects contribute in a
way which approximately cancels. As such, attempts to capture kinetic effects
in fluid models should treat both phenomena simultaneously.

In this study, the changes to $\gamma_{e}$ and $q_{\parallel,e}$ are in
contrast to the behaviour at equilibrium found using the PIC code BIT1 in
[42]. There, $\gamma_{e}$ is found to be well-approximated by the classical
value, and $q_{\parallel,e}$ is a non-monotonic function of collisionality,
which is contradicted by (9). There are significant differences in the
simulations carried out in [42], in particular that plasma-neutral
interactions were neglected and that only attached regimes were studied.
Furthermore, the differences in $\gamma_{e}$ seen here are of a similar
magnitude to those seen with the KIPP code in [6].

It should also be noted that this investigation has been done for equilibrium
plasma conditions. For the sheath boundary condition in particular, much
stronger kinetic effects may be present in transient regimes as shown in [3,
42], albeit for short durations relative to inter-ELM equilibria.

## Acknowledgements and supporting data

This work was part funded by the UK Engineering and Physical Science Research
Council (EPSRC) (grant number EP/W006839/1) and the UK Atomic Energy Authority
(UKAEA). It has also been improved through informal discussions with
colleagues at UKAEA and elsewhere, in particular David Moulton and Mike
Kryjak.

The simulations were carried out using the Imperial College London High
Performance Computing resource. The simulation data used in this study, along
with the analysis scripts, can be found at https://doi.org/10.14469/hpc/10979.

## References

  * [1] S I Braginskii.  Transport Processes in Plasmas.  Reviews of Plasma Physics, pages 205–311, 1965.
  * [2] P. C. Stangeby.  The Plasma Boundary of Magnetic Fusion Devices.  Plasma Physics and Controlled Fusion, 43(2):223–224, 2001.
  * [3] S. Mijin, F. Militello, S. Newton, J. Omotani, and R. J. Kingham.  Kinetic and fluid simulations of parallel electron transport during equilibria and transients in the scrape-off layer.  Plasma Physics and Controlled Fusion, 62(9):095004, 2020.
  * [4] M. Wigram, C. P. Ridgers, B. Dudson, J. P. Brodrick, and J. T. Omotani.  Incorporating nonlocal parallel thermal transport in 1D ITER SOL modelling.  2020\.
  * [5] Menglong Zhao, A. V. Chankin, and D. P. Coster.  Implementation of an inelastic collision operator into KIPP-SOLPS coupling and its effects on electron parallel transport in the scrape-off layer plasmas.  Contributions to Plasma Physics, 59(7):1–13, 2019.
  * [6] A. V. Chankin, G. Corrigan, and A. E. Jaervinen.  Assessment of the strength of kinetic effects of parallel electron transport in the SOL and divertor of JET high radiative H-mode plasmas using EDGE2D-EIRENE and KIPP codes.  Plasma Physics and Controlled Fusion, 60(11), 2018.
  * [7] J. P. Brodrick, R. J. Kingham, M. M. Marinak, M. V. Patel, A. V. Chankin, J. T. Omotani, M. V. Umansky, D. Del Sorbo, B. Dudson, J. T. Parker, G. D. Kerbel, M. Sherlock, and C. P. Ridgers.  Testing nonlocal models of electron thermal conduction for magnetic and inertial confinement fusion applications.  Physics of Plasmas, 24(9), sep 2017.
  * [8] D. Tskhakaya.  Kinetic Modelling of the Detached Divertor Plasma.  (633053), 2015.
  * [9] J. T. Omotani and B. D. Dudson.  Non-local approach to kinetic effects on parallel transport in fluid models of the scrape-off layer.  Plasma Physics and Controlled Fusion, 55(5), 2013.
  * [10] O. V. Batishchev, M. M. Shoucri, A. A. Batishcheva, and I. P. Shkarofsky.  Fully kinetic simulation of coupled plasma and neutral particles in scrape-off layer plasmas of fusion devices.  Journal of Plasma Physics, 61(2):347–364, 1999.
  * [11] Z. Abou‐Assaleh, M. Petravic, R. Vesey, J. P. Matte, and T. W. Johnston.  Non‐Local Transport in a Tokamak Plasma Divertor with Recycling.  Contributions to Plasma Physics, 34(2-3):175–179, 1994.
  * [12] S. Mijin, F. Militello, S. Newton, J. Omotani, and R. J. Kingham.  Kinetic effects in parallel electron energy transport channels in the Scrape-Off Layer.  Technical report, 2019.
  * [13] D. Power, S. Mijin, F. Militello, and R. J. Kingham.  Ion–electron energy transfer in kinetic and fluid modelling of the tokamak scrape-off layer.  The European Physical Journal Plus, 136(11):1–13, 2021.
  * [14] S. Mijin, A. Antony, F. Militello, and R. J. Kingham.  SOL-KiT—Fully implicit code for kinetic simulation of parallel electron transport in the tokamak Scrape-Off Layer.  Computer Physics Communications, 258:107600, 2021.
  * [15] I. P. Shkarofsky, M. P. Bachynski, and T. W. Johnston.  The Particle Kinetics of Plasmas.  Reading, Mass.; Dordrecht printed, 1966.
  * [16] Lyman Spitzer and Richard Härm.  Transport phenomena in a completely ionized gas.  Physical Review, 89(5):977–981, 1953.
  * [17] R. K. Janev, D Reiter, and U Samm.  Collision Processes in Low-Temperature Hydrogen Plasmas.  2003\.
  * [18] A. Kramida, Y. Ralchenko, and J. Reader.  Nist atomic spectra database, 2020.
  * [19] Richard J. Procassini and Dana A. Knoll.  Kinetically motivated boundary conditions for fluid models of scrape-off layer transport.  Journal of Nuclear Materials, 1992.
  * [20] M. Blommaert, W. Dekeyser, N. Horsten, P. Börner, and M. Baelmans.  Implementation of a consistent fluid-neutral model in SOLPS-ITER and benchmark with EIRENE.  Contributions to Plasma Physics, 58(6-8):718–724, 2018.
  * [21] B. Dudson, J. Allen, T. Body, B. Chapman, C. Lau, L. Townley, D. Moulton, J. Harrison, and B. Lipschultz.  The role of particle, energy and momentum losses in 1D simulations of divertor detachment.  Plasma Physics and Controlled Fusion, 61(6), 2019.
  * [22] Eric Todd Meier.  Modeling Plasmas with Strong Anisotropy, Neutral Fluid Effects and Open Boundaries.  (August 2011):213, 2011.
  * [23] N. Horsten, W. Dekeyser, G. Samaey, and M. Baelmans.  Assessment of fluid neutral models for a detached ITER case.  Nuclear Materials and Energy, 12:869–875, 2017.
  * [24] P. Helander, S. I. Krasheninnikov, and P. J. Catto.  Fluid equations for a partially ionized plasma.  Physics of Plasmas, 1(10):3174–3180, 1994.
  * [25] N. Horsten, W. Dekeyser, G. Samaey, and M. Baelmans.  Comparison of fluid neutral models for one-dimensional plasma edge modeling with a finite volume solution of the Boltzmann equation.  Physics of Plasmas, 23(1), 2016.
  * [26] Irina Veselova, Elizaveta Kaveeva, Vladimir Rozhansky, Ilya Senichenkov, Anastasia Poletaeva, Richard A. Pitts, and Xavier Bonnin.  SOLPS-ITER drift modelling of ITER burning plasmas with narrow near-SOL heat flux channels.  Nuclear Materials and Energy, 26:100870, 2021.
  * [27] Giulio Rubino, R. Ambrosino, G. Calabrò, V. Pericoli Ridolfini, and B. Viola.  Comparative analysis of the SOL plasma in DEMO using EDGE2D/EIRENE and TECXY codes.  Nuclear Materials and Energy, 12:864–868, 2017.
  * [28] Michael R K Wigram.  Modelling tokamak power exhaust and scrape-off-layer thermal transport in high-power fusion devices.  PhD thesis, 2019.
  * [29] M. Zhao, A. V. Chankin, and D. P. Coster.  Kinetic simulations of electron heat flux in the scrape-off layer.  Nuclear Materials and Energy, 12:819–824, 2017.
  * [30] M. Kryjak, B. Dudson, D. Power, S. Mijin, and C. Ridgers.  Sensitivity of scrape-off layer codes to modelling approaches.  48th EPS Conference on Plasma Physics, EPS 2022, pages 4–7, 2022.
  * [31] Xian Zhu Tang and Zehua Guo.  Sheath energy transmission in a collisional plasma with collisionless sheath.  Physics of Plasmas, 22(10), 2015.
  * [32] Yanzeng Zhang, Yuzhi Li, Bhuvana Srinivasan, and Xian-Zhu Tang.  Resolving the mystery of electron perpendicular temperature spike in the plasma sheath.  033504(February), 2022.
  * [33] Ivona Vasileska and Leon Kos.  Time-Dependent Boundary Conditions During ELMs in ITER Plasma.  Journal of Fusion Energy, 39(5):212–220, 2020.
  * [34] D. Tskhakaya, S. Jachmich, T. Eich, and W. Fundamenski.  Interpretation of divertor Langmuir probe measurements during the ELMs at JET.  Journal of Nuclear Materials, 415(1 SUPPL):S860–S864, 2011.
  * [35] D. Del Sorbo, J. L. Feugeas, Ph Nicolaï, M. Olazabal-Loumé, B. Dubroca, S. Guisset, M. Touati, and V. Tikhonchuk.  Reduced entropic model for studies of multidimensional nonlocal transport in high-energy-density plasmas.  Physics of Plasmas, 22(8), 2015.
  * [36] G. P. Schurtz, P. D. Nicolaï, and M. Busquet.  A nonlocal electron conduction model for multidimensional radiation hydrodynamics codes.  Physics of Plasmas, 7(10):4238–4249, 2000.
  * [37] Jeong Young Ji, Eric D. Held, and Carl R. Sovinec.  Moment approach to deriving parallel heat flow for general collisionality.  Physics of Plasmas, 16(2), 2009.
  * [38] D. Tskhakaya, F. Subba, X. Bonnin, D. P. Coster, W. Fundamenski, and R. A. Pitts.  On kinetic effects during parallel transport in the SOL.  Contributions to Plasma Physics, 48(1-3):89–93, 2008.
  * [39] P. C. Stangeby.  Basic physical processes and reduced models for plasma detachment.  Plasma Physics and Controlled Fusion, 60(4), mar 2018.
  * [40] L. Lengyel.  Analysis of Radiating Plasma Boundary Layers.  Technical report, IPP, 1981.
  * [41] A. Q. Kuang, S. Ballinger, D. Brunner, J. Canik, A. J. Creely, T. Gray, M. Greenwald, J. W. Hughes, J. Irby, B. Labombard, B. Lipschultz, J. D. Lore, M. L. Reinke, J. L. Terry, M. Umansky, D. G. Whyte, and S. Wukitch.  Divertor heat flux challenge and mitigation in SPARC.  Journal of Plasma Physics, (2020), 2020.
  * [42] D. Tskhakaya, K. Matyash, R. Schneider, and F. Taccogna.  The particle-in-cell method.  Contributions to Plasma Physics, 47(8-9):563–594, 2007.